Structure and Stability of Chemically Modified DNA Bases: Quantum Chemical Calculations on 16 Isomers of Diphosphocytosine

We studied for the first time 16 tautomers/rotamers of diphosphocytosine by four computational methods. Some of these tautomers/rotamers are isoenergetic although they have different structures. High-level electron correlation MP2 and MP4(SDQ) ab initio methods and density functional methods employing a B3LYP and the new M06-2X functional were used to study the structure and relative stability of 16 tautomers/rotamers of diphosphocytosine. The dienol tautomers of diphosphocytosine are shown to be much more stable than the keto-enol and diketo forms. The tautomers/rotamers stability could be ranked as PC3=PC12<PC2=PC11<PC1<PC10<PC8<PC9<PC15<PC16<PC6∼PC7<PC13<PC4∼PC14<PC5. This stability order was discussed in the light of stereo and electronic factors. Solvation effect has been modeled in a high dielectric solvent, water using the polarized continuum model (PCM). Consideration of the solvent causes some reordering of the relative stability of diphosphocytosine tautomers: PC3∼PC12∼PC2∼PC11<PC1<PC10<PC8<PC9<PC15∼PC16<PC13<PC6∼ PC7∼PC14<PC4∼PC5.


Computational Methods
All quantum chemical calculations were performed using the Gaussian 09 suite of programs [69].Geometry optimizations for all tautomers have been performed using ab initio density functional theory (DFT) at the B3LYP functional [70][71][72] and the new M06-2X functional in conjunction with the 6-31G(d) basis set.For each stationary point, we carried out a force constant harmonic frequency calculation at the same levels to characterize their nature as minima or transition states and to correct energies for zero-point energy and thermal contribution.The frequencies are scaled by a factor of 0.98.The vibrational modes were animated using the ChemCraft program [73].Natural charges were computed within full Natural Bond Orbital Analysis, using NBO implemented in Gaussian 09 [69].The vibrational modes were examined by using the ChemCraft program [73].Partial charge distributions were calculated using the natural population analysis (NPA) method [74].Solvation effect has been modeled in water using the polarized continuum model (PCM) of Tomasi et al. [75][76][77][78] at the B3LYP/6-31+G(d,p) level for all tautomers/rotamers.

Geometrical Aspects of Diphosphocytosine Tautomeric
Forms.Geometry of the 16 tautomers and rotamers was optimized at four different levels.Optimized structures of the 16 diphosphocytosine tautomers (at MP4(SDQ)/6-31+G(d)) are presented in Figure 1.For the notations of the tautomers and rotamers of diphosphocytosine, we use PC1, PC2, PC3, . . .16 for the sixteen forms.All structures with their notations and atomic numbering are shown in Figure 1.Structure geometrical parameters of all the 16 tautomers are listed in Table 1.Bond length changes among tautomers show the same trend as expected from formal bond orders.In tautomers PC2, PC3, PC11, and PC12, all CC and CP bond lengths are in the range characteristic for aromatic rings.C1-C2 bond length changes among tautomers depending on the amount of -bond character the tautomer has.In tautomers PC2, PC3, PC11, and PC12 this bond distance is about 1.40 Å while in the rest of the tautomers this distance increased as shown in Table 1.In all the dienol forms, the P3-C2 and P3-C4 and also P5-C6 and P5-C4 distances are nearly similar while those for the other keto-enol and diketo forms show considerable differences of about 0.147 Å.This reflects the extent of conjugation and the increasing aromatic character of the six-membered ring.An inspection of the dihedral angles in Table 1 proves that all the tautomers/rotamers structures are nonplanar.It is reported that [79] results of electron correlation calculations (MP2/6-31G * ) give nonplanar geometries for several nucleotide bases specifically, for cytosine, which showed a significant pyramidalization of the amino group from MP2/6-31G * calculations.Using enormous computer resources, a complete MP2/6-31G * * frequency calculation was carried out [80] for cytosine, confirming that the planar structure is a saddle point.Bond angles in the ring show quite dramatic changes between tautomers and rotamers, by up to more than 10 ∘ .Owing to the lack of polarization functions on hydrogens, geometry parameters of the latter are expected to show larger differences, but these less important parameters were not listed in.
In the dienol form of diphosphocytosine, there are three double bonds giving it the six -electrons needed to gain its aromatic character as in case of benzene.The phosphorous atom has an additional lone pair of electrons but not involved in the conjugated system.Once the keto form is formed in the diphosphocytosine, one resonance form is lost as the phosphorus atom is strongly pyramidal in the keto forms and its lone pairs are no longer available for conjugation.The extent of pyramidalization emerges from the sum of bond angles around the phosphorus atom.An extra destabilization occurs upon forming the second keto tautomer.Accordingly, diketo form is accompanied by destruction of conjugated electrons system in the ring of the diphosphocytosine.In case of cytosine, on the other hand, formation of keto or diketo forms still keeps planarity at the nitrogen atoms where they still able to resonate with the adjacent double bonds.Aromaticity is lost in both cytosine and diphosphocytosine, but the number of resonance structures decreased in the case of the latter system.

Relative Energies and Thermodynamic Parameters of
Tautomers.Since our major objective of this study was to obtain information on the relative stabilities of diphosphocytosine tautomers and rotamers, we tried to go as high as possible, with the level of theory and use both DFT and high correlation ab initio MP2 and MP4(SDQ) methods.The results on the relative stabilities of diphosphocytosine tautomers and rotamers using different computational levels are listed in Table 2.For a better view, this is also shown in Figure 2. The results of the four computational methods are consistent with slight differences.The 16 diphosphocytosine tautomers and rotamers could be classified, on the basis of total relative energy, into three groups: (1) low-energy tautomers group containing four tautomers namely, PC2, PC3, PC11, and PC12.These tautomers: are those satisfying the aromaticity of the pyrimidine ring (the benzene-like structure), (2) moderate-energy tautomers group containing also four tautomers, PC1, PC8, PC9, and PC10.(3) high energy tautomers group containing the rest eight tautomers, PC4, PC5, PC6, PC7, PC13, PC14, PC15, and PC16.DFT methods predict PC2 to be the most stable lowest-energy tautomer.High correlation ab initio methods, on the other hand, predict the most stable is the hydroxy-amino form PC3. Energy difference is less than 0.50 kcal/mol.The structural difference between the two tautomers/rotamers is the proximity of the positively charged H to the positively charged phosphorus atom.The same argument is applied to PC11 and PC12, where the positively charged H is repulsively interacting with P3 (0.503, Table 3) in PC11 and P5 (0.531, Table 3) in PC12.Previous studies suggest that the canonical amino-oxo structure is the most stable one [81].PC5 is the most unstable highest-energy tautomer as predicted by  all methods.This could be due to the destruction of the aromaticity and also attributed to the repulsive interaction of the positively charged H and the great positive charge accommodated on P3 (0.657, Table 3).As could be seen from Table 2 that including zero-point energy and thermal energy correction has a considerable effect on the molecular energy.The corrected energy is lower than the uncorrected by about 4.0 kcal/mol.Also, considering these corrections causes reordering of the tautomers stability.Figure 3 shows the MP2-corrected for zero-point and thermal energy versus uncorrected one for all the 16 diphosphocytosine tautomers.Based on the corrected energy, both B3LYP and M06-2X have  the same order as On the other hand, the relative stability order as predicted by ab initio Moller Plesset MP2 and MP4(SDQ However, MP4(SDQ) method shows slight differences like PC15 ∼ PC16 and PC4 < PC14 ∼ PC5. Figure 4 presents tautomers' order by DFT and MP2 methods.Enthalpy and free energy are consistent with the total energy as can be detected from Table 2.One of the electronic factors that affect the stability of the tautomers is the repulsion between the neighboring -OH and -PH or -CH groups, while the   other factor is the attraction between the O-H or P-H bonds and the lone pair of the neighboring phosphorous or oxygen atoms.Another factor is the attractive interaction between the O-H or P-H bonds, on one side, and the lone pair(s) of the adjacent P or O atoms, on the other one [82].
The amino N accommodates the highest negative charge (up to −0.904); then the O atom comes in the second place (up to −0.778).Positive charges are accommodated on the P and H atoms. Very little negative charge on H bonded to P atoms as appeared in PC4, PC6, and PC8.The values of the NPA charges at hydrogen atoms of -OH groups are larger than the values at the hydrogen atom of the -CH group (see Table 3).As a result, the repulsion between the neighboring -OH and -CH groups is stronger than the repulsion of the neighboring -PH and -CH groups.

Solvation Study.
One of the most crucial factors determining the tautomer distribution in the biological media is the environment.Solvent preferentially stabilizes structures with localized charges.Solvation effect has been explored in water ( = 78.54)using the polarized continuum model (PCM) at the B3LYP/6-31+G(d) level for all the 16 tautomers.Table 4 collects the relative energies for different tautomers/rotamers (kcal/mol) at B3LYP/6-31+G(d) in water.In water, the most stable diphosphocytosine tautomers are PC3 and PC12 and the most unstable tautomers/rotamers are PC4 and PC5.As seen from Table 4, the relative stabilities in water go in the order This order is slightly different from that predicted at the same level in gas phase ( = 1).For example, PC3 ∼ PC12 ∼ PC2 ∼ PC11.This result is easily understandable on the basis of the calculated dipole moment of these tautomers.Polar tautomers/rotamers are favored in polar solvents.
The NPA charges computed for tautomers/rotamers geometry optimized in water are given in Table 5.
The NPA charges computed for tautomers/rotamers geometry optimized in water are generally smaller than those computed for tautomers/rotamers optimized in gas phase.This is understandable as a role of the high dielectric constant solvent, water.The amino N accommodates the highest negative charge (up to −0.846), then the O atom comes in the second place (up to −0.707).Positive charges are accommodated on the P and H atoms. Very little negative charge on H bonded to P atoms as appeared in PC4, PC6, and PC8.The values of the NPA charges at hydrogen atoms of -OH groups are larger than the values at the hydrogen atom of the -CH group (see Table 5).As a result, the repulsion between the neighboring -OH and -CH groups is stronger than the repulsion of the neighboring -PH and -CH groups.

Summary
On the basis of the present theoretical study, we may confirm the following conclusions.
(i) We studied for the first time 16 tautomers/rotamers of diphosphocytosine by four computational methods.Some of these tautomers/rotamers are isoenergetic although they have different structure.
(ii) From the total energy point of view, the 16 tautomers/rotamers are categorized into three main groups.These are low-energy difference group, moderate-energy difference group, and high-energy difference group.
(iv) Our systematic calculations show that the energy difference between components of a rotamer pair is very stable, independent of the level of theory.This leads to the important conclusion that besides the dominant PC3 form, its pair PC2 may lie not higher than ∼0.4 kcal mol −1 ; therefore, any experimental study should consider this possibility when evaluating spectroscopic experimental studies.

ISRN Physical Chemistry
(vi) Charge is distributed on the ring atoms; −ve charges are on N and O, while +ve charges on P and H.

Figure 2 :
Figure 2: Relative energy of different diphosphocytosine tautomers at different levels.

Figure 3 :
Figure 3: MP2 corrected for zero-point and thermal energy versus uncorrected one for all the 16 diphosphocytosine tautomers.

Table 2 :
Thermodynamic parameters (kcal/mol) a for diphosphocytosine tautomers/rotamers at different levels in gas phase.
a Corrected for zero-point energy and thermal energy; b uncorrected energy.
a Uncorrected energy; b corrected for zero-point and thermal energy.