Stability of Nucleic Acid Base Pairs in Organic Solvents: Molecular

Feb 16, 2007 - Structure and dynamics of some macrocyclic pyrimidine derivatives. A. V. Kozlov , V. E. Semenov , A. S. Mikhailov , A. V. Il'yasov , V...
17 downloads 0 Views 1MB Size
J. Phys. Chem. B 2007, 111, 2591-2609

2591

Stability of Nucleic Acid Base Pairs in Organic Solvents: Molecular Dynamics, Molecular Dynamics/Quenching, and Correlated Ab Initio Study Lucie Zendlova´ , Pavel Hobza, and Martin Kabela´ cˇ * Institute of Organic Chemistry and Biochemistry, Academy of Sciences of the Czech Republic, and Center for Biomolecules and Complex Molecular Systems, FlemingoVo na´ m. 2, 166 10 Prague 6, Czech Republic ReceiVed: August 22, 2006; In Final Form: December 21, 2006

The dynamic structure and potential energy surface of adenine‚‚‚thymine and guanine‚‚‚cytosine base pairs and their methylated analogues interacting with a small number (from 1 to 16 molecules) of organic solvents (methanol, dimethylsulfoxide, and chloroform) were investigated by various theoretical approaches starting from simple empirical methods employing the Cornell et al. force field to highly accurate ab initio quantum chemical calculations (MP2 and particularly CCSD(T) methods). After the simple molecular dynamics simulation, the molecular dynamics in combination with quenching technique was also used. The molecular dynamics simulations presented here have confirmed previous experimental and theoretical results from the bulk solvents showing that, whereas in chloroform the base pairs create hydrogen-bonded structures, in methanol, stacked structures are preferred. While methanol (like water) can stabilize the stacked structures of the base pairs by a higher number of hydrogen bonds than is possible in hydrogen-bonded pairs, the chloroform molecule lacks such a property, and the hydrogen-bonded structures are preferred in this solvent. The large volume of the dimethylsulfoxide molecule is an obstacle for the creation of very stable hydrogen-bonded and stacked systems, and a preference for T-shaped structures, especially for complexes of methylated adenine‚‚‚thymine base pairs, was observed. These results provide clear evidence that the preference of either the stacked or the hydrogen-bonded structures of the base pairs in the solvent is not determined only by bulk properties or the solvent polarity but rather by specific interactions of the base pair with a small number of the solvent molecules. These conclusions obtained at the empirical level were verified also by high-level ab initio correlated calculations.

Introduction Noncovalent interactions play an important role in such biological systems as DNA and proteins. The structures of these biopolymers are determined not only by the presence of covalent bonds in the proteins or sugar-phosphate backbone (in DNA) but chiefly by noncovalent interactions, based on the balance between the hydrogen bonds (H-bonds), which are of electrostatic origin, and van der Waals interactions mainly directed by London dispersion as an attractive force.1,2 Specific H-bonds are necessary for correct base pairing in DNA, fidelity of replication, transcription processes, and generally for the formation of the secondary structures of biopolymers. Nonspecific dispersion forces, however, are responsible for basestacking interactions in the DNA, for the interaction of ligand with DNA, as well as for the interaction of aromatic side chains of amino acids with peptide bonds. The role of the environment (solvent and metal ions) cannot be underestimated, however. While the metal cations can neutralize the negatively charged phosphate groups of nucleic acids, water is essential for the stability of the double helix of DNA, because the helical form is unstable in the gas phase and the addition of nonaqueous solvents lowers its melting temperature.3 Typical denaturation agents of DNA, urea, formamide, methanol (CH3OH), and dimethylsulfoxide (DMSO), are used to lower DNA duplex thermostability during electrophoresis or polymerase chain reaction.4 * Author to whom correspondence should be addressed. Fax: (+420) 220 410 320. E-mail: [email protected].

The experimentally known very limited solubility of nucleic acid bases in water5,6 was confirmed theoretically by the calculation of solvation free energies and partition coefficients of nucleic acid bases in chloroform (CHCl3) and water, using both the explicit7-10 and the implicit11,12 model of solvents. Problems with the experimental studies of H-bonding in the presence of water arise from the facts that H-bonded water can absorb in the same IR region as nucleic acid bases and that in the NMR technique the protons of the solvent can rapidly exchange with amino protons of the bases, which leads to the preference for organic solvents in the experiments. It was ascertained experimentally (using IR, NMR, and calorimetric titration techniques) that isolated nucleic acid bases in a vacuum and in nonpolar solvents, such as carbon tetrachloride (CCl4) and CHCl3, associate by hydrogen bonding,13-16 whereas in aqueous solution bases form stacked (S) complexes.17,18 Katz and Penman,19 Shoup et al.,20 and Newmark et al.21 have found by proton magnetic resonance studies that guanosine and cytosine specifically interact in DMSO by H-bonds. These findings have been confirmed by free energy/molecular dynamics empirical potential studies by Kollmann et al., who calculated association free energies of base pairs in a vacuum and water22,23 using the Cornell et al. force field,24 and Monte Carlo simulations of base pairs in CHCl3 by Pranata et al.,25,26 in CCl4 by Pohorille et al.,27,28 and in DMSO by Danilov et al.29,30 Whereas all of the above-mentioned theoretical papers are dedicated to the behavior of nucleic acid bases or base pairs in

10.1021/jp065418j CCC: $37.00 © 2007 American Chemical Society Published on Web 02/16/2007

2592 J. Phys. Chem. B, Vol. 111, No. 10, 2007 the bulk solvent, much less interest has been devoted to the exact description of the solvation pattern around the base pairs, and it focused exclusively on hydration. The hydrated hydrogenbonded (HB) base pairs are almost the sole points of interests of the authors.31-33 The pioneering work describing the hydration of S base pairs was done by Claverie et al.34 Sivanesan et al. observed that for guanine‚‚‚cytosine35,36 and cytosine‚‚‚ cytosine37 base pairs the S base pairs hydrate better than the HB base pairs. The first hydration shell of the S pair can, according to these authors, accommodate 5-6 water molecules while the HB pair only 4-5. Nevertheless, only single-point correlated calculations based on Hartree-Fock and density functional theory optimized geometries performed with small basis sets were used. Our previous molecular dynamics studies devoted to specific microhydration of base pairs38,39 showed that a gradual increase of the hydration number of the adenine‚ ‚‚thymine base pair resulted in a transition from planar HB structures to the S ones. Furthermore, already in the presence of two water molecules, the population of S structures was higher than that of the HB ones. The same conclusions can be applied to the guanine‚‚‚cytosine base pair, but this transition requires a higher number of water molecules.38 These findings were also confirmed by high-level ab initio calculations.40,41 The aim of this study is to examine the dynamic structure and the potential and free energy surfaces of adenine‚‚‚thymine and guanine‚‚‚cytosine base pairs and their methylated analogues exposed to a small number of molecules of organic solvents that are widely used by experimentalists. Various theoretical approaches starting from the simple empirical methods (molecular dynamics, molecular dynamics/quenching technique) to highly accurate ab initio quantum chemical calculations, which are available nowadays for studies of medium-sized biological systems (the MP2 and CCSD(T) methods), will be utilized. Accurate results from the ab initio calculations would also serve as benchmarks to verify calculations at a lower level of theory. Methods Molecular Dynamics Techniques. The molecular dynamics (MD) trajectories of the adenine‚‚‚thymine (AT) and guanine‚‚‚cytosine (GC) base pairs and their N-methylated analogues (mAmT and mGmC) (for schematic drawings of the structures and numbering of the bases see Figures 1aS and 1bS in the Supporting Information) with a small number (1, 2, 4, 8, or 16) of molecules of organic solvents, CH3OH, CHCl3, and DMSO, were determined using the 10 ns NVE (constant number of particles, volume, and energy) MD simulations with the Cornell et al.24 empirical force field. The total kinetic energy corresponded to an average temperature of 298 K. The standard parameters were used for the bases, whereas for the solvents we employed force field parameters and scaled RESP charges calculated at the HF/6-31G* level introduced by Fox and Kollman.42 The scaled charges were applied as they provide better agreement42 with the physical and chemical properties of the solvents. Molecular Dynamics/Quenching (MD/Q) Techniques. The potential energy surfaces (PESs) of nonmethylated and methylated adenine‚‚‚thymine and guanine‚‚‚cytosine base pairs in a vacuum and with one and two molecules of CH3OH, CHCl3, and DMSO were investigated using the MD/Q method with the Cornell et al.24 empirical force field with the same parameters and under the same conditions of the MD part as mentioned in the previous section. The MD run was interrupted every 1 ps, the temperature was changed to 0 K, and the structure of the complex was optimized using the conjugate gradient method.

Zendlova´ et al. TABLE 1: Total Populations (in %) of Structures in a Given Class of Nonmethylated and Methylated Adenine‚‚‚Thymine and Guanine‚‚‚Cytosine Complexes in a Vacuum Obtained by MD Techniques Employing the Cornell et al. Force Fielda orientation of basesb

AT

mAmT

GC

mGmC

p np t s ns separated

36 34 28 1 1 0

11 12 21 7 49 0

48 37 15 0 0 0

44 33 23 0 0 0

aThese data were taken from ref 38. b Classes differ by the mutual orientation of bases and the distances between the centers of masses of bases: p means planar H-bonded (the distance between the centers of masses of bases ranges between 5.0 and 10.0 Å, and the deviation of the base planes lies between 0° and 10°); np means nonplanar H-bonded (the same distance as for the first type, the deviation of base planes lies between 10° and 45°); t means T-shaped (the deviation of base planes lies between 45° and 90°); s means planar stacked (the distance between center of masses of bases is smaller than 5.0 Å, and the deviation of the base planes lies between 0° and 10°); ns means nonplanar stacked (the distance is the same as that for a planar stacked structure, and the deviation of the base planes lies between 10° and 45°).

The minimized structures were saved, and the MD run continued from the point where it had been interrupted. The total length of each simulation was 250 ns. Separation of Structures into Classes Based on the Orientations of Bases. As we have recently shown,38 the S and HB structures can easily be separated by measuring the distances between the centers of the masses of the bases. The common distance for the HB structures is 6 Å, whereas for the S structures it is 3.5 Å. The T structures are characterized by almost perpendicular orientation between the bases. In our study we distinguish five basic structural types: (i) planar HB; (ii) nonplanar HB (nHB); (iii) T; (iv) planar S; (v) nonplanar S (nS); and additionally (vi) structures without any direct contact between the bases (see Table 1 for more detailed information about the separation criteria). All of these geometrical separations were performed using CARNAL, a program for trajectory analysis that is a part of the AMBER program package and our own script already used in our previous studies.40,41 Ab Initio Calculations. The geometries of the global minima of the planar HB, S, and also T structures obtained from the MD/Q method were further reoptimized at the correlated ab initio resolution of identity (RI) MP2 level of theory employing the cc-pVDZ basis set. For complexes with CHCl3, we also added the global minima of the nHB structures, because they had been found to be more stable than the planar ones. To avoid locating transition states, we changed the planar geometry of all of the bases’ amino groups (corresponding to the minimum in the Cornell et al. force field24) to a slightly nonplanar one (corresponding to the minimum at the RI-MP2 level). The stabilization energies were calculated at the RI-MP2 level of theory with various basis sets, cc-pVDZ, TZVPP, aug-cc-pVDZ, and aug-cc-pVTZ, to determine the sensitivity of the results to the basis set used. The last two basis sets were used to compute the stabilization energies at the complete basis set limit at this level of calculations using the two-point Helgaker’s extrapolation.43 All ab initio calculations were performed by the TURBOMOLE 5.6 program package.44 For monosolvated nonmethylated base pairs, we calculated also the CCSD(T) interaction energies with the 6-31G*(0.25) basis set to evaluate any error of the RI-MP2 method. All of the CCSD(T) calculations were carried out using the MOLPRO program.45

Nucleic Acid Base Pairs in Organic Solvents

J. Phys. Chem. B, Vol. 111, No. 10, 2007 2593

Figure 1. Histograms showing the population distribution of nonmethylated and methylated adenine‚‚‚thymine complexes with 0-16 molecules of solvent (CH3OH, DMSO, and CHCl3). Note that non-negligible populations were observed for nonmethylated and methylated AT complexes with 8 and 16 molecules of DMSO (distance between centers of masses of bases larger than 9 Å).

2594 J. Phys. Chem. B, Vol. 111, No. 10, 2007

Zendlova´ et al.

Figure 2. Histograms showing the population distribution of nonmethylated and methylated guanine‚‚‚cytosine complexes with 0-16 molecules of solvent (CH3OH, DMSO, and CHCl3). Note that non-negligible populations were also observed for nonmethylated and methylated GC complexes with 8 and 16 molecules of DMSO (distance between centers of masses of bases larger than 9 Å).

Nucleic Acid Base Pairs in Organic Solvents

J. Phys. Chem. B, Vol. 111, No. 10, 2007 2595

TABLE 2: Total Populations of Structures in a Given Class of Nonmethylated and Methylated Adenine‚‚‚Thymine Complexes in the Presence of 1-16 Molecules of Solvent (CH3OH, DMSO, and CHCl3) Obtained by MD Techniques Employing the Cornell et al. Force Field orientation of basesa/ number of solventb

H2O

CH3OH

CH3OH

DMSO

orientation of basesa/ number of solventb

H2O

CH3OH

CHCl3

DMSO

p np t s ns separated

AT 1 mol 24 6 28 52 35 33 4 2 9 7 0 0

6 59 33 0 1 0

3 35 41 4 17 0

p np t s ns separated

mAmT 1 mol 2 1 3 10 1 17 8 31 77 41 0 0

4 38 40 7 12 0

1 15 37 20 27 0

p np t s ns separated

AT 2 mol 16 3 18 31 1 33 11 10 24 24 0 0

7 64 28 0 0 0

1 15 35 9 40 0

p np t s ns separated

mAmT 2 mol 1 3 1 31 6 33 7 10 85 24 0 0

5 43 40 5 7 0

2 23 44 12 20 0

p np t s ns separated

AT 4 mol 5 1 6 11 17 21 19 23 53 44 0 0

7 65 28 0 0 0

1 17 43 6 33 0

p np t s ns separated

mAmT 4 mol 1 0 1 6 4 10 6 40 88 44 0 0

5 53 37 2 3 0

2 20 38 16 25 0

AT 8 mol p np t s ns separated

1 2 8 21 68 0

0 7 12 30 51 0

8 66 26 0 0 0

1 20 48 6 24 0

p np t s ns separated

mAmT 8 mol 2 0 2 5 5 10 7 39 84 45 0 0

6 55 36 1 2 0

2 22 42 13 20 1

p np t s ns separated

AT 16 mol 2 0 2 2 10 5 20 85 65 7 1 0

4 47 43 1 5 0

1 23 57 3 13 3

p np t s ns separated

mAmT 16 mol 2 0 3 7 7 13 6 35 81 44 1 0

3 33 44 6 14 0

1 21 46 9 16 7

a Classes differ by the mutual orientation of bases and the distances between the centers of masses of bases: p means planar H-bonded (the distance between the centers of masses of bases ranges between 5.0 and 10.0 Å, and the deviation of the base planes lies between 0° and 10°); np means nonplanar H-bonded (the same distance as for the first type, the deviation of base planes lies between 10° and 45°); t means T-shaped (the deviation of base planes lies between 45° and 90°); s means planar stacked (the distance between center of masses of bases is smaller than 5.0 Å, and the deviation of the base planes lies between 0° and 10°); ns means nonplanar stacked (the distance is the same as that for a planar stacked structure, and the deviation of the base planes lies between 10° and 45°). b “1 mol” means 1 molecule of solvent, “2 mol” means 2 molecules of solvent, etc.

Results and Discussion MD of the Base Pairs Microsolvated by Organic Molecules. In Figures 1 and 2 we present histograms showing the distribution of the centers of the masses of the bases and thus the ratio between the HB, T, and S motifs during the simulation in a vacuum and various solvents (see Tables 1-3 for the exact populations of each motif). It was established for all of the systems studied that both the HB and the S structures are nonplanar (about 10-20°) rather than fully planar at room temperature. This deviation from planarity is allowed by easily energetically accessible buckle and propeller modes in the base pairs,46 and this effect also can be observed in the crystal structures of DNA. In Tables 4 and 5 we present the distribution of solvent molecules around the atoms of nucleic acid bases that are either H-bond donors or acceptors. We used the distance criterion to determine whether the solvent molecule binds to selected atom of the base. As a contact we accepted only distances shorter than 3.4 Å between the heavy atoms of bases and the heavy

atoms of solvents (O in the case of CH3OH and DMSO and C in the case of CHCl3). We can conclude that in the presence of one molecule of the solvent there exist preferential binding sites for binding of the solvent. When we increase the number of the solvent molecules, the distribution of the solvent around the base pair becomes more homogeneous. Finally in the presence of 16 solvent molecules, corresponding roughly to the complete first solvation shell, all possible H-bond donors and acceptors of bases are employed in the interaction with the solvent, and the equal distribution of solvent around each atom of the bases was found. MD Simulations in CH3OH. We obtained almost identical results as for the MD simulation of base pairs in the presence of few water molecules38 showing a substantial similarity between the CH3OH and water solvents. A modest decrease of permittivity of the solvent seems to have only a small effect on the population of the HB and nHB structures. The abundances of the HB + nHB structures with CH3OH observed were only by a few percent higher than in the case of the population of the HB + nHB structures with water.

2596 J. Phys. Chem. B, Vol. 111, No. 10, 2007

Zendlova´ et al.

TABLE 3: Orientation of Bases, Total Populations of Structures in a Given Class of Nonmethylated and Methylated Guanine‚‚‚Cytosine Complexes in the Presence of 1-16 Molecules of Solvent (CH3OH, DMSO, and CHCl3) Obtained by MD Techniques Employing the Cornell et al. Force Field orientation of basesa/ number of solventb

H2O

CH3OH

CHCl3

DMSO

orientation of basesa/ number of solventb

H2O

CH3OH

CHCl3

DMSO

p np t s ns separated

GC 1 mol 43 13 34 73 22 13 0 0 1 0 0 0

12 74 14 0 0 0

10 68 22 0 0 0

p np t s ns separated

mGmC 1 mol 38 10 36 71 25 19 0 0 1 1 0 0

11 74 15 0 0 0

9 66 22 1 2 0

p np t s ns separated

GC 2 mol 36 9 33 65 28 24 1 0 2 2 0 0

12 73 15 0 0 0

9 66 24 0 1 0

p np t s ns separated

mGmC 2 mol 31 9 4 68 28 21 2 1 3 2 0 0

11 75 14 0 0 0

9 66 22 1 3 0

p np t s ns separated

GC 4 mol 30 8 29 60 2 28 3 1 5 4 1 0

12 75 13 0 0 0

7 54 33 1 5 0

p np t s ns separated

mGmC 4 mol 25 8 28 61 32 24 5 2 9 5 1 0

11 76 13 0 0 0

7 55 28 3 7 0

p np t s ns separated

GC 8 mol 12 4 15 38 30 36 10 5 32 17 1 0

12 75 13 0 0 0

6 47 41 1 5 0

p np t s ns separated

mGmC 8 mol 18 5 22 43 31 29 9 7 18 16 2 0

12 76 12 0 0 0

6 51 32 2 7 1

p np t s ns separated

GC 16 mol 5 2 7 21 27 36 12 11 43 30 4 0

7 60 32 0 1 0

3 34 55 1 5 2

p np t s ns separated

mGmC 16 mol 9 2 12 25 20 30 17 13 38 30 4 0

7 63 29 0 1 0

3 35 45 3 10 4

a

Classes differ by the mutual orientation of bases and the distances between the centers of masses of bases: p means planar H-bonded (the distance between the centers of masses of bases ranges between 5.0 and 10.0 Å, and the deviation of the base planes lies between 0° and 10°); np means nonplanar H-bonded (the same distance as for the first type, the deviation of base planes lies between 10° and 45°); t means T-shaped (the deviation of base planes lies between 45° and 90°); s means planar stacked (the distance between center of masses of bases is smaller than 5.0 Å, and the deviation of the base planes lies between 0° and 10°); ns means nonplanar stacked (the distance is the same as that for a planar stacked structure, and the deviation of the base planes lies between 10° and 45°). b “1 mol” means 1 molecule of solvent, “2 mol” means 2 molecules of solvent, etc.

The decisive factor for the explanation of the similarity of the results is the presence of the OH group in CH3OH as well as in water allowing the solvent to create the same binding patterns in both solvents. As in the case of complexes with water, we can conclude that an increase in the number of CH3OH molecules undoubtedly leads to a preference for S structures for the GC and AT complexes. More CH3OH molecules are necessary for the domination of S structures in the case of GC and mGmC systems than for AT and mAmT complexes, which is due to stronger H-bonding in the GC pair. In addition, owing to the absence of the second hydrogen, the CH3OH molecules do not create such a dense stabilizing network of H-bonds in S structures as in complexes with water. For both the methylated and the nonmethylated AT pair we found that the most preferred binding sites for one molecule of the solvent are to the N7 atom and amino group of adenine (see HB AT-1CH3OH in Figure 3) corresponding often to the simultaneous interaction of the solvent with both of these atoms. We used standard atom labeling for purine and pyrimidine bases (see also Figure 1aS in the Supporting Information). For the

nonmethylated GC pair the interaction of the hydrogen atom of the hydroxyl group of methanol with oxygen atoms of bases (O6 of guanine and O2 of cytosine, see T GC-1CH3OH in Figure 6) is preferred. The interaction with the N1 atom of cytosine (see HB GC1-CH3OH structure in Figure 6) is the third most important binding motif. Methylation of the GC pair leads to the preference of the O6 atom of guanine and the amino group of cytosine for binding of this solvent (see HB mGmC-2CH3OH in Figure 6). MD Simulations in DMSO. Since DMSO is a larger molecule than water or CH3OH and it contains a sulfoxyl group instead of a hydroxyl one, different binding patterns of this solvent with bases in comparison with water can be expected. Adding DMSO to the AT and mAmT base pairs leads to the preference for T structures, with S structures becoming the second most abundant ones. Conversely, in the cases of GC and mGmC complexes we can notice that HB complexes are dominant (corresponding to the Watson-Crick (WC) arrangement) without the influence of the increasing number of solvent molecules, thus verifying

Nucleic Acid Base Pairs in Organic Solvents

J. Phys. Chem. B, Vol. 111, No. 10, 2007 2597

TABLE 4: Frequency (in %) of Contacts of Each Possible Donor or Acceptor Atom of Nonmethylated and Methylated Adenine‚‚‚Thymine Bases with Molecules of Solvent during MD Simulationa adenine

TABLE 5: Frequency (in %) of Contacts of Each Possible Donor or Acceptor Atom of Nonmethylated and Methylated Guanine‚‚‚Cytosine Bases with Molecules of Solvent during MD Simulationa

thymine

guanine H-donors

H-donors H-acceptors H-donors H-acceptors complexes

N6

AT-1CH3OH AT-2CH3OH AT-4CH3OH AT-8CH3OH AT-16CH3OH

20 12 11 12 12

mAmT-1CH3OH mAmT-2CH3OH mAmT-4CH3OH mAmT-8CH3OH mAmT-16CH3OH

27 24 15 16 15

AT-1DMSO AT-2DMSO AT-4DMSO AT-8DMSO AT-16DMSO

4 12 23 30 25

mAmT-1DMSO mAmT-2DMSO mAmT-4DMSO mAmT-8DMSO mAmT-16DMSO

71 56 61 55 51

AT-1CHCl3 AT-2CHCl3 AT-4CHCl3 AT-8CHCl3 AT-16CHCl3

3 2 15 15 11

mAmT-1CHCl3 mAmT-2CHCl3 mAmT-4CHCl3 mAmT-8CHCl3 mAmT-16CHCl3

6 6 14 13 14

N9 N1 N3 N7 N1 10 15 13 11 13

2 3 5 6 10

10 9 10 10 10

20 11 11 11 11

13 17 13 11 11

4 3 17 5 4 17 6 6 17 12 10 14 14 14 14 30 39 31 25 25

62 40 30 26 25

N3

O2

O4

complexes

3 11 9 9 11

8 14 13 7 14 3 12 9 12 13 8 13 11 11 11 7 7 20 7 6 21 17 18 12 16 18 12 14 15 14

2 2 10 11 11

N4

N3

O2

2 4 6 5 6 4 8 7 4 6 4 7 10 6 6 7 11 9 6 6 10 10 10.0 10 10

11 13 14 13 10.

2 4 4 6 10

18 18 16 16 10

21 24 26 21 13

5 6 7 10 12

11 11 12 15 13

11 10 10 10 10

9 9 9 9 10

10 10 10 10 10

13 13 13 13 13

12 12 12 12 12

13 13 12 12 12

3 5 5 7 11

15 20 20 17 11

7 8 12 15 11

7 12 12 10 14

25 20 24 20 14

17 18 20 18 15

mGmC-1CH3OH 3 5 mGmC-2CH3OH 4 6 mGmC-4CH3OH 5 7 mGmC-8CH3OH 8 12 mGmC-16CH3OH 12 13

29 44 39 45 49 3 3 10 9 11

32 31 9 9 11

22 27 14 13 12

6 4 11 12 14

32 28 14 13 14

22 28 14 16 15

a

As a contact we accepted only distances shorter than 3.4 Å between the heavy atoms of the bases and the heavy atoms of the solvents (O in the case of CH3OH and DMSO, and C in the case of CHCl3). Interaction with carbon atoms of the bases was not considered due to a very low count of contacts.

the previous experimental conclusions.19-21 A broadening of the histogram in the presence of 8 or 16 molecules of the solvent with distances larger than 8 Å shows that also some structures of the base pairs separated by the solvent can be present. Since the DMSO interacts with bases almost exclusively via its oxygen atom, we focused on the distribution of the solvent molecules solely around H-bond donors of the bases. For the nonmethylated AT pair we found that one molecule of DMSO interacts mostly with the N1 atom of thymine (more than 60% of contacts during the simulation, see HB AT-1DMSO in Figure 4) and the N9 atom of adenine (30%, see HB AT2DMSO in Figure 4). Since both of these atoms are blocked after methylation, the most preferential binding site for the methylated AT pair becomes the amino group of adenine (70%, see S mGmC-1DMSO in Figure 4). For the nonmethylated GC pair there are three binding sites that are almost equally important: N9 of guanine, the amino group of cytosine (see HB GC-2DMSO in Figure 7), and N1 of cytosine (see HB GC-1DMSO in Figure 7). Due to the blocking of the first two named sites by methylation in the methylated GC pair (similarly to the AT pair) the interaction of one molecule of the solvent with the amino group of cytosine (see HB mGmC-1DMSO in Figure 7) is the most frequent (80%).

H-acceptors H-donors H-acceptors

N1 N2 N9 N3 N7 O6 N1

GC-1CH3OH GC-2CH3OH GC-4CH3OH GC-8CH3OH GC-16CH3OH

4 9 16 19 25

cytosine

3 3 3 3 12

14 11 10 9 12

27 21 20 16 10

19 13 13 10 10

38 35 30 23 13

GC-1DMSO GC-2DMSO GC-4DMSO GC-8DMSO GC-16DMSO

3 3 34 3 4 34 6 9 32 9 13 28 19 19 21

32 34 30 27 20

mGmC-1DMSO mGmC-2DMSO mGmC-4DMSO mGmC-8DMSO mGmC-16DMSO

8 11 19 23 32

GC-1CHCl3 GC-2CHCl3 GC-4CHCl3 GC-8CHCl3 GC-16CHCl3

8 9 12 9 9 11 8 9 12 8 10 11 10 10 10

9 10 10 10 10

11 9 11 9 11 9 11 9 10 10

mGmC-1CHCL3 mGmC-2CHCl3 mGmC-4CHCl3 mGmC-8CHCl3 mGmC-16CHCl3

11 11 11 11 12

13 13 13 13 13

14 14 14 14 13

9 18 31 37 34

12 12 13 13 13

28 25 23 23 21 83 71 50 40 34

12 12 12 12 12

12 12 12 12 10

a As a contact we accepted only distances shorter than 3.4 Å between the heavy atoms of the bases and the heavy atoms of the solvents (O in the case of CH3OH and DMSO, and C in the case of CHCl3). Interaction with carbon atoms of the bases was not considered due to a very low count of contacts.

MD Simulations in CHCl3. For all of the systems studied, we can conclude that the presence of CHCl3 leads to the preference of the HB complexes, confirming thus both the theoretical25,26 and the experimental13-16 results in the bulk solvent. We observed that the population of S structures decreases with the increasing number of CHCl3 molecules. This is in disagreement with previously accepted theory that the increasing permittivity of the surrounding leads to the preference for the S structures (relative permittivity in vacuum of 1 vs relative permittivity in CHCl3 of 4.8). This discrepancy can be explained because the polarity of the solvent is not a decisive factor for the preference of stacking or H-bonding between the bases, but rather specific interactions between the solvent and the bases are critical for this preference. Similarly like CH3OH, the CHCl3 molecule can interact with the atoms that are either H-bond donors or H-bond acceptors. For both the methylated and the nonmethylated AT pairs we found that CHCl3 prefers interaction via its hydrogen atom rather than by its chlorine atoms. The interaction with oxygen atoms of thymine (see HB AT-1CH3Cl in Figure 7) was observed more often than the interaction of the solvent with nitrogen atoms of adenine (see 49p AT-1CH3Cl in Figure 4aS in the Supporting Information). Surprisingly, a different situation was observed for both the methylated and the nonmethylated GC

2598 J. Phys. Chem. B, Vol. 111, No. 10, 2007

Zendlova´ et al.

Figure 3. Nonmethylated (upper two rows) and methylated (lower two rows) adenine‚‚‚thymine complexes with one or two molecules of CH3OH with a ball-and-stick representation of the most stable HB, T, and S structures optimized at the ab initio RI-MP2/cc-pVDZ level. The description below the structures shows the binding motif between bases; interaction energy in kcal/mol obtained as interaction energy extrapolated to the infinite basis set (deformation energies included, see also Tables 9 and 10). Oxygen atoms are red, nitrogen atoms are blue, carbon atoms are green, and hydrogen atoms are white. The H-bonds are depicted using dashed purple lines. For an easier comparison of the different structures the position of adenine is kept in the same orientation except the structures where such an orientation does not allow a clear view on the spatial arrangement of the bases.

Nucleic Acid Base Pairs in Organic Solvents

J. Phys. Chem. B, Vol. 111, No. 10, 2007 2599

Figure 4. Nonmethylated (upper two rows) and methylated (lower two rows) adenine‚‚‚thymine complexes with one or two molecules of DMSO with a ball-and-stick representation of the most stable HB, T, and S structures optimized at the ab initio RI-MP2/cc-pVDZ level. For a more detailed description see Figure 3. The sulfur atom is yellow.

pairs, where no preferential binding site was found already in the presence of one molecule of the solvent. Potential and Free Energy Surfaces of Base Pairs in the Presence of One or Two Molecules of Organic Solvent. Total numbers of structures for each system studied in the presence

of organic solvents and in a vacuum and their populations determined by the MD/Q method are shown in Tables 6-8. The four most stable HB, T, S as well as nS and nHB structures and their relative stabilities can be found in the Supporting Information section. Besides the location of the structural motifs,

2600 J. Phys. Chem. B, Vol. 111, No. 10, 2007

Zendlova´ et al.

Figure 5. Nonmethylated (upper two rows) and methylated (lower two rows) adenine‚‚‚thymine complexes with one or two molecules of CHCl3 with a ball-and-stick representation of the most stable HB, T, and S structures optimized at the ab initio RI-MP2/cc-pVDZ level. For a more detailed description see Figure 3. The chlorine atom is light green.

we also evaluated the populations of the investigated structures (Tables 6-8) and in this way considered the role of entropy. We found that the addition of one molecule of any solvent makes the PES much more complicated (having several hundreds of structures) than that for isolated base pairs (less than 20 structures);47 the presence of a second molecule of solvent leads to a further increase in the number of structural motifs by several times (thousands of structures). All systems studied with the presence of one or two DMSO molecules offer the largest variety of structural motifs when compared with those of the remaining two solvents (CH3OH and CHCl3) as well as with complexes with water.40,41 As could be expected, the minimization procedure leads in most cases to planarization of the complexes of the base pairs, which means that unlike with MD the planar S and HB complexes now become more populated than the nonplanar ones.

We received similar results for both the nonmethylated and the methylated AT and GC systems in the presence of CH3OH as for the hydrated AT and GC systems. We observed that the HB AT and HB GC structures were the most populated in the presence of one CH3OH molecule (40% of AT-1CH3OH and 68% of GC-1CH3OH). After the addition of a second CH3OH molecule, the population of the HB AT system drops, whereas the population of the HB GC systems does not change significantly. As shown38,40 for hydrated AT systems, the methylation of the base strongly influences the population of the HB AT motifs, contrary to the hydrated GC systems, where this effect on the HB population is negligible. According to this finding, the populations of the HB AT and mAmT structures in the presence of one or two CH3OH molecules upon methylation significantly decrease, and the S structures become the most populated (69% of mAmT-1CH3OH and 66% of mAmT-2CH3OH). In the case of methylated GC systems, the

Nucleic Acid Base Pairs in Organic Solvents

J. Phys. Chem. B, Vol. 111, No. 10, 2007 2601

Figure 6. Nonmethylated (upper two rows) and methylated (lower two rows) guanine‚‚‚cytosine complexes with one or two molecules of CH3OH with a ball-and-stick representation of the most stable HB, T, and S structures optimized at the ab initio RI-MP2/cc-pVDZ level. For an easier comparison of the different structures the position of guanine is kept in the same orientation except the structures where such an orientation does not allow a clear view on the spatial arrangement of the bases. For a more detailed description see Figure 3.

HB motifs remain the most abundant (55% of mGmC-1CH3OH and 48% of mGmC-2CH3OH), and the populations of the S structures do not exceed 20%. Although the DMSO molecule is a polar molecule quite like water or CH3OH, we observed significantly different populations

of HB, S, and T structures at the free energy surfaces (FESs) of AT and GC systems in the presence of this molecule in comparison with the above-mentioned solvents. It can be explained not only by the different chemical composition of DMSO but also by its size, which is larger than that of water

2602 J. Phys. Chem. B, Vol. 111, No. 10, 2007

Zendlova´ et al.

Figure 7. Nonmethylated (upper two rows) and methylated (lower two rows) guanine‚‚‚cytosine complexes with one or two molecules of DMSO with a ball-and-stick representation of the most stable HB, T, and S optimized at the ab initio RI-MP2/cc-pVDZ level. For a more detailed description see Figures 3 and 6. The sulfur atom is yellow.

or CH3OH. Like in the MD simulation, the T structures of the AT and GC systems were observed as the most populated at FESs, with the exception of that of mAmT-1DMSO and GC1DMSO, where the S (54%) and HB (44%) motifs were the most populated. The CHCl3 molecule is nonpolar in contrast to the polar CH3OH and DMSO solvents. For complexes containing only one or two molecules of CHCl3, we found the distribution of structural motifs closer to a vacuum than to the bulk solvent.

For the AT systems in the presence of one or two CHCl3 molecules, the most populated structures found were the T ones for the AT and the S ones for the mAmT complexes. The nonplanar HB structures were found as the most populated at all of the FESs of the GC systems with the exception of mGmC-2CHCl3, where the T motifs were the most populated. However, the trends corresponding to the situation in the bulk solvent (a decrease of abundance of the S structures and an increase in the HB ones) are already perceivable during a

Nucleic Acid Base Pairs in Organic Solvents

J. Phys. Chem. B, Vol. 111, No. 10, 2007 2603

TABLE 6: Total Populations of Structures in a Given Class of Nonmethylated and Methylated Adenine‚‚‚Thymine and Guanine‚‚‚Cytosine Complexes in a Vacuum Obtained by Molecular Dynamics/Quenching (MD/Q) Techniques Employing the Cornell et al. Force Field orientation of basesa

no.b

AT pop.c

p np t s ns separated total

9 1 10 6 1 0 27

57 0 23 19 1 0 100

mAmT no.b pop.c

no.b

GC pop.c

5 0 7 11 2 0 25

10 1 11 1 1 0 24

62 8 29 1 0 0 100

10 0 7 8 0 0 100

mGmC no.b pop.c 6 1 5 4 0 0 16

58 0 25 16 0 0 100

a Classes differ by the mutual orientation of bases and the distances between the centers of masses of bases: p means planar H-bonded (the distance between the centers of masses of bases ranges between 5.0 and 10.0 Å, and the deviation of the base planes lies between 0° and 10°); np means nonplanar H-bonded (the same distance as for the first type, the deviation of base planes lies between 10° and 45°); t means T-shaped (the deviation of base planes lies between 45° and 90°); s means planar stacked (the distance between center of masses of bases is smaller than 5.0 Å, and the deviation of the base planes lies between 0° and 10°); ns means nonplanar stacked (the distance is the same as that for a planar stacked structure, and the deviation of the base planes lies between 10° and 45°). b Total number of structures obtained by the MD/Q technique in the given class. c Total population (in %) of all structures of the given class obtained by the MD/Q technique.

transition from isolated base pairs in a vacuum to complexes with one or two molecules of the solvent. Geometrical Description of the Interaction between the Bases and between the Bases and the Solvent. The H-bond patterns between the bases found for the most stable HB complexes in a vacuum47 and with water40 were compared with those in complexes with CH3OH, DMSO, and CHCl3 molecule(s). We observed that they were in almost all cases identical and not very sensitive to the solvent. The binding patterns of the CH3OH molecule are reminiscent of those of the water molecule, which is a result of the chemical similarity of these two compounds. The water molecule itself can create H-bonds via all of its three atoms, whereas the CH3OH molecule can employ only the hydroxyl group for Hbonding. The smaller number of H-bonded centers in CH3OH does not become evident in HB structures, where the water molecule also usually uses only two of its atoms for H-bonds. However, it plays a significant role in T and S structures, where the water molecule usually exploits all three atoms for the H-bond network, connecting the bases by a solvent bridge. When CH3OH interacts with guanine, similar binding patterns of the solvent to bases were detected as had been observed in the work of Delchev and Mikosch.48 For DMSO the interaction with the bases in most stable complexes is mediated exclusively by the oxygen atom of the solvent. It concerns also the S complexes, where the oxygen atom is usually involved in bifurcated H-bonds connecting the hydrogen atoms of the two bases (see the S AT-1DMSO structure in Figure 4). While for the base pairs’ complexes with two water molecules (and partially also with CH3OH), we observed the tendency of the solvent to clusterize (interact with each other); for DMSO such an effect is excluded for steric reasons. We discovered three different binding motifs for the CHCl3 molecule in the HB complexes. In the first one the interaction is mediated by the H atom of the solvent (see HB mAmT1CHCl3 in Figure 5), and in the second, which is present only

TABLE 7: Total Populations of Structures in a Given Class of Nonmethylated and Methylated Adenine‚‚‚Thymine Complexes in the Presence of One or Two Molecules of Solvent (H2O, CH3OH, DMSO, and CHCl3) Obtained by MD/Q Techniques Employing the Cornell et al. Force Field orientation H2O CH3OH CHCl3 DMSO of basesa/ number of solventb no.c pop.d no.c pop.d no.c pop.d no.c pop.d p np t s ns separated total

67 38 115 67 20 30 337

31 11 30 26 2 0 100

AT 1 mol 93 40 77 12 275 21 164 24 55 2 18 0 682 100

80 47 330 187 34 41 719

17 13 52 16 1 1 100

121 201 609 317 129 25 1402

22 11 34 22 10 0 100

p np t s ns separated total

281 259 671 402 106 26 1745

33 10 24 31 2 0 100

AT 2 mol 384 17 582 14 1262 33 805 32 172 3 50 1 3255 100

281 806 2173 1053 117 90 4520

10 30 49 11 1 0 100

401 991 2984 1117 834 156 6483

8 17 44 14 15 2 100

p np T s ns separated total

33 31 86 96 18 9 273

mAmT 1 mol 6 29 12 3 61 3 8 261 12 78 253 69 5 55 4 0 12 0 100 671 100

21 39 382 358 28 16 844

1 5 40 54 1 0 100

46 157 620 574 107 19 1523

4 7 33 52 3 1 100

p np t s ns separated total

102 138 435 614 110 19 1418

mAmT 2 mol 4 138 7 4 464 7 12 1572 15 75 2076 66 6 374 5 0 37 0 100 4661 100

80 654 3008 2219 120 237 6318

2 9 41 47 1 1 100

213 988 3703 2734 440 279 8357

2 11 42 39 5 1 100

a Classes differ by the mutual orientation of bases and the distances between the centers of masses of bases: p means planar H-bonded (the distance between the centers of masses of bases ranges between 5.0 and 10.0 Å, and the deviation of the base planes lies between 0° and 10°); np means nonplanar H-bonded (the same distance as for the first type, the deviation of base planes lies between 10° and 45°); t means T-shaped (the deviation of base planes lies between 45° and 90°); s means planar stacked (the distance between center of masses of bases is smaller than 5.0 Å, and the deviation of the base planes lies between 0° and 10°); ns means nonplanar stacked (the distance is the same as that for a planar stacked structure, and the deviation of the base planes lies between 10° and 45°). b “1 mol” means 1 molecule of solvent, “2 mol” means 2 molecules of solvent, etc. c Total number of structures obtained by the MD/Q technique in the given class. d Total population (in %) of all structures of the given class obtained by the MD/Q technique.

rarely but is more stable than the first one, both H and Cl atoms of CHCl3 are employed in binding to bases (see the HB AT1CHCl3 structure in Figure 5). In the third one, which is energetically preferred, the whole CHCl3 molecule is localized above the planes of the bases, which leads to a significant distortion of the base pair from planarity (see the nHB AT1CHCl3 structure in Figure 5). The distortion can be compensated by the addition of a second molecule of the solvent when the first molecule is situated above and the second below the planes of the bases. Nonetheless, such complexes are less stable than the complexes with both molecules of CHCl3 situated above the planes of the bases. No secondary effect of the stabilization of S base pairs (like the creation of a “clip” H-bond network connecting two bases together by a water molecule40) was observed, which seems to be the main reason (rather than the

2604 J. Phys. Chem. B, Vol. 111, No. 10, 2007

Zendlova´ et al.

TABLE 8: Total Populations of Structures in a Given Class of Nonmethylated and Methylated Guanine‚‚‚Cytosine Complexes in the Presence of One or Two Molecules of Solvent (H2O, CH3OH, DMSO, and CHCl3) Obtained by MD/Q Techniques Employing the Cornell et al. Force Field orientation H2O CH3OH CHCl3 DMSO of basesa/ number of solventb no.c pop.d no.c pop.d no.c pop.d no.c pop.d p np t s ns separated total

81 35 126 30 16 0 292

68 11 18 2 2 0 100

GC 1 mol 77 68 75 10 240 16 44 3 36 2 5 0 477 100

81 49 231 37 30 15 443

20 46 33 1 0 0 100

91 147 394 63 54 0 749

44 18 28 5 5 0 100

p np t s ns separated total

257 259 625 204 114 0 1460

37 18 31 9 6 0 100

GC 2 mol 448 67 751 13 1527 12 391 6 266 3 14 0 3397 100

283 695 1734 174 131 59 3076

19 58 22 0 0 0 100

505 1433 3679 524 697 90 6928

15 25 45 6 8 1 100

p np t s ns separated total

39 26 84 53 12 0 217

mGmC 1 mol 70 47 55 7 38 4 14 193 19 8 126 19 1 30 2 0 4 0 100 438 100

33 27 183 125 8 9 385

6 45 38 10 1 0 100

40 129 376 199 80 4 828

19 24 27 20 9 0 100

p np t s ns separated total

128 180 438 285 90 0 1125

mGmC 2 mol 54 253 48 21 412 11 13 1408 20 9 988 17 5 230 3 0 4 0 100 3295 100

120 400 1708 746 69 56 3099

17 0 30 3 0 0 100

297 1363 3884 1619 865 132 8160

7 23 39 19 10 1 100

a Classes differ by the mutual orientation of bases and the distances between the centers of masses of bases: p means planar H-bonded (the distance between the centers of masses of bases ranges between 5.0 and 10.0 Å, and the deviation of the base planes lies between 0° and 10°); np means nonplanar H-bonded (the same distance as for the first type, the deviation of base planes lies between 10° and 45°); t means T-shaped (the deviation of base planes lies between 45° and 90°); s means planar stacked (the distance between center of masses of bases is smaller than 5.0 Å, and the deviation of the base planes lies between 0° and 10°); ns means nonplanar stacked (the distance is the same as that for a planar stacked structure, and the deviation of the base planes lies between 10° and 45°). b “1 mol” means 1 molecule of solvent, “2 mol” means 2 molecules of solvent, etc. c Total number of structures obtained by the MD/Q technique in the given class. d Total population (in %) of all structures of the given class obtained by the MD/Q technique.

polarity of the solvent) for the preference of the H-bonded structures in the bulk CHCl3. Ab Initio Calculations and Their Comparison with the Empirical Ones. The most stable HB, T, and S structures of the adenine‚‚‚thymine, 9-methyladenine‚‚‚1-methylthymine, guanine‚‚‚cytosine, and 9-methylguanine‚‚‚1-methylcytosine complexes with one or two molecules of CH3OH, DMSO, or CHCl3 and the nHB structures (with CHCl3 molecules only) obtained from the MD/Q calculations were selected for a further ab initio study. Only the global minima of each motif were considered, because the empirical force field provides results for the nucleic acid base pairs that are very close to the ab initio ones.40,41,49 All of these structures are shown in Figures 3-8. A Comparison of Empirical and Ab Initio Results. We found that the empirical results reproduce the accurate quantum

chemical calculations very well even in the presence of organic solvents. The structural motif is preserved during the ab initio optimization, and the root-mean-square (rms) deviations between optimized empirical and ab initio geometries lie between 0.2 and 0.6 Å. The differences between the bond length of the ab initio and empirical H-bonds in the HB motifs are not greater than 0.1 Å in most cases. The deviations between the planes of the bases (dP) of the AT planar HB structures show in most cases only a minor change during the ab initio optimization. The greatest differences (12-30°) were observed between the ab initio and empirical structures for the methylated structures with two molecules of solvents. Higher rms deviations as well as a more significant difference between the ab initio and the empirical dPs were observed for the GC complexes. The dPs are in a range of 10-30° for all of the GC complexes. The higher nonplanarity observed in the GC complexes originates from the fact that at the ab initio level the amino group of guanine is more nonplanar than that of adenine. The values of the torsion angles of the hydrogen of the amino groups of the optimized isolated guanine and adenine bases at the MP2/6-311G(2df,p) level of theory are -13.3°/ 39.2° and -15.3°/16.5°, respectively.50 These numbers seem, however, to be slightly overestimated in comparison with recent results employing larger basis sets.51 The empirical Cornell et al. force field yields planar amino groups. It means that the participation of the amino group of guanine in the H-bond pattern and further reoptimization at the ab initio level leads to nonplanarization of the GC WC motif, similar to what was observed for the isolated GC pair by Danilov and Kurita.52 Even from the energetic point of view, there is a good agreement between the empirical and ab initio data. The empirical force field correctly detects the global minima at the PESs as well as the relative differences in stability between the different arrangements of the base pairs in all of the solvents considered. The empirical stabilization energies lie between the cc-pVDZ and the TZVPP ab initio values, which means that they are still underestimated with respect to the complete basis set (CBS) limit. We are aware of the fact that in our molecular dynamics and molecular mechanics simulations we use the force field without a polarizability term. This can bring an inaccuracy to the results especially in DMSO solvent, where the role of polarizability is more important than for the other two organic solvents tested in our study. However, currently used force fields covering the polarizability term do not offer much better results than the nonpolarizable force fields. Further, comparing the empirical and ab initio results for complexes containing the DMSO solvent and complexes with other two solvents, we did not find a larger discrepancy either in qualitative agreement (the average rms deviations of ab initio and empirical optimized structures in DMSO solvent of 0.3 Å does not exceed the rms deviations in other solvents) or in quantitative agreement (the absolute average difference (AAD) between the empirical and the ab initio interaction energies). We found only a small increase of the AADs by approximately 1 kcal/mol for complexes with DMSO solvent in comparison with other complexes. It shows that the role of the polarizability (at least in microsolvated systems) is not as important as it was expected to be. Energetic Comparison of the Global Minima in Different Arrangements of the Bases in Organic Solvents. The structures and the ab initio energies (at the CBS limit) of the global minima can be found in Figures 3-8. The energetic properties

Nucleic Acid Base Pairs in Organic Solvents

J. Phys. Chem. B, Vol. 111, No. 10, 2007 2605

Figure 8. Nonmethylated (upper two rows) and methylated (lower two rows) guanine‚‚‚cytosine complexes with one or two molecules of CHCl3 with a ball-and-stick representation of the most stable HB, T, and S optimized at the ab initio RI-MP2/cc-pVDZ level. For more detailed description see Figures 3 and 6. The chlorine atom is light green.

of the complexes described in more detail as well as their comparison with the isolated base pairs can be found in Tables 9 and 10. Complexes with CH3OH. Like in the case of complexes of AT with water, we found that while in the case of a monosolvated AT pair the global minimum corresponds to the HB structure, in the case of disolvated complexes the S structure is the most stable. The preference for the S structures is undoubted for all mAmT complexes. For GC and mGmC complexes with CH3OH (except for GC-1CH3OH, where the stabilization energy of HB is significantly higher than that of S), the stability of the most stable S and HB structures detected was comparable. In all of the cases, T structures are less stable than the HB and S ones. Complexes with DMSO. The complexes with DMSO exhibit the highest stabilization energies of all of the systems studied. Surprisingly enough, apart from the fact that the solvents have different properties, the stability order of the most stable HB, S, and T structures remains the same as in CH3OH. Even though the S structures are in most cases the most stable, they are less populated, and thus, their role in the global context is less important. Complexes with CHCl3. The weakest interaction between the base pair and the solvent was found for CHCl3. In all of the

cases we found that the slightly nonplanar HB arrangement with the molecule(s) of the solvent situated above (below) corresponds to the global minimum. The planar structure of the base pair with a H atom of CHCl3 situated in the plane of the base pair is usually less stable by 1-4 kcal/mol than the nonplanar one (see the nonmethylated and methylated HB AT and GC structures in the presence of one or two molecules of CHCl3 in Figures 5 and 8). Also the most stable T and S structures are less stable. We also compared the stabilization energies of the AT, mAmT, CG, and mGmC complexes in a vacuum with the stabilization energies of those complexes in the presence of one or two molecules of a solvent to determine how the solvent contributed to the system’s stability. The highest increases in stability (more than 30 kcal/mol) were observed for both the nonmethylated and the methylated S structures with two molecules of CH3OH due to the mutual cooperation of these two molecules of the solvent, which can simultaneously interact with the bases. The size of a DMSO molecule has a crucial effect on the high stability of complexes with this solvent. The CHCl3 molecule is not capable of creating a network of H-bonds, as a result of which the increase of the stability in the complexes with CHCl3 was the smallest. On average (if we subtract the interaction between the solvent molecules), each molecule of

2606 J. Phys. Chem. B, Vol. 111, No. 10, 2007

Zendlova´ et al.

TABLE 9: Empirical and Ab Initio Interaction Energies of the Most Stable HB, T, and S Structures of the Nonmethylated and Methylated Adenine‚‚‚Thymine Complexes in a Vacuum and with One or Two Molecules of Solvent (CH3OH, DMSO, and CHCl3) ∆Eamberb

∆Edefc

HB (1) T (6) S (7)

-15.50 -11.40 -11.40

1.88 1.12 0.96

HB (1) T (9) S (2)

-13.30 -7.02 -13.20

3.42 2.61 1.25

HB (1) T (13) S (8)

-23.34 -20.89 -21.24

3.27 4.25 2.16

HB (6) T (28) S (1)

-29.52 -27.53 -30.66

3.87 3.28 3.95

HB (66) T (18) S (1)

-19.26 -21.11 -23.29

2.55 1.78 1.92

HB (257) T (101) S (1)

-25.38 -26.32 -29.01

2.29 2.52 2.06

HB (2) T (22) S (5)

-25.49 -23.60 -24.91

3.72 2.86 5.31

HB (134) T (20) S (5)

-34.79 -36.75 -38.28

5.33 6.10 7.54

HB (127) T (5) S (1)

-19.18 -22.19 -23.03

2.01 2.01 1.90

HB (506) T (33) S (1)

-29.39 -31.99 -33.76

3.33 4.42 3.63

HB (18) nHB (3) T (1) S (85)

-19.30 -21.02 -21.13 -16.49

2.07 1.82 3.22 1.06

HB (21) nHB (3) T (1) S (932)

-23.62 -24.11 -24.26 -20.56

1.70 2.33 2.08 1.31

HB (113) nHB (3) T (1) S (23)

-16.92 -21.23 -21.36 -18.68

2.08 1.70 1.80 0.73

HB (6) nHB (10) T (1) S (82)

-24.58 -24.12 -25.00 -22.98

1.68 1.11 2.01 1.19

motifa

∆EDztot d

∆ETztote

AT vacuum -15.92 -18.20 -7.29 -10.76 -6.57 -11.92 mAmT vacuum -10.39 -12.98 -4.51 -6.86 -8.48 -13.98 AT-1CH3OH -24.58 -29.17 -20.77 -25.17 -16.25 -23.56 AT-2CH3OH -31.43 -36.90 -26.69 -34.31 -28.54 -39.33 mAmT-1CH3OH -19.33 -24.07 -17.67 -23.38 -16.83 -25.64 mAmT-2CH3OH -25.65 -31.72 -23.14 -32.33 -25.34 -37.49 AT-1DMSO -26.85 -31.87 -19.62 -26.64 -6.80 -23.34 AT-2DMSO -36.15 -43.49 -36.64 -45.98 -34.75 -46.86 mAmT-1DMSO -17.70 -23.98 -17.33 -25.57 -18.48 -27.02 mAmT-2DMSO -26.97 -36.74 -25.26 -38.65 -28.04 -40.55 AT-1CHCl3 -20.23 -24.12 -18.46 -24.95 -5.25 -21.22 -12.33 -19.83 AT-2CHCl3 -18.75 -29.52 -21.47 -31.63 -22.22 -31.77 -17.34 -28.13 mAmT-1CHCl3 -17.43 -21.15 -16.06 -29.96 -15.94 -23.68 -13.97 -22.66 mAmT-2CHCl3 -19.54 -30.72 -21.93 -32.58 -21.24 -32.31 -19.12 -31.66

∆EaDztotf

∆EaTz totg

∆ECBStoth

-17.63 -10.83 -12.35

-18.84 -11.87 -13.65

-19.32 -12.42 -14.19

-12.74 -7.00 -15.35

-13.87 -7.69 -16.96

-14.33 -7.98 -17.63

-28.03 -24.23 -23.26

-30.41 -26.33 -25.66

-31.38 -27.20 -26.66

-33.42 -29.06 -38.17

-36.00 -31.74 -42.16

-37.06 -32.86 -43.82

-23.22 -23.06 -25.67

-25.48 -25.11 -28.21

-26.41 -25.96 -29.26

-31.13 -32.01 -37.05

-33.83 -35.02 -40.92

-34.94 -36.25 -42.54

-31.22 -26.55 -24.02

-33.39 -28.75 -25.50

-34.30 -29.70 -26.16

-43.12 -45.93 -47.07

-45.95 -49.06 -50.82

-47.17 -50.43 -52.41

-23.95 -25.94 -24.53

-26.05 -28.19 -29.83

-26.93 -29.13 -33.70

-37.25 -40.01 -41.02

-40.32 -43.58 -44.54

-41.63 -45.10 -46.05

-23.72 -24.82 -21.32 -20.56

-25.57 -27.11 -23.18 -22.60

-26.33 -28.06 -23.95 -23.47

-30.08 -29.13 -32.04 -29.42

-33.27 -41.79 -35.28 -32.60

-34.59 -47.50 -36.63 -33.95

-20.96 -21.74 -24.00 -23.53

-22.72 -26.54 -26.45 -25.82

-23.45 -28.62 -27.48 -26.79

-29.04 -30.23 -33.25 -32.43

-34.57 -37.33 -36.69 -36.51

-35.91 -40.44 -38.14 -38.24

a Binding motif between the bases, for a description see the text. Absolute energetic order of structures obtained by MD/Q techniques can be found in the parenthesis. b Empirical interaction energy in kcal/mol calculated as a difference of the energy of the complexes and subunits using the Cornell et al. force field. c Total deformation energy in kcal/mol obtained as a sum of deformation energies of subunits calculated at the RIMP2/cc-pVDZ level. d Basis set superposition error (BSSE)-corrected interaction energies (deformation energies included) obtained at the RI-MP2/ cc-pVDZ level (in kcal/mol) calculated by following equations: ∆Eint ) Ecomplex - (Esubunit1 + Esubunit2 + Esolvent1) for systems with one molecule of solvent and ∆Eint ) Ecomplex - (Esubunit1 + Esubunit2 + Esolvent1+ Esolvent2) for systems with two molecules of solvent. e BSSE-corrected interaction energies (deformation energies included) obtained at the RI-MP2/TZVPP level (in kcal/mol) calculated by the equations in footnote d. f BSSEcorrected interaction energies (deformation energies included) obtained at the RI-MP2/aug-cc-pVDZ level (in kcal/mol) calculated by the equations in footnote d. g BSSE-corrected interaction energies (deformation energies included) obtained at the RI-MP2/aug-cc-pVTZ level (in kcal/mol) calculated by the equations in footnote d. h Total interaction energy in kcal/mol (deformation energies included) obtained by RI-MP2 calculation and extrapolated to the infinite basis set.

Nucleic Acid Base Pairs in Organic Solvents

J. Phys. Chem. B, Vol. 111, No. 10, 2007 2607

TABLE 10: Empirical and Ab Initio Interaction Energies of the Most Stable HB, T, and S Structures of the Nonmethylated and Methylated Guanine‚‚‚Cytosine Complexes with One or Two Molecules of Solvent (CH3OH, DMSO, and CHCl3) motifa

∆Eamberb

∆Edefc

∆EDztotd

HB (1) T (6) S (8)

-24.60 -13.70 -12.43

3.59 2.12 2.66

-21.30 -14.17 -9.40

HB (1) T (4) S (7)

-25.20 -17.30 -13.60

3.66 1.88 1.30

HB (1) T (12) S (65)

-31.54 -28.73 -22.63

5.19 5.61 2.87

HB (1) T (10) S (5)

-39.74 -36.91 -37.49

6.49 5.49 5.43

HB (1) T (12) S (15)

-29.19 -26.50 -25.56

4.26 4.01 4.59

HB (2) T (10) S (4)

-40.18 -37.32 -38.03

5.90 5.45 5.16

HB (1) T (12) S (83)

-33.41 -28.99 -25.78

5.52 5.51 5.00

HB (1) T (21) S (25)

-46.03 -43.88 -43.86

6.70 6.59 8.54

HB (1) T (30) S (7)

-33.05 -26.73 -28.22

5.43 5.08 5.46

HB (1) T (15) S (8)

-43.73 -41.70 -42.03

8.43 7.26 7.85

HB (7) nHB (1) T (17) S (51)

-27.58 -30.62 -25.86 -20.14

4.91 5.08 4.13 1.65

HB (4) nHB (1) T (105) S (574)

-34.75 -35.36 -30.51 -24.61

4.53 5.60 4.82 2.27

HB (9) nHB (1) T (12) S (21)

-26.38 -32.37 -25.81 -22.70

4.13 3.15 3.99 1.61

HB (4) nHB (1) T (129) S (280)

-35.15 -35.73 -29.02 -25.97

4.62 2.03 4.77 1.91

GC

∆ETztote

∆EaDztotf

∆EaTztotg

∆ECBStoth

-24.59 -16.78 -15.69

-25.69 -16.15 -16.11

-26.57 -17.44 -17.69

-26.94 -17.97 -18.35

-25.50 -16.01 -16.52

-27.16 -17.05 -18.07

-27.84 -7.48 -18.72

-34.82 -31.68 -27.13

-37.63 -34.17 -29.77

-38.79 -35.19 -30.87

-45.42 -43.20 -44.98

-48.91 -46.72 -49.91

-50.34 -48.16 -51.92

-32.82 -30.78 -30.61

-35.06 -33.07 -33.88

-35.98 -33.88 -35.23

-44.73 -43.30 -44.39

-48.20 -46.85 -48.38

-49.62 -48.35 -50.02

-37.52 -34.37 -31.62

-40.19 -36.54 -34.41

-41.31 -37.47 -35.65

-51.46 -49.26 -52.52

-55.34 -51.93 -56.41

-56.95 -53.21 -58.15

-37.26 -34.25 -35.25

-39.95 -36.83 -37.97

-41.07 -37.94 -39.13

-44.70 -51.93 -51.53

-52.86 -55.22 -55.54

-56.23 -58.81 -60.84

-30.92 -31.74 -32.49 -25.66

-33.38 -34.58 -35.17 -28.03

-34.38 -35.76 -36.28 -29.03

-40.02 -42.29 -40.22 -34.64

-43.87 -46.07 -47.35 -37.96

-45.46 -47.64 -50.67 -39.35

-30.21 -33.87 -27.68 -26.91

-32.90 -36.67 -30.05 -31.80

-34.01 -37.83 -31.05 -31.85

-40.56 -46.31 -37.00 -36.37

-44.38 -50.11 -40.69 -39.80

-45.96 -51.68 -42.23 -41.24

mGmC -21.07 -24.36 -14.36 -18.00 -9.86 -16.12 GC-1CH3OH -30.38 -36.25 -26.72 -32.74 -19.67 -27.47 GC-2CH3OH -40.20 -45.68 -36.64 -44.48 -34.93 -47.09 mGmC-1CH3OH -28.44 -33.56 -26.26 -31.64 -22.66 -30.90 mGmC-2CH3OH -38.85 -46.18 -36.56 -44.56 -34.17 -45.07 GC-1DMSO -32.17 -38.40 -28.39 -34.59 -22.26 -31.18 GC-2DMSO -44.60 -52.91 -38.82 -48.88 -38.89 -52.10 mGmC-1DMSO -31.78 -37.99 -24.91 -33.99 -25.66 -34.96 mGmC-2DMSO -41.43 -50.26 -36.06 -49.26 -37.28 -51.35 GC-1CHCl3 -26.36 -31.70 -24.35 -32.05 -25.88 -32.83 -17.70 -25.07 GC-2CHCl3 -28.11 -39.94 -31.14 -41.91 -29.41 -40.05 -23.06 -33.76 mGmC-1CHCl3 -22.79 -30.32 -26.30 -33.53 -19.14 -27.23 -18.39 -26.27 mGmC-2CHCl3 -28.34 -40.40 -34.71 -46.85 -23.05 -35.86 -22.07 -34.82

a Binding motif between the bases, for a description see the text. Absolute energetic order of structures obtained by MD/Q techniques can be found in the parenthesis. b Empirical interaction energy in kcal/mol calculated as a difference of the energy of the complexes and subunits using the Cornell et al. force field. c Total deformation energy in kcal/mol obtained as a sum of deformation energies of subunits calculated at the RIMP2/cc-pVDZ level. d Basis set superposition error (BSSE)-corrected interaction energies (deformation energies included) obtained at the RI-MP2/ cc-pVDZ level (in kcal/mol) calculated by following equations: ∆Eint ) Ecomplex - (Esubunit1 + Esubunit2 + Esolvent1) for systems with one molecule of solvent and ∆Eint ) Ecomplex - (Esubunit1 + Esubunit2 + Esolvent1+ Esolvent2) for systems with two molecules of solvent. e BSSE-corrected interaction energies (deformation energies included) obtained at the RI-MP2/TZVPP level (in kcal/mol) calculated by the equations in footnote d. f BSSEcorrected interaction energies (deformation energies included) obtained at the RI-MP2/aug-cc-pVDZ level (in kcal/mol) calculated by the equations in footnote d. g BSSE-corrected interaction energies (deformation energies included) obtained at the RI-MP2/aug-cc-pVTZ level (in kcal/mol) calculated by the equations in footnote d. h Total interaction energy in kcal/mol (deformation energies included) obtained by RI-MP2 calculation and extrapolated to the infinite basis set.

2608 J. Phys. Chem. B, Vol. 111, No. 10, 2007 TABLE 11: Interaction Energies of the Most Stable Nonmethylated Adenine‚‚‚Thymine HB, T, and S Structures with One Molecule of Solvent (CH3OH, DMSO, and CHCl3) Extrapolated to the Infinite Basis Set

Zendlova´ et al. TABLE 12: Interaction Energies of the Most Stable Nonmethylated Guanine‚‚‚Cytosine HB, T, and S Structures with One Molecule of Solvent (CH3OH, DMSO, and CHCl3) Extrapolated to the Infinite Basis Set

motifa

∆ECBStotb

∆ECCSD(T) ∆EMP2c

∆ECBStot + [∆ECCSD(T) - ∆EMP2]d

motifa

∆ECBStotb

∆ECCSD(T) ∆EMP2c

∆ECBStot + [∆ECCSD(T) - ∆EMP2]d

HB (1) T (13) S (8)

-31.38 -27.20 -26.66

AT-1CH3OH 0.97 1.38 2.86

-30.41 -23.34 -23.80

HB (1) T (12) S (65)

-38.79 -35.19 -30.87

GC-1CH3OH -0.26 -0.16 1.08

-39.05 -35.35 -29.79

HB (2) T (22) S (5)

-34.30 -29.70 -26.16

AT-1DMSO 0.60 1.71 2.31

-33.70 -27.99 -23.85

HB (1) T (12) S (83)

-41.31 -37.47 -35.65

GC-1DMSO 0.18 0.40 2.93

-41.13 -37.07 -32.72

HB (18) nHB (3) T (1) S (85)

-26.33 -28.06 -23.95 -23.47

AT-1CHCl3 0.63 2.50 2.02 3.53

-25.70 -25.56 -21.93 -19.94

HB (7) nHB (1) T (17) S (51)

-34.38 -35.76 -36.28 -29.03

GC-1CHCl3 -0.24 1.09 0.50 2.01

-34.62 -34.67 -35.78 -27.02

a Binding motif between the bases, for a description see the text. Absolute energetic order of structures obtained by MD/Q techniques can be found in the parenthesis. b Total interaction energy in kcal/mol (deformation energies included) obtained by RI-MP2 calculation and extrapolated to the infinite basis set. c (∆ECCSD(T) - ∆EMP2) correction term obtained as the difference between CCSD(T)/6-31G*//RI-MP2/ cc-pVDZ and MP2/6-31G*//RI-MP2/cc-pVDZ interaction energies. d Final interaction energy in kcal/mol obtained as the interaction energy extrapolated to the infinite basis set (deformation energies included) and the (∆ECCSD(T)- ∆EMP2) correction term.

CH3OH contributes to the stability by approximately 12 kcal/ mol, DMSO by 15 kcal/mol, and CHCl3 by 9 kcal/mol. Effect of the Basis Set Size on the Reliability of the Ab Initio Interaction Energies. The interaction energies calculated with the various basis sets (cc-pVDZ, TZVPP, aug-cc-pVDZ, and aug-cc-pVTZ) as well as extrapolated to the CBS limit (a two-point extrapolation to the complete basis set limit was made from the aug-cc-pVDZ and the aug-cc-pVTZ energies) can be found in Tables 6 and 7. The stabilization energies obtained with the cc-pVDZ basis set are systematically much smaller than those obtained with the TZVPP basis set. The differences vary from 5 to 13 kcal/ mol depending on the type and number of molecules of the solvent. The variance is proportional to the system size; the smallest one was observed for nonmethylated bases in the presence of one molecule of methanol, while the largest one was observed for methylated base pairs interacting with two molecules of DMSO. It is thus evident that the DZ basis set is insufficient for a correct and accurate description of the stabilization energies. The differences between the TZVPP interaction energies and the values at the CBS limit ranged from 2 to 5 kcal/mol independently of the structural arrangement of the base pairs. It should be emphasized that the stability order of the structures calculated using the TZVPP basis set remains unchanged when passing to the CBS limit, which makes this basis set suitable for use in routine calculations of similar complexes. The aug-cc-pVDZ basis set provides results almost identical with those from the TZVPP basis set. The most exact (but also the most computationally demanding) interaction energies calculated employing the aug-cc-pVTZ basis set are typically only by 1-2 kcal/mol smaller than those at the CBS limit, and the difference does not exceed 3 kcal/mol even in the largest systems studied. Very similar dependencies on the size of the basis set were obtained also for the base pairs in the presence of water.41

a Binding motif between the bases, for a description see the text. Absolute energetic order of structures obtained by MD/Q techniques can be found in the parenthesis. b Total interaction energy in kcal/mol (deformation energies included) obtained by RI-MP2 calculation and extrapolated to the infinite basis set. c (∆ECCSD(T) - ∆EMP2) correction term obtained as the difference between CCSD(T)/6-31G*//RI-MP2/ cc-pVDZ and MP2/6-31G*//RI-MP2/cc-pVDZ interaction energies. d Final interaction energy in kcal/mol obtained as the interaction energy extrapolated to the infinite basis set (deformation energies included) and the (∆ECCSD(T)- ∆EMP2) correction term.

The (∆ECCSD(T) - ∆EMP2) correction terms determined solely for systems with one molecule of the solvent employing the 6-31G(d) (0.25) basis set are shown in Tables 11 and 12. The use of a rather small basis set is dictated by demanding CCSD(T) calculations. However, this basis set yields reliable results, which arises from the fact that the (∆ECCSD(T) - ∆EMP2) correction term is almost independent of the basis set size.53 As expected, the (∆ECCSD(T) - ∆EMP2) correction term is small (less than 1 kcal/mol in absolute value) for the HB and T systems independently of the type of the solvent. The (∆ECCSD(T) - ∆EMP2) correction term is positive and rather large (2-3 kcal/ mol) for the S structures, which also agrees with our previous results,41,53 showing thus a slight overestimation of the stability of S complexes by the MP2 method. Conclusions The Cornell et al. empirical force field properly describes the behavior of the nucleic acid base pairs not only in the vacuum and water but also in the various organic solvents. The MD simulations presented here performed with a small number of the solvent molecules have verified previous experimental and theoretical results from the bulk solvents, thus proving that whereas in CHCl3 the base pairs create H-bonded structures in water or CH3OH the S structures are preferred. The interaction of the base pairs with DMSO is less transparent, because, probably due to the extended size of the solvent, the structures preferred are mostly the T structures, at least in the case of AT and mAmT complexes. A significant deviation of the base pairs from planarity at room temperature was observed during all MD simulations while minimization led to planarization of the base pairs. The MD/Q technique in complexes with the presence of one or two molecules of the solvent verified the trends obtained from MD simulations. These results provide clear evidence for the fact that the preference of the S or HB structures of the base pairs in the solvent is determined not only by polarity and bulk properties of the solvent but primarily by specific interac-

Nucleic Acid Base Pairs in Organic Solvents tions of a small number of solvent molecules situated in the first solvation shell. The PESs of solvated systems are very complicated and cover from several hundreds to several thousands of structures. We revealed completely different interactions between the bases and the solvents studied. Whereas CH3OH (as well as water) can stabilize the S structures of the base pairs by a higher number of H-bonds than is possible in hydrogen-bonded pairs, the CHCl3 molecule lacks such a property, and the HB structures with molecule(s) of solvent situated above or below the base pair are preferred. The DMSO molecule is unique by its dimensions in comparison with other solvents, and the T structures are the most abundant ones. In terms of obtaining accurate results, the accurate ab initio correlated calculations confirmed the empirical results. It is, however, necessary to use extended basis sets (at least TZVPP or aug-cc-pVDZ) and consider the (∆ECCSD(T) - ∆EMP2) correction term. Acknowledgment. This work was a part of Project No. Z40550506 of the Institute of Organic Chemistry and Biochemistry, the Academy of Sciences of the Czech Republic, and was supported by Grants Nos. LC512 (the Ministry of Education of the Czech Republic), 203/05/009 (P.H.), 203/05/H001 (L.Z.), and KJB400550518 (M.K.) from the Grant Agency of the Czech Republic. Supporting Information Available: Standard numbering of atoms of adenine (A), thymine (T), guanine (G), and cytosine (C) and their methylated analogues and the four most stable motifs of the HB, nHB, T , S, and nS structures of nonmethylated and methylated adenine‚‚‚thymine and guanine‚‚‚cytosine complexes with one or two CH3OH, DMSO, and CHCl3 molecules obtained by MD/Q techniques. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Hobza, P.; Zahradnı´k, R.; Mu¨ller-Dethlefs, K. Collect. Czech. Chem. Commun. 2006, 71 (4), 443-531. (2) Mu¨ller-Dethlefs, K.; Hobza, P. Chem. ReV. 2000, 100 (1), 143167. (3) Saenger, W. Principes of Nucleic Acid Structure; SpringerVerlag: New York, 1984. (4) Chester, N.; Marshak, D. R. Anal. Biochem. 1993, 209 (2), 284290. (5) Herskovits, T. T.; Harrington, J. P. Biochemistry 1972, 11 (25), 4800-4811. (6) Cullis, P. M.; Wolfenden, R. Biochemistry 1981, 20 (11), 30243028. (7) Miller, J. L.; Kollman, P. A. J. Phys. Chem. 1996, 100 (20), 85878594. (8) Orozco, M.; Colominas, C.; Luque, F. J. Chem. Phys. 1996, 209 (1), 19-29. (9) Eksterowicz, J. E.; Miller, J. L.; Kollman, P. A. J. Phys. Chem. B 1997, 101 (50), 10971-10975. (10) Amutha, R.; Subramanian, V.; Nair, B. U. Chem. Phys. Lett. 2001, 335 (5-6), 489-495. (11) Giesen, D. J.; Chambers, C. C.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. B 1997, 101 (25), 5084-5088. (12) Li, J. B.; Cramer, C. J.; Truhlar, D. G. Biophys. Chem. 1999, 78 (1-2), 147-155. (13) Kyogoku, Y.; Lord, R. C.; Rich, A. J. Am. Chem. Soc. 1967, 89 (3), 496-504. (14) Binford, J.; Hollowary, D. M. J. Mol. Biol. 1968, 31, 91. (15) Miller, J. H.; Sobell, H. M. J. Mol. Biol. 1967, 24 (2), 345-350. (16) Sartorius, J.; Schneider, H. J. Angew. Chem., Int. Ed. Engl. 1996, 35 (21), 1446-1452.

J. Phys. Chem. B, Vol. 111, No. 10, 2007 2609 (17) Schweizer, M. P.; Broom, A. D.; Ts’o, P. O.; Hollis, D. P. J. Am. Chem. Soc. 1968, 90 (4), 1042-1055. (18) Ts’o, P. O. P.; Melvin, I. S.; Olson, A. C. J. Am. Chem. Soc. 85, 1963, 1289. (19) Katz, L.; Penman, S. J. Mol. Biol. 1966, 15 (1), 220-231. (20) Shoup, R. R.; Miles, H. T.; Becker, E. D. Biochem. Biophys. Res. Commun. 1966, 23 (2), 194-201. (21) Newmark, R. A.; Cantor, C. R. J. Am. Chem. Soc. 1968, 90 (18), 5010-5017. (22) Cieplak, P.; Kollman, P. A. J. Am. Chem. Soc. 1988, 110 (12), 3734-3739. (23) Dang, L. X.; Kollman, P. A. J. Am. Chem. Soc. 1990, 112 (2), 503-507. (24) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117 (19), 5179-5197. (25) Pranata, J.; Jorgensen, W. L. Tetrahedron 1991, 47 (14-15), 24912501. (26) Pranata, J.; Wierschke, S. G.; Jorgensen, W. L. J. Am. Chem. Soc. 1991, 113 (8), 2810-2819. (27) Pohorille, A.; Pratt, L. R.; Burt, S. K.; MacElroy, R. D. J. Biomol. Struct. Dyn. 1984, 1 (5), 1257-1280. (28) Pohorille, A.; Burt, S. K.; MacElroy, R. D. J. Am. Chem. Soc. 1984, 106 (2), 402-409. (29) Danilov, V. I.; Slyusarchuk, O. N.; Poltev, V. I.; Alderfer, J. L. J. Biomol. Struct. Dyn. 1997, 15 (2), 347-355. (30) Danilov, V. I.; Zheltovsky, N. V.; Slyusarchuk, O. N.; Poltev, V. I.; Alderfer, J. L. J. Biomol. Struct. Dyn. 1997, 15 (1), 69-80. (31) Zhanpeisov, N. U.; Leszczynski, J. J. Mol. Struct. (THEOCHEM) 1999, 487 (1-2), 107-115. (32) Zhanpeisov, N. U.; Leszczynski, J. J. Phys. Chem. A 1998, 102 (30), 6167-6172. (33) Moroni, F.; Famulari, A.; Raimondi, M. J. Phys. Chem. A 2001, 105 (7), 1169-1174. (34) Langlet, J.; Giessnerprettre, C.; Pullman, B.; Claverie, P.; Piazzola, D. Int. J. Quantum Chem. 1980, 18 (2), 421-437. (35) Sivanesan, D.; Babu, K.; Gadre, S. R.; Subramanian, V.; Ramasami, T. J. Phys. Chem. A 2000, 104 (46), 10887-10894. (36) Sivanesan, D.; Sumathi, I.; Welsh, W. J. Chem. Phys. Lett. 2003, 367 (3-4), 351-360. (37) Subramanian, V.; Sivanesan, D.; Ramasami, T. Chem. Phys. Lett. 1998, 290 (1-3), 189-192. (38) Kabela´cˇ, M.; Hobza, P. Chem.sEur. J. 2001, 7 (10), 2067-2074. (39) Kabela´cˇ, M.; Ryja´cˇek, F.; Hobza, P. Phys. Chem. Chem. Phys. 2000, 2 (21), 4906-4909. (40) Kabela´cˇ, M.; Zendlova´, L.; R ˇ eha, D.; Hobza, P. J. Phys. Chem. B 2005, 109 (24), 12206-12213. (41) Zendlova´, L.; Hobza, P.; Kabela´cˇ, M. ChemPhysChem 2006, 7 (2), 439-447. (42) Fox, T.; Kollman, P. A. J. Phys. Chem. B 1998, 102 (41), 80708079. (43) Halkier, A.; Helgaker, T.; Jørgensen, P.; Klopper, W.; Koch, H.; Olsen, J.; Wilson, A. K. Chem. Phys. Lett. 1998, 286 (3-4), 243-252. (44) Ahlrichs, R.; Bar, M.; Haser, M.; Horn, H.; Kolmel, C. Chem. Phys. Lett. 1989, 162 (3), 165-169. (45) Werner, H. J.; Eckert, F.; Hampel, C.; Heltzer, G.; Korona, T.; Lindh, R.; Lloyd, A. W.; McNicholas, S. J.; Manby, F. R.; Meyer, W.; Mura, M. E.; Knowles, P. J.; Nicklass, A.; Palmieri, P.; Pitzer, R.; Rauhut, G.; Schu¨tz, M.; Stoll, H.; Stone, A. J.; Tarroni, R.; Thorsteinsson, T.; Amos, R. D.; Bernhardsson, A.; Berning, A.; Celani, P.; Cooper, D. L.; Deegan, M. J. O.; Dobbyn, A. J. MOLPRO, version 2002.6. http://www.molpro.net. (46) Sˇ pirko, V.; Sˇ poner, J.; Hobza, P. J. Chem. Phys. 1997, 106 (4), 1472-1479. (47) Kabela´cˇ, M.; Hobza, P. J. Phys. Chem. B 2001, 105 (24), 58045817. (48) Delchev, V. B.; Mikosch, H. Monatsh. Chem. 2004, 135 (11), 1373-1387. (49) Hobza, P.; Huba´lek, F.; Kabela´cˇ, M.; Mejzlik, P.; Sˇ poner, J.; Vondra´sˇek, J. Chem. Phys. Lett. 1996, 257 (1-2), 31-35. (50) Bludsky´, O.; Sˇ poner, J.; Leszczynski, J.; Sˇ pirko, V.; Hobza, P. J. Chem. Phys. 1996, 105 (24), 11042-11050. (51) Wang, S.; Schaefer, H. F., III. J. Chem. Phys. 2006, 124 (4), 1-8. (52) Kurita, N.; Danilov, V. I.; Anisimov, V. M. Chem. Phys. Lett. 2005, 404 (1-3), 164-170. (53) Jurecˇka, P.; Hobza, P. J. Am. Chem. Soc. 2003, 125 (50), 1560815613.