PROTEINS: Structure, Function, and Genetics 11:142-152 (1991)
Molecular Dynamics Simulations of the Cytochrome c,-Rubredoxin Complex From Desulfovibrio vulgaris David E. Stewart' and John E. Wamplel.2 'University Computing and Networking Services, Specialized Systems Support, Computer Services Annex, and 'Department of Biochemistry, Boyd Graduate Studies Research Center, The University of Georgia, Athens, Georgia 30602
ABSTRACT Molecular dynamics simulations have been carried out on the complex formed between the tetraheme cytochrome c3 and the iron protein rubredoxin from the sulfate-reducing bacterium Desulfovibrio vulguris. These simulations were performed both with explicit solvent water molecules included, and without solvent molecules using a distancedependent dielectric constant to approximate the screening effects of solvent. The results of both simulations are strikingly different, indicating that the representation of environmental effects is important in such simulations. F o r example, a striking adaptation of the two proteins seen in the nonsolvated simulation is not seen when explicit solvent water is included; in fact, the complex appears to become weaker in the solvated simulation. Nonetheless, the iron-iron distance decreases more significantly in the solvated simulation than in the nonsolvated simulation. It was found that in both cases molecular dynamics optimized the structures further than energy minimization alone. Key words: AMBER, electrostatics, proteinprotein interactions, intermolecular bonding, electron transfer, redox proteins INTRODUCTION The importance of biological electron transfer in the energy metabolism of living systems has made it the subject of intense research. One aspect critical for a full understanding of this process is the manner in which electron transfer proteins interact. Based on spectroscopic evidence and the analysis of available X-ray structures we have proposed models for electron transfer complexes of small proteins from the sulfate-reducing bacterium Desulfovibrio vulgaris. The supporting evidence for these models has been discussed In this report one of these models has been refined using molecular dynamics (MD) simulations. Cytochrome c3 from Desulfovibrio is unique among the c-type cytochromes in that there are four hemes per polypeptide chain of approximately 110 amino acid residues. The hemes exhibit distinct re0
1991 WILEY-LISS, INC.
dox properties with negative redox potentials as low as -340 mV.5 In the two known crystal structures627 the hemes are well separated (average Fe-to-Fe distance of 14.5 A) and are arranged with a particularly obvious lack of planar alignment (average heme-to-heme angle of 72"). Rubredoxins from these same organisms are much smaller, with between 50 and 55 residues, and contain one nonheme iron atom ligated by four cysteine sulfhydryl groups. The physiological function of rubredoxin has not yet been determined,' but it has been shown t o act as an effective redox partner for cytochrome c3 in ~ i t r oRu.~ bredoxins exhibit redox potentials in the range of 0 to -50 mV.l0 The physiological function of cytochrome c3 is believed to be as an electron carrier for the enzyme hydrogenase.' Rubredoxin can be rapidly reduced in vitro by cytochrome c3 in the presence of hydroge n a ~ e Thus, .~ while no physiological significance can be claimed for the proposed complex, such models can serve to aid investigation of electron transfer. The surface of rubredoxin in the region of its redox active iron atom is uniformly negatively ~ h a r g e d . ~ , ~ Although each of the four hemes of cytochrome c3 is exposed t o the surface, only one (heme 1) is surrounded by a complementary uniformly positive charged It has been shown that many redox proteins interact via electrostatic interact i o n ~ . ' ' - ~On ~ the basis of the pattern of apparent electrostatic complementarity, along with steric and geometric considerations, this heme has been proposed to be the site a t which rubredoxin and flavodoxin interact.lP3 NMR data indicate that 1:l complexes are formed between the interacting proteins and that spectral lines associated with the same heme group are perturbed in both cases.2 Computational studies have been undertaken pre-
Received October 1, 1990; revision accepted February 20, 1991. Address reprint requests to David E. Stewart, The University of Georgia, University Computing and Networking Services, Specialized Services Support, Computer Services Annex, Athens, GA 30602.
143
DYNAMICS OF THE RUBREDOXIN-CYTOCHROME c3 COMPLEX 1
"
"
"
"
'
I
-05-150 -
-"""I m
x
m
-3500
L
al
-
c w c
L W
-3700
0
.d
3
c-'
rn
c
-166001
0 0 L 3
B
u 0
a
-150
W c,
a,
c
-16800
-250 -350
-17000 -17200
-5501 ,
-17400 10
40
70
100
,
,
10
viously on the interactions of redox protein^.'^-^^ A molecular dynamics study similar to that reported here was carried out by Wendoloski et al.15 on the interaction between cytochromes c and b,. Again, while this interaction has no physiological significance-the pairing modeled involved proteins from different organisms-an in vitro interaction had previously been demonstrated.l8 In the cytochrome c-cytochrome b, study it was observed that the total potential energy and the heme-to-heme distance decreased during the simulation. This observation indicated that, in addition to providing details about movements of atoms, the dynamics simulation optimized the structure of the complex more fully than energy minimization alone. Molecular dynamics has been shown to be effective in optimizing structures to a greater extent than energy minimization in other cases as well.19-22The advantage of molecular dynamics is its inherent ability to explore a larger region of the conformational space available to the molecule or complex than energy minimization since the kinetic energy of the atoms in MD simulations allows local energy barriers to be crossed.23 In the case of the cytochrome c-cytochrome b, model,15 a striking change occurred in the interface between the two proteins as the simulation progressed. The heme groups turned from a edge-on configuration and, along with a phenylalanine resi-
,
,
40
Time (ps) Fig. 1. Change in total AMBER potential energy of the complex. (A) During the nonsolvated simulation; ( 8 ) during the solvated simulation.
,
, 70
,
,
'I1 100
Time (ps) Fig. 2. Nonbonded interaction energies. (A) During the nonsolvated simulation; (B) during the solvated simulation.
due of cytochrome c, a stack of aromatic rings was formed. MD simulations may be performed using a distance-dependent dielectric constant to approximate the screening of electrostatic interactions by solvent, or with explicit water molecules included to represent the aqueous environment of the protein. The merits of each of these approaches have been discussed previously. One of the implications of the cytochrome c-cytochrome b, sir nu la ti on^'^ was that the two approaches can lead to similar results in the interface region of a protein-protein complex. In this study, we compare the results obtained from MD simulations using a distance-dependent dielectric constant to those in which explicit water molecules were included, focusing particularly on the interface of the complex. 15724--26
METHODS The program AMBER 3.0 UCSFZ7was used for the molecular mechanics and dynamics calculations and for the analysis of the results. Minimization and MD calculations were performed on the University of Georgia Cyber 205 supercomputer, and analysis and display of the results were performed on a Silicon Graphics Iris workstation. The version of AMBER 3.0 used was optimized to run on the Cyber 205. It has been found that a 100-fold increase in computational speed can be realized by vectorization
144
D.E. STEWART AND J.E. WAMPLER
-41ooLl
-1
A
-4300 -4500 x Is)
L
-4700
W
c
c
W
c 0 -+
0
-4900
c, 0
m
L
-1600
rn
c,
c
H
a
-1800 -2000
2.1
1.7
.gr . i I
-2200 -2400
.5
10
40
70
100
10
Fig. 3. Water-protein interaction energies during the solvated simulation. (A) rubredoxin; (B) cytochrome c,.
0
=I
191 L
70
100
Time (ps)
T i m e (ps)
1,
40
1
Structure # Fig. 4. (A) Distance (in A) between the centers of mass of the two proteins during solvated simulation (solid line, square symbols) and during the nonsolvated simulation (dashed line, circles). (B) Radius of gyration (in A) for the complex during the solvated simulation (solid line) and the nonsolvated simulation (dashed line).
Fig. 5. Root-mean-squared (rms) deviation for all atoms from the initial (0.1 psec) structure. (A) During the nonsolvated simulation; (B) during the solvated simulation.
of the program on such high performance computer~.~,~' The starting structure has been described previously.2 The coordinates of the proteins were obtained from the Brookhaven Protein Data Bank.29,30 The structures used in the complex were rubredoxin from D. vulgaris3' and cytochrome c3 from D. vulgaris strain M i y a ~ a k i . ~ ~ In the nonsolvated simulation the energy of the model-built complex was minimized until its rootmean-square (rms) gradient of the energy was less than 0.025 kcal/(mol A). A distance-dependent dielectric was used with a nonbonded cutoff of 9.5 A. Initially 125 cycles of constrained energy minimization were performed to gradually relax steric conflicts without introducing large distortions in the structure. This was followed by unconstrained conjugate gradient minimization until the limit of convergence was reached. A united-atom representation of the hydrogen atoms was used, with only the hydrogen atoms on polar atoms explicitly included,33 resulting in a total of 1,661 atoms (269 hydrogen 1,392 heavy atoms) for the complex. The energy-minimized structure was used as the starting coordinates for the MD run. Initial velocities assigned to the minimized structure correspond to a mean temperature of 50 K. The heating phase consisted of 30 K temperature increments applied every 200 fsec to a final temperature of 300 K. Equilib-
+
DYNAMICS OF THE RUBREDOXIN-CYTOCHROME c3 COMPLEX
145
A
5 ps.
30 ps.
B
5 ps.
30 ps.
Fig. 6. The adaptive changes in the two proteins over the course of the nonsolvated simulation are illustrated by two structures from the simulation shown with orthogonal, stereo views. In both cases the upper structures is that at 5 psec and the lower structure is from 30 psec. (A) View along binding cleft showing
symmetry of the rubredoxin (on left) with the iron center represented by a dot for the iron atom. (B) Orthogonal view to A showing how the beta sheets at the top of the two proteins are drawn together by formation of ionic and hydrogen bonds as listed in Table IA.
rium at 300 K was ensured by repeatedly reassigning random velocities from the appropriate Maxwellian distribution function every 200 fsec for a period of 2.0 psec. The system was then allowed to stabilize for a period of 4.0 psec after which data collection was begun for 65 psec. In the solvated simulation a 6 A shell of 1001 TIP3P water molecules was added to the model-built complex.34 A constant dielectric of E = 1 and nonbonded cutoff of 9.5 A was used. This system was
minimized for 3,500 cycles. The minimized structure was heated from 0 K to 300 K in 30 K increments applied every 200 fsec. The system was equilibrated at 300 K for 2.0 psec, followed by a stabilization period of 10 psec. Data collection continued for 100 psec. In both simulations a 0.001 psec time step was used. Bond lengths were fixed at their optimal values using the SHAKE algorithm.35 Structures were collected every 0.1 psec (100 steps). Trajectories and parameters for analysis were ob-
146
D.E. STEWART AND J.E. WAMPLER
A
100 ps.
B
Fig. 7. Orthogonal, stereo views of the initial (1 psec) and final (100 psec) structures from the solvated simulation. (A) View along binding cleft showing symmetry of the rubredoxin (on left) with the iron center represented by a dot for the iron atom. (B) Orthogonal view to A. In both cases little overall change is seen and this result is characteristic of the entire simulation.
tained using the AMBER analysis module. Most manual manipulations, display, and plotting were done using the program SYBYL (Tripos Associates, Inc., St. Louis, MO). Animation of the dynamics simulations was carried out using QUANTA (Polygen Corporation, Waltham, MA). Plotting and manipulations of the data sets representing the trajectories were done using SPECOS SA.36
RESULTS The General Character of the Complex The changes in the potential energy during the simulations of the complexes are shown i n Figure 1. The conformations of the complexes from the MD simulations were energetically more favorable than those from energy minimization. The potential energy of the complex rapidly becomes more favorable
DYNAMICS OF THE RUBREDOXIN-CYTOCHROME c3 COMPLEX
147
TABLE I. Status of Interface Bonds During the MD Simulation Intermolecular bond ( c , -+ rd) residue (atom*) K 57 (NZ) - E 17 (OE) K 58 (NZ) - E 17 (OE) K 58 (NZ) - D 19 (OD) K 95 (NZ) - E 50 (OE) KlOl (NZ) - E 12 (OE) K 15 (NZ) - E 17 (OE) K 94 (NZ) - E 50 (OE) D 71 (OD) - C 42 (N) K 57 (NZ) - E 12 (0) K 60 (NZ) - P 40 (0) K 72 (NZ) - V 41 (0) K 94 (NZ) - T 7 (0) K 94 (NZ) - V 8 (0) ClOO (SG) - E 12 (N) ClOO (SG) - Y 11 (N) K 58 (NZ) - G 18 (0) K 58 (NZ) - Y 11 (OH)
Y 66 (OH) - Y 11 (N) K 52 (NZ) - V 8 (0)
Present in Present min’d struct. at end Comments IA. During the Nonsolvated Simulation X X Maintained throughout X X Maintained throughout X X Maintained throughout X X Maintained throughout X X Maintained throughout X Formed between 30 and 35 psec X Formed between 25 and 30 psec X Formed during stabilization period X Formed between 25 and 30 psec X Broken between 0 and 5 psec Broken during stabilization period X Broken between 25 and 30 psec X X Broken between 0 and 5 psec Broken between 20 and 25 psec X Broken between 20 and 25 psec X Formed during stabilization period X X Formed during equilibration period Broken during stabilization period Reformed between 10 and 15 psec X Formed between 20 and 25 psec X Formed between 0 and 5 psec IB. During the Solvated Simulation
K 57 (NZ) - E 17 (OE) K 58 (NZ) - D 19 (OD) K 95 (NZ) - E 50 (OE) KlOl (NZ) - E 12 (OE) K 58 (NZ) - T 21 (OG) K 60 (NZ) - P 40 (0) K 72 (NZ) - V 41 (0) K 94 (NZ) - T 7 (0) K 94 (NZ) - V 8 (0)
X X X X X X X X X
X
Broken during equilibration Maintained throughout Broken during equilibration Broken during equilibration Broken during first 5 psec of stabilization Broken during equilibration Broken during first 5 psec of stabilization Broken during equilibration Broken between 30 and 35 psec
*Residues identified by standard one letter code. Number represents position in protein sequence (from N-terminal). Atom types are as follows: NZ, epsilon amino of lysine; OE and OD, terminal carboxyl group oxygen of glutamic and aspartic acid residues; N, main chain amide hydrogen; 0, main chain amide carbonyl; SG, terminal sulfhydryl of cysteine.
in approximately the first 35 psec of the data collection period in the nonsolvated simulation, while such a change is not seen in the solvated complex until after approximately 50 psec. It is apparent from the figure that the potential energy is much more negative in the solvated simulation than in the nonsolvated simulation. This arises from the inclusion of the water molecules and to a larger extent from the value of the dielectric constant used: in the nonsolvated simulation the charges within the protein are shielded due to the distance-dependent dielectric; in the solvated simulation the dielectric constant is 1 and no such explicit shielding occurs. The strength of the complex may be evaluated by calculating the nonbonded interaction energies (electrostatic, van der Waals, and hydrogen bonding) between rubredoxin and cytochrome cg. These are shown in Figure 2. The hydrogen bonding and van der Waals interaction energies are similar in both simulations; however, there is a significant difference in the electrostatic interaction energies. This
again arises mainly from the dielectric constant used. The total interaction energy in the nonsolvated simulation generally becomes more favorable during the simulation, particularly in the 30-35 psec time frame. Conversely, in the solvated simulation, no such trend is observed, instead, the rubredoxin-cytochrome cg interaction energy becomes less favorable until the latter portion of the simulation (90100 psec). The driving force for the unfavorable change is apparently the interactions between the solvent molecules and the protein (Fig. 3). The total interaction energy between the three “groups” of molecules-rubredoxin, cytochrome c3, and waterbecomes more favorable during the simulation due to the contributions from the solvent-protein interactions. The specific interactions (see below) seem to indicate that the two proteins are less directly involved as the solvated simulation progresses. Two gross measures of this are the distance between the centers of masses of the two proteins (Fig. 4A) and the
148
D.E. STEWART AND J.E. WAMPLER
radius of gyration of the complex (Fig. 4B). In both simulations the distances between the centers of mass change comparably through approximately 30 psec at which time this distance begins to increase steadily in the solvated simulation, indicating that the two proteins are separating. The complex radius of gyration decreases initially only in the nonsolvated simulation, probably due to the strong electrostatic interactions arising from the use of a distance-dependent dielectric constant. No such change is seen in the solvated simulation. In Figure 5, the rms deviations of the atoms during the simulations are shown. In both simulations, a maximum rms deviation of slightly less than 3 A is observed, however, the time required for the simulations to reach this value is different. The rms deviation increases most rapidly in the nonsolvated simulation, due to the absence of water molecules to damp the motions of the atoms. In this case, the rms deviation appears to stabilize after approximately 35 psec suggesting formation of a stable conformation. The final structures from each simulation are significantly different, with an rms deviation of 4.4 A between the 65 psec nonsolvated structure and the 100 psec solvated structure for all atoms; this decreases to 4.2 A if the hydrogen atoms are omitted from the calculation. The rms deviation between rubredoxin in the two final structures is 3.4 A; for cytochrome c3 this is 3.5 A (heavy atoms only). The lower rms deviation for the individual proteins as compared to the complexes suggests that there is a significant difference in the orientation of the molecules in the complex between the solvated and nonsolvated structures. In a qualitative sense the overall interpretation of these results is illustrated with two structures from both simulations as shown in Figures 6 and 7. Under the drive of the electrostatic effect, the two proteins change shape during the nonsolvated simulation (Fig. 6). The net result is that both the binding depression on the cytochrome and the shape of the rubredoxin become more symmetrical. On the large scale most of these changes represent “sliding” motions of beta sheet portions of both proteins. During the solvated simulation (Fig. 7), these changes do not take place and the shapes of both proteins stay more or less constant. The intercalation of water and loss of electrostatic bonding (see below) prevent these changes from taking place over the time period of the simulation.
Interactions in the Interface In the nonsolvated simulation, the minimized structure has 5 ionic and 5 hydrogen bonds, while the 0.1 psec data collection structure has 5 ionic and
I (u
U
13 12 11 10
10
40
70
100
Time (ps) Fig. 8. Iron-iron distances between the iron center of rubredoxin and the iron of heme 1 of cytochrome c, (distance in A). (A) During the nonsolvated simulation; (B) during the solvated simulation.
tion as discussed above. It has been suggested that an inordinate number of such bonds may be formed due to the lack of inclusion of explicit ~ o l v e n t . ~ ~ , ~ ~ This criticism is confirmed by the number of intermolecular ionic and hydrogen bonds found for the solvated structures: there are 4 ionic and 6 hydrogen bonds in the minimized structure, 2 ionic and 2 hydrogen bonds in the 0.1 psec data collection structure, and 1 ionic bond in the 100 psec structure. These changes are caused by the intercalation of water molecules into the interface between the two molecules disrupting the intermolecular bonds. The fate of these bonds during the simulations is summarized in Table I. Electron transport is dependent on the orientation and distance between redox group^.^^-*^ One criterion to judge the optimization of the structures resulting from these simulations, then, would be the distance between prosthetic groups as a function of time, with the later structures expected to have a smaller iron-to-heme distance. In the nonsolvated simulation, the iron-to-heme distance decreased slightly (Fig. 8). A more significant decrease was seen in the solvated simulation (Fig. €9, with the distance being much smaller a t the end of this simulation than in the nonsolvated case.
7 hydrogen bonds, and the 65 psec structure has 7
A striking rearrangement is seen in the organiza-
ionic and 6 hydrogen bonds. This is reflected in the increasing electrostatic energy during the simula-
tion of residues in the interface between the two proteins in the nonsolvated simulation. This is
DYNAMICS OF THE RUBREDOXIN-CYTOCHROME
149
COMPLEX
f
5 ps.
Y
30 ps.
Y
Y
a Y
(
Fig. 9. Interactions in the vicinity of heme 1 of cytochrome c, at the interface between the two proteins. Residues from rubredoxin are shown with the dithered line. The rearrangement of ionic and hydrogen bonds summarized in Table IA is accompanied by changes in mobility and orientation of the aromatic residues (see
Fig. 10). Residues shown for cytochrome c, are heme 1, histidine 70 and 106, tyrosine 66, lysines 15,57, and 58,and aspartate 56. Residues of rubredoxin shown are tyrosine 1 I , glutamate 17, and aspartate 19.
shown in the structures of Figure 9. To examine the fluctuations of the sidechains of the aromatic residues, the value of the dihedral angle x2 (C,-C,C,-C,1) was calculated for each of the 0.1 psec data collection structures. The trajectories of dihedral angle rotation of the residues involved (Fig. 10) showed a concerted change at 30 psec and obvious coupling between the aromatics as revealed by the changes in amplitude of their variations. This change is driven by a rearrangement of hydrogen and ionic bonding partners which occurs in a concerted fashion in a few picoseconds. The change is completed by 30 psec. The general conformation and configuration of the two molecules before this change is illustrated by the structure a t 5 psec shown in Figure 9 (top). Similarly, after the transition has occurred, the structure (illustrated by the 30 psec structure in Fig. 9, bottom) remains fairly constant for the remaining 35 psec of the simulation. It is interesting to note that while this transition occurs under the inf hence of a rearrangement of electrostatic interactions, that the major changes involved are interactions be-
tween residues within each of the proteins rather than between them. Thus, the ionic bonds of the heme propionate residue with both the main chain amide and €-amino of the nearby lysine (number 15) are changed into a n ionic bridge involving lysine 15, the heme propionate, and lysine 57. Similarly, tyrosine 11or rubredoxin flips over so that it bonds to the opposite oxygen of the carboxyl of glutamate 17 leaving the lower oxygen (as seen in the figure) free to also bind to lysine 57 of the cytochrome. The trajectories of the aromatic residues during the solvated simulation (Fig. 11)reveal no such coupled changes and the ionic and hydrogen bonds between the proteins are, for the most part, broken (Table I).
DISCUSSION While it has been shown that the exclusion of solvent molecules in MD calculations may result in the formation of improper hydrogen bond^'^,^^ and that the atomic motions exhibit artificially large amplitUdes,15,26,43 a widely held opinion is that simulations with a distance-dependent dielectric approxi-
150
D.E. STEWART AND J.E. WAMPLER
180
I
L TTyr y r 11 120 (ru) 60 I
i i
O?
k
-60 - 6 0 r TTyr y r 66 I
-120
t-I
(c,) (c3) ' Y '
-180-'
'
I
'
I
10
0
'
I
I
30
20
'
40
I
I
50
60
-
Time (ps) Fig. 10. Trajectories of dihedral angle xz changes in interface during the nonsolvated simulation. Angles for tyrosine 11 are shown offset ( + 90").
-180k'
'
I
'
I
0
I
10
'
I
20
I
30
I
40
'
50
I
I
60
I
70
I
80
I ' - I 90 100
Tine (ps) Fig. 11. Trajectories of dihedral angle x2 of same residues shown in Figure 10 during the solvated simulation. Angles for tyrosine 66 are shown offset (-330").
mation may still be instrumental in understanding the dynamic behavior of proteins. In the previous study of the cytochrome c-cytochrome b, complex,15 fully solvated and nonsolvated simulations lead to
very similar changes within the interface between the proteins, but on a much slower time scale for the solvated simulation. In addition, since the atomic motions are not damped by the solvent, a larger por-
DYNAMICS OF THE RUBREDOXIN-CYTOCHROME c3 COMPLEX
tion of conformational space may be searched by such sir nu la ti on^.'^.^^ Our results are less clearly in support of this idea that nonsolvated and solvated simulations lead to similar changes only on a different time scale. In both simulations, the resulting structure had a significantly lower potential energy than did the minimized structure. However, the interaction energies between the two proteins exhibited opposite trends. In the nonsolvated simulation, the interaction energy became more favorable, while in the solvated simulation the converse was true until the 90 psec time period. This unfavorable change in the interaction energy appears to be driven by favorable interactions between the proteins and the water molecules. The obvious explanation for this difference is intercalation of water molecules into the interface between the two proteins, solvating the ionic bonds. The consequences of these changes are summarized in Table I. The resulting difference is clear in the stereo plots of Figures 6 and 7. Without solvent, the proteins adapt to each other in a n intriguing way forming a fairly symmetrical “lock-and-key’’ type of complementary fit. With solvent they remain independent entities. The suggestion that the water molecules intercalate into the interface of the complex is supported not only by the breaking of the intermolecular bonds (and associated decrease in electrostatic interaction energy) but also by the change in the radius of gyration (R,). The R, of the nonsolvated complex decreases quickly, as the number of intermolecular bonds increases, while the R, of the solvated complex remains stable, increasing slightly in the latter half of the simulation. Overall, both the R, and the distance between centers of mass for the solvated complex increase. Extrapolating these two trajectories to longer time suggests that the two molecules are drifting apart although the iron-iron distance is actually decreasing. The comparison between early and late structures for both simulations is shown in Figures 6 and 7. The changes seen are in essentially opposite directions. Without solvent, rubredoxin flattens and adapts to the surface of the cytochrome. In the view shown, this involves a counterclockwise rotation of rubredoxin relative to cytochrome c,. With solvent, rubredoxin becomes slightly more spherical and appears t o rotate slightly clockwise in this view. Molecular dynamics simulations of the proposed cytochrome c,-rubredoxin complex reveal that the complex is stable on the picosecond time scale under the influences of the AMBER force field. Spectroscopic evidence’ supports both the specificity and the close hemeinonheme iron interactions modeled. However, the manually docked proteins from our previous work1-, are obviously not in a n optimum configuration since both simulations lead to lower
151
energies. The two simulations lead to different results. In the nonsolvated case the proteins “roll” together and change shape in an adaptive way. In the solvated simulation they “roll” in the opposite direction on trajectories which may eventually lead to disruption of the complex. It is possible that the nonsolvated simulation probes areas of conformational space that are not accessible during the 100 psec time period of the solvated simulation. The interesting changes in the interface between the two proteins might be similarly accessible. One approach to investigation of this question would be to solvate the terminal structure from the nonsolvated simulation and see if the aromatic stacking and electrostatic interaction (ionic and hydrogen bonds) are disrupted in this case. Within the limits of available computer time, this study and similar simulations of the cytochrome c,-flavodoxin complex2 are currently being pursued. A wealth of data comes from these simulations. Every residue in both proteins undergoes changes in relative orientation and in its interactions with others. A complete analysis of the changes in the aromatic residues within these structures (manuscript in preparation) reveals that coupling and stacking between residues like those shown in Figure 9 are common.
ACKNOWLEDGMENTS The authors thank Dr. Walter McRae (University Computing and Networking Services) and the University of Georgia Research Foundation for their support of this work. Vectorization of major portions of the AMBER suite of programs for operation on the Cyber 205 was performed by Dr. Steve Gallion. This research was supported in part by University Computing and Networking Services, by grants from the University of Georgia Research Foundation, and the National Institutes of Health (Grant GM 41482).
REFERENCES 1 Stewart, D.E., LeGall, J., Moura, L., Moura, J., Peck, H.D., Jr., Xavier, A.V., Weiner, P.K., Wampler, J.E. A hypothetical model of the flavodoxin-tetraheme cvtochrome c , complex of sulfate-reducing bacteria. Biochkmistry 27:54442450. 1987. 2. Stewak,D.E., LeGall, J., Moura, I., Moura, J., Peck, H.D., Jr., Xavier, A.V., Weiner, P.K., Wampler, J.E. Molecular modeling and NMR studies of the rubredoxin-cytochrome cg complex of sulfate-reducing bacteria. Eur. J . Biochem 185:695-700, 1989. 3. Stewart, D.E. The structure, interactions, and dynamics of electron transport proteins from Desulfouibrzo. A molecular modeling and computational study. Ph.D. Dissertation, The University of Georgia, 1989. 4. Wampler, J.E., Stewart, D.E., Gallion, S.L. Molecular dynamics simulations of proteins and protein-protein complexes. In: “Computer Simulation Studies in Condensed Matter Physics 11.” Landau, D.P., Mon, K.K., Schuttler, H.G. (eds.). Berlin: Springer-Verlag, 1990: 68-84. 5. Yagi, T., Inokuchi, H., Kimura, K. Cytochrome cj, a tetrahemoprotein electron carrier found in sulfate-reducing bacteria. Acc. Chem. Res. 16:2-7, 1983. 6. Higuchi, Y., Kusunoki, M., Matsuura, Y., Yasuoka, N.,
152
7.
8. 9.
10. 11. 12.
13. 14.
15. 16.
17. 18. 19.
20. 21. 22. 23. 24. 25. 26.
D.E. STEWART AND J.E. WAMPLER
Kakudo, M. Refined structure of cytochrome c3 at 1.8 A. J . Mol. Biol. 172:109-139, 1984. Haser, R., Pierrot, M., Frey, M., Payan, F., Astier, J.P., Bruschi, M., LeGall, J . Structure and sequence of the multiheme cytochrome c3. Nature (London) 282:806-810, 1979. LeGall, J., Fauque, G. Dissimilatory reduction of sulfur compounds. In: “Biology of Anaerobic Microorganisms.” Zehnder, A.J.B. (ed.). New York: Wiley, 1988: 587-639. Bell, G.R., Lee, J.P., Peck, H.D., Jr., LeGall, J . Reactivity of Desulfovibrio gigas hydrogenase toward artificial and natural electron donors or acceptors. Biochimie 60:315320,1978. Moura, I., Moura, J.J.G., Santos, M.H., Xavier, A.V., LeGall, J . Redox studies on rubredoxins from sulfate and sulfur reducing bacteria. FEBS Lett. 107:419-421, 1979. Salemme, F.R. A hypothetical structure for a n intermolecular electron transfer complex of cytochromes c and b,. J . Mol. Biol. 102:563-568, 1976. Smith, H.T., Staudenmayer, N., Millet, F. Use of specific lysine modifications to locate the reaction site of cytochrome c with cytochrome c oxidase. Biochemistry 16: 4971-4974, 1977. Mauk, M.R., Reid, L.S., Mauk, A.G. Spectrophotometric analysis of the interaction between cytochrome b, and cytochrome c. Biochemistry 21:1843-1846, 1982. Erman, J.E., Vitello, L.B. The binding of cytochrome c peroxidase and ferricytochrome c. A spectrophotometric determination of the equilibrium association constant as a function of ionic strength. J . Biol. Chem. 255:6224-6227, 1980. Wendoloski, J.J., Matthew, J.B., Weber, P.C., Salemme, F.R. Molecular dynamics of a cytochrome c-cytochrome b, electron transfer complex. Science 238794-797, 1987. Northrup, S.H., Boles, J.O., Reynolds, J.C.L. Electrostatic effects in the Brownian dynamics of association and orientation of heme proteins. J . Phys. Chem. 915991-5998, 1987. Northrup, S.H., Boles, J.O., Reynolds, J.C.L. Brownian dynamics of cytochrome c and cytochrome c peroxidase association. Science 241:67-70, 1988. Mauk, M.R., Reid, L.S., Mauk, A.G. Spectrophotometric analysis of the interaction between cytochrome b, and cytochrome c. Biochemistry 21:1843-1846, 1982. Brooks, B.R., Bruccoleri, R.E., Olafson, B.D., States, D.J., Swaminathan, S., Karplus, M. CHARMM: A program for macromolecular energy, minimization, and dynamics calculations. J . Comp. Chem. 4:187-217, 1983. Levitt, M. Protein folding by restrained energy minimization and molecular dynamics. J . Mol. Biol. 170:723-764, 1983. Elber, R., Karplus, M. Multiple conformational states of proteins: A molecular dynamics analysis of myoglobin. Science 235318-321, 1987. Brunger, A.T., Kuriyan, J., Karplus, M. Crystallographic R factor refinement by molecular dynamics. Science 235: 458-460, 1987. van Gunsteren, W.F. The role of computer simulation techniques in protein engineering. Protein Eng. 25-13, 1988. Brooks, C.L., Karplus, M. Solvent effects on protein motion and protein effects on solvent motion. J . Mol. Biol. 208:159-181, 1989. Wendoloski, J.J., Matthew, J.B. Molecular dynamics effects on protein electrostatics. Proteins 5313-321, 1989. Levitt, M., Sharon, R. Accurate simulation of protein dy-
27.
28.
29.
30.
31.
32. 33.
34.
35. 36. 37.
38.
39.
40. 41. 42. 43.
namics in solution. Proc. Natl. Acad. Sci. U.S.A. 8575577561,1988. Singh, U.C., Weiner, P.K., Caldwell, J.W., Kollman, P.A. AMBER (UCSF version 3.0). Department of Pharmaceutical Chemistry, University of California, San Francisco, 1986. Gallion, S.L., Levy, R.M., Weiner, P.K., Hirata, F. Implementation of a macromolecular mechanics program on a CYBER 205 supercomputer. Comp. Chem. 10:165-173, 1986. Bernstein, F.C., Koetzle, T.F., Williams, G.J.B., Meyer, E.F., Jr., Brice, M.D., Rodgers, J.R., Kennard, O., Shimanouchi, T., Tasumi, M. The Protein Data Bank: A computer-based archival file for macromolecular structures. J . Mol. Biol. 112535-542, 1977. Abola, E.E., Bernstein, F.C., Bryant, S.H., Koetzle, T.F., Weng, J. Protein Data Bank. In: “Crystallographic Databases-Information Content, Software Systems, Scientific Applications.” Allen, F.H., Bergerhoff, G., Sievers, R. (eds.). Bonn: Data Commission of the International Union of Crystallography, 1987: 107-132. Adman, E.T., Sieker, L.C., Jensen, L.H., Bruschi, M. LeGall, J . A structural model of rubredoxin from Desulfovibrio vulgaris at 2 Angstrom resolution. J. Mol. Biol. 112:113-120, 1977. Higuchi, Y., Kusunoki, M., Matsuura, Y., Yasuoka, N., Kakudo, M. Refined structure of cytochrome c3 at 1.8 Angstrom resolution. J. Mol. Bio. 172:109-139, 1984. Weiner, S.J., Kollman, P.A., Case, D.A., Singh, U.C., Ghio, C., Alagona, G., Profeta, S., Weiner, P. A new force field for molecular mechanical simulation of nucleic acids and proteins. J. Am. Chem. SOC.106:765-784, 1984. Jorgensen, W.L., Chandrasekhar, J., Madura, J., Impey, R.W., Klein, M.L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79:926935, 1983. van Gunsteren, W.F., Berendsen, H.J.C. Algorithms for macromolecular dynamics and constraint dynamics. Mol. Phys. 341311-1327, 1977. Wampler, J.E. SPECOS SA (The Stand Alone Version) User’s Guide, Instrument Software Services, Arnoldsville, GA, 1989. Levitt, M., Sharon, R. Simulating protein dynamics in solution: Bovine pancreatic trypsin inhibitor. In: “Crystallography in Molecular Biology.” Moras, D., Drenth, J., Strandberg, B., Suck, D., Wilson, K. (eds.). New York: Plenum, 1985 197-205. Mauk, A.G., Scott, R.A., Grav. H.B. Distance of electron transfer to and from metalliprotein redox sites in reactions with inorganic complexes. J . Am. Chem. Soc. 102: 4360-4363, 1980. Makinen, M.W., Schichman, S.A., Hill, S.C., Gray, H.B. Heme-heme orientation and electron transfer kinetic behavior of multisite oxidation-reduction enzymes. Science 222:929-931,1983. Siders, P., Cave, R.J., Marcus, R.A. A model for orientation effects in electron-transfer reactions. J. Chem. Phys. 81: 5613-5624,1984. Scott, R.A., Mauk, A.G., Gray, H.B. Experimental approaches to studying biological electron transfer. J. Chem. Ed. 62:932-938, 1985. Cave, R.J., Siders, P. Marcus, R.A. Mutual orientation effects on electron transfer between porphyrins. J . Phys. Chem. 90:1436-1444,1986. Chen, L.X.-Q., Engh, R.A. Brunger, A.T., Nguyen, D.T., Karplus, M., Fleming, G.R. Dynamics studies of apoazurin of Alcaligenes denitrificans. Biochemistry 27:6908-6921, 1988.
Molecular Dynamics Simulations of the Cytochrome c,-Rubredoxin Complex From Desulfovibrio vulgaris David E. Stewart' and John E. Wamplel.2 'University Computing and Networking Services, Specialized Systems Support, Computer Services Annex, and 'Department of Biochemistry, Boyd Graduate Studies Research Center, The University of Georgia, Athens, Georgia 30602
ABSTRACT Molecular dynamics simulations have been carried out on the complex formed between the tetraheme cytochrome c3 and the iron protein rubredoxin from the sulfate-reducing bacterium Desulfovibrio vulguris. These simulations were performed both with explicit solvent water molecules included, and without solvent molecules using a distancedependent dielectric constant to approximate the screening effects of solvent. The results of both simulations are strikingly different, indicating that the representation of environmental effects is important in such simulations. F o r example, a striking adaptation of the two proteins seen in the nonsolvated simulation is not seen when explicit solvent water is included; in fact, the complex appears to become weaker in the solvated simulation. Nonetheless, the iron-iron distance decreases more significantly in the solvated simulation than in the nonsolvated simulation. It was found that in both cases molecular dynamics optimized the structures further than energy minimization alone. Key words: AMBER, electrostatics, proteinprotein interactions, intermolecular bonding, electron transfer, redox proteins INTRODUCTION The importance of biological electron transfer in the energy metabolism of living systems has made it the subject of intense research. One aspect critical for a full understanding of this process is the manner in which electron transfer proteins interact. Based on spectroscopic evidence and the analysis of available X-ray structures we have proposed models for electron transfer complexes of small proteins from the sulfate-reducing bacterium Desulfovibrio vulgaris. The supporting evidence for these models has been discussed In this report one of these models has been refined using molecular dynamics (MD) simulations. Cytochrome c3 from Desulfovibrio is unique among the c-type cytochromes in that there are four hemes per polypeptide chain of approximately 110 amino acid residues. The hemes exhibit distinct re0
1991 WILEY-LISS, INC.
dox properties with negative redox potentials as low as -340 mV.5 In the two known crystal structures627 the hemes are well separated (average Fe-to-Fe distance of 14.5 A) and are arranged with a particularly obvious lack of planar alignment (average heme-to-heme angle of 72"). Rubredoxins from these same organisms are much smaller, with between 50 and 55 residues, and contain one nonheme iron atom ligated by four cysteine sulfhydryl groups. The physiological function of rubredoxin has not yet been determined,' but it has been shown t o act as an effective redox partner for cytochrome c3 in ~ i t r oRu.~ bredoxins exhibit redox potentials in the range of 0 to -50 mV.l0 The physiological function of cytochrome c3 is believed to be as an electron carrier for the enzyme hydrogenase.' Rubredoxin can be rapidly reduced in vitro by cytochrome c3 in the presence of hydroge n a ~ e Thus, .~ while no physiological significance can be claimed for the proposed complex, such models can serve to aid investigation of electron transfer. The surface of rubredoxin in the region of its redox active iron atom is uniformly negatively ~ h a r g e d . ~ , ~ Although each of the four hemes of cytochrome c3 is exposed t o the surface, only one (heme 1) is surrounded by a complementary uniformly positive charged It has been shown that many redox proteins interact via electrostatic interact i o n ~ . ' ' - ~On ~ the basis of the pattern of apparent electrostatic complementarity, along with steric and geometric considerations, this heme has been proposed to be the site a t which rubredoxin and flavodoxin interact.lP3 NMR data indicate that 1:l complexes are formed between the interacting proteins and that spectral lines associated with the same heme group are perturbed in both cases.2 Computational studies have been undertaken pre-
Received October 1, 1990; revision accepted February 20, 1991. Address reprint requests to David E. Stewart, The University of Georgia, University Computing and Networking Services, Specialized Services Support, Computer Services Annex, Athens, GA 30602.
143
DYNAMICS OF THE RUBREDOXIN-CYTOCHROME c3 COMPLEX 1
"
"
"
"
'
I
-05-150 -
-"""I m
x
m
-3500
L
al
-
c w c
L W
-3700
0
.d
3
c-'
rn
c
-166001
0 0 L 3
B
u 0
a
-150
W c,
a,
c
-16800
-250 -350
-17000 -17200
-5501 ,
-17400 10
40
70
100
,
,
10
viously on the interactions of redox protein^.'^-^^ A molecular dynamics study similar to that reported here was carried out by Wendoloski et al.15 on the interaction between cytochromes c and b,. Again, while this interaction has no physiological significance-the pairing modeled involved proteins from different organisms-an in vitro interaction had previously been demonstrated.l8 In the cytochrome c-cytochrome b, study it was observed that the total potential energy and the heme-to-heme distance decreased during the simulation. This observation indicated that, in addition to providing details about movements of atoms, the dynamics simulation optimized the structure of the complex more fully than energy minimization alone. Molecular dynamics has been shown to be effective in optimizing structures to a greater extent than energy minimization in other cases as well.19-22The advantage of molecular dynamics is its inherent ability to explore a larger region of the conformational space available to the molecule or complex than energy minimization since the kinetic energy of the atoms in MD simulations allows local energy barriers to be crossed.23 In the case of the cytochrome c-cytochrome b, model,15 a striking change occurred in the interface between the two proteins as the simulation progressed. The heme groups turned from a edge-on configuration and, along with a phenylalanine resi-
,
,
40
Time (ps) Fig. 1. Change in total AMBER potential energy of the complex. (A) During the nonsolvated simulation; ( 8 ) during the solvated simulation.
,
, 70
,
,
'I1 100
Time (ps) Fig. 2. Nonbonded interaction energies. (A) During the nonsolvated simulation; (B) during the solvated simulation.
due of cytochrome c, a stack of aromatic rings was formed. MD simulations may be performed using a distance-dependent dielectric constant to approximate the screening of electrostatic interactions by solvent, or with explicit water molecules included to represent the aqueous environment of the protein. The merits of each of these approaches have been discussed previously. One of the implications of the cytochrome c-cytochrome b, sir nu la ti on^'^ was that the two approaches can lead to similar results in the interface region of a protein-protein complex. In this study, we compare the results obtained from MD simulations using a distance-dependent dielectric constant to those in which explicit water molecules were included, focusing particularly on the interface of the complex. 15724--26
METHODS The program AMBER 3.0 UCSFZ7was used for the molecular mechanics and dynamics calculations and for the analysis of the results. Minimization and MD calculations were performed on the University of Georgia Cyber 205 supercomputer, and analysis and display of the results were performed on a Silicon Graphics Iris workstation. The version of AMBER 3.0 used was optimized to run on the Cyber 205. It has been found that a 100-fold increase in computational speed can be realized by vectorization
144
D.E. STEWART AND J.E. WAMPLER
-41ooLl
-1
A
-4300 -4500 x Is)
L
-4700
W
c
c
W
c 0 -+
0
-4900
c, 0
m
L
-1600
rn
c,
c
H
a
-1800 -2000
2.1
1.7
.gr . i I
-2200 -2400
.5
10
40
70
100
10
Fig. 3. Water-protein interaction energies during the solvated simulation. (A) rubredoxin; (B) cytochrome c,.
0
=I
191 L
70
100
Time (ps)
T i m e (ps)
1,
40
1
Structure # Fig. 4. (A) Distance (in A) between the centers of mass of the two proteins during solvated simulation (solid line, square symbols) and during the nonsolvated simulation (dashed line, circles). (B) Radius of gyration (in A) for the complex during the solvated simulation (solid line) and the nonsolvated simulation (dashed line).
Fig. 5. Root-mean-squared (rms) deviation for all atoms from the initial (0.1 psec) structure. (A) During the nonsolvated simulation; (B) during the solvated simulation.
of the program on such high performance computer~.~,~' The starting structure has been described previously.2 The coordinates of the proteins were obtained from the Brookhaven Protein Data Bank.29,30 The structures used in the complex were rubredoxin from D. vulgaris3' and cytochrome c3 from D. vulgaris strain M i y a ~ a k i . ~ ~ In the nonsolvated simulation the energy of the model-built complex was minimized until its rootmean-square (rms) gradient of the energy was less than 0.025 kcal/(mol A). A distance-dependent dielectric was used with a nonbonded cutoff of 9.5 A. Initially 125 cycles of constrained energy minimization were performed to gradually relax steric conflicts without introducing large distortions in the structure. This was followed by unconstrained conjugate gradient minimization until the limit of convergence was reached. A united-atom representation of the hydrogen atoms was used, with only the hydrogen atoms on polar atoms explicitly included,33 resulting in a total of 1,661 atoms (269 hydrogen 1,392 heavy atoms) for the complex. The energy-minimized structure was used as the starting coordinates for the MD run. Initial velocities assigned to the minimized structure correspond to a mean temperature of 50 K. The heating phase consisted of 30 K temperature increments applied every 200 fsec to a final temperature of 300 K. Equilib-
+
DYNAMICS OF THE RUBREDOXIN-CYTOCHROME c3 COMPLEX
145
A
5 ps.
30 ps.
B
5 ps.
30 ps.
Fig. 6. The adaptive changes in the two proteins over the course of the nonsolvated simulation are illustrated by two structures from the simulation shown with orthogonal, stereo views. In both cases the upper structures is that at 5 psec and the lower structure is from 30 psec. (A) View along binding cleft showing
symmetry of the rubredoxin (on left) with the iron center represented by a dot for the iron atom. (B) Orthogonal view to A showing how the beta sheets at the top of the two proteins are drawn together by formation of ionic and hydrogen bonds as listed in Table IA.
rium at 300 K was ensured by repeatedly reassigning random velocities from the appropriate Maxwellian distribution function every 200 fsec for a period of 2.0 psec. The system was then allowed to stabilize for a period of 4.0 psec after which data collection was begun for 65 psec. In the solvated simulation a 6 A shell of 1001 TIP3P water molecules was added to the model-built complex.34 A constant dielectric of E = 1 and nonbonded cutoff of 9.5 A was used. This system was
minimized for 3,500 cycles. The minimized structure was heated from 0 K to 300 K in 30 K increments applied every 200 fsec. The system was equilibrated at 300 K for 2.0 psec, followed by a stabilization period of 10 psec. Data collection continued for 100 psec. In both simulations a 0.001 psec time step was used. Bond lengths were fixed at their optimal values using the SHAKE algorithm.35 Structures were collected every 0.1 psec (100 steps). Trajectories and parameters for analysis were ob-
146
D.E. STEWART AND J.E. WAMPLER
A
100 ps.
B
Fig. 7. Orthogonal, stereo views of the initial (1 psec) and final (100 psec) structures from the solvated simulation. (A) View along binding cleft showing symmetry of the rubredoxin (on left) with the iron center represented by a dot for the iron atom. (B) Orthogonal view to A. In both cases little overall change is seen and this result is characteristic of the entire simulation.
tained using the AMBER analysis module. Most manual manipulations, display, and plotting were done using the program SYBYL (Tripos Associates, Inc., St. Louis, MO). Animation of the dynamics simulations was carried out using QUANTA (Polygen Corporation, Waltham, MA). Plotting and manipulations of the data sets representing the trajectories were done using SPECOS SA.36
RESULTS The General Character of the Complex The changes in the potential energy during the simulations of the complexes are shown i n Figure 1. The conformations of the complexes from the MD simulations were energetically more favorable than those from energy minimization. The potential energy of the complex rapidly becomes more favorable
DYNAMICS OF THE RUBREDOXIN-CYTOCHROME c3 COMPLEX
147
TABLE I. Status of Interface Bonds During the MD Simulation Intermolecular bond ( c , -+ rd) residue (atom*) K 57 (NZ) - E 17 (OE) K 58 (NZ) - E 17 (OE) K 58 (NZ) - D 19 (OD) K 95 (NZ) - E 50 (OE) KlOl (NZ) - E 12 (OE) K 15 (NZ) - E 17 (OE) K 94 (NZ) - E 50 (OE) D 71 (OD) - C 42 (N) K 57 (NZ) - E 12 (0) K 60 (NZ) - P 40 (0) K 72 (NZ) - V 41 (0) K 94 (NZ) - T 7 (0) K 94 (NZ) - V 8 (0) ClOO (SG) - E 12 (N) ClOO (SG) - Y 11 (N) K 58 (NZ) - G 18 (0) K 58 (NZ) - Y 11 (OH)
Y 66 (OH) - Y 11 (N) K 52 (NZ) - V 8 (0)
Present in Present min’d struct. at end Comments IA. During the Nonsolvated Simulation X X Maintained throughout X X Maintained throughout X X Maintained throughout X X Maintained throughout X X Maintained throughout X Formed between 30 and 35 psec X Formed between 25 and 30 psec X Formed during stabilization period X Formed between 25 and 30 psec X Broken between 0 and 5 psec Broken during stabilization period X Broken between 25 and 30 psec X X Broken between 0 and 5 psec Broken between 20 and 25 psec X Broken between 20 and 25 psec X Formed during stabilization period X X Formed during equilibration period Broken during stabilization period Reformed between 10 and 15 psec X Formed between 20 and 25 psec X Formed between 0 and 5 psec IB. During the Solvated Simulation
K 57 (NZ) - E 17 (OE) K 58 (NZ) - D 19 (OD) K 95 (NZ) - E 50 (OE) KlOl (NZ) - E 12 (OE) K 58 (NZ) - T 21 (OG) K 60 (NZ) - P 40 (0) K 72 (NZ) - V 41 (0) K 94 (NZ) - T 7 (0) K 94 (NZ) - V 8 (0)
X X X X X X X X X
X
Broken during equilibration Maintained throughout Broken during equilibration Broken during equilibration Broken during first 5 psec of stabilization Broken during equilibration Broken during first 5 psec of stabilization Broken during equilibration Broken between 30 and 35 psec
*Residues identified by standard one letter code. Number represents position in protein sequence (from N-terminal). Atom types are as follows: NZ, epsilon amino of lysine; OE and OD, terminal carboxyl group oxygen of glutamic and aspartic acid residues; N, main chain amide hydrogen; 0, main chain amide carbonyl; SG, terminal sulfhydryl of cysteine.
in approximately the first 35 psec of the data collection period in the nonsolvated simulation, while such a change is not seen in the solvated complex until after approximately 50 psec. It is apparent from the figure that the potential energy is much more negative in the solvated simulation than in the nonsolvated simulation. This arises from the inclusion of the water molecules and to a larger extent from the value of the dielectric constant used: in the nonsolvated simulation the charges within the protein are shielded due to the distance-dependent dielectric; in the solvated simulation the dielectric constant is 1 and no such explicit shielding occurs. The strength of the complex may be evaluated by calculating the nonbonded interaction energies (electrostatic, van der Waals, and hydrogen bonding) between rubredoxin and cytochrome cg. These are shown in Figure 2. The hydrogen bonding and van der Waals interaction energies are similar in both simulations; however, there is a significant difference in the electrostatic interaction energies. This
again arises mainly from the dielectric constant used. The total interaction energy in the nonsolvated simulation generally becomes more favorable during the simulation, particularly in the 30-35 psec time frame. Conversely, in the solvated simulation, no such trend is observed, instead, the rubredoxin-cytochrome cg interaction energy becomes less favorable until the latter portion of the simulation (90100 psec). The driving force for the unfavorable change is apparently the interactions between the solvent molecules and the protein (Fig. 3). The total interaction energy between the three “groups” of molecules-rubredoxin, cytochrome c3, and waterbecomes more favorable during the simulation due to the contributions from the solvent-protein interactions. The specific interactions (see below) seem to indicate that the two proteins are less directly involved as the solvated simulation progresses. Two gross measures of this are the distance between the centers of masses of the two proteins (Fig. 4A) and the
148
D.E. STEWART AND J.E. WAMPLER
radius of gyration of the complex (Fig. 4B). In both simulations the distances between the centers of mass change comparably through approximately 30 psec at which time this distance begins to increase steadily in the solvated simulation, indicating that the two proteins are separating. The complex radius of gyration decreases initially only in the nonsolvated simulation, probably due to the strong electrostatic interactions arising from the use of a distance-dependent dielectric constant. No such change is seen in the solvated simulation. In Figure 5, the rms deviations of the atoms during the simulations are shown. In both simulations, a maximum rms deviation of slightly less than 3 A is observed, however, the time required for the simulations to reach this value is different. The rms deviation increases most rapidly in the nonsolvated simulation, due to the absence of water molecules to damp the motions of the atoms. In this case, the rms deviation appears to stabilize after approximately 35 psec suggesting formation of a stable conformation. The final structures from each simulation are significantly different, with an rms deviation of 4.4 A between the 65 psec nonsolvated structure and the 100 psec solvated structure for all atoms; this decreases to 4.2 A if the hydrogen atoms are omitted from the calculation. The rms deviation between rubredoxin in the two final structures is 3.4 A; for cytochrome c3 this is 3.5 A (heavy atoms only). The lower rms deviation for the individual proteins as compared to the complexes suggests that there is a significant difference in the orientation of the molecules in the complex between the solvated and nonsolvated structures. In a qualitative sense the overall interpretation of these results is illustrated with two structures from both simulations as shown in Figures 6 and 7. Under the drive of the electrostatic effect, the two proteins change shape during the nonsolvated simulation (Fig. 6). The net result is that both the binding depression on the cytochrome and the shape of the rubredoxin become more symmetrical. On the large scale most of these changes represent “sliding” motions of beta sheet portions of both proteins. During the solvated simulation (Fig. 7), these changes do not take place and the shapes of both proteins stay more or less constant. The intercalation of water and loss of electrostatic bonding (see below) prevent these changes from taking place over the time period of the simulation.
Interactions in the Interface In the nonsolvated simulation, the minimized structure has 5 ionic and 5 hydrogen bonds, while the 0.1 psec data collection structure has 5 ionic and
I (u
U
13 12 11 10
10
40
70
100
Time (ps) Fig. 8. Iron-iron distances between the iron center of rubredoxin and the iron of heme 1 of cytochrome c, (distance in A). (A) During the nonsolvated simulation; (B) during the solvated simulation.
tion as discussed above. It has been suggested that an inordinate number of such bonds may be formed due to the lack of inclusion of explicit ~ o l v e n t . ~ ~ , ~ ~ This criticism is confirmed by the number of intermolecular ionic and hydrogen bonds found for the solvated structures: there are 4 ionic and 6 hydrogen bonds in the minimized structure, 2 ionic and 2 hydrogen bonds in the 0.1 psec data collection structure, and 1 ionic bond in the 100 psec structure. These changes are caused by the intercalation of water molecules into the interface between the two molecules disrupting the intermolecular bonds. The fate of these bonds during the simulations is summarized in Table I. Electron transport is dependent on the orientation and distance between redox group^.^^-*^ One criterion to judge the optimization of the structures resulting from these simulations, then, would be the distance between prosthetic groups as a function of time, with the later structures expected to have a smaller iron-to-heme distance. In the nonsolvated simulation, the iron-to-heme distance decreased slightly (Fig. 8). A more significant decrease was seen in the solvated simulation (Fig. €9, with the distance being much smaller a t the end of this simulation than in the nonsolvated case.
7 hydrogen bonds, and the 65 psec structure has 7
A striking rearrangement is seen in the organiza-
ionic and 6 hydrogen bonds. This is reflected in the increasing electrostatic energy during the simula-
tion of residues in the interface between the two proteins in the nonsolvated simulation. This is
DYNAMICS OF THE RUBREDOXIN-CYTOCHROME
149
COMPLEX
f
5 ps.
Y
30 ps.
Y
Y
a Y
(
Fig. 9. Interactions in the vicinity of heme 1 of cytochrome c, at the interface between the two proteins. Residues from rubredoxin are shown with the dithered line. The rearrangement of ionic and hydrogen bonds summarized in Table IA is accompanied by changes in mobility and orientation of the aromatic residues (see
Fig. 10). Residues shown for cytochrome c, are heme 1, histidine 70 and 106, tyrosine 66, lysines 15,57, and 58,and aspartate 56. Residues of rubredoxin shown are tyrosine 1 I , glutamate 17, and aspartate 19.
shown in the structures of Figure 9. To examine the fluctuations of the sidechains of the aromatic residues, the value of the dihedral angle x2 (C,-C,C,-C,1) was calculated for each of the 0.1 psec data collection structures. The trajectories of dihedral angle rotation of the residues involved (Fig. 10) showed a concerted change at 30 psec and obvious coupling between the aromatics as revealed by the changes in amplitude of their variations. This change is driven by a rearrangement of hydrogen and ionic bonding partners which occurs in a concerted fashion in a few picoseconds. The change is completed by 30 psec. The general conformation and configuration of the two molecules before this change is illustrated by the structure a t 5 psec shown in Figure 9 (top). Similarly, after the transition has occurred, the structure (illustrated by the 30 psec structure in Fig. 9, bottom) remains fairly constant for the remaining 35 psec of the simulation. It is interesting to note that while this transition occurs under the inf hence of a rearrangement of electrostatic interactions, that the major changes involved are interactions be-
tween residues within each of the proteins rather than between them. Thus, the ionic bonds of the heme propionate residue with both the main chain amide and €-amino of the nearby lysine (number 15) are changed into a n ionic bridge involving lysine 15, the heme propionate, and lysine 57. Similarly, tyrosine 11or rubredoxin flips over so that it bonds to the opposite oxygen of the carboxyl of glutamate 17 leaving the lower oxygen (as seen in the figure) free to also bind to lysine 57 of the cytochrome. The trajectories of the aromatic residues during the solvated simulation (Fig. 11)reveal no such coupled changes and the ionic and hydrogen bonds between the proteins are, for the most part, broken (Table I).
DISCUSSION While it has been shown that the exclusion of solvent molecules in MD calculations may result in the formation of improper hydrogen bond^'^,^^ and that the atomic motions exhibit artificially large amplitUdes,15,26,43 a widely held opinion is that simulations with a distance-dependent dielectric approxi-
150
D.E. STEWART AND J.E. WAMPLER
180
I
L TTyr y r 11 120 (ru) 60 I
i i
O?
k
-60 - 6 0 r TTyr y r 66 I
-120
t-I
(c,) (c3) ' Y '
-180-'
'
I
'
I
10
0
'
I
I
30
20
'
40
I
I
50
60
-
Time (ps) Fig. 10. Trajectories of dihedral angle xz changes in interface during the nonsolvated simulation. Angles for tyrosine 11 are shown offset ( + 90").
-180k'
'
I
'
I
0
I
10
'
I
20
I
30
I
40
'
50
I
I
60
I
70
I
80
I ' - I 90 100
Tine (ps) Fig. 11. Trajectories of dihedral angle x2 of same residues shown in Figure 10 during the solvated simulation. Angles for tyrosine 66 are shown offset (-330").
mation may still be instrumental in understanding the dynamic behavior of proteins. In the previous study of the cytochrome c-cytochrome b, complex,15 fully solvated and nonsolvated simulations lead to
very similar changes within the interface between the proteins, but on a much slower time scale for the solvated simulation. In addition, since the atomic motions are not damped by the solvent, a larger por-
DYNAMICS OF THE RUBREDOXIN-CYTOCHROME c3 COMPLEX
tion of conformational space may be searched by such sir nu la ti on^.'^.^^ Our results are less clearly in support of this idea that nonsolvated and solvated simulations lead to similar changes only on a different time scale. In both simulations, the resulting structure had a significantly lower potential energy than did the minimized structure. However, the interaction energies between the two proteins exhibited opposite trends. In the nonsolvated simulation, the interaction energy became more favorable, while in the solvated simulation the converse was true until the 90 psec time period. This unfavorable change in the interaction energy appears to be driven by favorable interactions between the proteins and the water molecules. The obvious explanation for this difference is intercalation of water molecules into the interface between the two proteins, solvating the ionic bonds. The consequences of these changes are summarized in Table I. The resulting difference is clear in the stereo plots of Figures 6 and 7. Without solvent, the proteins adapt to each other in a n intriguing way forming a fairly symmetrical “lock-and-key’’ type of complementary fit. With solvent they remain independent entities. The suggestion that the water molecules intercalate into the interface of the complex is supported not only by the breaking of the intermolecular bonds (and associated decrease in electrostatic interaction energy) but also by the change in the radius of gyration (R,). The R, of the nonsolvated complex decreases quickly, as the number of intermolecular bonds increases, while the R, of the solvated complex remains stable, increasing slightly in the latter half of the simulation. Overall, both the R, and the distance between centers of mass for the solvated complex increase. Extrapolating these two trajectories to longer time suggests that the two molecules are drifting apart although the iron-iron distance is actually decreasing. The comparison between early and late structures for both simulations is shown in Figures 6 and 7. The changes seen are in essentially opposite directions. Without solvent, rubredoxin flattens and adapts to the surface of the cytochrome. In the view shown, this involves a counterclockwise rotation of rubredoxin relative to cytochrome c,. With solvent, rubredoxin becomes slightly more spherical and appears t o rotate slightly clockwise in this view. Molecular dynamics simulations of the proposed cytochrome c,-rubredoxin complex reveal that the complex is stable on the picosecond time scale under the influences of the AMBER force field. Spectroscopic evidence’ supports both the specificity and the close hemeinonheme iron interactions modeled. However, the manually docked proteins from our previous work1-, are obviously not in a n optimum configuration since both simulations lead to lower
151
energies. The two simulations lead to different results. In the nonsolvated case the proteins “roll” together and change shape in an adaptive way. In the solvated simulation they “roll” in the opposite direction on trajectories which may eventually lead to disruption of the complex. It is possible that the nonsolvated simulation probes areas of conformational space that are not accessible during the 100 psec time period of the solvated simulation. The interesting changes in the interface between the two proteins might be similarly accessible. One approach to investigation of this question would be to solvate the terminal structure from the nonsolvated simulation and see if the aromatic stacking and electrostatic interaction (ionic and hydrogen bonds) are disrupted in this case. Within the limits of available computer time, this study and similar simulations of the cytochrome c,-flavodoxin complex2 are currently being pursued. A wealth of data comes from these simulations. Every residue in both proteins undergoes changes in relative orientation and in its interactions with others. A complete analysis of the changes in the aromatic residues within these structures (manuscript in preparation) reveals that coupling and stacking between residues like those shown in Figure 9 are common.
ACKNOWLEDGMENTS The authors thank Dr. Walter McRae (University Computing and Networking Services) and the University of Georgia Research Foundation for their support of this work. Vectorization of major portions of the AMBER suite of programs for operation on the Cyber 205 was performed by Dr. Steve Gallion. This research was supported in part by University Computing and Networking Services, by grants from the University of Georgia Research Foundation, and the National Institutes of Health (Grant GM 41482).
REFERENCES 1 Stewart, D.E., LeGall, J., Moura, L., Moura, J., Peck, H.D., Jr., Xavier, A.V., Weiner, P.K., Wampler, J.E. A hypothetical model of the flavodoxin-tetraheme cvtochrome c , complex of sulfate-reducing bacteria. Biochkmistry 27:54442450. 1987. 2. Stewak,D.E., LeGall, J., Moura, I., Moura, J., Peck, H.D., Jr., Xavier, A.V., Weiner, P.K., Wampler, J.E. Molecular modeling and NMR studies of the rubredoxin-cytochrome cg complex of sulfate-reducing bacteria. Eur. J . Biochem 185:695-700, 1989. 3. Stewart, D.E. The structure, interactions, and dynamics of electron transport proteins from Desulfouibrzo. A molecular modeling and computational study. Ph.D. Dissertation, The University of Georgia, 1989. 4. Wampler, J.E., Stewart, D.E., Gallion, S.L. Molecular dynamics simulations of proteins and protein-protein complexes. In: “Computer Simulation Studies in Condensed Matter Physics 11.” Landau, D.P., Mon, K.K., Schuttler, H.G. (eds.). Berlin: Springer-Verlag, 1990: 68-84. 5. Yagi, T., Inokuchi, H., Kimura, K. Cytochrome cj, a tetrahemoprotein electron carrier found in sulfate-reducing bacteria. Acc. Chem. Res. 16:2-7, 1983. 6. Higuchi, Y., Kusunoki, M., Matsuura, Y., Yasuoka, N.,
152
7.
8. 9.
10. 11. 12.
13. 14.
15. 16.
17. 18. 19.
20. 21. 22. 23. 24. 25. 26.
D.E. STEWART AND J.E. WAMPLER
Kakudo, M. Refined structure of cytochrome c3 at 1.8 A. J . Mol. Biol. 172:109-139, 1984. Haser, R., Pierrot, M., Frey, M., Payan, F., Astier, J.P., Bruschi, M., LeGall, J . Structure and sequence of the multiheme cytochrome c3. Nature (London) 282:806-810, 1979. LeGall, J., Fauque, G. Dissimilatory reduction of sulfur compounds. In: “Biology of Anaerobic Microorganisms.” Zehnder, A.J.B. (ed.). New York: Wiley, 1988: 587-639. Bell, G.R., Lee, J.P., Peck, H.D., Jr., LeGall, J . Reactivity of Desulfovibrio gigas hydrogenase toward artificial and natural electron donors or acceptors. Biochimie 60:315320,1978. Moura, I., Moura, J.J.G., Santos, M.H., Xavier, A.V., LeGall, J . Redox studies on rubredoxins from sulfate and sulfur reducing bacteria. FEBS Lett. 107:419-421, 1979. Salemme, F.R. A hypothetical structure for a n intermolecular electron transfer complex of cytochromes c and b,. J . Mol. Biol. 102:563-568, 1976. Smith, H.T., Staudenmayer, N., Millet, F. Use of specific lysine modifications to locate the reaction site of cytochrome c with cytochrome c oxidase. Biochemistry 16: 4971-4974, 1977. Mauk, M.R., Reid, L.S., Mauk, A.G. Spectrophotometric analysis of the interaction between cytochrome b, and cytochrome c. Biochemistry 21:1843-1846, 1982. Erman, J.E., Vitello, L.B. The binding of cytochrome c peroxidase and ferricytochrome c. A spectrophotometric determination of the equilibrium association constant as a function of ionic strength. J . Biol. Chem. 255:6224-6227, 1980. Wendoloski, J.J., Matthew, J.B., Weber, P.C., Salemme, F.R. Molecular dynamics of a cytochrome c-cytochrome b, electron transfer complex. Science 238794-797, 1987. Northrup, S.H., Boles, J.O., Reynolds, J.C.L. Electrostatic effects in the Brownian dynamics of association and orientation of heme proteins. J . Phys. Chem. 915991-5998, 1987. Northrup, S.H., Boles, J.O., Reynolds, J.C.L. Brownian dynamics of cytochrome c and cytochrome c peroxidase association. Science 241:67-70, 1988. Mauk, M.R., Reid, L.S., Mauk, A.G. Spectrophotometric analysis of the interaction between cytochrome b, and cytochrome c. Biochemistry 21:1843-1846, 1982. Brooks, B.R., Bruccoleri, R.E., Olafson, B.D., States, D.J., Swaminathan, S., Karplus, M. CHARMM: A program for macromolecular energy, minimization, and dynamics calculations. J . Comp. Chem. 4:187-217, 1983. Levitt, M. Protein folding by restrained energy minimization and molecular dynamics. J . Mol. Biol. 170:723-764, 1983. Elber, R., Karplus, M. Multiple conformational states of proteins: A molecular dynamics analysis of myoglobin. Science 235318-321, 1987. Brunger, A.T., Kuriyan, J., Karplus, M. Crystallographic R factor refinement by molecular dynamics. Science 235: 458-460, 1987. van Gunsteren, W.F. The role of computer simulation techniques in protein engineering. Protein Eng. 25-13, 1988. Brooks, C.L., Karplus, M. Solvent effects on protein motion and protein effects on solvent motion. J . Mol. Biol. 208:159-181, 1989. Wendoloski, J.J., Matthew, J.B. Molecular dynamics effects on protein electrostatics. Proteins 5313-321, 1989. Levitt, M., Sharon, R. Accurate simulation of protein dy-
27.
28.
29.
30.
31.
32. 33.
34.
35. 36. 37.
38.
39.
40. 41. 42. 43.
namics in solution. Proc. Natl. Acad. Sci. U.S.A. 8575577561,1988. Singh, U.C., Weiner, P.K., Caldwell, J.W., Kollman, P.A. AMBER (UCSF version 3.0). Department of Pharmaceutical Chemistry, University of California, San Francisco, 1986. Gallion, S.L., Levy, R.M., Weiner, P.K., Hirata, F. Implementation of a macromolecular mechanics program on a CYBER 205 supercomputer. Comp. Chem. 10:165-173, 1986. Bernstein, F.C., Koetzle, T.F., Williams, G.J.B., Meyer, E.F., Jr., Brice, M.D., Rodgers, J.R., Kennard, O., Shimanouchi, T., Tasumi, M. The Protein Data Bank: A computer-based archival file for macromolecular structures. J . Mol. Biol. 112535-542, 1977. Abola, E.E., Bernstein, F.C., Bryant, S.H., Koetzle, T.F., Weng, J. Protein Data Bank. In: “Crystallographic Databases-Information Content, Software Systems, Scientific Applications.” Allen, F.H., Bergerhoff, G., Sievers, R. (eds.). Bonn: Data Commission of the International Union of Crystallography, 1987: 107-132. Adman, E.T., Sieker, L.C., Jensen, L.H., Bruschi, M. LeGall, J . A structural model of rubredoxin from Desulfovibrio vulgaris at 2 Angstrom resolution. J. Mol. Biol. 112:113-120, 1977. Higuchi, Y., Kusunoki, M., Matsuura, Y., Yasuoka, N., Kakudo, M. Refined structure of cytochrome c3 at 1.8 Angstrom resolution. J. Mol. Bio. 172:109-139, 1984. Weiner, S.J., Kollman, P.A., Case, D.A., Singh, U.C., Ghio, C., Alagona, G., Profeta, S., Weiner, P. A new force field for molecular mechanical simulation of nucleic acids and proteins. J. Am. Chem. SOC.106:765-784, 1984. Jorgensen, W.L., Chandrasekhar, J., Madura, J., Impey, R.W., Klein, M.L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79:926935, 1983. van Gunsteren, W.F., Berendsen, H.J.C. Algorithms for macromolecular dynamics and constraint dynamics. Mol. Phys. 341311-1327, 1977. Wampler, J.E. SPECOS SA (The Stand Alone Version) User’s Guide, Instrument Software Services, Arnoldsville, GA, 1989. Levitt, M., Sharon, R. Simulating protein dynamics in solution: Bovine pancreatic trypsin inhibitor. In: “Crystallography in Molecular Biology.” Moras, D., Drenth, J., Strandberg, B., Suck, D., Wilson, K. (eds.). New York: Plenum, 1985 197-205. Mauk, A.G., Scott, R.A., Grav. H.B. Distance of electron transfer to and from metalliprotein redox sites in reactions with inorganic complexes. J . Am. Chem. Soc. 102: 4360-4363, 1980. Makinen, M.W., Schichman, S.A., Hill, S.C., Gray, H.B. Heme-heme orientation and electron transfer kinetic behavior of multisite oxidation-reduction enzymes. Science 222:929-931,1983. Siders, P., Cave, R.J., Marcus, R.A. A model for orientation effects in electron-transfer reactions. J. Chem. Phys. 81: 5613-5624,1984. Scott, R.A., Mauk, A.G., Gray, H.B. Experimental approaches to studying biological electron transfer. J. Chem. Ed. 62:932-938, 1985. Cave, R.J., Siders, P. Marcus, R.A. Mutual orientation effects on electron transfer between porphyrins. J . Phys. Chem. 90:1436-1444,1986. Chen, L.X.-Q., Engh, R.A. Brunger, A.T., Nguyen, D.T., Karplus, M., Fleming, G.R. Dynamics studies of apoazurin of Alcaligenes denitrificans. Biochemistry 27:6908-6921, 1988.
