Calculations of Magnesium−Nucleic Acid Site Binding in Solution

Juneau, K.; Podell, E.; Harrington, D. J.; Cech, T. R. Structure 2001, 9, 221. [Crossref] ...... Anton S. Petrov , Chad R. Bernier , Chiaolong Hsiao ,...
1 downloads 0 Views 236KB Size
Download PDF
6072

J. Phys. Chem. B 2004, 108, 6072-6081

Calculations of Magnesium-Nucleic Acid Site Binding in Solution Anton S. Petrov,* George R. Pack,* and Gene Lamm Department of Chemistry, UniVersity of LouisVille, LouisVille, Kentucky 40292 ReceiVed: NoVember 18, 2003; In Final Form: February 27, 2004

The interaction of nucleic acids with metal cations, particularly magnesium, plays an important role in stabilizing their tertiary structure. An accurate quantitative description of these interactions requires high-level quantum chemical calculations. Because metal-nucleic acid interactions occur in water, solvation effects are important. To explore the influence of solvent we calculated the interaction energies of guanine and the dimethyl phosphate anion with a hexahydrated magnesium ion, which represent the main binding modes of magnesium and nucleic acids. We used density functional theory within the CPCM model and compared these results with data obtained for these complexes using the GBSA and PB approaches with the AMBER ff99 force field. The comparison indicates that the empirical force field significantly underestimates the total interaction energy and that the main source of error is an inaccurate description of the gas-phase interaction term. The presence of charge transfer and polarization contributions in the latter leads to a contraction of the complex geometry relative to that predicted by the classical potential. This lack increases the error in the nuclear-centeredcharge approximation of the electrostatic energy. A method to calculate total interaction energies applicable for large systems in solution is proposed.

1. Introduction Macromolecular structure, which is a determinant of functional properties, is stabilized by an array of small species including water molecules and metal ions. For example, the binding of divalent metal ions not only results in changes in the local conformations of proteins and nucleic acids, but also determines their tertiary structures1 and can play an important role in catalysis.2 Because of this, the prediction of the binding affinities of these ions should be helpful in understanding macromolecular folding and in explaining structure-function relationships in these systems. In this work, we address the specific binding of magnesium to nucleic acids by using conventional classical and quantum mechanical approaches to calculate the total interaction energy in solution. A divalent cation may coordinate directly to a macromolecular binding site or indirectly through water molecules comprising the solvation shell of the ion. Site-specific binding of magnesium to nucleotide bases, particularly guanine, has been observed in certain crystal structures.1,3,4 Mobile counterions within a few angstroms from the macromolecular surface are strongly influenced by the electrostatic field of the polyelectrolyte and form a layer of “condensed” ions.5 For site-specific binding, the total binding affinity depends on the balance of several factors: the electrostatic interaction between opposite charges, nonelectrostatic forces involved in the formation of H bonds, the entropy of binding or displacing water, electrostatic screening by other ions in the system, and changes in the polarizability (the local dielectric coefficient) of the environment. Computational simulations of macromolecular systems generally rely on molecular dynamics (MD)6,7 or Monte Carlo (MC)8 methods based on empirical pairwise-additive force fields. Methods combining the generalized Born approach with Poisson-Boltzmann theory9 have been used to estimate the binding energy between cobalt(III) hexammine and RNA.10 * Address correspondence to these authors. E-mail: gwise.louisville.edu. E-mail: [email protected].

aspetr01@

However, a major drawback of these methods is a failure to account for the nonelectrostatic charge transfer and polarization contributions that stabilize the complex.11,12 This often leads to total interaction energies that underestimate the actual values10 as well as to a displacement of the metal cation position compared to crystal structure data.13 The importance of magnesium in determining nucleic acid structure and the fact that its interaction with guanine and phosphate are not adequately described by electrostatics alone led us to consider these systems in detail. A large number of studies on relatively small systems that model nucleic acid binding sites have been performed by using correlated ab initio methods. The results of these investigations have been extensively reviewed.14,15 Most of the calculations for these systems have been performed in a gas-phase environment.11,16-19 Although these studies provided much helpful information about competition between different sitespecific binding modes (inner- vs outer-sphere coordination) revealing the degree of nonadditivity of many-body interactions and the strength of nonelectrostatic interactions between monomers, the effect of solvent on interactions between biomolecules is often neglected. Recent advances in the development of quantum continuumsolvent models, e.g., the polarized continuum model (PCM)20 and its modifications,21-23 and the increased power of computational facilities have allowed the treatment of small systems of nucleotide base-metal complexes in which the presence of solvent was taken into account.12,24 These models have successfully been used to investigate the stability of charged25,26 and neutral27,28 H-bonded complexes as well as to estimate the hydration energies of various metal cations.29 There have also been a few publications describing first-principles-based calculations of the effective electrostatic potential of hydrated cations in solution.30,31 A major conclusion of these studies is that the accurate treatment of binding interactions in a condensed phase demands that solvent effects be included.

10.1021/jp037517s CCC: $27.50 © 2004 American Chemical Society Published on Web 04/17/2004

Calculations of Mg-Nucleic Acid Site Binding in Solution In the present study we focus on the interaction between a hexahydrated magnesium ion and its two primary specific binding sites in DNA and RNA molecules: the nucleotide base guanine (G) and the charged backbone, modeled as a tg-dimethyl phosphate anion (DMP-). The goal of this paper is to present an accurate but computationally inexpensive approach to estimate the total interaction energy that can be extended to systems of significantly larger size. To achieve this goal, the results of extensive quantum chemical calculations were compared with simpler techniques employing empirical analytical potentials. This revealed a possible source of error in the latter methods pointing out that solvent-influenced binding interactions in macromolecular systems should be treated with a higher degree of accuracy than is normally considered. This may be important in the prediction of the tertiary structure of nucleic acids, especially RNAs, which can be highly stabilized by magnesium ions.1 2. Computational Details Gas-phase calculations were performed by using density functional theory (DFT) at the B3LYP level, which combines the GGA exchange three-parameter hybrid functional developed by Becke32 and the correlation functional of Lee-Yang-Parr,33 as implemented in the Gaussian 98 suite of programs.34 Geometries of complexes and monomers were obtained by using the double-ζ split valence 6-31G(d,p) basis set with one set of polarization functions on all atoms. Single-point energy calculations were done with the SCF ) tight option at the same level of theory. All interaction energies were calculated by treating the complexes as G-Mg(H2O)62+ and DMP--Mg(H2O)62+ dimers. The basis set superposition error (BSSE) in the dimercentered basis set was obtained by applying the counterpoise procedure of Boys and Bernardi.35 The deformation of monomers was taken into account in the majority of the calculations according to a scheme proposed by van Duijneveldt-van de Rijdt and van Duijneveldt,36 unless otherwise stated. The interaction energy was partitioned by using the scheme proposed by Kitaura and Morokuma37,38 at the HF/6-31G(d) level, using the B3LYP/ 6-31G(d,p) geometries of the gas phase-optimized complexes; this decomposition was performed with the GAMESS package.39 To describe the solvent continuum we used the COSMOPCM (CPCM) approach,23,40 which does not require an expensive coupled iterative procedure within an SCF cycle. According to this model the polarization of a solvent with dielectric coefficient  induces a solute charge distribution described by a perturbed Hamiltonian:

H ˆ )H ˆ 0 + Vˆ

(1)

where H ˆ 0 is the Hamiltonian of the isolated solute in a vacuum and the operator Vˆ describes solute-solvent interactions. A perturbation Vˆ induced by a solvent modifies an eigenfunction Ψ0 of the Hamiltonian H ˆ 0 resulting in an eigenfunction Ψ. This yields the energy of a species in the gas phase

ˆ 0 |Ψ0〉 E)1 ) 〈Ψ0|H

(2)

 ˆ 0 + Vˆ /2|Ψ〉 + En.e. E ) 〈Ψ|H

(3)

and in solution

where the second term in eq 3 is the nonelectrostatic contribution responsible for the formation of the solvent-inaccessible cavity. The molecular cavity was formed by interlocking spheres

J. Phys. Chem. B, Vol. 108, No. 19, 2004 6073 centered on solute atoms or atomic groups with radii taken from the United Atom topological model41 and scaled by a factor of 1.2. Several values for the dielectric coefficient were used, ranging from that of bulk water (w ) 78.4) to lower values simulating less polar solvents (in the range 5-40). All calculations based on empirical potentials were carried out with the AMBER 7.0 suite of programs,42 using the parameters of the ff99 force field,43 most of which have been originally derived for the ff94 of Cornell et al.44 The only parameter changes from the ff99 force field were the RESP charges for the hexahydrated magnesium ion and for the guanine and DMP monomers, which were derived at the HF/6-31G(d) level according to the method of Kollman and co-workers.44,45 The necessary RESP charges for the monomers in the study were derived at the HF/6-31G(d) level. The solvent was taken into account by applying the generalized Born surface area approach (GBSA)46,47 with the modified set of Bondi radii as implemented in the Sander module of AMBER. The energy of a species in solution was obtained as

E ) INT + ES + VDW + GBSA

(4)

where INT is an energy term that accounts for the internal degrees of freedom (e.g., bond bending and stretching), ES is electrostatic interaction energy calculated by Coulomb’s law and with the set of nuclear-centered RESP charges, VDW is the van der Waals interaction energy obtained with the LennardJones potential, and GBSA is the contribution that accounts for the presence of a solvent within the generalized Born surface area approach and includes electrostatic and nonelectrostatic terms. For GBSA geometries, 1000 cycles of minimization were performed by applying the steepest decent algorithm for the first 10 cycles followed by the conjugate gradient procedure. The Poisson-Boltzmann (PB) calculations were performed with the finite difference method as implemented in the UHBD program.48 The potential at the boundary of a coarse grid (60 × 60 × 60) with grid spacing 1.5 Å was set to zero and a focusing procedure was performed at the center of mass by using a finer grid (0.3 Å) with the same number of grid points. Dielectric coefficients of 1 and 78.4 were used to represent the solute interior and exterior, respectively. A magnesium ion radius of 0.86 Å and the set of PARSE radii49 to describe the nucleic acid moieties and water molecules as well as the set of nuclear centered RESP charges were used in the PB calculations. To provide all computational details in a simple manner we modify the standard quantum mechanical notation of describing single-point energy and geometry calculations by adding in brackets the method of calculating the solvation energy (if applicable). For example, ff99[GBSA]//B3LYP/6-31G(d,p) stands for the calculation of a single-point energy performed with the ff99 force-field and the GBSA approach with a gasphase optimized geometry obtained at the B3LYP/6-31G(d,p) level. 3. Results and Discussion 3.1. Direct Calculation of the Total Interaction Energy. The total interaction energy (IE) between monomers A and B that form complex AB in a solvent with dielectric coefficient  can be calculated as

- E,def IEAB ) EAB - EA - EB - E,def A B

(5)

where Ei are the energies of complex AB and its monomers

6074 J. Phys. Chem. B, Vol. 108, No. 19, 2004

Petrov et al.

Figure 1. Equilibrium geometries of (a) G-Mg(H2O)62+ and (b) DMP--Mg(H2O)62+ obtained (I) by applying 1000 steps of minimization within the GBSA model of AMBER with the ff99 force field, (II) within the framework of the CPCM model with  ) 78.4 at the same level of theory, and (III) with the optimization in the gas phase at the B3LYP/6-31G(d,p) level. Hydrogen bond distances are shown in angstroms.

and E,def account for the deformation of monomers in the i complex compared to their isolated states. The absolute binding free energy is found by adding zero-point and entropic contributions to the total interaction energy. These energies are often difficult to obtain accurately but some work in this area has been done.10,50 IE values obtained within the framework of the CPCM model were calculated from eqs 3 and 5 and corrected for BSSE obtained from gas-phase calculations.51 Similarly, the force field-based IE values were calculated by combining eqs 4 and 5. 3.1.1. Classical Calculations. The specific interaction of magnesium with proteins and nucleic acids has proven to be very difficult to model with classical force fields. As a typical example, we present calculations using the ff99 force field in AMBER. The interaction energies calculated for G-Mg(H2O)62+ and DMP--Mg(H2O)62+ at the ff99[GBSA] level with dielectric coefficient w were found to be 2.4 and -0.5 kcal/mol, respectively. These complexes are shown in Figure 1 (structures Ia and Ib, respectively). Similar calculations performed on these complexes in the gas phase at the ff99//ff99[GBSA] level resulted in interaction energies of -32.1 and -160.4 kcal/mol. These results indicate that the presence of a high dielectric solvent strongly destabilizes the interaction between these monomers, yielding IE values near zero. Such a destabilization effect has been found for both charged and uncharged systems.27 There are a couple of simple explanations of this dielectric phenomenon. First, on the basis of electrostatic Born theory,52 the solvation of the separate Mg(H2O)62+ and guanine species is more favorable than solvation of the complex, which retains

the divalent charge but has a larger radius. This effect is more pronounced for the DMP--Mg(H2O)62+ system in which the complex has not only an enlarged cavity radius but also a decreased net charge in comparison with the separated components. Second, the decrease in stability of these complexes in solution is also due to the nonelectrostatic contribution of cavity creation. The formation of the complex leads to a reduction in the solvent-accessible surface area.28 When the complex is formed, polar groups of the monomers [G, DMP-, and Mg(H2O)62+] become solvent inaccessible, thus decreasing the free energy of solvation. 3.1.2. Quantum Calculations. To test the accuracy of classical results we also performed calculations using the quantum CPCM approach at the B3LYP/6-31G(d,p)[CPCM] level for G-Mg(H2O)62+ and DMP--Mg(H2O)62+ (structures IIa and IIb in Figure 1). These complexes are stabilized by formation of hydrogen bonds in which electrostatic interactions are augmented by the nonelectrostatic charge transfer and polarization contributions.53 Total interaction energies as a function of solvent dielectric coefficient are shown in Figure 2 (open circles). Although the absolute value of the interaction energy for both complexes decreases with increasing dielectric coefficient, this change is not simply proportional to 1/, as would be expected from a calculation of the interaction energy of two point charges immersed in a dielectric continuum. By performing CPCM calculations we found that extrapolation of the interaction energy to a solvent with infinite dielectric coefficient yields nonzero values: -19.2 kcal/mol and -27.3 kcal/mol for G-Mg(H2O)62+ and DMP--Mg(H2O)62+, respectively, as shown

Calculations of Mg-Nucleic Acid Site Binding in Solution

J. Phys. Chem. B, Vol. 108, No. 19, 2004 6075

Figure 2. Gas-phase energy (b) and solvation free energy (9) contributions to the total interaction energy (O) calculated for (a) G-Mg(H2O)62+ and (b) DMP--Mg(H2O)62+ as a function of solvent dielectric coefficient. Calculations were performed at the B3LYP/6-31G(d,p)[CPCM] level.

TABLE 1: Geometrical Parameters for the Complexes Shown in Figure 1a distance [Å] Mg-X Mg-O a

Ia

IIa

IIIa

Ib

IIb

IIIb

5.643 2.04-2.10

4.732 2.05-2.09

4.734 2.06-2.11

4.996 2.04-2.10

4.039 2.06-2.08

4.034 2.07-2.13

X denotes the C5 atom of guanine or the P atom of DMP-. A description of these complexes is given in the caption to Figure 1.

in Figure 2. Although individual values for each component of the IE cannot be determined within CPCM calculation, previous studies have shown that solvent dielectric effects reduce the electrostatic component of the interaction energy but leave the nonelectrostatic terms substantially unaffected.54,55 Therefore, the extrapolated IEs are expected to be primarily nonelectrostatic and are seen to contribute significantly to the total interaction energy. Although the value of the dielectric constant at the surface of nucleic acids lies in the range 30-50,56 in the rest of this study we shall consider two limiting cases [ ) 1 (gas phase) and 78.4 (bulk water)]. We assume that our main conclusions are still valid for intermediate values of the dielectric constant, since we have shown (Figure 2) that the total value of the interaction energy is almost independent of dielectric constant over the range 40-80. Below we investigate in detail why the classical approach underestimates the IE values and try to overcome the main disadvantage of the quantum CPCM approach when applied to large systems. 3.1.3. Gas-Phase Geometry Approximation. Optimization performed in the CPCM model is considerably more computationally expensive than that performed in the gas phase for the medium-sized complexes considered here and is currently impractical for large systems. Because of this we explored the effect of avoiding the optimization step by comparing the interaction energy with and without geometry optimization of the complex.57 The interaction energy for G-Mg(H2O)62+ and DMP--Mg(H2O)62+ at the B3LYP/6-31G(d,p)[CPCM]//B3LYP/ 6-31G(d,p) level were found to be -16.9 and -22.7 kcal/mol, respectively. The geometries of these complexes optimized in the gas phase at the B3LYP/6-31G(d,p) level are shown in Figure 1 (structures IIIa and IIIb). This figure also gives the hydrogen bond distances at the optimized geometries; additional parameters for these complexes are given in Table 1. By comparing the distance data in Table 1, it can be seen that the presence of the reaction field affects the nuclear positions of G-Mg(H2O)62+ (columns IIa and IIIa) and DMP-Mg(H2O)62+ (columns IIb and IIIb). Although the hydrogenbonding pattern remained the same, the presence of solvent leads to an increase in the calculated hydrogen bond distances for both complexes compared to their gas-phase analogues. A

Figure 3. Thermodynamic cycle used to calculate the total interaction energy IEAB between species A and B in solution. (a) The cycle consists of the transfer of interacting monomers A and B from the solvent into the gas phase (with free energy ∆Gdesol), calculation of the gas-phase interaction energy between these monomers (IEAB) 1), and solvation of complex AB (with free energy ∆Gsol). (b) Computational methods used to calculate the energy of each step.

similar trend has also been seen for other H-bonded systems formed by charged monomers.28 The lengthening of the H bonds stems from the attenuation of the electrostatic attraction. This is consistent with the observation that the elongation is significantly larger for the DMP--Mg(H2O)62+ “zwitterion” complex (0.14-0.15 Å) than for G-Mg(H2O)62+ (0.01-0.04 Å) due to the stronger electrostatic interactions in the former system. Optimization performed in the dielectric continuum (w ) 78.4) induces a slight shortening of the Mg-O distances in both complexes: these values were found to be in the range 2.052.09 Å for the G-Mg(H2O)62+ complex and 2.06-2.08 Å for DMP--Mg(H2O)62+, while their gas-phase analogues spanned 2.06-2.11 and 2.07-2.13 Å, respectively. A comparison of the geometries obtained by quantum methods with those found by the ff99[GBSA] approach (structures Ia and Ib in Figure 1) reveals that the classical treatment leads not only to elongation of the H bonds but also to a change in the H-bonding patterns, resulting in a significant increase in the Mg-C5 and Mg-P distances, as compared to values obtained by ab initio methods (Table 1). These calculations suggest that the hydrogen bonds in this complex are significantly different from the ones for which this force field was parametrized. 3.2. Thermodynamic Cycle Approach. CPCM calculations of the total interaction energy are quite computationally expensive and, furthermore, are limited to complexes of relatively small size. This difficulty can be worked around by applying a thermodynamic cycle approach (Figure 3),58 accord-

6076 J. Phys. Chem. B, Vol. 108, No. 19, 2004

Petrov et al.

Figure 4. Final positions of Mg(H2O)62+ in (a) G-Mg(H2O)62+ and (b) DMP--Mg(H2O)62+ optimized at the B3LYP/6-31G(d,p) level. Arrows indicate the total displacement of the magnesium ion from its original equilibrium position in each complex.

ing to which eq 5 can be rewritten in the form desol IEAB ) IE)1 + ∆Gdesol + ∆Gsol AB + ∆GA B AB

(6)

where the first term represents a gas-phase (or nonsolvent) is the free energy required to transfer a interaction, ∆Gdesol i monomer i from solvent to the gas phase, and ∆Gsol AB is the free energy associated with the transfer of complex AB from the gas phase to solvent. Combining the free energy terms gives

IEAB ) IE)1 AB + ∆∆Gsol

(7)

The conceptual advantage of the thermodynamic cycle approach is that it allows the total IE to be explicitly expressed as a sum of two terms, a gas-phase interaction energy and a solvation free energy difference. A computational advantage is that by choosing different methods for each term, for example, quantum versus classical, one may reduce the cost of the simulations and be able to treat larger systems. However, if different methods are used for the interaction and solvation energies, then use of a single geometry may not be consistent with the method used to calculate an energy. We have therefore investigated the effect of different methods for calculating the energy as well as the choice of geometry. 3.2.1. Gas-Phase Interaction Energy. Because the total interaction energy in a solvent can be represented according to eq 7, we decomposed the IE obtained at the B3LYP/631G(d,p)[CPCM] level into nonsolvent and solute-solvent contributions. This decomposition for G-Mg(H2O)62+ and DMP--Mg(H2O)62+ is shown in Figure 2 (filled circles and squares). Although the presence of solvent attenuates the interaction between complexes, it actually increases the stabilizing effect of nonsolvent interactions. This effect is more pronounced for G-Mg(H2O)62+ than for DMP--Mg(H2O)62+ and increases with the dielectric coefficient of the solvent. For example, the nonsolvent interactions calculated in water were found to be stronger by 11.5 kcal/mol for the former system and by 6.8 kcal/mol for the latter, compared to the gas phase. This additional stabilization stems from adjustments in the geometry of the complexes in response to the reaction field. To explore the effect of neglecting the costly CPCM optimization, we performed two sets of calculations to determine how the interaction energy for G-Mg(H2O)62+ and DMP-Mg(H2O)62+ changed as the distance between the magnesium ion and guanine or DMP- increased from the optimized value. In the first set, the calculations were performed at the B3LYP/ 6-31G(d,p)[CPCM] level, while in the second set, gas-phase geometries [at B3LYP/6-31G(d,p)] were used. For both sets all geometric coordinates were initially optimized and equilibrium positions were found without applying any constraints. To obtain the interaction energy as a function of distance between monomers, the positions of all atoms in each complex were

re-optimized but with the Mg-C5 and Mg-P distances fixed (at a prescribed value) and the angles for Mg-C5-C4 (for guanine-magnesium) and Mg-P-O (for DMP--magnesium) constrained to their equilibrium values. As the monomers are moved apart, the resulting COSMO cavity consists of two overlapping monomer cavities, each composed of interlocking spheres centered on individual monomer atoms, as described in Section 2. During the gas-phase calculations, as the magnesium ion was moved away from the binding site, water molecules originally H-bonded to guanine or to DMP- tended to preserve these hydrogen bonds instead of remaining with the magnesium ion. Therefore, to maintain the octahedral coordination of magnesium, additional constraints were applied to keep all Mg-O distances fixed at the equilibrium value. This illustrates that the proper balance of forces acting on a complex cannot be reproduced by gas-phase simulations. A similar set of constraints was also applied during the B3LYP/6-31G(d,p)[CPCM] calculations for comparison. Unfortunately, when the magnesium ion was somewhat distant (∼2 Å) from the equilibrium position for G-Mg(H2O)62+ and DMP--Mg(H2O)62+, further optimization of the complexes in the reaction field was not possible, presumably due to the decreasing overlap between monomer cavities. This prevented consideration of these interactions for separation distances larger than 2.4 and 1.5 Å, respectively. The final optimized gas-phase structures and the displacement vector of the magnesium ion are shown in Figure 4. Figure 5 shows the nonsolvent component ( ) 1) of the interaction energy as a function of the displacement of the magnesium ion from its original position in the fully optimized G-Mg(H2O)62+ and DMP--Mg(H2O)62+ complexes. Throughout the range where the B3LYP/6-31G(d,p)[CPCM] calculations could be performed, the nonsolvent interactions calculated by using the gas-phase geometries are less stable by 10-11.5 [6.7-14.2] kcal/mol for G-Mg(H2O)62+ [DMP--Mg(H2O)62+] compared to the complex optimized in the presence of the reaction field. This demonstrates that the extra stabilization that the reaction field confers to the nonsolvent interaction extends several angstroms from the ion binding site. Because the use of gas-phase B3LYP/6-31G(d,p) geometries seemed to be a reasonable approximation in the thermodynamic cycle approach, we also tested the ability of the ff99 force field to reproduce the nonsolvent interaction energy terms based on the previously obtained geometries of the complexes and individual components at this level of theory. The results of these calculations are also shown in Figure 5. These data indicate that the ff99//B3LYP/6-31G(d,p) calculations significantly underestimate the gas-phase interaction energy at close distances. As the hydrated magnesium ions approach their equilibrium positions, there is a great difference between the classical and quantum mechanical interaction energies. The differences between the energies obtained for G-Mg(H2O)62+ and DMP--

Calculations of Mg-Nucleic Acid Site Binding in Solution

J. Phys. Chem. B, Vol. 108, No. 19, 2004 6077

Figure 5. Nonelectrostatic interactions as a function of distance from the equilibrium position of Mg in (a) G-Mg(H2O)62+ and (b) DMP-Mg(H2O)62+ calculated by the following methods: B3LYP/6-31G(d,p)[ CPCM] (2); B3LYP/6-31G(d,p)[CPCM]//B3LYP/6-31G(d,p) (9); and ff99[GBSA]//B3LYP/6-31G(d,p) (b).

Figure 6. Morokuma decomposition performed at the HF/6-31G(d) level (solid lines) and AMBER ff99 force field (dashed lines). Energy terms as a function of distance from the equilibrium position of Mg for (a) G-Mg(H2O)62+ and (b) DMP--Mg(H2O)62+ with B3LYP/6-31G(d,p) geometries: total energy (b), ES (4), CT (]), POL (O), EX (0), and MIX (g); see the text for energy term descriptions.

Mg(H2O)62+ at the equilibrium positions (R ) 0) were 29.5 and 28.0 kcal/mol, respectively. However, at distances larger than 2.5 Å the nonsolvent contributions to the interaction energy calculated with the B3LYP/6-31G(d,p) and AMBER ff99 force field methods agree. 3.2.2. Gas-Phase Energy Decomposition. To understand the approximations responsible for the large underestimation of the gas-phase interaction potential by the classical force field, we decomposed the energies calculated by each of the methods into electrostatic and nonelectrostatic terms. For the gas-phase ab initio calculations, a Morokuma decomposition analysis was performed at the HF/6-31G(d) level taking into account BSSE. This scheme partitions the gas-phase interaction energy into electrostatic (ES), polarization (POL), charge transfer (CT), exchange-correlation (EX), and mixing (MIX) contributions. However, any deformation energy is neglected in this scheme, that is, the geometry of each component is taken to be that in the complex. The energy contributions were compared with similar terms calculated according to the classical force field. To be consistent with the Morokuma decomposition scheme, we recalculated the interaction potentials of the complexes using the AMBER ff99 force field, neglecting the deformational energy of the components. Those energies consist of a Couloumb electrostatic term approximated by nuclear-centered point-charge interactions, a van der Waals energy obtained from a Lennard-Jones potential, and an internal energy term that accounts for bond stretching and bending. The results of the quantum and classical decompositions are presented in Figure 6. The first observation is that the electrostatic component is significantly underestimated by the

classical approximation. This term of the gas-phase interaction energy calculated by the AMBER ff99 force field is too small by 8.5 and 32.3 kcal/mol for G-Mg(H2O)62+ and DMP-Mg(H2O)62+, respectively, in comparison with the ab initio calculations, a relative error of about 11% for both systems. As noted above, this inconsistency is observed within 2.5 Å from the equilibrium distance. These results indicate that neglect of higher multipoles of the atomic charge distributions can lead to significant errors in the molecular electrostatic potential in the vicinity of the molecular surface.59,60 To improve the electrostatic term in the classical force field, higher order terms61,62 as well as nonnuclear-centered charges63 have been suggested and may be necessary for these complexes. A second observation stemming from the comparison in Figure 6 is that polarization, charge transfer, and exchange repulsion each make significant contributions to the interaction energy. In many hydrogen-bonded complexes, these terms, while individually significant, tend to cancel out. The importance of the CT and POL terms was found to be more pronounced for G-Mg(H2O)62+. Figure 6a shows that the magnitude of the total interaction energy is greater than the Coulomb term, revealing the stabilizing role of these terms in this system. In the case of DMP--Mg(H2O)62+, the stabilizing effect of these terms is less significant (Figure 6b) and is less in magnitude than the repulsion terms, making the electrostatic energy greater than the total stabilization. This indicates that there is a large repulsion overlap of the closed-shell electron clouds of the interacting species. Overall, the classical methodology of treating a tightly bound complex in terms of the properties of the individual monomers

6078 J. Phys. Chem. B, Vol. 108, No. 19, 2004

Petrov et al.

TABLE 2: Solvation Energies ∆Gsol and ∆∆Gsol [kcal/mol] in a Solvent with E ) 78.4 Calculated for G-Mg(H2O)62+ and DMP--Mg(H2O)62+ and Their Monomersa Ib

IIc

IIId

IVe

Vf

∆Gsol[G-Mg(H2O)62+] ∆Gsol[G] ∆Gsol[Mg(H2O)62+] ∆∆Gsol[G-Mg(H2O)62+]

-166.13 -23.29 -202.15 59.31

-158.09 -21.36 -191.01 54.28

-261.18 -26.20 -287.54 52.56

-288.09 -25.20 -298.35 35.46

-332.18 -42.53 -347.72 58.07

∆Gsol[DMP--Mg(H2O)62+] ∆Gsol[DMP-] ∆Gsol[Mg(H2O)62+] ∆∆Gsol[DMP--Mg(H2O)62+]

-69.83 -60.01 -202.15 192.33

-55.86 -58.39 -191.01 193.54

-160.60 -74.00 -287.54 200.94

-213.21 -74.80 -298.35 159.93

-255.56 -95.67 -347.72 187.83

a

The energies were calculated by using the methods (see the text for notation) in footnotes b-e. b B3LYP/6-31G(d,p)[CPCM]. c B3LYP/631G(d,p)[CPCM]//B3LYP/6-31G(d,p). d ff99[GBSA]//B3LYP/6-31G(d,p). e B3LYP/6-31G(d,p)[GBSA]//B3LYP/6-31G(d,p). f B3LYP/6-31G(d,p)[PB]// B3LYP/6-31G(d,p).

Figure 7. Solvation free energy differences as a function of distance from the equilibrium positions of Mg for (a) G-Mg(H2O)62+ and (b) DMP-Mg(H2O)62+ calculated by the following methods: B3LYP/6-31G(d,p)[CPCM] (2); B3LYP/6-31G(d,p)[CPCM]//B3LYP/6-31G(d,p) (9); and ff99[GBSA]//B3LYP/6-31G(d,p) (b).

fails for magnesium binding to DNA. The use of nuclearcentered charges and incomplete cancellation of polarization and charge-transfer effects by exchange-repulsion contributes to the inaccuracy. This inaccuracy in the calculation of gasphase energies is the principal source of error in computing the aqueous phase interaction energy between the hydrated magnesium ion and specific sites on nucleic acids. The reason for this is that the nonsolvent interaction term enters the expression of the total interaction energy without cancellation, as seen by eq 7. 3.2.3. SolVation Free Energy Differences. The second term in eq 7 accounts for the presence of solvent and tends to decrease the stability of a complex. Because of the difficulty in applying quantum calculations to large biochemical systems, this section discusses some alternative ways to calculate solvation energies of such systems. Solvation energies calculated at the B3LYP/6-31G(d,p)[CPCM] level for G-Mg(H2O)62+ and DMP--Mg(H2O)62+ and their monomers are given in Table 2 (column I). The destabilizing effect of the solvent on complex formation (i.e., the difference between the solvation energies of the complexes and the monomers) was found to be 59.3 and 192.3 kcal/mol, respectively. The geometrical rearrangement caused by solvation alters the calculated solvation energy. These energies calculated at the B3LYP/6-31G(d,p)[CPCM]//B3LYP/6-31G(d,p) level (gas-phase geometries) were a few kilocalories per mole less negative than their analogues obtained from the calculations with the CPCM-optimized structures (Table 2, column II). Values for ∆∆Gsol (the second term of eq 7) were evaluated with both approaches. These data are plotted in Figure 7 as a function of the distance from the equilibrium position of the magnesium ion. The difference between solvation free energies calculated by using B3LYP/6-31G(d,p)[CPCM] and B3LYP/6-

31G(d,p)[CPCM]//B3LYP/6-31G(d,p) was 5.0 and 1.2 kcal/mol for G-Mg(H2O)62+ and DMP--Mg(H2O)62+, respectively. As an alternative to the quantum mechanical CPCM theory we used the GBSA method as implemented in the AMBER package for the calculation of solvation free energies. Because this method is based upon empirical parameters, it is much faster than the CPCM method and allows treatment of large biologically relevant systems.10 The solvation energies of both complexes and their monomers calculated at the ff99[GBSA]// B3LYP/6-31G(d,p) level are presented in Table 2 (column III). A comparison of these energies with those obtained within the CPCM model shows them to be significantly higher, especially for DMP--Mg(H2O)62+. However, the values of ∆∆Gsol (presented in Table 2 and Figure 7) reveal only a few kilocalories per mole discrepancy with the CPCM data obtained for the same set of geometries. This clearly illustrates the cancellation of errors occurring in the calculation of the solvation energy of the complexes and the corresponding monomers and testifies to the applicability of this approach in estimating solvation free energies. Another observation drawn from Figure 7 is that ∆∆Gsol monotonically decreases as a function of the distance from the equilibrium position of the magnesium ion. Another classical method for calculating solvation free energies, Poisson-Boltzmann (PB) theory, was also tested. Data in Table 2 (columns III and V) indicate that there is also a significant discrepancy in the values of the solvation energies obtained by the GBSA and PB approaches and with the same geometry of the complexes. However, the same data reveal that ∆∆Gsol calculated by these two methods differs by only 7.5 and -13.1 kcal/mol for G-Mg(H2O)62+ and DMP--Mg(H2O)62+, respectively. At the same time we need to point out that, although the GBSA approach is less computationally expensive and is applicable to molecular dynamics simulations, this method

Calculations of Mg-Nucleic Acid Site Binding in Solution

J. Phys. Chem. B, Vol. 108, No. 19, 2004 6079

Figure 8. Total interaction energies ( ) 78.4) as a function of distance from the equilibrium positions of Mg in (a) G-Mg(H2O)62+ and (b) DMP--Mg(H2O)62+ calculated by the following methods: B3LYP/6-31G(d,p)[CPCM] (2); B3LYP/6-31G(d,p)[CPCM]//B3LYP/6-31G(d,p) (9); and ff99[GBSA]//B3LYP/6-31G(d,p) (b). Values indicated with an open circle (O) were obtained by the combined B3LYP/6-31G(d,p)[GBSA]// B3LYP/6-31G(d,p) approach.

TABLE 3: Total Interaction Energies [kcal/mol] in a Solvent with E ) 78.4 Calculated for G-Mg(H2O)62+ and DMP--Mg(H2O)62+ and Their Monomersa complex 2+

G-Mg(H2O)6 DMP--Mg(H2O)62+

Ib

IIc

IIId

IVe

Vf

VIg

2.4 -0.5

9.8 16.4

-20.2 -30.3

-16.9 -22.8

-18.6 -16.2

-13.1 -28.4

a The energies in column iii shown in bold are considered to be the most reliable. The energies were calculated by using the methods (see the text for notation) in footnotes b-g. b ff99[GBSA]//ff99[GBSA]. c ff99[GBSA]//B3LYP/6-31G(d,p). d B3LYP/6-31G(d,p)[CPCM]//B3LYP/ 6-31G(d,p)[CPCM]. e B3LYP/6-31G(d,p)[CPCM]//B3LYP/631G(d,p). f B3LYP/6-31G(d,p)[GBSA]//B3LYP/6-31G(d,p). g B3LYP/ 6-31G(d,p)[PB]//B3LYP/6-31G(d,p).

fails to reproduce the solvation energies of large polymeric molecules (proteins and nucleic acids).64 Because of this the choice of the PB method is more attractive in calculations of the binding energies of large metal-nucleic acid complexes. We also found that minimization of the complexes and their monomers using the ff99[GBSA] approach leads to a significant change in both the nonsolvent and solvation energies compared to previously discussed data obtained using ab initio equilibrium geometries: the nonsolvent energy was found to be -32.1 (-160.4) kcal/mol (Figure 5) and the solvation free energy was 35.4 (159.9) kcal/mol for G-Mg(H2O)62+ (DMP--Mg(H2O)62+) (Figure 7). This clearly illustrates that the final values of the interaction energy between two monomers would be affected not only by the methods chosen for the calculations but also by the geometries used. 3.2.4. Total Interaction Energy in a SolVent. The effect of solvation on the total interaction energy is summarized in Figure 8, which shows the energy as a function of distance from the equilibrium position of magnesium ions, and in Table 3 for the equilibrium geometries of the complexes. The energies were calculated according to eq 7 as a sum of the nonsolvent and solvation energies for G-Mg(H2O)62+ and DMP--Mg(H2O)62+ in a dielectric continuum ( ) 78.4). The energies for the systems optimized at the B3LYP/6-31G(d,p)[CPCM] level were -20.2 and -30.3 kcal/mol, respectively. The simplified approach that skips the CPCM-optimization routine and uses gasphase geometries at the B3LYP/6-31G(d,p)[CPCM]//B3LYP/ 6-31G(d,p) level leads to an underestimation of the energies by only 3.3 and 7.5 kcal/mol. The small difference between these two approaches allows us to conclude that the use of gasphase geometries, which is computationally less expensive, is a reasonable approximation in calculating total interaction energies in solution. This is in agreement with comments by

Klicic´ et al. concerning similar calculations of pKa values of organic compounds.57 The classical ff99(GBSA)//B3LYP/6-31G(d,p) method fails to predict a negative interaction energy for either complex. However, the calculations discussed above showed that the differences in solvation energies calculated within the GBSA and PB approximations with use of the B3LYP/6-31G(d,p) geometries are in a good agreement with the CPCM results. Therefore, the main source of error in calculating the total energy comes from an underestimation of the gas-phase interaction term calculated with the AMBER force field. Although relaxation of the complexes and their monomers performed within the fully classical ff99[GBSA] approach lowers the interaction energies to 2.4 and -0.5 kcal/mol for G-Mg(H2O)62+ and DMP-Mg(H2O)62+, respectively, these values are still in poor agreement with the CPCM data. Figure 8 also indicates that the interactions in solution are of short range and exist only within a distance of 2.5-3.0 Å between monomers. Because the ff99 force field is unable to accurately reproduce either the energetic or geometrical aspects of the site-specific binding interactions studied here, we propose an alternative method to estimate the interactions between metal ions and macromolecules. The nonsolvent interaction between an ion and its binding site should be calculated with ab initio methods, while solvation contributions can be estimated by using either the generalized Born or Poisson-Boltzmann approaches. In Figure 8 we also plot the total interaction energy for G-Mg(H2O)62+ and DMP--Mg(H2O)62+ as a function of distance from the equilibrium position of the magnesium ion obtained by the combined B3LYP/6-31G(d,p)[GBSA]//B3LYP/ 6-31G(d,p) method. The results of this approach and its alternative, the B3LYP/6-31G(d,p)[PB]//B3LYP/6-31G(d,p) method, for the complexes at the equilibrium geometries are presented in Table 3 (columns V and VI). Although these combinations of methods still underestimate the total interaction energies compared to the B3LYP/6-31G(d,p)[CPCM] calculations, they still provide reasonable results, are better than the classical methods, and may be applied to large systems. 4. Conclusions The presence of solvent significantly destabilizes the interactions in G-Mg(H2O)62+ and DMP--Mg(H2O)62+. We have shown that to calculate the total interaction energy for these systems, accurate calculation of nonsolvent and solvation contributions is important. While ab initio CPCM calculations are sufficient for small systems, they are not applicable to large biochemical complexes. The simplified approach to calculate

6080 J. Phys. Chem. B, Vol. 108, No. 19, 2004 solvation energies with gas-phase geometries within the CPCM method underestimates their strength by only a few kilocalories per mole for the systems studied here. Force field-based methods, such as AMBER, that are often used to describe the interaction between macromolecules have difficulty describing site-specific interactions with solvated magnesium ions. In particular, we found that the major drawback of the AMBER ff99 force field is an inaccurate description of the electrostatic interaction between monomers at close distances of approach. Another drawback for this system studied here is incomplete cancellation of the attractive (POL and CT) terms by the repulsion (EX). Both problems stem from a classical treatment of the extended electronic clouds of the monomers. While this leads to a significant underestimation of the total interaction energy, the ES and POL components might be improved by inclusion of higher order electrostatic terms and by accounting for polarizability. However, the fact that the charges of the separate monomers are different from those in the complex makes it difficult to see how a purely classical approach could accurately represent the binding energies. The good agreement between the calculation of solvation energy differences with the Poisson-Boltzmann or generalized Born approach with the CPCM data suggests that the proposed method can readily be extended to study the binding between cations and polynucleotides in the presence of electrolyte. This could be done by using a scheme computationally similar to the Bashford-Karplus (BK) approach to macromolecular pKa calculations.65,66 For example, consider the site-specific binding of a hydrated magnesium ion to the N7 atom of a guanine base in a RNA segment. One would first identify the “titratable site” as a G-Mg(H2O)62+ complex and calculate the free energy of desolvation used for the separate guanine and magnesium monomers according to the PB theory. The interaction energy for complex formation in the gas phase could then be evaluated by using DFT or some other quantum mechanical approach, followed by a second PB calculation of the free energy of solvation of the guanine-magnesium complex. The free energy cycle implicit in this method, identical with that shown in Figure 3, mirrors that of Bashford and Karplus. The only computational distinction is that in the BK approach to pKa calculations, the pKa of the model compound is known (although it could, in principle, be computed), whereas in site-binding complex formation, as discussed here, the free energy for complex formation (e.g., the pKMg for N7 of guanine) is not known a priori and must be calculated explicitly. Finally, our method may be applied to a wider class of systems than cation-polynucleotide binding. For example, if a simple quantum mechanical gas-phase calculation of a particular ligand-substrate complex suggests that the effects of charge transfer and polarization might make significant contributions to the binding energy, then the approach used here should be of value. Acknowledgment. This work was supported in part by NIH grant GM 29079-20. The authors also gratefully acknowledge the University of Kentucky HP Superdome and the Advanced Biomedical Computing Center (ABCC) at the NCI-Frederick for allocation of computer time. Supporting Information Available: Listing of RESP charges derived at the HF/6-31G(d) level for guanine, tgDMP-, and Mg(H2O)62+ monomers. This material is available free of charge via the Internet at http://pubs.acs.org.

Petrov et al. References and Notes (1) Juneau, K.; Podell, E.; Harrington, D. J.; Cech, T. R. Structure 2001, 9, 221. (2) Gorman, M. A.; Morera, S.; Rothwell, D. G.; delaFortelle, E.; Mol, C. D.; Tainer, J. A.; Hickson, I. D.; Freemont, P. S. EMBO J. 1997, 16, 6548. (3) Shui, X. Q.; McFail-Isom, L.; Hu, G. G.; Williams, L. D. Biochemistry 1998, 37, 8341. (4) Ennifar, E.; Walter, P.; Dumas, P. Nucleic Acids Res. 2003, 31, 2671. (5) Manning, G. S. Quart. ReV. Biophys. 1978, 11, 179. (6) Dornberger, U.; Spackova, N.; Walter, A.; Gollmick, F. A.; Sponer, J.; Fritzsche, H. J. Biomol. Struct. Dyn. 2001, 19, 159. (7) Spackova, N.; Cheatham, T. E.; Ryjacek, F.; Lankas, F.; van Meervelt, L.; Hobza, P.; Sponer, J. J. Am. Chem. Soc. 2003, 125, 1759. (8) Essex, J. W.; Severance, D. L.; TiradoRives, J.; Jorgensen, W. L. J. Phys. Chem. B 1997, 101, 9663. (9) Massova, I.; Kollman, P. A. Perspect. Drug DiscoV. 2000, 18, 113. (10) Tsui, V.; Case, D. A. J. Phys. Chem. B 2001, 105, 11314. (11) Petrov, A. S.; Lamm, G.; Pack, G. R. J. Phys. Chem. B 2002, 106, 3294. (12) Rulisek, L.; Sponer, J. J. Phys. Chem. B 2003, 107, 1913. (13) Haider, S.; Parkinson, G. N.; Neidle, S. J. Mol. Biol. 2002, 320, 189. (14) Sponer, J.; Leszczynski, J.; Hobza, P. Biopolymers 2001, 61, 3. (15) Siegbahn, P. E. M.; Blomberg, M. R. A. Annu. ReV. Phys. Chem. 1999, 50, 221. (16) Sponer, J.; Sabat, M.; Burda, J. V.; Leszczynski, J.; Hobza, P.; Lippert, B. J. Biol. Inorg. Chem. 1999, 4, 537. (17) Sponer, J.; Burda, J. V.; Leszczynski, J.; Hobza, P. J. Biomol. Struct. Dyn. 1999, 17, 61. (18) Sponer, J. E.; Glahe, F.; Leszczynski, J.; Lippert, B.; Sponer, J. J. Phys. Chem. B 2001, 105, 12171. (19) Munoz, J.; Sponer, J.; Hobza, P.; Orozco, M.; Luque, F. J. J. Phys. Chem. B 2001, 105, 6051. (20) Cossi, M.; Barone, V.; Cammi, R.; Tomasi, J. Chem. Phys. Lett. 1996, 255, 327. (21) Cances, E.; Mennucci, B.; Tomasi, J. J. Chem. Phys. 1997, 107, 3032. (22) Foresman, J. B.; Keith, T. A.; Wiberg, K. B.; Snoonian, J.; Frisch, M. J. J. Phys. Chem. 1996, 100, 16098. (23) Barone, V.; Cossi, M. J. Phys. Chem. A 1998, 102, 1995. (24) Gresh, N.; Sponer, J. E.; Spackova, N.; Leszczynski, J.; Sponer, J. J. Phys. Chem. B 2003, 107, 8669. (25) Tunega, D.; Haberhauer, G.; Gerzabek, M.; Lischka, H. J. Phys. Chem. A 2000, 104, 6824. (26) Aquino, A. J. A.; Tunega, D.; Haberhauer, G.; Gerzabek, M. H.; Lischka, H. Phys. Chem. Chem. Phys. 2001, 3, 1979. (27) Nakabayashi, T.; Sato, H.; Hirata, F.; Nishi, N. J. Phys. Chem. A 2001, 105, 245. (28) Aquino, A. J. A.; Tunega, D.; Haberhauer, G.; Gerzabek, M. H.; Lischka, H. J. Phys. Chem. A 2002, 106, 1862. (29) Martinez, J. M.; Pappalardo, R. R.; Marcos, E. S. J. Am. Chem. Soc. 1999, 121, 3175. (30) Wasserman, E.; Rustad, J. R.; Xantheas, S. S. J. Chem. Phys. 1997, 106, 9769. (31) Floris, F. M.; Martinez, J. M.; Tomasi, J. J. Chem. Phys. 2002, 116, 5448. (32) Becke, A. D. Phys. ReV. A 1988, 38, 3098. (33) Lee, C. T.; Yang, W. T.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (34) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A., Jr.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Salvador, P.; Dannenberg, J. J.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Baboul, A. J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Andres, J. L.; Gonzalez, C.; Head-Gordon, M.; Replogle, E. S.; Pople, J. A. Gaussian 98, Revision A.11; Gaussian, Inc.: Pittsburgh, PA, 2001. (35) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553. (36) van Duijneveldt-van de Rijdt, J. G. C. M.; van Duijneveldt, F. B. Ab Initio Methods Applied to Hydrogen Bonded Systems. In Theoretical Treatments of Hydrogen Bonding; Hadzi, D., Ed.; Wiley: Chichester, UK, 1997; pp 13-47. (37) Morokuma, K. J. Chem. Phys. 1971, 55, 1236. (38) Kitaura, K.; Morokuma, K. Int. J. Quantum Chem. 1976, 10, 325.

Calculations of Mg-Nucleic Acid Site Binding in Solution (39) 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. J.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. J. Comput. Chem. 1993, 14, 1347. (40) Andzelm, J.; Kolmel, C.; Klamt, A. J. Chem. Phys. 1995, 103, 9312. (41) Barone, V.; Cossi, M.; Tomasi, J. J. Chem. Phys. 1997, 107, 3210. (42) Case, D. A.; Pearlman, D. A.; Cadwell, J. W.; Cheatham, T. E., III; Wang, J.; Ross, W. S.; Simmerling, C. L.; Darden, T. A.; Merz, K. M.; Stanton, R. V.; Cheng, A. L.; Vincent, J. J.; Crowley, M.; Tsui, V.; Gohkle, H.; Radmer, R. J.; Duan, Y.; Pitera, J.; Massova, I.; Seibel, G. L.; Singh, U. C.; Weiner, P. K.; Kollman, P. A. AMBER 7; University of California: San Francisco, CA, 2002. (43) Wang, J. M.; Cieplak, P.; Kollman, P. A. J. Comput. Chem. 2000, 21, 1049. (44) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Kollman, P. A. J. Am. Chem. Soc. 1993, 115, 9620. (45) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269. (46) Tsui, V.; Case, D. A. Biopolymers 2000, 56, 275. (47) Weiser, J.; Shenkin, P. S.; Still, W. C. J. Comput. Chem. 1999, 20, 217. (48) Davis, M. E.; Madura, J. D.; Luty, B. A.; McCammon, J. A. Comput. Phys. Commun. 1991, 62, 187. (49) Sitkoff, D.; Sharp, K. A.; Honig, B. J. Phys. Chem. 1994, 98, 1978. (50) Srinivasan, J.; Cheatham, T. E.; Cieplak, P.; Kollman, P. A.; Case, D. A. J. Am. Chem. Soc. 1998, 120, 9401.

J. Phys. Chem. B, Vol. 108, No. 19, 2004 6081 (51) Lundell, J.; Latajka, Z. Chem. Phys. 2001, 263, 221. (52) Born, M. Z. Phys. 1920, 1, 45. (53) Scheiner, S. Hydrogen bonding: a theoretical perspectiVe; Oxford University Press: New York, 1997. (54) Burda, J. V.; Sponer, J.; Leszczynski, J. J. Biol. Inorg. Chem. 2000, 5, 178. (55) Sponer, J. E.; Leszczynski, J.; Glahe, F.; Lippert, B.; Sponer, J. Inorg. Chem. 2001, 40, 3269. (56) Lamm, G.; Pack, G. R. J. Phys. Chem. B 1997, 101, 959. (57) Klicic, J. J.; Friesner, R. A.; Liu, S. Y.; Guida, W. C. J. Phys. Chem. A 2002, 106, 1327. (58) Jorgensen, W. L.; Buckner, J. K.; Boudon, S.; TiradoRives, J. J. Chem. Phys. 1988, 89, 3742. (59) Pack, G. R.; Wang, H.; Rein, R. Chem. Phys. Lett. 1972, 17, 381. (60) Jensen, F. Introduction to computational chemistry; John Wiley: New York, 1999. (61) Rein, R.; Rabinowitz, J. R.; Swissler, T. J. J. Theor. Biol. 1972, 34, 215. (62) Koch, U.; Egert, E. J. Comput. Chem. 1995, 16, 937. (63) Aleman, C.; Orozco, M.; Luque, F. J. Chem. Phys. 1994, 189, 573. (64) Onufriev, A.; Bashford, D.; Case, D. A. J. Phys. Chem. B 2000, 104, 3712. (65) Bashford, D.; Karplus, M. Biochemistry 1990, 29, 10219. (66) Bashford, D.; Karplus, M. J. Phys. Chem. 1991, 95, 9556.