HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Citation Map Services Email this article to a friend Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Da Silva, F. L. B. Articles by Penfold, R. Search for Related Content PubMed PubMed Citation Articles by Da Silva, F. L. B. Articles by Penfold, R. Social Bookmarking                   What's this?

A critical investigation of the Tanford-Kirkwood shceme by means of Monte Carlo simulations

¹ Theoretical Chemistry, Chemical Centre, Lund University, S-221 00 Lund, Sweden
² Institute of Food Research, Norwich Research Park, Colney, Norwich NR4 7UA, United Kindgom
³ Group of Biomolecular Physics, Department of Physics, Faculty of Science, UNESP/Bauru, 17033-360 Bauru, São Paulo, Brazil

Reprint requests to: Fernando Luís B. da Silva, Theoretical Chemistry, Chemical Centre, Lund University, POB 124, S-221 00 Lund, Sweden; e-mail: fernando{at}signe.teokem.lu.se; fax: 46-46-222-4543.

(RECEIVED October 10, 2000; FINAL REVISION March 15, 2001; ACCEPTED April 23, 2001)

Article and publication are at http://www.proteinscience.org/cgi/doi/10.1101/.

	Abstract

TOP Abstract Introduction Model systems Results and Discussion Conclusion Materials and methods References

 
Monte Carlo simulations are used to assess the adequacy of the Tanford-Kirkwood prescription for electrostatic interactions in macromolecules. Within a continuum dielectric framework, the approach accurately describes salt screening of electrostatic interactions for moderately charged systems consistent with common proteins at physiological conditions. The limitations of the Debye-Hückel theory, which forms the statistical mechanical basis for the Tanford-Kirkwood result, become apparent for highly charged systems. It is shown, both by an analysis of the Debye-Hückel theory and by numerical simulations, that the difference in dielectric permittivity between macromolecule and surrounding solvent does not play a significant role for salt effects if the macromolecule is highly charged. By comparison to experimental data, the continuum dielectric model (combined with either an approximate effective Hamiltonian as in the Tanford-Kirkwood treatment or with exact Monte Carlo simulations) satisfactorily predicts the effects of charge mutation on metal ion binding constants, but only if the macromolecule and solvent are assigned the same or similar permittivities.

Keywords: Electrostatic interactions; Debye-Hückel; low dielectric cavity; computer simulations; continuum model; proteins model

	Introduction

TOP Abstract Introduction Model systems Results and Discussion Conclusion Materials and methods References

 
It seems reasonable to believe that the binding of a charged ligand to a protein or some other macromolecule will depend, at least in part, on purely electrostatic interactions. All chemically interesting interactions for molecular systems ultimately stem from an electrostatic origin, of course, as only Coulombic potential terms enter the Schrödinger equation. Very often, however, to render a many-body problem tractable, an effective Hamiltonian is introduced where large classes of microscopic interactions are replaced by a suitable average and not accounted for explicitly. Partitioning an effective interaction into an electrostatic (ionic) term, which dominates at long range, and other terms of covalent character with much shorter range is, then, a useful heuristic and computational device. Consider, for example, the titration of a residue side chain in a protein. The free energy change for proton transfer will be broadly similar to that of the corresponding process in the isolated amino acid, but further augmented by an electrostatic term that depends on the overall charge and charge distribution of the protein. In such cases, the binding constant shift arising solely from electrostatic interactions can amount to several pK units.

The role of electrostatic interactions in protein structure and function has long attracted interest. In one of the first theoretical attempts to deal with the problem, Tanford and Kirkwood assumed that the long-range part could be described within a dielectric continuum model (Tanford and Kirkwood 1957). This effective Hamiltonian model has been invoked extensively to study the interaction between charged ligands and proteins, membranes, and other macromolecules (Hill 1955, 1956; Warwicker and Watson 1982; Warshel et al. 1984; Harvey 1986; Bashford et al. 1988) (the reference list given is by no means exhaustive). The original Tanford-Kirkwood (TK) scheme considers the macromolecule as a sphere of low dielectric permittivity different from the surrounding (aqueous) solution of high permittivity. With respect to the rigid macromolecular reference frame supporting fixed charged groups, the salt ions and other charged ligands are supposed to be in relative motion throughout the solvent medium. Tanford and Kirkwood avoided the formidable statistical thermodynamics problem this motion implies by again appealing to the effective interaction construct, eliminating explicit reference to the mobile particles and introducing the Debye-Hückel potential.

Despite the evident success of the TK result, it suffers from two significant limitations: (1) closed form analytical solutions are only available in simple geometric configurations (usually spherical); and (2) nonlinear effects and explicit ion–ion correlations are ignored.

At the cost of computational simplicity, the first constraint has been removed in a number of calculations by several groups (Warwicker and Watson 1982; Bashford et al. 1988; Bashford and Karplus 1990; Beroza et al. 1991; Juffer et al. 1999). The second limitation has only been partially tested by replacing the DH equation with the nonlinear Poisson-Boltzmann equation (Fushiki et al. 1991), though the description remains at the mean-field level ignoring ion–ion correlations. Monte Carlo (MC) simulations have addressed the full statistical mechanical problem but then only for a model with a uniform dielectric permittivity (Fushiki et al. 1991; Kesvatera et al. 1994; Penfold et al. 1998). An additional simplification of the TK prescription, which may prove unrealistic at low ionic strength, lies in the assumption of infinitesimally small macromolecule concentration (Linse et al. 1995). In NMR studies of proteins the typical protein concentration is around 1 mM, whereas the concentration of accompanying counterions can be an order of magnitude larger or more.

The original application of the TK scheme dealt with titrating groups in proteins and how the pK_a for an ionizable amino acid is affected by the rest of the protein. The quantity calculated is the shift in pK_a between the "free" amino acid and the same residue in the protein. It became clear at an early stage that the scheme invoking a low dielectric interior for the protein did not properly reflect the difference in solvation energy between the free amino acid and when it is part of a protein (Tanford and Roxby 1972; Shire et al. 1974; Warshel et al. 1984; Papazyan and Warshel 1997). It seems, however, as there is no consensus of how to best tackle the problem (King et al. 1991; Warshel and Åqvist 1991; Antonsiewicz et al. 1994a, b; Simonson and Perahia 1995; Simonson and Brooks 1996; Löffler et al. 1997; Sham et al. 1997; Baptista and Soares, unpubl.) without invoking additional empirical parameters describing solvation and water penetration (Dwyer et al. 2000). A simple alternative is to assume a uniform dielectric response throughout the solution, including the macromolecule, equal to that of pure water. For titrating groups at or close to the protein surface this latter ansatz seems to give better agreement with experiment (Spassov and Bashford 1998; Kesvatera et al. 1999, 2001).

Another application of the TK scheme is to the binding of a charged ligand to a protein at different conditions, that is, mutation of charged amino acids and/or addition of inert salt. The pK for such a process does not invoke the problematic solvation term mentioned above, and we will use the binding of divalent cations to a negatively charged protein as a model throughout this work. By using MC simulation to resolve the statistical mechanical problem, we will present a critical investigation of the dielectric continuum model originally studied by Tanford and Kirkwood. The simulations will be performed for a model including the dielectric discontinuity and the mobile species in the solution will comprise of salt particles as well as additional counterions in order to maintain electroneutrality.

	Model systems

TOP Abstract Introduction Model systems Results and Discussion Conclusion Materials and methods References

 
We investigate a model system mimicking a globular protein with a net charge, which is supposed to bind divalent ions. It could, for example, be a caricature of calbindin D_9k (Linse et al. 1988), and consists of a rigid sphere with radius R_p and a centrally located charge of valency Z. A surrounding electrolyte solution, including neutralizing counterions, is described by the so-called restricted primitive model (Levesque et al. 1986). Each uniformly charged mobile ion k with valency z_k is treated as a hard sphere of radius R_k; the solvent degrees of freedom are characterized only by an average dielectric permittivity _s that is also assigned to the ionic interior. Internal structural details of the biopolymer are neglected so that the protein interior is represented by a concentric spherical cavity of radius R_d (referred to as the dielectric radius) with a continuum dielectric permittivity _p < _s). Typically, R_d < R_p though this is not a constraint of the model. In any case, max{R_d,R_p} defines an exclusion boundary impenetrable by mobile salt particles and counterions. The entire system is placed in a spherical cell of radius R_c determined by the protein concentration (see Fig. 1).

View larger version (26K): [in this window] [in a new window]   Fig. 1. Schematic representation of the model systems (see also Table 1). A spherical protein in a spherical cell with radii Rp and RC, respectively. The protein interior with a low dielectric permittivity is shown as a shaded region of radius Rd < Rp. A central charge of valency Z and two binding sites marked as black dots.

Classical techniques suggest a straightforward development of the electrostatic potential at any point (because of a fixed charge distribution and a dielectric discontinuity) in an orthogonal polynomial basis (Böttcher 1973), that is, the multipole expansion. This power series is not ideally suited for computer simulations, though, as the number of terms required for satisfactory convergence increases sharply when sources lie very close to the dielectric interface. Nevertheless, if the charges are more than 0.5Å from this boundary, it is feasible to achieve a precision of about 10^-10, with an acceptable number of terms, that is, <1000. Given the elementary charge e and the vacuum permittivity ₀, two ionised sites of valency z_i, z_j (either fixed or mobile) located at radial positions r_i, r_j and subtending an angle _ij, contribute an electrostatic potential energy given by,

(1)

(2)

(3)

where r²_ij = r_i² + r_j² - 2r_ir_j cos _ij if the spatial separation, P_n is the Legendre polynomial of order n, and n_max is the maximum number of terms included in the sum. The first term in Equations 1–3 is the direct Coulomb interaction, whereas the second describes the reaction field accounting for induced polarization charge at the dielectric discontinuity. In the absence of a dielectric inhomogeneity (_p = _s), it is easy to verify that Equations 1–3 reduce to the ordinary Coulomb potential,

(4)

with the characteristic divergence at small separation (r_ij 0). On the other hand, for a nonuniform dielectric response, the reaction field term converges in this limit to yield a nonvanishing contribution, that is, the so-called self-image energy. Together with electrostatic interactions, the simulation Hamiltonian also includes a short-ranged hard-core overlap restriction among the mobile ions (thereby preventing Coulomb collapse in the configurational Markov chain),

(5)

as well as a one-body external field that accounts for the protein excluded volume and imposes the cell boundary constraint acting as a hard wall,

(6)

The full configurational energy of the system then becomes,

(7)

where N = N_c + N_s + N_p is the total number of charges comprising N_c mobile counterions, N_s mobile added salt ions, and N_p fixed protein charges.

The calbindin-like model, introduced above, with a central net charge and two binding sites located at the usual positions, will be denoted SC8 (Spherical model with a Central charge of -8). Three other variations of this model, listed in Table 1, have been considered. A mutation that neutralizes a charged amino acid of the wild type will, in this model, simply lead to a change of the protein net charge, that is, we will refer to this "mutant" as SC7. Slightly more realistic is the Spherical model with a central charge of -7 and a Peripheral charge of -1 close to the protein surface, referred to as SP8 (see Table 1 for further details). One purpose of this work is to investigate the limitations of the DH approximation in a biophysical context. We have, therefore, also included a model protein that is otherwise identical to SC8 but with a variable central charge (ranging from -4 to -24). The most highly charged model system with Z = -24, SC24, has been studied in greater detail.

	Results and Discussion

TOP Abstract Introduction Model systems Results and Discussion Conclusion Materials and methods References

The binding constant shift pK is related to the free energy cost for the protein to incorporate a mobile ion from the solution, see Equations 9 and 10. To calculate pK, therefore, estimates of µ^ex for both free and bound ions are required. The free ion contribution µ_F^ex is usually minor, but can become significant in some circumstances. Figure 2 illustrates the dependence of µ_F^ex on salt concentration c_s for a protein solution at experimental conditions (c_p = 0.02 mM) typically encountered. A moderately charged protein (model SC8) has only a small effect on the free ion chemical potential at low c_s, whereas it is efficiently screened at higher salt concentrations where µ_F^ex closely follows the ionic excess chemical potential of a bulk electrolyte solution. Moreover, without a very substantial net protein charge, µ_F^ex is entirely insensitive to the dielectric properties of the protein interior, as Figure 2 demonstrates, where any shifts on exchanging a solvent filled cavity for a vacuum are below the statistical noise level, regardless of salt content. On the other hand, a highly charged protein (model SC24) dramatically lowers µ_F^ex by several kT in dilute electrolyte. Added salt will screen the strong electrostatic field from the protein, however, effecting an increase in the excess chemical potential with increasing c_s, contrary to the behaviour in a bulk electrolyte solution at least up to physiological concentrations. A dielectric discontinuity now contributes a small but noticeable effect, amounting to a few tenths of kT, as displayed in Figure 2.

View larger version (13K): [in this window] [in a new window]   Fig. 2. The reduced excess chemical potential for a free divalent cation in a 20 µM protein solution is indicated as a function of monovalent salt concentration. The (superimposed) solid lines refer to the SC8 model and the broken lines to the SC24 model. Lines marked by solid circles are calculated with a dielectric boundary radius Rd = 14Å; lines without symbols refer to Rd = 0Å. The thin broken line shows the excess chemical potential of a divalent ion in a bulk salt solution.

View larger version (12K): [in this window] [in a new window]   Fig. 3. The reduced excess chemical potential for a divalent cation as a function of protein concentration with a homogeneous dielectric response (dielectric radius Rd = 0). The solid lines refer to the SC8 model and the broken lines to the SC24 model. Lines marked by circles denote a monovalent simple electrolyte concentration of cs = 100 mM; lines without symbols are for cs = 1 mM. The arrows indicate the value of ßµFex in the corresponding bulk salt solutions. The upper arrow indicates the system at cs = 1 mM (ßµFex = -0.156).

To mimic the calcium binding property of calbindin D_9k, the protein model SC8 binds two divalent ions. The stoichiometric ion binding constant varies with the amount of salt present and provides a suitable test of theoretical predictions for the shift pK with respect to c_s. Figure 4a shows a very significant shift of 3–4 pK units on increasing c_s from 1 mM (reference state) to 500 mM. The simple TK calculation performs very satisfactorily in predicting pK, disagreeing with MC simulation results by less than a few tenths of a pK unit. More significantly, it appears that the presence of a low dielectric protein interior has practically no influence on pK. With a dielectric radius R_d = 12Å, pK is 3.5, whereas a shift of 3.8 is obtained within a uniform dielectric model (see Fig. 4b). Increasing R_d further to 18Å reduces the shift down to 3.2 but because the protein radius is 14Å this corresponds to suggesting that the contact water layer, about 4Å thick, is both impenetrable to salt ions and unpolarizable with a relative dielectric permittivity of unity. Experimental studies have measured pK 4.6 for aqueous solutions of calbindin D_9k on raising the salt concentration from 1 mM to 500 mM (Kesvatera et al. 1994). Although model SC8, with the net charge contracted to a single point and no higher electric moments, is surely too crude a description for any real protein, the salt dependence of cation binding constant shifts for a model with eight randomly placed negative charges is almost indistinguishable from the results obtained with model SC8, as illustrated in Figure 4c.

View larger version (38K): [in this window] [in a new window]   Fig. 4. (a) Total shift in ion binding constant pK, owing to solution ionic strength, is plotted as a function of salt concentration cs for the SC8 model. A reference state at cs = 1 mM is chosen and the protein is allowed to bind two divalent cations. Results from MC simulations and TK calculations are compared as solid and broken lines, respectively. Lines marked by solid circles correspond to a dielectric boundary radius Rd = 18Å and those marked by open squares denote Rd = 12Å; unmarked lines refer to a homogeneous dielectric response (Rd = 0Å). (b) The shift of pK arising from the presence of an interior vacuum pK = pK (Rd = 0Å) - pK (Rd = 14Å) is demonstrated as a function of salt concentration cs. In the solvent region, the relative dielectric permittivity is s = 78.7. Solid line shows MC data and broken line indicates TK results. (c) Total shift in ion binding constant pK, owing to solution ionic strength, is plotted as a function of salt concentration cs for a spherical protein model with eight randomly placed negative charges. Results from MC simulations and TK calculations are compared as solid line and solid circles, respectively. The protein dielectric boundary radius was Rd = 12Å in all calculations.

View larger version (32K): [in this window] [in a new window]   Fig. 5. (a) Total shift in ion binding constant pK, owing to solution ionic strength, is plotted as a function of salt concentration cs for the SC24 model. Other details are as for Fig. 4A. (b) The difference of excess chemical potentials ßµBex = ßµBex(cs) - ßµBex(cs = 1 mM) for a divalent cation in a bulk monovalent electrolyte solution is plotted as a function of salt concentration cs.

View larger version (11K): [in this window] [in a new window]   Fig. 6. Total shift in ion binding constant pK, owing to solution ionic strength, is plotted as a function of the protein net charge (Z). A reference state at salt concentration cs = 1 mM is chosen and the final salt concentration is 500 mM. Results from MC simulations and TK calculations are compared as solid and broken lines, respectively. Lines marked by solid circles correspond to a dielectric boundary radius Rd = 12Å; unmarked lines refer to a homogeneous dielectric response (Rd = 0Å).

Turning now to the consequences of a dielectric boundary, even for weakly charged proteins, enlarging the interior region of low permittivity has only a marginal effect on the salt shift of the binding constant as depicted in Figure 7 for the SC8 model. Although both binding sites remain in the highly polarizable region, a small reduction in pK is observed on expanding the dielectric radius R_d from zero to about 10Å. Further increasing R_d to incorporate the binding sites within the dielectric boundary has no effect at all on the calculated shift. When R_d extends beyond the protein radius (R_p = 14Å) a decrease of pK is again seen. The latter is a trivial effect caused primarily by salt exclusion, as mobile ions are very strongly repelled from the low dielectric region. For the highly charged model SC24, see Figure 8, MC data is almost constant in this scale and TK completely overpredicts the shifts, although it behaves in a qualitatively similar manner to TK results for model SC8.

View larger version (12K): [in this window] [in a new window]   Fig. 7. Total shift in ion binding constant pK, owing to solution ionic strength, is plotted as a function of the dielectric radius Rd for the SC8 model. A reference state at salt concentration cs = 1 mM is chosen and the protein (of radius Rp = 14Å) is allowed to bind two divalent cations. Results from MC simulations and TK calculations are compared as solid and broken lines, respectively.

View larger version (16K): [in this window] [in a new window]   Fig 8. Total shift in ion binding constant pK, owing to solution ionic strength, is plotted as a function of the dielectric radius Rd for the SC24 model. A reference state at salt concentration cs = 1 mM is chosen and the protein (of radius Rp = 14Å) is allowed to bind two divalent cations. The final salt concentration is 500 mM. Results from MC simulations and TK calculations are compared as open and solid symbols, respectively, joined by straight line segments to guide the eye.

To assess the impact on ligand binding equilibria of perturbations in protein charge and charge distribution (realized experimentally by site directed mutagenesis of ionizable residues), a closely related model is also considered. By suggesting the replacement of, say, a carboxylic acid group by the corresponding amine, model SC7 represents a "mutant" protein obtained in either of two ways. First, the SC7 variant of model SC8 is formed by simply increasing the valency of the central charge from -8 to -7. Second, the peripheral charge of model SP8 is neutralized again decreasing the overall charge to -7 units. Figure 9 shows the change in pK of divalent cation binding attributable to both these mutations as a function of dielectric radius. In both cases, the shift pK is very nearly independent of R_d as long as the binding sites are in the high dielectric region. Moreover, the calculated value of pK is also in reasonable agreement with experimental observations (Kesvatera et al. 1994) of calbindin, where a single charge mutation typically alters the calcium binding constant with 0.5–1.5 pK units at low salt and protein concentration (Svensson et al. 1990). Once the low dielectric cavity incorporates the binding sites and the point of mutation, however, pK climbs steeply and more so for the SP8 variant. When R_d matches the protein radius (R_d = 14Å), clearly unphysical binding constant shift around 25 pK units is apparent. Because the simulation results are faithfully tracked by the TK predictions over the full range of R_d, this failure must rest with inadequacies of the model rather than theoretical approximations.

View larger version (10K): [in this window] [in a new window]   Fig. 9. Total shift in ion binding constant pK, owing to a single charge mutation, is plotted as a function of the dielectric radius Rd for the SC8 and SP8 models. The salt concentration is cs = 1 mM and the protein concentration is cp = 0.02 mM. MC results are shown with solid (SC8) and broken (SP8) lines; TK data are denoted by solid circles.

	Conclusion

TOP Abstract Introduction Model systems Results and Discussion Conclusion Materials and methods References

 
The TK prescription based on the DH theory is an excellent approximation for studies of the binding of charged ligands to macromolecules. For sufficiently charged systems the linearization in the DH equation is unjustified and one must invoke a more accurate theory. Whether or not the nonlinear Poisson-Boltzmann theory is sufficient will depend on the system. For example, the description of highly charged aggregates, such as DNA in the presence of multivalent ions, where ion–ion correlations become important requires a correlated theory or numerical simulations (Guldbrand et al. 1986; Mel'nikov et al. 1999). The TK theory also becomes less reliable at high protein and low salt concentration, mainly because it approximates the excess chemical potential of the ligand with its value in the corresponding bulk salt solution. Typically this will happen when N_cc_p c_s.

The relative permittivity of the macromolecule is usually assumed to be much smaller than the surrounding aqueous solvent. For a highly charged macromolecule this is immaterial, when studying effects caused by salt concentration, as the leading term in the free energy difference is proportional to the net charge and the dielectric effect only enters in higher order terms. This can be demonstrated from the TK theory but is also found in MC simulations.

In a more general sense, any type of protein/solvent model that proposes an effective interaction potential by partitioning space into a finite number of dielectric continuna with different permittivities can be traced back to the intuitive design of Tanford and Kirkwood. Of course, the original TK analysis also simplifies the geometry of the ensuing electrostatic problem, but the essential idea of dielectrically distinct "background" regions remains common to more sophisticated models. The major weakness of applying such a scheme is the lack of an a priori prescription for locating the intervening dielectric boundaries and a method for estimating the dielectric response of the macromolecule. Charged ligand binding equilibria are highly sensitive to this aspect of electrostatic interactions, hence it may not always be possible to generate biologically relevant predictions without an independent method to sensibly determine these parameters. This is most clearly seen when studying the effect of charge mutations on the ligand binding, where models invoking a low dielectric response for the protein in many cases fail completely.

	Materials and methods

TOP Abstract Introduction Model systems Results and Discussion Conclusion Materials and methods References

 
Numerical simulations
Although the protein remains fixed at the center of the cell, the surrounding counterions and monovalent salt ions sample equilibrium configurations from the canonical (NVT) ensemble generated by the standard Metropolis MC algorithm (Metropolis et al. 1953). At high bulk protein concentrations, the cell radius is determined straightforwardly to obtain the appropriate volume per macromolecule, whereas for low protein and high salt concentrations the number of mobile particles required in the simulation rapidly becomes prohibitive. Tests have indicated, though, that the configuration-dependent part of the ligand chemical potential at the binding sites is rather insensitive to the cell size and that 50–100 salt pairs are typically sufficient to ensure convergence on the size-consistent limit. The excess chemical potential of the free calcium ion µ_F^ex is, however, greatly influenced by the macromolecular electrostatic field so that a naive reduction of the cell radius R_c would yield an erroneous result. Nevertheless, in low concentration, the protein will have only a marginal influence on µ_F^ex, particularly at high salt concentration where the electrostatic interactions are efficiently screened. Under these conditions it becomes more appropriate to calculate µ_F^ex from a bulk salt solution at the corresponding concentration. The number of counterions, N_c, was set to comply with overall electroneutrality, the temperature was fixed at T = 298.15K and the solvent dielectric constant was taken to be _s = 78.7.

The protein radius was chosen to be equal to R_p = 14Å whereas each mobile ion was assigned a common radius of R_k = R_s = 2Å for all k. The low dielectric response of the protein was explored in the extreme limit by introducing a spherical (radius R_d) vacuum cavity (_p = 1) centered on the cell origin. Commonly, the dielectric permittivity of a protein is taken to be comparable with that of a pure hydrocarbon phase, though there is no clear consensus on this issue (Warshel and Aqvist 1991; Antonsiewicz et al. 1994b; Simonson and Perahia 1995; Simonoson and Brooks 1996; Löffler et al. 1997; Papazyan and Washel 1997; Warshel and Papazyan 1998; Sham et al. 1998; Warwicker 1999; Baptista and Soares, unpubl.). To focus on qualitative physical mechanisms rather than quantitative prediction, we adopt the vacuum value here to exaggerate any effects arising from the dielectric discontinuity, although the consequences of high protein permittivity are examined in the following section. Note that the penetration of a real ion into a region of weak dielectric response is associated with a high energy, thereby implying a vanishingly low probability for the process. This is accounted for in the simulations by treating the dielectric discontinuity as an impenetrable boundary to the mobile ions. Normally we will set R_p > R_d so that the protein excluded volume overrides the dielectric depletion mechanism, though we admit the general condition R_s + max{R_p, R_d} defining the mobile ion exclusion sphere, see Equation 5. Values of R_d in the range of 0–18Å were studied, with the salt concentration c_s between 1 mM and 500 mM. For both equilibration and production runs, 10⁵ trial configurations per particle were generated.

Direct evaluation of the explicit TK equations (Tanford and Kirkwood 1957) was also carried out for the same model systems, with the limitation that all the protein charges and binding sites must lie inside the dielectric cavity. Moreover, these calculations were confirmed by numerical finite differece solutions obtained using the MEAD package (Bashford and Gerwert 1992).

Free energy calculations
When an ion binds to a protein, the changes in the free energy can be separated into two parts,

(8)

where the first term represents the electrostatic interactions and the second accounts for all the remaining contributions (structural changes, etc.). If G^rest is assumed to be independent of the salt concentration or charge mutations, then the logarithmic change in the stoichiometric binding constant is given by

(9)

where, for all the cases and models examined here, G^el_ref identifies a reference electrostatic free energy for the system at c_s = 1 mM with respect to a standard state and ß = 1/kT (k is the Boltzmann constant). In terms of the excess chemical potentials of protein (P), protein + ligand (P + L) and free ions (F), the electrostatic free energy can be expressed as

(10)

where µ_B^ex is the excess chemical potential of the bound ion(ligand). If several ligands are involved in the binding process, than µ_F^ex should, besides the excess chemical potential of each bound ligand, also include the interaction between the ligands. Similarly, µ_F^ex should contain the excess chemical potential of the free ligands.

Widom's (1963) test particle insertion method is an effective way to evaluate chemical potentials in continuum simulations. When applying the Widom technique to Coulomb fluids, where the explicit particle density is not too high, it should be realized that insertion of a single ion means that the total system becomes nonelectroneutral, possibly introducing spurious systematic errors. Most of this error can be corrected for by rescaling some or all of the charges in the system in such a way that the total system including the ghost particle becomes electroneutral. Such a scheme, referred to as the modified Widom technique, has been implemented and described in detail elsewhere (Svensson and Woodward 1988). The chemical potential of a free calcium ion is obtained by random particle insertions over the entire system, whereas µ_B^ex is calculated from an insertion directly into the binding site. This is done at every fifth particle move throughout the simulation to obtain a typical statistical error on the excess chemical potential around 0.01 - 0.03 kT.

TK analysis
Instead of directly computing the excess chemical potential at a binding site as an ensemble average, it is conveniently expressed in terms of the total electrostatic free energy change on detaching a ligand from the protein,

(11)

Consistent with the linear response analysis of Tanford and Kirkwood, the chemical potential for a free charged hard sphere ion with valency z and radius R_s is approximated using the DH theory of strong electrolytes,

(12)

where is the inverse DH screening length, which is proportional to the square root of the salt concentration. Following Tanford and Kirkwood, the electrostatic free energy of the neutral macromolecule is set to zero, and G^el for a protein containing N_p fixed charges becomes

(13)

with,

(14)

(15)

(16)

The solution ionic strength enters only through the variable x = (R_p + R_s) and the auxiliary functions,

(17)

described in an early paper by Kirkwood (1934). Notwithstanding the thorough discussion of results 13–16 in the original work by Tanford and Kirkwood (1957), it is instructive here to consider the three terms of Equation 13 in some detail. The direct Coulombic interaction between the charges of the protein is accounted for by the A_ijs and is independent of salt concentration, whereas the B_ijs are just the reaction field components (see Equation 2) contributing the effects from induced polarization charge at the dielectric boundary between protein and solvent. In the case of uniform dielectric permittivity, that is, _s = _p, all the B_ijs will vanish. The C_ij terms describe the effect of mobile counterions and added salt. Now suppose a ligand of valency z_l binds to the protein at radial position r_l subtending angles _kl with all other charges k, and consider the shift in its excess chemical potential on altering the salt concentration from x₁ to x₂, so that,

(18)

It is easy to see that only the C_ij terms contribute to µ^ex_B. Writing the net charge of the apo protein Z = ^Np_i=1 z_i and retaining only the leading terms of the multipole expansion (16) yield,

(19)

With the additional geometric simplification R_p + R_s = R_d             , then g(x) takes the form,

(20)

Further, noting that K₀ = 1 and K₂(x) = 1 + x + x²/3, the limiting cases of interest are obtained,

(21)

(22)

so that the excess chemical potential of the bound ion µ_B^ex finally becomes,

(23)

For a weakly charged protein and a simple ligand (both Z and z_l small), the monopole term may well be dominated by "dipolar" contributions describing the first moment of the charge distribution with m_i = z_ir_i/(R_p + R_s) and M_l = ^Np_i=1 m_i cos _il. The ionic strength dependent factor of the dipole term is also sensitive to the location of the dielectric boundary as well as the ratio of permittivities and will be approximately doubled on changing from a high to a low dielectric interior.

	Acknowledgments

 
This work has been supported by the Conselho Nacional de Desenvolvimento Científico e Technológico (CNPq/Brazil and FAPESP/Brazil), whom F.L.B.D.S. wishes to thank. We also thank Don Bashford for the access to MEAD. It is also a pleasure to acknowledge fruitful discussions with Henry S. Ashbaugh and Herman J.C. Berendsen.

The publication costs of this article were defrayed in part by payment of page charges. This article must therefore be hereby marked "advertisement" in accordance with 18 USC section 1734 solely to indicate this fact.

	References

TOP Abstract Introduction Model systems Results and Discussion Conclusion Materials and methods References

 
Antonsiewicz, J., McCammon, J.A., and Gilson, M.K. 1994a. The determinants of pK_as in proteins. Biochemistry 35: 7819–7833.

Antonsiewicz, J., McCammon, J.A., and Gilson, M.K. 1994b. Prediction of pH-dependence properties of proteins. J. Mol. Biol. 238: 415–436.[CrossRef][Medline]

Bashford, D. and Gerwert, K. 1992. Electrostatic calculations of the pK_a values of ionizable groups in bacteriorhodopsin. J. Mol. Biol. 224: 473–486.[CrossRef][Medline]

Bashford, D. and Karplus, M. 1990. pK_a's of ionizable groups in proteins: Atomic detail from a continuum electrostatic model. Biochemistry 29: 10219–10225.[CrossRef][Medline]

Bashford, D., Karplus, M., and Canters, G.W. 1988. Electrostatic effects of charge pertubations introduced by metal oxidation in proteins—a theoretical analysis. J. Mol. Biol. 203: 507–510.[CrossRef][Medline]

Beroza, P., Fredkin, D.R., Okamura, M.Y., and Feher, G. 1991. Protonation of interacting residues in a protein by a Monte Carlo method: Application to lysozyme and the photosynthetic reaction center of rhodobacter sphaeroides. Proc. Natl. Acad. Sci. 88: 5804–5808.[Abstract/Free Full Text]

Böttcher, C.J.F. 1973. Theory of Electric Polarization. Elsevier, Amsterdam, The Netherlands.

Fushiki, M., Svensson, B., Jönsson, B., and Woodward, C.E. 1991. Electrostatic interactions in protein solution—A comparison between Poisson-Boltzmann and Monte Carlo calculations. Biopolymers 31: 1149–1158.[CrossRef][Medline]

Guldbrand, L., Nilsson, L., and Nordenskiöld, L. 1986. A Monte Carlo simulation study of electrostatic forces between hexagonally packed DNA double helices. J. Chem. Phys. 85: 6686–6698.[CrossRef]

Harvey, S.C. 1989. Treatment of electrostatic effects in marcomolecular modeling. Proteins: Struct., Func. Genet. 5: 78–92.

Hill, T.L. 1955. Approximate calculations of the electrostatic free energy of nucleic acids and other cylindrical macromolecules. Arch. Biochem. Biophys. 57: 229–239.[CrossRef]

Hill, T.L. 1956. Influence of electrolyte on effective dielectric constants in enzymes, proteins and other molecules. J. Chem. Phys. 60: 253–255.[CrossRef]

Honig, B. and Nicholls, A. 1995. Classical electrostatics in biology and chemistry. Science 268: 1144–1149.[Abstract/Free Full Text]

Kesvatera, T., Jönsson, B., Thulin, E., and Linse, S. 1994. Binding of Ca² to Calbindin D_9K: Structural stability and function at high salt concentration. Biochemistry 33: 14170–14176.[CrossRef][Medline]

Kesvatera, T., Jönsson, B., Thulin, E., and Linse, S. 1999. Ionization behaviour of acidic residues in calbindin D_9k. Proteins: Struc., Func. Genet. 37: 106–115.

Kesvatera, T., Jönsson, B., Thulin, E., and Linse, S. 2001. Focusing the electrostatic potential at EF-hands of Calbindin D_9K. Titration of acidic residues.

King, G., Lee, F.S., and Warshel, A. 1991. Microscopic simulations of macroscopic dielectric constants of solvated proteins. J. Chem. Phys. 95: 4366–4377.[CrossRef]

Kirkwood, J.G. 1934. Theory of solutions of molecules containing widely separated charges with special applications to zwitterions. J. Chem. Phys. 2: 351–361.[CrossRef]

Levesque, D., Weis, J.J., and Hansen, J.P. 1986. Simulation of classical fluids. Monte Carlo Methods Stat. Phys. 5: 47–119.

Linse, S., Brodin, P., Johansson, C., Thulin, E., Grundström, T., and Forsen, S. 1988. The role of protein surface changes in ion binding. Nature 335: 651–652.[CrossRef][Medline]

Linse, S., Johansson, C., Brodin, P., Grundström, T., Drakenberg, T., and Forsen, S. 1991. Electrostatic contribution to the binding of calcium in calbindin D_9K. Biochemistry 30: 154–162.[CrossRef][Medline]

Linse, S., Jönsson, B., and Chazin, W.J. 1995. The effect of protein concentration on ion binding. Proc. Natl. Acad. Sci. 92: 4748–4752.[Abstract/Free Full Text]

Löffler, G., Screiber, H., and Steinhauser, O. 1997. Calculation of the dielectric properties of a protein and its solvent: Theory and a case study. J. Mol. Biol. 270: 520–534.[CrossRef][Medline]

Mel'nikov, S., Lindman, B., Khan, M.O., and Jönsson, B. 1999. Phase behaviour of a single DNA in mixed solvents. J. Am. Chem. Soc. 121: 1130–1136.[CrossRef]

Metropolis, N.A., Rosenbluth, A.W., Rosenbluth, M.N., Teller, A., and Teller, E. 1953. Equation of state calculations by fast computing machines. J. Chem. Phys. 21: 1087–1097.[CrossRef]

Papazyan, A. and Warshel, A. 1997. Continuum and dipole-lattice models of solvation. J. Phys. Chem. B 101: 11254–11264.[CrossRef]

Penfold, R., Warwicker, J., and Jönsson, B. 1998. Electrostatic models for calcium binding proteins. J. Phys. Chem. B 108: 8599–8610.[CrossRef]

Sham, Y.Y., Chu, Z.T., and Warshel, A. 1997. Consistent calculations of pK_a's of ionizable residues in proteins: Semi-microscopic and microscopic approaches. J. Phys. Chem. B 101: 4458–4472.[CrossRef]

Sham, Y.Y., Muegge, I., and Warshel, A. 1998. The effect of protein relaxation on charge-charge interactions and dielectric constants of proteins. Biophys. J. 74: 1744–1753.[Abstract/Free Full Text]

Shire, S.J., Hanania, G.I.H., and Gurd, F.R.N. 1974. Electrostatic effects in myoglobin. Hydrogen ion equilibria in sperm whale Ferrimyoglobin. Biochemistry 13: 2967–2974.[CrossRef][Medline]

Simonson, T. and Brooks III, C.L. 1996. Charge screening and the dielectric constant of proteins: Insights from molecular dynamics. J. Am. Chem. Soc. 118: 8452–8458.[CrossRef]

Simonson, T. and Perahia, D. 1995. Microscopic dielectric properties of cytochrome c from molecular dynamics simulations in aqueous solution. J. Am. Chem. Soc. 117: 7987–8000.[CrossRef]

Spassov, V. and Bashford, D. 1998. Electrostatic coupling to pH-titrating sites as a source of cooperativity in protein-ligand binding. Prot. Sci. 7: 2012–2025.[Abstract]

Svensson, B.R. and Woodward, C.E. 1988. Widom's method for uniform and non-uniform electrolyte solutions. Mol. Phys. 64: 247–259.[CrossRef]

Svensson, B., Jönsson, B., and Woodward, C.E. 1990. Electrostatic contributions of the binding of Ca²⁺ in calbindin mutants. A Monte Carlo study. Biophys. Chem. 38: 179–183.

Svensson, B., Jönsson, B., Woodward, C.E., and Linse, S. 1991. Ion binding properties of calbindin D_9k—A Monte Carlo simulation study. Biochemistry 30: 5209–5217.[CrossRef][Medline]

Tanford, C. and Kirkwood, J.G. 1957. Theory of protein titration curves. I. General equations for impenetrable spheres. J. Am. Chem. Soc. 79: 5333–5347.[CrossRef]

Tanford, C. and Roxby, R. 1972. Interpretation of protein titration curves. Application to lysozyme. Biochemistry 11: 2192–2198.[CrossRef][Medline]

Warshel, A. and Aqvist, J. 1991. Electrostatic energy and macromolecular function. Annu. Rev. Biophys. Biophys. Chem. 20: 267–298.[CrossRef][Medline]

Warshel, A. and Papazyan, A. 1998. Electrostatic effects in macromolecules: Fundamental concepts and practical modeling. Curr. Opin. Struct. Biol. 8: 211–217.[CrossRef][Medline]

Warshel, A., Russel, S.T., and Churg, A.K. 1994. Macroscopic models for studies of electrostatic interactions in proteins: Limitations and applicability. Proc. Natl. Acad. Sci. 81: 4785–4789.

Warwicker, J. 1999. Simplified methods for pK_a and acid pH-dependent stability estimation in proteins: Removing dielectric and counterion boundaries. Prot. Sci. 8: 418–425.[Abstract]

Warwicker, J. and Watson, H.C. 1982. Calculation of the electric potential in the active site cleft due to -helix dipoles. J. Mol. Biol. 157: 671–679.[CrossRef][Medline]

Widom, B. 1963. Some topics in the theory of fluids. J. Chem. Phys. 39: 2808–2812.[CrossRef]

CiteULike    Connotea    Del.icio.us    Digg    Reddit    Technorati    What's this?