DFT Study on the Conformational and Vibrational Properties of 3 ’-Deoxycytidine and Its Analogues

Impact of saturation of the sugar C(3’)=C(4’) bond of acytidine nucleoside derivative, 3’, 4’-didehydro-3’-deoxycytidine (I) is revealed usingsimulated vibrational spectra, with respect to 3’-deoxycytidine. The density functional theory based calculations found that the C(3’)=C(4’) double bondrestricts the sugar flexibility and affects the sugar-base intramolecular hydrogen bond network. For 3’-dC, two minimal energy conformers are identified on the potential energy surface. The first conformer (IIA) takes anO4’-endo and gauche-gauche (gg) orientation whereas the second (IIB) has a C3’-exo and gauche-trans (gt) orientation in gas phase. The two conformers which have been observed previously in the crystal structure are separated by a low energy barrier. The carbon double bond in the sugar moiety of I confines the pseudorotation of the pentagon ring to be significantly flatter than that of the sugar ring in IIA and IIB. The simulated vibrational spectra, both in gas phase and in solutions, report such structural caused spectral changes with red and blue vibrational frequency shifts.The firsthydration shell of 3’-dC has been also investigated applying hybrid QM/MM molecular dynamic simulation.


Introduction
Cytidine analogues are an important class of chemotherapeutic agents used mainly as anti-tumoror anti-viralagents (Galmarini, Jordheim, & Dumontet, 2003;Schinazi et al., 2002).They may interfere with certain nuclear enzymes, acting as enzyme inhibitors, or be incorporated in a growing nucleic acid chain causing cell death (Galmarini, Mackey, & Dumontet, 2002).On the other hand, resistance to current anticancer and antiviral agents are emerging and wide-spreading, which requires insight understanding of the molecular behaviour of the drugs including their detailed structures as DNA building blocks (Fernández-Calotti, Colomer, & Pastor-Anglada, 2011;Clercq, 2004).Previous studies have reported that anti-HIV activity of the nucleoside antibiotics drugs depends upon the ribose conformer (Torrance et al., 2001), usually with an unsaturated sugar ring and lack of the C2'-OH or C3'-OH groups.It has been also shown that conformational changes may impact the anti-viral or the anti-tumor activity (Rodriguez et al., 1994;Wang, 2000).This may somehow guide the process of drug design by, for example, the synthesis of conformationally locked nucleoside analogues to facilitate the interaction with their target enzymes (Rodriguez et al., 1994;Wang, 2000).
For studying the solvent effect, several solvents of various polarities, benzene ( = 2.27), tetrahydrofuran (THF,  = 7.43), tetrahydrothiophene-s,s-dioxide (THT,  = 43.96)and water ( = 78.36),are employed to investigate the changes in different solvation media and applying the implicit solvent model Polarized Continuum Model (PCM).We have also further investigated the first hydration shell of 3'-dC using explicit solvent MD simulations applying the hybrid QM/MM method and utilizing the DFTB method for the QM subsystem.Results obtained by this method range from excellent to good agreement with other highly accurate wavefunction and/or DFT based models.DFTB geometrical data are in particular with excellent agreement with higher level calculations (Yang et al., 2007).The method also provides some improvement via the introduction of third order dispersion as well as other corrections (Yang et al., 2008;Dolgonos et al., 2009;Gaus, Cui, & Elstner, 2011;Elstner et al., 2001).These modifications allow DFTB to perform very well for proton affinities calculations, DNA bases stacking and several more (Yang et al., 2007;Kaminski et al., 2012).
For hydration shell calculation, the optimized nucleoside was immersed in a cubic box of water with a buffer region from the boundary of 20 Å. QM/MM dynamics is performed using the dispersion corrected DFTB Hamiltonian for the QM system composed of 3'-dC only.The MM subsystem was treated using the flexible SPC/Fw water model as implemented in AMBER and applying the ff99SB force field (te Velde et al.; Wu, Tepper, & Voth, 2006;Lindorff-Larsen et al., 2010;Case et al., 2005;Salomon-Ferrer, Case, & Walker, 2012).
The system was first minimized then heated to 300 K in the NVT ensemble for 2 psthen relaxed classically for 2 ns in the NPT ensemble.Again the system is relaxed in the NPT ensemble but applying the DFTB Hamiltonian for the nucleoside only subsystem and for 750 ps and extended by another 1 ns production simulation.In all stages, periodic boundary conditions are applied and a small time step of 0.5 fs was utilized to integrate the equations of motion.For long range electrostatic interactions, particle mesh Ewald (PME) (Darden, York, & Pedersen, 1993) was utilized and truncated after 15 Å.The DFTB Hamiltonian is known to give excellent structural results for nucleosides (Kubař et al., 2007).

Geometric Properties
Studying geometrical properties of nucleosides and their analogues is important, as this may affect the interactions with their macromolecular targets (Harte Jr. et al., 1991;Choi et al., 2003).This may have certain implications on some fields such as drug design.Table 1 compares the globally minimized geometric parameters of I and II with available results obtained from experiment and previous studies (Wang, 2007;Selvam, Chen, & Wang, 2010;Karthe et al., 1997).The results are in good agreement with the available computational data and experimental crystal structures.As commonly observed in most pyrimidine nucleosides and their analogues, the minimal structures are generally belonging to the anti category of conformation, which refers to the orientation of the sugar ring relative to the base moiety, and is defined by the  angle (= O(4′)-C(1′)-N(1)-C( 2)).The Potential Energy Surface (PES) of each nucleoside has been also explored for other conformational structures that may also be populated.For I, the second most stable conformer was belonging to the syn category of conformers ( = 59.90°).For II, two favourable conformational structures have been identified with a total energy difference of 5.675 kJ•mol -1 (Figure 1b).Indeed, the X-ray crystal structure of 3'-dC has already identified the existence of two conformers, the first one was belonging to the gauche-gauche (gg) category and the second one was belonging to the gauche-trans (gt) category of conformations.The gauche and trans notations refer to the orientation of the exocyclic C(5')-O(5') arm, and this is defined by two angles; the ω (O(5′)-C(4′)-C(5′)-O(5′)) and the (C(3′)-C(4′)-C(5′)-O(5′)) angles.In the experimental paper, the authors used the letters A and B to define the gg and the gt conformational structures, respectively, and we will follow the same convention, so, IIA is the gg conformer and IIB is the gt conformer.For I, we will limit our discussion to the globally identifiedanti conformer.For comparison purposes, we will use bracket to denote the IIB properties.
As shown in Table 1, the double bond C(3')=C(4')of 3',4'-D3C (I) makes the sugar moiety more rigid and flat compared to the sugar of II.For example, I exhibits an unusual C1'-endo conformation with the pseudorotationalphase angle P of 322.89° and a puckering amplitude  max of 10.41 whereas the saturated derivative, 3′-dC (II), possesses an O4'-endo (C3'-exo) conformations with the pseudorotational phase of 73.12° (182.91°) and the amplitude of  max of 38.23 (27.05).
For II and as mentioned before, the X-ray crystal structure study identified two conformers to co-exist in the same crystal, A and B. In the present study, the two conformers are indeed located, with the gg conformer (IIA) being more stable by 5.675 kJmol -1 than the gt conformers (IIB) (Table 1, Figure 1b).The X-ray structure was obtained from the Cambridge Structural Database (CSD) (van, 2006).B3LYP results are cross-checked against MP2 results at the same level of theory (6-311++G**).Each conformer is stabilized by a number of intramolecular H-bonds (see above) as shown in Figure 1b.Some discrepancies between the experimentally determined crystal structures and the gas phase structures are expected and may be linked to the presence of intermolecular interactions and crystal packing forces.As noted previously (Altona & Sundaralingam, 1972), such packing forces may overcome certain geometrical disadvantages and/or shift the equilibrium geometry in gas phase and/or solution to another form that has a lower incidence of occurrence in crystal.For example, the experimental  angle of the gg conformer is given by -142.4° which is almost 34° (ignoring signs) less than the corresponding value in gas phase).
Both conformers (IIA and IIB) possess a pure anti configuration in gas phase/crystal structure with the angle of about -176° and -174.14° in gas phase and -142.4° and -172.10° in the crystal structure, respectively.However, the discrepancy in the stabilization energy (‫|‬ΔTE‫|‬ of 5.675 kJ•mol -1 ) is attributed to the differences of the orientation of the exocyclic arm (C 4' -C 5' ) and is estimated using the  angle.In contrast to the IIA conformer which has a gauche conformation with a angle of 54.01° in gas phase and 60.7° in crystal, IIB possess a trans configuration with a angle of 179.06° in gas phase and 173.70° in crystal and is less stabilized.For the ω angle, both molecules have a comparable values and lie within the same category (gauche).The use of solvent does not affect either the geometries or the stability orders of I or both the II conformers.For more information see supplementary materials (S1 and S2 tables).

Other Molecular Properties
Hirshfeld charges (Q H ) are important indicators of atomic site dependent activities.Figure 1a  All the nitrogen and oxygen sites exhibit negative charges, which serve as electron donor to attract H atom to form H-bonds.
For the two 3'-dC conformers, it appears that no significant changes are taking place and the most affected atomic site by this local conformational change is the O(5') site, which possess a |0.009 a.u.| more negative charge on IIA than in IIB.2012), and the more active this site the easier the activation and hence, the more potent the nucleoside.There's not much difference in the activity between both 3'-dC conformers.

Vibrational Spectroscopy
Table 2 compares simulated vibrational frequencies of I and IIA in vacuum using the B3LYP/6-311++G**model.In this table, the vibrational frequencies of IIA previouslycalculated using the B3LYP/6-31G* model by (Selvam, Chen, & Wang, 2010) are also listed.It has been demonstrated that when combined with a moderate basis set such as 6-31G*, the B3LYP functional could provide accurate vibrational frequencies of 2'-dC with experimental measurement.All wave numbers are scaled by 0.9613 (Selvam, Chen, & Wang, 2010).The calculated vibrational frequencies of IIA using the same B3LYP functional but larger basis set agree well with those of (Selvam, Chen, & Wang, 2010).The small discrepancies between the two models may attribute to the different basis set that is used.For example, the C(2')-H stretch mode is given by 2911.38 cm -1 in the present study whereas this stretch vibration is given by 2866.90 cm -1 in (Selvam, Chen, & Wang, 2010).As a result, this model is applied to calculate the vibrational frequencies.A scale of 0.9613 is applied.Lorentzian broadening shape function with a full width at half maximum (FWHM) of 10 cm -1 and 4 cm -1 , respectively, are applied to the main and insert spectra

Solvent Effects on the IR Spectra
Explicit solvent-solute interaction is the most appropriate way to study the solvent-solute interactions and to track the solute property changes upon introduction of solvation (Ganesan et al., 2012).For studying the first hydration shell, for example, explicit water molecules must be included as we will see in the next section (3'-dC solvation shell).The computational cost of the explicit solvent models, however, is often very high.Alternatively, implicit solvent models are a cost effective way to study the solvent effects.This model studies the effects of solvent on molecular properties rather than the explicit solvent-solute interactions.Several solvent models, such as the generalized-Born model (GB model), the Poison-Boltzmann model (PB model) (Bashford & Case, 2000), the conductor like Screening MOdel (COSMO) (Klamt & Schuurmann, 1993) and the polarisable continuum model (PCM) have been previously described in several studies (Roux & Simonson, 1999;Feig & Brooks-Iii, 2004;Mennucci et al., 2002;Tomasi, Mennucci, & Cammi, 2005).Like most implicit solvent models, the PCM model ignores direct short range solvent-solute explicit interactions, which may be a limitation if the explicit interaction is necessary, namely for solute-solvent complexes exhibiting hydrogen bonds (Alagona & Ghio, 2006).The red-shifts or blue shifts showed in the IR spectra of I and IIA in solvents could provide useful information to reveal their intramolecular interactions as found by (Selvam, Chen, & Wang, 2010).That is, the nucleosides in each solvent are optimized again before their IR spectral calculations.Four solvents of various polarities are selected in the study including benzene ( = 2.27), tetrahydrofuran (THF,  = 7.43), tetrahydrothiophene-s,s-dioxide (THT,  = 43.96)and water ( = 78.36).Table 3 presents some selected IR spectral shifts of I and IIA in those solvents with respect to the gas phase spectra in the IR region of  > 1500 cm -1 .Solvent causes mostly red-shift I and IIA and the amplitudes of the shifts are proportional to the dielectric constant of the solvents.For I, small red-shift (|| < 20 cm -1 ) occurs in the region of > 3460 cm -1 .In the region of  < 1700 cm -1 , the red-shift of the C=O vibration is the largest in all solvents, .The C(3)=C(4) double bond exhibits a shift switching from red shift in non(low)-polar solvents to blue shift in polar solvents: the C=C stretch vibration slightly red shifts (< 0) in non-polar solvents such as benzene but switches sign in polar solvent such as water (> 0).The -NH 2 asymmetric band at 3590.44 cm -1 red shifts -4 cm -1 in benzene ( = 2.27), -10.29 cm -1 in THF ( = 7.43), -13.91 cm -1 in THT ( = 43.96)and -14.32 cm -1 in water ( = 78.36).The largest red-shift is the C=O stretch vibration at 1692.86 cm -1 , which red shifts -36.37 cm -1 , -68.51 cm -1 , -82.21 cm -1 , and -83.28 cm -1 in benzene, THF, THT and water solvents, respectively.
As shown in Table 3, the C=O stretch vibration of IIA exhibits the largest red-shift in solvent and the shift increases as the dielectric constant becomes larger, the shift is given by -27 cm -1 , -60 cm -1 , -66.43 cm -1 and -67.00 cm -1 in benzene, THF, THT and water solvents, respectively, with respect to their gas phase vibration.

3'-Dc Solvation Shell
Figure 5(a1) represents the final trajectory snapshot of 3'-dC together with the nearest 10 water molecules.We have also included a confirmation that is similar to the gas phase/implicit solvent globally optimized conformation in Figure 5(a2).Figure 5b displays the radial pair distribution functions (RDFs) between this atom (O 2 ) and any water hydrogen (HW) within the first hydration shell.The figure also includes the RDF data between (O 5' )H 5' and any water oxygen (OW) within the first hydration shell.From a geometrical point of view, it could have been very convenient to turn on a pure QM Hamiltonian, but, this is very challenging for several reasons.First; the computational expense associated with such treatment is very high limiting the simulation to the ps time scale.Second; for such simulations at this short time scale, it's very unlikely to effectively sample the conformational space and more complex sampling techniques such as Replica Exchange Molecular Dynamics (REMD) may also be necessary (Rosta & Hummer, 2009;García & Sanbonmatsu, 2002;Sugita & Okamoto, 1999).Research in this direction is currently in progress.
As shown in Figure 5(a,b), the O 2 atom is tetra-coordinated directly with four water molecules by the aid of four intramolecular H-bonds.This is also explained by the very strong peak at 1.45 Å (attributed to O 2 HW 1 and O 2 HW 2 ) as shown in the RDF figure (Figure 5b).For H 5' OW, the RDF plot shows a much lower (almost 1/4) peak at around 1.75 Å indicating mono-coordination.

Conclusions
The present study reveals how the properties and vibrational spectra of the cytidine nucleoside derivative, 3',4'-didehydro-3'-deoxycytidine (I) from those of 3'-deoxycytidine (II), as the saturation of the C(3')=C(4') bond in the sugar ring.The study also identified the presence of two low energy conformers and linked to the two experimentally resolved conformations, IIA and IIB.Such structural alternation impacts on their intramolecular H-bonding network, such as O(2')-HO(2)=C(2), which has been indicated by their IR spectral shifts.It is found that the double bond leads to hydrogen bond related IR spectral red-shifted in solvents.The polarisable continuum model (PCM) confirms our previous founding that the amplitudes of the red-shifts increase as the solvent dielectric constants (Selvam, Chen, & Wang, 2010).Studying the first hydration shell reveals that water molecules have higher preference for the carbonyl oxygen which is tetra-coordinated with four water molecules.The present study further indicates using the Fukui function that the double bond can be used to activate or deactivate sugar carbon sites by chemical modification, which provides useful information for drug design.
Figure 1.(a) The 2D and 3D structures with the Hirshfeld charges in parentheses of 3',4'-didehydro-3'-deoxycytidine (3',4'-D3C, i.e., I) and the two conformers of 3'-deoxycytidine (IIA (gg) and IIB (gt)) and (b) the two 3'-dC conformers (IIA and IIB)showing the calculated total energy difference also gives the atomic Hirshfeld charges on the atoms other than hydrogens calculated using the LB94/et-pVQZ model.In general, for the non-hydrogen atoms, nitrogens and oxygens are negatively charged whereas carbons are either negatively (when not directly bond with N or O) or positively charges (when directly bond with N or O), depending on their positions in the nucleoside.Saturation of the C(3')=C(4') in I doesn't affect the charge distribution of all sites of the nucleoside evenly: the apparent Q H changes happen almost locally in the vicinity of C(3')=C(4') of sugar except for O(2) of the keto oxygen.For example, apparent changes Q H are on the C(3'), C(4'), C(5'), O(2), O(2') and O(4') sites when the C(3')=C(4') become C(3')-C(4').

Figure 2 .
Figure 2. Comparison of the atomic site based condensed Fukui function of I, IIA and IIB

Figure 3 .
Figure 3.Comparison of thesimulated IR spectra of I (lower panel) and IIA (upper panel) in vacuum.A scale of 0.9613 is applied.Lorentzian broadening shape function with a full width at half maximum (FWHM) of 10 cm -1and 4 cm -1 , respectively, are applied to the main and insert spectra

Figure 4 Figure 4 .
Figure4(a,b) illustrates the significant IR spectral shift associate with the C(2)=O(2), O(2')=H and O(5')-H in I and IIA, respectively.Significant red shift of the C=O vibration is seen in both nucleosides, which is enhanced by polar solvents.The trends somehow obeys an exponential relation as discussed earlier by(Selvam, Chen, & (a1) First hydration shell of the final QM/MM MD trajectory snapshot of 3'-dC, (a2) first hydration shell of an intermediate QM/MM MD trajectory snapshot of 3'-dC that is similar to the optimized gas phase/implicit solvent conformation, potential solute-solvent H-bonds are shown in dashed lines and (b) RDF plots for between two 3'-dC atomic sites and the corresponding solvent atoms; the O 2 HW and the O(5')HOW

Table 2 .
Comparison of selected simulated IR frequencies and intensities of I and IIA (in vacuum)

Table 3 .
IR spectral frequency shifts I and IIA in various solvents with respect to vacuum a a Based on the B3LYP/6-311++G** level of theory and applying the PCM model for solvent effect, scaled by 0.9613, ∆v values are coloured according to the type of shift, red for red shifts and blue for blue shifts.