9554
J. Phys. Chem. B 2000, 104, 9554-9563
Are Many-Body Effects Important in Protein Folding? Arjan van der Vaart,† Badry D. Bursulaya,‡ Charles L. Brooks, III,‡ and Kenneth M. Merz, Jr.*,† Department of Chemistry, The PennsylVania State UniVersity, 152 DaVey Laboratory, UniVersity Park, PennsylVania 16802, and Department of Molecular Biology, The Scripps Research Institute, 10550 North Torrey Pines Road, La Jolla, California 92037 ReceiVed: March 29, 2000; In Final Form: July 10, 2000
In this article we investigate the importance of many-body and nonclassical effects, such as polarization and charge transfer, on the folding of the betanova protein. Our calculations show that these effects are crucial in stabilization of the system. Moreover, both polarization and charge transfer significantly alter the charge distribution of the system. Our detailed study shows that these fluctuations in charge are solely dependent on the local environment and not on the overall fold of the protein. Moreover, the contributions of polarization and charge transfer are roughly constant during the protein folding process. This means that the folding driving force is largely determined by the electrostatic energy. Our findings indicate that the folding of betanova can be accurately described by effective two-body potentials, despite the absence of explicit polarization and charge transfer in these models.
Introduction Recently, the folding of several small proteins and polypeptides has been studied by all atom simulations.1-10 These simulations yielded not only experimentally available parameters, such as folding free energies and folding intermediates, but also gained detailed insights into processes for which experimental resolution is still too low. Indeed, in computer simulations folding pathways can be followed at a time resolution of several femtoseconds under controlled conditions (such as temperature, pressure, and pH). Despite the potential power of simulation methods, several fundamental questions on the applicability of computer simulation to the protein folding problem still need to be addressed. For example, is the potential function accurate enough for both native and unfolded states? Does the simulated pathway (or reaction coordinate) correspond to the experimental pathway (or reaction coordinate)? How can proper sampling of phase space be ensured along the (un)folding reaction coordinate? In this article we will address the accuracy of the potential function for protein folding simulations. Thus far, all reported folding simulations were performed with effective two-body potentials.11-14 These potentials contain bond stretching, angle bending, torsional rotation, Coulomb, and Lennard-Jones type interactions. For computational convenience, only interactions between pairs of atoms are included, and many-body interactions such as polarization are neglected. Also, the charge distribution of the system is assumed to be constant. Quantum mechanical calculations have shown that the true nature of interactions in small hydrogen-bonded clusters is significantly different than effective two-body potentials suggest. For instance, many-body effects such as polarization are estimated to contribute up to 30% of the total interaction energy in small water clusters.15 For other clusters, polarization accounts * To whom correspondence should be addressed. † The Pennsylvania State University. ‡ The Scripps Institute.
for about 10% of the total interaction energy.16,17 Another interaction not included in the effective two-body potentials is charge transfer (CT), which is intermolecular charge flow. Estimates of this energy contribution vary between 2018-20 and 60%17,21-23 of the total interaction energy in small hydrogenbonded clusters. In addition to their energetic importance, both polarization and CT significantly impact the charge distribution of the system of interest.21,22,24-26 For example, in a water cluster polarization increases the average charge of the oxygen atoms by 0.06 electrons through intramolecular charge flow.17 CT has less impact on the charge distribution of the system than does polarization. In most hydrogen-bonded clusters, between 0.01 and 0.1 electron is transferred between the donor and acceptor21,22,24,26 and only the charges of the atoms directly involved in the hydrogen bond are affected by CT.24,25 Recently, the importance of polarization and CT was shown for large biomolecular systems. Polarization and CT drastically altered the charge distribution of solvated major cold shock protein A, resulting in the net transfer of two electrons from the protein to the first solvation layer.27,28 Moreover, both polarization and CT were shown to be crucial in stabilization of the system.28 Despite the absence of explicit polarization and CT interactions, effective two-body potentials have been very successful in the theoretical description of biomolecular systems.11,12 Because effective potentials are fitted to experiment, correct structural and thermodynamic data are generally obtained. This success suggests that for many systems, the effect of polarization and CT can be included in an averaged way through the proper choice of parameters.13,29 The possibility that polarization and CT can be treated in an averaged way indicates that polarization and CT are either constant, or that fluctuations in these effects average to zero. However, quantum mechanical calculations show that polarization and CT are very sensitive to the environment (i.e., atom types, and bond distances and angles).17,19-22 This suggests that polarization and CT may behave quite differently for systems that undergo large conformational changes (e.g., in protein folding studies), as opposed to systems
10.1021/jp001193f CCC: $19.00 © 2000 American Chemical Society Published on Web 09/14/2000
Many-Body Effects and Protein Folding in which the conformation stays roughly the same (e.g., solvated species at equilibrium). Investigation of the role of polarization and CT on the folding of proteins would require an all quantum treatment of the folding pathway. Given the large computational expense of quantum mechanical calculations, this is clearly not feasible. Therefore, we will follow a different approach in which we will perform fully quantum mechanical calculations on snapshots from a classical folding simulation. In this way, we can quantify the instantaneous effects of polarization and CT on the energetics and charge distribution of the system. However, because the folding trajectory was obtained from an effective two-body potential, which excludes polarization and CT, the effect of polarization and CT on the protein folding pathway cannot be directly assessed. These effects can only be inferred from comparison of the charge distribution and energetics of the protein at different stages of the folding process. We will study the effect of polarization and CT on the energetics and charge distribution of the folding of betanova.30 Betanova is a small hydrophilic protein (20 residues), the native state of which is formed by three β-sheets.30 We will investigate fluctuations in charge during the protein folding process and link these fluctuations to changes in the environment. The change in environment during the protein folding process consists of the change in the packing of water around the protein and the change in conformation of the protein. Therefore, we will first analyze the effect of water on the charge distribution of the protein. We will decompose the effect of solvation into CT and polarization effects and discuss in detail how CT and polarization affect the charge distribution of the system. Then, we will analyze the intramolecular charge rearrangement (CR) of betanova. Intramolecular CR is the flow of charge through space from one residue of the protein to another, which is very sensitive to the conformation of the protein. With these studies we will be able to assess the importance of many-body effects on the folding of betanova and the quality of effective twobody potentials for protein folding simulations. Methods Snapshots of the protein folding process were obtained from a molecular dynamics simulation of betanova in 2831 TIP3P31 water molecules using the CHARMM32 force field. Details of the simulation are described elsewhere.10 Five snapshots were selected from the native state, and from the partially unfolded states in which 55, 45, 30, 15, and 10% of the number of native contacts were retained (see Figure 1). Only the native state contained three fully formed β-strands. Snapshots with 30 and 45% of the native contacts contained β1 and β2, while β2 and β3 were present in the snapshots with 55% of the native contacts. The native contacts were obtained from an analysis of two native molecular dynamics (MD) runs10 and were defined as contacts between residues with a probability higher than 56%, or with hydrogen bonds with a probability higher than 66%.10 Missing hydrogen atoms were added using the protein database within the AMBER package.33 The total number of atoms was 8804 for all snapshots and the formal charge of betanova was +3 electrons. All calculations were performed with the linear scaling Divide and Conquer (D&C) algorithm,34-36 using the AM137 and PM338,39 Hamiltonians as implemented in the DivCon 99 program.40 The D&C method divides the total system into subsystems that overlap.34-36 Subsystems were generated using the cluster-based subsetting protocol of Dixon36 with one water molecule or amino acid per core, an inner buffer layer of 4.5 Å, and an outer buffer layer of 2.0 Å. A cutoff of 8.0 Å was
J. Phys. Chem. B, Vol. 104, No. 40, 2000 9555
Figure 1. Structure of betanova. The black areas indicate the regions capable of β-sheet formation. For clarity, only the side chains of charged groups are shown. The first residue is positioned at the bottom of each structure. (A) Ten percent of native contacts; (B) 15% of native contacts; (C) 30% of native contacts; (D) 45% of native contacts; (E) 55% of native contacts; (F) native state.
applied for the off-diagonal elements for the Fock, 1-electron, and density matrixes. All reported charges are CM2.41 Interaction energies between protein and water were decomposed using the Divide and Conquer interaction energy decomposition,17 to yield electrostatic, polarization and charge-transfer components:
Eint ) Ees + Epol + ECT
(1)
Ees ) E[(Fp,Fw),r,P∞,(p:p; w:w)] E[(Fp,Fw),∞,P∞,(p:p; w:w)] (2) Epol ) E[(Fp,Fw),r,Pr,(p:p;w:w)] E[(Fp,Fw),r,P∞,(p:p; w:w)] (3) ECT ) E[F,r,Pr′,(p:pw; w:pw)] E[(Fp,Fw),r,Pr,(p:p; w:w)] (4) Here F is the Fermi energy, r indicates the protein-solvent separation, and Pr the density matrix at this separation. Subsetting is indicated by (p:p; w:w) if no overlap between protein and water subsystems exists, (p:pw; w:pw) indicates that overlap between protein and water subsystems is allowed.28 To avoid artificial polarization, P∞ was constructed from the density matrix of the protein in a vacuum and the water system in which the protein was replaced by a continuum of dielectric constant 80.28,42 The effect of polarization on the charge distribution of the system was obtained from Pr - P∞, the effect of charge transfer from Pr′ - Pr. These charge effects will be reported as negative numbers when electrons were gained and positive
9556 J. Phys. Chem. B, Vol. 104, No. 40, 2000
van der Vaart et al.
numbers when electrons were lost. Polarization of water in the field of the protein was calculated by keeping the density matrix elements for the protein at their vacuum values during the calculation of E[(Fp,Fw),r,Pr,(p:p; w:w)]. A similar procedure was used to calculate the polarization of betanova in the field of the water molecules. The intramolecular charge rearrangement was obtained by allowing each protein subsystem (i.e., amino acid residues) to only include atoms as far as 2 residues along the N and C-terminal directions of the amino acid chain, followed by complete relaxation of the system:
Eintra (vac) ) E[Fp,(p:p)] - E[Fp,(pi:pi, pi(1, pi(2)] (5) where pi indicates the ith amino acid of the protein. We chose the value of 2, because we found that for an extended linear chain of betanova (or any other polypeptide in the extended conformation), E[Fp,(pi:pi, ..., pi(n)] ) E[Fp,(p:p)] for n g 2. E[Fp,(pi:pi, pi(1, pi(2)] is the energy of a polarizable chain in which charge can only flow indirectly through the backbone, but not through space. Relaxation of the polarizable chain in solution yields intramolecular and intermolecular charge flow. The intramolecular CR can be obtained by subtracting the intermolecular charge-transfer energy (ECT) from the relaxation energy of the polarizable chain in solution:
Eintra (sol) ) (E[F,r,Pr′,(p:pw; w:pw)] E[(Fp, Fw),r,Pr′′,(pi:pi, pi(1, pi(2; w:w)]) - ECT ) E[(Fp,Fw),r,Pr,(p:p; w:w)] E[(F , F ),r,Pr′′,(pi:pi, pi(1, pi(2; w:w)] (6) p
w
Although AM1 and PM3 show the same trends in all calculated properties, PM3 generally shows larger atomic charge effects and associated energies. This is probably caused by the difference in hydrogen bonding properties for AM1 and PM3.43,44 Because we found that AM1 is a better predictor of the amount of charge transferred in small hydrogen-bonded clusters,26 we will focus here on the AM1 results. The reader is referred to the Supporting Information for a listing of the PM3 results. Results Analysis of the atomic charges of betanova revealed that all charges could be grouped into a total of 15 classes based on the identity and average charge of the atom (Figure 2). We observed six carbon, two hydrogen, four nitrogen, and three different types of oxygen atoms. The average charge of these classes is solely dependent on the atom types that directly bond with the atom of interest, and the bond order. For example, the average charge of a carbon atom bonded to a hydroxyl group is 0.1 electron, while the average charge of a carbonyl carbon is 0.5 electron. Proximity to electron withdrawing groups also influences the atomic charges, e.g., the charge of the oxygen atom in the Tyr hydroxyl group is 0.05 electron more negative than in the Ser and Thr hydroxyl groups. More importantly, the charges are very similar for all folding stages. For both the vacuum and the solvated system, the average and fluctuation in charge of a certain folding stage nearly equals those of any other folding stage. This means that the atomic charges do not depend on the folding stage, and are basically transferable from the one folding stage to another. It also indicates that the combined effect of polarization and CT, which is observed as fluctuations in atomic charges, is similar for all folding stages. These fluctuations in charge are between 0.03 and 0.1 electron
Figure 2. Charge fluctuations in betanova. Shown are all atomic AM1/ CM2 charges, for all snapshots. The atomic charges are grouped into 15 classes, indicated by the letters a-o. These classes are (with the atom of interest underlined): (a) C)C (Cγ Trp); (b) CH (all carbon atoms not included in a-g); (c) C-NH3+; (d) C-OH, C-NH (Arg, Trp, and all CR); (e) COO-, C)O (Gln, Asn, and all peptide CdO); (f) -C+-(NH2)2; (g) CH (all hydrogen atoms not included in h); (h) OH, NH3+, NH2+, NH2, NH (Trp, Arg, and all peptide NsH); (i) NH2, NH Arg; (j) NH2+, NH peptide; (k) NH3+; (l) NH Trp; (m) COO-; (n) C)O (Gln, Asn, and all peptide (CdO), OH (Ser, Thr); (o) OH Tyr. (A) C and H charges for the protein in a vacuum. (B) C and H charges for the protein in solution. (C) N and O charges for the protein in a vacuum. (D) N and O charges for the protein in solution.
per atom for each folding stage. However, the amino acids in betanova are solvent exposed, so it is not possible to definitively state that this observation (i.e., atom type charges are a constant at all folding stages) is transferable to larger globular proteins. Table 1 summarizes the average charge and fluctuations in charge for all the carbon, hydrogen, nitrogen, and oxygen classes. The reader is referred to the Supporting Information for a complete listing of the averages and fluctuations in charge for all stages of folding. Figure 2 and Table 1 show that solvation has a large effect on the oxygen charges. All oxygen atoms become ∼0.05e more negative upon solvation, and the charge fluctuations of the carbonyl and hydroxyl oxygen atoms decrease 6-fold. Solvation has less impact on other atoms, most nitrogen atoms become slightly less negative, hydrogen atoms gain small amounts of charge, and carbon atoms either gain or lose small amounts of charge. The net effect of solvation on the charge distribution of betanova is a transfer of negative charge from the water to the protein. Table 2 shows that the average amount of charge transferred is ∼0.20 electron, while the observed range of charge transferred is very similar for most snapshots. The amount of charge transferred varies from 0.117 to 0.287 electron (see Supporting Information), which is a considerable amount of charge, considering that betanova consists of only 20 amino acids. To further investigate the effect of solvation, we decomposed the influence of water on the charge distribution of betanova into CT and polarization effects. The polarization effects describe intramolecular charge flow (i.e., charge flow within betanova and the solvent), while CT effects describe intermolecular charge flow across the protein-water interface. In Figure 3 the effect of CT on the charge of each protein residue is plotted. The bars in this figure indicate the observed range of the charge effect for each amino acid at each folding stage and bold bars indicate the charge effect on the charged residues. Negative values indicate a gain in electrons, while positive values indicate a loss of electrons. Figure 3 shows that the
Many-Body Effects and Protein Folding
J. Phys. Chem. B, Vol. 104, No. 40, 2000 9557
TABLE 1: Average Atomic Chargesa atom
figure
vacuum
solution
atom
figure
vacuum
solution
C1b C2c C3d C4e C5f C6g H1h H2i
2a 2b 2c 2d 2e 2f 2g 2h
-0.211 ( 0.039 -0.129 ( 0.020 -0.019 ( 0.011 0.102 ( 0.001 0.507 ( 0.000 0.755 ( 0.003 0.095 ( 0.042 0.414 ( 0.029
-0.205 ( 0.007 -0.133 ( 0.018 0.014 ( 0.005 0.101 ( 0.019 0.525 ( 0.001 0.749 ( 0.032 0.098 ( 0.006 0.423 ( 0.010
N1j N2k N3l N4m O1n O2o O3p
2i 2j 2k 2l 2m 2n 2o
-0.826 ( 0.017 -0.691 ( 0.024 -0.657 ( 0.014 -0.532 ( 0.014 -0.646 ( 0.008 -0.505 ( 0.033 -0.457 ( 0.007
-0.814 ( 0.017 -0.690 ( 0.025 -0.647 ( 0.005 -0.538 ( 0.037 -0.679 ( 0.013 -0.553 ( 0.005 -0.500 ( 0.010
a AM1/CM2, in electron. Shown are the averages and standard deviations over all snapshots. b C)C (Cγ Trp). c CH (all carbon atoms not included in C1 and C3-C6). d C-NH3+. e C-OH, C-NH (Arg, Trp, and all CR). f COO-, C)O (Gln, Asn, and all peptide). g -C+-(NH2)2. h CH (all hydrogen atoms not included in H2). i OH, NH3+, NH2+, NH2, NH (Trp, Arg, and all peptide). j NH2, NH Arg. k NH2+, NH peptide. l NH3+. m NH Trp. n COO-. o C)O (Gln, Asn, and all peptide), OH (Ser, Thr). p OH Tyr.
TABLE 2: Net Amount of Charge Transferred at the Protein-Water Interface N (%)a
∆qb
10 15 30 45 55 100 total
0.228 ( 0.013 0.231 ( 0.066 0.209 ( 0.050 0.232 ( 0.061 0.220 ( 0.023 0.194 ( 0.050 0.219 ( 0.046
a Fraction of native contacts. b Average and standard deviation for AM1/CM2, in electron. Betanova accepts a net negative charge.
Figure 3. The effect of charge transfer on the charge distribution of betanova. These effects are induced by the solvent. Charges are listed by residue. The bars indicate the observed range of charge effects, bold bars indicate the charged residues. Negative values indicate a gain of electrons, positive values indicate a loss of electrons. Shown are CM2 charges, obtained using the AM1 Hamiltonian. (A) 10% of native contacts; (B) 15% of native contacts; (C) 30% of native contacts; (D) 45% of native contacts; (E) 55% of native contacts; (F) native state.
charged residues are responsible for most of the charge transferred from water to protein. Positively charged groups gain electrons and negatively charged groups lose electrons upon CT. Table 2 and Figure 3 suggest that CT is relatively independent of the local environment, since the total amounts of CT and CT per residue are very similar for all frames. At the atomic level, we found that the hydrogen atoms of the ammonium groups of Lys 9, 15 and Arg 1, the hydrogen atoms of the guanidinium groups of Arg 1 and 20, and the oxygen atoms of the carboxylate groups of Glu 18 and Arg 20 had the largest CT effects. Indeed, the CT effect on these atoms was about an order of magnitude larger than the CT effect of the other atoms of betanova. The average CT effect of each
ammonium group is -0.02 ( 0.005e and for each guanidinium group -0.013 ( 0.004e per hydrogen. The average CT effect for each carboxylate is +0.018 ( 0.005e per oxygen. The fluctuations, here given as standard deviations, are caused by differences in hydrogen bonding of these atoms with the water layer. Hydrogen bonds with short bond lengths enhance CT, while longer bond lengths reduce CT. The CT effect on Arg 1 and Arg 20 indicates that the CT of each “functional” group is to a great extent independent of the CT effect on other groups. In other words, the CT effect on Arg 1 is roughly equal to the sum of the CT on the ammonium group and guanidinium group, while the CT on Arg 20 is roughly equal to the sum of the CT on the guanidinium group and carboxylate group. The average difference between the total CT on the residue and the sum of the CT on the guanidinium, ammonium, and carboxylate groups is -0.03 ( 0.005e for Arg 1 and -0.01 ( 0.006e for Arg 20. Arg 20 is an interesting residue, because the guanidinium group gains electrons, while the carboxylate group loses electrons upon CT. On average, a gain of electrons is to be expected, because there are 4 guanidinium hydrogen atoms and only 2 carboxylate oxygen atoms. This is exactly what was observed for most frames. There are only two frames in which the total CT on Arg 20 is positive: one frame of the native state where CT is 0.046 electron and one frame with 45% of the native contacts where CT is 0.001 electron. Analysis showed that this is due to bad hydrogen bonding contacts between water molecules and the guanidinium groups: in these two frames the hydrogen bonds are long and have small angles. In both cases water was prevented from hydrogen bonding by the proximity of the carboxyl group of Glu 18 to the guanidinium group of Arg 20. In the former snapshot a salt bridge was formed between Glu 18 and Arg 20 (2.0 Å apart), while in the latter snapshot a much weaker salt bridge was observed (3.9 Å apart). There are other frames in which close contacts between Glu 18 and Arg 20 exist, but in these frames the guanidinium group is still able to hydrogen bond to surrounding water molecules. For these cases we found that the CT effect is very similar to the CT observed in frames without the Glu 18-Arg 20 contact. Analysis of the effect of CT on the charge of the water molecules showed that charge was transferred to the water molecules that were directly hydrogen bonded to the protein. Most of the water molecules that have CT effects smaller than -0.004 or larger than 0.004 electron were hydrogen bonded to the charged groups, while the remainder were hydrogen bonded to hydrophilic groups or the amide and carbonyl groups of the backbone. Water molecules that were hydrogen bonded to negatively charged groups or the carbonyl backbone gain electrons, while water molecules that were hydrogen bonded to positively charged groups or the amide backbone lose
9558 J. Phys. Chem. B, Vol. 104, No. 40, 2000
van der Vaart et al. TABLE 3: Charge Transfer Parametersa
a
Figure 4. Charge transfer across the protein-water interface. Shown are the amounts of charge accepted by the water molecules via hydrogen bonds with the protein. Negative values indicate a gain in electrons, while positive numbers indicate a loss in electrons for the water molecules. (A) CT in hydrogen bonds of water with hydrogen atoms of the protein. (B) CT in hydrogen bonds of water with oxygen atoms of the protein, excluding carboxylate oxygen atoms. (C) CT in hydrogen bonds of water with carboxylate oxygen atoms of the protein.
electrons. Water molecules that were hydrogen bonded to the oxygen atoms of the Ser, Asn, Gln, Thr, and Tyr side chains gain electrons, while water molecules that were hydrogen bonded to the hydrogen atoms of these side chains lose electrons upon CT. Charge transferred from the protein to water stays localized in the first solvation shell and does not smear out over the entire water network. The observation that CT only occurs between hydrogen bond donors and acceptors (i.e., CT is strongly dependent on simple geometrical factors and it is largely independent of CT of other groups) suggests that CT can be described by a simple additive model. To investigate this notion, we plotted the amount of charge transferred to water atoms in hydrogen bonds with the protein as function of the hydrogen bond length and angle (see Figure 4). To prevent any ambiguity, we only included hydrogen bonds for which the water molecule was not involved in any other hydrogen bonds with the protein. If a water molecule was hydrogen bonded to more than one protein atom, then all hydrogen bonds involving that particular protein atom were also skipped. In Figure 4 the data for all these hydrogen bonds are shown, for all snapshots at once. This figure clearly shows that CT is a simple exponential function of the hydrogen bond distance. Surprisingly, CT is not dependent on the hydrogen bond angle. All hydrogen bonds could be grouped into three classes, based on the rate of decay of CT. These classes are hydrogen bonds of the solvent with protein hydrogen atoms (Figure 4a), with carboxylate oxygen atoms (Figure 4c), and with all other oxygen atoms (Figure 4b). The amount of charge transferred could be described by: H-bonds
∆qwrp ) -
∑p
H-bonds
∆qprw )
∑ w
Ape-Bprpw
Ape-Bprpw
CT to water
CT to protein
(7)
where p indicates a protein atom involved in a hydrogen bond with the solvent, w a water molecule involved in a hydrogen bond with the protein, and rpw is the hydrogen bond distance. Ap and Bp are coefficients that only depend on the class of hydrogen bond (Table 3). This model worked very well to describe all CT effects in the betanova snapshots. The standard deviation between the predicted and observed amount of CT over all betanova atoms averaged over all snapshots was 0.0009e, and 0.0032e if only the atoms involved in hydrogen bonding with the solvent were considered. It is important to note that only one class of protein hydrogen atoms was needed, even though the amount of charge transferred in hydrogen bonds between the solvent and am-
protein atom type
A (e)
B (Å-1)
hydrogen carboxylate oxygen other oxygen
-20.5151 0.5785 4.1024
4.0217 2.3987 3.7038
Obtained by simulated annealing.
monium or guanidinium groups is very different from hydrogen bonds with hydroxyl or peptide groups. The main reason for this involves the observation that hydrogen bonds with ammonium and guanidinium groups are generally shorter, and more hydrogen bonds are formed between the solvent and these groups than between other groups involving hydrogen atoms. This means that the distinction in CT of hydrogen atoms follows from differences in hydrogen bonding preference. The observation that the total amount of CT and the amount of CT per residue stay fairly constant can then be explained by the observation of low fluctuations in hydrogen bonding geometry. All residues, especially those of the charged groups (with large CT effects) have a significant amount of exposure to the surrounding solvent during the folding process and hydrogen bond lengths with the solvent show only small fluctuations. Thus, most of the groups that undergo CT are solvent exposed and the observation of fairly constant CT along the unfolding coordinate is not too surprising. Future studies on proteins where there are a significant number of buried residues may show larger fluctuations in CT. When applied to another protein, the model did not work so well. For major cold shock protein A (CspA),27,28 the mean standard deviation over all protein atoms was 0.0029e, and 0.0088e for atoms involved in hydrogen bonds with the solvent only. This was mainly caused by the low quality of A and B parameters for carboxylate oxygen atoms (the mean standard deviation of the carboxylate oxygen atoms was 0.0318e; for the hydrogen atoms, it was 0.0035e). Figure 4 shows that this is due to a lack of carboxylate oxygen data: most water molecules that hydrogen bonded to these atoms also bonded to other protein atoms and were subsequently excluded from the fit. Because CspA contains 9 carboxylate groups (betanova contains only 2), deficiencies in the carboxylate fit are more pronounced in the CspA system. From these data is clear that much more research is needed for an accurate assessment of the CT parameters. However, for the argument pursued in this article, the functional form of CT is more important than the actual value of the parameters. Upon solvation, polarization affects the charge distribution of betanova much more than CT. This is visualized in Figure 5 where the effect of polarization on the charge distribution of betanova is shown for all residues. Figure 5 shows that negatively charged residues generally gain electrons, while positively charged residues generally lose electrons upon polarization. However, polarization was not limited to the charged residues. Polarization affects almost all residues and the effect varies substantially during the folding process. The large fluctuation of the polarization effect clearly indicates the strong dependence of polarization on the local environment. To gain a further understanding of polarization, we analyzed the effect of polarization on the atomic charge distribution. As Figure 5 suggested, we found that many atoms are highly polarizable (g0.02 electron). In general, all oxygen atoms, the hydroxyl- and nitrogen-bound hydrogen atoms, and the carbonyl carbon atoms show very high polarization effects. Oxygen polarization effects are ∼-0.05 electron, carbonyl carbon effects are ∼0.02 electron, and hydroxyl and nitrogen bound hydrogen
Many-Body Effects and Protein Folding
J. Phys. Chem. B, Vol. 104, No. 40, 2000 9559
Figure 5. The effect of polarization on the charge distribution of betanova. These effects are induced by the solvent. Charges are listed by residue. The bars indicate the observed range of charge effects, bold bars indicate the charged residues. Negative values indicate a gain of electrons; positive values indicate a loss of electrons. Shown are CM2 charges, obtained using the AM1 Hamiltonian. (A) 10 percent of native contacts; (B) 15% of native contacts; (C) 30% of native contacts; (D) 45% of native contacts; (E) 55% of native contacts; (F) native state.
effects are ∼0.02 electron. The polarization effect for hydroxyl hydrogens is larger than for other hydrogen atoms (∼0.025 electron); for the oxygen atoms the oxygen carbonyls of the Gln and Asn side chains show the largest polarization effect (∼-0.07 electron). Fluctuations in the polarization effect are fairly large and in many cases fluctuations are of the same order of magnitude as the average (Figure 6). However, we found that for polarization, large average charge effects do not necessarily mean large fluctuations. This behavior is entirely different from the CT effect, where fluctuations are much smaller and the averages and fluctuations are closely related (Figure 6). This observation is another indication that polarization is a many-body effect and cannot simply be described in terms of atom types and hydrogen bonding qualities. A closer look at the polarization effect on the amide groups of Asn and Gln, and the CβH2 group of the Arg residues, indicated that the polarization of the hydrogen atoms attached to these groups is highly asymmetrical (Table 4). In many cases, polarization of one of the hydrogen atoms contained in these groups was an order of magnitude larger than for the other hydrogen atom. Sometimes, the polarization effect was of an opposite sign for the pair of hydrogen atoms. The difference in this behavior can easily be explained by the differences in the local chemical environment. Analysis showed that water molecules close to negatively polarized solute atoms will be
Figure 6. Charge fluctuations upon solvation of betanova. Shown are the average charge effect 〈q〉 and standard deviation σq for all atoms. Averages are taken over all 30 snapshots, calculated by AM1/CM2. Negative values indicate the gain of electrons; positive values indicate the loss of electrons. Atoms with large charge effects are printed in bold. (A) Charge-transfer effect; (B) polarization effect.
positively polarized, while water molecules will be negatively polarized, when close to positively polarized solute atoms. Upon solvation, water molecules close to the carbonyl oxygen atom of the amide groups will be positively polarized. The hydrogen atom on the same side of the carbonyl oxygen atom will be close to the positively polarized water molecules and will therefore tend to be negatively polarized. Polarization of the hydrogen atom on the same side as the carbonyl oxygen atom will therefore be smaller (less positive or even negative) than polarization of the other hydrogen atom. This was indeed the case in most frames for the Asn and Gln amide side chains. However, proximity of the hydrogen atom to other groups can change this picture. For instance, in several frames the hydrogen atom on the same side of the carbonyl oxygen atom was close to the guanidinium group of Arg 1. In that case the water molecules between the amide hydrogen and guanidinium group will be negatively polarized and the hydrogen atom will be positively polarized. The positive polarization of the hydrogen atom and the proximity of negatively polarized water molecules
TABLE 4: Asymmetrical Polarization Effectsa residue Gln 6 Asn 7 Asn 12 Asn 13
atomb CONHH CONHH CONHH CONHH CONHH CONHH CONHH CONHH
∆qpolc 0.0216 ( 0.009 0.0060 ( 0.007 0.0206 ( 0.008 0.0082 ( 0.006 0.0150 ( 0.012 0.0082 ( 0.011 0.0243 ( 0.008 0.0116 ( 0.008
residue
atomb
∆qpolc
Arg 1
H-Cβ-H
0.0004 ( 0.020 -0.0257 ( 0.016 0.0198 ( 0.027 -0.009 ( 0.016
Arg 20
H-Cβ-H H-Cβ-H H-Cβ-H
a Shown are the polarization effects for two identical atoms of some groups. Note that a positive polarization effect means that electrons are lost upon polarization, a negative polarization effect indicates the gain of electrons. b Polarization effect of the underlined atom. c AM1/CM2 charge effect, in electron. Shown are the averages and standard deviations of the particular atom over all frames.
9560 J. Phys. Chem. B, Vol. 104, No. 40, 2000
van der Vaart et al.
Figure 7. Interactions at the betanova-water interface using the AM1 Hamiltonian. N indicates the percentage of native contacts present. (A) Interaction energies and decomposition of the interaction energies into electrostatic, polarization, and charge-transfer components. Shown are the energies for all snapshots at each stage of the protein folding process. (B) Average interaction energies.
will also have some influence on the carbonyl oxygen atom of the amide group. Indeed, we found that this oxygen atom was less negatively polarized in those cases. For the β-hydrogen atoms of Arg 1, we found that the hydrogen with the largest, positive polarization effects were in close proximity to the ammonium group of Arg 1. For Arg 20, the sign of polarization was mostly controlled by proximity to the carboxyl group or guanidinium group of Arg 20. The very simple approach described above worked well to explain other polarization effects at the betanova-water interface as well. The sign of polarization was generally determined by chemical identity, while the magnitude was found to be controlled by the environment. Proximity to positively polarized groups makes the polarization effect of the atom of interest more positive or less negative, while proximity to negatively polarized groups makes it less positive or more negative. In this picture, water molecules were essential in communicating the charge distribution of two closely positioned groups. In classical terms, the energy of the system is minimized by enhancement of the dipoles of the water molecules and each of the closely positioned groups. Because the charge distribution of two close amino acid residues becomes either more positive or more negative, this process introduces some repulsion between the residues. However, this repulsion is more or less neutralized by the screening of the water molecules. After establishing the effect on the charge distribution of the system, we will now turn to the energetic effects of CT and polarization. In Figure 7A the electrostatic, polarization, CT, and interaction energies between the protein and water molecules are plotted for all snapshots. On the x-axis, the snapshots are listed by the number of native contacts (N - 5 snapshots per folding stage). The vertical dotted lines mark the boundaries between the folding stages. The averages are graphed for each of the different states in Figure 7B. These figures show that CT is the main contributor to the interaction energy. On average CT accounts for 75% of the interaction energy, polarization 15%, and electrostatics 10% (see Supporting Information). CT and polarization were always stabilizing and electrostatics were only slightly destabilizing in two snapshots where the protein retains 55% of the native contacts. Polarization was a constant contributor to the interaction energy, while CT slightly decreases during the protein folding process. Hence, the interaction energy curve shows roughly the same shape as the electrostatic energy curve. This means that the interaction between protein and water can be approximated by the electrostatic energy augmented by a nearly constant bias due to polarization and CT during the folding process. The importance of the fluctuations in the electrostatic energy certainly warrant a detailed investigation into the reasons these
Figure 8. Polarization at the protein-water interface. Shown are the polarization of betanova in the field of water (Epol(p)), the polarization of water in the field of betanova (Epol(w)), the sum of Epol(p) and Epol(w), and the total polarization energy Epol. These energies are listed for all snapshots at each stage of the protein fold. (A) Polarization energies; (B) polarization energies as percentage of the total polarization energy.
fluctuations occur. However, the complexity of the proteinwater interface severely impedes efforts along these lines. For example, no obvious structural reasons were found for the increase in electrostatic repulsion at the 55% folding stage. Although this stage is the only stage in which the β2 and β3 sheets are formed at the same time, we did not find close contacts between charged groups which would have increased the electrostatic energy. Moreover, the structures of the third and fourth snapshot of the 55% folding stage are remarkably similar, even though the difference in electrostatic energy is ∼170 kcal/mol. A detailed comparison of all the individual Coulomb integrals along the folding trajectory would certainly shed light on this important issue, but the sheer volume of integrals prevented us from doing such an analysis at this point. Polarization of the protein is more or less independent of the polarization of the water molecules, and vice versa. This can be seen from Figure 8, where the polarization energy of the protein in the field of the water and the polarization energy of water in the field of the protein are plotted for all snapshots. These energies were obtained by optimizing the charge distribution of one of the systems (water or protein), while maintaining the charge distribution of the other system at its vacuum value. This is different from the total polarization energy where the charge distribution of the entire system is optimized. Figure 8 shows that the sum of the two in-field polarization energies follows the total polarization energy closely. Indeed, the sum of the in-field polarization energies amounts to ∼90% of the total polarization energy for all frames. This means that only ∼10 kcal/mol was gained by coupling the protein and water systems. The charge distribution of the in-field polarized systems also closely resembles the charge distribution of the fully polarized system. For the protein, the maximum increase of charge per atom on the fully versus the in-field polarized system was only ∼0.01 electron. Not surprisingly, atoms that have high polarization effects see a greater increase of charge than atoms that have smaller polarization effects. Overall, the maximum increase of charge was ∼10% per atom. For water, the change in charge was even smaller: maximally 0.001 electron per water molecule. Moreover, only the first solvation layer was affected. Note that for both water and the protein, polarization effects were generally enhanced in going from in-field to full polarization. This means that negative polarization effects become more negative and positive polarization effects more positive.
Many-Body Effects and Protein Folding
Figure 9. Intramolecular charge transfer. (A) Intramolecular chargetransfer energies for the protein in solution (Eintra(sol)) and vacuum (Eintra(vac)) and the difference between these terms (∆Eintra). (B) Correlation between the number of close contacts and the intramolecular charge transfer energy in solution. Only contacts between atoms of residues at least 2 residues upstream or downstream are included, contacts are e 2.5 Å.
We will now discuss the influence of the conformation of the protein on the charge distribution of the system, by looking at the intramolecular CR. Figure 9A shows that the intramolecular CR in solution and vacuum are nearly identical. The maximum difference between the two is 0.9 kcal/mol, while for most frames the difference is 0.0 kcal/mol. Intramolecular CR accounts for between 0 and 19% of the total energy of the protein in a vacuum (see Supporting Information). We note that for PM3 this percentage is between 0 and 28%, because the intramolecular CR of PM3 is twice as large as that of AM1. Although intramolecular CR stabilizes the protein significantly in a vacuum, comparison with the protein-water interaction energies (Figure 7) shows that it is not as important in solution. In solution, charge transfer between protein and water is at least 25 times more stabilizing than intramolecular CR. Intramolecular CR is a very localized and short-range effect. This can be seen from Figure 9B where intramolecular CR energies are plotted versus the number of protein-protein contacts shorter than 2.5 Å. The observed correlation clearly suggests that intramolecular CR only affects atoms that are very close in space. Analysis of the effect of intramolecular CR on the charge distribution of the system supports this observation. Figure 10 shows the effect of intramolecular CR on the charge distribution of the protein for a snapshot at which the protein has 30% of the native contacts. In this figure, the effect of intramolecular CR is plotted for the charge of each protein atom. A negative change in charge indicates a gain of electrons; a positive sign indicates a loss of electrons. Figure 10 shows that the effect of intramolecular CR is very similar in a vacuum (Figure 10A) and solution (Figure 10B). Although large intramolecular CR charge effects were observed for the same atoms in a vacuum and solution, the magnitude of the effect differs slightly, usually between 0 and 10% (but larger deviations up to ∼40% were observed for some hydrogen bonds, see below). Further analysis indicated that large intramolecular CR effects occur exclusively for atoms that are close in space. Moreover,
J. Phys. Chem. B, Vol. 104, No. 40, 2000 9561
Figure 10. Intramolecular charge rearrangement for a snapshot with 30% of the native contacts. Hydrogen bonds are indicated by the letters a-c, other contacts are indicated by d-g, and the atom of interest is underlined according to: (a) 33 CONH Gly 2-195 CONH Asn 13 (2.09 Å, 148.6°); (b) 45 NH Trp 3-277 COO- Glu 18 (2.06 Å, 169.4°); (c) 59 CONH Ser 4-179 CONH Thr 11 (2.00 Å, 136.2°); (d) 14 Hδ Arg 1-82 Hγ Val 5 (1.67 Å); (e) 37 HR Trp 3-183 HR Asn 12 (1.98 Å); (f) 43 H (aromatic) Trp 3-151 Hβ Tyr 10 (1.84 Å); (g) 48 H (aromatic) Trp 3-271 Hβ Glu 18 (1.90 Å). (A) Vacuum; (B) solution.
all close contacts were observed as peaks in the intramolecular CR charge plots. In Figure 10 all peaks have been assigned to close contacts through the use of the letters a-g. Contacts a-c are hydrogen bonds and these are assigned in Figure 10A. All other contacts are H-H contacts and are assigned in Figure 10B. Figure 10 shows that intramolecular charge transfer only occurs through hydrogen bonds, while other close contacts are merely polarized. Hydrogen bond donor and acceptor atoms have opposite signs for the charge effect (Figure 10A), while other atoms that are close to each other have identical signs for the charge effect (Figure 10B). In hydrogen bonds, the hydrogen donor gains electrons, while the hydrogen acceptor loses electrons. Note that this is also observed for the protein-water hydrogen bonds. The effective amount of charge transferred in the intramolecular hydrogen bonds can be deduced from the peak heights in the intramolecular CR plots. The amount of charge transferred depends on the quality of the hydrogen bond: shorter bond lengths and large bond angles promote charge transfer, whereas longer bond lengths and small bond angles decrease the amount of charge transferred. In the O‚‚‚H hydrogen bond between the peptide side chains of Ser 4 and Thr 11, only ∼0.001e is transferred (Figure 10A, pair c), whereas in the O‚‚‚H hydrogen bond between the backbone amide of Trp 3 and the carboxylate of Glu 18, ∼0.009e is transferred (Figure 10A, pair b). Although the hydrogen bond lengths are very similar (2.00 and 2.06 Å, respectively), the bond angle for the Ser 4‚‚‚Thr 11 pair is 136°, whereas the bond angle for the Trp 3‚‚‚Glu 18 pair is 169°. In solution, the amount of charge transferred in the Trp 3‚‚‚Glu 18 pair decreases to ∼0.005e, because upon solvation a water molecule hydrogen bonds to the Glu 18 carboxylate group (with a bond length of 1.95 Å and an angle of 169°, this water molecule accepts -0.007e). The Gly 2‚‚‚Asn 13 and Ser 4‚‚‚Thr 11 hydrogen bonds are not affected by solvation, because no water molecules directly hydrogen bond to these groups. The contact distance and charge effect are also correlated for the non-hydrogen bond close contacts. However, charge effects not only occur for the atoms directly involved in the close
9562 J. Phys. Chem. B, Vol. 104, No. 40, 2000 contact (as in hydrogen bonds), but also for atoms where the close contacts are covalently bonded. In Figure 10B all atoms directly preceding the atoms involved in the close contact (i.e., atom 13 and 81 for contact d, atom 36 and 182 for contact e, etc.) show charge effects about equal in magnitude, but opposite in sign of the close contact atom. These are the charge effects on the carbon atoms to which the hydrogen atoms of the close contacts are covalently bound. Because the sign of the charge effect was identical for two non-hydrogen bond atoms that were in close contact, no charge flows between these atoms. These atoms were merely polarized, meaning that they either both gain or both lose electrons. They obtain this charge from the atom to which they covalently bond, causing the charge effect of this covalently bound atom to be opposite from that of the atom involved in the close contact. This observation makes sense because by altering the charge distribution in a non-hydrogen bond close contact, the electron-electron repulsion is reduced, which thereby reduces the destabilizing tendency of the close contact. Further analysis showed that similar results were obtained for all other frames as well. Because large charge effects were only observed for atoms that were in close contact with each other, water generally does not have a large effect on the observed charges. Only when water directly hydrogen bonds to the atoms involved were large deviations observed (10% or up). For most close contacts, this is a rare event. This explains why the observed energies of intramolecular CR in a vacuum and solution were observed to be so similar (Figure 9A). Discussion Our calculations showed the importance of many-body and CT effects on the energetics and charge distribution of betanova during the protein folding process. Polarization and CT are the main stabilizing factors in solvation of the protein (respectively 15.5 and 74.1% of the total interaction energy, see Supporting Information), while electrostatics are less important. CT significantly impacts the charge of the charged groups, while polarization significantly affects the charge of the oxygen atoms, carbonyl carbon- and nitrogen-bound atoms, and hydroxyl hydrogen atoms. Total fluctuations in charge are up to 0.1 electron per atom. These findings suggest that polarization and CT should be included in a proper description of the protein folding process. However, we also found that both the averages and fluctuations in charge are very similar for all stages of folding. The average charge distribution is not dependent on the overall fold of the protein, which makes the average atomic charges (and also the fluctuations in charge) portable from one folding stage to the other. The fundamental reason for this observation is the local character of CT and polarization. Our detailed analysis showed that CT only occurs in hydrogen bonds and can be described by a simple exponential function of the hydrogen bonding distance. Polarization was shown to affect almost all atoms (hydrogen-bonded and non-hydrogen-bonded atoms), but polarization effects could be rationalized by considering the chemical identity of the atom and changes in the local environment (proximity of negatively and positively polarized atoms, etc.). For betanova, which is a hydrophilic protein without a large solvent-excluded core, changes in local environment are very similar at each folding stage. Therefore, fluctuations in charge will also be similar at each folding stage. The similarity of the average charge and fluctuations in charge during the protein fold suggest that the charge distribution of betanova may be approximated by the average charge distribu-
van der Vaart et al. tion. This approximation becomes increasingly more accurate as the number of configurations is increased. Application of this approximation should therefore be limited to averaged properties; more instantaneous properties will require a fluctuating charge model. Moreover, our results also suggests that it may be straightforward to classify all atoms in proteins according to the classes given in Figure 2. However, as noted earlier, to generalize this conclusion will require the examination of systems which have clearly defined core regions. It is our expectation that there will be significant differences between buried and exposed atoms, but this has still to be verified. In terms of the energy, we observed that both polarization and CT stay roughly constant during the folding process. Although CT governs the interactions between protein and water, the evolution of protein-water interactions roughly follows the change in electrostatic interactions upon folding. Therefore, the interactions between protein and water can be approximated by the electrostatic energy, increased by a constant due to polarization and CT. This means that for classical simulations, correct description of electrostatic energy is essential, whereas polarization and CT energies may be averaged through proper parametrization. Our analysis suggests that a well-parametrized effective twobody potential is able to properly describe the folding of betanova. The energetic effects of polarization and CT can be included through proper parametrization, while the effects of fluctuating charges can be eliminated by sufficiently long sampling. We note that an effective two-body potential would suffice here, since changes in the local environment are very similar at all stages of folding. Care is needed in the description of systems with radical changes of local environment, e.g., in the folding of proteins with a large hydrophobic core, and the folding of membrane proteins. Application of a fluctuating charge model or application of multiple atomic charge sets may be necessary for these systems. We realize that more sophisticated ab initio or density functional theory (DFT) methods are necessary for a more detailed quantitative analysis of the observed charge effects. However, previous studies showed that semiempirical models correctly predict the amount of charge transferred in small hydrogen-bonded clusters.26,45 Moreover, the use of the CM2 charge model,41 which has also been consistently parametrized for ab initio and DFT calculations,41 provides us with highquality charge distributions. Hence, we believe that our calculations will be in qualitative agreement with high-level ab initio methods, while quantitative agreement on the charge distributions may have been reached through use of the CM2 model. Conclusion Our calculations showed that many-body and CT effects are very important in the folding of betanova. CT is the main contributor to the protein-water interaction energy and both polarization and CT significantly impact the charge distribution of the system. CT only affects the atoms in a hydrogen bond, whereas polarization greatly affects oxygen atoms, carbonyl carbon- and nitrogen bound-atoms, and hydroxyl hydrogen atoms. Total fluctuations in charge are up to 0.1 electron per atom. Despite the absence of explicit polarization and CT interactions, effective two-body potentials can properly describe the folding process of betanova. Both CT and polarization only depend on the local environment and not on the overall fold of the protein. Because fluctuations in local environment are similar for betanova, averages and fluctuations in charge are very similar
Many-Body Effects and Protein Folding at each stage of the folding process. Therefore, the charge distribution of betanova may be modeled by the average charge distribution, an approximation that becomes increasingly more accurate as sampling is increased. The contributions of polarization and CT to the interaction energies are nearly constant during the protein folding process. This means that protein folding is driven by electrostatic interactions, which can be properly modeled in effective two-body potentials. More accurate descriptions of the protein folding process will require fluctuating charge models that better represent the instantaneous properties of the system. This may become important for systems in which the local environment changes radically, e.g., in the folding of proteins with large hydrophobic cores. We note that classical polarizable models already exist,46-50 and work is under way to accurately model CT.51,52 Acknowledgment. We thank the National Center for Supercomputer Applications for generous allocations of supercomputer time. This research was supported by the Department of Energy (DE-FG02-96ER622670). Supporting Information Available: Average charges and fluctuations for the carbon, hydrogen, nitrogen, and oxygen classes for all folding stages. Tables with the interaction, electrostatic, polarization, charge transfer, and intramolecular charge rearrangement energies for the AM1 and PM3 Hamiltonians and a listing of the net amount of charge transferred between water molecules and the protein are included as well. This information is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Daura, X.; Jaun, B.; Seeback, D.; van Gunsteren, W. F.; Mark, A. E. J. Mol. Biol. 1998, 280, 925-932. (2) Daura, X.; Gademann, K.; Jaun, B.; Seebach, D.; van Gunsteren, W. F.; Mark, A. E. Angew. Chem., Int. Ed. Engl. 1999, 38, 236-240. (3) Daura, X.; van Gunsteren, W. F.; Mark, A. E. Proteins 1999, 34, 269-280. (4) Chipot, C.; Pohorille, A. J. Am. Chem. Soc. 1998, 120, 1191211924. (5) Bartels, C.; Schaefer, M.; Karplus, M. Theor. Chem. Acc. 1996, 101, 62-66. (6) Boczko, E. M.; Brooks, C. L., III. Science 1995, 269, 393-396. (7) Sheinerman, F. B.; Brooks, C. L., III. J. Mol. Biol. 1998, 278, 439456. (8) Sheinerman, F. B.; Brooks, C. L., III. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 1562-1567. (9) Duan, Y.; Kollman, P. A. Science 1998, 282, 740-744. (10) Bursulaya, B. D.; Brooks, C. L., III. J. Am. Chem. Soc. 1999, 121, 9947-9951. (11) Pettitt, B. M.; Smith, J. C. Comput. Phys. Comm. 1995, 91, 1-344. (12) van Gunsteren, W. F.; Berendsen, H. J. C. Angew. Chem., Int. Ed. Engl. 1990, 29, 992-1023. (13) Sprik, M. EffectiVe Pair Potentials and Beyond; Allen, M. P., Tildesley, D. J., Eds.; NATO ASI Series C., Vol 397; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1993; pp 211-259. (14) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987.
J. Phys. Chem. B, Vol. 104, No. 40, 2000 9563 (15) Milet, A.; Moszynski, R.; Wormer, P. E. S.; van der Avoird, A. J. Phys. Chem. A 1999, 103, 6811-6819. (16) Gao, J.; Xia, X. Science 1992, 258, 631-635. (17) van der Vaart, A.; Merz, K. M., Jr. J. Phys. Chem. A 1999, 103, 3321-3329. (18) Kitaura, K.; Morokuma, K. Int. J. Quantum Chem. 1976, 10, 325340. (19) Stone, A. J. Chem. Phys. Lett. 1993, 211, 101-109. (20) Stevens, W. J.; Fink, W. H. Chem. Phys. Lett. 1987, 139, 15-22. (21) Reed, A. E.; Weinhold, F.; Curtiss, L. A.; Pochatko, D. J. J. Chem. Phys. 1986, 84, 5687-5705. (22) Reed, A. E.; Curtiss, L. A.; Weinhold, F. Chem. ReV. 1988, 88, 899-926. (23) Weinhold, F. J. Mol. Struct. 1997, 399, 181-197. (24) Umeyama, H.; Morokuma, K. J. Am. Chem. Soc. 1977, 99, 13161332. (25) Morokuma, K. Acc. Chem. Res. 1977, 10, 294-300. (26) van der Vaart, A.; Merz, K. M., Jr. Int. J. Quantum Chem. 2000, 77, 27-43. (27) Nadig, G.; van Zant, L. C.; Dixon, S. L.; Merz, K. M., Jr. J. Am. Chem. Soc. 1998, 120, 5593-5594. (28) van der Vaart, A.; Merz, K. M., Jr. J. Am. Chem. Soc. 1999, 121, 9182-9190. (29) Rizzo, R. C.; Jorgensen, W. L. J. Am. Chem. Soc. 1999, 121, 48274836. (30) Kortemme, T.; Ramirez-Alvarado, M.; Serrano, L. Science 1998, 281, 253-256. (31) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926-935. (32) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187-217. (33) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179-5197. (34) Yang, W.; Lee, T.-S. J. Chem. Phys. 1995, 103, 5674-5678. (35) Dixon, S. L.; Merz, K. M., Jr. J. Chem. Phys. 1996, 104, 66436649. (36) Dixon, S. L.; Merz, K. M., Jr. J. Chem. Phys. 1997, 107, 879893. (37) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. J. Am. Chem. Soc. 1985, 107, 3902-3909. (38) Stewart, J. J. P. J. Comput. Chem. 1989, 10, 209-220. (39) Stewart, J. J. P. J. Comput. Chem. 1989, 10, 221-264. (40) Dixon, S. L.; van der Vaart, A.; Gogonea, V.; Vincent, J. J.; Brothers, E. N.; Sua´rez, D.; Westerhoff, L. M.; Merz, K. M., Jr. DiVCon 99; The Pennsylvania State University, 1999. (41) Li, J.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. A. 1998, 102, 1820-1831. (42) Gogonea, V.; Merz, K. M., Jr. J. Phys. Chem. A 1999, 103, 51715188. (43) Zheng, Y. J.; Merz, K. M. Jr. J. Comput. Chem. 1992, 13, 11511169. (44) Dannenberg, J. J. J. Mol. Struct. 1997, 401, 279-286. (45) Pearl, G. M.; Zerner, M. C. J. Am. Chem. Soc. 1999, 121, 399404. (46) Rick, S. W.; Stuart, S. J.; Berne, B. J. J. Chem. Phys. 1994, 101, 6141-6156. (47) Caldwell, J.; Dang, L. X.; Kollman, P. A. J. Am. Chem. Soc. 1990, 112, 9144-9147. (48) Applequist, J.; Carl, J. R.; Fung, K.-K. J. Am. Chem. Soc. 1972, 94, 2952-2960. (49) Kuwajima, S.; Warshel, A. J. Phys. Chem. 1990, 94, 460-466. (50) Wallqvist, A.; Ahlstrom, P.; Karlstrom, G. J. Phys. Chem. 1990, 94, 1649-1656. (51) Gresh, N. J. Chim. Phys. 1997, 94, 1365-1416. (52) van der Vaart, A.; Merz, K. M., Jr. To be published.
