J . Phys. Chem. 1988, 92, 7134-7141
7134
occurring at the electrode/electrolyte interface. X-ray absorption spectroscopy has been shown to be a viable tool for the in situ determination of the concentrations of high-Z (>20) elements at electrode/electrolyte interfaces.
Acknowledgment. Our work was generously supported by the Materials Science Center at Cornel1 University, the National
Science Foundation, and the Office of Naval Research. H.D.A. also acknowledges support by the Presidential Young Investigator Award Program at the NSF and the A.P. Sloan Foundation. The authors acknowledge the assistance of G. M. Bommarito and M. J. Albarelli in collecting the data. Registry No. I-, 20461-54-5; I*, 7553-56-2; Pt, 7440-06-4.
Ionic Strength Dependence of Enzyme-Substrate Interactions. Monte Carlo and Poisson-Bottzmann Results for Superoxide Dismutase Russell J. Bacquet? J. Andrew McCammon,* Department of Chemistry, University of Houston, Houston, Texas 77004
and Stuart A. Allison Department of Chemistry, Georgia State University, Atlanta, Georgia 30303 (Received: March 23, 1988; In Final Form: June 1, 1988)
Monte Carlo (MC) simulations have been carried out for simplified models of the enzyme superoxide dismutase at infinite dilution in several aqueous salt solutions. Comparison is made to results obtained from numerical algorithms based on Poisson-Boltzmann (PB) theory. The impact of approximations inherent in the PB approach is found to be small. The diffusion-limited rate constant for enzymesubstrate association is computed from the MC data. The results are qualitatively in accord with experimental data. Better agreement can be obtained by using more detailed models, but this requires the use of PB rather than MC due to computational requirements. The major conclusion of this study is therefore that the approximations inherent in the PB theory are acceptable in the present application.
1. Introduction
Experimental'*2and t h e o r e t i ~ a l ~works - ' ~ have demonstrated that the charge distribution in an enzyme is important not only in the final orientation and docking of substrate but also during its diffusive approach to the active site from an initially large separation. Accurate calculation of the electrostatic interaction energy is difficult, even for simple model systems, since the nonuniform distribution of salt ions contributes and the dielectric responses of both solvent and protein must be considered. Given the potential energy surface and assuming diffusion control of the reaction, a rate constant k can be computed.6 The ionic strength (I)dependence of k is of particular interest. Although the addition of salt inevitably reduces the interaction between two charged species, it is not always immediately clear how the reaction rate will be affected, since some moments of their charge distributions may facilitate reaction while others impede it. For example, their net charges may have the same sign while a dipole or quadrupole moment is favorably oriented with respect to an active site. This is the case for the reaction of superoxide dismutase (SOD) with superoxide radical, 02-.With net charges of -4 and -1, respectively, there exists an overall repulsion, but the potential energy near the two active sites of the enzyme is lowered by its large and strategically placed quadrupole moment.' Experimentally, the reaction rate declines as salt is added's2 (see Figure 9), implying that reduction of the steering quadrupole moment is of more consequence than reduction of the net charge interaction. However, it has also been suggested' that direct interference by salt in the catalytic mechanism may explain some of the rate dependence. Theoretical study of a simple model for this system, using an approximate treatment of salt effects, yielded a rate constant which initially increased but then declined as the salt addition continued.'
Recently, more realistic enzyme and dielectric models have been e ~ a m i n e d ~ J * with * - ' ~ salt effects determined by a numerical form of the linearized Poisson-Boltzmann equation. The increased realism of these studies constitutes a major advance. Nevertheless, several approximations are made in this approach and it is of interest to determine their impact. This can be done only for simple models, where alternate and more accurate methods can be applied for comparison. Here we reexamine the ionic strength dependence of the SOD/OF interaction energy surface and the rate constant computed from it. A simple model of these reactants in various solutions of aqueous NaCl is studied. Enzyme shape, charge distribution, and dielectric response are treated in a crude fashion, but, within that context, salt effects are accurately obtained by Monte Carlo simulation. For comparison, and to explore the approximations involved, Poisson-Boltzmann and Debye-Huckel computations for this model are performed. The ionic strength dependence of the rate constant is then determined by means of a Brownian dynamics trajectory method. The model is detailed in section 2 of this paper, and the application of the various (1) Cudd, A.; Fridovich, I. J. Biol. Chem. 1982, 257, 11443. (2) Argese, E.; Viglino, P.; Rotilio, G.; Scarpa, M.; Rigo, A. Biochemistry 1987, 26, 3224. (3) McCammon, J. A.; Harvey, S. C. Dynamics of Proteins and Nucleic Acids; Cambridge University: Cambridge, 1987. (4) McCarnmon, J. A. Science ( Wushington, D.C.) 1987, 238,486. ( 5 ) Northrup, S . H.; Luton, J. A.; Boles, J. 0.; Reynolds, J. C. L. J. Cornput.-Aided Mol. Des. 1988, I , 291. (6) Northrup, S . H.; Allison. S. A.; McCammon, J. A. J. Chem. Phys. 1984, 80, 1517. (7) Allison, S. A.; McCammon, J. A. J. Phys. Chem. 1985, 89, 1072. (8) Warwicker, J.; Watson, H. C. J . Mol. Biol. 19112, 157, 671. (9) Sharp, K.; Fine, R.; Honig, B. H. Science (Washington, D.C.) 1987, 236, 1460. (10) Allison, S. A.; Baquet, R. J.; McCammon, J. A. Biopolymers 1988, 27. 251. (1 1) Northrup, S. H.; Boles,J. 0.; Reynolds, J. C. L. J . Phys. Chem. 1987, 91, 5991. (12) Klapper, I.; Hagstrom, R.; Fine, R.; Sharp, K.; Honig, B. Proteinr 1986, 1, 41. I~
'Current address: Agouron Pharmaceuticals, Inc., 11025 North Torrey Pines Road, P.O. Box 12209, La Jolla, CA 92037.
~
0022-3654/88/2092-7134$01.50/00 1988 American Chemical Society
Ionic Strength Dependence of Superoxide Dismutase
The Journal of Physical Chemistry, Vol. 92, No. 25, 1988 7135
techniques is described in section 3. Results are presented and discussed in section 4. Section 5 contains a summary and concluding remarks.
2. Model This work deals with the calculation by several different methods of potential energies of mean force for salt ions and substrate particles near an enzyme. In this section we describe the model system, beginning with the solvent, dielectric, and enzyme. These features are the same regardless of which method is being applied. Next we discuss the different ways in which salt is handled. Finally, the additional details of the model needed for the rate calculations are presented. Explicit inclusion of water molecules in the simulation of a salt solution is generally not feasible due to the large number required (e.g., -250000 for a 0.009 M solution modeled with 100 ions). While the dielectric properties of solvent are important, its molecularity has little significance for our purposes. Thus, it is adequate to treat water as a dielectric continuum with c = 78.54. For simplicity, all electrostatic interactions are reduced by this factor although the dielectric response of the protein interior is much less than that of solvent. The error can be severe for interior protein charges but is less serious for surface protein charges interacting with free ions or charged ligands in solution. The role of clefts in focusing the electrostatic potential* is also missed when t is uniform throughout the system. The use of a low internal dielectric constant would provide results more appropriate for comparison with experiment. However, the added computational expense would not enhance the comparison of theoretical methods which is the primary concern of this work. To the extent possible for a very simple model, we seek to reproduce the shape, size, and charge distribution of the enzyme superoxide dismutase (SOD). Two views of SOD based on its detailed structure obtained from X-ray cry~tallographyl~ appear in Figure 1. Contours of the potential energy surface for the dielectric model just described are shown based on 2196 assigned partial charges in the protein and no salt in the solution. Although roughly spherical, SOD is somewhat elongated in one of its dimensions. There are numerous surface irregularities including channels leading to the catalytic metal ions. The electrostatic force on the substrate is attractive near the active sites and generally repulsive elsewhere. The model protein, also shown in Figure 1, is a smooth sphere of radius 30.0 8, containing three collinear point charges: one of -23.38 proton units at the center and two of 9.69 units located 10 8,above and below. This charge set approximately reproduces the monopole, dipole, and quadrupole terms associated with all the charged groups in the detailed X-ray structure. The development of a similar model has been thoroughly described.I4 The energy contours of Figure 1 clearly show that this model is only an extremely idealized version of the real system. Nevertheless, we shall refer to it as 3SOD and will compare theoretical results based on it with the experimental data for the enzyme. For the 3SOD model we define a polar axis which passes through the three charges. The equatorial region then corresponds to the plane with a z coordinate of zero. Calculations were also done for an enzyme model with the same size and shape, but with a single charge of -4 located at its center. This model, termed lSOD, represents a hypothetical SOD-like molecule with its charges uniformly distributed, so that no steering force is present. In the Monte Carlo simulation Na' and C1- are modeled by charged particles with sizes defined by a repulsive exponential core potential of the form
\ \
- - - - - _ _ _ -.' ..-. . --____--
/'
.
0
0
I
I
, -
' . 0
0 .
- e - - -
/
0
/
- 1
'
0
/
\
0 * - - - -
0
I
1 -
- e - - _
...... \' \
.'- - - - _ _ - -.-
-
\
- - - e - -
0
-----. .
0
'
\
\
.- _ _ _ - - -.'' ..'.;--- - _ _ - - *
'\
/'
0
I
C
-
\
/
/
I
Figure 1. Electrostatic potential energy contours for two views of detailed SOD (a and b) and for the 3SOD model (c). Solid lines are 0.25kT (-0.15 kcal/mol) apart and are given for values between 2.0- and -1.OkT. The outermost solid contour is at 0.25kT; the energy increases as the two large lobes are traversed and decreases as the two small lobes are entered. To show the saddle points more clearly, dashed contours are given at 0.30-, 0 . 3 5 , 0.40-, and 0.45kT.
TABLE I: Core Potentials A , erg
particle Na' c1protein
Na+ 3.75 x 1 0 4 3 3.0 x 10-13 3.0 X lo-"
c13.0 x 1 0 4 3 2.25 x 10-13 3.0 x 10-13
r*, A 0.95 1.81
27.0
j and Aij,ri*, and rj* are parameters given in Table I. The ion-ion
where rij is the center-to-center distance separating species i and (13) Tainer, J. A.; Getzoff, E. D.; Beem, K. M.; Richardson, J. S.; Richardson, D. C. J . Mol. Biol. 1982, 160, 181. (14) Allison, S. A,; Ganti, G.; McCammon, J. A. Biopolymers 1985, 24, 1323.
interactions are based on studies of simple electrolyte solutionsIs while the enzyme-ion interactions are reasonable but rather arbitrary choices. The equilibrium distribution of ions near the (1 5 ) Friedman, H. L.; Smitherman, A,; DeSantis, R. J. J . Solution Chem. 1973, 2, 59.
Bacquet et al.
7136 The Journal of Physical Chemistry, Vol. 92, No. 25, 1988 TABLE 11: Simulation Parameters' ionic strength radius of simulation or salt concn,b M region, A 0.0005
0.009 0.104 0.507
340 130 70 60
max. length no. of Na',
C1- ions 52, 52, 84, 242,
48 48 80 238
no. of
of trial
move, A
no. of trial movesc
configurations analyzed
124 42 17 12
18 X lo6 14 X lo6 11.808 X lo6 8.64 X lo6
180 000 140 000 72 000 18000
'Parameters given here were used in the simulationsof both lSOD and 3SOD. bThis is the local ionic strength or salt concentrationat a point far from the enzyme but not immediately adjacent to the boundary of the simulation cell. CApproximately70% of trial moves were accepted. enzyme is determined primarily by Coulombic forces with the precise nature of the core interactions having a minimal impact. In the Poisson-Boltzmann and Debye-Huckel calculations, salt ions are not included as discrete particles. Instead, the net density of ionic charge at each point is computed from the total electrostatic potential at that point. Since the influence of ionic charge at all other points must be included and is not initially known, it is clear that the salt distribution must be determined by an iterative procedure. Superoxide radical has the same charge and approximately the same size as CI-. Due to this similarity, it is unnecessary to include one or more particles representing substrate in the Monte Carlo simulations or the Poisson-Boltzmann calculations. To a very good approximation, the potential of mean force for 02-is equivalent to that for C1-. The model system adopted for the BD calculations is somewhat different. The equilibrium ion distribution, dielectric properties of water, and enzyme charges are represented by an effective potential energy surface obtained from the Monte Carlo work. The enzyme shape and size are unchanged, but reactive patches are defined as those portions of the surface within 10' of the polar axis. 0, is a sphere with a hydrodynamic radius of 2.05 A, based on ab initio studies16 of O2-(H2O),, and a diffusion constant defined by stick boundary conditions. These values are improvements over those used in earlier s t ~ d i e s .They ~ result in a factor of 2 change in the mutual diffusion constant, DM,which leads to a roughly comparable change in the computed rate constants.
3. Methods In this work an approximate potential of mean force between a simplified enzyme and a particle is computed as a function of salt concentration. The particle represents either the substrate or a salt anion, reflecting their real similarity for the enzyme/ O2-/NaC1 system under consideration. Due to the simplicity and symmetry of the model system, computer simulation of the salt distribution is feasible; Le., a moderate amount of computer time will produce an acceptably low level of statistical noise. N o other approximations are made in the course of obtaining a solution by this method. Although the real strength of the procedure is its applicability to more realistic and complex systems, the finite difference Poisson-Boltzmann algorithm is also applied to the same simple protein model. Errors due to the inherent approximations of this formulation can in principle be detected by comparison with an accurate simulation of the same model. Finally, a Brownian dynamics trajectory method3" is used to compute diffusion-controlled reaction rates from the potentials of mean force. A brief description of each of these theoretical approaches and details of their implementation are provided in this section. 3.1. Monte Carlo Simulations. Equilibrium ion distributions were determined by using Metropolis Monte Carlo computer sim~lations'~ in the canonical ensemble at a temperature of 298.15 K. No periodic boundary conditions (PBC) were imposed. Salt and protein were confined to a large spherical region. The long range of the Coulombic interaction coupled with the inhomogeneous ion distribution makes PBC with a minimum image or (16) Lopez, J. P.; Albright, T. A,; McCammon, J. A. Chem. Phys. Lett. 1986, 125, 454.
spherical cutoff inappropriate. Ewald ~ u m m a t i o n is ' ~ feasible, but the chosen method is fully adequate when the simulation region is sufficiently large and the ion distribution near the hard outer boundary is disregarded. Parameters describing the simulations are listed in Table 11. Data were collected on a spherical coordinate grid with averaging over $ (the angle around the polar axis) and with +z equivalent to -z, since deviations from such symmetry can only be due to noise. For the 3SOD simulations, the ranges of r (30 A to cell boundary) and 8 (0 to a/2 rad) were covered by -55 and 19 unevenly spaced points, respectively. A relatively coarse grid was used both far from the protein, where the salt concentration gradient is small, and near the polar axis, where the $-averaged region of space is small. The ionic strength given in Table I1 for each simulation is obtained from the local salt concentrations 10-15 A from the boundary of the simulation cell where interactions with the enzyme have been screened to zero and no spurious boundary effects are present. It therefore corresponds to the bulk ionic strength of a solution that contains enzyme at infinite dilution and is not equivalent to the average ionic strength within the simulation cell. 3.2. Nonlinear Least-Squares Fits. It is desirable to have a functional representation of the Monte Carlo results. Using standard IMSL routines, a least-squares fit of the data to the function
was obtained for each simulation. q, e, c , and /3 = (kT)-' have their conventional meanings with c = 78.54, T = 298.15 K,and q = -1. N is the number of model protein charges while qp and rp are their magnitudes and positions as given in section 2. cp and up are the variables optimized in the fitting procedure. In effect, the protein charges are scaled and screened so that the potential field due only to them closely matches that calculated for the entire enzyme/salt system. Equation 2 is similar to eq 3, discussed in the next section, which contains the Debye K and an ion exclusion radius, but its use here is justified only because it does provide a fully satisfactory smooth curve through the M C points (see Figures 4-6). Optimum values for cp and up are listed in Table 111. One measure of the quality of the global fit is the sum of squared residuals given in Table 111. However, this quantity also reflects and is likely dominated by the amount of noise in the data points. 3.3. Poisson-Boltzmann and DebyeHiickel Calculations. Poisson-Boltzmann (PB) and Debye-Huckel (DH) theory have often been used to estimate the distribution of salt near a charged biomolecule. PB neglects correlations in position among the small ions, and D H makes the further approximation of linearizing an exponential factor (so that the DH method is also referred to as linearized PB). We have numerically computed D H and PB results for both the lSOD and 3SOD protein models using a finite-difference algorithm provided by J. Warwicker.18 In this treatment the boundary defining the protein interior, from which salt is excluded, is constructed on a Cartesian grid which extends over a large region. The potential is calculated at all grid points (1 7) Valleau, J. P.; Whittington, S. G . In Modern Theoretical Chemistry; Berne, B. J., Ed.; Plenum: New York, 1977; Vol. 5A. (1 8) Current address: Department of Molecular Biophysics and Biochemistry, Yale University, New Haven, CT 06511.
Ionic Strength Dependence of Superoxide Dismutase
The Journal of Physical Chemistry, Vol. 92, No. 25, 1988 7137
TABLE 111: Least-Squares Fits salt solution concn. M kn.. A-' 0.007 0.0005 0.031 0.009 0.106 0.104 0.234 0.507 0.007 0.0005 0.009
0.031
0.104
0.106
0.507
0.234
charges and coordinates 4. proton units X.Y.Z. A -4 -4 -4 -4 -23.38 9.69 -23.38 9.69 -23.38 9.69 -23.38 9.69
from the finite-difference relation. Local salt ion concentrations are determined at grid points exterior to the protein based on either a linear or an exponential response to the potential. Potentials are recalculated and the salt distribution is adjusted in an iterative fashion until convergence is achieved. We have examined only the 0.104 M bulk salt concentration, using a grid with 87 points in each dimension, spaced at 1-25-A intervals. This provided a water layer on all sides of the protein with a thickness of 23 A. Boundary values for the potential must be specified at the edges of the Cartesian grid. Rather than using zeros or an estimate based on the simplest type of D H treatment, we have obtained boundary values by interpolating from the results of a prior calculation that used a longer, coarser grid. Such a grid is inaccurate near the enzyme but gives excellent results at intermediate distances for test cases we have studied. For the calculations reported in this paper the long grid had a spacing of 2 A and consisted of 87 points in each dimension. This provided a water layer thickness of 56 A. In addition to providing boundary values, the preliminary result was used as an initial guess for the solution on the fine grid. A truncated Taylor series expansion was used for the interpolations. As validation of the numerical procedure and in order to show the ionic strength dependence of the discrepancy between M C and DH, we also exhibit the analytic D H solution for the lSOD model
where (I is the distance of closest approach between an ion and the enzyme, K is the reciprocal Debye screening length which is related to the ionic strength, and the other symbols have their conventional meanings. For our case, a = 31.25 A and q1q2= 4. 3.4. Brownian Dynamics Trajectory Method. Given a force field for enzyme-substrate interactions, a diffusion-limited rate constant for their association can be computed by recently developed methods.% The procedure has been carefully described,M and only a brief sketch is given here. The rate constant is expressed in terms of three principal quantities. The first is kD(b),the rate at which substrate diffuses to within a distance b of the enzyme. If b is chosen in such a way that the potential of mean force, V(r), felt by the substrate at that distance is approximately centrosymmetric, then kD(b)is given by
where DM is the diffusion constant for relative motion of the particles and k T is the product of Boltzmann's constant and the temperature. The second quantity is a, the probability that substrate initially at a distance 6 from the enzyme will react rather than diffuse away to the larger distance q. Reaction is usually defined as contact between specified portions of the two species, although other orientational requirements can be imposed. q is chosen sufficiently
fitting parameters c, A-' a, 8, sum of 0.0104 9.03 0.0492 17.21 0.1771 23.16 0.8135 28.28 0.0053 13.99 0.0046 16.29 0.0240 1.81 -1.37 0.0214 0.1134 12.19 0.0979 6.22 0.3388 20.39 13.71 0.3207
quality of fit sauared residuals no. of mints 0.0084 58 0.0161 53 0.0093 43 0.0030 43 0.9607 836 0.5745
779
0.4723
646
0.6849
646
large compared to b to ensure that the collision probability is a significant fraction of its maximum value. a is measured by repeatedly initiating substrate trajectories and recording whether reaction with enzyme or diffusion to the cutoff distance q occurs in each case. The substrate trajectory is determined by its initial position, the forces exerted on it by the enzyme as modulated by solvent and salt, and a random force representing displacements due to solvent collisions. The last quantity needed in the rate constant expression is Q, the probability that trajectories terminated at a distance q would eventually have reached the distance b rather than diffusing away to infinity. It is given by
It has been shown that the rate constant, k, for a diffusionlimited bimolecular reaction is given by6
4. Results We begin this section by illustrating with smoothed Monte Carlo results the effect of salt on the SOD/O, potential energy surface. We then display the raw M C data, examining its statistical accuracy, the adequacy of the simulation cell dimensions, and the suitability of the chosen fitting function. For the uniformly charged enzyme model, analytic results obtained via the DebyeHuckel approximation are given for each ionic strength. Numerical Debye-Huckel and Poisson-Boltzmann data are then presented for both enzyme models at one salt concentration, and the impact of the approximations made in these theories is explored. Finally, we report rate constants obtained from Brownian dynamics calculations based on the Monte Carlo potential surfaces. The ionic strength dependence of the theoretical results is discussed and contrasted with experimental data. 4.1. Salt Dependence of Steering. That substrate will diffuse toward the active sites of SOD due to its quadrupole moment is clear from Figure 1, which displays the potential energy surface in the absence of counterions or added salt. The precise ionic strength dependence of this surface is not trivially evident, although a gradual flattening is expected. Figures 2 and 3 are based on Monte Carlo simulations of SOD in 5 X lo4 and 0.104 M salt solutions, respectively. Due to the symmetry of the model, it is sufficient to show only one quadrant of a plane that contains the polar axis. At low I , substrate is strongly repelled from most of the enzyme surface. However, approaches along an axis through the active sites face only a modest barrier and then reach an attractive region. At high Z,repulsion from the equatorial domain has weakened and the electrostatic barrier near the axial positions has disappeared. While the overall repulsion between the reacting species has decreased, there is also less steering toward the active sites. Clearly, the behavior of the rate constant depends on the balance between these opposing effects as the salt concentration changes.
7138 The Journal of Physical Chemistry, Vol. 92, No. 25, 1988
Bacquet et al.
A
I
I
I
I
1
70.
90.
110.
130.
0.7
0.e
0.5
zg I
1
Faure 2. Potential energy of 02-for positions on a quadrant of a plane passing through the two active sites of SOD (the only unique region due to the symmetry of the 3SOD model). The displayed surface corresponds to the quadrant farthest from the viewer. Values are given for enzyme-substratecontact (30.0 A) at two points, the active site and midway between active sites. Mesh lines are 3.0 8, apart and are displayed oui to 135.0 A.
0.4
0.3
a
0.2
0.1
0.0
.104M NaCl
30.
50.
DISTANCE (A)
Figure 4. Potential energy of mean force for substrate (or CI-) vs distance
from lSOD model center. Monte Carlo simulation data points are given for four NaCl concentrations: 0.0005M (A),0.009 M ( O ) , 0.104 M (0), 0.507 M (X). Solid lines are fits of eq 2 to these data. Dotted lines are analytic Debye-Hiickel results obtained from eq 3.
Figure 3. Same as Figure 2, but at higher salt concentration.
4.2. Simulation Potentials of Mean Force. The Monte Carlo results for Owcl-(r) were obtained for both the lSOD (uniformly charged) and 3SOD models over a wide range of salt concentrations. Full angular averaging leads to very little noise in the lSOD simulations given in Figure 4. The 3SOD simulations were run for an equal number of passes (see Table 111), but the data can be averaged over only one of the two spherical angles. Figure 5 shows @w(r) for 8 = 90 and 4 arbitrary, what we have termed the equatorial approach to the SOD model. Statistics are still quite good except for the highest concentration where the large number of particles makes a shorter and more noisy simulation necessary. The noise is quite severe for the axial approach (8 = 0) shown in Figure 6, since axial positions correspond to much smaller regions than equatorial positions. A preferential sampling scheme, in which particles near the polar axis have a larger probability of undergoing a trial move, would have partly alleviated this problem but was not used in the present simulations. Nevertheless, the results obtained are sufficient for meaningful estimates of the potential surface even near the active sites. Another concern is that the simulation cell be sufficiently large so that w(r) decays to essentially zero before reaching its boundaries. The final 10-15 A of data is removed from consideration to avoid spurious effects due to the hard outer boundary of the simulation cell. The remaining points shown in Figures 4-6 indicate that adequate cell
radii were chosen in all cases. The 5 X M data are not given in their entirety. For that concentration the grid extended to 340 A, which was fully satisfactory. Display and use of the M C data are greatly facilitated if a functional representation is obtained. As described in section 3, a sum of Yukawa potentials was chosen, with scaled charges and screening parameters determined by nonlinear least-squares fitting and listed in Table 111. The adequacy of this function is evident from Figures 4-6, where good fits are obtained at all positions and concentrations. Note that the fit is a global one and is not necessarily the best local fit to the data at a single angle. Also, no simulation points are shown near contact for the lowest concentration where the fit indicates an attractive region. However, the raw data, which were collected on a finer grid, support the presence of such a feature. The behavior of the potential energy surface near the active site is clear from the smooth fits to the M C data in Figure 6. The surface drops as the salt concentration moves from 0.0005 to 0.009 to 0.104, but the surface rises as the salt molarity increases to 0.507. Such an increase is not surprising since in the infinite salt limit the potential energy is zero everywhere. Near the active site there is considerable noise in the simulations, and the differences among the various curves are small. 4.3. Comparison of Theoretical Results. Recently, finitedifference forms of the Poisson-Boltzmann and Debye-Hiickel equations have been applied to models involving realistic shapes and charge distributions, as well as a spatially varying dielectric respon~e!*~,*-~~ The ability to solve the PB and DH equations for highly detailed systems is an exciting one; however, the errors incurred by neglect of ion-ion correlation and linearization must be considered. There also exists the possibility of purely numerical errors due to the choice of boundary values and grid parameters.
r-i
I
I
I
--IT--
I
0.5
i
0.4
4
0.3
i I
y
4
r
The Journal of Physical Chemistry, Vol. 92, No. 25, 1988 7139
Ionic Strength Dependence of Superoxide Dismutase
1
A
~
5 Y
\
4
0.1
-0.2 -O'l
0.0
30.
50.
70.
eo.
110.
130.
f'
30.
50.
110.
130.
DISTANCE (A)
DISTANCE (A)
Figure 5. Simulation data points and fits to eq 2 for the equatorial approach to the 3SOD model (see section 2). Four NaCl concentrations are given: 0.0005 M (A and ---), 0.009 M ( 0 and ---), 0.104 M (0 and -), and 0.507 M (X and
BO.
70.
Figure 6. Same as Figure 5 , but for the polar approach (see section 2).
.e-).
Here we apply the finite-difference algorithm to the simple models for which M C simulations have been performed, in order to assess the impact of the theoretical approximations involved. It is evident from Figure 7 that, for a 0.104 M salt solution, the numerical D H and PB results agree extremely well with each other and with the curve computed from eq 3, the analytic D H expression for the lSOD model. Thus, the numerical procedure appears satisfactory and, furthermore, linearization has little effect for this ionic strength. For the 0.507 M case displayed in the inset to Figure 7, D H is noticeably inferior to PB, as expected. The 3SOD model results depicted for two angles in Figure 8 also confirm the accuracy of the numerical D H procedure. Judged by the Monte Carlo data, the D H and PB curves in Figures 7 and 8 overestimate w(r). This is the typical consequence of neglecting ion-ion correlation when the correlations are dominated by Coulombic interactions. The opposite effect can sometimes occur due to the fact that ions are treated as point charges with no volume in their interactions with each other (although their size is considered for ion-enzyme interactions). However, this excluded-volume effect is negligible except for very large local ion concentrations and/or ion sizes. The concentration dependence of the discrepancy between M C and DH is displayed in Figure 4 for the ISOD model. As a percentage of the M C result, the difference increases with ionic strength, although its absolute magnitude does not change greatly. 4.4. Rate Constants. The simple model used here is well-suited for the comparison of DH, PB, and M C determinations of the potential of mean force between SOD and its substrate. However, it is only marginally more realistic than some earlier models.' The M C treatment of salt effects is an improvement, but severe approximations (Le., cylindrical symmetry, a spherical enzyme surface, and c = 78.54 inside as well as outside the protein) are still present. The realism of these results can be judged indirectly
0.8
I
0.1
0.0
30.
40.
50.
60.
70.
DISTANCE (A)
Figure 7. Comparison of theoretical results for the lSOD model in 0.104 M (main figure) and 0.507 M (inset) electrolyte solutions. Solid lines are fits of MC data to eq 2, dotted lines are analytic DH results from eq 3, dashed line is the numerical DH result (given only for the 0.104 M case), and dot/dash line is the numerical PB result.
by computing rate constants from them and comparing with experiment. Although significantly better agreement is not expected, it is of interest to determine the extent to which accuracy changes as a result of adopting the present model. Figure 9 displays experimental and theoretical rate constants for the dismutation of 02-by SOD as functions of the ionic strength (I)of the solution. Three different salts were used in the experiments, but very similar results are expected since the ionic strength rather than the particular identity or valence of the salt ions is the crucial quantity. In fact, there is extremely good
Bacquet et al.
7140 The Journal of Physical Chemistry, Vol. 92, No. 25, 1988 \
I
I
I
0.3
0.2
F v
I
o.l O.O
a -0.1
-0.2
eo.
50.
40.
30.
70.
DISTANCE (A)
Figure 8. Comparison of theoretical results for the 3SOD model in 0.104 M electrolyte solution. Upper set of curves is for the equatorial approach, and the lower set is for the polar approach (see section 2). Solid lines are fits of MC data of eq 2, dashed lines are numerical DH results, and dot/dash lines are numerical PB results.
10.0
I---.---*. .......................
8.
' '
,* -x-
".,
+-*-)cex
)c
e.6 8.0
t
4
1 -3.0
-2.0
-1.0
0.0
108 I
Figure 9. Ionic strength dependence of the rate constant. Solid lines are experimental measurements: NaCl from ref 2 (O),Na3PO4from ref 1 (a),and NaC104from ref 2 (-). Dashed lines are computed from simple charge and shape models of SOD: MC data of this paper (+), ref 7 (X), and ref 19 (0). Dotted lines are computed from highly detailed SOD models: ref 9 (A)and ref 10 (#). All curves are merely straight lines connecting measured or computed data points. Zero ionic strength points are plotted on the left border. Units of k are s-l M-I. The data from ref 19 are based on e = 2 for the interior.
agreement for NaCl and Na3FQ4 as measured by pulse radiolysis.' The NaC104 dataZobtained by a polarographic technique differ somewhat, but agreement is sufficient to indicate the qualitative and semiquantitative validity of these measurements. Error bars or related comments were not provided with any of these experimental data sets. We note that the experimental data originally reported in terms of molar concentrations of components have been converted to solution ionic strengths incorrectly in some previous studies? the correct conversion originally given by Allison et a1.I0 has been used here. The theoretical rate constants computed in this study by the methods of section 3.4 are also displayed in Figure 9 along with the results of an earlier study7 which used a more approximate treatment of salt effects. These computed data sets are similar to each other in shape but differ by a factor of approximately 2.
The quantitative difference is due to an improved value of D M , the mutual diffusion constant for SOD and 0;. The change in D M is due to a new value for the hydrodynamic radius of O2-I6 and a change from slip to stick boundary conditions. D M scales the theoretical result and was approximately halved to 0.128 A2/ps. Of more significance is the qualitative similarity of the two curves. One must conclude that the different treatment of salt effects had little impact. The calculated rate constants based on these simple models are only roughly consistent with the experimental data. The experiments do not cover the range of very low I where an increase in the rate constant is computed. As I continues to grow, both theory and experiment indicate a decrease in k but the computed k's are too high and do not begin falling early enough or with sufficient steepness. The only other theoretical study19of a simple SOD model gave results very similar in shape to ours (an increase in k followed by a decrease with small slope), but the rates were higher by a factor of 2-3. D M was not reported and may partly explain the discrepancy. In that model a Debye-Huckel treatment of salt effects was employed and two dielectric values were considered for the protein interior. The realistic value of t = 2 gave rates 4560% higher than the convenient approximation of t = 80. Thus, the use of a better t for the protein interior in the study reported here is unlikely to improve agreement with experiment. Our current results probably benefit from some cancellation among errors. As discussed earlier, two theoretical studies of SOD involving highly detailed descriptions of enzyme shape and electrostatic field and sophisticated treatments of dielectric and salt effects have recently appeared.9,10One studylo produced a rate curve similar at moderate I to the one obtained in the present study but with a more rapid decrease in k in excellent agreement with experiment. That study did not examine extremely small values of I. The other study9 found a much higher and flatter profile for k vs I. The rate of zero ionic strength was computed and exceeded slightly the rate at I = 0.01, in marked contrast to the findings of all three studies involving simple enzyme models. 5. Conclusions This work dealt with the calculation by several different methods of potentials of mean force for salt ions and substrate particles near a model enzyme possessing high symmetry. Monte Carlo simulations and various forms of the Poisson-Boltzmann equation were applied and their feasibility and accuracy assessed. Brownian dynamics methods were then used to obtain the ionic strength dependence of the diffusion-controlled rate of enzyme/substrate reaction. The computed rates were compared with other theoretical results and with experimental data. We conclude that the potential of mean force, w(r),on which substrate diffuses in an enzyme-salt system can be conveniently determined by Monte Carlo (MC) simulation methods when symmetry is present in the model. Although the model is an approximate one, no further approximations are made in order to obtain a solution. Only a statistical uncertainty is present which can be made arbitrarily small by lengthening the simulation. Although better accuracy is available in principle from MC, the Debye-Hiickel (DH) formulation often used for simple models does not differ greatly and is much easier and faster to apply. DH slightly overestimates w(r) for all cases considered here. With careful choices of boundary conditions, grid size, etc., we find excellent agreement between analytic and numerical forms of the Debye-Hiickel equation. Furthermore, at least up to roughly physiological salt concentrations, we find no difference between numerical DebyeHiickel and Poisson-Boltzmann results. These facts, taken with the reasonable agreement with MC, support the validity of studies that apply numerical D H algorithms to highly detailed model^.^^^-^*'^ This finding is perhaps the most important conclusion of the present study. Rates computed from the M C results are not qualitatively different from those reported earlier,7 although a better hydro(19) Head-Gordon,
T.;Brooks, C. L. J . Phys. Chem. 1987, 91, 3342.
J . Phys. Chem. 1988, 92, 7141-7147 dynamic radius for 02-and the use of stick rather than slip boundary conditions have lowered the rates, moving them closer to experimental data. The computed k vs I curve initially rises as I grows from small values, where experimental results are unavailable, and then begins to fall for moderate to high I, in rough accord with experiment. However, the computed decline appears insufficiently steep and does not begin soon enough. The two published theoretical studies based on detailed models of SOD differ from each other and from our simple model results. Relative to the experimental data, one appears much better than the M C
7141
results while the other is perhaps less satisfactory.
Acknowledgment. J.A.M. is the recipient of the G. H. Hitchings Award from the Burroughs Wellcome Fund. S.A.A. is the recipient of an NSF Presidential Young Investigator Award, a Camille and Henry Dreyfus Grant for Newly Appointed Faculty in Chemistry, and a Cottrell Grant from Research Corporation. This work has been supported in part at the University of Houston by the N I H and the Robert A. Welch Foundation. Registry No. SOD,9054-89-1; 02-,11062-77-4.
Thermotropic Phase Transitions in Bis(n-tetradecylammonium) Tetrachlorocadmate( I I ) and Some Homologous Compounds K. J. Schenk and G . Chapuis* UniversitB de Lausanne, Institut de Cristallographie, BSP- Dorigny, CH-1015 Lausanne, Switzerland (Received: April 30, 1987; In Final Form: February 29, 1988)
The room-temperature structure, determined by single-crystalX-ray diffraction, of the title compound is presented. Its four high-temperature phases, identified by means of powder X-ray diffraction and thermal analysis, are also described. Calorimetric results of some bis(n-alkylammonium) tetrachlorocadmates(I1) are interpreted regarding the number (four for n t 13 and and (CloH21NH3)2CdC14 two for 7 I n 5 12) and sequence of phase transitions. From the structures of (C14H29NH3)2CdC14 at room temperature, it is suggested that the number and sequence of phase transformations,observed in the (CfiwlNH3)2CdCb family, depends on the orientations of the n-alkylarnmonium chains in the ground state. For n t 13, the chains adjacent to the CdC14 layers are parallel to one direction and four-phase transitions can be observed, whereas for 7 I n I 12, the chains form a herringbonelike pattern, and two transitions can be detected. Available spectroscopic data combined with structural knowledge are finally woven into a model of the phase transitions. It involves essentially the creation, annihilation, and interconversionof G forms and kinks, the activation of ammonium and methyl rotation, the disordering of the ammonium groups in several potential wells, and the static and dynamical puckering of the CdC14 polyanion.
Introduction The investigation of the bis(n-alkylammonium) tetrahalometallates(II), (CnHZn+lNH3)2MX4 ( M stands for a divalent transition-metal atom, and X for a halogen atom) has greatly contributed to our understanding of phase transitions in layer structures and, by analogy, in two-dimensional lipid bilayers as found in bi~mernbranes.'*'~The advances in synthesis along with the ease of controlling various structural parameters (metal, halogen, number of carbon atoms in the n-alkylammonium ion) have made them ideal objects for studies by spectroscopy, calorimetry, diffraction, and a variety of other techniques. (1) Needham, G. F.; Willett, R. D.; Franzen, H. F. J . Phys. Chem. 1984, 88, 674. (2) Mokhlisse, R.; Chanh, N . B.; Couzi, M.; Haget, Y.; Hauw, C.; Meresse A. J . Phys. C Solid State Phys. 1984, 1 ?, 233. (3) Depmeier, W. Acta Crystallogr.,Sect. B 1981,837, 330; J. Solid State Chem. 1979, 29, 15. (4) Muralt, P.; Kind, R.; Blinc, R.; PekS, B. Phys. Reu. Lett. 1982, 49, 1019. (5) Doudin, B.; Chapuis, G., Acta Crystallogr., Sect. B, in press. (6) Stewart, J. M.; Kundell, F. A.; Baldwin, J. C. X-RAY 72. Technical Report TR-192 of the Computing Center, University of Maryland; version by D. Schwarzenbach, Universite de Lausanne. (7) Cromer, D. T.; Mann, J. D. Acta Crysrallogr.,Sect. A 1968, A24, 321. (8) Cromer, D. T.; Liberman, D. J. Chem. Phys. 1970,53, 1981. (9) Schwarzenbach, D. Abstracts of the 4th European Crystallographic Meeting, Oxford, 1977. (10) Kind, R.; Plesko, S.; Arend, H.; Blinc, R.; ieks, B.; Seliger, J.; Loiar, B.; Slak, J.; Levstik, A.; FilipiE, C.; kagar, V.; Lahajnar, G.; Milia, F.; Chapuis, G. J . Chem. Phys. 1979, 71, 2118. (1 1) Schenk, K. J. ThBse, Universite de Lausanne, Switzerland, 1984. (12) Klyne, W.; Prelog, V. Experientia 1960, 16, 521. (13) Ziiiiiga, F. J.; Chapuis, G. Acta Crystallogr.,Sect. B 1983,839,620. (14) Taupin, D. J . Appl. Cryst. 1968,1, 178, and the routine based on this report, version of February 1, 1974. (IS) Nagle, J. F. J . Chem. Phys. 1973, 58, 252. Nagle, J. F. Ann. Reu. Phys. Chem. 1980, 31, 157.
0022-3654/88/2092-7141$01 S O / O
The crystal structures of (CnH2*lNH3)2 CdC14 (CnCd in this paper) consist of layers of corner-sharing CdCls octahedra alternating with compact layers formed out of n-alkylammonium ions. These ordered layers are transformed into various disordered smectic layers in successive phase transition prior to melting. These thermotropic transformations have been shown to involve disordering of the aliphatic chains, in particular the creation and interconversion of G forms or kinks, disordering of the ammonium groups, activation of ammonium- and methyl-group rotations, and finally static and dynamical puckering of the CdC14 polyanion. This last phenomenon plays an important part in quite challenging diffuse scattering patterns* and even incommensurately modulated
phase^.^-^ The subtle interplay of the effects mentioned above results in quite intriguing phase transformations encompassing continuous, discontinuous, reconstructive, displacive, and order-disorder effects. Their properties produce transition schemes sufficiently complicated to make the CnCd powerful model systems for biological membranes.' Both systems undergo structural phase transitions at comparable temperatures and compositions, and it is tempting to correlate' some of the transition mechanisms deduced from model systems with those observed in lipid bilayers. In this context, the crystal structure of C14Cd at room temperature will be presented along with a characterization of its high-temperature phases as deduced from differential scanning calorimetry (DSC) and X-ray diffraction (XRD) experiments on powders. In addition, the sequence of phases in C14Cd will be compared with those observed in some homologous members of the layer perovskite family. In this report we shall adhere to the convention regarding the labeling of the phases established by previous ~ o r k e r s : ' the *~~~~~~~ (16) White, M. A. J . Chem. Phys. 1984, 81, 6100.
0 1988 American Chemical Society
