Vibrational spectroscopy simulation of solvation effects on a G-quadruplex

Abstract It is commonly believed that solvation effects on the vibrational properties of a solute are easily accounted for by simple rules of thumbs, that is, solvating a polar molecule in a polar medium has the only effect of red shifting all its spectroscopical features and, similarly, solvating a polar molecule in a nonpolar medium has the opposite effect. In this work, we use theoretical vibrational spectroscopy at quasi-classical and quantum approximate semiclassical level to gain atomistic insights about solvent–solute interactions for 2′-deoxyguanosine and the G-quadruplex. We employ the quasi-classical trajectory method to include full anharmonicity into our calculated spectra, and then introduce quantum nuclear effects by means of divide-and-conquer semiclassical spectroscopy calculations. Solvation is treated explicitly leading to a good reproducibility of the available experimental data and reliable predictions when an experimental reference is missing. Communicated by Ramaswamy H. Sarma


Introduction
The nowadays well-known double helix structure of DNA is widely regarded as one of the most important discoveries of the last century.This kind of structure is characterized by a well-defined pattern of interactions between pairs of nucleobases and it is able to preserve the entire genetic information of living creatures.It is also now known that some nucleobases are able to interact with each other in a different way than the one described by Watson and Crick (WC) (Hoogsteen, 1963).To emphasize this concept, some computational studies demonstrated that guanine seems particularly prone to interact with other guanine molecules (Kelly et al., 2005;Sieranski, 2020;Sponer et al., 2013).In many living organisms we are indeed able to observe this behavior, for instance in the guanine rich sequences found in the telomeric and promoter regions of DNA (Lipps & Rhodes, 2009).It is likely that these sequences fold into a noncanonical double-helix WC structure, where four guanines interacting through Hoogsteen-type bonding are arranged in a square planar geometry.This type of geometry that is formed by one, two or four DNA strands, is also known as G-quadruplexes.Factors like the sequence of nucleobases, presence of certain ligands and other environmental effects determine the way the quadruplex will fold.For this reason, many shapes of this kind of molecule can be found.Recently, a lot of attention has been payed by the scientific community to these kinds of systems due to their potential functional role in gene regulation and DNA repairing mechanisms from lesions (Guiblet et al., 2021;Rhodes & Lipps, 2015).
The interest in G-quadruplexes comes mainly from their potential applicability as targets of cancer therapeutics (Adam et al., 2002;Cogoi & Xodo, 2006;Hansel-Hertsch et al., 2017;Mohanty et al., 2021;Rankin et al., 2005), but also for their possible role in virus replication activity for viruses like HIV and SARS-CoV-2 (Liu et al., 2022;Musumeci et al., 2015;Perrone et al., 2013;Ruggiero et al., 2021).This is due to the fact that G-quadruplexes serve as kind of "knot" in the DNA replication process signaling the stop of transcription.This means that the development of ligands stabilizing the quadruplexes could be used to induce the apoptosis of target cells (Neidle, 2016).Many efforts have been put also in the in vitro detection of specific G-quadruplexes by means of UV spectroscopy and Circular Dichroism (CD) (Mergny et al., 1998;Paramasivan et al., 2007).Other techniques used for the analysis of the quadruplexes in water solutions, which are also good candidates for the in vivo detection, are Raman spectroscopy (RS) methods such as surface enhanced Raman spectroscopy (SERS), polarized Raman scattering and UV resonant Raman (UVRR) (Di Fonzo et al., 2020;Friedman & Terentis, 2016;Pagba et al., 2010).While all of these approaches focus mainly on the effects of topology change, change of the templating cation and signal identification, none of them seem to point towards solvent effects, which are fundamental to be taken into consideration before moving to the in vivo analysis.
In this work, our main goal is to understand the effects of solvation on the vibrational features of the G-quadruplex.In particular, our work focuses on the differences that occur between gas phase spectroscopic simulations and simulations carried out in a box of solvent.The latter should also allow us to demonstrate that use of an appropriate explicit solvent is fundamental to catch the correct shifts caused by system-solvent interactions.
Vibrational spectroscopy is a useful tool often employed by experimental chemists and physicists to understand the structure and dynamics of molecules.Techniques such as the Fourier transform IR (FTIR) or more recently the IR multi photon dissociation (IRMPD) are able to give many important insights about molecules in gas phase, solid state and, in some cases, even liquid state (Asami et al., 2009;Choi & Miller, 2006;Fraschetti et al., 2014;Ivanov, 2010;Ivanov et al., 2015;2019;Mondragon-Sanchez et al., 2004;Oh et al., 2005;Setvin et al., 2015;Wu & McMahon, 2007;Xu et al., 2011).However, their applicability to biological systems still remains problematic.This is due to the fact that the absorption of liquid water in the IR is very intense, causing the spectrum to be entirely dominated by the solvent features instead of those of the system.RS is ideal for this kind of applications, since liquid water has a low Raman activity.For these reasons RS has found large application to both chemistry, biochemistry and also the biotechnological field (Friedman & Terentis, 2016;Gualerzi et al., 2021).However, while all these techniques are able to give useful information about structure, dynamics and presence of a given molecular species, they are not able to catch the atomic details for any of these.
In the last years computer simulations have served as a tool to interpret and predict experiments.The work in our group has been mostly focused on the development and application of path-integral-derived semiclassical techniques for the simulation of vibrational spectroscopy (Conte & Ceotto, 2020) at the (approximate) quantum mechanical level.In particular, very recently, we have also shown that application of our methods to large systems is possible even adopting empirical force fields (Gabas et al., 2020(Gabas et al., , 2022)).
In this article, we focus on biological systems.We start with a simple system like 2 0 -deoxyguanosine, which we study both in gas and solvated phase.We adopt the information gained from this first set of simulations to show how solvation acts differently when the same molecule is found in a quadruplex arrangement, and compare our results with experimental data.

Theoretical and computational details
The reported simulations are characterized by more than 20,000 degrees of freedom.Therefore, the only viable way to describe the potential energy surface (PES) of the system is the adoption of classical force fields (FF).In a previous work presented by three of us (Gabas et al., 2020), we benchmarked several FF against ab initio data, and the polarizable FF AMOEBABIO18 included in the Tinker 8.6.1 software suite emerged as the most effective one (Rackers et al., 2018;Zhang et al., 2018).An important point to make while using FF applied to spectroscopic simulations is that as the quality of the PES description is limited, so it is expected to be the accuracy with respect to experimental data.However, a qualitative and insightful description of the spectroscopic features is still achievable.For every studied system, a preliminary geometry optimization is carried out in gasphase conditions, with a quasi-Newton algorithm.Then, the system is inserted in a solvent box with periodic boundary conditions (PBC), and particle mesh Ewald (PME) scheme for the long range interactions is adopted.Choice of the dimension and shape of the box is specific for the system and will be specified later in this section.Once the molecule is inserted in the PBC box, a second optimization with the steepest descent method is carried out.The molecular dynamics (MD) simulations are propagated for 3000 steps with a time step of 0.2 fs by using the Beeman integrator (Beeman, 1976).A particular set of initial mass-scaled conditions is used for the NVE trajectories p j ð0Þ ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi � hx j ð2n j þ 1Þ p q j ð0Þ ¼ q j, eq , ( where the initial position for the jth normal mode is the one found for the equilibrium configuration, the initial linear momentum is obtained from the harmonic oscillator (HO) approximation, x j is the harmonic frequency for the jth mode and n j is a positive integer number which represents the quantum of excitation for the jth mode.With the initial conditions shown in Eq. ( 1), it is possible to run trajectories at an energy close to the harmonic zero point energy (ZPE).By adding a quantum of excitation on a given mode, the energy of the trajectory will be closer to the energy of the quantum state corresponding to the excitation in that mode.

Quasiclassical trajectory
Most of the spectroscopic simulations we present are made at the quasi classical trajectory (QCT) level, where the initial conditions of the trajectory are quantized (Eq. 1) and the evolution of the system is purely classical.In this framework, the molecular power spectrum is calculated for each vibrational mode j as the Fourier transform (FT) of the classical momentum (or velocity) autocorrelation function.
In QCT spectra, nuclear quantum effects are not accounted for, meaning that the spectral peaks correspond only to classical vibrational frequencies.Our approach with this method is to run a single trajectory with a quantum of excitation in the mode of interest.Usually, a trajectory per mode of interest is sufficient to obtain a clean and accurate enough spectrum, but as it will be shown in the results section, for some modes a wider area of the phase space needs to be sampled to achieve a better qualitative description of the spectroscopic feature.In the latter instance, we use a set of trajectories with randomized initial mass-scaled conditions starting from those in Eq. (1).p j ð0Þ ¼ À ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi � hx j ð2n j þ 1Þ p sinðh j Þ q j ð0Þ ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where h j is a random uniformly distributed phase between 0 and 2p, n j is set equal to 1 for the modes of interest and equal to 0 for all other modes.The resulting spectrum when several MD runs are performed is simply the average of the single spectra, that is, where I ðiÞ j ðxÞ is the intensity at frequency x for the jth normal mode given by the i-esime trajectory.

Semiclassical method
Semiclassical initial value representation (SCIVR) has proved to be a practical and accurate tool for adding quantum effects to molecular dynamics (Begusic et al., 2022;Begusic & Vanicek, 2020, 2021;Bonnet, 2020Bonnet, , 2021;;Church et al., 2017Church et al., , 2018;;Church & Ananth, 2019;Grossmann, 2006;Heller, 1981;Kay, 2006;Larregaray & Bonnet, 2021;Malpathak et al., 2022;Miller, 1974Miller, , 2001Miller, , 2005;;Pollak, 2007;Vanicek & Begusic, 2021).Several semiclassical simulations are carried out also by means of the multiple coherent divide-and-conquer semiclassical initial value representation method (MC-DC SCIVR; Aieta et al., 2020;Bertaina et al., 2019;Ceotto et al., 2009Ceotto et al., , 2017;;Di Liberto et al., 2018;Gabas et al., 2018Gabas et al., , 2019;;Gandolfi et al., 2020).With this method, it is possible to compute quantum mechanical vibrational power spectra of high dimensional systems upon projection of the dynamical quantities onto lower dimensional vibrational subspaces.In this case the working formula is where ðp i ð0Þ, qi ð0ÞÞ is the initial phase space configuration projected on the subspace for the ith trajectory, St is the projected classical action evaluated along the trajectory, Ẽ is a parameter and /t is the phase of the projected prefactor, which is defined as ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi To evaluate the sub blocks of the monodromy matrix which appear in / t it is necessary to calculate the Hessian matrix at every time step of the MD simulation.This represents a huge bottleneck for the whole method.The problem of Hessian calculations has been tackled with different methods in our group (Ceotto et al., 2013;Conte et al., 2019;Gandolfi & Ceotto, 2021).Ab initio methods applied to molecules with more than 50 atoms can be seriously problematic under this aspect, while when using an FF the Hessian matrix can be computed at every step without approximations even for molecules with more than 100 atoms.For systems bigger than 1000 atoms, we opt to propagate along the trajectory only the elements of the Hessian related to the normal modes belonging to the subspace of interest.These elements of the Hessian matrix are approximated by differentiation of the forces calculated along the dynamics (Cazzaniga et al., 2020(Cazzaniga et al., , 2022)).With this method, we are able to focus only on the normal modes centered on the molecule and eliminate from the spectrum those coming from the solvent molecules.For SC calculations, a single trajectory has been used upon setting to zero the initial momentum of all the normal modes with harmonic frequencies lower than 300 cm À 1 , and providing all other normal modes with the basic zero-point energy contribution (i.e., n ¼ 0 in Eq. 1).

Computational details
As anticipated in the introduction, we focus our attention on the study of the 2 0 -deoxyguanosine, and a tetrad of 2 0 -deoxyguanosine (Figure 1) which from now on will be simply referred to as the 'quadruplex'.Guess geometries for both systems are extrapolated by removing the necessary residues from a DFT optimized structure published in a paper by Sponer et al. (2013).Then, minimization is performed for a single deoxyguanosine inserted in a 30-side cubic box filled with water (2702 atoms), and the quadruplex inserted in a rectangular box with a 50 � 50 2 base, and 30 height (7517 atoms).For both systems, long range interactions are considered using PME (Darden et al., 1993;Essmann et al., 1995).As a comparison with experimental data, we considered the Raman spectrum in water of the thrombin binding aptamer (TBA) by Pagba et al. (2010).In this case, we made the assumption that the quadruplex considered in our simulations can be thought of as a minimal sub-unit of the full TBA.Other experiments with other types of G-quadruplexes exist as well, but we think that since TBA is one of the smallest systems that can be found in literature it could better resemble the behavior of the minimal unit that we are considering.We are not currently aware of how much the presence of a rigid backbone like the one present in the TBA can influence our discussion, but this factor will be taken into consideration in the result discussion part.The vibrational normal modes considered in this work, for both studied systems are the asymmetric and symmetric stretches of the C2 0 H 2 group of the sugar moiety, the C6 ¼ O stretch, and the 2 0 -deoxyguanosine (dG) ring deformation mode.

0 -Deoxyguanosine
In the first set of QCT simulations performed on the 2 0 -deoxyguanosine molecule, we observe that solvation acts differently on different normal modes.Specifically, it is possible to observe that for the C ¼ O stretch (Figure 2, panel a), multiple vibrational features are present in both solvated and gas phase systems.In order to assign the C ¼ O stretch to the correct peak, a more detailed analysis has been conducted.For both systems, two additional trajectories are run.In one, all the initial momenta are set to the harmonic ZPE (ZPE trajectory).In the other one the initial momentum of the C ¼ O stretching mode is set to zero, while all the other modes have the same initial momenta as in the ZPE trajectory (deactivated trajectory).By doing so, in the gas phase spectrum (Figure s1 panel a in the supplementary information file) we observe for the deactivated trajectory a single peak much more intense than the one observed in the original trajectory, which we recall was characterized by a harmonic quantum of excitation in the stretch and the corresponding harmonic zero-point energy contribution in all other modes.The position of this peak is � 30 cm À 1 red-shifted with respect to the most intense peak in Figure 2a.In the ZPE trajectory the same feature of the spectrum in Figure 2a is observed, even if it is slightly blue-shifted and with a change in the intensity ratio between the most intense peak and the secondary one.In particular, it seems that the secondary one gains in intensity with these set of initial conditions, while the main one loses intensity.Based on these observations we can conclude that the most intense peak can be assigned to the C ¼ O stretch because in the deactivated trajectory it provides an almost harmonic signal.Moving on to the analysis of the solvated system (Figure s1 panel b in the SI file), we observe that both the deactivated and ZPE spectra are much less intense than the one in Figure 2a.In the spectrum generated by the deactivated trajectory, a single peak is again observed.It is positioned almost at the same frequency of the secondary peak  on the right of the main one.A similar situation is observed with the ZPE trajectory, with the main difference that the observed feature is enlarged with a shoulder near the position of the most intense peak of the original spectrum (Figure 2a).Again, these two observations led us to the conclusion that the C ¼ O stretch is the one positioned at 1553 cm À 1 , which corresponds to the most intense one.With this we can conclude that a red shift is observed when solvent molecules are added, on the contrary of the harmonic estimates.
By a simple visual analysis of the MD simulation whereby the spectrum has been calculated, it is possible to observe that this part of the molecule is strongly interacting with the surrounding water molecules.We think that this interaction with water is what causes the observed red shift of the QCT signal with respect to the gas phase case.Another peculiarity that can be observed from the spectrum in Figure 2, is that while the QCT signal is reported to be strongly harmonic in gas phase, in the water box simulation the harmonic estimate is 115 cm À 1 far from the QCT signal.Probably this anharmonicity derives from the fact that in the equilibrium conformation, from which the harmonic estimate is made, the C ¼ O group sense some of the interactions with the solvent as repulsive.This could explain also why the harmonic estimates of the CO stretch in the water box system are blue shifted with respect to the harmonic estimate made in vacuo.This result demonstrates the importance of anharmonicity and dynamical effects when simulating biomolecules in solution.
The opposite effect was observed for the ring deformation mode (dG RDM).In this case solvation seems to blue shift both the QCT and harmonic estimates with respect to the signal observed in gas phase (Figure 2

panel b).
DG RDM is a vibrational normal mode delocalized on all atoms of the guanine.For this reason, a simple visual analysis of the MD simulation may be not sufficient to understand why the spectrum in Figure 2 (panel b) features blue shifted QCT and harmonic frequencies with respect to the corresponding ones in the gas phase simulation.By calculating an average (Levine et al., 2011) of the radial pair distribution function (RDF) between all the atoms of the guanine ring and all that ones of the water molecules (blue solid line in Figure 3), it is possible to observe that water tends to occupy outer solvation shells, except for a little peak near 2, which can be interpreted as the presence of labile short range interactions.This probably means that the interactions between water and the whole guanine ring tend to be unfavorable.By observing the RDF for C6 ¼ O (green solid line in Figure 3) atoms, it can also be found the proof that during the MD simulation this part of the molecule is strongly interacting with water.In particular, this can be deduced by the fact that a high peak is present at distances shorter than the first solvation shell.
The last set of normal modes that we analyze for the 2 0 -deoxyguanosine is localized on the C2 0 H 2 group on the sugar moiety of the molecule.
By observing the spectrum of the C2 0 H 2 modes (Figure 2 panel c), it is expected to be difficult to extract much information.This is due to the fact that in the QCT gas phase spectrum associated to the trajectory of the symmetric stretch there is also a signal related to the asymmetric stretch, while for a water solvated system the signals are split but do not provide a clear fingerprint.For the gas phase system, we relate the presence of several peaks to the fact that a single trajectory for these two modes is not sufficient to sample a wide enough phase-space volume around the minimum geometry.
For this reason, we move to a multiple trajectory approach, with initial conditions as stated in Eq. ( 3).In particular, we try this method with, 10, 100 and 1000 trajectories for the gas phase system, for all the normal modes to be studied.
By observing the spectra in Figure 4, it is possible to see that by increasing the number of trajectories from 10 to 1000 a better resolution of the separation between the QCT asymmetric and symmetric stretches for C2 0 H is obtained.On the other hand, for the C ¼ O stretch and the ring deformation mode (Figure s2 in the SI file) a reduced number of trajectories is already sufficient to provide a properly converged signal.This could be due to the fact that normal modes like the C ¼ O stretches and dG RDM are usually rather harmonic, meaning that a more restricted volume of phase space can be sampled to recover the correct spetroscopic feature.Given these results, we are now ready to apply our multitrajectory approach to water solvated systems.We do this by employing 100 trajectories given the unaffordable overhead of adopting 1000 trajectories for such systems.By comparing gas-phase and 'in solution' spectra (Figure 5) we are eventually able to assess the effects of solvation on the C2 0 H 2 normal modes.We find that results are very similar in the two different phases.This makes us believe that water molecules have no sizable interactions with the C2 0 H 2 group.This can be verified by looking at the RDF plots for the two trajectories in question (Figure 3 black and red lines).In these cases, it can be observed that water tends to be organized in solvation shells without any sign of short range interaction.

Quadruplex
In our vibrational analysis of the single 2 0 -deoxyguanosine molecule we observe that the effects of solvation are quite complex and not easily detectable by adopting simple rules of thumb.Moving on to a more complex system like the quadruplex (Figure 1b) will further prove this concept true.In particular, in a system like the quadruplex it is possible to ascertain that a well-defined pattern of interactions between the four 2 0 -deoxyguanosine residues is present.We will confirm if the water molecules are able to interact mainly with the external part of the system, while the core of the system is partially hindered to interact with the solvent.More than that, in a structure like the one of the quadruplex, we expect that the vibrational frequencies of the modes under investigation stick to a symmetry related order.The whole supramolecolar system which composes the quadruplex (Figure 1) should ideally belong to the D 4h symmetry group.Therefore, in an ideal scenario, the normal modes which we are considering should reveal a 4-fold degeneracy.In reality, there are small deviations from this symmetric arrangement and so symmetry-related frequencies are expected to be close to each other but not exactly degenerate.
In the experimental spectrum of TBA (Pagba et al., 2010), the signals of the C2 0 H 2 stretches are assigned to a broad low intensity feature.The main peak at 2952 cm À 1 is related to the asymmetric stretch, while the shoulder at 2890 cm À 1 to the symmetric stretch.The C ¼ O stretches are experimentally assigned to the intense and broad peak centered at 1663 cm À 1 , while the dG RDM mode is assigned to an intense narrow peak at 1582 cm À 1 : We start our vibrational analysis of the quadruplex by looking at the C2 0 H 2 normal modes.As it was observed on the previous case, a single trajectory per normal mode is not sufficient to extract fully reliable information (Figure s3 in the SI file).
For this reason, we calculate the spectrum by running 100 trajectories for both the gas and solvated system.With this approach we are able to reproduce the experimental feature.On the other hand, our results do not have a strong quantitative agreement with respect to the experimental data, resulting in a mean absolute error (MAE) of 107 cm À 1 : As anticipated in the Introduction, this is mainly due to the force field we need to employ to describe the PES of the system.
On the basis of the calculated spectra in Figure 6, it is also easy to see that even in this case the solvation of the system has nearly no effect on the C2 0 H 2 signal positions.In the case of a single 2 0 -deoxyguanosine, this part of the molecule seems to have almost no interaction with the water  molecules and for the quadruplex we think that this should not change much.In TBA, the tetrad of guanine is surrounded by a rigid backbone and the region near the C2 0 is excluded from interacting with water as well, and, for this reason, we think that considering a model where this part of the molecule is weakly interacting with water should not be a rough approximation.
Moving on to lower frequencies, we analyze the C ¼ O stretches.In our optimized structure of the quadruplex, water molecules are able to interact with the C ¼ O group of the quadruplex only from above or below and this makes this part of the molecule less prone to solvent interaction.
By observing the spectrum in Figure 7, it is clear that in this case the addition of solvent molecules to the simulation has the effect to blue shift the frequencies with respect to the gas-phase simulation.The opposite effect is observed for the single 2 0 -deoxyguanosine molecule, which means that solvation has a different effect on the vibrational features of the quadruplex.We use again the RDF as a tool to better understand our results (Figure s4 in the SI file).In this case, we find out that the shape of the function is much similar to the one of the dG RDM mode reported in Figure 3 for the single 2 0 -deoxyguanosine.This means that the interaction of the C ¼ O group with water is weak and unfavorable.The agreement with the experiment for the C ¼ O stretches is good, and the QCT signals for these modes computed for the quadruplex are also blue shifted with respect to the one computed for the single 2 0 -deoxyguanosine.This shift, as reported in the literature (Krafft et al., 2002), is a distinctive sign of the presence of the 2 0 -deoxyguanosine in a G-quadruplex arrangement instead of a single/double strand configuration and can also be observed in our computed data (Table 1).
The last set of normal modes that are left to be examined are the dG RDM modes.Our observations in this case are mostly aligned with what was observed with the single 2 0 -deoxyguanosine, that is, this set of modes are strongly harmonic, and solvation makes both the QCT and harmonic estimates blue shift.The main difference lies in the reason of the blue shift, since the water molecules are not able to get as close to the guanosine rings as in 2 0 -deoxyguanosine.
Likewise the case of the C ¼ O stretches of the quadruplex, in the dG RDM spectrum (Figure 8) we observe that both the harmonic and QCT signals are in good agreement with experimental data.The main difference between the harmonic estimates and the QCT frequencies is that the latter are nearly degenerate like the experimental ones.As anticipated, our simulation concerns only a minimal unit of the full TBA.This means that compensation of errors could be a factor in our last set of results, but it has to be noted that, like in the full structure, also the core part of the quadruplex has low capability of interaction with water molecules.Therefore, it is possible to suggest that the effect of the presence of a rigid backbone is to preserve the core part of the G-quadruplex, which means that a structure without it opens the possibility on longer time scales to see a dissociation of the supramolecular system.Following this argument, we think that the role of the water molecules in the quadruplex simulations is to constrain the four fragments of 2 0 -deoxyguanosine in the quadruplex geometry.The geometry of the system restrained in the quadruplex Table 1.QCT frequencies of the studied modes compared between the single 2 0 -deoxyguanosine and the quadruplex and to experimental data of TBA from ref. Pagba et al. (2010) form can lead to repulsive interactions with the other 2 0 -deoxyguanosines, so that it should be possible to explain the observation of blue shifts when solvent molecules are added.

Semiclassical spectra
So far, we have considered only a classical nuclear framework to sample the Born-Oppenheimer surface of the system.In this framework, the energy of the trajectory is ideally fixed during the dynamics and only the accessible configurations of the phase space will be sampled.We have employed short-time MD simulations and the sampled phase space in some cases could be insufficient.Propagating longer in time the MD simulations seems to be an appealing and easy solution, but this kind of approach brings in additional issues since more and more energy will be transferred to low frequency modes and even rigid rotations due to numerical inaccuracies and energy leakage.This would lead to a less accurate spectrum.One solution we find out to be helpful is to run more trajectories with same energy, as shown previously.Another solution is to consider an approximate quantum mechanical evolution of the system, as in the MC-DC-SCIVR method (Ceotto et al., 2009(Ceotto et al., , 2017;;Di Liberto et al., 2018).In this case, a Gaussian wavepacket is propagated in the phase space, meaning that even by means of a single trajectory a wider area with respect to a classical trajectory is sampled.Since the two low frequency modes are strongly harmonic and in good agreement with the expected results, we focus again our attention on the symmetric and asymmetric C2 0 H 2 stretches.We start again our analysis by considering the 2 0 -deoxyguanosine molecule.For the gas phase calculation to partition the system into subspaces we use the average Hessian on all the Hessians along the trajectory (Ceotto et al., 2017).For the water box system the average Hessian is calculated   starting from fifteen Hessians uniformly distributed along the trajectory.For the gas phase system we employ a bidimensional subspace in which the asymmetric and symmetric stretches are present.This means that these two normal modes are strongly coupled along the dynamics, and could also explain why only a large number of trajectories is able to reproduce a separation of the two signals.
The MC-DC SCIVR spectrum of 2 0 -deoxyguanosine in gas phase (Figure 9 panel a) shows two signals that are nearly degenerate, which, similarly to the QCT for the same molecule, is due to the fact that asymmetric/symmetric modes are strongly coupled to each other.Moving on to the water box system, we choose two monodimensional subspaces for both modes.By looking at the spectrum (Figure 9 panel b) for the asymmetric mode a double peak is present, with the most intense peak still degenerate with the one of the symmetric mode (upper part of panel b).In this case by changing the initial conditions as done with 2 0 -deoxyguanosine, we were able to assign the less intense peak in Figure 9b (black line) to the asymmetric stretch mode.We notice that by using MC-DC SCIVR to calculate the spectrum we obtain a better resolution with respect to the QCT counterpart.As far as the peak around 3350, we think that is a combination band originated by the C2 0 H 2 symmetric stretch.
Finally, we move to apply the semiclassical approach to the quadruplex employing the same protocol adopted for 2 0deoxyguanosine.
By analyzing the spectra in Figure 10, we conclude that for the gas phase there is an improvement moving from QCT (Figure s3 in the SI file) to SC single trajectory calculations.In particular, it is evident that most of the asymmetric and symmetric signals seem to be split.One problem that arise in this case is that the spectra are not always reporting a single or few peaks, due to the presence of resonances and other vibrational states strongly coupled to the C2 0 H 2 stretches, which hampers the assignment of the signals.The same is true also for the water box system, in which it is possible to observe that the experimental feature seems to be reproduced reasonably.On a quantitative level the experimental values in this case are about 100 cm À 1 apart from the theoretical estimates, which seems to be fair since the employment of 100 trajectories gave nearly the same results for this set of modes.

Conclusions
With our analysis, we were able to show how solvation acts differently on different vibrational features and compare them between the molecular and supramolecular arrangement.Solvation cannot be accounted for only by considering simple scaling factors, but it is often necessary to consider every normal mode and molecule explicitly employing a high level spectroscopy method.To corroborate our results, we employed the radial pair function, which is able to assess how much a certain part of the molecule is able to interact with the solvent during the trajectory.FF do not provide highly precise descriptions of the PES and, for this reason, the use of FF in the theoretical spectroscopy field should be limited to a qualitative comparison rather than a quantitative one.In order to combine the possibility to deal with large size molecular systems and maintain an accurate description of the PES, QM/MM methods should be employed.Relative to this matter, a series of promising QM/MM methods are the ones that implement the mutual polarization of the QM and MM part (Lipparini & Mennucci, 2021;Macaluso et al., 2021).In this work, we also showed that for some kind of normal modes a single trajectory with the 'correct' initial conditions, can lead to a good qualitative agreement with respect to the experiment.Instead for other normal modes a distribution of trajectory run in the same energy shell is needed to increase the agreement with the experiment.By using FF it is possible to run a high number of trajectories without worrying too much about computational overheads.The last point in the conclusions, comes from the fact that we were also able to observe that a quantum mechanical description of the dynamics obtained by a single trajectory, is able to give us results similar to what it was obtained by calculating the spectra by means of multiple trajectories.

Figure 1 .
Figure 1.(a) 2 0 -Deoxyguanosine equilibrium geometry in the water box; (b) quadruplex minimum geometry in the water box.

Figure 2 .
Figure 2. 2 0 -Deoxyguanosine QCT spectra with a single trajectory per normal mode.Panel (a) C ¼ O stretch, (b) ring deformation and (c) asymmetric and symmetric stretches of C2 0 H 2 : Vertical dashed lines indicate the harmonic estimates.

Figure 3 .
Figure 3. 2 0 -Deoxyguanosine radial distribution function (RDF) for the trajectory associated to the analyzed normal modes calculated between the group of atoms of interest and water.

Figure 6 .
Figure 6.Gas phase (a) and solvated (b) C2 0 H 2 stretches QCT spectra calculated by means of 100 trajectories for the quadruplex.The profile of the experimental feature is presented as a broad black line, where the asymmetric stretching is highlighted in orange and the asymmetric one in teal.

Figure 7 .
Figure 7. Quadruplex C6O stretches QCT spectrum (solid colored lines) and the relative harmonic estimates (vertical dotted lines) for the gas-phase and solvated quadruplex.

Figure 8 .
Figure 8. Quadruplex dGRDM QCT spectrum (solid colored lines) and the relative harmonic estimates (vertical dotted lines) for the ring deformation modes.The green band identifies the experiment.

Figure 9 .
Figure 9. 2 0 -Deoxyguanosine semiclassical spectra of the asymmmetric and symmetric stretching of the C2 0 H 2 group in gas (left) and water solution (right).

Figure 10 .
Figure 10.Semiclassical spectra of the quadruplex in gas phase (a) and water solution (b), for the asymmetric and symmetric stretches of C2 0 H 2 : .