J. Phys. Chem. 1994, 98, 10554-10557
10554
Simulation of Superoxide-Superoxide Dismutase Association Rate for Six Natural Variants. Comparison with the Experimental Catalytic Rate Alessandro Sergi* and Mauro Ferrario Dipartimento di Fisica, Universith di Messina, Italy
Fabio Polticelli Dipartimento di Biologia, Universith di Roma “Tor Vergata ”, Italy
Peter O’Neill MRC, Radiobiology Unit, Chilton, Didcot, Oxon OX1 1 ORD, U.K.
Alessandro Desideri Dipartimento di Chimica Organica e Biologica, Universith di Messina, Italy Received: May 16, 1994@
A Brownian dynamics simulation method has been implemented to study the superoxide-superoxide dismutase association reaction. Electrostatic potential and forces are calculated solving the linearized finite difference Poisson-Boltzmann equation. The accuracy of the algorithm has been tested carrying out simulations at different ionic strength values on bovine Cu,Zn superoxide dismutase (SOD), whose three-dimensional structure is known at 2 A resolution, and comparing the calculated association rates with the experimentally determined catalytic rates. Application of the algorithm to six Cu,Zn SOD variants shows that simulations well reproduce the experimental catalytic rate and demonstrates that the diffusion of the superoxide anion to the active site copper is the rate-limiting step in the catalytic process of Cu,Zn superoxide dismutase.
Introduction Cu,Zn superoxide dismutases (Cu,Zn SODS) are a class of enzymes that catalyze the dismutation of the toxic superoxide radical (OZ-)to molecular oxygen and hydrogen per0xide.l The mechanism of action extensively studied on the bovine enzyme, for which a detailed three-dimensional structure exists? has been shown to involve two half cycles of reduction and oxidation of the catalytically active copper ion which constitutes the redox enter.^ For this reaction an enzyme-substrate recognition process based on electrostatic interactions has been hypothesized by chemical modification of positively charged residue^^,^ and by the negative dependence of the enzyme activity on alkaline pH and increasing ionic strength.6 Brownian dynamics simulations have confirmed that electrostatic forces play an important role in steering superoxide into the active site,7t8 being able to reproduce the trend of the ionic strength dependence experimentally observed on the native and chemically modified bovine e n ~ y m e . ~ The J ~ electrostatically driven enzyme-substrate recognition process has been suggested to be unique for all the Cu,Zn SODS on the basis of the similarity of the electrostatic potential distribution” and of the catalytic rate constant value.12 These results suggest that the rate-limiting step of the catalytic reaction is the diffusion process of the incoming superoxide toward the copper active center. In fact, the dismutation reaction can be described by the following scheme:I3
SOD
+ 0,-
ki
0,--SOD
k2
products
(1)
where kl is the second-order association rate, k-1 is the firstorder enzyme-substrate dissociation constant, and k2 is the frst-
* Address for correspondence: Dipartimento di Fisica, Universith di Messina, Sezione Teorica, CP 50, 1-98166 S. Agata di Messina, Italy. Abstract published in Advance ACS Abstracts, September 1, 1994. @
0022-365419412098-10554$04.50/0
order catalytic rate constant. These quantities are related in the Michaelis-Menten model for the enzyme kinetics by the following equation: kf = k21k,,, = klk2/(k-,
+ k2)
(2)
where kf is the second-order enzyme rate constant and k, is the Michaelis-Menten constant. The catalytic process can be defined as diffusion limited, i.e. kf = kl, if k2 >> k-1. Recently the experimental catalytic rate constant has been measured for several natural variants, and in all cases a comparable value has been found. 12,14 In this paper we present the results of Brownian dynamics simulations to calculate the kl association rate constant between 0 2 - and Cu,Zn SOD. In particular, such an approach has been applied to a set of Cu,Zn SODs, having known X-raydetermined or computer-modeled three-dimensional structure, for which the catalytic rate has been experimentally measured. 12,14,15 The results obtained indicate that the calculated kl is very similar for all the systems investigated, in agreement with the experimental data, confirming that the catalytic mechanism of all Cu,Zn SODs is diffusion controlled.
Computational Methods The atomic coordinates of the bovine and human Cu,Zn SOD have been obtained from the Brookhaven Protein Data Bank,16 the coordinates of the yeast Cu,Zn SOD have been kindly provided by Prof. M. Bolognesi, and those of spinach leaves Cu,Zn SOD have been kindly provided by Prof. Y. Kitagawa. The three-dimensional structure of the sheep and the shark enzymes have been computer modeled using the bovine enzyme as a temp1ate.l1J7 0 1994 American Chemical Society
J. Phys. Chem., Vol. 98, No. 41, 1994 10555
Superoxide-Superoxide Dismutase Association Rate
Electrostatic Potential Calculation. Electrostatic potentials around the proteins have been calculated by the program Delphi uses a distributed by Biosym T e c h n o l ~ g i e s .This ~ ~ ~program ~~ continuum macroscopic approach to solve the finite difference Poisson-Boltzmann equation (FDPB) around macromolecules with arbitrary shape and polarity. In detail, the molecules have been mapped into five different cubic grids of 81 x 81 x 81 site points, spaced by 1.0, 1.5, 2.0,4.0, and 8.0 A, respectively; the FDPB has then been solved at each lattice point, using the results obtained with the bigger grids to accurately determine the boundary conditions of the finer ones. The five different potential grids have been used to generate a 81 x 81 x 81 x 5 map from which forces have been calculated. Values of 4 and 78.5 have been used for the dielectric constant of the protein interior and the solvent, respectively. The ion exclusion radius, i.e. the closest approach of solvated ions to the protein surface, was, depending on the simulations, 1 or 2 A. The proteins have been analyzed at pH = 7, considering all lysines and arginines protonated and all carboxyl groups dissociated; a charge of +le, where e is the proton charge, was assigned to the €-nitrogen of each lysine, +0Se to each guanidinium nitrogen, -0.5 to each carboxylate oxygen, and +2e to both the copper and zinc atoms. For all of the six proteins, a charge of -le equally distributed on the two nitrogen atoms of the side chain was given to His61 (numbering of the residues refers to the bovine enzyme), while a charge of + l e was assigned to the imidazole ring of His41 in the mammalian proteins.20 His19 in the bovine enzyme was considered half protonated, and a charge of +0.25e was assigned to each imidazole nitrogen, while His129 in the spinach enzyme was considered fully protonated because of the salt bridge which links this residue to a glutamic acid.21 Brownian Dynamics. The diffusion of superoxide relative to SOD has been evaluated typically on an ensenble of 10 000 different trajectories. The trajectory propagation has been described by the equation
r(t
+ At,) = r(t) + AtpF(t)k,T + s
However, the probability P that the diffusing particle hits the target must be corrected since calculations are carried out in a limited space and there is a finite probability that a trajectory escaped at R2 reenters the simulation space at R1.I0 This correction can be analytically calculated only if the flux at the starting radius R1 has no angular dependence or, in other words, if for r > R1 the potential either vanishes or can be considered centrosymmetric. In the former case, when R1 is large enough to ensure that the potential U % 0 for r > RI, kd(R1) = 4nDR1 and the corrected probability PC is given by the equation
P, =
(3)
where r is the position vector, D is the diffusion constant of superoxide relative to SOD (0.1283 ii/ps,10,22k~ is the Boltzmann constant, T = 298 K, Atg is the time step, F(t) is the electrostatic force applied on superoxide, and s is a vector whose components are independent Gaussian numbers determined by the properties (Si) = 0
each direction has been obtained interpolating over the two faces of the cube perpendicular to this direction and formed by the eight grid points. Trajectory Propagation. A trajectory is typically initiated at R1 = 61.5 and continues until it reacts by entering a sphere of radius 4 A centered on the copper atomlo or escapes by reaching a sphere of R2 = 150 A centered on the protein center. When the diffusing particle collides with a nonreactive area of the protein surface (i.e. outside the metal active site), the particle is “reflected’ starting from the position occupied before collision occurred and translated with an appropriate random vector.* The aim of this procedure is to prevent a particle close to the protein surface from colliding against the protein because of a local strong attractive force, therefore remaining trapped in that position for an unnaturally long time without either reacting with the copper or escaping. Rate Calculation. The association rate kl between superoxide and superoxide dismutase is calculated by multiplying the probability P that a trajectory starting at R1 hits the target by the rate h(R1) of diffusion of particles to the sphere of radius R1:23
P 1 - (1 - P)Rl/R,
(7)
In the latter, the potential outside the sphere of radius R1 is estimated analytically by assuming a centrosymmetric field and representing the protein as a spherical ion of finite radius.24The correction is given by the equation D
(4)
and
where (s;)
= 2DAtg ( i = x , y , z )
(5)
The time step value is directly related to the distance of the superoxide from the center of the protein; in particular, a different time step is used in each grid, its value being 0.125 ps for the finer grid (grid index g = 1) and increasing proportionally with the increase in grid spacing. Calculation of Forces. Electrostatic force components at each grid point have been calculated at the start of the simulation by central differentiation of the potential U(r) over two neighboring points in each direction of the space. In this way, three 79 x 79 x 79 x 5 grids, giving the Cartesian components of the electrostatic field, have been obtained; these grids have been used to calculate by a linear interpolation procedure the forces acting on superoxide at each time step. In detail, the position of the diffusing particle has been used to determine the grid index g, which indicates which field grid has to be used, and the eight nearest neighbor grid points. The force in
and c = q1q2e2exp(kr,)/ek,T( 1
+ kr,)
being ke the Boltzmann constant, T the absolute temperature, k the Debye screening factor, q1q2e2 the product of the charges on superoxide and SOD, ro the radius of SOD represented as a spherical ion, and E the dielectric constant of the solvent. In the calculation reported in this paper, electrostatic potential and forces are accurately calculated by using the PoissonBoltzmann equation throughout the simulation space. Analysis of the calculated electric fields shows that anisotropies in the potential extend to large distances, and the assumption discussed
Sergi et al.
10556 J. Phys. Chem., Vol. 98, No. 41, 1994
TABLE 1: Dependence of the Simulated Association Rate (kl) of Bovine Cu,Zn Superoxide Dismutase on Trajectory Starting (R1)and Terminating (Rz)Radii
1.06 1.09 1.09 1.08 1.12 0.99 1.00
41.5 51.5 61.5 71.5 81.5 101.5 121.5 131.5 201.5
1.07 1.08
t
1.14 1.12
0'85 0.751..
0
1.14 1.06
( x lo-''
M-'
S-')
RdA)
sheep
yeast
41.5 51.5 61.5 71.5 81.5
0.99 1.02 1.03 0.95 1.07
0.81 0.86 0.91 0.96 0.91
TABLE 3: Experimental (kf)and Simulated (kl) Association Rates for Six Different Cu,Zn Superoxide Dismutasesf ox
human yeast sheep
spinach shark
1
.
0.04
.
.
0.08
. . 0.12
I0.05 0.16
Ionic Strength (M)
TABLE 2: Dewndence of the Simulated Association Rate (kl) of Sheep aid Yeast Cu,Zn Superoxide Dismutase on the Trajectory Starting (R1)Radius kl
.
kf
%
ki
%
0.39b 0.25' 0.33b 0.32b 0.39d 0.39e
100 64 85 82 100 100
1.09 0.83 0.91 1.03 0.97 1.22
100 76 83 94 89 112
Catalytic rate constants determined at pH 8 and p = 0.02 M. O'Neil et al., 1988. Getzoff et al., 1992. Polticelli et al.,unpublished results. e Calabrese et al., 1989. f Rates in units of x lo-'' M-I s-l. in the preceding paragraph will lead to a correction, given by eq 8, which is not more exact than the one given by eq 7. Therefore, since electrostatic potential values for r > R1 are in any case small, we used eq 7 to calculate the corrected probability and the rate constant of the encounter between superoxide and SOD. Moreover, carrying out simulations with different values of R1 and R2, we estimated the average relative error in such a procedure, which resulted to be in all cases within the intrinsic statistical variance of the trajectory method (see Tables 1 and 2).
Results and Discussion Test of the Accuracy of the Algorithm. In order to verify the validity of the assumption that V ( r ) 0 for r > R1, the encounter probability between the SOD active site and superoxide has been calculated for three Cu,Zn SODs, having different net charge and surface potential," varying R1 and/or R2.
Table 1 reports the results of the simulations carried out on the bovine Cu,Zn SOD (PI = 5.2, net charge -4 at pH = 7) varying R1, at two different R2 values to quantify the error in the calculation of the association rate constant of the encounter between superoxide and bovine SOD. The average relative error is negligible, being smaller than 3%. As already noticed by previous w o r k e r ~ , ~the J ~ absolute ~ ~ * value of the simulated rates is 2-3 times higher than the experimental ones12 (see Table 3). Some authors attribute this result to the static representation of the macromolecule. A combined approach in which Brown-
Figure 1. Comparison between experimental6catalytic rate kf (0)of bovine Cu,Zn SOD and the simulated association kl with an exclusion radius of 1 (0)or 2 (0)A. The data are reported as a function of ionic strength and refer to different scales because of their difference in the absolute value. ian dynamics and molecular dynamics simulations are coupled has been shown to reduce by 50% the association rate calculated by Brownian dynamics per se.25 Table 2 shows the results obtained on the sheep (PI = 8.0, net charge 0) and yeast (PI = 4.5,net charge -8) Cu,Zn SODs, fixing the R2 value at 150 8, and varying R1; also for these enzymes the average relative error is less than 3% and therefore independent from the particular charge distribution of the molecule analyzed. On the basis of the calculations reported in Tables 1 and 2 we decided that values of R1 = 61.5 A and R2 = 150 A were the optimal choice for obtaining results with the smallest computational effort compatible with the required accuracy. The efficiency of our algorithm has been tested by comparing the calculated ionic strength dependence of the association rate constant of the bovine Cu,Zn SOD with the experimentally measured catalytic rate dependence. The simulated rates, obtained by using two different ion exclusion radius values (1 and 2 A) for the calculation of the electrostatic field, are reported in Figure 1, together with catalytic rate values experimentally obtained by increasing the ionic strength by varying the amount of NaC104 in the assay mixture.6 The results are reported in two different scales because of the above mentioned difference in the absolute values of the experimental and simulated rates. The experimental trend is reproduced by both calculations. The best agreement is obtained using an ion exclusion radius of 2 A, very similar to the radius of the salt used in the determination of the association rate: indicating that the algorithm can reproduce quite well subtle properties of the experimental system, such as the dependence of the association rate on the average radius of the ions in solution. Simulation of the Association Reaction. Bovine Cu,Zn SOD has been shown to display an asymmetric distribution of the potential which has been p r o p o ~ e d ~to~be, ~the ~ ,driving ~~ force guiding the substrate toward the active site. This enzyme is known in fine detail in both its three-dimensional structure2 and kinetic parameters.28 In particular, the reaction rate of bovine Cu,Zn SOD with superoxide is constantly high between pH 5 and pH 10 at low salt concentrations, but decreases with increasing ionic strength and P H . ~This is the expected behavior for the case of a positively charged area attracting the negatively charged substrate into the active site, with the rest of the protein surface being insensitive to ionic modulation because of its constant negative charge above pH 5. However, it was found12 that the same catalytic efficiency and the same pattern of pH and ionic strength dependence are displayed by Cu,Zn superoxide dismutases having more positive or more negative protein
Superoxide-Superoxide Dismutase Association Rate net charge than the bovine enzyme. In fact, an asymmetric distribution of the electrostatic potential has been found for all Cu,Zn SODS investigated.” On the basis of such considerations, the algorithm described in this paper has been applied to the calculation of the association rate constant between superoxide and superoxide dismutase for six different Cu,Zn SODS for which the catalytic rate has been accurately measured, namely, the bovine, the ovine, the human, the yeast, the spinach, and the shark enzymes. The bovine, the human, the yeast, and the spinach Cu,Zn SODS were selected because their three-dimensional structure is known in fine detail by X-ray ~ r y s t a l l o g r a p h y ; ~moreover, ~ ~ ~ - ~ ~ ~ ~the * yeast and spinach enzymes display interesting changes in the charge distribution around the active site, the former lacking Lys 120 and the latter Lys 120, Lys 134, and Glu 131, all residues which have been supposed to play an essential role in steering 0 2 - to the catalytically active ~ o p p e r . ~ The ~ , *selection ~ of the ovine enzyme was dictated by the availability of detailed kinetic datal2 and by the fact that it displays a very high homology (97%)31 with respect to the bovine enzyme while strikingly differs from the acidic bovine protein as far as the charge parameter is concerned; finally, the shark enzyme, whose structure has been computer modeled,17 was chosen because in this enzyme the same net charge as the bovine SOD is obtained with a different global charge distribution. The simulated association rate constants for the six enzymes are reported in Table 3. In all cases the algorithm well reproduces the pattern of the catalytic rate constants experimentally determined by pulse radiolysis,12as evidenced by expressing the experimental and simulated rates as a percentage of that of the bovine enzyme, although, as previously stated, the absolute calculated values are 2-3 times higher than the experimental ones. In particular, the simulated rate of the human enzyme is the lowest, as is also experimentally measured. Inspection of the threedimensional structure and charge distribution of this enzyme does not reveal any evident property that can explain such a result. A 12% deviation between calculated and experimental rates is found in the sheep and the shark enzymes. This result may be attributed to the fact that the three-dimensional structure of these two enzymes has been computer modeled and for this reason is less reliable than those obtained by X-ray diffraction. A very good agreement between the simulated and the experimental value is found for the yeast enzyme, for which the threedimensional structure has been determined at 2.5 8, r e s o l ~ t i o n . ~ ~ In the case of the spinach enzyme, for which the deviation between calculated and experimental values is 11%, the threedimensional structure is also known by X-ray crystallography at 2 8, resolution. On the other hand, in this enzyme, a large rearrangement of the charge distribution around the active site takes place,21which may cause in turn variations of the pKa of charged protein groups with a subsequent uncertainty in the assignment of their net charge. An NMR study on the spinach Cu,Zn SOD is under way in our laboratory in order to correctly assign net charges in this enzyme. Despite the small differences observed between simulated and experimental values, the results reported in Table 3 indicate that Brownian dynamics association rates nicely reproduce the trend of the experimental catalytic rates, c o n f i i n g that the association reaction is the rate-limiting step in the reaction mechanisms of all Cu,Zn superoxide dismutases. A reliable estimate of the catalytic rate can then be obtained by computing the enzyme-
J. Phys. Chem., Vol. 98, No. 41, 1994 10557 substrate association rate by Brownian dynamics simulation. In particular, the results reported indicate that Brownian dynamics simulation can be quantitatively used to predict the properties of these enzymatic systems on a relative basis, since the absolute values do not match the experimental ones, as previously observed by other workers.8,22
Acknowledgment. This research was partially supported by the Italian Research Council via the C.N.R. P.F. “Sistemi Infoxmatici e Calcolo Parallelo” and by the EC program “Human Capital and Mobility”. Dr. Fabio Polticelli is a “Human Capital and Mobility” fellowship holder. References and Notes (1) Bannister, J.; Bannister, W.; Rotilio, G. CRC Biochem. 1987, 22, 111. (2) Tainer, J.; Getzoff, E.; Beem, K.; Richardson, J.; Richardson, D. J. Mol. Biol. 1982, 160, 181. (3) McAdam. M.: Fielden. E.; Lavelle. F.; Calabrese, L.; Cocco, D.; Rotilio, G. Biochem. J. 1977, 167, 271. (4) Cocco, D.: Rossi, L.; Barra, D.; Bossa, F.; Rotilio, G. FEBS Lett. 1982, 150, 303. (5) Cudd, A.; Fridovich, I. J . Biol. Chem. 1982, 257, 11443. (6) Argese, E.; Viglino, P.; Rotilio, G.; Scarpa, M.; Rigo, A. Biochemistry 1987, 26, 3224. (7) Allison. S.: McCammon, J. J . Phvs. Chem. 1985, 89, 1072. (8) Sharp, K.; Fine, R.; Schulten, K.;Honig, B. J . Phys. Chem. 1987, 91, 3624. (9) Sharp, K.; Fine, R.; Honig, B. Science 1987, 236, 1460. (10) Allison, S.; Bacquet, R.; McCammon, J. Biopolymers 1988, 27, 251. (1 1) Desideri, A.; Falconi, M.; Polticelli, F.; Bolognesi, M.; Djinovic, K.; Rotilio, G. J . Mol. Biol. 1992, 223, 337. (12) O’Neill, P.; Davies, S.; Calabrese, L.; Capo, C.; Marmocchi, F.; Natoli, G.; Rotilio, G. Biochem. 1988, 251, 41. (13) Wong, Y.; Clark, T.; Shen, J.; McCammon, J. Mol. Simul. 1993, 10, 277. (14) Calabrese, L.; Polticelli, F.; O’Neill, P.; Galtieri, A,; Barra, D.; SchininB, E.; Bossa, F. FEBS Lett. 1989, 250, 49. (15) Getzoff, E.; Cabelli, D.; Fisher, C.; Parge, H.; Viezzoli, M.; Banci, L.; Hallewell, R. Nature (London) 1992, 358, 347. (16) Bemstein, F.; Koetzle, T.; Williams, G.; Meyer, E., Jr.; Brice, M.; Rodgers, J.; Kennard, 0.;Shimanouchi, T.; Tasumi, M. J. Mol. Biol. 1977, 112, 535. (17) Polticelli, F.; Falconi, M.; O’Neill, P.; Petruzzelli, R.; Galtieri, A.; Lania, A,; Calabrese, L.; Rotilio, G.; Desideri, A. Arch. Biochem. Biophys., in press. (18) Klapper, I.; Hagstrom, R.; Fine, R.; Sharp, K.; Honig, B. Proteins 1986, I, 47. (19) Gilson, M.; Sharp, K.; Honig, B. J. Comput. Chem. 1987, 9, 327. (20) Cass, A.; Allen, H.; Hill, 0.; Smith, B.; Bannister, J.; Bannister, W. Biochemistry 1977, 16, 3061. (21) Kitagawa, Y.; Tanka, N.; Hata, Y.; Kusonoki, M.; Lee, G.; Katsube, Y.; Asada, K.; Aibara, S.; Morita, Y. J . Biochem. 1991,109,477. (22) Sines, J.; Allison, S.; McCammon, J. Biochemistry 1990,29,9403. (23) Northrup, S.; Allison, S.; McCammon, J. J . Chem. Phys. 1984,80, 1517. (24) Northup, S.; Hynes, J. J. Chem. Phys. 1979, 71, 871. (25) Luty, B.; Amrani, S. E.; McCammon, J. J . Am. Chem. SOC. 1993, 115, 11874. (26) Koppenol, W. In In Oxygen and Oxy-radicals in Chemistry & Biolom: Rodeers. M. A. J.. Powers. E. L.. Eds.: Academic Press: New Yorkri981; p’p 671-674. (27) Getzoff, E.; Tainer, J.: Weiner, P.; Kollman, P.; Richardson, J.; Rich&dson, D. Nature (London) 1983, 306, 287. (28) Fielden, E.; Roberts, P. B.; Bray, R.; Lowe, D.; Mautner, G.; Rotilio, G.; Calabrese, L. Biochem. J . 1974, 139, 49. (29) Parge, H.; Hallewell, R.; Tainer, J. Proc. Natl. Acad. Sci. U.S.A. 1992, 89, 6109. (30) Djinovic, K.; Gatti, G.; Coda, A.; Antolini, L.; Pelosi, G.; Desideri, A.; Falconi, M.; Marmocchi, F.; Rotilio, G.; Bolognesi, M. J . Mol. Biol. 1992, 225, 791. (31) SchininB, M. E.; Barra, D.; Gentilomo, S.; Bossa, F.; Capo, C.; Rotilio, G.; Calabrese, L. FEBS Lett. 1986, 7, 207.
