- Email: [email protected]
Molecular dynamics simulation of superoxide interacting with superoxide dismutase
Chemical Physics 158 (1991) 191-198 North-Holland
Molecular dynamics simulation of superoxide interacting with superoxide dismutase Jian Shen and J. Andrew McCammon Lkpartment of Chemistry, University ofHouston, Houston, TX 77204-5641. USA Received 20 May 1991
Molecular dynamics simulations have been used to study the equilibrium distribution of a superoxide substrate molecule (0~ ) in the channel to the active site of the enzyme Cu, Zn superoxide dismutase (SOD). The results are used to consider aspects of the kinetics of SOD.
1. Intro&r&on
The enzyme Cu, Zn superoxide dismutase (SOD) is important in nature, and potentially in pharmaceutical and industrial applications, as a scavenger of superoxide (0, ) [ l-41. SOD catalyzes the destruo tion of 0~ with a rate constant that is high and inversely dependent on solvent viscosity [ 5-7 1. While this suggest that the reaction is diffusion controlled, variations of the rate constant with conservative mutations and substrate concentration suggest that the enzyme operates near the limit of diffusion control, or that the viscosity dependence arises from a mechanistic step other than the initial diffusional encounter of 0, with the enzyme [ 8- 10 1. Even if SOD operates with a conventional diffusion-controlled mechanism in some situations, there could be other situations in which events other than initial encounters with substrate influence the reaction kinetics. For example, the rate of transport of 02 to the enzyme could be increased by modifying the distribution of electrostatic charges on its surface, or by convective effects such as may arise from stirring of the solution. If the transport becomes fast enough, the kinetics of the enzyme will be controlled by mechanistic steps that follow the initial encounter [ 111. In the present paper, we describe preliminary results of molecular dynamics simulations of SOD interacting with 0~. Previous Brownian dynamics 0301-0104/91/$03.50
Fig 1. Schematic view of the active site channe.l in SOD. The approximate locations are shown for two residues that seem to play a particularly important role in the motion ofQ- down the channel.
simulations have described the difisional encounter of 0, with SOD [ 12-161. The present work is intended to provide a framework for understanding the events that follow these initial encounters. Ultimately, we hope to develop a description of the motion of O,- down the 10 A channel from the surface of SOD to the active-site Cu (fig. 1 ), and the subsequent electron and proton transfer steps. This information will be combined with that on the initial encounters to provide a more complete kinetic model. In addition to allowing an analysis of SOD kinetics beyond the diffusion-controlled regime, this model may display features that are important in other en-
0 1991 ElseGer Science Publishers B.V. All rights reserved.
192
J. Shen, J.A. McCammon /&peroxide interacting with superoxide dismutase
zymes, ion channels, etc., in which ions are transported through a molecular channel. Two simulations are described here. The first is with water but no 0, in the SOD channel. This is a refinement of a simulation described previously, in that an improved potential function is used and the simulation covers a longer time [ 17 1. The analysis is limited to features relevant to the interaction with 0,. The second simulation includes the 0, and focuses on the equilibrium distribution of 0, in the channel.
2. Methods The simulation of SOD with water in the channel used methods that are largely identical with those described previously [ 17 1. Here, we only outline the methods and identify the changes that were made. As before, the system consists of protein and solvent water in a 26 A-radius sphere centered on the Cu of the “O-subunit” of SOD. All atoms in the outermost 6 A shell except for the water hydrogens were harmonically restrained, to the X-ray coordinates in the case of protein atoms [ 18 1, and to dynamically relaxed positions from a larger model of the hydrated protein in the case of water molecules. The whole system contains 1629 atoms from SOD (including hydrogens that may participate in hydrogen bonding) and 1780 water molecules. The GROMOS package and potential function were used for all simulations [ 191, but with the following modifications to the potential function used before. The electrostatic charges on the metal atoms and their ligands, previously estimated by an electronegativity equalization argument, were replaced by charges obtained from ab initio calculations [20]. The new model reflects somewhat less charge transfer from the ligands to the metals, with the charges on Cu and Zn increasing from 0.3 for both, to 0.65 and 0.50, respectively. Also, the SPC model for water was replaced by the SPC/E model [ 2 11. As before, the nonbonded interactions were cut off at 8.0 A. (But a longer distance was used when the Oj’ was included, see below. ) After the system was constructed, it was equilibrated by molecular dynamics with reassignment of all velocities from 300 K Maxwellian distributions every 0.5 ps. The time step was 2 fs. The solvent was dynamically re-
laxed for 18 ps, and then the solvent and protein were dynamically relaxed for 15 ps. The latter run was continued for 8 ps with velocity reassignment every 1.Ops. Then the system was coupled to a 300 K bath and simulated for 75 ps; only the last 60 ps were analyzed. The calculations with 0, present started from the final structure of the above simulation. The electrostatic cutoff distance was increased to 20.0 A. The 0, was represented by potential functions that are described elsewhere [22]. A water molecule in the middle of the channel was graduahy transformed into 0, during a 20 ps slow growth process. The remainder of the calculations consisted of a series of umbrella sampling simulations with harmonic restraints on the distance of one oxygen in 02 from the Cu [ Ill. In some simulations, a simultaneous restraint on the 0, distance from Oy of Thr 135 was used to insure overlap of the sampled configurations of 0, in three dimensions. For most simulations, the force constant for the restraint was 9.6 kcal mol-* Ae2. Each umbrella simulation was preceded by a 10.5 ps equilibration run consisting successively of 0.5 ps simulation, Maxwellian reassignment at 300 K every 0.2 ps for 1 ps, Maxwellian reassignment at 300 K every 1.O ps for 4.0 ps, and then 5.0 ps of additional simulation. Then each simulation was continued for 20 ps during which coordinates were saved every 10 Is for further analysis. Selected simulations were extended for as much as an additional 60 ps. In the first umbrella simulation, the minimum of the restraint potential was 6.5 A from the Cu. In subsequent simulations, this minimum was typically shifted in increments of 0.5 A to sample regions closer to and farther from the Cu. The configurations of 0, from each simulation were graphically examined to check qualitatively on the uniformity of their distribution and their overlap with previous distributions. In a few regions where the average force on 05 was relatively large, increments smaller than 0.5 A and restraint force constants twice as large as usual were needed to get adequate sampling. Also, in one region between Thr 135 and Arg 141, it was necessary to use two restraining potentials because the apparent path preferred by 0, deviated substantially from the radial coordinate from Cu. By using a second restraining potential based on the distance from Or:,, it was possible to generate a set of umbrella sampling distribu-
J. Shen, J.A. McCammon /Superoxide interacting with superoxide dismutase
tions that had adequate overlap in three dimensions. The total number of umbrella sampling runs used here is 21, corresponding to a total simulation time of 220.5 ps (equilibration) plus 640 ps (data collection). Each of the umbrella simulations provided a local estimate of the probability distribution p(r) of one of the atoms of 0, in some part of the SOD channel [ 111. To develop a continuous description of p( r) throughout the channel, the following procedure was used. The coordinates of the reference oxygen atom of 0, were binned with a grid spacing of 0.5 A. The resulting distribution p*(r) is biased by the restraining potential u(r). This bias was removed by using the relationship p(r) =~*(r) exp t -b(r)
(1)
1.
The local probability distributions obtained in this way still differ by constant factors that are unknown. Starting with the distribution pI (r) in the most densely sampled region, a spatially overlapping distribution pZ(r) was selected and scaled to be consistent with the first distribution by finding the scaling factor S that satisfies 1 [PI (0 -%(i)
1[n*(h(O
12=0*
(2)
Here, the sum is over the grid elements and ni( i) indicates the number of times the tagged oxygen was found in grid element i during the simulation that produced the probability distribution i. The factor involving the nj( i) is designed to increase the weight of grid points with the greatest joint sampling in the evaluation of S. The two probability distributions are then merged into an extended distribution through p(i)
=
pl(ih
(i) +@2(On2(i) n,(i)+n2(i)
’
(3)
where the ni( i) are used to increase the weight of the most heavily sampled local distribution in forming the average for the extended distribution. Eqs. (2) and ( 3 ) are used to merge the remaining local probability distributions successively into an extended distribution of increasing range. Finally, the potential of mean force w(r) for the tagged oxygen is evaluated from w(r)=-kTlnp(r)+C,
(4)
193
where kT is Boltzmann’s constant times absolute temperature, and C is an arbitrary constant [ 111. For the above procedure to work, there must be a continuous chain of substantially overlapping local distributions from the original umbrella sampling runs. That this is true in the present case is illustrated in fig. 2. Also, the relative values of the potential of mean force in different regions of space must be reasonably insensitive to the particular sequence chosen for merging local probability distributions. In the present case, the choice of different sequences leads to differences of less than 1 kcal/mol for points separated by about 3 A, but may exceed 1 kcal/mol for widely separated points. This is adequate for the semiquantitative analysis presented here, which focuses on local barriers and wells in the potential of mean force.
3. I&salts and discus~on Although the initial simulation did not include 02, there are a number of results that are relevant to the interaction of SOD with its substrate. We focus on those features here; many additional results from the initial simulation are described elsewhere [ 23 1. In the second part of this section, we present equilibrium results of the simulation that included 02. A first concern in any simulation of a protein is whether the structure of the model remains reasonably close to that of the X-ray structure. In the present case, the essential features of the X-ray structure appear to be preserved. The average coordinates of the unrestrained atoms of the protein during the final 10 ps of the simulation have an overall rms deviation from the X-ray structure of 1.4 A. The deviation is largest for residues at the surface of the protein, but is reduced for residues that are more than 15 A from the Cu due to the restraints applied in the outer shell. The overall average rms amplitude of atomic fluctuation about the average coordinates during the final 10 ps is 0.5 A, rising to nearly 1.0 A at the protein surface. These results are quite similar to those seen in other simulations of hydrated proteins [ 111. Again, suppression of the fluctuations is apparent for residues that are more than 15 A from the Cu. This should not markedly affect the analysis of channel dynam-
194
J. Shen, J.A. McCammon /&peroxide interacting with superoxide dismutase
Fig. 2. Stereo view of the sampling density inside the channel of SOD. The total number of 02 configurations recorded in each bin (with edge lengths of 0.5 A) for all ofthe umbrella sampling runs is contoured at levels of 5 (dashed line) and 50 (solid line). A representative instantaneous structure from the middle of the calculations is superposed. The residues Thr 135 and Arg 14 1 arc in front and in back of the contours, respectively.
its, since the channel only extends about 10 A from the Cu. Turning to more specific features of the channel and the active site, several substantial displacements from the X-ray structure are seen at the opening of the channel at the protein surface. The backbone of residues 136 through 139 in one rim of the channel shifts by about 3 A, possibly due in part to the insertion of a bridging water molecule into what had been a direct hydrogen bond between the N of Lys 120 and the 0 of Ala 138. In the other rim, the backbone of residues 56 through 59 shifts by about 2 A.The net effect is a slight narrowing of the channel mouth. The amplitudes of fluctuation in the backbone of the channel rims are about 1.2 A during a period of 60 ps. Among the sidechains at the outer part of the channel, the most notable displacement is in Lys 134, whose NC atom moves from a position in tbe channel and 7.6 A from the Cu, to the external solvent at 13.6 A from the Cu. At the other end of the channel, near the active site, the displacements and fluctuations are smaller, with typical magnitudes of 1.O A and 0.5 A, respectively., The structure of the metal-ligand system is generally similar to that in the X-ray structure. The four histidine nitrogens that are coordinates to the Cu are slightly displaced toward a planar structure, in part perhaps because the solvent water in the fifth coordination site moves in toward the Cu to strengthen its interaction with the metal (see below).
The solid angle at the Cu formed by the surface that has these four nitrogens as vertices is reduced from 4.8 sr in the X-ray structure to 4.5 sr, corresponding to a greater exposure of Cu to solvent in the channel. This solid angle fluctuates between about 3.0 and 5.5 sr, with a strong frequency component near 250 cm-‘. These fluctuations are correlated with fluctuations in the distance from Cu to the oxygen of the bound water; the latter distance ranges from 2.1 to 2.9 A during the last 30 ps of the simulation. The solvent in the channel has a pronounced structure near the Cu, as is evident in the radial distribution function shown in fig. 3. The water molecule bound to the Cu is quite stable. It does not exchange at all during the last 30 ps of the simulation, but rather oscillates around its average distance of 2.4 A from the Cu. The existence of a strongly bound water molecule is consistent with magnetic resonance evidence for such a molecule with a residence time between 4x 10m6 s and 10e8 s with a Cu-H distance of 2.8 kO.4 A [24]. Assuming the water molecule has its dipole pointing radially away from the Cu (as found in the simulation), the distance between the Cu and the water oxygen would be about 2.1 A,within the range found here. The radial distribution function also shows a second binding site for a water molecule between 3.5 and 4.5 A from the Cu. The water molecules that occupy this site do exchange several times during the simulation with other water mole-
J. Shen, J.A. McCammon /Superoxide interacting with superoxide dismutase
0.
3.
6.
9.
12.
15.
18.
r, A
Fig. 3. Radial distribution functions for distances between Cu and the oxygen atoms (heavy lines) and hydrogen atoms (thin lines) of water molecules in the channel. The dashed lines are the cumulative numbers of oxygen or hydrogen atoms obtained by integrating the radial distribution functions. These data are from the final 30 ps of the simulation without 0~ present.
cules in the channel. The number of water molecules in the channel increases steadily with distance beyond 5 A from the Cu. On average, there are three, five, and ten water molecules within 6,7, and 8 A of Cu. It is also useful for future reference to consider the hydration of two important polar sidechains in the channel, namely Thr 135 and Arg 14 1. The radial distribution functions ofwater molecules around key atoms in these sidechains have been computed. By defining the first hydration shell in terms of the first minimum in each radial distribution function, the numbers of water molecules in these shells were computed as 3.9, 6.1, and 1.3 for IV&, Nl&, and O&, respectively. The average residence times of water molecules in these shells were slightly less than 1 ps, but certain water molecules remain bound for periods as long as 30 ps. The major result of the simulations with 02 is the potential of mean force (pmf ) for 0, displacement in the SOD channel. It is not sufficient to use a single radial coordinate from Cu to construct this pmf, because a small variation in such a pmf could mask the occurrence of free energy barriers for non-radial displacements along the preferred pathways for 02 motion. The pmf is therefore represented as a function w(r) of the coordinates of one of the oxygen atoms in 0,. The pmf is shown in figs. 4 and 5 from two orthogonal directions. In fig. 4, the viewpoint is the
195
same as in fig. 2, with the channel running from the region of the Cu up to the left and away from the viewer. In fig. 5, the channel runs from the Cu up to the left but toward the viewer. As in fig. 2, a single instantaneous structure of the protein is included in the figures. A detailed examination of the pmf indicates that the motion of 0, from the surface of SOD to the vicinity of the Cu is not a simple “downhill” process in terms of the free energy change. One possible path involves the following steps. The 0, initially binds between the backbone rim segments formed by residues 56 through 59 and residues 138 through 141. This corresponds to the small low energy patch in the upper left of fgs. 4a and 5a. From there, the O,crosses a free energy barrier of about 1 kcal/mol to move into the broader, deeper minimum between the sidechains of Thr 135 and Arg 14 1. Finally, the O,crosses a barrier of about 2 to 3 kcal/mol to displace the water molecule bound to the Cu. The minimum free energies in the initial binding site and in the Cu binding site are, respectively, about 1 kcal/mol and 2 kcal/mol higher than in the primary minimum between the Thr 135 and Arg 14 1 sidechains. Although the pmf does not decrease monotonically in going from the mouth of the channel to the Cu, the variations in free energy are not large. This fact, together with the significant number of local minima indicated by figs. 4a and 5a, suggests that substrate molecules are likely to follow a variety of paths that differ in their details in moving through the SOD channel. The dynamics of these motions are also likely to be quite complicated, involving the rearrangement of hydrogen bonds among the substrate, water molecules and the channel residues. These rearrangements will be coupled to the fluctuations of the protein itself Indications of these complications are evident in very preliminary analyses of the average interactions of 02 with water molecules and protein atoms for different positions of Ot in the channel, and of trajectories in which the 02 is allowed to wander in the channel with no restraining potentials. The former analyses suggest that the initial binding to SOD may involve a direct interaction between 0, and the guanidinium moiety of Arg 14 1; the latter appears to move toward the mouth of the channel and to replace water molecules in the first hydration shell of 0,. In the broad minimum of the pmf be-
196
J. Shen, J.A. McCammon /Superoxide interacting with superoxide dismutase
(4
Fig. 4. Stereo views of the 05 potential of mean force (pmf ), viewed from the same direction as in fig. 2. The channel runs from the region of the Cu up to the left and away from the viewer. The contour levels are (a) 4 kcal/mol, and (b) 6 kcal/mol.
tween the Thr 135 and Arg 14 1 sidechains, the hydration number of O,- may decrease somewhat further, although the 02 may have bridging water molecules to one or both of the Thr 135 and Arg 14 1 sidechains. In one unrestrained dynamics trajectory, one can see an 02 associating with the Arg 141 sidechain and then moving to form a fairly long-lived ( > 5 ps) structure with bridging waters to one NH2 group of Arg 141 and a direct interaction with the OH of Thr 135. Much additional work must be done to analyze such details, however. Finally, it is of interest to consider the general implications of the present results in terms of the kinetics of SOD. Suppose that the steps in the SOD catalyzed reaction can be represented as SOD+06
kl eSOD//02
ka =SOD/O1-...
ing site, respectively. For simplicity, we assume for now that k3 is large enough to be neglected in estimating the overall rate. This is suggested by the fact that the dismutation of 01 by Cu2+ in water is diffusion controlled [ 25 1, although more analysis of the situation in the enzyme is required. Using Brownian dynamics to estimate k,, and using the volume sampled by 0, at the outer part of the channel to estimate k, /k_ ,, rough values for k, and k_ , are 6 x 1O9 M-’ s-’ and 1.6~ 10” s-‘, respectively. If the limiting step in moving through the channel is crossing a free energy barrier of 3 kcal/mol, the largest that k2 could be is 4.0 x 10’ 1s- I. The effective forward rate constant, k,k2 +kz
lir=k_,
b
’
(6)
,
k-1
(5) where // and / indicate structures in which the 02 is at the outer end of the channel and at the Cu bind-
would then be 0.7k, =4x lo9 M-’ s-‘, in the range of the experimental results, and the reaction would essentially be diffusion controlled. The channel crossing rate constant k2 is, however, very likely to be reduced somewhat by frictional effects and perhaps
J. Shen, J.A. McCammon /Superoxide interacting with superoxide dismutase
197
(4
M
Fig. 5. Stereo views ofthe 0~ pmffrom a direction orthogonal to that in f@. 4. The channel runs from the Cu up to the left but toward the viewer. The contour levels are (a) 4 kcal/mol, and (b) 6 kcal/mol.
by gating effects due to the motion of the protein [ 111. These issues are being analyzed in current work. Other issues that should be analyzed in the future include the effects of including electronic polarizability in the simulations, and the nature of the electron and proton-transfer steps at the end of the reaction.
Acknowledgement This work has been supported in part by NIH, NSF, the Robert A. Welch Foundation, the Texas Advanced Research Program, and the National Center for Super-computing Applications. JAM is the recipient of the George H. Hitchings Award from the Burroughs Wellcome Fund. We also thank Wilfred van Gunsteren for providing the GROMOS program.
References [ 1] I. Fridovich, in: Biochemical and Medical Aspects of Active Oxygen, eds. 0. Hayaishi and K. Asada (Univ. Park Press, Baltimore, 1977) pp. 3-12,171. [2] J.A. Fee, in: Metal Ion Activation of Dioxygen, ed. T.G. Spiro (Wiley, New York, 1980) p. 209. [3]R.A.GreenwaldandG.Cahen,eds,OxyRadicalsandtheir Scavenger Systems, Vol. II, Cellular and Medical Aspects (Elsevier, Amsterdam, 1983). [4]J.R. Stewart, W.H. Merill and W.H. Frist, in Oxygen RadicaB in Biology and Medicine, eds. M.G. Simic, K.A. Taylor, J.F. Ward and C. von Sonntag (Plenum Press, New York, 1988) p. 957. (5lE.M. Fielden, P.B. Roberts, R.C. Bray, D.J. Lowe, G.N. Mautner, G. Rotilio and L. Calabrese, Biochem. J. 139 (1974) 49: [6] A. Cudd and I. Fridovich, J. Biol. Chem. 257 (1982) 11443. [ 71 E. Argese, P. Viglino, G. Rotilio, M. Scarpa and A. Rigo, Biochemistry 26 (1987) 3224. [8] C. Bull and J.A. Fee, J. Am. Chem. Sot. 107 (1985) 3295. [9] J.A. FeeandC.Bull, J.Biol.Chem. 261 (1986) 13QOO. [ IO] L. Banci, I. Bertini, C. Luchinat and R.A. Hallewell, J. Am. Chem.Soc. 110 (1988) 3629. [ 111 J.A. McCammon and S.C. Harvey, Dynamics of Proteins and Nucleic Acids (Cambridge Univ. Press, Cambridge, 1987).
198
J. Shen, J.A. McCammon /Superoxide interacting with superoxide dismutase
[ 12 ] S.A. Allison and J.A. McCammon, J. Phys. Chem. 89 (1985) 1072. [ 131 K. Sharp, R. Fine and B. Honig, Science 236 (1987) 1460. [ 14) T. HeadGordonand C.L. Brooks, J. Phys. them. 91 (1987) 3342. [ 15 ] S.A. Allison, R.J. Banquet and J.A. McCammon, Biopolymers 27 (1988) 251. [ 161 J.J. Sines, S.A. Allison and J.A. McCammon, Biochemistry 29 ( 1990) 9403. [ 171 J. Shen, S. Subramaniam, C.F. Wongand J.A. McCammon, Biopolymers 28 ( 1989) 2085. [ 181 J.A. Tainer, E.D. Getzoff, KM. Beem, J.S. Richardson and D.C. Richardson, J. Mol. Biol. 160 (1982) 181.
[ 191 H.J.C. Beret&en and W.F. van Gunsteren, personal communication. [20] J. Shen, S. Subramaniam, C.F. Won& T. Albright and J.A. McCammon, J. Comput. Chem. 11 ( 1990) 346. [ 2 1 ] H.J.C. Berendsen, J.R. Gigera and T.P. Straatsma, J. Phys. Chem. 91 (1987) 6269. [ 22 ] J. Shen, C.F. Wong and J.A. McCammon, J. Comput. Chem. II (1990) 1003. [23] J. Shen, Ph.D. Thesis (University of Houston, 1990). [24] B.P.Gaber,R.D.Brown,S.H.KoenigandJ.A.Fee,Biochim. Biophys. Acta 271 (1972) 1. [ 251 J. Rabani, D. Klug-Roth and J. Lilic, J. Phys. Chem. 77 (1973) 1169.
Molecular dynamics simulation of superoxide interacting with superoxide dismutase Jian Shen and J. Andrew McCammon Lkpartment of Chemistry, University ofHouston, Houston, TX 77204-5641. USA Received 20 May 1991
Molecular dynamics simulations have been used to study the equilibrium distribution of a superoxide substrate molecule (0~ ) in the channel to the active site of the enzyme Cu, Zn superoxide dismutase (SOD). The results are used to consider aspects of the kinetics of SOD.
1. Intro&r&on
The enzyme Cu, Zn superoxide dismutase (SOD) is important in nature, and potentially in pharmaceutical and industrial applications, as a scavenger of superoxide (0, ) [ l-41. SOD catalyzes the destruo tion of 0~ with a rate constant that is high and inversely dependent on solvent viscosity [ 5-7 1. While this suggest that the reaction is diffusion controlled, variations of the rate constant with conservative mutations and substrate concentration suggest that the enzyme operates near the limit of diffusion control, or that the viscosity dependence arises from a mechanistic step other than the initial diffusional encounter of 0, with the enzyme [ 8- 10 1. Even if SOD operates with a conventional diffusion-controlled mechanism in some situations, there could be other situations in which events other than initial encounters with substrate influence the reaction kinetics. For example, the rate of transport of 02 to the enzyme could be increased by modifying the distribution of electrostatic charges on its surface, or by convective effects such as may arise from stirring of the solution. If the transport becomes fast enough, the kinetics of the enzyme will be controlled by mechanistic steps that follow the initial encounter [ 111. In the present paper, we describe preliminary results of molecular dynamics simulations of SOD interacting with 0~. Previous Brownian dynamics 0301-0104/91/$03.50
Fig 1. Schematic view of the active site channe.l in SOD. The approximate locations are shown for two residues that seem to play a particularly important role in the motion ofQ- down the channel.
simulations have described the difisional encounter of 0, with SOD [ 12-161. The present work is intended to provide a framework for understanding the events that follow these initial encounters. Ultimately, we hope to develop a description of the motion of O,- down the 10 A channel from the surface of SOD to the active-site Cu (fig. 1 ), and the subsequent electron and proton transfer steps. This information will be combined with that on the initial encounters to provide a more complete kinetic model. In addition to allowing an analysis of SOD kinetics beyond the diffusion-controlled regime, this model may display features that are important in other en-
0 1991 ElseGer Science Publishers B.V. All rights reserved.
192
J. Shen, J.A. McCammon /&peroxide interacting with superoxide dismutase
zymes, ion channels, etc., in which ions are transported through a molecular channel. Two simulations are described here. The first is with water but no 0, in the SOD channel. This is a refinement of a simulation described previously, in that an improved potential function is used and the simulation covers a longer time [ 17 1. The analysis is limited to features relevant to the interaction with 0,. The second simulation includes the 0, and focuses on the equilibrium distribution of 0, in the channel.
2. Methods The simulation of SOD with water in the channel used methods that are largely identical with those described previously [ 17 1. Here, we only outline the methods and identify the changes that were made. As before, the system consists of protein and solvent water in a 26 A-radius sphere centered on the Cu of the “O-subunit” of SOD. All atoms in the outermost 6 A shell except for the water hydrogens were harmonically restrained, to the X-ray coordinates in the case of protein atoms [ 18 1, and to dynamically relaxed positions from a larger model of the hydrated protein in the case of water molecules. The whole system contains 1629 atoms from SOD (including hydrogens that may participate in hydrogen bonding) and 1780 water molecules. The GROMOS package and potential function were used for all simulations [ 191, but with the following modifications to the potential function used before. The electrostatic charges on the metal atoms and their ligands, previously estimated by an electronegativity equalization argument, were replaced by charges obtained from ab initio calculations [20]. The new model reflects somewhat less charge transfer from the ligands to the metals, with the charges on Cu and Zn increasing from 0.3 for both, to 0.65 and 0.50, respectively. Also, the SPC model for water was replaced by the SPC/E model [ 2 11. As before, the nonbonded interactions were cut off at 8.0 A. (But a longer distance was used when the Oj’ was included, see below. ) After the system was constructed, it was equilibrated by molecular dynamics with reassignment of all velocities from 300 K Maxwellian distributions every 0.5 ps. The time step was 2 fs. The solvent was dynamically re-
laxed for 18 ps, and then the solvent and protein were dynamically relaxed for 15 ps. The latter run was continued for 8 ps with velocity reassignment every 1.Ops. Then the system was coupled to a 300 K bath and simulated for 75 ps; only the last 60 ps were analyzed. The calculations with 0, present started from the final structure of the above simulation. The electrostatic cutoff distance was increased to 20.0 A. The 0, was represented by potential functions that are described elsewhere [22]. A water molecule in the middle of the channel was graduahy transformed into 0, during a 20 ps slow growth process. The remainder of the calculations consisted of a series of umbrella sampling simulations with harmonic restraints on the distance of one oxygen in 02 from the Cu [ Ill. In some simulations, a simultaneous restraint on the 0, distance from Oy of Thr 135 was used to insure overlap of the sampled configurations of 0, in three dimensions. For most simulations, the force constant for the restraint was 9.6 kcal mol-* Ae2. Each umbrella simulation was preceded by a 10.5 ps equilibration run consisting successively of 0.5 ps simulation, Maxwellian reassignment at 300 K every 0.2 ps for 1 ps, Maxwellian reassignment at 300 K every 1.O ps for 4.0 ps, and then 5.0 ps of additional simulation. Then each simulation was continued for 20 ps during which coordinates were saved every 10 Is for further analysis. Selected simulations were extended for as much as an additional 60 ps. In the first umbrella simulation, the minimum of the restraint potential was 6.5 A from the Cu. In subsequent simulations, this minimum was typically shifted in increments of 0.5 A to sample regions closer to and farther from the Cu. The configurations of 0, from each simulation were graphically examined to check qualitatively on the uniformity of their distribution and their overlap with previous distributions. In a few regions where the average force on 05 was relatively large, increments smaller than 0.5 A and restraint force constants twice as large as usual were needed to get adequate sampling. Also, in one region between Thr 135 and Arg 141, it was necessary to use two restraining potentials because the apparent path preferred by 0, deviated substantially from the radial coordinate from Cu. By using a second restraining potential based on the distance from Or:,, it was possible to generate a set of umbrella sampling distribu-
J. Shen, J.A. McCammon /Superoxide interacting with superoxide dismutase
tions that had adequate overlap in three dimensions. The total number of umbrella sampling runs used here is 21, corresponding to a total simulation time of 220.5 ps (equilibration) plus 640 ps (data collection). Each of the umbrella simulations provided a local estimate of the probability distribution p(r) of one of the atoms of 0, in some part of the SOD channel [ 111. To develop a continuous description of p( r) throughout the channel, the following procedure was used. The coordinates of the reference oxygen atom of 0, were binned with a grid spacing of 0.5 A. The resulting distribution p*(r) is biased by the restraining potential u(r). This bias was removed by using the relationship p(r) =~*(r) exp t -b(r)
(1)
1.
The local probability distributions obtained in this way still differ by constant factors that are unknown. Starting with the distribution pI (r) in the most densely sampled region, a spatially overlapping distribution pZ(r) was selected and scaled to be consistent with the first distribution by finding the scaling factor S that satisfies 1 [PI (0 -%(i)
1[n*(h(O
12=0*
(2)
Here, the sum is over the grid elements and ni( i) indicates the number of times the tagged oxygen was found in grid element i during the simulation that produced the probability distribution i. The factor involving the nj( i) is designed to increase the weight of grid points with the greatest joint sampling in the evaluation of S. The two probability distributions are then merged into an extended distribution through p(i)
=
pl(ih
(i) +@2(On2(i) n,(i)+n2(i)
’
(3)
where the ni( i) are used to increase the weight of the most heavily sampled local distribution in forming the average for the extended distribution. Eqs. (2) and ( 3 ) are used to merge the remaining local probability distributions successively into an extended distribution of increasing range. Finally, the potential of mean force w(r) for the tagged oxygen is evaluated from w(r)=-kTlnp(r)+C,
(4)
193
where kT is Boltzmann’s constant times absolute temperature, and C is an arbitrary constant [ 111. For the above procedure to work, there must be a continuous chain of substantially overlapping local distributions from the original umbrella sampling runs. That this is true in the present case is illustrated in fig. 2. Also, the relative values of the potential of mean force in different regions of space must be reasonably insensitive to the particular sequence chosen for merging local probability distributions. In the present case, the choice of different sequences leads to differences of less than 1 kcal/mol for points separated by about 3 A, but may exceed 1 kcal/mol for widely separated points. This is adequate for the semiquantitative analysis presented here, which focuses on local barriers and wells in the potential of mean force.
3. I&salts and discus~on Although the initial simulation did not include 02, there are a number of results that are relevant to the interaction of SOD with its substrate. We focus on those features here; many additional results from the initial simulation are described elsewhere [ 23 1. In the second part of this section, we present equilibrium results of the simulation that included 02. A first concern in any simulation of a protein is whether the structure of the model remains reasonably close to that of the X-ray structure. In the present case, the essential features of the X-ray structure appear to be preserved. The average coordinates of the unrestrained atoms of the protein during the final 10 ps of the simulation have an overall rms deviation from the X-ray structure of 1.4 A. The deviation is largest for residues at the surface of the protein, but is reduced for residues that are more than 15 A from the Cu due to the restraints applied in the outer shell. The overall average rms amplitude of atomic fluctuation about the average coordinates during the final 10 ps is 0.5 A, rising to nearly 1.0 A at the protein surface. These results are quite similar to those seen in other simulations of hydrated proteins [ 111. Again, suppression of the fluctuations is apparent for residues that are more than 15 A from the Cu. This should not markedly affect the analysis of channel dynam-
194
J. Shen, J.A. McCammon /&peroxide interacting with superoxide dismutase
Fig. 2. Stereo view of the sampling density inside the channel of SOD. The total number of 02 configurations recorded in each bin (with edge lengths of 0.5 A) for all ofthe umbrella sampling runs is contoured at levels of 5 (dashed line) and 50 (solid line). A representative instantaneous structure from the middle of the calculations is superposed. The residues Thr 135 and Arg 14 1 arc in front and in back of the contours, respectively.
its, since the channel only extends about 10 A from the Cu. Turning to more specific features of the channel and the active site, several substantial displacements from the X-ray structure are seen at the opening of the channel at the protein surface. The backbone of residues 136 through 139 in one rim of the channel shifts by about 3 A, possibly due in part to the insertion of a bridging water molecule into what had been a direct hydrogen bond between the N of Lys 120 and the 0 of Ala 138. In the other rim, the backbone of residues 56 through 59 shifts by about 2 A.The net effect is a slight narrowing of the channel mouth. The amplitudes of fluctuation in the backbone of the channel rims are about 1.2 A during a period of 60 ps. Among the sidechains at the outer part of the channel, the most notable displacement is in Lys 134, whose NC atom moves from a position in tbe channel and 7.6 A from the Cu, to the external solvent at 13.6 A from the Cu. At the other end of the channel, near the active site, the displacements and fluctuations are smaller, with typical magnitudes of 1.O A and 0.5 A, respectively., The structure of the metal-ligand system is generally similar to that in the X-ray structure. The four histidine nitrogens that are coordinates to the Cu are slightly displaced toward a planar structure, in part perhaps because the solvent water in the fifth coordination site moves in toward the Cu to strengthen its interaction with the metal (see below).
The solid angle at the Cu formed by the surface that has these four nitrogens as vertices is reduced from 4.8 sr in the X-ray structure to 4.5 sr, corresponding to a greater exposure of Cu to solvent in the channel. This solid angle fluctuates between about 3.0 and 5.5 sr, with a strong frequency component near 250 cm-‘. These fluctuations are correlated with fluctuations in the distance from Cu to the oxygen of the bound water; the latter distance ranges from 2.1 to 2.9 A during the last 30 ps of the simulation. The solvent in the channel has a pronounced structure near the Cu, as is evident in the radial distribution function shown in fig. 3. The water molecule bound to the Cu is quite stable. It does not exchange at all during the last 30 ps of the simulation, but rather oscillates around its average distance of 2.4 A from the Cu. The existence of a strongly bound water molecule is consistent with magnetic resonance evidence for such a molecule with a residence time between 4x 10m6 s and 10e8 s with a Cu-H distance of 2.8 kO.4 A [24]. Assuming the water molecule has its dipole pointing radially away from the Cu (as found in the simulation), the distance between the Cu and the water oxygen would be about 2.1 A,within the range found here. The radial distribution function also shows a second binding site for a water molecule between 3.5 and 4.5 A from the Cu. The water molecules that occupy this site do exchange several times during the simulation with other water mole-
J. Shen, J.A. McCammon /Superoxide interacting with superoxide dismutase
0.
3.
6.
9.
12.
15.
18.
r, A
Fig. 3. Radial distribution functions for distances between Cu and the oxygen atoms (heavy lines) and hydrogen atoms (thin lines) of water molecules in the channel. The dashed lines are the cumulative numbers of oxygen or hydrogen atoms obtained by integrating the radial distribution functions. These data are from the final 30 ps of the simulation without 0~ present.
cules in the channel. The number of water molecules in the channel increases steadily with distance beyond 5 A from the Cu. On average, there are three, five, and ten water molecules within 6,7, and 8 A of Cu. It is also useful for future reference to consider the hydration of two important polar sidechains in the channel, namely Thr 135 and Arg 14 1. The radial distribution functions ofwater molecules around key atoms in these sidechains have been computed. By defining the first hydration shell in terms of the first minimum in each radial distribution function, the numbers of water molecules in these shells were computed as 3.9, 6.1, and 1.3 for IV&, Nl&, and O&, respectively. The average residence times of water molecules in these shells were slightly less than 1 ps, but certain water molecules remain bound for periods as long as 30 ps. The major result of the simulations with 02 is the potential of mean force (pmf ) for 0, displacement in the SOD channel. It is not sufficient to use a single radial coordinate from Cu to construct this pmf, because a small variation in such a pmf could mask the occurrence of free energy barriers for non-radial displacements along the preferred pathways for 02 motion. The pmf is therefore represented as a function w(r) of the coordinates of one of the oxygen atoms in 0,. The pmf is shown in figs. 4 and 5 from two orthogonal directions. In fig. 4, the viewpoint is the
195
same as in fig. 2, with the channel running from the region of the Cu up to the left and away from the viewer. In fig. 5, the channel runs from the Cu up to the left but toward the viewer. As in fig. 2, a single instantaneous structure of the protein is included in the figures. A detailed examination of the pmf indicates that the motion of 0, from the surface of SOD to the vicinity of the Cu is not a simple “downhill” process in terms of the free energy change. One possible path involves the following steps. The 0, initially binds between the backbone rim segments formed by residues 56 through 59 and residues 138 through 141. This corresponds to the small low energy patch in the upper left of fgs. 4a and 5a. From there, the O,crosses a free energy barrier of about 1 kcal/mol to move into the broader, deeper minimum between the sidechains of Thr 135 and Arg 14 1. Finally, the O,crosses a barrier of about 2 to 3 kcal/mol to displace the water molecule bound to the Cu. The minimum free energies in the initial binding site and in the Cu binding site are, respectively, about 1 kcal/mol and 2 kcal/mol higher than in the primary minimum between the Thr 135 and Arg 14 1 sidechains. Although the pmf does not decrease monotonically in going from the mouth of the channel to the Cu, the variations in free energy are not large. This fact, together with the significant number of local minima indicated by figs. 4a and 5a, suggests that substrate molecules are likely to follow a variety of paths that differ in their details in moving through the SOD channel. The dynamics of these motions are also likely to be quite complicated, involving the rearrangement of hydrogen bonds among the substrate, water molecules and the channel residues. These rearrangements will be coupled to the fluctuations of the protein itself Indications of these complications are evident in very preliminary analyses of the average interactions of 02 with water molecules and protein atoms for different positions of Ot in the channel, and of trajectories in which the 02 is allowed to wander in the channel with no restraining potentials. The former analyses suggest that the initial binding to SOD may involve a direct interaction between 0, and the guanidinium moiety of Arg 14 1; the latter appears to move toward the mouth of the channel and to replace water molecules in the first hydration shell of 0,. In the broad minimum of the pmf be-
196
J. Shen, J.A. McCammon /Superoxide interacting with superoxide dismutase
(4
Fig. 4. Stereo views of the 05 potential of mean force (pmf ), viewed from the same direction as in fig. 2. The channel runs from the region of the Cu up to the left and away from the viewer. The contour levels are (a) 4 kcal/mol, and (b) 6 kcal/mol.
tween the Thr 135 and Arg 14 1 sidechains, the hydration number of O,- may decrease somewhat further, although the 02 may have bridging water molecules to one or both of the Thr 135 and Arg 14 1 sidechains. In one unrestrained dynamics trajectory, one can see an 02 associating with the Arg 141 sidechain and then moving to form a fairly long-lived ( > 5 ps) structure with bridging waters to one NH2 group of Arg 141 and a direct interaction with the OH of Thr 135. Much additional work must be done to analyze such details, however. Finally, it is of interest to consider the general implications of the present results in terms of the kinetics of SOD. Suppose that the steps in the SOD catalyzed reaction can be represented as SOD+06
kl eSOD//02
ka =SOD/O1-...
ing site, respectively. For simplicity, we assume for now that k3 is large enough to be neglected in estimating the overall rate. This is suggested by the fact that the dismutation of 01 by Cu2+ in water is diffusion controlled [ 25 1, although more analysis of the situation in the enzyme is required. Using Brownian dynamics to estimate k,, and using the volume sampled by 0, at the outer part of the channel to estimate k, /k_ ,, rough values for k, and k_ , are 6 x 1O9 M-’ s-’ and 1.6~ 10” s-‘, respectively. If the limiting step in moving through the channel is crossing a free energy barrier of 3 kcal/mol, the largest that k2 could be is 4.0 x 10’ 1s- I. The effective forward rate constant, k,k2 +kz
lir=k_,
b
’
(6)
,
k-1
(5) where // and / indicate structures in which the 02 is at the outer end of the channel and at the Cu bind-
would then be 0.7k, =4x lo9 M-’ s-‘, in the range of the experimental results, and the reaction would essentially be diffusion controlled. The channel crossing rate constant k2 is, however, very likely to be reduced somewhat by frictional effects and perhaps
J. Shen, J.A. McCammon /Superoxide interacting with superoxide dismutase
197
(4
M
Fig. 5. Stereo views ofthe 0~ pmffrom a direction orthogonal to that in f@. 4. The channel runs from the Cu up to the left but toward the viewer. The contour levels are (a) 4 kcal/mol, and (b) 6 kcal/mol.
by gating effects due to the motion of the protein [ 111. These issues are being analyzed in current work. Other issues that should be analyzed in the future include the effects of including electronic polarizability in the simulations, and the nature of the electron and proton-transfer steps at the end of the reaction.
Acknowledgement This work has been supported in part by NIH, NSF, the Robert A. Welch Foundation, the Texas Advanced Research Program, and the National Center for Super-computing Applications. JAM is the recipient of the George H. Hitchings Award from the Burroughs Wellcome Fund. We also thank Wilfred van Gunsteren for providing the GROMOS program.
References [ 1] I. Fridovich, in: Biochemical and Medical Aspects of Active Oxygen, eds. 0. Hayaishi and K. Asada (Univ. Park Press, Baltimore, 1977) pp. 3-12,171. [2] J.A. Fee, in: Metal Ion Activation of Dioxygen, ed. T.G. Spiro (Wiley, New York, 1980) p. 209. [3]R.A.GreenwaldandG.Cahen,eds,OxyRadicalsandtheir Scavenger Systems, Vol. II, Cellular and Medical Aspects (Elsevier, Amsterdam, 1983). [4]J.R. Stewart, W.H. Merill and W.H. Frist, in Oxygen RadicaB in Biology and Medicine, eds. M.G. Simic, K.A. Taylor, J.F. Ward and C. von Sonntag (Plenum Press, New York, 1988) p. 957. (5lE.M. Fielden, P.B. Roberts, R.C. Bray, D.J. Lowe, G.N. Mautner, G. Rotilio and L. Calabrese, Biochem. J. 139 (1974) 49: [6] A. Cudd and I. Fridovich, J. Biol. Chem. 257 (1982) 11443. [ 71 E. Argese, P. Viglino, G. Rotilio, M. Scarpa and A. Rigo, Biochemistry 26 (1987) 3224. [8] C. Bull and J.A. Fee, J. Am. Chem. Sot. 107 (1985) 3295. [9] J.A. FeeandC.Bull, J.Biol.Chem. 261 (1986) 13QOO. [ IO] L. Banci, I. Bertini, C. Luchinat and R.A. Hallewell, J. Am. Chem.Soc. 110 (1988) 3629. [ 111 J.A. McCammon and S.C. Harvey, Dynamics of Proteins and Nucleic Acids (Cambridge Univ. Press, Cambridge, 1987).
198
J. Shen, J.A. McCammon /Superoxide interacting with superoxide dismutase
[ 12 ] S.A. Allison and J.A. McCammon, J. Phys. Chem. 89 (1985) 1072. [ 131 K. Sharp, R. Fine and B. Honig, Science 236 (1987) 1460. [ 14) T. HeadGordonand C.L. Brooks, J. Phys. them. 91 (1987) 3342. [ 15 ] S.A. Allison, R.J. Banquet and J.A. McCammon, Biopolymers 27 (1988) 251. [ 161 J.J. Sines, S.A. Allison and J.A. McCammon, Biochemistry 29 ( 1990) 9403. [ 171 J. Shen, S. Subramaniam, C.F. Wongand J.A. McCammon, Biopolymers 28 ( 1989) 2085. [ 181 J.A. Tainer, E.D. Getzoff, KM. Beem, J.S. Richardson and D.C. Richardson, J. Mol. Biol. 160 (1982) 181.
[ 191 H.J.C. Beret&en and W.F. van Gunsteren, personal communication. [20] J. Shen, S. Subramaniam, C.F. Won& T. Albright and J.A. McCammon, J. Comput. Chem. 11 ( 1990) 346. [ 2 1 ] H.J.C. Berendsen, J.R. Gigera and T.P. Straatsma, J. Phys. Chem. 91 (1987) 6269. [ 22 ] J. Shen, C.F. Wong and J.A. McCammon, J. Comput. Chem. II (1990) 1003. [23] J. Shen, Ph.D. Thesis (University of Houston, 1990). [24] B.P.Gaber,R.D.Brown,S.H.KoenigandJ.A.Fee,Biochim. Biophys. Acta 271 (1972) 1. [ 251 J. Rabani, D. Klug-Roth and J. Lilic, J. Phys. Chem. 77 (1973) 1169.