J. Phys. Chem. 1995, 99, 1614-1618
1614
Integral Equation Theory for Proteins: Application to Ca2+Binding in Calbindin Dgk Bo Svensson and C. E. Woodward* Department of Chemistry, University College Australian Defence Force Academy, University of New South Wales, ACT 2600, Australia Received: July 14, 1994; In Final Form: October 26, 1994@
The site-site Omstein-Zernike equation with the RISM closure is used to study the binding of calcium ions to the protein calbindin. A static, charged hard-sphere model is used for the protein, and the primitive model is used for the surrounding electrolyte solution. The free energy of calcium binding to the protein over a range of ambient electrolyte concentrations and mutations was determined from the random phase approximation. Theory is found to be in remarkably good agreement with computer simulation predictions. Extensions of the theory to include a molecular description of the solvent are discussed.
1. Introduction
The importance of electrostatic interactions in biological systems has now been well established thanks to a considerable number of experimental and theoretical studies focused on highly charged biomolecules.' One important class of such biomolecules are those proteins responsible for the binding of cationic species. This group includes a large number of calcium binding proteins which have been the subject of intensive experimental investigation.2 Quite recently, Monte Carlo simulations of the protein calbindin D9k demonstrated the importance of electrostatic interactions in determining the binding affinity of calcium ions (Ca2+) for a range of electrolyte concentrations and valencies and protein m~tations.~ Calbindin is a charged protein with a valency of -8 and is capable of binding two (Ca2+) very strongly. In the simulations the protein atoms were modeled as (partially) charged hard spheres, of appropriate radii, with positions consistent with the X-ray determined structure. The primitive model was used for the surrounding electrolyte solution; i.e., the salt ions were treated as charged hard spheres immersed in a dielectric continuum, with a dielectric constant appropriate for the solvent. The calculated binding constant shifts were in remarkably good agreement with experiments, despite the model lacking a dielectric discontinuity at the protein surface. Calculations were also performed using both the linearized (Debye-Huckel) and nonlinear Poisson-Boltzmann theories. In both these approximations, the asymmetric electrostatic potential is solved on a distributed grid of points. It was established that the Debye-Huckel (DH) and PoissonBoltzmann predictions for the Ca2+ binding constants agreed very well with simulations in the 1:l electrolytes. The fact that the DH theory performed so well was quite surprising, given that it is a linear approximation and that the surface charge density of calbindin is moderately high. Indeed, the separate ion-protein density profiles predicted by the DH are qualitatively poor. For example, the co-ion density close to the protein surface is predicted to be negative! However, the free energy for Ca2+ binding is an integral over the charge density. The long-ranged nature of the Coulomb interactions and the fact that the net charge density profile from the DH is reasonably accurate is a plausible explanation for its success. The DH theory was also quite accurate in 1:-2 electrolytes where the divalent ion was the co-ion to the protein. Such was not the @
Abstract published in Advance ACS Abstracts, December 1, 1994.
0022-365419512099- 1614$09.00/0
case for 2:- 1 electrolytes; however, the Poisson-Boltzmann approximation continued to perform well in this system. In this work will explore a different theoretical approach to study the binding of Ca2+ to calbindin. In particular, we shall apply the site-site Omstein-Zemike (SSOZ) equation. The SSOZ equation was developed by Chandler and co-workers to study molecular fluids and is an extension of a similar integral equation for monat~mics.~ As it deals specifically with sitesite rather than angular dependent, orientational, correlation functions, it provides a particularly convenient formalism for the calculation of molecular fluid properties. The so-called reference interaction site model RISM and extended RISMS approximations provide closures of the SSOZ equation based on similar forms from the theory of monatomics. Both of these approximations have been extensively applied to a range of problems. These include infinite dilution studies of molecular conformational equilibria in solution6 and electron solvation.' More recently, the SSOZ equation and the RISM approximation have been generalized to nonuniform molecular fluids and polymers.* Despite the well-known shortcoming of these closures in predicting long-ranged orientational correlations, they have still proved to be very useful tools in providing, sometimes quantitative, predictions for the structure and thermodynamics of molecular systemsSg Indeed, the RISM approximations are admirably suited for the study of large biomolecules in solution. This notwithstanding, such integral equation theories are still not widely used theoretical tools in this area. Part of the reason for this is undoubtedly the computational expense; some calculations would, for example, require the simultaneous solution of many hundred correlation functions. Improvements in computational hardware, e.g., the advent of fast dedicated machines, however, have made such calculations more viable. This has been illustrated in recent work by Kitao et al., who solved the extended RISM for the protein melittin immersed in water.'O The latter was modeled using the SPC potential. Those authors solved for the solvent structure around 436 protein atoms! The calculations presented here on calbindin will be on a simplified protein model with no more than 49 explicit sites. For the study of charged biomolecules it would seem that the extended RISM (eRISM) provides the most useful approximation. It is based on the hypemetted chain (HNC) approximation which has been shown to be quite successful for monatomic charged fluids. However, it is quite difficult to obtain properly converged solutions of the eRISM approximation. For example, Kitao et al. found in their iterative solutions 0 1995 American Chemical Society
Ca2+ Binding in Calbindin Dgk that convergence depended strongly on the initial guess. This clearly makes it difficult to determine whether a stable solution even exists in some cases. In this work we will use a RISM closure based on the mean spherical approximation (MSA). Using the nomenclature of Monson and Morris? we shall denote this the RISM-MSA. This approximation can be considered a type of extended Debye-Hiickel approach, which accounts for hard-core correlations. Given the success of the nonuniform DH theory in predicting Ca2+binding free energies in calbindin, we expect the RISM-MSA to be a potentially useful approximation, especially in 1:-1 and 1:-2 electrolytes. Furthermore, the RISM-MSA is very easy to solve and usually quite robust. Though we chose to solve the RISM by iterating on both the direct and total correlation functions at discrete quadrature, one can also express the RISM-MSA as a result of a variational problem, where one seeks the value of the direct correlation functions only inside the core. The functional to be minimized is the random phase approximation to the free energy,” which provides one with a natural expression for the free energy of binding.
2. Theory
2.1. A Simple Model for Calbindin in Solution. In a detailed calbindin model, used in earlier simulation work, all 722 atoms in the protein were explicitly treated being placed at fKed sites according to coordinates PDB3ICB of the Brookhaven Protein Data Bank. We believed this to be a reasonable model as ‘HNMR studies indicated that the solution structure of calbindin was very similar to its crystalline form. The crystalline structure of calbindin was originally determined by Szebenyi and Moffat” and has more recently been refined by Svensson et a1.I2 In the detailed model, all protein atoms were treated as uncharged or partially charged hard spheres. It was found that a suitably accurate representation of the charge density for calbindin could be obtained by assigning a partial charge of -0.5 to all carboxylic oxygens and a charge of +1 to all lysine nitrogens. In order to facilitate the numerical solutions in this work, we decided to simplify the model somewhat and reduce the number of protein atoms to be explicitly treated. In particular, the hydrophobic interior of the protein was modeled as a large hard sphere with a radius of 11.6 A. In addition to this core, we only treated charged atoms explicitly (i.e., carboxylic oxygens and lysine nitrogens); these were given the same charges as in the detailed model described above, but their radii were increased from about 2 to 6 hi in order to give the simplified model approximately the same excluded volume as the detailed one. That is, the expanded radius accounted for other neutral surface atoms not explicitly accounted for. The total number of explicit sites in this new model for “calbindin” was thus reduced to 49. Despite this simplification, it was found that the simulated Ca2+ binding free energies in the simplified and detailed model differed by at most 20%. The electrolyte solution into which the calbindin is immersed is approximated by the primitive model, wherein the salt ions were treated as charged hard spheres all with radius 2 hi. All particle interactions thus had the form
where qi is the (partial) charge of the ith particle which has diameter oi. The quantity r is the distance between particles. The water solvent is not accounted for explicitly and enters the calculation only via the relative dielectric constant, e. We chose
J. Phys. Chem., Vol. 99, No. 5, 1995 1615 = 78.7, which is the appropriate value for water at a temperature of 298 K. 2.2. Monte Carlo Simulations. Simulations were performed on the simplified calbindin model, described above, using the cell model. In this model, all the salt ions were confined to a spherical cell, with the rigid protein molecule placed at the center. In addition to the handful of counterions to the protein, the cell contained 200 salt ions. The cell radius was chosen to give the desired salt concentration. At the lowest salt concentration of 2 mM, the resulting protein concentration was similar to typical experimental values (about 20 pM). At higher salt concentrations the system was large enough to make finitevolume effects negligible. The free energy for Ca2+ binding is the difference between the excess chemical potential of the Ca2+ ions at the binding sites in the protein and in the bulk electrolyte solution. The bulk solution was simulated in a cube with periodic boundary conditions. The chemical potentials were determined by using a variation of the Widom test particle developed specifically for charged systems.13 A description of the simulation details can be found in previous work3 and will not be reiterated here. All simulations were carried out on an IBM RS6000 320, with a typical run averaging about half an hour of CPU time. 2.3. The Site-Site Omstein-Zernike Equations. The SSOZ equation has been extensively reviewed in the literature so it will suffice to present here a brief development relevant to our particular problem. As our model consists of a single protein molecule immersed in a primitive model electrolyte, the appropriate infinite dilution limit of the SSOZ equation is used to obtain the following relations for the protein-ion correlation functions in Fourier space, E
In this expression, the subscripts i a n d j denote the protein sites and the Greek subscripts denote the salt species. We have (3) yhere ea is the number density of the a t h salt species and hia(k) is the Fourier transform of the total correlation function, hia(r), for site i and species a , i.e.,
hia(k) = ( 2 / k ) l d rhia(r)rsin(2nkr)
(4)
Similarly,
where E&) is the Fourier transform of the direct correlation function cia(r). Finally, &j(k) is the Fourier transform of the intramolecular correlation function between sites i and j in the protein, d i j ( k )=
sin( 2nk1,) 2nk1,
(6)
with 10 the distance between the sites. The correlation functions for the salt species satisfy the usual Omstein-Zernike equations
fiM(k) = e M ( k ) +
eaY(k) fiyp(k>
(7)
Y
As Coulombic interactions are present, it is necessary to renormalize eq 2.14 We let the long-ranged part of the direct
1616 J. Phys. Chem., Vol. 99, No. 5, 1995
Svensson and Woodward 2.4. Chemical Potential. The random phase approximation (RPA) functional is given bf
correlation functions be denoted by
and
’I
ia
(9) with p = l / k ~ Twhere T is the temperature and k~ is Boltzmann’s constant. Thus, for example
tia(k)
-
&)ia(k)
for k
-o
(10)
Using the long-ranged contribution to the total site-solute correlation function
and
J
Taking the functional derivative of eq 19 with respect to Ci,(r) gives dZ/dCia(r) = -ea1/2gia(r) that is for r inside the core dZ/dCia(r) = 0
(20) (21)
Using eq 18 to approximate Cia(r) outside the core, eq 21 gives the well-known result that the RISM-MSA direct correlation function minimizes the RPA functional. Introducing a charging parameter A which multiples the long-ranged (electrostatic) part of the protein-ion potential, one obtains in the RISM-MSA approximation.
= Jdr D u k a ( r ) gka(A,r) @a
(22)
a,k
With similar definitions for the solute correlation functions, one is able to rewrite eqs 2 and 7 in renormalized form,
and hence
PAP,, = Z(A = 1) - Z(A = 0)
(23)
where Apel is the change in the chemical potential of the protein obtained by switching on the protein charges.
3. Results and
where
with a similar definition for the solute correlation functions. In order to solve these equations further, relations must be found between the direct and total correlation functions. For a bulk electrolyte it has been established that the hypemetted chain (HNC) appr~ximation’~ is reasonably accurate; this is given by
+
+
z$
