Conformational Analysis, Molecular Structure and Solid State Simulation of the Antiviral Drug Acyclovir (Zovirax) Using Density Functional Theory Methods

The five tautomers of the drug acyclovir (ACV) were determined and optimised at the MP2 and B3LYP quantum chemical levels of theory. The stability of the tautomers was correlated with different parameters. On the most stable tautomer N1 was carried out a comprehensive conformational analysis, and the whole conformational parameters (R, β, Φ, φ1, φ2, φ3, φ4, φ5) were studied as well as the NBO Natural atomic charges. The calculations were carried out with full relaxation of all geometrical parameters. The search located at least 78 stable structures within 8.5 kcal/mol electronic energy range of the global minimum, and classified in two groups according to the positive or negative value of the torsional angle φ1. In the nitrogen atoms and in the O2' and O5' oxygen atoms of the most stable conformer appear a higher reactivity than in the natural nucleoside deoxyguanosine. The solid state was simulated through a dimer and tetramer forms and the structural parameters were compared with the X-ray crystal data available. Several general conclusions were emphasized.


Introduction
Acyclovir (ACV, 9-(2-hydroxyethoxymethyl) guanine, Zovirax [1,2], Figure 1), is a guanine derivative possessing antiviral activity and commonly used in the treatment of herpes. It is a potent antiviral agent that is used as a highly specific inhibitor of herpes viruses (HSV) types 1 and 2 [3][4][5][6][7]. A series of new guanine base modified tricyclic analogues of ACV and ganciclovir were evaluated for activity against herpes simplex virus type 1 and 2, showing similar antiherpetic potency as the parent compounds ACV and ganciclovir [8]. The antiherpetic activity was found to be strongly dependent on the nature and esteric demands of the substituents in the 6 and/or 7 positions [9]. Most persons who are infected with human immunodeficiency virus type (HIV-1) are also infected with herpes simplex type2 (HSV-2), which is associated with increased plasma and genital levels of HIV-1. Thus, it has been reported that ACV inhibits HIV upon human herpes virus (HHV) coinfection in tissue cultures. This activity was found to be correlated with the phosphorylation of the parent drug to the monophosphate form mediated by HHV-encoded kinase [10]. Recent studies show that ACV decreases the HIV-1 RN and suppress both viruses in coinfected tissues [11], with great impact on HIV and HSV-2 seropositive patients [12,13]. The treatment of ACV-resistant HSV infections [7,14] and tests for release it at oral therapeutic levels have been also reported [15].
As it is shown above, ACV has been extensively studied from the pharmaceutical and medical point of view. However, few studies appear on its molecular structure, only at low level [16], and to our knowledge, there is no data on its conformational characteristics by Density Functional Methods (DFT) or ab initio quantum chemical methods. Thus, it is one of the tasks of the present manuscript. Conformations of some derivatives of ACV with biological activities have been reported [17], as well as studies of complexes of ACV with several metals [18,19].
The tautomeric study of DNA structural components has a great interest today with numerous research publications [20,21]. The tautomeric equilibrium in ACV between the keto and enol forms has been observed from the UV/VIS spectra [22]. This equilibrium depends on the polarity of the solvent, and therefore, in water solution the keto form prevails, while in methylene chlorides it is the enol one. Thus, another goal of the present manuscript is to study the possible tautomers of ACV, and the determination of their % populations at room temperature and at 273.15 K. The stability of these tautomers and its dependence of different parameters is another point analyzed here.
Finally, from our understanding, would be interesting to calculate the different conformational possibilities of ACV and compare the results with the natural nucleoside deoxyguanosine (dG). This is the last goal of the present manuscript. An accurate knowledge of the different conformers of ACV, its charge distribution, inter-and intra-molecular interactions, solid structure, and flexibility would be an important help for the interpretation of drug-target interactions, as well as to design new antivirus.
For this reason, the conformers of natural and analogues nucleosides have been analyzed by different authors. Now, an extensive theoretical study of the conformational preferences in ACV has been carried out with full relaxation of all geometric parameters, in an attempt to gain insights into molecular features responsible for activity. We will attempt to determine herein, if the various geometric features in ACV are correlated or interact with one another. We are also interested in whether alternative forms of hydrogen bonding make significant contributions to the conformational behavior of ACV.
conformers to confirm that they corresponded to local minima. All optimized structures showed only positive harmonic vibrations (local energy minima). Relative energies were obtained by adding zero-point vibrational energies (ZPEs) to the total energy. For the calculation of the ZPEs, the frequencies were retained unscaled. The ΔG values were sums of electronic and thermal Free Energies. The conformational equilibrium at 298.15 K was evaluated by means of the Boltzmann distribution formula exp(−ΔG/kT), where ΔG is the relative Gibbs energy.

Results and Discussion
ACV has five possible tautomers (Figure 2), that were fully optimized at different levels of computation, Table 1. The most stable one corresponds to N1, and thus we have focused the study only in this tautomer. The remaining forms were left for future research. In the last two column of this Table is shown the % population of the different tautomers at 298.15 K and 273.15 K. At room temperature the largest population corresponds to tautomer N1 (48.1%). The second population is due to OHC tautomer (37.7%) and the third one is to OHT (14.3%). Tautomers N3 and N7 have very little population, less than 0.05%.

N1
N3 N7 OHC OHT An analysis of the relative energies of these tautomers shows that they can be related to the dipole moment (μ) and to the torsional angle φ 1 , Figure 3. Thus, the least stable tautomer in the isolated state (N7) has the highest μ, i.e., it is the most stable in water solution. Other relations can be observed between ΔE and φ 1 . It is noted that tautomer N7 has the highest negative value of φ 1 angle.
Considering the structure of the chain, that likes the structure of the sugar in the nucleosides, and in accordance to previous works [45], other three structural parameters were defined to fix the chain position respect to the plane of the nucleobase: (vi) The vector R (N9···O5') which determines the distance of the OH group relative to the base; (vii) The angle β (C4-N9···O5') which defines the angle of the OH group relative to the base plane; and (viii) the angle Φ (C1'-N9···O5') which also determines the position of the OH group.

Conformers and Energetics
An extensive conformational study of tautomer N1 was carried out through a rotation of the exocyclic φ 1 , φ 2 , φ 3 , φ 4 and φ 5 torsional angles. A detailed collection of the most important conformational parameters of these optimized forms is included in Table 2. The conformers were classified according to the two ranges of rotation of φ 1 : conformers A with the φ 1 values negative and conformers B with φ 1 positive.
Two energy criteria were considered for each conformer: the electronic energy ΔE + ZPE correction, and the Gibbs energy ΔG. For the numbering of the conformers in each range of rotation of φ 1 was followed the ΔE + ZPE criterion. Calculations at different levels, as well as single point calculations at the MP2/6-31G(d,p)//B3LYP/6-31G(d,p) level were carried out to confirm the stability of the main conformers, Table 2 and Table 3. In general, the stability order remains, although several changes are observed. Thus, conformer B2 appears as the most stable one instead of A1 predicted by the DFT methods. The % population of the different conformers at 298.15 K (P 298. 15 ) and at 273.15 K (P 273.15 ) was calculated. It indicates that only confomers A1 (41.5%), B1 (36.1%), B2 (9.6%) and A2 (8.7%) have importance. The population of the remaining conformers is lower than 1% and thus they are not of interest. The temperature effect is not significant.
By the same methodology, the relative energies of the different conformers have been determined in related nucleosides [34,[46][47][48]. The global minimum calculated in these molecules by MP2/6-31G(d,p) was in accordance to that found by B3LYP/6-31G(d,p). Thus, our results by B3LYP can be considered acceptable.           The conformers differ in general very little in energy. Thus, in our calculations 78 optimized conformers were found within the electronic energy range ΔE = 0-8.5 kcal/mol, and Gibbs energy range ΔG = 0-6.6 kcal/mol with respect to the global minimum. This range of values of ΔG is smaller than that calculated in dG, in dT [50], 0-7 kcal/mol and in dU [31], 0-9 kcal/mol.
Only four conformers are found within the electronic energy range ΔE = 0-1.0 kcal/mol (by criterium of ΔE + ZPE), Table 2, with φ 1 −72°/−97° as g − and 70°/96° as g + by MP2. Among these conformers, B1 has the highest dipole moment 6.20 D although very close to that in A1, 6.18 D. These conformers are slightly favored in a polarizable environment with water.
Another seven conformers appear within the electronic energy range ΔE = 1.0-3.0 kcal/mol, three values were anti (φ 1 ca. −100°) and four syn (φ 1 ca. 90°). The anti structures are the expected forms in the natural nucleosides that form the nucleotides and polynucleotides in biological systems [33,35,[46][47][48]. The ratio anti/syn in conformers A is 0.8 in the low-energy group (<3 kcal/mol) but it increases up to 1.7 in the 3.0-8.5 kcal/mol range.
The global minimum by DFT methods corresponds to the conformer denoted as A1 ( Figure 4) and it appears stabilized by an intramolecular H-bond. The optimised bond lengths and natural NBO atomic charges on this conformer are collected in Figure 5. This global minimum by criterium of ΔE + ZPE agree well to that obtained by criterion of ΔG, but differs of that obtained by MP2/6-31G(d) and MP2/6-31G(d,p)//B3LYP/6-31G(d,p), Table 3 and Figure S3. It is because of the small difference in energy between both forms A1 and B2. This global minimum in the syn form by MP2 is in accordance to that obtained in other nucleosides [33][34][35] but differs of the anti form expected for the natural nucleosides and nucleotides in biological systems [43]. The second most stable conformer is B1 with values of φ 1 = 70°, φ 2 = −152° and φ 3 = 94° by MP2. Figure 4 shows the six best optimum conformers selected in the two ranges of φ 1 : three are A (A1-A3), and three are B (B1-B3). The values of the intramolecular H-bonds and the most important structural angles of each conformer are also included. Figures 6-9 show the distribution of the 78 optimised conformers according to their energies, exocyclic torsional angles ϕ 1 -ϕ 5 , the angles β, Φ, and the vector R. The best significant conformers are pointed in these figures.

Conformational Angle Analysis
An overall examination of the five exocyclic torsional angles and two bond angles, defining the conformational space in tautomer N1 of ACV, leads to conclude the following, nucleosides [51][52][53], and it is the form for biological activity [54]. Although, in many nucleoside analogues syn and anti forms have similar energy, however, the global minimum corresponds in general to the syn form, as in AZT [35]. In purine nucleosides there is relatively little restraint to rotation about the N9-C1' bond [55]. In ACV the high population of the anti forms and the wider range of φ 1 values are factors that facilitate its antiviral activity. The sugar ring of the natural nucleoside dG can be mimics by the similar exocyclic angles of ACV. The bond Φ angle also shows a regular distribution in the 10° ≤ Φ ≤ 83° range for conformers A, and 29° ≤ Φ ≤ 95° for conformers B. The value of this angle appears between 51° and 54° for the four best conformers with ΔE < 1 kcal/mol −1 . In general, the most stable conformers do not have large values of Φ, neither low values.
Finally, the R vector has a large range of values, between 2.905Å and 5.956 Å, and it determines the distance of the OH end of the chain from the base plane.
Comparing the bond lengths of ACV (conformer A1) with those of dG, it is observed that ACV has several bonds of the guanine moiety larger than the corresponding bonds in dG, while in the remaining bonds are smaller. The largest difference appears in C4-N3.

The Side Chain
The most important characteristic of the structure of ACV is the conformation of the side chain that is attached to N9, which is characterized by the exocyclic torsional angles φ 1 to φ 5 , by the bond angles β and Φ that determine the position of the chain from the base, and by the distance R. The bonds along this chain can be in trans respect to each other and thus it gives an extended zig-zag structure or, by contrast, one or several alternate bonds can be in gauche in such a way that it resembles at least partially a portion of the furanose moiety of dG [49,56]. In this last case, it is capable of adopting conformations resembling a portion of the pentose ring, a factor which undoubtedly plays an important role in their biological activities [57]. Thus, three types of structures appear: two that have the side chain partially folded (conformers A and B), and one structure which has all bonds in trans orientation and leading to an almost planar zig-zag arrangement [49].

Intramolecular H-Bonds
Several authors have studied the intramolecular H-bonds in related nucleosides, in special using AIM method [58,59]. However, the absence of the furanose ring in ACV reduces largely the number of possible H-bonds. Therefore, only two intramolecular H-bonds may be observed in the main conformers of ACV: (i) The hydroxyl hydrogen H5'(O5') and the guanine's position 3 nitrogen atom, H5'···N3, and (ii) The hydroxyl hydrogen H5'(O5') and the guanine's position 2 oxygen atom, H5'···O2'. Figure 4 shows these types of H-bonds in the six most stable conformers. H-bond (i) appears in conformers A1, B1, B2 and A2; while (ii) is observed in conformers B3 and A3 with values of 2.346 Å and 2.350 Å, respectively. H-bond (i) is stronger than (ii), and it gives a great stability to the structure. Thus, the conformers with the H-bond (i) are the most stable, ΔE < 1 kcal·mol −1 , while the conformers with the H-bond (ii) have higher energy, ΔE > 2 kcal·mol −1 .

Natural NBO Atomic Charges
The calculated values in conformer A1 appear collected in Figure 5. The largest negative charge corresponds to N2 and O5' atoms, −0.90 and −0.84 e, respectively by MP2 (where e is the charge of an electron). The next atom with large negative charge is N1 (ca. −0.7 e). The value of this charge is slightly higher than in dG, −0.69 e. The main effect of the NH 2 group is a remarkable increment in the positive charge on C2, 0.81 e, because of the high negative charge on the amino nitrogen, −0.90 e. Consequently, a noticeable increase of the negative charge on N3 is observed, −0.73 e (−0.70 e in dG). The electron-rich sites of the guanine moiety are involved in H-bonds, with N3 and O2' acting as single acceptors.
In the nitrogen N9 the negative charge (−0.52 e) is lower than on N1, N3 and N7, but slightly higher than in dG. It is because in ACV the bonding to the chain increases the negative charge on N9.
The value of the charge on N7 is important because in the neutral form of anti-tumour platinum drugs, the platinum atom has a strong preference for nitrogen N7 rather than for oxygen atoms of the base for its coordination [51]. Also, the N7 position in DNA is the most open to attack. N1 position is also important because when deprotonation of the weakly acid ACV occurs, the metal binding site changes to N1, which is the formally deprotonated site [60]. Raman spectra of related nucleosides in H 2 O indicate that the site of deprotonation in basic solutions is N1, while the site of protonation in acidic solutions is N3 and N7, the same sites that in its phosphorylated form. These results are useful for identification and characterization of its structure in natural occurring biopolymers [61].
The negative charge on the oxygen atoms is high, −0.65 e on O2' and −0.69 e on O6. The values in O2' and O5' are slightly higher than in dG, by B3LYP −0.603 e and −0.761 e, respectively. By contrast, in O6 the negative charge is lower than in dG, −0.592 e.
C6 is the atom with the highest positive charge, 0.829 e by MP2, in concordance to the high negative charge on O6, i.e., it is the most reactive. With a slightly lower value appears C2 (0.812 e), and with much lower values H5'(O5 ' ), H1(N1), H2(N2) and H2'(N2). The remaining hydrogen atoms have much less positive charge and they are less reactive.

Solid State Simulation
The structure of ACV (CSD code CEHTAK) in the solid state was determined by X-ray diffraction [53,[60][61][62][63]. Birnbaum et al. [49,62] obtained the crystals that belong to monoclinic space with three independent molecules (A, B and C), together with two water molecules [49]. Molecules A and B have similar conformation while that of molecule C is different, Table 4. In molecules A and B the angle ϕ 1 and ϕ 2 are in the preferred g − form, but in molecule C ϕ 2 is in trans. It means that molecules A and B showed a partially folded conformation of the side chain, while molecule C appears extended [62]. Molecules A and B correspond to conformer A5 in the isolated state, while molecule C is represented by conformer A34. Molecules A and B are ca. 3 kcal·mol −1 more stable than molecule C. The difference in energy between molecule A and B is small and it is due to the flexibility of the side chains. Table 4. A comparison of the most important structural parameters in the dimer and tetramer forms calculated at the B3LYP/6-31G(d,p) level in acyclovir molecule with those in the crystal. The torsional angles are in degrees and R in Å. We have simulated this arrangement through a dimer and tetramer forms, Figure 10. The calculated value of the inter-and intramolecular H-bonds distances and the total energy is included in the Figure. The calculated bond lengths, bond angles and torsional angles are collected in Tables S1 and S2. The simulated structure is tightly bonded but not so tight as it the reported in the crystal. There are four donor sites: H1(N1), H2(N2), H2'(N2) and H5', and five possible acceptors: O6, N3, N7, O2' and O5'.  Comparing the calculated bond lengths and angles of molecule A in the dimer with the X-ray data [62], Tables S2 and S3, it is observed the good agreement in the values of our simulated structure. The greater differences appear in C4=C5 (0.025 Å), C2-N2 (0.023 Å), C3'-C4' and N1-C6 (0.016 Å), and the difference is almost null in C6=O, C5-N7, N7=C8 and C4-N9. A slight larger error is obtained when the comparison theory-experiment is carried out with the molecule B of the dimer. This error is slight reduced when the comparison is carried out with the simulated molecule A of the tetramer. In the angles the difference is almost null, with similar values of all the angles. Only several angles in the side chain show some noticeable differences, i.e., C3'-C4'-O5', 4.7°.

Solid form
Other four anhydrous forms of ACV and a new hydrate have been characterized by X-ray diffraction, with significant differences in the intermolecular H-bonding networks among the ACV forms [63][64][65]. A study of the solution forms of ACV shows that ACV can exist as polymorphic and pseudopolymorphic solvates [66]. We have simulated other dimer, trimer and pentamer forms with different conformers and collected in Figure S4.

Conclusions
In the present work we have shown a comprehensive compendium of the possible conformers in tautomer N1 of ACV. The geometries and values of the properties presented here appear to the most accurate to date. The most important findings of the present work are the following: (1) Five tautomers of ACV were identified and fully optimized. At room temperature only tautomer N1 (48.1%), OHC (37.7%) and OHT (14.3%) have a noticeable population. It is very small in tautomers N3 and N7, less than 0.05%. (2) The relative energies of the five tautomers appear related to the dipole moment and to the torsional angle φ 1 . The least stable tautomer N3 in the isolated state has the highest μ and thus, it is the most favoured in a polarisable environment with water. (3) In the isolated state the most stable tautomer is N1 by both B3LYP and MP2 methods. In this tautomer, and through a rotation of φ 1 , φ 2 , φ 3 , φ 4 and φ 5 exocyclic torsional angles, 78 optimized stable conformers were identified, two syn and two anti failing into the 0-1 kcal·mol −1 ΔE + ZPE energy range. (4) The calculated most stable conformer by all DFT levels corresponds to A1, while by MP2 is B2. In the nitrogen atoms and in the O2' and O5' oxygen atoms of conformer A1 appear a higher reactivity than in the corresponding natural nucleoside deoxyguanosine. (5) The distribution of all the conformers according to the ranges of stability of the characteristic torsional angles was established. The values obtained indicate the flexible nature of ACV, which is higher than dG. An increase of the stability appears when the side chain is near to the purine base, with a value of R that fails into 3.925-2.892 Å range, and an angle Φ close to 54°. (6) Only two intramolecular H-bonds may be observed in the main conformers of ACV, in contrast to the six H-bond types calculated in dG. It leads to a flexibility higher in ACV than in dG. (7) The solid state was simulated through a dimer and tetramer forms. An excellent agreement with the X-ray crystal data was obtained, which indicates the good accuracy of the theoretical methods used.