Effective Charges for Macromolecules in Solvent - The Journal of

3868

J. Phys. Chem. 1996, 100, 3868-3878

Effective Charges for Macromolecules in Solvent R. R. Gabdoulline and R. C. Wade* European Molecular Biology Laboratory, Meyerhofstrasse 1, D-69012 Heidelberg, Germany ReceiVed: October 13, 1995X

In order to study protein-protein and protein-DNA association, the electrostatic forces and interaction free energies for two macromolecules at different mutual orientations and separations need to be evaluated. Realistic systems typically consist of thousands of atomic charges in an environment with a nonuniform dielectric permittivity and a solvent of nonzero ionic strength. Consequently, accurate evaluations of electrostatic forces and energies for such systems are only computationally feasible for a limited number of macromolecule positions. Here, we show that, by representing each molecule by a small number of effective charges in a uniform dielectric, intermolecular electrostatic interactions can be calculated simply and with high accuracy. The effective charges are derived by fitting them to reproduce the molecular electrostatic potential calculated by numerical solution of the finite-difference linearized Poisson-Boltzmann equation. The derived charges are expected to be useful in applications such as the simulation of the diffusional encounter of proteins where they will provide improved accuracy over the commonly-used test charge approximation.

Introduction Electrostatic forces and interaction free energies between biological macromolecules can be calculated by numerical solution of the Poisson-Boltzmann equation. The finite difference technique1-4 may be used in which the PoissonBoltzmann equation is solved by representing three-dimensional space by a finite cubic lattice. An alternative approach is the boundary element method,5-8 in which the problem is reduced to two-dimensional integral equations on the macromolecular surfaces. These methods can give reliable results for proteins consisting of several hundred amino acid residues, but they require considerable computational resources. This means that, in practice, the evaluation of electrostatic interaction free energies by these methods can only be carried out for a limited number of positions of the macromolecules. However, for the Brownian dynamics study of the association of two macromolecules, for example, interaction forces need to be evaluated for many different arrangements of the molecules and such computationally demanding calculations of the forces are not feasible. In applications such as the simulation of the diffusional encounter of proteins, this problem is usually handled by evaluating interaction forces using the electrostatic potential of one molecule (considered rigid) and representing the other by test charges placed in the field of the first molecule (see, for example, refs 9-11). Namely, the second molecule is considered as a set of point charges immersed in the solvent with the charge values equal to the partial atomic charges or, more simply, the net charges of charged groups of atoms. Here, we show that a more rigorous evaluation with much better accuracy may be achieved by representing the second molecule by a set of effective charges in a uniform dielectric medium. By fitting the effective charges to reproduce the external molecular potential calculated by solution of the Poisson-Boltzmann equation, effective charges that model the rescaling of charges that results from the reaction field due to the nonuniform dielectric environment of the molecule are derived. When the potential of a molecule is considered to be * To whom correspondence should be addressed. X Abstract published in AdVance ACS Abstracts, February 1, 1996.

0022-3654/96/20100-3868$12.00/0

due to its partial atomic charges and a set of image charges created by the dielectric boundary, the image charges determining the external potential are all situated in the interior of the molecule. Thus, the external potential can be well represented by effective charges positioned inside the molecule. In contrast, the internal potential of a molecule is dependent on external image charges and cannot be well represented by internal effective charges. The effective charges are thus derived by fitting the macromolecule’s external electrostatic potential. This is done in a manner similar to that by which partial atomic charges12,13 are derived from electron densities calculated by quantum mechanical methods. An advantage of using effective charges compared to a multipole representation of the molecule is that they can be positioned to reflect the shape and charge distribution of the molecule. The accuracy of the representation can be improved by increasing the number of effective charge sites or by assigning higher order multipoles, e.g. dipoles, to each of the effective charge sites. Such a ‘multicenter multipole expansion’ is especially suitable for nonspherical objects such as macromolecules, and it makes it possible to describe the potentials and interactions in the regions where a ‘one-center multipole expansion’ does not converge. As we shall see, the use of a small number of monopole terms results in sufficient accuracy for medium-sized proteins. The applicability of the ‘multicenter multipole expansion’ formalism for determining partial charges from ab initio quantum calculations has already been shown.14 Truncation to monopole terms can be sufficiently accurate, as is discussed by Kollman and Dill.15 A related expansion for macromolecules is given by the inducible multipole solvation model (IMPS),16,17 where multipoles are ascribed to each atom from theoretical considerations. The model is applicable to zero ionic strength, and the electrostatic properties of small molecules can be calculated faster than by finite difference methods. Large systems require comparable computing times because of the need to evaluate pairwise quantities.16 The effective charge model presented here is not intended to describe intramolecular interactions but is instead suitable for describing intermolecular electrostatic interactions simply and accurately. It is straightforward to apply the effective charge model to the case of ionic solvent. © 1996 American Chemical Society

Effective Charges for Macromolecules in Solvent

J. Phys. Chem., Vol. 100, No. 9, 1996 3869

In the following theory sections, we show how effective charges fit to a molecular potential can be used to evaluate those terms of the electrostatic interaction free energy which correspond to the interaction of the charges of one molecule with the potential of the other, with account taken of the reaction field from its own dielectric inhomogeneity. The errors in the evaluation of these terms are determined by the errors in the potential fit. The potential fit procedure used is applicable to the derivation of a large number of effective charges and can be made insensitive to errors caused by the large number of parameters involved. We then describe a test of the accuracy of the effective charge representation with a model molecule having 15 atoms. The results of finite difference calculations of the interaction free energies between two identical model molecules at various orientations and displacements are compared with those from effective charge and test charge representations. In the following section, we apply the effective charge representation to proteins. First, a protein with a simplified charge distribution is studied. Namely, only the partial atomic charges on the “head groups” of the formally charged residues have nonzero values. The term “head group” is used here to denote the carboxylate groups of Asp and Glu residues, the amino group of Lys, and the guanido group of Arg residues. The effective charges derived from fitting the potentials of different conformations of the molecule are compared. The subsequent section deals with full partial atomic charge models of two different proteins having zero and nonzero net charge. Interaction free energies computed by finite differences and by using effective and test charges are compared. Appendix A contains the details of the finite difference computations. Appendix B describes a simplified method for computing the cavity insertion terms of the interaction free energy. These correspond to the interaction of the charges on one molecule with the dielectric inhomogeneity of the other and cannot be represented by effective charges. In Appendix C, we derive the effective charges for a system of charges in a spherical cavity, for which simple analytical estimates are available. Theory In section A, we show which part of the full electrostatic interaction free energy between two molecules effective charges can be used to compute. We show that effective charges can account for more of the total interaction free energy than use of the standard test charge approximation. They do not, however, account for cavity insertion terms, which can be estimated as shown in Appendix B, or for some higher order terms of the interaction free energy whose magnitude we estimate and show to be relatively unimportant. In section B, we describe the method by which effective charges can be derived to reproduce the electrostatic potential surrounding a molecule. We also apply principal component analysis to examine the contributions of individual effective charge sites to the electrostatic potential. A. Electrostatic Interaction Free Energies. The electrostatic potential Φ(r) at a point r in a system of nonuniform dielectric (r) containing molecules with a fixed charge density F(r) surrounded by ionic solution with a Debye-Hu¨ckel parameter k(r) can be described by the linearized PoissonBoltzmann equation

-∇(r)∇Φ(r) + κ(r)‚Φ(r) ) 4πF(r)

(1)

where κ(r) ) (r)k2(r). The electrostatic free energy of the system E is given by the integral

E)

1 3 ∫d r Φ(r) F(r) 2

for which we will use the following notation:

1 E ) (F,Φ) 2 We will consider the charge density F(r) to be given by partial atomic point charges at the nuclei of the molecules. The dielectric permittivity (r) is considered to be constant at s in solvent and to have a different value p inside the molecules. The Debye-Hu¨ckel parameter k(r) is zero inside the molecules and has a value ks in the solvent unless the solvent is nonionic, in which case k(r) is zero everywhere. For a system with two non-overlapping interacting molecules, the following decompositions can be unambiguously defined:

F(r) ) F1(r) + F2(r) (r) ) s + 1(r) + 2(r) κ(r) ) κs + κ1(r) + κ2(r) where Fi(r) is the charge distribution for molecule i, i ) 1, 2; i(r) ) p - s and κi(r) ) -κs if r is inside molecule i; and i(r) ) 0 and κi(r) ) 0 elsewhere. The more complicated case of two overlapping molecules is discussed briefly at the end of this section, along with the definition of overlap. The electrostatic interaction free energy, ∆E12, between two rigid solute molecules (1 and 2) immersed in ionic solvent can be calculated by numerical solution of the finite difference (FD) linearized Poisson-Boltzmann (LPB) equation as

∆E12 ) E12 - E1 - E2 where E12 ()E) is the free energy of the system with molecules 1 and 2 present and Ei is the free energy of the same system when only molecule i is present. Ei ) (Fi,Φi(0)), where Φi(0) is the potential of molecule i in the absence of the other molecule. The interaction free energy, ∆E12, needs to be computed at different intermolecular separations and orientations. The electrostatic potential of the system may be decomposed exactly into two terms Φ1 and Φ2, given as solutions to

-∇(r)∇Φi(r) + κ(r)‚Φi(r) ) 4πFi(r), i ) 1, 2 Thus, Φi is the potential defined by the charges of molecule i in the presence of the dielectric inhomogeneities of both molecules. The full electrostatic interaction free energy ∆E12 can then be decomposed into four components that are dependent on Φ1 and Φ2:

∆E12 ) 1/2((F1,Φ1) - (F1,Φ1(0))) + /2((F2,Φ2) - (F2,Φ2(0))) + 1/2(F2,Φ1) + 1/2(F1,Φ2)

1

) ∆E1-c + ∆Ec-2 + 1/2∆E1-2 + 1/2∆E2-1

(2)

The first two terms, which we label ∆E1-c and ∆Ec-2 with the subscript c represent “cavity”, arise from the effect of the dielectric cavity of one molecule on the potential of the other. The last two terms, 1/2∆E1-2 and 1/2∆E2-1, correspond to the interaction of each molecule’s charges with the potential of the other in the presence of both molecular cavities. They are equal, and therefore, we will use ∆E1-2 ) 1/2∆E1-2 + 1/2∆E2-1 to represent their sum. The effective charges Fieff for molecule i are defined as point charges in uniform solvent, giving the same potential Φi(0) outside molecule i (to within the fitting accuracy) as the partial

3870 J. Phys. Chem., Vol. 100, No. 9, 1996 atomic charges in the inhomogeneous dielectric system consisting of molecule i immersed in solvent. Let us first consider the case when molecule 1 is a set of charges in a dielectric cavity but molecule 2 is simply a set of charges immersed directly in solvent rather than in a dielectric cavity. This case may be realized when the atoms of molecule 2 have vanishing radii. The first term, ∆E1-c, is then 0. The second term, ∆Ec-2, accounts for the effect of the reaction field of the interior of molecule 1 on Φ2. It reflects the decrease in solvation of molecule 2’s charges due to the presence of the low dielectric cavity of molecule 1, and it depends only on the shape and size of molecule 1 (not on its charges). The third term, 1/2∆E1-2, is given exactly by the test charge approximation in this special case, i.e. by multiplying molecule 2’s charges by the value of the potential of molecule 1 at their location. The fourth term, 1/2∆E2-1, is equal to the third. The use of the test charge approximation for molecule 1 to evaluate the fourth term would, however, make it different from the third term, since the test charge approximation treats molecule 1’s charges as if they were immersed in uniform solvent rather than a low dielectric cavity and consequently results in the incorrect assumption that Φ2 ) Φ2(0). On the other hand, if we replace molecule 1 by a set of effective charges in uniform solvent, the symmetry of the third and fourth terms is restored. The use of effective charges instead of molecule 1 correctly describes the interaction of molecule 2’s charges with the potential of molecule 1 (term 3s1/2∆E1-2), as well as the interaction of molecule 1’s charges with the potential of molecule 2 (term 4s1/2∆E2-1). In general, the simplification above is not applicable; i.e., molecule 2 generally has a different dielectric constant from the solvent. This means that term 1 (∆E1-c) describing the dependence of the self interaction of molecule 1’s charges on the dielectric inhomogeneity of molecule 2 will be nonzero. Term 2 (∆Ec-2) will be more complex because of the presence of two dielectric cavities. In addition, the representation of both molecules by effective charges will not provide an exact description of terms 3 and 4 (∆E1-2) anymore. Higher order components of these terms (estimated below) which vanish rapidly with increasing intermolecular separation are not accounted for in the effective charge representation. Thus, for the general case, the electrostatic interaction free energy will be evaluated as the sum of two cavity insertion terms (1 and 2 in eq 2), which will be estimated as described in Appendix B, plus the interaction free energy between either the effective charges on both molecules or the effective charges of one molecule with the potential (Φi(0)) of the other (calculated with only one molecule present). The errors due to the neglected components in this estimate of the electrostatic interaction free energy can be examined by considering Green’s functions. Let us define the Green’s functions Gs12, Gs1, Gs2, and Gs for operators (4π)-1‚(-∇∇ + κ) with ) s + 1 + 2 and κ ) κs + κ1 + κ2 for Gs12; ) s + i and κ ) κs + κi for Gsi with i ) 1, 2; and ) s and κ ) κs for Gs. In these notations

Φ1 ) Gs12F1 and Φ2 ) Gs12F2 where GF ) ∫d3r′ G(r-r′)‚F(r′). The kernel of the Green’s function Gs, Gs(r-r′), may be written analytically as

Gs(r-r′) )

e-ks|r-r′| s|r - r′|

The effective charges for molecule i, Fieff, are derived such that, to within the fitting accuracy,

Gabdoulline and Wade

GsiFi = GsFieff

(3)

outside the molecule i, i ) 1, 2. The estimation of the terms in eq 2 makes use of the approximations

Gs12F1 ≈ Gs2F1eff and Gs12F2 ≈ Gs1F2eff

(4)

In the special case discussed above when molecule 2 has no dielectric cavity, the relations in eq 4 are equalities not approximations. In general, however, they are approximations and here we estimate the dependence of the errors due to these approximations on intermolecular separation distance. Gs2F1eff and Gs1F2eff can be estimated as described in Appendix B. Using relations 3 and 4, terms 1 and 2 in eq 2 are given by

(F1,Gs12F1) - (F1,Gs1F1) ≈ (F1eff,Gs2F1eff) - (F1eff,GsF1eff) (F2,Gs12F2) - (F2,Gs2F2) ≈ (F2eff,Gs1F2eff) - (F2eff,GsF2eff) (5) Term 3 in eq 2 can be expressed as follows by using relations 3 and 4 above and the property of Gs2 and Gs to be self-adjoint:

(F2,Gs12F1) ≈ (F2,Gs2F1eff) ) (Gs2F2,F1eff) = (GsF2eff,F1eff) ) (F2eff,GsF1eff) = (F2eff,Gs1F1) Term 4 in eq 2 is given by an equivalent expression. The importance of the approximations invoked above can be evaluated by using a series representation of the Green’s functions ∞

Gsi ) [∑(GsLi)n]Gs n)0

where Li ) (4π)-1‚(∇i∇ - κi). This expansion is obtained when the equation

[-∇(s + i)∇ + (κs + κi)]φ ) 4πF is solved iteratively as

(4π)-1[-∇s∇ + κs]φm ) F + Liφm-1, m ) 0, 1, ... and φ-1 ) 0 Then the potential φm is given by the sum of terms with n ) 0, ..., m in the expansion ∞

GsiF ) [∑(GsLi)n]GsF n)0

This representation expands the potential GsiF (created by charges F in the presence of the dielectric cavity of molecule i) by iteratively perturbing the potential due to F in uniform solvent (GsF) to that in the presence of the dielectric inhomogeneity. An expansion for Gs12 is, similarly, ∞

Gs12 ) [∑(Gs(L1 + L2))n]Gs

(6)

n)0

Expressing approximation 4 in series form: ∞

∞

n)0

n)0

Gs12 ≈ (∑(GsL2)n)(∑(GsL1)n)Gs

(7)


J. Phys. Chem., Vol. 100, No. 9, 1996 3871

We see that the first energy term neglected when term 3 of eq 2 is computed is (F2,GsL1GsL2GsF1) although higher order terms of eq 6 are evaluated in expression 7. The first neglected term reads as potential GsF1 results from charges F1 in uniform solvent; this potential induces image charges at molecule 2’s dielectric inhomogeneity L2GsF1. The potential due to these induced charges is GsL2GsF1; this potential induces image charges at molecule 1 (L1GsL2GsF1). The interaction of these induced charges with F2 gives the neglected term. Let us consider the case of nonionic solvent, large intermolecular separations, and molecules having nonzero net charge. Charges induced on molecule 1 by the potential from molecule 2 may be represented by a dipole, the moment of which is proportional to the electric field at the center of molecule 1. Then the neglected contributions in terms 3 and 4 of eq 2 may be estimated to have a leading order of (r12)-7, where r12 is the intermolecular distance parameter, defined, for example, as r12 ) |Rc1 - Rc2|/(R1 + R2), where Rci is the coordinate of the ith molecule center and Ri is the largest center-to-atom distance for molecule i. The leading term of ∆Ec-2 and ∆E1-c is of order (r12)-4, and the first neglected term has order (r12)-10. In the simple test charge approximation, ∞

Gs12 ≈ Gs1 ) (∑(GsL1)n)Gs n)0

and, consequently, the first neglected term of ∆E1-2 of eq 2 is (F2,GsL2GsF1) with an order of (r12)-2. We now consider molecular overlap. Difficulties in calculating the cavity insertion terms (∆E1-c and ∆Ec-2) occur when the dielectric cavity cannot be ascribed to a single molecule. This obviously occurs when two molecules are in van der Waals contact. This situation can also occur when the molecules do not touch and the solvent accessible molecular surface is used to define the dielectric boundary. When there is no space for the solvent probe between the two molecules, a dielectric region appears which belongs to the interface rather than to any one of the molecules. It is appropriate to define this as a limit for the use of the formalism described here. It requires the molecular surfaces to be separated by more than twice the probe radius used in calculating the solvent accessible molecular surface to define the dielectric boundary. B. Potential Fitting. We define the effective charge distribution Feff for a molecule as the point charge distribution n qjeff‚δ(r - rj) in uniform solvent giving the same Feff ) ∑j)1 potential Φ(0) outside the molecule (to within the fitting accuracy) as the partial atomic charges in the inhomogeneous dielectric system consisting of the molecule immersed in solvent. The molecular electrostatic potential of the molecule Φ(0)(r) may be calculated by numerically solving the finite-difference Poisson-Boltzmann equation taking account of the partial charges on all the atoms, the inhomogeneous dielectric, and a surrounding ionic solution. Effective point charges may be derived to describe the electrostatic interactions of a macromolecule by fitting them to reproduce its molecular electrostatic potential. This can be done in the same way as partial atomic charges are derived from electron densities.18 n effective charges qjeff at positions rj are obtained by minimizing a fitting functional:

J ) ∫d3r |Φ(0)(r) - ∑qjeff‚Fj(r)|2 Fj(r) )

e-ks|r-rj| s‚|r - rj|

(8)

Integration is carried out over the region of space relevant to intermolecular interactions. The solution is given by a linear system of equations

Aqeff ) b

(9)

where A is a matrix with elements aij,

aij ) ∫d3r Fi(r)‚Fj(r) qeff ) (q1eff, q2eff, ..., qneff) and

b ) (b1, b2, ..., bn), bi ) ∫d3r Fi(r)‚Φ(0)(r) The fitting accuracy can be measured by

χ ) 1 - ∫d3r |Φ(0)(r) - ∑qjeff‚Fj(r)|2/∫d3r |Φ(0)(r)|2

(10)

It has a maximum value of 1 (100%) when the potential is fitted by itself and is 0 when the fitting potential is zero everywhere. The contributions of the individual effective charges to the potential are not orthogonal, but orthogonal (or generalized) contributions can be constructed from linear combinations of effective charges. These linear combinations may be obtained by diagonalization of the matrix A. We use the orthogonal matrix X (X-1 ) XT)19 with elements xij, so that X-1AX ) diag(λ1, λ2, ..., λn). Then the generalized charges q˜ i (which are linear combinations of the charges on specific sites) with q˜ i ) ∑j xji‚qjeff form an orthogonal set of vectors, and the fitting problem is just

λi‚q˜ i ) b˜ i

(11)

where b˜ i ) ∑k xki‚bk. The functional J has a simple form:

J ) J0 + ∑λi‚(q˜ i - b˜ i/λi)2

(12)

i

with J0 ) ∫d3r |Φ(0)(r)|2 - ∑ib˜ 2i /λi as a minimum value. This representation is useful to distinguish the generalized charges that are more important for the representation of the electrostatic potential from the relatively unimportant generalized charges. The smaller the eigenvalue, the less important the value of the corresponding generalized charge is. However, when effective charges are constructed by inverse transformation from the generalized charges, it is apparent that all of the generalized charges q˜ i can make comparable contributions to the values of the effective charges qjeff. Thus, very small changes in the electrostatic potential to be fit can result in sizable changes in some effective charges. This situation is analogous to the case of partial atomic charge determination12,13 for small molecules when the charges on some atoms appear to make only small contributions to the molecular potential. If the fitting is carried out in a spherical shell region at a large distance from the center of the molecule, effective charges are obtained that are equivalent to a one-center multipole expansion of the atomic charges and their image polarization charges resulting from the dielectric discontinuity. Model Molecule We have studied the accuracy of the interaction energies calculated with test and effective charges using the small model molecule shown in Figure 1. The molecule has four unitpositive (marked by +) and three unit-negative (-) atoms; the remaining eight atoms are neutral. The atoms all have radii of 1.2 Å and are located at a distance of 2 Å from the central atom. Their coordinates are given in Table 1. The molecular

3872 J. Phys. Chem., Vol. 100, No. 9, 1996


Figure 1. Model molecule used for comparison of different energy calculation methods. There are seven charged atoms (four having charge +1 e, three having charge -1 e). The dot surface shows the dielectric boundary between the solvent and the interior of the molecule.

TABLE 1: Effective Charges Derived for the Model Molecule Shown in Figure 1a atom

Qatom

Qeff

x

y

z

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 net molecular charge

+1 +1 0 0 0 0 -1 +1 0 0 -1 0 -1 +1 0 +1

1.29557 1.53335 0.0 0.0 0.0 0.0 -1.58333 1.51519 0.0 0.0 -1.56678 0.0 -1.63584 1.44167 0.0 0.99983

0.000 1.155 -1.155 -1.155 1.155 2.000 0.000 -2.000 0.000 0.000 1.155 -1.155 -1.155 1.155 0.000

0.000 1.155 1.155 -1.155 -1.155 0.000 2.000 0.000 -2.000 0.000 1.155 1.155 -1.155 -1.155 0.000

2.000 1.155 1.155 1.155 1.155 0.000 0.000 0.000 0.000 0.000 -1.155 -1.155 -1.155 -1.155 -2.000

Figure 2. Electrostatic interaction free energy between two identical model molecules at different center-to-center separations along the x coordinate. The shifted molecule was rotated by 90° about the y axis before translating it along the x axis. Energies ∆E1-2 are calculated with the following: FD method (diamonds); effective charge model (full lines); test charges (dashes). Full interaction free energies ∆E12 are also shown (squares).

a Actual charges (Q atom), effective charges (Qeff), and Cartesian coordinates (x, y, z) are given for each atom. Fitting accuracy for the electrostatic potential is χ ) 0.9999.

potential was calculated using the UHBD program as described in Appendix A. The potential of the molecule in the spherical region R ) 7-13 Å, where R is the distance from the central atom, was fitted by seven effective charges at the positions of the charged atoms. The resulting charges are listed in Table 1. The potential is fitted with an accuracy of χ ) 0.9999. The use of a nine term expansion of multipoles on the center of the molecule (using a monopole and three dipole and five quadrupole terms) gives χ ) 0.9996. The difference is negligible when the potential is fitted at larger distances (here, R > 13 Å). Ten times worse accuracy (χ ) 0.997 with nine multipoles vs χ ) 0.9998 with seven effective charges) is seen for fitting at shorter distances R ) 7-8 Å. The effective charges are considerably larger in magnitude than 1.0 (1.3-1.6), but the sum of the charges is +1 to high accuracy. This is a consequence of fitting the potential in a region enclosing the charges: the monopole terms defined by the net charge are dominant in the integrand of the fitting functional and are, consequently, fitted first. Figures 2-5 show the electrostatic interaction free energies of two identical model molecules at different separations and orientations obtained by five different methods. One of the molecules is moved away from the other (Figs. 2 and 3) or rotated at a center-to-center distance of 10 Å (Figs. 4 and 5). In the calculations, presented in Figure 2, the shifted molecule was first rotated by 90°. Squares show the full electrostatic interaction free energy ∆E12 of the system of two charged molecules calculated by numerical solution of the FD LPB equation. Diamonds show the energy ∆E1-2 ) ∆E12 - ∆E1-c - ∆Ec-2 (see eq 2). The free energy ∆E1-2 differs from the full interaction free energy by terms corresponding to the introduction of the dielectric

Figure 3. Electrostatic interaction free energy as in Figure 2, but for model molecules separated along the z axis.

cavity of each molecule into the electric field of the other. These terms do not appear in either the effective or test charge approximations of the free energy. The interaction of a low dielectric cavity with a charge can easily be estimated for this model system in the same way as in ref 16 by keeping only the dipole contributions of the induced polarization and replacing the low dielectric cavity by a single low dielectric sphere. These terms are positive and, for the model system, are on the order of 0.05 kcal/mol, i.e. less than the difference between using effective and test charges but more than the difference between using effective charges and the finite difference solution. A more detailed description of a simplified method for accounting for these terms is given in Appendix B. Test charge free energies (dash) were calculated by multiplying the atomic charges (+1 or -1) of the second molecule by the value of the first molecule’s potential at the position of the charges. The effective charge free energies (full line) were calculated similarly using potential fitted charges on the second molecule. The agreement between the free energies calculated with the accurate (finite difference) method and the effective charges is good at all separation ranges (see Figures 2 and 3). The discontinuity at the center-to-center separation of 28 Å occurs because at this separation one of the charges is outside the region


Figure 4. Electrostatic interaction free energy between two identical model molecules when one of the model molecules is shifted with respect to the other by 10 Å along the x axis and then rotated gradually about the y axis. Note that angle ) 90° corresponds to the same arrangement as x ) 10 Å in Figure 2. Energy calculation methods are as in Figure 2.

Figure 5. Electrostatic interaction free energy between two identical model molecules when one of the model molecules is shifted with respect to the other by 10 Å along the y axis and then rotated gradually about the z axis. Energy calculation methods are as in Figure 2.

where the potential was calculated. The discontinuity at 7 Å is caused by overlap of the two molecules. The results for different orientations (Figures 4 and 5) show that, at the separation investigated here, the effective charges give an interaction free energy with an error comparable to that of the finite difference method. The test charge model, in general, underestimates the amplitude of angular variations of interaction energies. Of course, the picture seen here is dependent on the model chosen for computation. However, as we will see later, the underestimation observed here should always be seen for cases where the interaction is not solely dependent on monopole terms. The same relative features were seen when the calculations were carried out with ionic solvent. Applications A. Potential Fit. The model molecule results show that electrostatic interaction free energies can be calculated accurately with effective charges placed at the positions of all the charged atoms. The small size and rather spherical shape of the test molecule are not important for the accuracy of the evaluated interaction free energy ∆E12. This is demonstrated by application of the method to two larger and more realistic protein

J. Phys. Chem., Vol. 100, No. 9, 1996 3873

Figure 6. Distance of each atom of IL-4 from the center of geometry of the molecule. Diamonds show the distances for the 37 effective charge sites positioned at the formally charged residues.

systems. We have carried out the fitting of the potential for two proteins: interleukin-4 (IL-4), having 129 residues and a net charge of +7 e (assuming standard protonation states for the residues at pH 7), and interleukin-2 (IL-2), having 134 residues and a net charge of 0 e; both of the molecules are far from spherical in shape. Figure 6 shows the distance of all 1316 atoms of IL-4 from the center of geometry of the molecule; the diamonds indicate the distance from the center for the 37 effective charge sites used below. Protein with a Simplified Charge Distribution. The first computations were carried out on a simplified model of IL-4 in which only the head group atoms of charged residues were assigned nonzero charges (240 charged atoms). The potential was calculated by numerical solution of the FD LPB equation. The effective charges were positioned at the geometric centers of the 35 head groups, and their values were derived by fitting the potential over distances 30-45 Å from the center of the molecule. A high accuracy of χ ) 0.998 was obtained. The diagonlization of the fitting matrix A in eq 9 shows that the errors are only doubled when only the half of the generalized charges with the largest eigenvalues are used. However, the generalized charges corresponding to the lower eigenvalues make approximately the same contribution to the effective charge values, so the best fit effective charges vary considerably according to how many generalized charges are used. This problem is also encountered in the derivation of partial atomic charges from ab initio quantum mechanical electron densities and may be handled similarly12,13,18 by introducing constraints that force the effective charges to have reasonable values. Some differences are, however, evident. The net charge on the molecule and the lower order multipole moments should not be artificially constrained because they appear to be the main contributors to the fitted potential and are determined by the higher eigenvalue generalized charges. Analysis of the lower order moments of the effective charges shows that the scaling of the atomic charge distribution (due to polarization screening on the molecular surface) is different for monopole and higher order terms. The generalized charge corresponding to the highest eigenvalue is the same when derived from the initial atomic charge distribution on the charge sites and when derived from potential fitting. This simply means that the net polarization charge is exactly the same as for a point charge in a solvent, leading to scaling of the charge by 1. However, the dipole and higher moments are screened less, so that the dipole moment of a set of charges inside a single low dielectric region (and the polarization charges caused by them) is about 1.5 times larger than if these charges were individually immersed into the

3874 J. Phys. Chem., Vol. 100, No. 9, 1996

Figure 7. Correlation plot of the generalized charges derived from the initial partial atomic charges vs the potential derived generalized charges. The charges are calculated for 17 different conformations of IL-4 with a simplified model in which partial atomic charges are only assigned to the head group atoms of formally charged residues. The numbering of the generalized charges is such that the corresponding eigenvalues are in descending order. Generalized charge 1 (represented as a square in the correlation plot) corresponds to monopole, generalized charges 2-4 (diamonds), to dipole, and generalized charges 5-9 (crosses), to quadrupole terms in the usual expansion of a charge distribution.

Figure 8. Set of effective charges obtained by minimizing the MEPfitting functional for the simplified model of IL-4. The dotted line indicates the (sum of) partial atomic charges at the corresponding effective charge sites.

solvent. This can be shown analytically when the low dielectric region is a sphere (see Appendix C). Figure 7 shows that, for the model IL-4 system, the average scaling factor for dipole and quadrupole terms is, on the basis of the scaling factors for the generalized charges for 17 different conformations of IL-4, about 1.6. Figure 8 shows the set of effective charges calculated by minimizing the fitting functional (eq 8). They show considerable conformational dependence and can be very large in magnitude. The presence of combinations of effective charges that are unimportant for potential fitting permits the derivation of effective charges that are less sensitive to molecular conformation. One way to achieve this is to force the site charges to have values close to the net charges of the atomic groups replaced by the site (here, (1). The regularization procedure we used to do this involves simply assigning corresponding linear combinations of atomic charges to the unimportant (low eigenvalue) generalized charges; this results in more conformationally invariant effective charges. Another approach is to compute multiple-conformation molecular


Figure 9. Set of effective charges obtained after the regularization procedure described in the text for 17 conformations of the simplified model of IL-4. The six generalized charges with the lowest eigenvalues were set to the values derived from partial atomic charges. The multiple conformation MEP-derived charges are shown by the dotted line.

electrostatic potential (MEP)-derived charges following ref 20. We did this by applying the same weighting terms to each of 17 different conformations of IL-4 obtained from different structure determinations and molecular dynamics simulations (see Appendix A). Figure 9 shows the set of effective charges after application of the regularization procedure and the effective charges derived from the multiple-conformation MEP. The regularization involved simply equating the six generalized charges corresponding to the lowest eigenvalues to the values calculated from partial atomic charges. The correlation between the multiple-conformation MEP-derived charges (broken lines) and the set of regularized charges is reasonably good. There are several sources of uncertainty in the generalized charges corresponding to the low eigenvalues (and, consequently, in the values of the effective charges). Here we discuss two of them. (1) The errors due to the operations of matrix inversion are small: only two or three of the generalized charges with the smallest eigenvalues were affected (at most by a factor of 1.1) when the computations were carried out to single precision (throughout the study, we used double precision). Bigger errors were found for larger matrices when more effective charge sites were used. For example, on inversion of a matrix derived for 124 effective charge sites (used for IL-2 in the next section) the 15 smallest eigenvalues differ by more than 10% for single and double precision calculations. These errors make a negligible contribution to the potential fit, but their contribution to the effective charge values may be large. (2) The effective charges are dependent on the positions of the charge sites. The clearest cause of large effective charge variations is the close proximity of two charge sites. Assigning effective charges of opposite sign to close charge sites produces a small dipole. The contribution to the potential of such effective charges is small, and consequently, large effective charges can be assigned with little effect on fitting accuracy. Quantitatively, this singularity is reflected in the lowest eigenvalues of the matrix A. Roughly speaking, a numerical procedure with precision δ (defined here as a maximal δ satisfying the ‘computer’ equality 1 + δ ) 1) gives correct values for those eigenvalues which are bigger than δ times the largest eigenvalue. The number of eigenvalues within this limit is determined by the positions of the effective charges. The regularized effective charges described above as well as the multiple-conformation MEP-derived effective charges are, not surprisingly, less sensitive to these errors. Figure 10 shows


Figure 10. Logarithm of the eigenvalues of the fitting matrix A obtained for the simplified model of IL-4. The bold line corresponds to eigenvalues of the matrix obtained by averaging the A matrices corresponding to specific conformations of IL-4; the other lines show the eigenvalues of matrices derived from eight individual conformations of IL-4.

J. Phys. Chem., Vol. 100, No. 9, 1996 3875

Figure 12. Generalized charges for different conformations of the simplified model of IL-4. The bold line corresponds to multipleconformation charges derived from 17 conformations.

Figure 13. Two interacting IL-2 molecules. The electrostatic free energies were calculated for different orientations of the right-hand molecule, which was rotated about the axis perpendicular to the page and passing through the molecular center. The center-to-center separation was 48 Å.

Figure 11. Right-hand sides (b˜ i) of eq 11 for eight conformations of the simplified model of IL-4. The bold line corresponds to the RHS of the multiple-conformation equation (derived from 17 conformations).

the eigenvalues of the fitting matrix for eight individual conformations and for the average matrix for 17 conformations of IL-4. The spectrum of the average matrix is more uniform (flat). At the same time, the averaging of the right-hand sides bi in eq 9 leads to more uniformly distributed b˜ i in eq 11 (see Figure 11). Consequently, there is no noticeable change in the distribution of the generalized charges (solutions of eq 11), as is seen from their plot in Figure 12. One important change is, however, evident. The last few generalized charges, with the smallest eigenvalues, are not as large in the average matrix as they are for the individual conformations. Protein with Full Partial Atomic Charge Distribution. A second set of potential calculations was carried out with charges on all atoms of the protein. The assignment of effective charges to the positions of the 35 charged residues and the two termini to fit the potential from all the protein’s partial atomic charges results in reasonable accuracy with significant simplification of the system. In the worst case, when the fit is made over the region 2 atomic radii from any atom of molecule, the fit accuracy is χ ) 0.97. As we will see in the next section, this fit accuracy allows interaction free energies to be calculated with good accuracy. B. Electrostatic Interaction Free Energies. Below, we present three examples of the use of effective charges to

calculate electrostatic interaction free energies between two identical protein molecules. Interleukin 2. The system consisted of two IL-2 molecules in water. One of the molecules was kept fixed while the other was rotated. The rotation axis was perpendicular to the plane of Figure 13, which shows the mutual positions of the two IL-2 molecules. During rotation, the minimum distance between the atoms of the two different molecules varied from 5-12 Å. Figure 14 shows the interaction free energy ∆E1-2 ) ∆E12 ∆E1-c - ∆Ec-2 for two IL-2 molecules with 48 Å separation between their centers. The finite difference free energies have relatively small errors, as is seen from a comparison of results from calculations with different grid spacings. The finite difference free energies were computed with grid spacings 1.0, 0.8, and 0.65 Å; the effective charges were obtained by fitting the potential with grid spacing 0.8 Å. Using 32 test charges positioned at the centers of the head groups of the charged residues gives a free energy estimate of approximately the same quality as using test charges on all atoms and, as noted for the model molecule, results in underestimation of changes in interaction free energy with changes in orientation. However, 32 effective charges on the charged residues give, on average, a much better fit of the free energy. Further imporvement in fitting was obtained by assigning 26 more effective charge sites on the side chains of polar residues. Essentially no further improvement was gained by the use of 124 effective charge sites on the side chains of the residues. The cavity insertion terms ∆E1-c + ∆Ec-2 can be calculated as shown in Appendix B to obtain the full interaction free energy ∆E12. Interleukin 4. The same features may be seen when the interaction free energies are calculated for two IL-4 molecules

3876 J. Phys. Chem., Vol. 100, No. 9, 1996

Figure 14. Electrostatic interaction free energy vs rotation angle for two IL-2 molecules. The crosses, squares, and diamonds show free energies (∆E1-2 ) ∆E12 - ∆E1-c - ∆Ec-2) calculated by the finite difference method at grid spacings of 1.0, 0.8, and 0.65 Å, respectively. Full, dot, and dash lines show ∆E1-2 obtained by using, respectively, 32, 58, and 124 effective charges. Dash-dot and bold dash-dot lines show the energies obtained using 32 and 1022 test charges. The full interaction energy ∆E12 can be obtained by adding the cavity insertion terms shown in Figure 17 to ∆E1-2.

Figure 15. Electrostatic interaction free energy vs rotation angle for two IL-4 molecules. Diamonds show free energies ∆E1-2 calculated by the finite difference method at grid spacing of 0.8 Å. Dashed line shows ∆E1-2 computed with 112 effective charges. Dot and dashdot lines show energies obtained using 37 and 1049 test charges, respectively. Full line represents the full interaction free energy ∆E12 obtained by adding the cavity insertion terms as estimated in Appendix B to the ∆E1-2 computed with the effective charges. The full interaction free energy (∆E12) calculated by the finite difference method is shown by squares.

immersed in water at a center-to-center separation of 38 Å (Figure 15). It should be noted that the dielectric surface was defined by the solvent accessible molecular surface in this case and that at angles 210-270° the dielectric cavities of the two molecules were joined; i.e., there was no solvent separation of the two molecules. This does not, however, affect the accuracy of the estimated free energies ∆E1-2 (eq 2), although the estimate of the full interaction energy ∆E12 has larger errors. The same computations carried out for two IL-4 molecules in ionic solvent showed the same relative features. Figure 16 shows the free energies ∆E1-2 for the same orientations of two IL-4 molecules as above but calculated with ionic solvent with a Debye-Hu¨ckel parameter ks ) 0.125 Å-1. Discussion We have studied the electrostatic interactions between two proteins at intermediate separation distances. At larger separa-


Figure 16. Electrostatic interaction free energy vs rotation angle for two IL-4 molecules but for the case of ionic solvent (with DebyeHu¨ckel parameter ks ) 0.125 Å-1). Key as in Figure 15.

tion distances, the electrostatic interaction between charged molecules is mainly governed by monopole terms corresponding simply to the interaction of two point charges in a medium of dielectric permittivity s which can be well described by simple test charge models. However, intermolecular charge-dipole interactions are wrongly described by the simple test charge model (see Appendix C), although at large distances these interactions themselves may be small and hence unimportant. At small separation distances (less than about 3 Å separation between the molecular surfaces of two proteins), a detailed atomic description of solvent is necessary because the interactions are affected by the specific solvent atoms at the molecular interface. The macroscopic description of solvent becomes sensitive to the probe radius used to define the dielectric boundary. Within this limit of macroscopic theory, potential fit effective charges may be derived to describe the electrostatic interaction free energies with high accuracy. At very large distances, the effective charge representation gives the same accuracy as the representation of molecules by one-center multipoles with the same number of parameters. At closer distances, the effective charges may be chosen to reflect the molecular shape. Thus they can describe the potential near the surface of highly nonspherical molecules, where one-center multipoles fail. To achieve accuracy comparable to that of finite difference calculations, a relatively large number of effective charges and, therefore, fitting parameters may be necessary. The eigenvalue analysis of the fitting procedure implemented here resolves problems arising due to numerical uncertainties in effective charge values. The effective charge representation is an extension of the testcharge representation of the protein, when all the charges are scaled by the same factor of 1.0. The analysis here shows that this scaling is exact for the monopole term of the atomic charge distribution of the molecule. The scaling factor for the higher order terms is different, having an average value between 1.0 and 2.0. This is to be compared to two limiting cases of dielectric boundaries: spherical and planar. The potential from a charge q at the center of a dielectric sphere of permittivity p is q/(sr) outside the sphere, where s is the permittivity of the surrounding medium; thus, the scaling factor is 1.0. The potential from the charge q in a dielectric half-space of permittivity p is 2q/(s + p)r (see ref 21), when it is measured in the other half-space of permittivity s; thus, the scaling factor is 2s/(s + p) ≈ 2 when s is much larger than p. In the case of a complex shaped molecule, the dielectric boundary is intermediate between these two cases. However, the first scaling

Effective Charges for Macromolecules in Solvent is always exact for the monopole term of the charge distribution if the molecular surface is of closed form: the net charge can be shown by Gauss’s law to be scaled in exactly the same way as a charge at the center of the sphere. The dipoles and higher order multipoles are less screened by the dielectric boundary than the monopoles. For the separations and molecules investigated here, the screening factor for dipoles and quadrupoles is ca 1.6, and roughly this charge scaling may be expected for molecules of similar size having zero net charge. However, the factor depends on the shape of the molecule. (The same computations on a water molecule, using radii of 1.2 and 1.48 Å for hydrogens and oxygen, respectively, gave a factor of 1.40-1.45.) When the net charge of the molecule is not 0, the effective charges are governed by two different scaling factors for monopole and multipole charges. Resulting effective charges are, on average, larger than 1.0, as is seen for the effective charges of the model molecule (Table 1) and IL-4 (Figure 9). Once computed, effective charges may be used to calculate electrostatic forces in Brownian dynamics simulations and electrostatic interaction free energies in docking studies with the same computational effort as test charges. The time for the computation of effective charges is mainly defined by the matrix inversion operation. For 124 effective charge sites on IL-4, this operation took approximately 1/6th the time necessary to compute the potential of the molecule on a 150 × 150 × 150 grid. Acknowledgment. We thank V. Helms, A. Juffer, and V. Lounnas for critical reading of the manuscript and T. Mu¨ller and H. Oschkinat for providing IL-4 coordinates. Appendix A. Calculation Details The UHBD program4,22 was used to calculate electrostatic interaction free energies and finite difference molecular potentials. A 150 × 150 × 150 grid was used with a spacing of 0.4 Å (model molecule); 0.8 Å (IL-4); and 0.65, 0.8, and 1.0 Å (IL-2). In the last case, the effective charges were fit to the potential with a 0.8 Å spacing grid, but all three spacings were used for finite difference calculations of energy. The OPLS parameter set23 was used to assign partial charges and radii to atoms. An interior dielectric constant of 2 and an exterior dielectric constant of 78 were used. Ionic strengths of 0 and 145 mM were used. The latter value corresponds to a DebyeHu¨ckel parameter ks in eq 1 of 0.125 Å-1. The assignment of the solvent dielectric boundary was based on the van der Waals molecular surface (model molecule, IL-2) and the solvent accessible molecular surface of a probe with radius 1.4 Å (IL4). Dielectric boundary smoothing24 was used in all finite difference calculations. Molecular dynamics simulations of IL-4 were carried out with the GROMOS package25 using, as a starting point, X-ray coordinates from ref 26. NMR coordinates were supplied by T. Mu¨ller et al.27 These two sets of coordinates represent different conformations of IL-4, having approximately a 2 Å root mean square deviations for CR atoms. The 17 conformations used for effective charge derivations include the X-ray and NMR coordinates and 15 coordinates taken from the molecular dynamics simulations at 10 ps intervals. The conformation of IL-2 was a snapshot obtained after 30 ps of molecular dynamics simulation starting from X-ray coordinates from ref 28.

J. Phys. Chem., Vol. 100, No. 9, 1996 3877

Figure 17. Cavity insertion free energy terms for the interaction of two IL-2 molecules. The symbols show values from finite difference calculations, while the lines show values evaluated with formula B1. The upper line is the sum of the lower two, which give the cavity term for each molecule (namely, E1-c and Ec-2).

of a low dielectric cavity associated with another molecule, a simplified form of the IMPS method16 is used. The cavity is modeled by a number of spherical cavities. The change in free energy due to each spherical cavity is evaluated, and these changes are summed. This implies that the sphere-sphere interactions are neglected. We use spheres of the same radius. The best results were obtained using spheres located at the geometric centers of the side chain of each residue of the protein associated with the low dielectric cavity. The radius, ak, of the spheres was adjusted to fit the energies at one orientation of IL-2, and ak ) 3.25 Å. This radius appears to be suitable for any other orientation and, moreover, also for IL-4. It agrees well with the hydrodynamic radii of amino acid residues (see, for the example, ref 30). To calculate the free energy changes, only dipole terms of the induced polarization on each sphere are taken into account. The interaction free energy is then given by the sum

∆E1-c )

s - p

qiqj‚(rik,rjk) ‚a3k s(2s + p) ijk |rik|3‚|rjk|3 ‚∑

(B1)

The triple sum runs over all the charges i and j ) 1, 2, ..., N (twice) and the inducible dipole sites k ) 1, 2, ..., K. It can naturally be decomposed into consecutive double sums by computing all 3K induced dipoles due to N charges and then computing the reaction-field energy on each charge due to each dipole. The computing time is proportional to NK. Consequently, the evaluation of the forces and energies needs a similar amount of computational time to the evaluation of chargecharge pairwise interactions. The formula is valid for nonionic solvent. Figure 17 shows the terms ∆E1-c and ∆Ec-2, estimated from formula B1, for the IL-2 interaction energy calculations. In these computations, the center-to-center separation between the two IL-2 molecules was 48 Å, and the closest atom-to-atom separation of the molecules varied from 5 to 13 Å. The sum ∆E1-c + ∆Ec-2 varies much less than ∆E1-2 and makes a small contribution to the full free energy in this case. The accuracy of the approximations above and possible corrections are described in ref 16.

Appendix B. Cavity Insertion Terms

Appendix C. Effective Charge Representation of the Electrostatic Potential Due to a Point Charge in a Spherical Cavity

In order to compute the change in free energy of the charges associated with one molecule in a solvent due to the presence

Let us consider a sphere of radius R located at the origin, having a dielectric permittivity p, and surrounded by a medium

3878 J. Phys. Chem., Vol. 100, No. 9, 1996


of permittivity s. A charge q, located at point d ) (0, 0, d) inside the sphere produces a potential φ(r,θ) outside the sphere given by the following expansion (the derivations are similar to those in ref 29):

q

∞

2n + 1

dn

φ(r,θ) ) ‚∑ ‚ ‚Pn(cos θ) (C1) s n)0(n + 1) + np/s rn+1 This expansion implicity includes the radius of the sphere and is only valid when d < R and r > R. Let us compare this expansion of the potential with the expansion of the potential of a single effective charge qeff located at the same point d in a uniform dielectric s:

φfit(r,θ) )

qeff s‚|r - d|

)

qeff ∞ dn ‚∑ ‚Pn(cos θ) s n)0rn+1

(C2)

The presence of the variable coefficient C(n) ) 2n + 1/[(n + 1) + np/s], changing monotonically with n, makes a unique assignment of qeff impossible. However, there are several limiting cases, when a representation with a single charge in a uniform dielectric is accurate. (1) The charge is at the center of the sphere: only the n ) 0 term of expansion (C1) is nonzero, C(0) ) 1, and qeff ) q. (2) The potential at large distances (r . d, r > R) is considered. Then the n ) 0 term is dominant with C(0) ) 1 and qeff ) q. (3) There are two opposite charges at d+ ) (0, 0, d) and d- ) (0, 0, -d). Then we have two expansions like C1 for each of the charges, and they should be approximated by two expansions like C2. The monopole, n ) 0, term of one charge’s expansion cancels the monopole term of the other’s. Expansion C1 effectively starts from n ) 1 with C(1) ) 3/(2 + p/s) ≈ 1.5 if p , s. The effective charges to fit the potential at large distances are qeff ≈ ( 3/(2 + p/s)‚q ≈ (1.5q. This gives a scaling factor of 1.5 for a dipole within the low permittivity dielectric. Quadrupoles and higher multipoles will be screened even less. (4) Charge q and r are close to the boundary of the sphere (but still on opposite sides: d < R < r). Higher order terms then determine the expansion, and the coefficient C should be taken at its limiting value C(∞) ) 2/(1 + p/s) ≈ 2.0. The corresponding effective charge is qeff ≈ 2/(1 + p/s)‚q ≈ 2.0q, which is the same as for a planar dielectric boundary (see ref 21).

References and Notes (1) Warwicker, J.; Watson, H. C. J. Mol. Biol. 1982, 157, 671. (2) Klapper, I.; Hagstrom, R.; Fine, R.; Sharp, K.; van Keulen, A. M.; van der Ploeg, A.; Berendsen, H. J. C. Proteins: Struct. Funct. Genet. 1986, 1, 47. (3) Davis, M. E.; McCammon, J. A. J. Comput. Chem. 1990, 11, 401. (4) Madura, J. D.; Briggs, J. M.; Wade, R. C.; Davis, M. E.; Luty, B. A.; Ilin, A.; Antosiewicz, J.; Gilson, M. K.; Bagheri, B.; Scott, L. R.; McCammon, J. A. Comput. Phys. Commun. 91, 57-95. (5) Zauhar, R. J.; Morgan, R. S. J. Comput. Chem. 1990, 11, 603. (6) Juffer, A. H.; Botta, E. F. F.; van Keulen, B. A. M.; van der Ploeg, A. J. Comput. Phys. 1991, 97, 144. (7) Yoon, B. J.; Lenhoff, A. M. J. Phys. Chem. 1992, 96, 3130. (8) Zhou, H.-X. Biophys. J. 1993, 65, 955. (9) Northrup, S. H.; Boles, J. O.; Reynolds, J. C. L. Science 1988, 241, 67. (10) Luty, B. A.; Wade, R. C.; Madura, J. D.; Davis, M. E.; Briggs, J. M.; McCammon, J. A. J. Phys. Chem. 1993, 97, 233. (11) Kozack, R. E.; d’Mello, M. J.; Subramaniam, S. Biophys. J. 1995, 68, 807. (12) Stouch, T. R.; Williams, D. E. J. Comput. Chem. 1993, 14, 858. (13) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269. (14) Buckingham, A. D.; Flower, P. W. J. Chem. Phys. 1983, 79, 6426. (15) Kollman, P. A.; Dill, K. A. J. Biomol. Struct. Dyn. 1991, 8, 1103. (16) Davis, M. E. J. Chem. Phys. 1994, 100, 5149. (17) Davis, M. E. J. Chem. Phys. 1994, 101, 3414. (18) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Kollman, P. A. J. Am. Chem. Soc. 1993, 115, 9620. (19) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical recipes; Cambridge University Press: Cambridge, U.K., 1992; Chapter 11. (20) Reynolds, C. A.; Essex, J. W.; Richards, W. G. J. Am. Chem. Soc. 1992, 114, 9075. (21) Landau, L. D.; Lifshitz, E. M. Electrodynamics of continuous media; Pergamon: New York, 1963; Chapter 1. (22) Davis, M. E.; Madura, J. D.; Luty, B. A.; McCammon, J. A. Comput. Phys. Commun. 1990, 62, 187. (23) Jorgensen, W.; Tirado-Rives, J. J. Am. Chem. Soc. 1988, 110, 1657. (24) Mohan, V.; Davis, M. E.; McCammon, J. A.; Pettitt, B. M. J. Phys. Chem. 1992, 96, 6428. (25) van Gunsteren, W. F.; Berendsen, H. J. C.; Hermans, J.; Hol, W. G. J.; Postma, J. P. M. Proc. Natl. Acad. Sci. U.S.A. 1983, 80, 4315. (26) Wlodawer, A.; Pavlovsky, A.; Gustchina, A. FEBS Lett. 1992, 309, 59. (27) Mu¨ller, T.; Sebald, W.; Oschkinat, H. J. Mol. Biol. 1994, 237, 423. (28) Brandhuber, B. J.; Boone, T.; Kenney, W. C.; McKay, D. B. Science 1987, 238, 1707. (29) Jackson, J. D. Classical electrostatics, 2nd ed.; Wiley: New York, 1975; Chapter 4. (30) McCammon, J. A.; Northrup, S. H.; Karplus, M.; Levy, R. M. Biopolymers 1980, 19, 2033.

JP953109F

Effective Charges for Macromolecules in Solvent - The Journal of

Recommend Documents