Complexes of Pentahydrated Zn2+ with Guanine, Adenine, and the

J. Phys. Chem. B 1999, 103, 11415-11427

11415

Complexes of Pentahydrated Zn2+ with Guanine, Adenine, and the Guanine-Cytosine and Adenine-Thymine Base Pairs. Structures and Energies Characterized by Polarizable Molecular Mechanics and ab Initio Calculations Nohad Gresh*,† Unite´ de Pharmacochimie Mole´ culaire et Structurale, U266 INSERM, UMR 8600 CNRS, U.F.R. des Sciences Pharmaceutiques et Biologiques, 4, aVenue de l’ObserVatoire, 75270 Paris Cedex 06, France

Jıˇrı´ Sˇ poner‡ J. HeyroVsky´ Institute of Physical Chemistry, Academy of Sciences of the Czech Republic, DolejsˇkoVa 3, 182 23 Prague, Czech Republic, and Institute of Biophysics, Academy of Sciences of the Czech Republic, Kra´ loVopolska´ 135, 612 65 Brno, Czech Republic ReceiVed: June 25, 1999; In Final Form: October 19, 1999

We investigate the binding of Zn2+ and pentahydrated Zn2+ to guanine and adenine and to the WatsonCrick G-C and A-T base pairs using the SIBFA molecular mechanics and ab initio Hartree-Fock and MP2, as well as density functional theory approaches. The ab initio computations use four basis sets: a 4-31G+(2d) coreless effective potential, the 6-31G* basis set combined with a pseudopotential description of the zinc cation, and the 6-311G** and 6-31G+(2d,2p) basis sets. For the G and G-C complexes, our computations show that the previously published structures were not global minima on the potential energy surface. The calculations demonstrate a very significant flexibility in the modes of binding of hydrated Zn2+ to nucleobases. This cation can adopt coordination numbers in the 4-6 range and bind to either N7 or O6 while the binding energies vary by small amounts. The SIBFA computations can reproduce the values of the ab initio binding energies using the CEP 4-31G+(2d) with an accuracy of 3% and can correctly account for the significant cooperative character (-15 kcal/mol) of Zn2+ binding to the guanine-cytosine base pair. The present results show this procedure to be adequate for computations on the complexes of divalent cations with large oligonucleotides, which cannot be carried out by ab initio calculations.

Introduction The interaction of Zn2+ with nucleic acids can significantly influence their structure and function. Like all metal cations, Zn2+ can interact electrostatically with negatively charged phosphate groups (reviewed in ref 1). However, in contrast to Mg2+ which has a similar ionic radius, Zn2+ often binds to nucleobases, mainly to the N7 position of purines or other nitrogenous sites.2,3 The difference between Zn2+ and Mg2+ is due to their different electronic structure. Quantum-chemical calculations have convincingly demonstrated that Zn2+ prefers to interact with nitrogen sites of nucleobases due to the attractive 3d electron-lone pair interaction4,5 which is absent for Mg2+. This picture is consistent with crystal database studies showing an abundant binding of Zn2+ to nitrogen atoms.6 Its binding to nucleobases can destabilize DNA and cause local denaturation of the double helix. On the other hand, some three-dimensional DNA architectures can be stabilized by the presence of this cation. These include intramolecular triplexes,7 hairpins,8 homopurine duplexes,9 and formation of kinks.10 Zn2+ can also influence formation of protein-DNA complexes, although here its influence on the protein rather than DNA structure can be essential.11

In the present paper, we demonstrate that the polarizable SIBFA molecular mechanics force field12 is capable of reproducing all major aspects of interactions between Zn2+ and DNA bases and base pairs. Therefore, this procedure can be used in studies of complexes between metal cations and large DNA fragments which would be prohibitively costly by ab initio computations. In the present study, several new aspects of Zn2+ binding to DNA bases were revealed with the aid of polarizable molecular mechanics. It should be noted that development of nonadditive force fields which correctly describe complexes involving metal cations is important since the currently widely used pair-additive force fields are in general not suitable for such treatments. Thus, it was recently demonstrated that the deficiency in the force field description of guanine-Na+ interactions leads to clear artifacts of H-bonding of bases in the course of nanosecond scale molecular dynamics simulations of DNA quadruplexes.13 Errors in pairwise additive potentials, upon handling the complexes of nucleic acids with divalent transition metal cations such as Zn2+, would be by an order of magnitude larger than those for monovalent alkali metal cations, rendering pairwise additive force fields useless. Procedures

†

Tel: (33) 01-45.73.96.89. Fax: (33) 01-43.26.69.18. E-mail: gresh@ bisance.citi2.fr. ‡ Tel: 420-2-6605-2056. Fax: 420-2-858-2307. E-mail: sponer@indy. jh-inst.cas.cz.

1. Ab Initio Computations. We use the following basis sets: (1) the coreless effective potential 4-31G set derived by Stevens et al.,14 supplemented on the heavy atoms by two

10.1021/jp9921351 CCC: $18.00 © 1999 American Chemical Society Published on Web 12/07/1999

Gresh and Sˇ poner

11416 J. Phys. Chem. B, Vol. 103, No. 51, 1999 diffuse, uncontracted 3d orbitals. This is the basis set initially used for the calibration of the SIBFA procedure; (2) the 6-31G* basis set which was previously used in several studies devoted to the binding of divalent cations to the nucleobases;4,5,15 (3) the 6-311G** basis set; and (4) the 6-311 + G(2d,2p) basis set. For Zn2+, two 4f diffuse functions are included, and the overall contraction is (5111111111/51111/311/11). This will enable us to evaluate both energy and configurational dependencies of the results upon the size of the basis as well as to compare the trends of the MP216 calculations to the density functional theory ones. Computations using basis sets 2-4 are carried out with the Gaussian 94 package.17 Because of the high cost of the computations, the MP2 computations with basis set 4 were limited to five out of the 11 investigated guanine complexes. The energy decompositions using the CEP 4-31G+(2d) basis set are carried out for three representative complexes using the restricted variational space approximation method (RVS) of Stevens and Fink18 using the Gamess package.19 This analysis also provides values for the basis set superposition error (BSSE)20 at the HF level, as computed with the virtual orbital space.21 No BSSE corrections are done at the correlated level. 2. DFT Computations. The DFT computations are done with the B3LYP functional,22 using basis sets 3 and 4, and the Gaussian 94 package. DFT computations were very recently reported by one of us upon investigating of the binding of hydrated IIa and IIb group metal cations to the thioguaninecytosine base pair.15 3. SIBFA Computations. The SIBFA intermolecular interaction energy is computed as a sum of five contributions (see ref 12 for details):

∆E ) EMTP + Erep + Epol + Ect + Edisp EMTP denotes the electrostatic (multipolar) energy contribution, computed with distributed multipoles derived from the ab initio SCF wave function of the constitutive fragments. The multipoles (up to quadrupoles) are distributed on the atoms and bond barycenters using the procedure developed by Vigne´Maeder and Claverie.23 This procedure enables a close reproduction of the ab initio electrostatic potential around biologically relevant molecules.23 It was used shortly after it was available, during the inception of SIBFA,24 and retained in the subsequent applications and developments of SIBFA. Erep is the short-range repulsion energy, computed as a sum of bond-bond, bond-lone pair, and lone pair-lone pair interactions. Epol is the polarization energy contribution, computed with distributed, anisotropic polarizabilities on the individual molecules. The polarizabilities are distributed on the centroids of the localized orbitals (heteroatom lone pairs and bond barycenters) using a procedure due to Garmer and Stevens.25 A Gaussian screening of the polarizing field is used. The field polarizing each molecule is computed with the permanent multipoles and the induced dipoles of all the other molecules in an iterative fashion. In addition to Epol, we calculate also Epol*, which is the polarization energy without inclusion of the contribution of the induced dipoles to the electrostatic field. Denoting by E° and E, respectively, the complete fields due to the permanent multipoles and the permanent multipoles + induced dipoles, and by RP the polarizability of centroid P

which denotes a lone pair or a bond barycenter, the expression of Epol(P) is

Epol(P) ) -0.5

∑i EP0(i)∑j RP(i,j) EP(j)

Ect is the charge-transfer energy contribution, in which a coupling with the polarization is accounted for,12c and Edisp is the dispersion energy contribution. The distributed multipoles and polarizabilities of formate and water are derived from ab initio calculations with the CEP 4-31G+(2d) basis set. Standard geometries were used, derived from X-ray crystal structures of DNA26 and used in previous simulations of oligonucleotides and their complexes with drugs27 and divalent cations.12 Part of the starting positions for SIBFA energy minimization was taken from the previous HF/6-31G* ab initio studies.15 Additional novel structures were generated with the help of molecular graphics. For these the criteria are those of preferential Zn2+ binding to either N7 or O6 as well as the imposition of a starting coordination number for this cation in the 4-6 range. The energy minimizations using SIBFA are done with the Merlin minimizer28 and bear on the intermolecular variables. The energy-minimized positions are used for single-point ab initio energy computations using the CEP 4-31G+(2d) basis set and basis sets 3 and 4. The solvation energies ∆Hsolv are computed using the Langlet et al. Continuum procedure29 interfaced in the SIBFA software.30 The second set of ab initio calculations are done at the SCF/ 6-31G* level while the Zn2+ cation is described by Christiansen relativistic pseudopotential. Starting from the SIBFA energyminimized structures, full optimization is performed including all intramolecular and intermolecular degrees of freedom as in our previous studies.4,5,15 It results in final structures closely similar to the starting ones, except in three instances indicated below. However, in contrast to the previous papers, the complexation energies reported in this paper, and using this basis set as well as basis sets 3 and 4, are obtained as the energy difference between the complex and isolated monomers relaxed in the monomer-centered basis set. In other words, these values include the deformation energies of monomers but they are not corrected for the basis set superposition error. It does not mean that the basis set superposition error is negligible for the present system.5b,15 Nevertheless, for the purpose of the present paper and for obtaining the relative values the basis set superposition error correction is not inevitable. Further, the quantitative accuracy of the calculations is limited also by other factors such as size of the basis set, neglect of electron correlation effects, and uncorrected basis set superposition error during the optimization. On the other hand, we wish to underline that concerning the size of the basis set for the presently studied cluster the 6-31G* basis set provides quite reasonable results. The actual performance of the 6-31G* basis set is for large clusters much better for large clusters than one could assume on the basis of studies carried out on very small systems. For the complexes investigated here, this point is illustrated by the comparisons reported between the numerical values of ∆EHF with this basis and those with extended basis sets 3 and 4. Results and Discussion A. Binding of Pentahydrated Zn2+ to Guanine and Adenine. The results are reported in Tables 1A,B and 2A,B. Tables 1A and 2A list for the complexes of guanine and adenine,

Complexes of Pentahydrated Zn2+

J. Phys. Chem. B, Vol. 103, No. 51, 1999 11417

TABLE 1: Binding Energies (kcal/mol) of Pentahydrated Zn2+ to Guanine (See Text for Definitions) A. Calculations at the Uncorrelated Level (1) Direct Binding Modes a ab initio E1 Epol Epol* Ect E2 E1 + E2 ∆E(SIBFA) ∆EHFa/ BSSE ∆EHFb ∆EHFc ∆EHFd ∆EHFe ∆EHFf

b SIBFA

ab initio

c SIBFA

ab initio

-227.7 -235.0 -135.5 -107.9 (-139.7) -26.9 -19.0 -162.4 -385.1 -361.9

-240.9 -249.6 -123.7 -99.3 (-126.0) -25.3 -16.9 -149.0 -389.7 -365.8

-356.4

-363.0

-361.5 -355.1 -378.5 -381.0 -384.9

-367.4 -363.4 -387.5 -386.6 -387.0

d SIBFA

ab initio

e SIBFA

ab initio

f SIBFA

ab initio

g SIBFA

ab initio

SIBFA

-209.1 -124.6 (-149.9) -23.3

-200.1 -134.8 (-153.5) -27.2

-189.0 -138.3 (-158.5) -27.5

-217.6 -121.3 (-149.6) -23.1

-208.1 -120.5 (-146.9) -20.0

-356.9

-362.1

-354.8

-362.0

-348.6

-358.2 -353.1 -375.0 -376.8 g

-360.4 -357.5 -377.8 -379.2 -387.9

-360.7 -356.6 -377.4 -379.0 h

-364.3 -358.5 -381.4 -384.4 -387.9

-359.3 -354.1 -374.9 -376.7 -386.6

(2) Through-Water Binding Modes h

E1 Epol Epol* Ect E2 E1 + E2 ∆E(SIBFA) ∆EHFa ∆EHFb ∆EHFc ∆EHFd ∆EHFe ∆EHFf

i SIBFA

-196.4 -138.1

-213.1 -114.5 (-133.9) -25.1

-204.3 -119.1 (-133.9) -25.8

-219.9 -106.7 (-131.1) -23.4

-352.6

-349.3

-350.0

-25.0 -163.1 -368.0

ab initio

j

ab initio

-346.9 -351.7 -346.8 -369.1 -371.1 -380.0

SIBFA

ab initio

-348.6 -345.0 -366.6 -367.2 -375.7

SIBFA

-346.6 -340.0 -363.2 -364.9

B. MP2 and DFT Calculations and SIBFA Calculations Including the Dispersion Energy (1) Direct Binding Modes a

b

c

d

e

f

g

ab initio SIBFA ab initio SIBFA ab initio SIBFA ab initio SIBFA ab initio SIBFA ab initio SIBFA ab initio SIBFA ∆Etot ∆EMP2i ∆EMP2jDEMP2b ∆EMP2k ∆EDFTj ∆EDFTk

-409.5 -407.9 -390.0 -412.4 -387.6 -423.1

-403.7 -403.3 -388.1 -411.0 -390.4 -424.9

-403.7 -400.8

-405.4 -395.5 -387.7 -407.2 -394.0 -426.0

-407.6 -388.6 -422.3

-400.0 -400.2 -390.3 -409.8 -395.7 -428.1

-408.9 -407.8 -392.9

-392.0 -401.0

-414.4 -394.2 -429.4

-406.5 -388.7 -421.4

(2) Through-Water Binding Modes h ab initio ∆Etot ∆EMP2i ∆EMP2j ∆EMP2k ∆EDFTj ∆EDFTk

i SIBFA

ab initio

-397.6

j SIBFA

ab initio

-392.3

SIBFA -394.3

-391.4

-385.6

-389.5

-401.3 -382.2 -417.0

-396.6 -379.7 -413.3

-396.2 -375.5 -410.4

a CEP 4-31G +(2d) basis set. Result from the RVS procedure including the BSSE correction with the virtual orbitals. b CEP 4-31G +(2d) basis set. Result from the SCF procedure without the BSSE correction. c 6-311 +G(2d,2p)** basis set. d 6-311G** basis set. e 6-31G* basis set. Singlepoint computation at the SIBFA-minimized position. f 6-31G* basis set. Energy-minimized position. g Energy-minimization led to structure f. h Energyminimization led to structure g. i MP2 calculation with the CEP 4-31G +(2d) basis set. j MP2 or DFT calculation with the 6-311 +G(2d,2p) basis set. k MP2 or DFT calculation with the 6-311G** basis set. i-k are single-point calculations at the SIBFA-energy-minimized position.

respectively, the SIBFA interaction energies, ∆E, and their components without the dispersion energy, and the corresponding values of the ab initio binding energies at the HartreeFock level. Tables 1B and 2B list for the G and A complexes the values of the total SIBFA binding energies now including

Edisp, ∆Etot, and the corresponding MP2 as well as DFT results. To facilitate the comparison with the SIBFA results, we list first the CEP 4-31G(2d) binding energies followed by the 6-31G-type basis sets in their order of decreasing size, namely 4, 3, and 2.

Gresh and Sˇ poner

11418 J. Phys. Chem. B, Vol. 103, No. 51, 1999 TABLE 2: Binding Energies (kcal/mol) of Pentahydrated Zn2+ to Adeninea A. Calculations at the Uncorrelated Level a ab initio E1 Epol Epol* Ect E2 ∆E(SIBFA) ∆EHFb ∆EHFc ∆EHFd ∆EHFe

b SIBFA

c

ab initio

SIBFA

-204.0 -97.7 (-130.5) -16.6

SIBFA

-187.6 -109.7 (-138.4) -20.6

-318.3 -321.0 -314.9 -338.6 -348.0

ab initio

ab initio

SIBFA

-185.8 -101.2 (-123.4) -23.5

-317.9 -322.8 -318.0 -340.8 -352.8

d -187.2 -100.3 (-122.0) -22.8

-310.5 -310.5 -306.4 -328.2 -340.0

-310.3 -308.8 -304.4 -326.3 -340.0

B. MP2 and DFT Calculations and SIBFA Calculations Including the Dispersion Energy a ab initio ∆Etot ∆EMP2f ∆EMP2h ∆EDFTg ∆EDFTh

b SIBFA

ab initio

-361.4 -368.2 -373.2 -363.7 -382.1

c SIBFA

ab initio

-359.9 -366.8 -374.8 -370.1 -388.0

d SIBFA

ab initio

-348.8 -349.0 -359.8 -357.5 -374.0

SIBFA -347.2

-349.7 -359.6 -355.8 -372.9

a a and b are direct binding modes of Zn2+ to adenine, and c and d are through-water binding modes b CEP 4-31G +(2d) basis set. Result from the SCF procedure without the BSSE correction. c 6-311 +G(2d,2p)** basis set. a-c are single-point calculations at the SIBFA energy-minimized position. d 6-311G** basis set. e 6-31G* basis set. Energy-minimized position. fMP2 calculation with the CEP 4-31G +(2d) basis set. g MP2 or DFT calculation with the 6-311 +G(2d,2p) basis set. h MP2 or DFT calculation with the 6-311G** basis set.

As in ref 31, only single-point computations are performed using the CEP 4-31G+(2d) and basis sets 3 and 4 with the SIBFA geometries. For the 6-31G* computations, two rows of results are provided: the first gives the ∆E(SCF) values using the same coordinates as resulting from the SIBFA energyminimization. The final row of data refers to the HF/6-31G* optimized structures. 1. Guanine. Seven representative direct-binding mode complexes, a-g, and three through-water ones, h-k, are investigated and their structure optimized using SIBFA. Zn2+ adopts a sixfold coordination in the first two direct binding modes a and b, as well as in the three through-water binding modes. Fivefold coordinated Zn2+ occurs in complexes c, f, and g, while in complexes d and e it adopts a coordination number of four. Single-point energy 6-31G* computations give binding energies that are larger than the CEP 4-31G+(2d) ones by virtually constant amounts, in the 18-20 kcal/mol range for majority of the guanine complexes. HF/6-31G* energy minimizations result into close energy equalization (to within 4 kcal/mol) between the structures with Zn2+ directly bound to G on the one hand, and those in which it is bound to G through water on the other hand. In one case a new minimum was obtained by the HF/631G* optimization and then used as a new starting point for SIBFA but retaining standard internal SIBFA geometries (complex d). In another case (complex c defined below) the HF/6-31G* gradient optimization led to another minimum, actually coincident with f. Due to the enormous flexibility of the hydration shells some differences of structures predicted by the SIBFA and HF/6-31G* optimizations must necessarily occur. It is due to the differences in energy description and optimization procedures. Note, that the HF/6-31G* optimizations completely relax the monomer geometries and one can assume that it provides more flexibility for interconversions among structures. We did not detect, however, any qualitative difference between predictions obtained by the two methods. (a) Direct Binding Modes (Figure 1a-g). In structure a, the metal cation is bound to N7, and two of its solvation waters act as hydrogen-bond donors to O6.

In b, Zn2+ is bound to O6, and one solvation water is a hydrogen-bond donor to N7. Both ∆E(SIBFA) and ∆EHF give this structure as the most stable direct-binding mode although the energy preference is less than 2% of the total binding energy. c is similar to a, with Zn2+ bound to N7, but has one solvation water in the second shell so that Zn2+ is fivefold coordinated. The outer-shell water W5 is an H-bond acceptor with respect to one inner-shell water W4 but with a very short (d ) 1.59 Å) O-H distance. As noted in our previous work devoted to formate complexes with pentahydrated Zn2+,31 this is indicative of strong local cooperative effects along the Zn2+-W4-W5 path. c is 9 kcal/mol less favorable than b following both SIBFA and SCF computations with the CEP 4-31G+(2d) set. Inclusion of Edisp in SIBFA nullifies this difference. Structure d initially derived by HF/6-31G* optimization is similar to b but has two outer-shell water molecules, so that Zn2+ is fourfold coordinated. The very short H-bond distance between water W1 and N7 (1.64 Å) is indicative of cooperative effects beween O6, Zn2+ and W1. Recent analysis showed cooperative effects to be enhanced when water molecules could act simultaneously as H-bond donors and acceptors to two congeners.32 This appears to be the case in Zn2+ complexes as well in Zn2+-bound water molecules which can behave simultaneously as H-bond donors to O6 or N7 sites of guanine. The SIBFA computations indicate that in terms of ∆E, d is less stable than b by only 3.7 kcal/mol out of 350. The corresponding SCF energy difference with the CEP 4-31G+(2d) basis set is somewhat larger (6.5 kcal/mol) but this represents only 2% of the total binding energy. Including the contribution of Edisp in SIBFA reverts the trend and yields a small preferential stabilization of d over b. Similar to d, structure e also has a tetracoordinated Zn2+, the cation now being bound to N7, whereas one coordinating water is an H-bond donor to O6, with a very short dO-H distance of 1.61 Å. f is similar to c: Zn2+ is bound to N7 and four water molecules, but the fifth, outer-shell molecule W5 is an H-bond

Complexes of Pentahydrated Zn2+

Figure 1. Binding of [Zn(H2O)5]2+ to guanine. Direct binding modes.

J. Phys. Chem. B, Vol. 103, No. 51, 1999 11419

Gresh and Sˇ poner

11420 J. Phys. Chem. B, Vol. 103, No. 51, 1999 acceptor with respect to two inner-shell waters W3 and W4 instead of only W4. This results in an about 6 kcal/mol energy gain, comparable to the ab initio computations. Finally, in g Zn2+ is bidentate bound to both O6 and N7. It retains only three inner-shell water molecules, two of which act as H-bond donors to O6, whereas the third acts as an H-bond donor to both outer-shell water molecules. The SIBFA computations give g as the least favorably bound complex (∆E ) -348.6 kcal/mol), whereas the CEP 4-31G+(2d) computations give a binding energy of -358.6 kcal/mol and a stability similar to that of c so that both c and g rank last on the scale. The most important numerical discrepancy with respect to the ab initio computation occurs with g. The relative error is of 3%, but can be anticipated to be reduced upon taking into account the BSSE correction which in the range of 3-6 kcal/mol out of 350, as observed for complexes a, g, and h. The energy difference between the best and the least-bound complexes is 17.1 kcal/ mol out of 350 kcal/mol. The ab initio HF computations give rise to smaller differences in the binding energies, of

d 4.1 4-fold O6bound

>

b 5.8 6-fold O6bound

)

c 5.8 5-fold N 7bound

>

e 9.5 4-fold N7bound

>

g 17.5 5-fold O6/N7

The first five complexes are characterized by significant variations in Zn2+ binding modes, as this cation adopts coordination numbers ranging from 4 to 6, and binds to either N7, or O6. They nevertheless differ by only 5.8 kcal/mol out of 400 kcal/mol, i.e., relative difference of less than 1.5%. This highlights the considerable versatility of Zn2+ binding modes. Completion of Zn2+ hydration shell by a sixth water molecule, as well as inclusion of bulk solvation effects, can be anticipated to further enhance such a variability. A somewhat altered ordering, favoring Zn2+ binding to O6 rather than to N7, is found in the absence of Edisp. This reflects the contribution of this

term into favoring the “softer” N site with respect to the “harder" O one for Zn2+ binding.12d The numerical agreement between ∆Etot and ∆EMP2 with the CEP 4-31G+(2d) basis set remains at

a

0.5 2.9

13.9 12.2

d 26 26.8

A similar ordering is also given by the energy-minimized 6-31G* structures: c δ∆E(SCF/6-31G*)

0.0

g

b

>

a 8.6

>

d 17.5

The ordering found in SIBFA is dictated by the E1 term. We note that among modes a-c, a with Zn2+ bound to N7 having the smallest E1 values, also has the largest E2 ones. In fact, the largest values of both Epol and Ect are those computed for through-water binding mode d. The Zn2+-induced distortions of the G-C base pair as in b and c, should be buttressed in actual double helices due to stacking with the neighboring base pairs, and the influence of the backbone. This will be examined in a forthcoming work.

Complexes of Pentahydrated Zn2+

J. Phys. Chem. B, Vol. 103, No. 51, 1999 11425

Figure 5. Binding of [Zn(H2O)5]2+ to the G-C base pair.

Analysis of nonadditivity effects was carried out in two ways. First, [Zn(H2O)5]2+ was considered as a single molecular entity. This “trimer” approach is simpler and has been used in the previous ab initio studies.4b,5b,15 The nonadditivity evaluated in this way provides directly the polarization enhancement (see above) of base pairing due to the cation binding. In the second approach, each solvation water and Zn2+ is considered as a separate entity. This is consistent with the fact that each water bears the same initial set of CEP 4-31G+(2d) permanent multipoles and polarizabilities: these only undergo the necessary reorientational changes following geometrical variations, and to them are added the induced dipoles specific to each. The first type of analysis shows all four complexes to retain a strong cooperative character. δEnadd is in the -7 to -11 kcal/mol range, being the smallest in the through-water binding mode. The SIBFA value is again in close agreement with the published ab initio result for complex a. Compared to the complexes with bare Zn2+ (δEnadd ) -15.4 kcal/mol), cooperativity is reduced by one-third in the direct binding modes. This reflects the reduction of the dicationic field due to the water molecules. δEnadd stems again essentially from Epol, but the relative weight of Ect has increased in modes b and c as compared to the bare Zn2+ complexes with G-C. The second type of analysis gives all complexes as being strongly anticooperative, and differentiates much more strongly between direct and through-water binding modes. Thus δEnadd amounts to 163-166 kcal/mol in a-c, and to 97 kcal/mol in

mode d. It is obviously dominated by anticooperative effects occurring within the cation coordination shell, leaving out a minor role for the involvement of cytosine. In other words, the repulsive nature of the total nonadditivity reflects the increased interligand repulsion due to polarization effects. The numerical values of δEnadd found for the direct binding mode are comparable to those recently derived upon studying the direct binding modes of [Zn(H2O)5]2+ to the formate anion which are in the range 160-190 kcal/mol.31 The onset of nonadditive effects is clearly a crucial feature, necessary to reduce the energy gap between direct and through-water Zn2+ complexes. Thus, in SIBFA the summed pairwise interactions in complex d is 95 kcal/mol smaller than in the best-bound complex a, while this difference is reduced to 27 kcal/mol when nonadditivity is taken into account. It will be very instructive, upon extending our approach to actual models of DNA and RNA oligonucleotides, to monitor the extent of additional modulation of this energy separation due to the oligonucleotide field on the one hand and to environmental factors on the other hand. A preliminary estimate of the latter type of effects was afforded by doing a single-point energy computation of the solvation effects of the four complexes using the Langlet-Claverie Continuum procedure.29 Thus the values of ∆Hsolv are -215.4, -217.1, -214.1, and -226.4 kcal/mol in modes a-d, respectively, preferentially favoring the through-water binding mode. This further reduces the energy gap between direct and through-water binding modes, from 27 to 13 kcal/mol out of the total of >650 kcal/mol.

11426 J. Phys. Chem. B, Vol. 103, No. 51, 1999

Figure 6. Binding of [Zn(H2O)5]2+ to the A-T base pair.

(b) [Zn(H2O)5]2+ Binding to A-T and A-A (Figure 6a-b). Similar to the bare cation case, the binding of pentahydrated Zn2+ to A-T and A-A base pairs does not produce any significant distortion of the base-pairing pattern. With [Zn(H2O)5]2+ considered as one single entity, nonadditivity is virtually null with the A-T base pair, whereas a weak cooperative effect (-1.4 kcal/mol) can be seen with the A-A base pair. When all eight interacting partners are considered each as a separate entity, a strong anticooperativity can be evidenced. δEnadd amounts to 174 kcal/mol, a value somewhat larger than in the complexes of [Zn(H2O)5]2+ with the G-C pair. Within the “trimer” approach, the SIBFA calculations have a close agreement with the ab initio nonadditivities, though the ab initio procedure shows null nonaddivity for AA pair and weak (-1.5 kcal/mol) for the AT base pair.5b However, for all base pairs studied SIBFA qualitatively reproduces the magnitude of the polarization enhancement of base pairing including the stark difference between adenine-containing and guanine-containing base pairs. Conclusions The present computations highlight the significant amount of variability in the modes of binding of pentahydrated Zn2+ to guanine and adenine. Thus, for guanine binding we have characterized structures in which the coordination number varied between 4 and 6, and N7 and O6 could be interchanged as Zn2+ ligands. Nevertheless, such variations resulted in changes of only 1-2% in the binding energies. A similar versatility in the through-water binding modes was also observed, with structures having Zn2+ bound to guanine through 2, 3, and 4 water molecules being very close energetically. In the prospect of large-scale computations on divalent cationbound oligonucleotides, an essential objective of the present

Gresh and Sˇ poner work was to benchmark the SIBFA polarizable molecular mechanics procedure and its ability to predict new structures. We found that the values of the binding energies ∆E(SIBFA) and ∆Etot(SIBFA) reproduce with an accuracy of 2-3% the corresponding HF and MP2 values using the CEP 4-31G+(2d) basis set. This basis set provides binding energies intermediate between those of the 6-311G** and 6-311G+(2d,2p) sets. ∆E(DFT) have a close numerical agreement with the corresponding ∆E(MP2) values, although in some cases minor inversions in the ordering of the bound complexes could be observed. HF/6-31G* energy minimizations allowing for the monomer relaxation reduce the energy differences between the sets of direct and through-water Zn2+ complexes, but yield in the general case a closely similar structure as the initial one. A match between ab initio and molecular mechanics results is meaningful and reproducible only if the individual components of the binding energies, namely E1 in first-order and Epol and Ect within E2 in second-order, are themselves correctly computed. This was borne out by energy decomposition of ∆E(HF) using the RVS procedure. It was observed that the individual components of the ab initio and SIBFA binding energies had much larger variations in the different complexes (of up to 60 kcal/mol for E1). It is in marked contrast with the total energies, whose amplitudes of variations were contained within 10 kcal/ mol out 400 kcal/mol. Another striking illustration of the independent character of E1 and E2 was given by the comparisons between inner- and outer-shell energies of the Zn2+ complexes with guanine. Thus E1 favored the inner-shell complexes, whereas E2 favored the outer-shell ones, in spite of the increased separation between Zn2+ and guanine. An important feature of the Zn2+ complexes with G-C base pair as contrasted to the A-T base pair resides in their substantial cooperativity, as previously demonstrated by the ab initio calculations of Burda et al.4a It amounts to -15 kcal/mol which significantly strengthens the three standard G-C hydrogen bonds, whereas it is only -2 kcal/mol for the A-T base pair. The SIBFA computations were able to account faithfully for these contrasted cooperative characters and resulted in numerical values of δEnadd fully consistent with the ab initio ones. Such computations furthermore predicted a novel distorted arrangement of the Zn-complexed G-C base pair to be stabilized. In it the O6(G)-N4(C) and N1(G)-N3(C) hydrogen bonds are disrupted, whereas O2(C) acts as a simultaneous H-bond acceptor from both N2(G) and N1(G). These features were confirmed by the SCF computations with both CEP 4-31G+(2d) and 6-31G* basis sets. The binding of pentahydrated Zn2+ to the G-C, A-T, and A-A base pairs was also investigated. This confirmed that several competing structures with the G-C base pair should have closely comparable energies, of which some involve a disruption of the N2-O2 and N1-N3 hydrogen bonds, and/or the insertion of a water molecule between G and C. Accounting for solvation effects by means of a Continuum reaction field procedure reduced the energy gap between direct and through-water Zn binding to guanine to 13 kcal/mol out of 650. The present computations showed the SIBFA procedure to be adequate to investigate the binding of Zn2+ to larger-sized oligonucleotides using representative starting conformations. In nucleic acids, the perturbations observed on the isolated base pairs should be buttressed due to stacking interactions with the neighboring base pairs, and to the limited accessibility of O6/ N7(G) or N6/N7(A) in the major groove. It will be highly instructive in this respect to compare the differential responses

Complexes of Pentahydrated Zn2+ of deoxy- versus ribooligonucleotides to the effects of Zn2+ and hydration. This will be the next step of our computations. Acknowledgment. The CEP 4-31G+(2d) Hartree-Fock and MP2 computations, and the DFT computations were performed on the computers of the Institut du De´veloppement et des Ressources en Informatique Scientifique (CNRS, Orsay, France) under project no. 990220, and of the Centre de Ressources Informatiques de Haute-Normandie (Mont-Saint-Aignan, France) under project no. 1998053. The CEP-4-31G+(2d) computations of the complexes of Zn2+ with the G-C base pair were performed thanks to Dr. Claude Giessner-Prettre (Laboratoire de Chimie The´orique, UMR 716 CNRS Universite´ Pierre-etMarie-Curie, Paris), whom we gratefully acknowledge. J.S. was supported by the following grants: A4040903 from IGA AS CR, VW-Stiftung I/74657 and 203/97/0029 from the GA CR. We also thank the Supercomputer Center Brno for a generous allotment of computer time. References and Notes (1) Klevickis, C.; Grisham, C. M. In Metal ions in Biological Systems; Sigel, A., Sigel, H., Eds.; Marcell Dekker: New York, 1996; Vol. 32, p 1. (2) Potaman, V. N.; Soyfer, V. N. J. Biomol. Struct. Dyn. 1994, 11, 1035 (3) Sabat, M.; Lippert, B. In Metal ions in Biological Systems; Sigel, A., Sigel, H., Eds.; Marcel Dekker: New York, 1996; Vol. 33, p 143. (4) (a) Burda, J. V.; Sponer, J.; Leszczynski, J.; Hobza, P. J. Phys. Chem. B 1997, 101, 9670. (b) Sponer, J.; Burda, J. V.; Sabat, M.; Leszczynski, J.; Hobza, P. J. Phys. Chem. A 1998, 102, 5951. (5) (a) Sponer, J.; Sabat, M.; Burda, J. V.; Doody, A. M.; Leszczynski, J.; Hobza, P. J.; Hobza, P. J. Biomol. Struct. Dyn. 1998, 16, 139. (b) Sponer, J.; Sabat, M.; Burda, J. V.; Leszczynski, J.; Hobza, P. J. Phys. Chem. B 1999, 103, 2528. (6) Bock, C. W.; Katz, K. A.; Glusker, J. P. J. Am. Chem. Soc. 1995, 117, 3754. (7) Beltran, R.; Martinez-Balbas A.; Bernues, J.; Bowater, R.; Azorin, F. J. Mol. Biol. 1993, 230, 966. (8) Matrinez-Balbas, A.; Azorin, F. Nucl. Acids Res. 1993, 21, 2557. (9) Ortiz-Lombardia M.; Eritja, R.; Azorin, F.; Kypr, J. Tejralova, I.; Vorlickova, M. Biochemistry 1995, 34, 14408. (10) Laundon, C.; Griffith, J. D. Biochemistry 1987, 26, 3759. (11) Palecek, E.; Brazdova, M.; Cernocka, H.; Vlk, D.; Brazda, V.; Vojtesek, B. Oncogene 1999, 18, 3617. (12) (a) Gresh, N.; Claverie, P.; Pullman, A. Int. J. Quantum Chem. 1986, 29, 101. (b) Gresh, N.; Leboeuf, M.; Salahub, D. R. In Modeling the Hydrogen Bond; ACS Symp. Ser.; Smith, D. A., Ed.; American Chemical Society: Washington, DC, 1994; No. 569, p 82. (c) Gresh, N. J. Comput. Chem. 1995, 16, 856. (d) Gresh N.; Garmer, D. R. J. Comput. Chem. 1996, 17, 1481. (e) Gresh, N. J. Chim.-Phys. 1997, 94, 1365. (f) Gresh, N. J. Phys. Chem. 1997, 101, 8680. (g) Garmer, D. R.; Gresh, N.; Roques B.-P.

J. Phys. Chem. B, Vol. 103, No. 51, 1999 11427 Proteins: Struct., Funct., Genet. 1998, 31, 42. (h) Gresh, N.; Parisel, O.; Giessner-Prettre, C. J. Mol. Struct. (THEOCHEM) 1999, 458, 27. (13) Spackova, N.; Berger, I.; Sponer, J. J. Am. Chem. Soc. 1999, 121, 5519. (14) Stevens, W. J.; Basch, H.; Krauss, M. J. Chem. Phys. 1984, 81, 6026 (15) Sponer, J.; Burda, J. V.; Leszczynski, J.; Hobza, P. J. Biomol. Struct. Dyn. 1999, 17, 61. (16) Pople, J. A.; Binkley, J. S.; Seeger, Int. J. Quantum Chem., Symp. 1976, 10, 1. (17) Gaussian 94, Revision E.2; Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Gill, P. M. W.; Johnson, B. G.; Robb, M. A.; Cheeseman, J. R.; Keith, T.; Petersson, G. A.; Montgomery, J. A.; Raghavachari, K.; Al-Laham, M. A.; Zakrzewski, V. G.; Ortiz, J. V.; Foresman, J. B.; Cioslowski, J.; Stefanov, B. B.; Nanayakkara, A.; Challacombe, M.; Peng, C. Y.; Ayala, P. Y.; Chen, W.; Wong, M. W.; Andres, J. L.; Replogle, E. S.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Binkley, J. S.; Defrees, D. J.; Baker, J.; Stewart, J. P.; Head-Gordon, M.; Gonzalez, C.; Pople, J. A. Gaussian, Inc.: Pittsburgh, PA, 1995. (18) Stevens, W. J.; Fink, W. H. Chem. Phys. Lett. 1987, 139, 15. (19) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; J. A. Montgomery Jr., J. A. J. Comput. Chem. 1993, 14, 1347. (20) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 533. (21) Sokalski, W. A.; Roszak, F.; Hariharan, P. C.; Kaufman J. J. Int. J. Quantum Chem. 1983, 23, 847. (b) Cammi, R.; Bonnacorsi, R.; Tomasi, J. Theor. Chim. Acta 1985, 68, 271. (22) (a) Lee, C., Yang, W., Parr, R. G. Phys. ReV. 1988, B37, 785. (b) Becke, A. D. J. Chem. Phys. 1993, 98, 1372. (23) Vigne´-Maeder, F.; Claverie, P. J. Chem. Phys. 1988, 88, 4934. (24) Gresh, N.; Claverie, P.; Pullman, A. Theor. Chim. Acta 1984, 20, 1. (25) Garmer, D. R.; Stevens, W. J. J. Phys. Chem. 1989, 93, 8263. (26) Arnott, S. R.; Chandrasekaran, D.; Birdsall, A.; Leslie, A.; Ratliff, R. Nature 1980, 283, 743. (27) Gresh, N.; Pullman, B.; Arcamone, F.; Menozzi, M.; Tonani, R. Mol. Pharmacol. 1990, 35, 251. (28) Evangelakis, G.; Rizos, J.; Lagaris, I.; Demetropoulos, G. N. Comput. Phys. Commun. 1987, 46, 401. (29) Langlet, J.; Claverie, P.; Caillet, J.; Pullman, A. J. Phys. Chem. 1988, 92, 1631. (30) (a) Langlet, J.; Gresh, N.; Giessner-Prettre, C. Biopolymers 1995, 36, 765. (b) Gresh, N.; Roques, B.-P. Biopolymers 1997, 41, 145. (31) Tiraboschi, G.; Roques, B.-P.; Gresh, N. J. Comput. Chem. 1999, 20, 1379. (32) Masella, M.; Gresh, N.; Flament, J.-P. J. Chem. Soc., Faraday 1998, 94, 2745. (33) (a) Sponer, J.; Hobza, P. J. Phys. Chem. 1994, 98, 3161. (b) Sponer, J.; Hobza, P. J. Am. Chem. Soc. 1994, 116, 716. (c) Bludsky, O.; Sponer, J.; Leszczynski, J.; Spirko, V.; Hobza, P., J. Chem. Phys. 1997, 105, 11042. (d) Sponer, J.; Florian, J.; Leszczynski, J.; Hobza, P. J. Biomol. Struct. Dyn. 1996, 13, 827. (e) Sponer, J.; Leszczynski, J.; Hobza, P., J. Biomol. Struct. Dyn 1996, 14, 117. (f) Luisi, B.; Orozco, M.; Sponer, J.; Luque, F. J.; Shakked, Z. J. Mol. Biol. 1998, 279, 1123.