Molecular thermodynamics for charged biomacromolecules

This article was originally published in a journal published by Elsevier, and the attached copy is provided by Elsevier for the author’s benefit and for the benefit of the author’s institution, for non-commercial research and educational use including without limitation use in instruction at your institution, sending it to specific colleagues that you know, and providing a copy to your institution’s administrator. All other uses, reproduction and distribution, including without limitation commercial reprints, selling or licensing copies or access, or posting on open internet sites, your personal or institution’s website or repository, are prohibited. For exceptions, permission may be sought for such use through Elsevier’s permissions site at: http://www.elsevier.com/locate/permissionusematerial

Fluid Phase Equilibria 241 (2006) 317–333

Molecular thermodynamics for charged biomacromolecules Jianzhong Wu ∗ , Dimitrios Morikis ∗ Department of Chemical and Environmental Engineering, University of California, Riverside, CA 92521-0444, United States Received 16 October 2005; received in revised form 28 December 2005; accepted 29 December 2005

Abstract Significant progress has been made over the past 25–30 years toward understanding the molecular basis of stability and biological functions of proteins and nucleic acids. Electrostatic calculations hold a prominent place among most advanced and popular computational methods to describe the influence of solution conditions on and the precise role played by individual amino acids or atoms in stability, folding and molecular recognition. In this article, we present a tutorial overview of coarse-grained and atomistic thermodynamic methods commonly used in modeling electrostatic interactions that are prevalent among biological processes including structural transitions of biomacromolecules and ligand–receptor associations. Illustrative examples are discussed on applications of these thermodynamics models to computational design of proteins with improved function and stability, to charge-transfer equilibria, to structure transitions of nucleic acids and proteins, to protein–ligand interactions, and to ion transport through lipid membranes. Some perspectives are also given on the future trends of computational modeling in biological thermodynamics. © 2006 Elsevier B.V. All rights reserved. Keywords: Computational biology; Electrostatics; Protein–protein interactions; Protein association; Biomacromolecular stability

1. Introduction In a departmental colloquium presented at the University of California at Riverside in November 2004, Professor John Prausnitz urged an audience of chemical and environmental engineers, the authors of this article in particular, to be like the bee instead of the ant or spider. The analogy comes from a quotation of Sir Francis Bacon discussed in the preface of the 2nd ed. of Molecular Thermodynamics of Fluid-Phase Equilibria [1]. Experimenters are like ants; they only collect and use. Theoreticians resemble spiders that make cobwebs out of their own substance. But bees gather materials from the flowers and digest them by their own power. Professor Prausnitz insists (and exemplifies by himself) that a molecular thermodynamicist should be like the bee, seeking useful solutions to real world problems by a delicate combination of rigorous statistical–mechanical theories with adequate semi-empirical approximations. In traditional chemical engineering, molecular thermodynamics is concerned primarily with a single problem, i.e., to

∗ Corresponding authors. Tel.: +1 951 827 2413; fax: +1 951 827 5696 (J. Wu), Tel.: +1 951 827 2696; fax: +1 951 827 5696 (D. Morikis). E-mail addresses: [email protected] (J. Wu), [email protected] (D. Morikis).

0378-3812/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.fluid.2005.12.044

find the compositions of two (or occasionally, more) coexisting phases at various combinations of known and unknown thermodynamic variables. Toward that end, concepts from molecular physics, statistical mechanics, and more recently molecular simulations, have been integrated into equations of state and excess Gibbs-energy models that provide effectual correlations/predictions of the thermodynamic properties of fluid mixtures, fugacity coefficients in particular. Because phaseequilibrium calculations are essential in optimization of separation processes that are routinely used in the petrochemical industry, the fugacity-coefficient models consist of a cornerstone of traditional chemical engineering. In recent years, however, conventional vapor–liquid and liquid–liquid equilibria have faded away from the center stage of chemical engineering, yielding to burgeoning fields mostly labeled with bio or nano. Instead of reactivity and phase behavior of petroleum fluids and natural gases, chemical engineers today are mostly concerned with the microscopic structures and thermophysical properties of complex fluids, i.e., solutions of polymers, surfactants, nanoparticles, and biomacromolecules, that can be deployed for synthesis of advanced materials, for design and delivery of therapeutic drugs, for fabrication of smart sensors, and for development of various environmentally friendly and energy-efficient chemical and biological products and processes. Applications of thermodynamics are not lim-

318

J. Wu, D. Morikis / Fluid Phase Equilibria 241 (2006) 317–333

ited to correlations and predictions of thermophysical properties and vapor–liquid phase diagrams. Indeed, thermodynamics also finds new applications in formulation of experimental systems, in control of self-assembled microscopic structures, and in prediction of materials performance. In this article, we discuss a few applications of molecular thermodynamics to charged biomacromolecules. Because the subject is so vast, no attempt has been made the discussion exhaustive. Instead, attention is given to representative biological systems where electrostatic interactions play a pivotal role in determining structural and thermodynamic properties. To introduce the basic concepts of electrostatic calculations, we begin with discussions of the coarse-grained models of biomacromolecules and predictions based on the Poisson–Boltzmann (PB) equation. We then discuss the limitations of the PB-based methods and recent progress toward improvement. Next, we introduce computational methods from an atomic or molecular perspective [2] and in certain instances, in combination with molecular dynamics (MD) simulations [3,4]. Some pedagogical examples are provided on applications of electrostatic calculations to charge-transfer equilibria, to structural transitions of proteins and nucleic acids, and to protein–protein or protein–ligand interactions. We conclude this article with some perspectives on the future computational research in biological thermodynamics. 2. Simple electrostatic models Under physiological conditions, virtually all biomacromolecules are charged and their properties and biological functions are closely related to the surface electrostatic potential and interactions with surrounding ions [5]. For example, nucleic acids bear a negative charge (PO4 − ) for each base, making them among the strongest natural polyelectrolytes. More than 20% of amino-acid residues in a globular protein are ionizable in an aqueous solution; these include contributions from the carboxyl groups in the side chains of aspartic and glutamic acids and from the basic groups located at lysine, arginine and histidine residues. The vast majority of biological membranes entail certain percentages of negatively charged phospholipids. Furthermore, a typical biological solution contains ample ions such as Na+ , K+ , Mg2+ , Cl− , CH3 COO− , PO4 − and derivatives of adenosine such as ATP and ADP. Electrostatic interactions are crucial to the structure and function of biomacromolecules, ranging from enzyme catalysis, ligand binding and the finetuning of redox potentials, to the stability of folded proteins and the translocon-mediated integration of transmembrane segments into the endoplasmic reticulum [6]. Coarse-grained models provide a first-step toward understanding electrostatic interactions in biological systems. In these highly simplified models, biomacromolecules are depicted as uniformly charged objects of simple geometry, namely proteins as spheres or spheroids, nucleic acids as cylinders, and lipid membranes as planar surfaces; small ions are represented by point charges or charged hard spheres, and the solvent by a dielectric continuum [7]. Within the “primitive” model, the electrostatic potential and ionic distributions near a biomacro-

molecular surface can be described by the Poisson–Boltzmann (PB) equation. Specifically, the Poisson equation connects the electrostatic potential ψ(r) to the local charge density Q(r): ∇ 2 ψ(r) = −

Q(r) , ε0 εr

(1)

where ε0 = 8.854 × 10−12 C2 J−1 m−1 is the permittivity of free space, and εr is the dielectric coefficient or relative permittivity of the solvent. The local charge density, Q(r) = i Zi eρi (r), is in turn determined by the Boltzmann law for the spatial distributions of ions: ρi (r) = ρi0 exp[−Zi eβψ(r)],

(2)

where ρi0 is the ion concentration at zero electrostatic potential, and β = 1/kB T. While the Poisson equation is exact in the framework of classical electrostatics, the Boltzmann law for the distribution of small ions is only approximate; it neglects the ion sizes, non-Coulomb interactions, and correlations in charge distributions [8,9]. In conjunction with an appropriate boundary condition for the surface electrostatic potential or charge density, the electrostatic potential ψ(r) and ionic distribution function ρi (r) can be solved from Eqs. (1) and (2) by using various numerical procedures [10]. Despite the simplicity of the PB equation, an analytical solution is limited to a few special circumstances, such as systems containing only one-type of ions or if the exponential function in the Boltzmann law is approximated by a linear expansion. In the latter case, the PB equation becomes ∇ 2 ψ(r) = κ2 ψ(r),

(3) 1/2

0 2 2 is the Debye-screening where κ = i ρi Zi e β/(ε0 εr ) parameter [11]. Table 1 presents asymptotic solutions to the linearized PB equation in planar, spherical, and cylindrical geometries [12]. In Appendix A, we give analytical solutions to Eq. (3) in more details and discuss important length scales affiliated with electrostatic interactions [7,13]. With the advent of modern computers, the PB equation and its variations have been scrutinized by extensive comparison with molecular simulations [14–16]. Regrettably, the comparison is often unsatisfactory because the PB equation ignores the excluded volumes of ions and correlations in ionic distributions. The situation is most contentious for systems containing multivalent counterions. For example, the PB equation is unable to describe attraction between similar charges and

Table 1 Asymptotic solutions of the linear PB equation in spherical, cylindrical, and planar geometries ϕ(x) Sphere Cylinder Plane

∼e−x /x √ ∼ π e−x / x ∼e−x

Here both the electrostatic potential and the surface distance are expressed in dimensionless units, ϕ = eψ/kB T and x = κr. The proportionality constants are determined by the boundary conditions.


319

molecular models are required to incorporate the atomic details of solvent and ions. Toward that end, DFT holds major advantages over the PB equation not only because it is able to account for ion size and density correlations that are ignored in the PB equation but more important, because it provides a systematic approach to incorporating the specific properties of ions and solvent molecules beyond the primitive model [30]. 3. Cell model

Fig. 1. Surface electrical potential of a charged particle vs. the charge density in an asymmetric electrolyte solution. Here C stands for the molar concentration of the electrolyte with divalent cations and monovalent anions; the spherical particle has a radius 1.5 nm and the diameter of small ions is 0.425 nm. The opened circles are results from Monte Carlo simulations [125,126]; the dashed, dotted, and solid curves are predictions of the PB equation, integral-equation theory, and DFT, respectively, all within the primitive model [25].

charge inversion2 that have been revealed by both molecular simulations and experiments [17–20]. In recent years, a number of more sophisticated statistical–mechanical methods have been proposed that yield much improved results in comparison with simulations [21,22]. The classical density functional theory (DFT) represents a powerful alternative to the Boltzmann equation for ionic distributions [23–27]. In addition to the electrostatic potential, DFT is able to account for non-Coulomb interactions and intermolecular correlations in terms of an excess Helmholtz energy Fex . In general, the distribution of mobile ions is given by [24]: ρi (r) = ρi0 exp[−Zi eβψ(r) − βF ex ].

(4)

Eq. (4) reduces to the Boltzmann law if the excess Helmholtz energy due to ion size and charge correlations is neglected. As a result, the difference between DFT and the PB equation is most significant for systems with strong electrostatic interactions that lead to charge localization and subsequently magnified excluded-volume effects [24]. Fig. 1 shows, for example, the surface electrostatic potential of a charged sphere as a function of the surface charge density predicted from the PB equation, from an integral-equation theory, and from DFT [25]. While the predictions of DFT agree well with the simulation results, the performance of the PB equation is unsatisfactory, in particular in the presence of divalent counterions. The discrete nature of water molecules and non-Coulombic interactions of small ions play an important role in biological specificities of ion-binding proteins, nucleic acids, ion channels, and biomembranes [28,29]. In such cases, more accurate 2

Charge inversion is a non-intuitive electrostatic phenomenon in which a macroion adsorbs an excess amount of counterions such that the surface electrostatic potential is opposite to that of the bare charge.

The cell model was originally proposed for predicting the osmotic coefficients of synthetic or natural polyelectrolyte solutions and the potentiometric titrations of weakly ionizable polyelectrolytes [31,32]. In recent years, it has also been used to describe structural transitions of biomacromolecules, particularly DNA duplex [33–37]. Within the cell model, a DNA molecule is depicted as an infinitely long, rigid rod of uniform surface charge density.3 A DNA solution is represented by uniformly distributed, cylindrical cells containing individual DNA molecules coaxially placed at the cell centers. The electrostatic potential and ionic distributions within each cell are described by the PB equation, which, in the cylindrical geometry, depends only on the radial distance: 1 ∂ ∂ϕ r = −4πlB ρi0 Zi e−Zi ϕ(r) , (5) r ∂r ∂r i

where ϕ = eψ/kB T is the reduced electrostatic potential, lB = e2 /(4πε0 εr kB T) is the Bjerrum length,4 and ρi0 is the number density of ion i at the cell boundary where a zero electrostatic potential is selected. As discussed later, ρi0 can be determined from the electrostatic potential and average concentration of DNA molecules in the solution. With the boundary conditions for the electrostatic potential at the DNA surface,5 ϕ (a) = 2ξ/a, and at the cell boundary ϕ (R) = 0, Eq. (5) can be formally integrated r r ϕ(r) = 2ξ ln(r/R) − 4πb0 Zi ln r ρi (r ) dr R a i

R r (6) + r ρi (r ) ln dr . R r As shown in Fig. 2, a is the radius of the cylindrical rod, b the cylindrical length per unit charge, i.e., the spacing between nearest-neighboring charges along the axis of the cylindrical rod. The parameter ξ = lB /b stands for a linear charge density in dimensionless units; it has a pronounced effect on the radial distribution of small ions and on the salt-dependence of the electrostatic potential. With the Boltzmann equation for ionic 3 At physiological conditions, the persistence length of a double-strand (ds)DNA is about 50 nm and that for a single-strand (ss)-DNA is about 1 nm. 4 The Bjerrum length provides a measure of the distance at which the electrostatic potential two elementary charges is equal to the thermal energy (kB T). For water at 25 ◦ C, the Bjerrum length is 0.714 nm. 5 Gauss’ law states ∂ψ/∂r| r=a = −ζ/(ε0 εr ) where ζ = −e/(2πab) represents the surface charge density for a negatively charged poly electrolyte.

320


Fig. 2. In the cell model, a DNA molecule is represented by an infinitely long, cylindrical rod immersed in a coaxial cylindrical cell containing point charges. Here a stands for the DNA radius, b is the axial distance between successive charges, and a2 /R2 is equal to the DNA volume fraction.

distributions, ρi (r) = ρi0 e−Zi ϕ(r) , the electrostatic potential can be numerically solved using an iteration method [38]. For a DNA solution free of salt, the PB equation can be solved analytically. In that case, the reduced electrostatic potential is given by [32]:

2 2 2 B ln(Ar) ϕ(r) = ln κc r sinh , (7) 2B2 where κc = 4πlB ρc0 Zc2 , subscript “c” denotes counterions (i.e., those ions with charges opposite to that of the DNA), and A and B are parameters obtained from 1 + B coth[B ln(AR)] = 0 and 1 + B coth[B ln(Aa)] = ξ. From the boundary condition ϕ(R) = 0, we find the concentration of counterions at the cell boundary: ρc0 1 − B2 = , ρ¯ c 2ξ

(8)

where ρ¯ c is the average concentration of the counterions in the salt-free DNA solution.

Fig. 3. The osmotic coefficient vs. DNA concentration predicted by the cell model (solid and dotted curves correspond to, respectively, calculations based on solid and hollow cylinder models) and from experiments (points). The dashed line is the limiting value of the osmotic coefficient predicted by the counterioncondensation theory [40].

If the osmotic pressure is measured against a salt buffer with concentration ρ0 , the salt concentration in the DNA phase can be estimated from the Donnan equilibrium6 : (1 − B2 )ρ¯ c + ρs ρs = ρ02 . (11) 2ξ In that case, the osmotic coefficient becomes (1−B2 )ρ¯ c 2ξ

+ 2(ρs − ρ0 )

3.1. Osmotic pressure

φ=

The osmotic pressure of a DNA solution can be approximated by the summation of that for a salt-free DNA and that for the simple electrolyte solution [32]. The former can be calculated from the osmotic pressure at the cell boundary where the electric field vanishes, i.e., Π/kB T = ρc0 = ρ¯ c (1 − B2 )/2ξ; the latter can be expressed in terms of the average salt concentration ρs and the osmotic coefficient φs for a simple electrolyte. The overall osmotic pressure is given by

At low DNA concentration and free of salt, a/R and subsequently B vanishes. The osmotic coefficient is reduced to φ∞ = 1/(2ξ), which is the same as that predicted by the counterion-condensation theory discussed in Appendix A [39]. Fig. 3 shows the osmotic coefficients predicted by the cell model in comparison with experimental data [40]. In these calculations, both the DNA radius (a = 1 nm) and the spacing between two nearest phosphate charges (b = 0.17 nm) are estimated from independent measurements. The solid lines shown in Fig. 3 are calculated from Eq. (12). The dashed lines, providing slightly improved results, are from a refined cell model where the DNA duplex chains are treated as hollow cylinders, taking into account

Π (1 − B2 )ρ¯ c = + 2ρs φs . kB T 2ξ

(9)

The osmotic coefficient is defined as the osmotic pressure of the DNA solution divided by that of an ideal solution: Π φ ≡ id ≈ Π

(1−B2 )ρ¯ c 2ξ

+ 2ρs φs

ρ¯ c + 2ρs

6

.

(10)

ρ¯ c

.

(12)

Donnan equilibrium requires an equal mean activity of the salt across the semi-permeable membrane. Eq. (11) is obtained by replacing the ionic activity with number density at the cell boundary.


321

Fig. 4. The melting-point temperature for the helix-coil transition of a double-strand T4 phage DNA as a function of NaCl concentration (CNa , mM) in presence of MgCl2 . Experimental points are marked for different ratios of the magnesium to DNA concentrations, CMg /CDNA with CDNA = 0.042 mM. The curves are calculated ˚ ah = 9.5 A, ˚ for ds-DNA and bc = 3.38 A, ˚ ac = 7.0 A ˚ for from the PB theory (A) and from the CC theory (B) [37]. The parameters used in PB theory are bh = 1.69 A, ˚ dMg = 3 A; ˚ in CC theory are bh = 1.69 A ˚ and ξ h = 4.2 for ds-DNA and bc = 4.5 A, ˚ ξ h = 1.6 for ss-DNA. ss-DNA, dNa = dCl = 2 A,

the major and minor grooves in the ␣-helical structure.7 Over a wide range of DNA concentrations, both the original and modified cell models predict the osmotic coefficients in good agreement with experimental data. 3.2. Stability of DNA duplex In development of new medicinal strategy to target nucleic acids, thermal denaturation experiments are commonly carried to access the influence of drugs on the stability of the doublestrand structure [41]. The denaturation curve reflects the double helix to single-strand equilibrium and is generally analyzed using a “two-state” model. Each nucleic acid is considered to be either totally in the double helix form or totally dissociated. The melting temperature (Tm ) is defined as the temperature at which half of the nucleic acids are in the single-strand form. Experimental studies show that the melting temperature is mainly determined by the charge densities of double-strand (ds) and single-strand (ss) DNA chains. Furthermore, this transition is highly sensitive to the salt concentration (Cs ), or more precisely, the mean activity (a± ) and charges of cations in the solution. Electrostatic models have been used to account for the variation of DNA duplex stability with respect to solution conditions, in particular the types and concentrations of salts. The Gibbs energy of the helix-coil transition for a ds-DNA can be effectively expressed in terms of the electrostatic and nonelectrostatic contributions:

Gm = G0 + GC ,

(13)

where G0 accounts for all non-electrostatic contributions, and

GC accounts for the electrostatic contributions. At the melt7

The electrostatic potential according to this refined cell model is given by ϕ0 (r) = 2 ln[κr cos(2 ln(r/RM )] − ln(2c1 ) for a < r < R and ϕi (r) = ϕ0 + 2 ln(l + c2 r2 ) for r < a. The integration constants c1 , RM , c2 and ϕ0 are obtained from the boundary conditions ϕ (0) =ϕ R) = 0, the Gauss law ϕ0 − ϕi = 2ξ/a and the continuity condition ϕi = ϕ0 .

ing point, ds-DNA chains are in equilibrium with the ss-DNA chains, i.e., Gm = 0. Because the non-electrostatic free energy is relatively insensitive to electrolyte conditions, G0 can be obtained from calorimetric data for denaturation of a ds-DNA at a particular salt condition (e.g., NaCl). With an empirical correlation for the non-electrostatic free energy, the cell model then provides a predictive tool for describing the variation of the melting point temperature due to addition of a salt. To a good approximation, one may ignore the change in volume due to the helix-coil transition, i.e., GC ≈ UC − T SC . The electrostatic energy and entropy for both ds- and ss-DNA polyions can be calculated from the electrostatic potential and ionic distributions within the cell model [37]:

UC ϕ(a) 1 R + ρi (r )Zi ϕ(r )2πr dr , (14) = bkB T 2b 2 a i

R SC ρi (r ) 0 ρi (r ) ln − ρ =− (r ) + ρ i i 2πr dr , 0 bkB ρ a i i (15) where the electrostatic potential is given by Eq. (6). Fig. 4 shows a comparison of the theoretical predictions with experiments when the transition occurs in the presence of magnesium ions. While the cell model slightly underpredicts the melting temperature, the opposite is true for the counterion-condensation theory (see Appendix A). Similar results have been obtained when Eqs. (13)–(15) are applied to DNA helix-coil transition in the presence of the natural polyamines such as putrescine2+ , spermidine3+ , spermine4+ , and their synthetic homologs with different spacing between the charged amino groups [42]. 4. Generalized Poisson–Boltzmann equation Computational modeling of biomacromolecules from an atomic or molecular perspective often starts with an atomic resolution of the biomacromolecular structure that can be deter-

322


mined by experimental or computational means.8 Each atom from a biomacromolecule is represented by a sphere embedded in a medium with the dielectric coefficient in the range of 2–20, depending on the computational protocol. The atomic radius is set to the van der Waals radius and when the atom belongs to a polar group, it carries a partial charge. As in simple electrostatic models, the aqueous environment is represented by the primitive model. Near the surface of a biomacromolecule, the electrostatic potential satisfies the generalized Poisson equation: ∇ · [ε0 εr (r)∇ψ(r)] = −Q(r),

(16)

where εr (r) stands for the local dielectric coefficient, the symbol “” denotes gradient when applied to a scalar and divergence when applied to a vector. The inhomogeneous dielectric coefficient εr (r) is fully determined by the 3D-structure of the biomacromolecule and its value in the protein interior is different from that of the solvent. Conversely, the local charge density Q(r) consists of contributions from the biomacromolecule and from small ions. Assuming that the biomacromolecule has a fixed charge distribution denoted by qM (r), we may express the generalized PB equation as ∇ · [ε(r)∇ψ(r)] = −qM (r) − Zi eΓi (r)ρi0 e−Zi eψ(r)kB T , i

(17)

where Γ i (r) represents the accessibility to ion i at position r, typically defined by the van der Waals radii of the surface atoms ˚ for monovalent and a probe sphere of fixed radius (e.g., of 2 A ions) representing the ion. The PB equation can be linearized if the reduced electrostatic potential is small. Nevertheless, it has been shown that in terms of practical applications, the difference between the linear and non-linear PB equations is not appreciable even at large value of the reduced electrostatic potential [43]. As a result, the linear PB equation is popular in modeling biomacromolecules and biological processes including those described below [12,43–47]. Numerical solutions of the generalized PB equation are computationally demanding due to the complex shapes and charge distributions of biomacromolecules. Typically, the values for charges and dielectric coefficients are assigned on a 3D-grid surrounding the biomacromolecule. The electrostatic potential for each grid point is then numerically solved by using space discretization methods [48,49]. The accuracy of the numerical solution depends on the scheme used for the space discretization and on the size of the computational grid. Popular software packages include APBS [50], UHBD [51], DELPHI [52], MEAD [53], WHAT IF [54], and CHARMM [55]; most of these algorithms are based on different versions of finite difference or finite element methods for the solution of the PB equation. 8 As of December 2005, there are 34,303 biomacromolecular structures deposited at the Protein Data Bank (http://www.rcsb.org/pdb/), most of which have been determined by X-ray diffraction and, to a lesser extent, by nuclear magnetic resonance (NMR) methods. In addition, the three-dimensional structures of biomacromolecules can be estimated from computer models using a template from a homologous experimental structure (structural homology method).

The electrostatic potentials of biomacromolecules can be intuitively visualized using software packages such as GRASP [56], MOLMOL [57], SWISS PDB VIEWER [58], and INSIGHT II (Accelrys Software Inc., San Diego, CA). Visualization of the electrostatic properties is a standard tool for analysis of protein structures. The electrostatic potential is projected on the biomolecular surface or alternatively, its spatial distribution is viewed in the form of iso-potential surfaces or lines at a fixed potential value (in units of kB T/e). Such presentations are popular in the publications of newly determined protein structures and complement other types of analyses for structure specificity and stability, and for binding either at a molecular level, such as surface complementarity, or at the individual amino acid level, such as van der Waals interactions, hydrogen bond, and salt bridge formations. Electrostatic calculations are helpful to analyze the electrostatic character of surfaces, volumes, and surrounding space for proteins and nucleic acids, biomacromolecular complexes and assemblies, and binding and active sites [44,50]. For example, a combined molecular dynamics (MD) and electrostatics study has provided a rational explanation for the high catalytic rate of the enzyme acetylcholinesterase [59]. It was shown that acetylcholinesterase generates a strong electrostatic field capable to attract the positively charged substrate acetylcholine to the active site; however, the high catalytic rate of the enzyme was inconsistent with the long and narrow gorge at the active site observed in the crystal structure. The MD/electrostatics study revealed an additional transient opening of a “back door” channel for substrate entry [59]. It was proposed that the dual substrate entry is responsible for the high catalytic rate of the enzyme, which could be tested by specific mutations. Another example is provided by the interaction of the immune system protein C3d with its B and T cell receptor CR2. Site-directed mutagenesis studies proposed a binding site based on loss or enhancement of the binding ability by selected mutants involving charged amino acids. The selection of charged amino acids for mutations was based on several earlier experimental studies which had shown pH and ionic strength dependence for the C3d–CR2 association. However, when the crystal structure of the C3d–CR2 complex became available, the binding site was found elsewhere. This ambiguity was resolved by electrostatic calculations which lead to the proposal of a general model for association driven by electrostatic macro-dipoles (recognition process) and short range interactions (binding process) [60]. This model provides a rational explanation of all available experimental data, including those which were previously thought to be contradictory. Fig. 5 shows the iso-potential contour surfaces for native C3d, two of its mutants with variable binding ability, and CR2. The single mutation K162A alters the balance in favor of the negative electrostatic potential, which correlates with the significant increase of binding ability observed in experiment. Conversely, the triple mutation D36A/E37A/E39A alters the balance in favor of the positive electrostatic potential, which correlates with the significant decrease of binding ability. Application of electrostatic calculations can also be illustrated with the analysis of the highly homologous pox viral proteins VCP (of the vaccinia virus) and SPICE (of the variola virus). The electrostatic potentials were


323

Fig. 5. Iso-potential surfaces at ±1kB T/e for native C3d (A) and its mutants K162A (B) and D36A/E37A/E39A (C). The macrodipole of C3d is responsible for the interaction with the excessively positive receptor CR2 [60]. The mutant with the increased activity has an increased negative electrostatic potential (B) and that with decreased activity has a decreased negative potential (C). Negative and positive potentials are colored in red and blue, respectively.

used to explain the up to 1000-fold difference in their activities against the immune system protein C3b [61,62]. These studies provided the theoretical basis for the design of several mutants of VCP/SPICE with predicted variable activities. A combination of MD simulations and electrostatic calculations for VCP/SPICE and selected mutants was used to demonstrate the effect of correlated inter-modular motions on the spatial distribution of the electrostatic potential [62]. Knowledge of the role of dynamics and electrostatics in VCP–C3b interaction is essential for the development of efficient vaccines against the evasion of the immune system by the pox viruses. Fig. 6 shows projections of calculated electrostatic potentials on the surface of VCP, using structure snapshots from an MD simulation. The figure shows the effect of dynamics (flexibility and mobility) in the variability of the surface electrostatic potential. The electrostatic potentials are also useful to study the regulation of protein–membrane association by ions, for proteins involved in signal transduction and vesicle trafficking. It was shown that Ca2+ -mediated nonspecific electrostatic interactions contribute to the association of C2 protein domains with different phospholipid membrane surfaces [63]. This is possible by the binding of Ca2+ to clusters of aspartic acids in a C2 loop which interacts with the membrane. Some efforts have been made to quantify electrostatic potentials by using a single parameter that provides an easy comparison of molecular interaction fields across homologous proteins or mutants. For example, similarity indices are useful to distinguish molecular properties, such as sequences or electrostatic potentials. The electrostatic similarity index (ESI) between pro-

teins ‘a’ and ‘b’ is defined as ESIa,b = 2(Pa , Pb )/(P2a + P2b ) where (Pa , Pb ), P2a , and P2b are the scalar products of the protein electrostatic potentials at the grid-points (i, j, k), within a specified “skin region” in the space surrounding the proteins, e.g., (Pa , Pb ) = i,j,k ψa (i, j, k)ψb (i, j, k) [64]. Pairwise comparisons of electrostatic similarity indices are possible using principal component analysis or clustering methods, as was the case for the classification of 104 proteins of the Pleckstrin homology domain family [65] and 33 homologous blue copper proteins [66]. A charge density probability function has been introduced to predict conservation of electrostatic properties and enzymatic function within four different enzyme families and one super-family [67]. Long-range electrostatic interactions are responsible for steering the substrate from the solvent toward the catalytic site and the short-range electrostatic interactions are responsible for local proton sharing and transfer during the catalytic process. Such analyses are useful to predict common functional properties of homologous proteins and phylogenetic trees of proteins. 5. Ionization of biomacromolecules Acid–base equilibria, ion-binding, and electron-transfer redox equilibria are important for the functions of biomacromolecules in many biological processes. All these chargetransfer equilibria are critically dependent on electrostatic interactions and can be modeled in essentially the same manner [68–77].

Fig. 6. Projections of electrostatic potentials on molecular surfaces for three MD snapshots of the VCP structure. VCP is made of four modules connected with short and flexible loops. The figure depicts the surface variation of the electrostatic potentials owed to the spatial mobility of the four VCP modules and the internal flexibility of each individual VCP module. The snapshots represent structures at 1, 5, and 8 ns (shown at the bottom panels in ribbon representation) [62]. Negative and positive potentials are colored in red and blue, respectively.

324


We consider first the dissociation of a proton from a single acidic group isolated in the solution (e.g., an aspartic acid): Ka

AH←→A− + H+ ,

(18)

where Ka is the dissociation constant, usually expressed as pKa −log Ka = 0.434 G/kB T, G is the free energy of protonation. The Henderson–Hasselbalch equation relates the pH with the pKa : pH = pKa + log

[A− ] [AH]

(19)

and shows that the pKa is the pH where the concentration of protonated state is equal to that of the deprotonated state. The fractional protonation or protonation probability, f = [AH]/([AH] + [A− ]), is related to pKa by ln

f = 2.303(pKa − pH). 1−f

(20)

We will proceed with our analysis for acid–base equilibria but similar schemes can be applied to ion-binding or electrontransfer redox reactions. For example, in a redox reaction: KET

Aox + e− ←→Ared ,

(21)

the redox potential is given by the Nernst equation: Φ = Φ0 +

kB T [Aox ] ln , F [Ared ]

(22)

where Φ0 = (kB T/F ) ln KET is the potential in which there is an equal concentration of oxidized and reduced states, F is the Faraday constant, and KET is the electron transfer equilibrium constant. Apparently, the Nernst equation is analogous to the Henderson–Hasselbalch equation. We now consider protonation of an ionizable group in a biomacromolecule that bears permanent charges and also other ionizable groups. The pKa for the isolated ionisable group is referred to as the model pKa or pKa0 , and its value in the biomacromolecular environment is referred to as the apparent app pKa or pKa . To a good approximation, pKa0 is the same as that for the protonation of a single acid in solution as discussed earlier. In other words, the values for pKa0 are known from experimental measurements for single ionizable groups (e.g., pKa0 = 4.0 for aspartic acid in water). Because of the permanent background charges (or partial charges), desolvation effects, and app other ionizable groups of the biomacromolecule, pKa is in 0 general different from pKa . The presence of other ionizable groups in a biomacroapp molecule makes the calculation of pKa complicated because at a given pH, the ionization states of neighboring ionizable groups are interdependent. To circumvent this problem, a pHindependent hypothetical quantity has been proposed, called intrinsic pKa or pKaintr [78]. The intrinsic pKa describes the ionization process of a specific ionizable group when all other ionizable groups are in their neutral states. The apparent pKa for an acid is related to pKaintr and a free energy reflecting the

Fig. 7. Thermodynamic cycle used to calculate the intrinsic pKa from experimentally known model pKa0 and electrostatic calculations. The subscripts “M” and “m” denote the ionizable group in the biomacromolecule and in the isolated state, respectively.

pH-dependent interactions of this ionizable group with all other ionizable groups of the biomacromolecule: pKaapp = pKaintr +

Ginter . 2.303kB T

(23)

Once pKaintr values are calculated, they are used to determine the pairwise interactions among all ionizable groups, Ginter , using statistical–mechanical methods as discussed below to account for all possible ionization states of the biomacromolecule. In app general, evaluation of pKa requires the titration curve of the system rather than direct use of Eq. (23). Fig. 7 shows a thermodynamic cycle that allows us to compute the pKaintr relative to pKa0 of the same group isolated in a solution. The upper horizontal process describes the dissociation of an isolated ionisable group (e.g., an aspartic acid in solution) in the high dielectric aqueous solution, whereas the lower horizontal process describes the dissociation of the same ionizable group in the biomacromolecular environment. The vertical processes describe “desolvation” effects, i.e., the transfer of the protonated and deprotonated states from an aqueous solution to a biomolecular environment. The lower horizontal process in Fig. 7 provides the pKaintr : pKaintr = pKa0 −

Z

Genv , 2.303kB T

(24)

where Z takes −1 for acids and +1 for bases. The term

Genv describes the interactions with the environment, i.e., solvent and background charges, for the ionizable group when all other ionizable groups are in their neutral states

Genv = Ga − G0a = GA− − GAH .

(25)

The free energy double difference (Fig. 7 and Eq. (25)) consists of two contributions:

Genv =

GB +

GC ,

(26)

B where

GB = GB − GB AH with Gi being due to the A− change of the dielectric environment for i = A− or AH, and similarly

GC is the corresponding change due to the interactions with the background charges. For either the protonated or deprotonated state, the single-difference term GB is the change in


reaction field (or Born) energy: Ze B B (ψ − ψm ), 2 M

GB =

(27)

where ψM and ψm are the reaction field potentials, owed to induced polarization of charges in the surrounding environment, at the position of the charged group. The term GC is the change in the Coulombic interactions with background charges: C C

GC = Ze(ψM − ψm ),

(28)

C and ψ C are the Coulombic potentials at the position where ψM m of the charged group. Because of the background charges, the changes of the electrostatic energies apply to both the protonated and deprotonated states of the ionizable group. In Eq. (23), Ginter is not calculated directly. Instead, it is estimated by the interaction potential between different ionic groups at all possible ionization states. The ionization microstates of a biomacromolecule with N ionizable groups can be described by n = 2N vectors, δn = {δn1 , δn2 , . . . , δnN }, with each element δni taking value 0 or 1 for the protonated or deprotonated ionizable group i, respectively. The free energy difference of the nth microstate from a fictitious microstate with all the ionizable groups in their neutral state is given by [77]: ⎧ ⎫ N ⎨ ⎬

Gn = 2.303kB TZi δni (pH − pKiint )+δni δnj Gij , ⎩ ⎭ i=1

325

A small protein of 100 ionizable amino acids can assume 2100 ionization states, which makes an exact calculation of Eqs. (29)–(33) computationally impossible. To overcome the 2N computational challenge, several approximations have been proposed to reduce the number of meaningful interactions, such as the reduced site approximation [74], the hybrid statistical mechanical/Tanford–Roxby approximation [77], the clustering method [75], and a Monte Carlo method [76]. The titration curve for each ionizable group is described by plotting the approximated qi as a function of pH and for the whole biomacromolecule by plotting the approximated Q as a function of pH. The apparent pKa value for ionizable group i is equal to the pH value at which the group i is 50% charged. The pKa values predicted by these methods have been successfully tested with experimental results [71–73,77,79–83]. app Analysis of pKa values has been used to identify electrostatic interactions which contribute to the structural stability of D/D PKA RII␣ [84], the catalytic mechanism of the protein GART [85], the coil-helix transition of GART [86], and the association of C3d with CR2 [60]. Fig. 8 illustrates the effect of app the protein environment on pKa values. It shows the shifts in app the pKa values of ionizable amino acids compared to model

1≤j 1/κ, the overall electrostatic interaction may be considered unimportant from a practical point of view. A conventional wisdom is that the electrostatic interactions between charged species are insignificant when their separation is beyond the Debye screening length. Eq. (A9) provides the essential basis of the DLVO11 theory for describing the electrostatic contribution to colloidal forces. A.2. Electric double layer When a charged plate is surrounded only by the counterions, the non-linear PB equation can be solved analytically. In this

331

case, the PB equation in reduced units is expressed as ∂2 ϕ = −4πlB ρc0 Zc e−Zc ϕ(r) , ∂r 2

(A10)

where r is the perpendicular distance from the surface, ρc0 the counterion density at the surface, and Zc is the counterion valence. The density of counterions at contact can be obtained from the contact-value theorem [22]: kB Tρc0 =

εr ε0 E2 , 2

(A11)

where E = −ψ (0) is the electric field at the surface. Using Gauss’ law ψ (r = 0) = σ/ε0 ε, we have ρc0 = 2πlB (σ/e)2 with σ standing for the surface charge density. With the boundary conditions ψ(0) = 0 and ψ (r = 0) = σ/ε0 ε, Eq. (A10) can be solved analytically: r 2 ln λGC +1 ϕ(r) = , (A12) Zc where λGC = (2πZc lB σ/e)−1 is the Gouy–Chapman length. Because the counterion–wall electrostatic energy is given by u(r) = kB Tr/λGC , the Gouy–Chapman length provides a measure of the distance at which the thermal energy equals the counterion–wall interaction energy. Substituting Eq. (A12) into the Boltzmann equation for counterion distribution gives ρ(r) =

ρc0 r

λGC

2 .

(A13)

+1

From Eq. (A13), one can find that the Gouy–Chapman length also corresponds to a distance from the charged surface that contains half of the counterions. A.3. Counterion-condensation (CC) theory The concept of counterion condensation stems from consideration of an infinitely long rod with a linear charge density ξ = lB /b surrounded by counterions of valence Zi . The system serves as a limiting condition for a salt-free polyelectrolyte at infinite dilution. The unscreened Coulomb interaction between the charged rod and a counterion at a distance r from the rod is r βu(r) = 2|Zi |ξ ln (A14) R where R is a reference position where the electrostatic potential is defined as zero. With the assumptions that the rod has negligible radius and all counterions are independent to each other, the partition function of each counterion within a cylindrical cell of radius R is given by

R

R 1 2 −u(r)/kB T Ω= e 2πr dr = 2−2|Z |ξ r 1−2|Zi |ξ dr. i πR2 0 R 0 (A15)

11

DLVO theory is named after Russian physicists B.V. Derjaguin and L.D. Landau and Dutch scientists E.J.W. Verwey and J.Th.G. Overbeek who constructed the theory independently.

Because the integration starts from zero, Ω diverges when |Zi |ξ ≥ 1. In order to avoid such divergence, Manning argued

332


that in this case, the ion atmosphere is unstable and some counterions must condense to the polyion such that its linear density is reduced to |Zi |ξ = 1 [39]. As a result, a fraction of the fixed charges, 1 − (|Zi |ξ)−1 , becomes completely neutralized by the counterions. Manning’s argument coincides with an independent work by Oosawa who proposed a two-state model for the distribution of counterions around a charged rod, i.e., a bound state and a free state [131]. In the PB equation, the distribution of counterions is continuous and there is no bound state. Although the instability condition is valid only for an infinitely thin line charge in the limit of zero concentration, the counterioncondensation (CC) theory presumes the prevalent applicability of the “limiting law”, i.e., |Zi |ξ ≤ 1 for all linear polyelectrolytes. It further assumes that the distributions of uncondensed ions and electrostatic potential can be calculated from the linear PB equation. Clearly the divergence of the partition function can be avoided if a more realistic model of polyelectrolyte is considered [132]. For this reason, the CC theory has been controversial from the outset and subjected to a number of modifications, in particular on the spatial distribution of condensed counterions and their division from the uncondensed ions [133–135]. Nevertheless, it remains a popular choice in practical applications not only because of its simplicity but also, probably more important, because of its good agreement with experiments for a number of important properties of polyelectrolytes solutions including those containing DNA. References [1] J.M. Prausnitz, R.N. Lichtenthaler, E.G.D. Azevedo, Molecular Thermodynamics of Fluid-Phase Equilibria, Prentice-Hall PTR, Upper Saddle River, NJ, 1999. [2] M.K. Gilson, Biophysics Textbook Online, 2000. [3] J.A. McCammon, S.C. Harvey, Dynamics of Proteins and Nucleic Acids, Cambridge University Press, Cambridge, UK, 1987. [4] C.L.I. Brooks, M. Karplus, B. Montgomery Pettitt, Protiens, A Theoretical Perspective of Dynamics, Structure, and Thermodynamics, John Wiley & Sons, New York, 1988. [5] R. Cotterill, Biophysics: An Introduction, Wiley, Chichester, 2002. [6] M.E. Davis, J.A. McCammon, Chem. Rev. 90 (1990) 509–521. [7] J.N. Israelachvili, Intermolecular and Surface Forces, Academic Press, London, 1992. [8] R.R. Netz, H. Orland, Eur. Phys. J. E 1 (2000) 203–214. [9] V. Vlachy, Annu. Rev. Phys. Chem. 50 (1999) 145–165. [10] J.E. Sader, D.Y.C. Chan, Langmuir 16 (2000) 324–331. [11] Y. Levin, Rep. Prog. Phys. 65 (2002) 1577–1632. [12] F. Fogolari, A. Brigo, H. Molinari, J. Mol. Recogn. 15 (2002) 377–392. [13] W.B. Russel, D.A. Saville, W.R. Schowalter, Colloidal Dispersions, Cambridge University Press, Cambridge, UK, 1999. [14] H. Boroudjerdi, R.R. Netz, J. Phys.: Conden. Matter 17 (2005) S1137–S1151. [15] A. Naji, S. Jungblut, A.G. Moreira, R.R. Netz, Phys. A: Stat. Mech. Appl. 352 (2005) 131–170. [16] J.P. Hansen, H. Lowen, Annu. Rev. Phys. Chem. 51 (2000) 209–242. [17] K. Besteman, M.A.G. Zevenbergen, H.A. Heering, S.G. Lemay, Phys. Rev. Lett. 93 (2004). [18] S. Ravindran, J.Z. Wu, Langmuir 20 (2004) 7333–7338. [19] M. Quesada-Perez, E. Gonzalez-Tovar, A. Martin-Molina, M. LozadaCassou, R. Hidalgo-Alvarez, Chem. Phys. Chem. 4 (2003) 235–248. [20] A.Y. Grosberg, T.T. Nguyen, B.I. Shklovskii, Rev. Modern Phys. 74 (2002) 329–345.

[21] C.W. Outhwaite, M. Molero, L.B. Bhuiyan, J. Chem. Soc., Faraday Trans. 89 (1993) 1315–1320. [22] D. Henderson, L. Blum, J. Chem. Phys. 69 (1978) 5441–5449. [23] Z.D. Li, J.Z. Wu, Macromol. Symp. 219 (2005) 51–57. [24] Z.D. Li, J.Z. Wu, Phys. Rev. E 70 (2004) (Art. No. 031109). [25] Y.X. Yu, J.Z. Wu, G.H. Gao, J. Chem. Phys. 120 (2004) 7223–7233. [26] Y.X. Yu, J.Z. Wu, G.H. Gao, Chin. J. Chem. Eng. 12 (2004) 688–695. [27] K. Wang, Y.X. Yu, G.H. Gao, Phys. Rev. E 70 (2004). [28] M. Bostrom, D.R.M. Williams, B.W. Ninham, Europhys. Lett. 63 (2003) 610–615. [29] M. Karplus, J.A. McCammon, Nat. Struct. Biol. 9 (2002) 646–652. [30] J.Z. Wu, AIChE J., 2006, in press. [31] R.M. Fuoss, S. Lifson, A. Katchalsky, PNAS 37 (1951) 579–589. [32] A. Katchalsky, Pure Appl. Chem. 26 (1971) 327–373. [33] S. Nilsson, L. Piculell, Macromolecules 24 (1991) 3804–3811. [34] S. Nilsson, L. Piculell, Macromolecules 23 (1990) 2776–2780. [35] S. Nilsson, L. Piculell, Macromolecules 22 (1989) 3011–3017. [36] S. Nilsson, L. Piculell, B. Jonsson, Macromolecules 22 (1989) 2367–2375. [37] N. Korolev, A.P. Lyubartsev, L. Nordenskiold, Biophys. J. 75 (1998) 3041–3056. [38] N. Korolev, A.P. Lyubartsev, A. Rupprecht, L. Nordenskiold, Biophys. J. 77 (1999) 2736–2749. [39] G.S. Manning, J. Chem. Phys. 51 (1969) 924–933. [40] P.L. Hansen, R. Podgornik, V.A. Parsegian, Phys. Rev. E 6402 (2001) (Art. No. 021907). [41] C.H. Spink, S.E. Wellman, Drug: Nucl. Acid Interactions 340 (2001) 193–211. [42] N. Korolev, A.P. Lyubartsev, L. Nordenskiold, Biophys. Chem. 104 (2003) 55–66. [43] F. Fogolari, P. Zuccato, G. Esposito, P. Viglino, Biophys. J. 76 (1999) 1–16. [44] B. Honig, A. Nicholls, Science 268 (1995) 1144–1149. [45] K.A. Sharp, B. Honig, J. Phys. Chem. 94 (1990) 7684–7692. [46] B. Honig, K. Sharp, A.S. Yang, J. Phys. Chem. 97 (1993) 1101– 1109. [47] J.D. Madura, M.E. Davis, M.K. Gilson, R.C. Wade, B.A. Luty, J.A. McCammon, in: D.B. Boyd (Ed.), Reviews in Computational Chemistry, vol. V, VCH Publishers, Inc., New York, 1994, pp. 229–267. [48] N.A. Baker, D. Sept, M.J. Hoist, J.A. McCammon, Ibm J. Res. Dev. 45 (2001) 427–438. [49] N.A. Baker, Meth. Enzymol. 383 (2004) 94–118. [50] N.A. Baker, D. Sept, S. Joseph, M.J. Hoist, J.A. McCammon, Proc. Natl. Acad. Sci. U.S.A. 98 (2001) 10037–10041. [51] J.D. Madura, J.M. Briggs, R.C. Wade, M.E. Davis, B.A. Luty, A. Ilin, J. Antosiewicz, M.K. Gilson, B. Bagheri, L.R. Scott, J.A. McCammon, Comput. Phys. Commun. 91 (1995) 57–95. [52] A. Nicholls, B. Honig, J. Comput. Chem. 12 (1991) 435–445. [53] D. Bashford, in: M. Tholburn (Ed.), An object-oriented programming suite for electrostatic effects in biological molecules, vol. 1343, Springer-Verlag, Berlin, 1997, pp. 233–240. [54] G. Vriend, J. Mol. Graph. 8 (1990) 52. [55] B.R. Brooks, R.E. Bruccoleri, B.D. Olafson, D.J. States, S. Swaminathan, M. Karplus, J. Comput. Chem. 4 (1983) 187–217. [56] A. Nicholls, K.A. Sharp, B. Honig, Proteins: Struct. Funct. Genet. 11 (1991) 281–296. [57] R. Koradi, M. Billeter, K. Wuthrich, J. Mol. Graph. 14 (1996) 51–55. [58] N. Guex, M.C. Peitsch, Electrophoresis 18 (1997) 2714–2723. [59] M.K. Gilson, T.P. Straatsma, J.A. McCammon, D.R. Ripoll, C.H. Faerman, P.H. Axelsen, I. Silman, J.L. Sussman, Science 263 (1994) 1276–1278. [60] D. Morikis, J.D. Lambris, J. Immunol. 172 (2004) 7537–7547. [61] G. Sfyroera, M. Katragadda, D. Morikis, S.N. Isaacs, J.D. Lambris, J. Immunol. 174 (2005) 2143–2151. [62] L. Zhang, D. Morikis, Biophys. J. (2006), in press. [63] D. Murray, B. Honig, Mol. Cell 9 (2002) 145–154. [64] R.C. Wade, R.R. Gabdoulline, F. De Rienzo, Int. J. Quant. Chem. 83 (2001) 122–127.

J. Wu, D. Morikis / Fluid Phase Equilibria 241 (2006) 317–333 [65] N. Blomberg, R.R. Gabdoulline, M. Nilges, R.C. Wade, Proteins: Struct. Funct. Genet. 37 (1999) 379–387. [66] F. De Rienzo, R.R. Gabdoulline, M.C. Menziani, R.C. Wade, Protein Sci. 9 (2000) 1439–1454. [67] D.R. Livesay, P. Jambeck, A. Rojnuckarin, S. Subramaniam, Biochemistry 42 (2003) 3464–3473. [68] K. Sharp, in: E. Di Cera (Ed.), Thermodynamics in Biology, Oxford University Press, New York, 2000. [69] G.M. Ullmann, E.W. Knapp, Eur. Biophys. J. Biophys. Lett. 28 (1999) 533–551. [70] G.M. Ullmann, J. Phys. Chem. B 104 (2000) 6293–6301. [71] J. Antosiewicz, J.A. McCammon, M.K. Gilson, J. Mol. Biol. 238 (1994) 415–436. [72] J. Antosiewicz, J.A. McCammon, M.K. Gilson, Biochemistry 35 (1996) 7819–7833. [73] D. Bashford, M. Karplus, Biochemistry 29 (1990) 10219–10225. [74] D. Bashford, M. Karplus, J. Phys. Chem. 95 (1991) 9556–9561. [75] M.K. Gilson, Proteins: Struct. Funct. Genet. 15 (1993) 266–282. [76] P. Beroza, D.R. Fredkin, M.Y. Okamura, G. Feher, Proc. Natl. Acad. Sci. U.S.A. 88 (1991) 5804–5808. [77] A.S. Yang, M.R. Gunner, R. Sampogna, K. Sharp, B. Honig, Proteins: Struct. Funct. Genet. 15 (1993) 252–265. [78] C. Tanford, J.G. Kirkwood, J. Am. Chem. Soc. 79 (1957) 5333–5339. [79] D. Bashford, D.A. Case, C. Dalvit, L. Tennant, P.E. Wright, Biochemistry 32 (1993) 8045–8056. [80] W.R. Forsyth, J.M. Antosiewiez, A.D. Robertson, Proteins: Struct. Funct. Genet. 48 (2002) 388–403. [81] R.E. Georgescu, E.G. Alexov, M.R. Gunner, Biophys. J. 83 (2002) 1731–1748. [82] J.E. Nielsen, J.A. McCammon, Protein Sci. 12 (2003) 1894–1901. [83] J.E. Nielsen, G. Vriend, Proteins: Struct. Funct. Genet. 43 (2001) 403–412. [84] D. Morikis, M. Roy, M.G. Newlon, J.D. Scott, P.A. Jennings, Eur. J. Biochem. 269 (2002) 2040–2051. [85] D. Morikis, A.H. Elcock, P.A. Jennings, J.A. McCammon, Protein Sci. 10 (2001) 2379–2392. [86] D. Morikis, A.H. Elcock, P.A. Jennings, J.A. McCammon, Protein Sci. 10 (2001) 2363–2378. [87] M.J. Ondrechen, J.G. Clifton, D. Ringe, Proc. Natl. Acad. Sci. U.S.A. 98 (2001) 12473–12478. [88] D. Bashford, K. Gerwert, J. Mol. Biol. 224 (1992) 473–486. [89] V.Z. Spassov, H. Luecke, K. Gerwert, D. Bashford, J. Mol. Biol. 312 (2001) 203–219. [90] R.V. Sampogna, A.S. Yang, B.H. Honig, Biophys. J. 66 (1994) A45. [91] A.S. Yang, B. Honig, J. Mol. Biol. 231 (1993) 459–474. [92] C. Tanford, Adv. Protein Chem. 24 (1970) 1–95. [93] J.A. Schellman, Biopolymers 14 (1975) 999–1018. [94] A.H. Elcock, J. Mol. Biol. 294 (1999) 1051–1062. [95] A.H. Elcock, J. Mol. Biol. 284 (1998) 489–502. [96] L. Xiao, B. Honig, J. Mol. Biol. 289 (1999) 1435–1444. [97] Y.N. Vorobjev, H.A. Scheraga, B. Honig, J. Phys. Chem. 99 (1995) 7180–7187. [98] A.S. Yang, B. Honig, J. Mol. Biol. 237 (1994) 602–614. [99] E. Alexov, Eur. J. Biochem. 271 (2004) 173–185. [100] T. Selzer, G. Schreiber, Proteins: Struct. Funct. Genet. 45 (2001) 190–198.

333

[101] T. Selzer, S. Albeck, G. Schreiber, Nat. Struct. Biol. 7 (2000) 537–541. [102] A.H. Elcock, D. Sept, J.A. McCammon, J. Phys. Chem. B 105 (2001) 1504–1518. [103] D. Sitkoff, K.A. Sharp, B. Honig, J. Phys. Chem. 98 (1994) 1978– 1988. [104] N. Froloff, A. Windemuth, B. Honig, Protein Sci. 6 (1997) 1293–1301. [105] Z.S. Hendsch, B. Tidor, Protein Sci. 8 (1999) 1381–1392. [106] V. Zoete, O. Michielin, M. Karplus, J. Comput. Aided Mol. Design 17 (2003) 861–880. [107] G. Archontis, T. Simonson, D. Moras, M. Karplus, J. Mol. Biol. 275 (1998) 823–846. [108] G. Archontis, T. Simonson, M. Karplus, J. Mol. Biol. 306 (2001) 307–327. [109] F.B. Sheinerman, B. Honig, J. Mol. Biol. 318 (2002) 161–177. [110] V.K. Misra, K.A. Sharp, R.A. Friedman, B. Honig, J. Mol. Biol. 238 (1994) 245–263. [111] V.K. Misra, J.L. Hecht, K.A. Sharp, R.A. Friedman, B. Honig, J. Mol. Biol. 238 (1994) 264–280. [112] A.H. Elcock, J.A. McCammon, J. Mol. Biol. 280 (1998) 731–748. [113] L. Zhang, B. Mallik, D. Morikis, submitted for publication. [114] D.A. Doyle, J.M. Cabral, R.A. Pfuetzner, A.L. Kuo, J.M. Gulbis, S.L. Cohen, B.T. Chait, R. MacKinnon, Science 280 (1998) 69–77. [115] B. Roux, Annu. Rev. Biophys. Biomol. Struct. 34 (2005) 153–171. [116] M. LeMasurier, L. Heginbotham, C. Miller, J. Gen. Physiol. 118 (2001) 303–313. [117] J. Aqvist, V. Luzhkov, Nature 404 (2000) 881–884. [118] D.J. Tobias, Curr. Opin. Struct. Biol. 11 (2001) 253–261; J.U. Bowie, Nature 438 (2005) 581–589. [119] G.D. Cymes, Y. Ni, C. Grosman, Nature 438 (2005) 975–980. [120] A.H. Elcock, Curr. Opin. Struct. Biol. 12 (2002) 154–160. [121] R. Mendez, R. Leplae, L. De Maria, S.J. Wodak, Proteins: Struct. Funct. Genet. 52 (2003) 51–67. [122] R. Mendez, R. Leplae, M.F. Lensink, S.J. Wodak, Proteins: Struct. Funct. Bioinform. 60 (2005) 150–169. [123] B. Parker, J. Sohn, G. Buhrman, P. Brown, K. Kristjansdottir, A. Safi, H. Edelsbrunner, W. Yang, J. Rudolph, Biochemistry 44 (2005) 16563–16573. [124] R.R. Gabdoulline, R.C. Wade, Curr. Opin. Struct. Biol. 12 (2002) 204–213. [125] L. Degreve, M. Lozada-Cassou, Mol. Phys. 86 (1995) 759–768. [126] L. Degreve, M. Lozada-Cassou, E. Sanchez, E. Gonzalez-Tovar, J. Chem. Phys. 98 (1993) 8905–8909. [127] B.J. Yoon, S. Kim, J. Colloid Interface Sci. 128 (1989) 275–288. [128] J.P. Hsu, B.T. Liu, J. Colloid Interface Sci. 178 (1996) 785–788. [129] R.E. Rice, J. Chem. Phys. 82 (1985) 4337–4340. [130] K.S. Pitzer, Acc. Chem. Res. 10 (1977) 371–377. [131] F. Oosawa, Polyelectrolytes, Marcel Dekker, New York, 1971. [132] M. Muthukumar, J. Chem. Phys. 120 (2004) 9343–9350. [133] U. Mohanty, B.W. Ninham, I. Oppenheim, Proc. Natl. Acad. Sci. U.S.A. 93 (1996) 4342–4344. [134] M. Deserno, C. Holm, S. May, Macromolecules 33 (2000) 199–206. [135] B. O’Shaughnessy, Q. Yang, Phys. Rev. Lett. 94 (2005) (Art. No. 048302). [136] M.T. Record Jr., W. Zhang, C.F. Anderson, Adv. Protein Chem. 51 (1998) 281–353. [137] D.E. Draper, RNA 10 (2004) 335–343.