Complementary Electrostatics for the Study of DNA Base-Pair

A Novel Theoretical Model for Molecular Recognition of Multiple-Site Interacting .... Journal of Quantum Chemistry 2009 109 (10.1002/qua.v109:10), 223...
0 downloads 0 Views 348KB Size
Download PDF
3298

J. Phys. Chem. B 1997, 101, 3298-3303

Complementary Electrostatics for the Study of DNA Base-Pair Interactions Shridhar R. Gadre* and Savita S. Pundlik Department of Chemistry, UniVersity of Pune, Pune-411 007, India ReceiVed: December 13, 1996; In Final Form: February 3, 1997X

The complementary “lock and key” patterns of weakly interacting molecules are explored by mapping the topography of respective molecular electrostatic potentials (MESP). A new model, viz., electrostatic potential for intermolecular complexation (EPIC), which incorporates such MESP features, has been employed for studying interactions between DNA base pairs. A wide variety of pairs of bases involving adenine (A), guanine (G), and cytosine (C) are chosen as test cases. The interaction energy within the EPIC model is expressed as (1/2){∑VA,iqB,i + ∑VB,iqA,i}, where VA,i is the MESP value due to A at the ith atom of molecule B and qB,i is the MESP-derived charge at this site. The interaction energies and geometrical parameters obtained by this model agree remarkably well with the corresponding literature values obtained by full geometry optimization at the ab initio level. Being intuitively appealing and simple in application, the EPIC model seems to have the potential of being a good predictive tool for investigating weak intermolecular interactions.

1. Introduction For unraveling the biological functioning of carriers of genetic information, it is necessary to understand their chemical structure.1 The latter is facilitated by a variety of experimental and theoretical tools such as X-ray crystallography,2 spectroscopy,3 and potential energy calculations.4 Hydrogen bonding between the DNA base pairs is an extremely important area of research since it is associated with the stabilization of the DNA molecule, with its mutations and interaction with drugs, proteins, etc. Information regarding energetic aspects of base-pair formation is available from gas phase experiments5 as well as from studies in nonpolar solvents.6 The interatomic distances in the H-bonded as well as stacked base pairs are known from the X-ray data.1 It has been pointed out7a that due to the inherent difficulties accompanying the experimental investigations, one has to supplement the latter with the appropriate theoretical treatments. Further, reliable information can only be obtained by employing the ab initio theory with a reasonable basis set (such as 6-31G**) incorporating polarization functions. The primary stage in any electronic structure calculations is that of geometry optimization. This demands a very high computational power for a reliable ab initio treatment of DNA base pairs. A noteworthy attempt7 in this direction has recently been made on a Cray Y-MP and other RISC-based machines, and details regarding energies and geometries of base pairs are now available. However, the properties of any new derivative of a given base or a different conformation of a base pair can be known only after going through such a computation-intensive job. It is quite impractical, as of today, to even think of ab initio geometry optimization for investigating the stacking of base pairs, albeit in the gas phase. A need is felt, therefore, for developing reasonably reliable yet computationally less expensive methods for weak biomolecular interactions. An obvious alternative for the study of large biological systems is offered by use of empirical or semiempirical procedures, as suggested earlier by Hobza and Zahradnik.8a Recently,8b they have compared the applicabilities of various empirical potential functions as well as semiempirical methods for obtaining binding energies of H-bonded DNA base pairs. X

Abstract published in AdVance ACS Abstracts, March 15, 1997.

S1089-5647(96)04064-3 CCC: $14.00

In general, however, these methods are found to be rather unsatisfactory. For example, the MM calculations9 can correcly predict the preference of WC (Watson-Crick) base pairing of d(CGTACG)2 with triostin A. Such calculations, however, cannot explain why DNA intercalations occur with a neighbor exclusion rule.9 The pair-potential-based MM2 method cannot correctly predict the energy associated with the H-bond.10 The relative ordering of energies of different tautomers of uracil is not correct within the semiempirical calculations.11 This problem can be partly circumvented by constructing better empirical parameters. Thus, the knowledge gained by investigating systems of limited size, with a high-level ab initio theory, may be used for reparametrization to obtain reliable empirical parameters. Although parametrization of all the component energy terms is quite a task,12 it is expected to be foolproof. The theory then should be able to correctly predict, for example, the relative energies of various binding sites in proteins, energetics for various conformations in molecules, and geometric dispositions of two interacting species. Unfortunately, it has been found to be difficult to form a universal set of parameters12 that could account for various possible modes of interaction in a given system. For example, with AMBER, a different set of parameters for atoms forming H-bonds was found to be essential.13 The parameters may not be transferable even on distinguishing all the atoms in biomolecules.14 Thus, even though the same type of CO‚‚‚NH bonds are present in G-G and T-T pairs, they differ greatly in stability (-30 and -16 kcal/mol, respectively). Also, it has been recently remarked7b that empirical potentials do not consider attractive interactions involving out of plane amino group hydrogens of DNA bases, resulting in a structure of guanine molecules that is 1.8 kcal/ mol higher in energy. It is therefore worthwhile considering other procedures that make use of the ab initio wave function for each of the participating species for assessing an interaction. This approach has been employed to varying extent in the construction of different models, for example, in the seminal work due to Buckingham and Fowler (B-F),15 in the molecular mechanics for clusters model (MMC) by Dykstra,16 and in the solvation potential method of Alhambra, Luque, and Orozco.17 Each one of these models has some advantages as well as limitations, which have been reviewed recently18 by Gadre and co-workers. There seems to be a general consensus that molecular recognition and weak bond formation between two © 1997 American Chemical Society

Study of DNA Base-Pair Interactions

J. Phys. Chem. B, Vol. 101, No. 16, 1997 3299

species are dominated by electrostatics. The electrostatic component is found19,20 to bear a linear relationship with the SCF binding energy. The former provides a correct ordering of energies in different chemical sites during protonation of three-membered-ring molecules.19 Zahradnik and Hobza14 had recommended the use of the B-F model for predicting structures of H-bonded complexes of moderate size. There have been notable successes of these models. For example, the structures of benzene-water21 and benzene-benzene22 dimers obtained by rotational spectroscopy agree with the molecular mechanics calculations (MMC) by Dykstra.16 Similarly, the B-F model is able to predict the experimental structures of many complexes23 as well as yield the interaction energies that are close to the ones obtained using ab initio calculations. Despite the success for small molecules, the B-F model may not turn out to be very suitable for probing the base-pair interactions since it involves rolling the entire van der Waals (vdW) surface of one molecule over the other. This procedure would be extremely time consuming especially for the interactions between two base pairs. This seems to be the apparent reason for not applying B-F model toward the study of biomolecules. Further, the electrostatic potentials in the B-F model are constructed using distributed multipole analysis24 (DMA), which itself is approximate in nature. Hence, MESP obtained by DMA may often lack the finer topographical features,25 which could be important for the description of binding patterns with other molecules. It is also observed that although the correlated DMA (CDMA) potentials26 offer good insights into the stabilization of various base pairs, some regions on the potential energy surface are not reproduced satisfactorily. The present form of the solvation potential model by Alhambra, Luque, and Orozco also suffers from a limitation, in that water has to be one of the interacting species. The MMC model has apparently not been tried out on complexes involving DNA bases, perhaps owing to the complexity of the systems. In this paper, we present the results of application of a new electrostatics-based model, electrostatic potential for intermolecular complexation (EPIC),18,32 to the investigations on basepair interactions. 2. The EPIC Model This model is based on full utilization of the spatial characteristics of the molecular electrostatic potential (MESP). The latter is a point property, which is defined within the quantum chemical framework as follows:

V(r) )

∑ A |R

ZA

A

- r|

∫|r′ - r|d3r′ F(r′)

-

In this expression, ZA stands for the positive nuclear charge of atom A located at RA in the molecule and F(r) is the molecular electronic charge density. Both these quantities are separately positive; hence, V(r) may be positive, negative, or zero. This engenders rich topographical features of V(r) in a region surrounding the given molecule. These can be brought out by rigorous topographical analysis27 involving location and identification of critical points (CPs). The latter are the points in space where all the first-order partial derivatives of the given scalar field, viz., MESP, vanish. The type of a CP is known by analyzing the eigenvalues of the matrix of second-order derivatives of V(r) (Hessian). If all the three eigenvalues are nonzero (rank r ) 3), it is a nondegenerate CP (r,s) having a signature, s, of +3, -3, +1, or -1. For example, (3,+3) represents a minimum, and (3,-3) stands for a maximum in V(r), while (3,+1) and (3,-1) denote saddle points.

MESP has long been recognized as the quantity that governs some classes of chemical reactions.28 Politzer and co-workers29 have used this scalar field extensively for describing various phenomena including electrophilic substitution reactions for substituted benzenes. The applications of MESP to various fields including reactivity, solvation, interactions, and structureactivity relationships have been reviewed by Na`ray-Szabo` and Ferenczy.30 It has also been recently pointed out31 that one should consider electronic moments higher than order one (dipole moments) in the analysis of noncovalent interactions. However, the entirety of topographical features of MESP is rarely exploited in the study of such reactions. To the best of the knowledge of the authors, such a topographical analysis has not been applied in the study of weak interactions, especially those between DNA bases, though, at times, a glancing reference to MESP is made. The new, EPIC model18,32 aims at predicting the interaction patterns between two species using MESP features of each of the subsystems. The CPs are located and characterized for the two monomers. Further, MESP-driven atom-centered point charges are obtained for these using the individual molecular wave function. The two molecules A and B are positioned so that the electrophilic atoms such as H in A are close to the charge concentration sites (minima) in B and vice versa. This is achieved by translating A or B with respect to the other so that a weighted sum of the functions f, ∑f(dij), dij being the distance between the ith electrophilic atom in A (or B) and the jth negative-valued CP in B (or A), is minimum. To avoid interpenetration of the two species, heavy atoms are embedded in the respective vdW spheres. After this initial positioning, the interaction energy expressed by E ) (1/2){∑VA,iqB,i + ∑VB,iqA,i} is minimized. This expression for the interaction energy is such that both A and B interchangeably play the roles of “lock” and “key”. E is minimized by keeping either A or B stationary and giving all the translational and rotational degrees of freedom to the other molecule, without allowing the crossing of vdW boundaries. A binary complex is thus obtained using complementary electrostatic features of the two molecules. 3. Method The base pairs containing adenine (A), guanine (G), and cytosine (C) are studied in the present work. These have been randomly selected, and there is no specific reason for excluding thymine (T) during the testing of the EPIC model. Geometries of A, G, and C, assumed to be planar, are optimized at the ab initio HF/6-31G* level using GAUSSIAN94.33 Single-point SCF calculations are carried out at the HF/6-31G++(d,p) level using UNIMOL.34 This basis set consists of diffuse polarization functions and is found to be quite suitable8 for the study of molecular interactions. Since electron correlation does not affect significantly the MESP CP characteristics35 for the level of basis set chosen here, post-SCF treatments are not performed. The SCF-level molecular wave functions are used for obtaining the positions and nature of MESP CPs of each molecule using the program UNIPROP.36 The MESP-driven atom-centered charges are obtained using the GRID37 program. Instead of excluding the double vdW surface, as is conventionally done, a more suitable selection of grid points is made. This is done to ensure that the essential electrostatic features of the base pairs are taken into account while fitting the charges. All this information is used for obtaining minimum energy structures for the bases involving A, G, and C, via a cyclic relaxation routine. Scheme 1 presents the steps involved in developing an EPIC model for a weak binary complex. Several runs (step 6, Scheme 1) with different initial guesses (step 5, Scheme 1) are performed so as to avoid local minima. The sketches of the dimer geometries

3300 J. Phys. Chem. B, Vol. 101, No. 16, 1997

Gadre and Pundlik

Figure 1. Structures and atom numbering for three DNA bases: (A) adenine, (B) guanine, and (C) cytosine. The asterisks indicate locations of MESP minima; the negative-valued (3,+1) saddles are marked by a + sign. (cf. Table 1 for the characteristics and MESP values at the CPs. The bond lengths and CP positions are not to scale.)

TABLE 1: Details of Negative MESP Valued Critical Points (CPs) in A, G, and C Marked in Figure 1 (See Text for Details) CP and naturea

distance (Å)b

MESP at CP (au)

adenine

(a) (3,+3) (b) (3,+3) (c) (3,+3) (d) (3,+3) (e) (3,+3)

1.26 1.50 1.50 1.25 1.26

-0.091 -0.027 -0.027 -0.101 -0.099

guanine

(a) (3,+3) (b) (3,+3) (c) (3,+1) (d) (3,+3)

1.24 1.25 1.28 1.29

-0.096 -0.117 -0.095 -0.119

(e) (3,+3) (f) (3,+1) (g) (3,+3) (h) (3,+3) (a) (3,+3) (b) (3,+3) (c) (3,+3) (f) (3,+1)

1.26 2.00 1.54 1.54 1.28 1.53 1.53 1.27

-0.077 -0.085 -0.003 -0.003 -0.115 -0.017 -0.017 -0.107

(d) (3,+3) (g) (3,+1) (e) (3,+3)

1.23 1.74 1.22

-0.130 -0.095 -0.108

molecule

cytosine

atomic chargesc (au) N1: -0.85, C2: 0.61, N3: -0.83 C4: 0.72, C5: -0.12, N6: -1.02 C6: 0.90, N7: -0.58, C8: 0.32 N9: -0.57, H2: 0.06, H61: 0.43 H62: 0.43, H8: 0.11, H9: 0.38 RMSD: 3.23 × 10-3 N1: -0.86, N2: -1.07, C2: 1.04 N3: -0.81, C4: 0.50, C5: -0.001 C6: 0.79, N7: -0.55, C8: 0.24 N9: -0.46, H1: 0.44, H21: 0.44 O6: -0.63 H22: 0.45, H8: 0.12, H9: 0.37 RMSD: 2.83 × 10-3 N1: -0.73, C2: 1.16, N3: -0.98 N4: -1.23, C4: 1.26, C5: -0.84 C6: 0.36, H1: 0.38, H41: 0.49 H42: 0.48, H5: 0.26, H6: 0.13 O2: -0.74 RMSD: 2.41 × 10-3

a cf. Figure 1 for labels of CP’s. b From the nearest heavy atom. c MESP-derived atomic charges derived by minimizing RMSD ) {∑i)1N(Vimodel - Viab)2/N}0.5, where Vmodel is the charge model MESP and Vab is the ab initio MESP calculated for N grid points.

SCHEME 1

4. Results and Discussion

and the pictures containing iso-valued surfaces of individual monomers are generated with the help of the graphics program UNIVIS38 on a PC-486.

A, G, and C are shown schematically in Figure 1, in which are also marked the locations of negative-valued MESP CPs. Their characteristics are summarized in Table 1. The off-plane minima corresponding to amino-nitrogen lone pairs are not shown in the figures. The positions of CPs represent the sites for hydrogen bonding, whose strength is determined to a great extent by the value of the MESP there. Thus, complexes involving N1 in adenine, N7 and O6 in guanine, and O2 in cytosine are expected to give rise to comparatively stronger binding. The MESP-derived charges for all the atoms in A, G, and C along with the RMSD (root mean square deviation) values are also displayed in Table 1. The iso-surfaces of MESP shown in Figure 2 indicate that regions of negative potential arising due to the nitrogen lonepair minima are quite localized in the case of adenine, whereas these are rather spread out in the case of guanine and cytosine. The picture18 for a six-membered ring with a “hamburger” formed by the negative MESP buns is no longer observed here. The full aromatic π-character is thus conspicuous by its absence, and a possible off-plane interaction would always involve the amino nitrogen. Although the lone pair associated with N6 in adenine has an MESP value comparable to the π-charge minima in benzene (-0.029 au at 6-31G++(d,p)), it is highly localized, as can be seen from Figure 2a. A weak electrostatic field

Study of DNA Base-Pair Interactions

J. Phys. Chem. B, Vol. 101, No. 16, 1997 3301

Figure 2. Iso-MESP surfaces of (a) adenine, (b) guanine, and (c) cytosine, with values of -0.02 au (spreadout region represented with a meshlike surface) and -0.07 au (darker regions).

(iso-surface corresponding to MESP value of -2 × 10-2 au in Figure 2b) is seen to be quite extensive in the case of guanine. In this molecule, the negative MESP corresponding to the N2 lone-pair is neglegibly small. In cytosine, the N4 lone pair minima are comparatively stronger (cf. Table 1). The negatively charged region off the plane of these molecules could have some effect upon the binding patterns; for example, this may be partly resposible for the nonplanarity in base pairs observed recently by Sponer et al.7c However, very high charge concentration sites lie in the molecular planes of these bases, which are primarily responsible for binding with the protons in another species. The three (3,+3) CPs corresponding to ring nitrogens N1, N3, and N7 in A are isolated from each other, and the charge conjugation effects are thus not observed here. In guanine, the charge delocalization in the molecular plane is manifested by the presence of (3,+1) saddles marked (c) and (e) in Figure 1. Similarly, in cytosine, the CPs (a), (g), (d), (f), and (e) form a highly negative region, favorable for proton binding. The base pairs are thus expected to have mostly planar configurations. This is in contrast to the binding features observed for benzene, wherein it is found that a cationic species is held firmly by the π-charge31 over the plane of the ring. Out of the 30 DNA base pairs studied by Sponer et al.7a at the ab initio level, we have randomly selected 10 as test cases with nomenclature identical to that followed in ref 7a. The interaction energies for AA1, AA2, AA3, GA1, GA2, GA3, GA4, GCWC, GC1, and CC as determined by the EPIC model are reported in Table 2. The corresponding geometries for some of the base pairs are depicted in Figure 3. In Table 2 are also included the interaction energies at the HF/6-31G** level in ref 7a for a comparison. We believe39 that the the basis set

TABLE 2: Interaction Energies (kcal/mol) Obtained for Various Base Pairs Using the EPIC Model (Emodel) and Those in Ref 7a (See Text for Details) base pair

Emodel

EHF17a

max. norm

AA1 AA2 AA3 GA1 GA2 GA3 GA4 GCWC GC1 CC

-10.1 -9.3 -8.3 -12.1 -7.9 -10.1 -8.9 -22.7 -10.6 -16.0

-8.7 -8.0 -6.9 -12.6 -7.5 -11.0 -8.8 -25.5 -12.7 -17.3

0.016 0.016 0.050 0.024 0.013 0.019 0.014 0.028 0.019 0.029

effects would not be large beyond the 6-31G* level. It is seen that the energies predicted by the present model are in fairly good agreement with the corresponding ab initio ones, the average error being 1.2 kcal/mol. The maximum percent deviation is observed in the case of AA3, which amounts to about 20%. The rank order of stability for multiple configurations of the given base pair is generally predicted correcly by the model. It turns out from the EPIC model that AA1 is the most stable complex formed by two adenine molecules. This fact, which is an outcome of very sophisticated computation, has a very simple explanation. In AA1, N1 atoms of both the adenines, which are accompanied by the strongest MESP minima, participate in hydrogen bonding. Similarly, AA3 is the most weakly bound among AA1, AA2, and AA3 since the CPs associated with the nitrogen lone pairs of the five-membered ring have comparatively less negative MESP values (cf. Table 1). A similar heuristic explanation may be offered in the case

3302 J. Phys. Chem. B, Vol. 101, No. 16, 1997

Gadre and Pundlik TABLE 3: Distance Parameters Involved in the Hydrogen Bonding in Various Base Pairs base pair AA1 AA2 (b)

(a)

AA3 GA1 GA2 GA3 GA4 GCWC

(c)

(d)

GC1 CC

(e)

Figure 3. Configurations for some of the DNA base pairs obtained using the EPIC model: (a) AA1, (b) AA3, (c) GA4, (d) GCWC, (e) CC.

of other base pairs. The strongest guanine-adenine complex is GA1, in which the deepest minima of both the moieties happen to be the binding sites. GA2 and GA4 are expected to be comparatively less stable since protons in adenine are trapped in shallower MESP regions of guanine, forming local minimum structures. GA1 is endowed with greater stability than GA3 because the proton-donating site for G is stronger with respect to the negative MESP value than the latter. The strongest base pair in Table 2 happens to be GCWC not only because there are three hydrogen bonds involved in the complex formation but also since the deepest MESP minima of both G and C are accessed by the protons of each other. In short, the binding patterns in a complex can be guessed from the knowledge of the MESP topography of the individual interacting moieties. It is thus not surprising that the rank order of various base pairs with respect to energetics is generally predicted correctly by the EPIC model. Table 3 gives values of distances involved in weak hydrogen bonding for the base pairs as obtained by this model. The agreement with the corresponding values reported in ref 7a is quite encouraging. Also encouraging is the fact that various intermolecular bonds have lengths that are within the range of values determined by X-ray methods.1 On comparing these bonding distances with the distances of MESP CPs from the nearest heavy atoms, it may be inferred that the distance between two closest heavy atoms associated with a weak hydrogen bond in a complex is approximately equal to the sum of the distance of the heavy atom from the CP and the length of a covalent bond involving that hydrogen. Scaling the radii of these hydrogen atoms while building a model plays a crucial role in obtaining the hydrogen bond distances. The construction of optimum base-pair configurations within the EPIC model is quite economical in computer time. For a

H bond

distances (Å) with EPIC model

distances (Å) in ref 7a

N6(H)‚‚‚N1 N1‚‚‚(H)N6 N1‚‚‚(H)N6 N6(H)‚‚‚N7 N6(H)‚‚‚N7 N7‚‚‚(H)N6 O6‚‚‚H(N6) N1(H)‚‚‚N1 N2(H)‚‚‚N7 N3‚‚‚H(N6) O6‚‚‚H(N6) N1(H)‚‚‚N7 N2(H)‚‚‚N1 N3‚‚‚H(N6) N2(H)‚‚‚O2 N1(H)‚‚‚N3 O6‚‚‚(H)N4 N2(H)‚‚‚N3 N3‚‚‚H(N4) N3‚‚‚(H)N4 N4(H)‚‚‚N3

3.18 3.18 3.13 3.10 3.10 3.10 3.03 3.16 3.15 3.25 3.05 3.18 3.16 3.30 3.08 3.18 3.16 3.17 3.67 3.18 3.18

3.16 3.16 3.17 3.16 3.18 3.18 2.95 3.19 3.12 3.24 2.94 3.24 3.10 3.19 3.02 3.04 2.92 3.01 3.22 3.05 3.05

single base pair, it takes a total of typically 10 h (including SCF calculations for the two molecules, their topographic analysis, and applying the EPIC model) on the HP9000/715 machine, which would be less than half an hour on a Cray. This time is indeed extremely small in comparison to the Cray timings (∼a few 100 h) required for a complete gradient-based optimization reported by Sponer et al.7a To check the suitability of the geometries obtained by the EPIC model as starting guesses for ab initio optimization, these have been subjected to single-point SCF calculations with the 6-31G++(d,p) basis set. The respective maximum norms of energy gradients (cf. Table 2) are observed to be typically on the order of 10-2. It is expected, therefore, that such starting guesses will lead to an enormous amount of saving of valuable computational time. This is a very important advantage, and one may think of complete optimization of a complex that involves a few hundred atomic basis functions. 5. Concluding Remarks A firm connection is established between the phenomenon of molecular recognition with the corresponding three-dimensional MESP distribution. The topography of the latter allows its complete description as opposed to the conventional methodology involving DMA or plotting of this scalar field on different discrete planes. The negative-valued CPs in MESP, obtained by topographic investigation, act as attractors for protons in another molecule, thus leading to a “lock and key” type binding of the two species, possibly each playing the dual role of “lock” as well as “key”. The extent of charge concentration at the CPs (MESP values) determines the strength of the hydrogen bond. The present EPIC model exploits these MESP features as well as the atom-centered point charges, to predict the orientations and the interaction energies of weak complexes. Thus, the model is developed solely on the basis of electrostatic properties of the two molecules. It is even possible to enhance the predictive capacity of the EPIC model by improving the MESP-driven charge model. The results of application of this model to DNA base pairs involving adenine, guanine, and cytosine have been presented. The agreement with the corresponding ab initio data suggests that the contributions essential to the interaction are included in this model. Also, an a priori and quick semiquantitative analysis of an interaction

Study of DNA Base-Pair Interactions can be performed by using such an electrostatic approach. The configurations generated by the model would act as very good starting guesses even if complete geometry optimization is desired at an ab initio level of theory. It is thus expected that the EPIC model would be of immense use in studies of DNA base trimers as well as base pair-base pair interactions. Such studies are currently being pursued in our laboratory. Acknowledgment. The authors are thankful to Dr. L. J. Bartolotti for help in carrying out single-point SCF calculations with gradient evaluation for DNA base pairs. The financial support from the Council for Scientific and Industrial Research (CSIR), New Delhi, in the form of a fellowship to S.S.P. as well as a grant to S.R.G. (80/(0022)/96/EMR-II) is greatfully acknowledged. References and Notes (1) Saenger, W. In Principles of Nucleic Acid Structures, Springer AdVanced Texts in Chemistry Series; Cantor, C. R., Ed.; Springer: New York, 1984. (2) (a) Recke, G. N.; Marsh, R. E. Acta Crystallogr. 1965, 20, 703. (b) Richards, F. M. J. Mol. Biol. 1965, 37, 225. (3) (a) Tinoco, I,. Jr.; Bustamente, C. In Methods in Stuctural Molecular Biology; Davies, D. B., Saenger, W., Danyluk, S. S., Eds.; Plenum: London, 1981; pp 269-305. (b) Peticolas, W. L. In Methods in Stuctural Molecular Biology; Davies, D. B., Saenger, W., Danyluk, S. S., Eds.; Plenum: London, 1981; pp 237-266. (4) (a) Scheraga, H. A. AdV. Phys. Org. Chem. 1968, 6, 103. (b) Broyde, S. B.; Wartell, R. M.; Stellman, S. D.; Hingerty, B.; Langridge, R. Biopolymers 1975, 14, 1597. (5) (a) Yanson, I. K.; Teplitsky, A. B.; Sukhodub, L. F. Biopolymers 1979, 18, 1149. (b) Dey, M.; Moritz, F.; Grotemeyer, J.; Schlag, E. W. J. Am. Chem. Soc. 1994, 116, 9211. (6) (a) Mayersky, A. A.; Sukhorukov, B. I. Nucl. Acids Res. 1980, 8, 3029. (b) Imoto, T. Biochim. Biophys. 1977, 475, 409. (7) (a) Sponer, J.; Leszczynski, J.; Hobza P. J. Phys. Chem. 1996, 100, 1965. (b) Sponer, J.; Hobza, P. Int. J. Quantum Chem. 1996, 57, 959. (c) Sponer, J.; Florian, J.; Hobza, P.; Leszczynski, J. J. Biomol. Struct. Dynam. 1996, 13, 827. (8) (a) Hobza, P.; Zahradnik, R. Chem. ReV. 1988, 68, 871. (b) Hobza, P.; Hubalek, F.; Kabelac, M.; Mejzlik, P.; Sponer, J.; Vondrasek, J. Chem. Phys. Lett. 1996, 257, 31. (9) Kollman, P. A.; Kenneth, M. M., Jr. Acc. Chem. Res. 1990, 23, 246. (10) Allinger, N. L.; Kok, R. A.; Imam, M. R. J. Comput. Chem. 1988, 9, 591. (11) Leszczynski J. J. Phys. Chem. 1992, 96, 1649.

J. Phys. Chem. B, Vol. 101, No. 16, 1997 3303 (12) Hobza, P.; Zahradnik, R. Int. J. Quantum Chem. 1992, 42, 581. (13) Weiner, P. K.; Kollman, P. A. J. Comput. Chem. 1981, 2, 287. (14) Zahradnik, R.; Hobza, P. Int. J. Quantum Chem. 1986, 29, 663. (15) (a) Buckingham, A. D.; Fowler, P. W. J. Chem. Phys. 1983, 79, 6426. (b) Buckingham, A. D.; Fowler, P. W. Can. J. Chem. 1985, 63, 2018. (16) Augspurger, J. D.; Dykstra, C. E.; Zweir, T. S. J. Phys. Chem. 1992, 96, 7252. (17) Alhambra, C.; Luque, F. J.; Orozco, M. J. Phys. Chem. 1995, 99, 3084. (18) Gadre, S. R.; Bhadane, P.; Pundlik, S. S.; Pingale, S. S. In Molecular Electrostatic Potentials: Theory and Applications; Murray, J. S., Sen, K. D., Eds.; Elsevier: Amsterdam, 1996. (19) Ghio, C.; Tomasi, J. Theor. Chim. Acta (Berlin) 1973, 30, 151. (20) Alagona, G.; Tani, A. J. Chem. Phys. 1981, 74, 3980. (21) Gutowsky, H. S.; Emilsson, T.; Arunan, E. J. Chem. Phys. 1993, 99, 4883. (22) Arunan, E.; Gutowsky, H. S. J. Chem. Phys. 1993, 98, 4294. (23) Hurst, G. J. B.; Fowler, P. W.; Stone, A. J.; Buckingham, A. D. Int. J. Quantum Chem. 1986, 29, 1223. (24) Stone, A. J. Chem. Phys. Lett. 1981, 83, 233. (25) Jug, K.; Kolle, C. Proc. Indian Acad. Sci. (Chem. Sci.) 1994, 106, 283. (26) Hobza, P.; Sponer, J.; Polasek, M. J. Am. Chem. Soc. 1995, 117, 792. (27) (a) Stewart, I. Sci. Am. 1991, 264, 123. (b) Bader, R. F. W.; NguyenDang, T. T.; Tal, Y. A. Rep. Phys. 1981, 44, 893. (28) Politzer, P.; Truhlar, D. G. Chemical Applications of Atomic and Molecular Electrostatic Potentials; Plenum: New York, 1981. (29) Politzer, P.; Abrahmsen, L.; Sjo¨berg, P. J. Am. Chem. Soc. 1984, 106, 855. (30) Na`ray-Szabo`, G.; Ferenczy, G. G. Chem. ReV. 1995, 95, 829. (31) Dougherty, D. A. Science 1996, 271, 163. (32) Gadre, S. R.; Pundlik, S. S.; Bartolotti, L. J. Unpublished work. (33) Frisch, M. J.; Trucks, G. W.; Head-Gordon, M.; Gill, P. M. V.; Wong, M. W.; Foresman, J. B.; Johnson, B. G.; Schlegel, H. B.; Rob, M. A.; Replogle, E. S.; Gomperts, R.; Andres, J. L.; Raghavachari, K.; Binkley, J. S.; Gonzales, C.; Martin, R. L.; Fox, D. J.; Defrees, D. J.; Baker, J.; Stewart J. J. P.; Pople, J. GAUSSIAN94; Gaussian Inc.: Pittsburgh, PA, 1994. (34) Limaye, A. C.; Gadre, S. R. UNIMOL. J. Chem. Phys. 1994, 100, 1303. (35) Kulkarni, S. A. Chem. Phys. Lett. 1996, 254, 268. (36) Shirsat R. N.; Bapat, S. V.; Gadre, S. R. UNIPROP. Chem. Phys. Lett. 1992, 200, 373. (37) Chipot, C. GRID: The FORTRAN program for fitting charges to molecular electrostatic potentials and fields; University de Nancy I: France, 1992. (38) Limaye, A. C.; Inamdar, P. V.; Dattawadkar, S. M.; Gadre, S. R. UNIVIS: a visualization package. J. Mol. Graph. 1996, 14. (39) Gadre, S. R.; Kulkarni, S. A.; Suresh, C. H.; Shrivastava, I. H. Chem. Phys. Lett. 1994, 239, 273.