INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Ubiquinol/cytochrome c oxidoreductase (cytochrome bc₁ complex) plays a central role in the electron transport chains of bacteria, mitochondria, and chloroplasts, converting redox free energy into a proton gradient that drives the cell's metabolism through ATP synthesis. The bc₁ complex is found in the plasma membrane of bacteria and in the inner mitochondrial membrane of eukaryotes. It catalyzes the oxidation of ubiquinol in the membrane, as well as the reduction of cytochrome c and the translocation of protons across the membrane (Gennis et al., 1993; Brandt and Trumpower, 1994; Brandt, 1997; Crofts and Berry, 1998). All bc₁ complexes contain three essential subunits to which prosthetic groups are bound: a cytochrome b with high- and low-potential hemes b_H and b_L, a "Rieske" iron-sulfur protein (ISP) containing an Fe₂S₂ cluster, and a cytochrome c₁ with another heme group.

Crystal structures of several mitochondrial bc₁ complexes with and without inhibitors bound at the Q_o and Q_i sites have become available recently (Xia et al., 1997; Zhang et al., 1998; Iwata et al., 1998). The structure of the bc₁ complex is a dimer, each monomer of which is composed of 10 or 11 different polypeptide subunits. Cytochrome b is located predominantly in the transmembrane region of the complex. Cytochrome c₁ comprises a water-soluble part and a C-terminal transmembrane -helix. The ISP subunit of each monomer consists of an N-terminal anchoring transmembrane -helix and a water-soluble head containing the Fe₂S₂ cluster, which participates in the functional cycle of the other monomer. The matrix side of the bc₁ complex is composed of several domains that account for about half of the mass of the complex, but are not involved in the catalysis.

The catalytic mechanism of the bc₁ complex illustrated in Fig. 1 involves two catalytic sites for oxidation (Q_o site) or reduction (Q_i site) of the quinones, and is known as the modified Q-cycle (Mitchell, 1976; Crofts, 1985). The complete cycle is initiated when a quinol QH₂ binds to the Q_o site, releases two protons into the intermembrane region, and transfers two electrons to two redox sites: one electron to the Fe₂S₂ cluster, from where it is transferred to cytochrome c₁ and subsequently consumed by cytochrome c₂; the second electron is transferred to cytochrome b_L, then to cytochrome b_H and consumed at the Q_i site. The overall reaction involves oxidation of two QH₂ molecules at the Q_o site that releases four protons into the intermembrane space, and the formation of QH₂ at the Q_i site that utilizes two protons from the matrix side of the membrane (N-side). The question arises how the bc₁ complex splits the electron pathways from the Q_o site to reduce two different redox centers, i.e., the Fe₂S₂ cluster and the heme b_L.

View larger version (50K): [in this window] [in a new window]   FIGURE 1   Catalytic mechanism of bc1 complex. Two cycles of QH2 oxidation at the Qo site are required to generate the two-electron reduction of Q to QH2 at the Qi site. The shaded area represents the membrane.

The available crystal structures of the bc₁ complexes suggest an answer to this question. In fact, in the earlier structure of the bc₁ complex from bovine heart mitochondria, the soluble part of the ISP subunit could not be resolved. This was attributed to a high mobility of the ISP in the crystal; the electron density was sufficient, however, to identify the Fe₂S₂ cluster at a position 31 Å from the iron of cytochrome c₁, and 27 Å away from the iron of heme b_L (Xia et al., 1997), close to the Q_o binding site. The ISP soluble head assumed a conformation with a similar position of the Fe₂S₂ cluster in the structure of the bc₁ complex from chicken heart mitochondria (Zhang et al., 1998) with an inhibitor (stigmatellin) bound at the Q_o site (proximal position). In this structure the His-161 of the ISP coordinating the Fe₂S₂ cluster formed a hydrogen bond to stigmatellin. However, in the structure of the native bc₁ complex from chicken (Zhang et al., 1998) the soluble part of the ISP domain was rotated with respect to its position in the structure with the inhibitor, as shown in Fig. 2. In this conformation (distal position) the Fe₂S₂ cluster was located 21.3 Å from cytochrome c₁. In P6₅22 crystals from bovine heart, the Fe₂S₂ cluster was found sufficiently close to heme c₁ to form a hydrogen bond to the propionate side chain of the heme (Zhang et al., 1998; Iwata et al., 1998). In another structure of the bc₁ complex from bovine heart from P6₅ crystals (Iwata et al., 1998), the ISP was observed to assume a third, "intermediate" position, in which the Fe₂S₂ cluster was located ~27.5 Å from cytochrome c₁ and 13 Å from the Q_o site, so that it would be impossible for His-161 to form a hydrogen bond to the quinone. The ability of the ISP soluble head to occupy different positions in the complex has been interpreted to imply that the quinol oxidation mechanism involves a substantial movement of the soluble head of the ISP between reaction domains in the cytochrome b and cytochrome c₁ subunits, while the transmembrane part of the ISP remains fixed, i.e., that when quinol binds to the Q_o site, the mobile ISP head moves into the proximal position, bringing the Fe₂S₂ cluster close to the Q_o site. After the Fe₂S₂ cluster is reduced by quinol, the ISP head moves into the distal position to reduce cytochrome c₁ (Zhang et al., 1998; Iwata et al., 1998; Kim et al., 1998; Crofts and co-workers, submitted for publication).

View larger version (47K): [in this window] [in a new window]   FIGURE 2   Rotation of the ISP mobile head in the bc1 complex upon binding of the inhibitor (stigmatellin). The cytochrome b, cytochrome c1, and ISP subunits are shown. (a) Arrangement of the bc1 complex subunits when stigmatellin is bound to the Qo site. (b) Arrangement of the bc1 complex subunits in the absence of the inhibitor.

One can investigate the movement of the ISP through molecular dynamics simulations. However, the expected time scale of the rotation (Crofts and co-workers, submitted for publication) is beyond the reach of molecular dynamics simulations, which are presently limited to time scales of a few nanoseconds for large proteins. Steered molecular dynamics (SMD) provides a means of overcoming this limitation by inducing the movement on the time scale accessible to molecular dynamics (Izrailev et al., 1997, 1998; Balsera et al., 1997; Isralewitz et al., 1997; Gullingsrud et al., 1999). In SMD simulations, time-dependent external forces are applied to a molecular system in order to accelerate the kinetics of the conformational change of interest.

In this paper we present a study of the mechanism of the ISP movement by SMD simulations on a time scale of 1 ns. For this purpose we have built the structure of a solvated bc₁ complex in the membrane bilayer. Modeling and simulating a complete solvated structure of the cytochrome bc₁ complex in a membrane bilayer posed extreme computational and methodological challenges due to the large size and complexity of the system. The structure of the bc₁ complex dimer alone comprises 61,180 atoms (including hydrogens). Building a bc₁ complex model suitable for the SMD simulations that included the solvent and a patch of membrane bilayer required manipulation of a few hundred thousand atoms and resulted in a system comprising 206,720 atoms, which was then reduced to a size of 91,061 atoms, of which 45,131 were actually subject to full dynamics.

The movement of the ISP between the configurations seen in the chicken complex can be described as a rotation by ~57° from the proximal to the distal position (Zhang et al., 1998). It was suggested that such a rotation could be accomplished through a one-dimensional constrained diffusion about the axis of rotation (Crofts and co-workers, submitted for publication). In this paper we test this suggestion by application of a torque to the soluble head of the ISP. The simulations reported below revealed that the conformations of the ISP observed in the crystal structures can indeed be connected through a smooth rotational path.

Simulations permitted an exploration of possible contacts of ISP with other subunits of the bc₁ complex along the rotational path. Analysis of these contacts, some of which were not present in the crystal structures, identified amino acid residues that may control the domain motion of ISP.

In the following section we report the point charge distributions calculated for different oxidation states of the redox centers of the enzyme, describe the building of the structure composed of the protein, the membrane bilayer, and the solvent, and define a 91,061 atom segment of the complete structure included in the SMD simulation. A comparison of the structure of the equilibrated model with the crystal structure, as well as the details of the trajectory followed by ISP, are presented in the Results section. The Discussion section summarizes the methodological advance achieved in this study, as well as our conclusions about the ISP motion and its functional relevance.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Movement of the rigid domains of the iron-sulfur protein

The ISP domain movement suggested by the crystal structures of the chicken cytochrome bc₁ complex was analyzed using the Hingefind algorithm (Wriggers and Schulten, 1997) implemented in a Tcl (Ousterhout, 1994) plug-in to the molecular visualization and analysis program VMD (Humphrey et al., 1996; Dalke and Schulten, 1997). The algorithm compares two known conformations of a protein and identifies connected regions that exhibit preserved packing within a specified tolerance of positional fluctuations. After dividing the protein into these relatively rigid domains, the algorithm determines effective rotation axes that characterize the movements of the identified domains relative to each other. The analysis was performed based on the coordinates of the -carbon atoms of the ISP from the structure of the native bc₁ complex and the structure with stigmatellin bound in the Q_o site of the complex (Zhang et al., 1998). The comparison of the two structures is illustrated in Fig. 3. The ISP movement was identified as a nearly rigid-body rotation of the water-soluble part of the ISP (comprising residues 73-196) by 56.11° with respect to the transmembrane helix of ISP (comprising residues 1-67). This result defined the exact axis along which a torque needed to be applied in an SMD simulation.

View larger version (35K): [in this window] [in a new window]   FIGURE 3   Comparison of the iron-sulfur protein backbone from the structure of the native bc1 complex (thick gray tube) and the structure with stigmatellin bound in the Qo site of the complex (thin black tube), done with the Hingefind algorithm (Wriggers and Schulten, 1997). The "hinge" region (residues 68-72) is shown as a thin gray line. The water-soluble parts of the protein (residues 73-196) are aligned; the transmembrane parts are separated by an angle of 56.11°. The arrow represents the effective axis of left-handed rotation from the initial to the final structure. In the protein the transmembrane helix is actually fixed and the soluble domain rotates (cf. Fig. 2).

Molecular dynamics

Molecular dynamics simulations reported in this paper were performed with the parallel molecular dynamics program NAMD (Nelson et al., 1996) using the CHARMM22 force field (MacKerell, Jr., et al., 1992, 1995, 1998; Schlenkrich et al., 1996). In all simulations a dielectric constant  = 1, and a cutoff of Coulomb forces with a switching function starting at 12 Å and reaching zero at a distance of 14 Å were assumed. All atoms, including hydrogens, were described explicitly. An integration time step of 1 fs was employed. All simulations were performed with the non-hydrogen atoms coupled to a thermal bath at 310 K, i.e., at a temperature 6.5 K above the gel  liquid crystal phase transition temperature of the dilauroyl-phosphatidylethanolamin (DLPE) lipid membrane (Blume, 1983) included in the model of the bc₁ complex (see below). The energy minimization was performed with the Powell minimization (Powell, 1977) feature of the program X-PLOR (Brünger, 1992) using the same set of parameters. The distribution of the point charges on the heme groups and the Fe₂S₂ cluster were calculated as described below.

Parametrization of the heme groups

There are three heme groups in the cytochrome bc₁ complex, which can be in either an oxidized or a reduced state, corresponding to the Fe³⁺ and Fe²⁺ oxidation states of the iron. In our simulations, all three hemes were assumed to be in the oxidized state and carry the same point charge distribution. The two carboxyl chains of the hemes were assumed to be deprotonated and carry a net charge of 2. The heme in the oxidized state carriers, therefore, a total charge of 1, while in the reduced state the total charge is 2. A charge distribution for the heme atoms is currently available in the CHARMM22 force field only for the reduced state, so for the purposes of this work the Merz-Kollman distributions of point charges (Singh and Kollman, 1984; Besler et al., 1990) on the oxidized and the reduced heme groups were calculated using the program GAUSSIAN-94 (Frisch et al., 1995) at the Hartree-Fock level with a 6-311G basis set. The coordinates of the atoms of a perfectly planar heme with deprotonated carboxyl chains were generated with the program QUANTA (MSI, 1994). A radius of 1.22 Å (Sutton, 1965) was used for the iron atom in the electrostatic potential fitting calculations.

The calculations of the charge distribution on the deprotonated heme in the oxidized state resulted in a highly nonuniform distribution with a large positive charge placed on one of the carboxyl oxygens. To achieve a more uniform charge distribution, the hydrogens protonating the carboxyl groups were added to the structure using scripting features of the molecular visualization and analysis program VMD (Humphrey et al., 1996; Dalke and Schulten, 1997), and the calculations were repeated. The four resulting charge distributions on reduced and oxidized, protonated, and deprotonated, hemes were symmetrized by averaging the charges on atoms in similar functional groups of the heme. For example, the charges obtained for all hydrogens of the methyl groups were averaged and assigned to these hydrogens. A diagram of the heme chemical structure is presented in Fig. 4.

View larger version (27K): [in this window] [in a new window]   FIGURE 4   Chemical structure of the heme group. Every atom is assigned a unique name. The first letter of the name corresponds to the chemical element of the atom.

Comparison of the resulting charge distribution on several groups of atoms for protonated oxidized heme to the corresponding distributions obtained for deprotonated and protonated hemes in the reduced state allowed us to identify the changes in the charge distribution of the protonated oxidized heme due to protonation and oxidation. If the total charge on a group did not change due to protonation, but changed as a function of the oxidation state, then the charge distribution from the oxidized protonated heme was adopted for the corresponding group of the oxidized deprotonated heme; if the total charge of a group did not change as a function of the oxidation state, but changed due to protonation, then the distribution from the reduced (deprotonated) heme was used. For example, the total average charge on the methyl groups of the reduced deprotonated, reduced protonated, and oxidized protonated hemes were 0.054161, 0.042947, and +0.002509, respectively. The first two numbers are close to each other and quite different from the third, so it can be concluded that the total charge of the methyl groups changes primarily due to oxidation rather than due to protonation. Therefore, the charge distribution with a total charge of +0.002509 was adopted for the methyl groups of the oxidized deprotonated heme.

This approach yielded a charge distribution for the oxidized deprotonated heme with a total charge of 1.0362. To correct the discrepancy from the expected value of 1, a charge of +0.0362 was added to the overall charge distribution. Since the (CH₂)₂COO groups of the oxidized heme can be expected to be slightly more positive than those of the reduced heme, the extra +0.0362 charge was distributed equally on these groups. Based on the differences in the charge distributions on the (CH₂)₂COOH groups of the oxidized and the reduced protonated hemes, a charge of +0.01 was added to the charge on the atoms CAA and CAD, and a charge of +0.0081 to the charge on the atoms CBA and CBD. The final charge distributions employed in the MD simulations are presented in Table 1 (see Fig. 4 for the atom name assignment).

It should be noted that the procedure described above is crude, but adequate for the purposes of MD simulations presented in this work. In cases where a precise description of the heme structure is required, more sophisticated methods such as geometry optimization should be used, and more structural information, e.g., the coordinates of the histidines coordinating the iron, should be included in the calculation of the charge distributions.

Parametrization of the Fe₂S₂ cluster

The Rieske Fe₂S₂ cluster in the bc₁ complex is coordinated by two histidine and two cystein residues, namely, by His-141, His-161, Cys-139, and Cys-158 of the ISP. When the Fe₂S₂ cluster is oxidized, both irons are in the Fe³⁺ oxidation state, whereas in the reduced cluster one iron is in the Fe²⁺ state. The overall charge of the cluster is 0 and 1 for the oxidized and reduced states, respectively. To calculate the charge distribution on the Fe₂S₂ cluster and the coordinating histidines and cysteins, the atomic coordinates from the 1.5 Å resolution structure of the water-soluble fragment of the ISP from bovine heart bc₁ complex (Iwata et al., 1996) were used [entry 1RIE in the Brookhaven Protein Data Bank (Bernstein et al., 1977)].

The coordinates of the hydrogen atoms of the ISP were generated using the HBUILD routine of X-PLOR (Brünger, 1992). For the purpose of generation of the hydrogen atoms, the irons and the sulfurs of the iron-sulfur cluster were modeled with point charges of +0.7 and 0.7, respectively. The coordinates of His-141, His-161, Cys-139, and Cys-158 were extracted from the resulting structure and used to model these residues as two 4-methylimidazole and two methylthiolate molecules to be included in the calculation of the charge distribution. The final configuration of the Fe₂S₂ cluster and of these molecules is represented schematically in Fig. 5.

View larger version (14K): [in this window] [in a new window]   FIGURE 5   Chemical structure of the Fe2S2 cluster and the side groups coordinating it. This structure was used in the calculations of the charge distribution. Atoms of every residue are assigned unique names. The first letter of the name corresponds to the chemical element of the atom.

The Merz-Kollman charge distribution for the reduced and oxidized states of the Fe₂S₂ cluster was calculated with a method similar to that used for the heme groups (see above) using the program GAUSSIAN-94 at the Hartree-Fock level with a 3-21G basis set. The obtained distributions were symmetrized by averaging the charges on identical atoms of the 4-methyl- imidazole and methylthiolate molecules, and on the sulfurs of the Fe₂S₂ cluster. The charges on the methyl group HB3 atoms, which are not present in the protein, were added to the charges on the corresponding CB atoms, in accordance with the method used to derive the charge distributions for histidine and cystein from the charge distributions of 4-methylimidazole and methylthiolate in the CHARMM22 force field. The final charge distributions used in the simulations are presented in Table 2 (see Fig. 5 for the atom name assignment).

The equilibrium bond length, angles, and torsional angles for the Fe₂S₂ cluster were derived from the coordinates of the water-soluble fragment of ISP (Iwata et al., 1996). The missing force constants for the Fe₂S₂ cluster and the His and Cys groups coordinating it were derived from existing parameters used for similar chemical structures in the CHARMM22 and CHARMM19 force fields, in particular, from the parameters available in CHARMM19 for prosthetic groups of the photosynthetic reaction center of Rhodopseudomonas viridis (Treutlein et al., 1988) and for the Met-liganded heme.

Preparation of the protein structure

The atomic coordinates of the cytochrome bc₁ complex from chicken heart mitochondria (Zhang et al., 1998) after refinement (Crofts and co-workers, submitted for publication) (entries 1BCC and 3BCC in the Brookhaven Protein Data Bank) were used. The coordinates of the hydrogen atoms for all protein subunits were generated with the HBUILD feature of X-PLOR. The structure of the bc₁ complex with stigmatellin bound at the Q_o site was used for building a model corresponding to the state of the complex in which the reduced Fe₂S₂ cluster is in the proximal position, and hemes b_L, b_H, and c₁ are oxidized by assigning appropriate charge distributions to the Fe₂S₂ cluster and the hemes, and removing stigmatellin from the model. Because the processes considered in this study take place in the intermembrane region of the complex, the protein domains located on the matrix side of the transmembrane region were not included in the simulations. The resulting structure, comprising 32,310 atoms, was refined by 500 steps of energy minimization of hydrogen atoms followed by 1500 steps of energy minimization of all atoms.

Placement of internal waters

The internal water molecules were placed in the protein following a procedure slightly modified from that implemented in the program package DOWSER (Zhang and Hermans, 1996). One of the initial steps of the procedure involves generating the coordinates of water oxygens on the molecular surface of the protein using the molecular surface calculation program MS (Connolly, 1983). Due to the large size of the bc₁ complex, the surface calculation for the whole bc₁ dimer failed, so the molecular surface was calculated separately for each of the two monomers. By using DOWSER, the candidate positions of the water molecules obtained were combined, sparsified, and classified as located on the solvent-exposed protein surface or buried inside the protein; 366 buried water molecules were identified by DOWSER as having low energy and separated from the rest of the system. Since in the calculations performed by DOWSER the non-polar hydrogen atoms, heme groups, and the Fe₂S₂ cluster were not taken into account, the hydrogens of the selected water molecules were rebuilt in the protein environment with full-atom representation using the HBUILD routine of X-PLOR. After an additional 500 steps of energy minimization of the water molecules within the fixed protein structure, 121 water molecules with interaction energies with the rest of the system smaller than 12 kcal/mol (a value used by DOWSER) were selected, and the rest of the water molecules removed. Finally, the remaining water molecules were subjected to 500 more steps of energy minimization while the rest of the system was fixed. This yielded a structure of the bc₁ complex with internal water molecules, which was used in further modeling.

Preparation of the lipid bilayer

To model the cytochrome bc₁ complex in a membrane bilayer, a membrane patch sufficiently large to be suitable for insertion of the complete dimer was generated. For this purpose, the structure of a solvated DLPE lipid bilayer modeled previously (Zhou and Schulten, 1995) was chosen as a starting point. This structure contained 101 lipids in each monolayer with a surface area of ~50 Ų/lipid, and 8108 water molecules. The bilayer was oriented with the membrane surface parallel to the xy plane and shaped as a cylinder of a 43 Å radius and 71.5 Å height, aligned with the z-axis. In order to generate a larger membrane patch, a 60 Å × 60 Å patch of lipids and water, composed of 72 lipids in each monolayer and 4847 water molecules, was cut out of the cylindrical patch. Since the structure did not contain polar hydrogens, they were generated using the program X-PLOR. The resulting structure comprising ~28,500 atoms was subjected to 500 steps of energy minimization with periodic boundary conditions to achieve periodicity of the patch. The system was then equilibrated for 100 ps with NAMD running on 36 300-MHz processors of a Cray T3E, which required 0.19 h/ps of simulation time.

The obtained equilibrated periodic structure of the DLPE bilayer was used to generate a 120 Å × 155 Å patch of membrane by multiplication of the 60 Å × 60 Å patch, resulting in a structure which consisted of 372 lipids in each monolayer and 25,162 water molecules, a total of ~147,500 atoms. The new patch was centered at the origin, with the longer side of the patch aligned with the y-axis. To restore periodicity in the y-direction all atoms of this structure, except those located within 10 Å of the faces of the patch perpendicular to the y-axis, were fixed and the structure was minimized for 500 steps with periodic boundary conditions.

Preparation of the solvent

Water for the solvation was generated from a 45 Å sphere of water molecules prepared as described previously (Bishop et al., 1997; Kosztin et al., 1997). First, a 50 Å × 50 Å × 35 Å box containing 3026 water molecules was selected from the sphere and subjected to 500 steps of energy minimization with periodic boundary conditions. The obtained periodic box of water molecules was further equilibrated for 100 ps. A 120 Å × 155 Å × 35 Å box of water molecules, comprising ~67,500 atoms, was generated from the resulting structure by duplication, and the periodicity was achieved by energy minimization of a 3 Å layer of water molecules on the borders of the box.

Solvation and placement of the protein into the lipid bilayer

Positioning of the protein in the membrane bilayer and solvent was carried out according to the following procedure. The x- and y-axes were chosen along the 120 Å and 155 Å sides of the created membrane patch, respectively, and the z-axis was aligned with the symmetry axis of the bc₁ complex dimer. The protein was then positioned in the center of the 120 Å × 155 Å rectangle in the x- and y-directions, and rotated about the z-axis to fit into the rectangle.

The center of the membrane bilayer in the z-direction was determined by calculation of the distributions of the positions of the lipid phosphate groups, and of the ends of the lipid hydrophobic tails. The latter distribution featured one peak centered in the middle between the two peaks of the distribution of the phosphate groups. The distribution of positions of the tryptophan residues along the z-axis in the protein exhibited two pronounced peaks located ~22 Å apart. The center between these two peaks was interpreted as the center of the transmembrane region of the bc₁ complex. The solvated membrane bilayer was then positioned along the z-axis so that the centers of the bilayer and of the transmembrane region of the protein were located at the same point. To completely cover the intermembrane region of the protein with solvent, the 35 Å layer of solvent water molecules prepared as described above was added to the system. This layer was positioned next to the boundary of the water molecules solvating the bilayer. Fig. 6 illustrates the relative position of the protein with respect to the membrane and to the layer of solvent water molecules.

View larger version (90K): [in this window] [in a new window]   FIGURE 6   Model of the structure of the cytochrome bc1 complex from chicken heart mitochondria with the iron-sulfur protein in the position proximal to the Qo site, inserted in the membrane bilayer. The core proteins and other domains of the bc1 complex on the matrix side of the membrane are not included. The z-axis is pointing vertically up. The protein backbone is represented by gray tubes; the Fe2S2 cluster and the heme groups are rendered through black spheres; the lipid molecules are drawn with lines, and the solvent water molecules are shown as dots. The shading of the solvent reflects the procedure of building the model. Water molecules that were used for solvation of the lipid bilayer are shown in light gray, while the water molecules added to the system to completely solvate the intermembrane domains of the bc1 complex after placing it into the membrane are shown in dark gray. The labels indicate naming conventions for three layers of water, as used in the text.

All water molecules that had the oxygen closer than 2.4 Å to the bc₁ complex were removed. Those of the remaining solvent water molecules, which were buried in the protein, were identified using DOWSER and also removed, so that they did not interfere with the internal water molecules placed earlier. Lipid molecules that had non-hydrogen atoms closer than 2.4 Å to the protein were separated and visually examined using VMD (Humphrey et al., 1996; Dalke and Schulten, 1997). Lipids that had a significant overlap with the protein were removed from the system, while those that had no or little contact with the protein were retained. By using several sets of 100 steps of energy minimization with periodic boundary conditions imposed in the x- and y-directions, manual editing of the lipid conformations using VMD, and removal of some of the lipid molecules, all unfavorable contacts between the protein and the lipids were eliminated. During the minimization all protein atoms and all water molecules, except for a 10 Å layer of water molecules on each side of the membrane, were fixed. Finally, a 6 Å layer of water molecules on the boundary between the 35 Å layer of solvent water (water I in Fig. 6) and the layer of water molecules solvating the membrane (water II in Fig. 6) was minimized for 100 steps, while the rest of the system was held fixed.

With exclusion of the core proteins and of other bc₁ complex domains located on the matrix side of the membrane, the resulting structure contained the bc₁ complex, 121 buried water molecules, 10,446 solvent water molecules on the matrix side of the membrane, 19,873 solvent water molecules on the intermembrane side, 288 lipid molecules forming the matrix side, and 283 lipid molecules forming the intermembrane side of the membrane, with a total of 206,720 atoms (see Fig. 6). The size of the system needed to be further reduced in order to carry out the 1 ns SMD simulation required for the present study.

Building the structure for SMD simulations

Since the focus of this investigation is the motion of the ISP domain of the bc₁ complex, we found it appropriate to further reduce the system size and only include in the simulation those parts of the system that are likely to affect the motion of the domain. The ISP of monomer I, which is involved in the redox reactions with the cytochromes b and c₁ of monomer II of the bc₁ complex, corresponding to chain E in the notation used in the 3BCC structure from the Protein Data Bank, was chosen as the subject of application of external forces in SMD simulations. Accordingly, the protein residues containing atoms within the electrostatic cutoff of 14 Å from the ISP of monomer I were identified. This set of residues was used to select the parts of the protein that were allowed to move in the simulations, while the rest of the protein was held fixed. These mobile parts included half of the transmembrane region, and all the soluble segments of cytochrome b, cytochrome c₁, subunit 7, and subunit 8 from monomer II that were located in the intermembrane region of the bc₁ complex, as well as residues of cytochromes b and c₁ from monomer I located on the interface between monomers. Thus, the mobile parts of the bc₁ complex consisted of residues 41-92, 122-188, 237-298, 332-357, and heme b_L of cytochrome b from monomer II; residues 1-212 and heme c₁ of cytochrome c₁ from monomer II; residues 49-196 and the Fe₂S₂ cluster of ISP from monomer I; residues 47-80 and 163-188 of cytochrome b from monomer I; residues 43-107 of cytochrome b from monomer I; residues 53-79 of subunit 7 from monomer II; and all residues of subunit 8 from monomer II. The water molecules solvating the lipid monolayer forming the matrix side of the membrane (water III layer in Fig. 6) were removed from the system. The remaining lipid molecules and water molecules from the water II layer (see Fig. 6) that did not have atoms separated by <20 Å from the mobile parts of the protein were also removed. Finally, all solvent water molecules from the water I layer that did not have any atoms within 15 Å of the mobile parts of the protein were deleted from the system. The lipid and water molecules that did not have any atoms within 15 Å of the mobile parts of the protein and all the lipids forming the matrix side of the membrane were fixed during the simulations in order to maintain the structure of the lipid bilayer remaining in the system. The resulting structure, shown in Fig. 7, included the protein, solvated by a 15 Å layer of mobile water and an additional 5 Å layer of fixed water, and lipid molecules, with a total of 91,061 atoms, 45,131 of which were fixed during the simulations.

View larger version (64K): [in this window] [in a new window]   FIGURE 7   The final simulated structure of the bc1 complex. Cytochrome b, cytochrome c1, subunits 7 and 8 of monomer II, and the ISP from monomer I are presented as black spheres; the rest of the protein is presented through gray spheres; the lipid molecules are shown as gray lines, and water molecules are shown as dots. (a) A view with the z-axis pointing vertically up. Water molecules in front of the protein are omitted. (b) A view along the z-axis. Water molecules are not shown.

The structure of the bc₁ complex obtained was equilibrated for 100 ps. During the first 50 ps of equilibration, all protein backbone atoms were fixed. To prevent water molecules from escaping from the system, it was surrounded by repulsive harmonic walls defined as the surface of a cylinder with a radius of 85 Å, a height of 170 Å, and faces parallel to the xy plane.

Setup of steered MD simulation

To induce the rotation of the mobile head of ISP about its rotation axis, as identified by the Hingefind algorithm and described above, external forces generating a torque applied on the ISP domain were exerted in SMD simulations on the -carbon atoms of residues 73-196 of ISP. The forces were implemented by harmonically restraining these atoms to restraint points located at the initial positions of the atoms and then rotating the restraint points about the rotation axis with a constant angular velocity . This corresponds to attaching harmonic springs to the atoms and pulling the ends of the springs along circular trajectories about the rotation axis. The force _i applied to the atom i by the spring is

(1)

where K is the force constant of the spring _i(t) is the position of the restrained atom, _i(t) is the position of the restraint point of this atom, and _i(0) = _i(0). The torque N applied to a restrained atom i with respect to the rotation axis is then

(2)

where is a unit vector of the rotation axis, and _i is a vector normal to the axis, pointing from the atom to the axis.

When the system has to overcome a potential energy barrier, e.g., break a hydrogen bond between ISP and another subunit, the attached springs stretch, which results in an increase of the applied force, and, therefore, the applied torque. After the hydrogen bond is broken, the springs relax and the applied torque decreases. Consequently, the profile of the applied torque exhibits peaks that correspond to transitions between the local potential energy minima along the path of rotation. These peaks are characteristic for all SMD simulations and carry the most essential information revealed (Grubmüller et al., 1996; Izrailev et al., 1997; Isralewitz et al., 1997; Stepaniants et al., 1997; Lu et al., 1998; Hermans et al., 1998; Kosztin et al., 1999; Wriggers and Schulten, 1999).

The force constant of the restraints was chosen to be K = 1 kcal/mol Ų  70 pN/Å, which corresponded to thermal position fluctuations of the restrained atoms of x =  ~ 0.8 Å. To cover the span of the rotation angle  from 0 to 56.11° within 1 ns, an angular velocity  = 0.05611°/ps was chosen. The simulation was carried out with NAMD running on 64 450-MHz processors of a Cray T3E and required ~0.18 h/ps of simulation time.

The SMD trajectory was analyzed by calculating the torque applied to each of the restrained -carbon atoms with respect to the assumed axis of rotation and the actual angle of rotation of each of these atoms every 2 ps. The applied torque and the rotation angles were averaged over all 126 restrained residues, as well as analyzed for individual residues.

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Equilibration reveals water channel

During the equilibration of the solvated structure of the bc₁ complex from chicken heart mitochondria placed in the DLPE bilayer, with stigmatellin removed from the Q_o site, the positions of several cytochrome b residues in the Q_o site changed as shown in Fig. 8. In the equilibrated structure, Lys-288 that in the crystal structure was hydrogen-bonded to Ser-152 had moved to form a hydrogen bond with the backbone oxygen of His-141 of ISP. Ser-152 changed its orientation toward the inside of the Q_o site. Tyr-279 shifted inside the Q_o site by ~2.5 Å. The Fe₂S₂ cluster, as well as His-141 and His-161 coordinating it, shifted accordingly in order to maintain a hydrogen bond with the backbone oxygen of Cys-160 of ISP found in the crystal structure. In this position the ring of Tyr-279 occupied some of the space taken by the head of stigmatellin in the crystal structure. Pro-271 moved inside the Q_o site by ~2.2 Å. The hydroxyl of Tyr-274, the carboxylate of Glu-272, and the nearby heme propionate reoriented so as to connect through hydrogen bonds to internal waters forming a water channel to the aqueous phase, as shown in Fig. 9. This hydrogen-bonded network was relatively stable during the subsequent simulation. Phe-129 rotated by ~120° to occupy the space that was occupied by the tail of stigmatellin in the crystal structure. The loop and the helix of cytochrome b composed of residues 156-166 moved by ~3 Å toward the center of the membrane bilayer, which might have been caused by the contacts made with the lipid molecules. This displacement has not been seen in the crystal structures available, and is discussed further below. The root-mean-square deviation of the positions of the backbone atoms in the crystal structure and the equilibrated structure was 1.77 Å.

View larger version (29K): [in this window] [in a new window]   FIGURE 8   Changes observed in the structure of the Qo binding site after equilibration of the bc1 complex structure with stigmatellin removed. The positions of the protein residues in the crystal structure with bound stigmatellin are shown in gray; positions after equilibration are shown in black. Stigmatellin and the Fe2S2 cluster are presented through thin lines. Dashed lines represent hydrogen bonds.

View larger version (29K): [in this window] [in a new window]   FIGURE 9   A water channel connecting Glu-272 and Tyr-274 located in the lobe of the Qo site to the solvent. Asn-249 participates in the hydrogen bond network with buried water molecules forming the channel. Protein residues and the heme bL are shown in black; water molecules are shown in gray; location of the bulk solvent is indicated by a thin dashed line.

Rotation of ISP

A 56° rotational motion of the ISP mobile head was realized in our SMD simulations during a period of 1 ns. The angle of rotation averaged over all residues experiencing the applied torque (residues 73-196 of ISP) exhibited very little fluctuation and depended linearly on time. The torque applied to the protein with respect to the rotation axis, averaged over residues 73-196, as a function of time and of the average rotation angle, is presented in Fig. 10. The applied average torque featured five pronounced peaks at ~50, 175, 320, 600, and 900 ps. These peaks corresponded to distinct structural changes of the bc₁ complex in the course of the ISP rotation.

View larger version (40K): [in this window] [in a new window]   FIGURE 10   Torque applied to residues 73-196 of ISP averaged over these residues as a function of time (bottom axis) and the rotation angle (top axis). The thick line represents a running average of the data with a window of 20 points.

To analyze the relative contributions of the restrained residues to each of these peaks, the torques applied to individual residues were averaged over every five consecutive residues and over time periods of 50 ps, yielding a set of values presented as a color map in Fig. 11. The color of each square represents the average value of the torque applied to the residues with numbers corresponding to this square within the respective 50 ps time period. Darker colors denote larger values of the torque. It can be concluded from Fig. 11, for example, that between 150 ps and 300 ps residues 113-118, 138-163, and 173-178 of ISP experienced the largest values of the applied torque. Another region with large applied torque was observed between 450 ps and 650 ps for residues 138-163 of the ISP, which form the binding pocket of the Fe₂S₂ cluster. The time intervals with the darker regions in Fig. 11 correlate with the time intervals of the peaks of the applied torques in Fig. 10.

View larger version (90K): [in this window] [in a new window]   FIGURE 11   Torque applied to residues 73-196 of the ISP averaged over groups of five residues and over 50-ps time intervals. The color of each square represents the average value of the torque applied to the residues with numbers corresponding to this square within the respective 50-ps time period. Darker colors denote larger values of the torque.

Analysis of the torques applied to individual residues showed that the contributions of each residue to the average torque were unevenly distributed. For example, due to formation of hydrogen bonds between the residues of ISP close to the Fe₂S₂ center and residues of cytochrome b, the torque applied to His-161, shown in Fig. 12 a, was significantly larger than that applied to ISP residues, such as Arg-126, which did not interact with other subunits, as shown in Fig. 12 b. Two peaks of the torque applied to His-161 observed in Fig. 12 a contributed to the peaks of the average torque at ~320 ps and 600 ps.

View larger version (30K): [in this window] [in a new window]   FIGURE 12   (a) Torque applied to His-161 of ISP during the simulation. (b) Torque applied to Arg-126 of ISP during the simulation.

Effect of initial torque

Initially, i.e., at t = 0, the Fe₂S₂ cluster was located in the position proximal to the Q_o site of cytochrome b with hydrogen bonds formed between the backbone oxygen of His-141 of ISP and Lys-288 of cytochrome b, and between the backbone oxygen of Cys-260 of ISP and Tyr-279 of cytochrome b as shown in Fig. 13 a. As the torque was applied to the system, the strain applied to the -sheet region of ISP (residues 73-78 and 182-196) adjacent to the "hinge," induced the rupture of a hydrogen bond between the backbone atoms of Leu-78 and Asp-80. This rupture, as well as the breaking of a hydrogen bond between Lys-104 of ISP and Glu-51 of subunit 8, resulted in the first peak of the applied torque observed in Fig. 10 at ~50 ps.

View larger version (74K): [in this window] [in a new window]   FIGURE 13   Snapshots of the configuration of the Fe2S2 cluster (atoms shown as gray spheres) and key protein residues during the induced rotation of ISP. His-141, Cys-160, and His-161 of ISP, Pro-267, Tyr-279, and Lys-288 of cytochrome b, and the carboxyl groups of heme c1 are shown in licorice. The loop composed of residues 263-268 of cytochrome b is presented as a tube segment. (a) Initial configuration, t = 0; (b) the cytochrome b loop has been shifted by ISP, t = 200 ps; (c) hydrogen bond between Tyr-179 and His-161 has been formed, t = 340 ps; (d) ISP lost most of the contacts with cytochrome b, t = 600 ps; (e) hydrogen bond between His-161 and Pro-267 has been formed, t = 750 ps; (f) final configuration, t = 1 ns.

Displacement of a loop of cytochrome b

At ~150 ps, the ISP encountered a loop of the cytochrome b composed of residues 263-268, and shifted it by ~2 Å. The shift can be seen from the decrease of the distance between the loop and heme c₁ in Fig. 13 b as compared to that in Fig. 13 a. After the loop was shifted, the applied torque decreased, resulting in a second peak at ~200 ps (see Fig. 10). Three hydrogen bonds formed by the ISP with the cytochrome c₁ subunits, which were broken at the same time, also contributed to the decrease of the applied torque: 1) a weak water-mediated hydrogen bond between Lys-90 of ISP and Arg-102 of cytochrome c₁ from monomer II; 2) a water-mediated hydrogen bond between Lys-85 of ISP and Arg-144 of cytochrome c₁ from monomer II; 3) a bond between the backbone oxygen of Ala-70 of the ISP and Lys-86 of cytochrome c₁ from monomer I. The bonds 1-3 were ruptured at 172, 174, and 192 ps, respectively.

Detachment of Cys-160 from Tyr-279 of cytochrome b

The largest applied torque was observed at ~320 ps, when the rotation reached ~18° (see (Fig. 10). At this time, the loops of ISP holding the Fe₂S₂ cluster, i.e., residues 139-148, 155-164, 172-178, and an adjacent region composed of residues 113-121 of ISP, experienced the majority of the torque applied to the system, as can be discerned from Fig. 11. During the time between 280 and 325 ps, several hydrogen bonds between the ISP and the other subunits were broken. The largest contribution to the torque was due to the need to rupture a hydrogen bond between Tyr-279 of the cytochrome b and the backbone oxygen of Cys-160 of ISP. After the ISP continued its rotation and the torque decreased, this bond was temporarily compensated by formation of a hydrogen bond between Tyr-279 of cytochrome b and His-161 of ISP, coordinating the Fe₂S₂ cluster as shown in Fig. 13 c. The latter hydrogen bond was subsequently ruptured at ~364 ps. Another contribution to the decrease of the applied torque came from the rupture of the hydrogen bonds between the backbone oxygen of Leu-142 of ISP and Asn-149 of cytochrome b, and between Lys-90 of ISP and Glu-76 of cytochrome c₁ from monomer I.

Formation of a hydrogen bond with Pro-267 of cytochrome b

During the time interval between 450 ps and 670 ps, the ISP rotated sufficiently to lose all its contacts to cytochrome b except the contacts made to the cytochrome b loop composed of residues 263-268 (see Fig. 13 d). In the conformation assumed by the ISP at ~600 ps these contacts, in particular, the contacts between the loop and the atoms forming a disulfide bond between Cys-144 and Cys-160 of ISP, were very similar to those observed in the crystal structure of the bc₁ complex from beef with ISP in the "intermediate" position (Iwata et al., 1998). The loss of these contacts resulted in a broad peak of the applied torque centered at ~600 ps. The largest contribution to the applied torque was again made by the residues surrounding the Fe₂S₂ cluster. The rotation was accompanied by breaking of a water-mediated hydrogen bond between His-141 of ISP and Lys-288 of cytochrome b, and the hydrogen bonds formed by Glu-108 and Asp-152 of ISP of Glu-49 and Glu-49 of subunit 8, respectively. At ~670 ps, His-161 of ISP formed a hydrogen bond with the backbone oxygen of Pro-267 of cytochrome b, as shown in Fig. 13 e. Breaking of this bond (see Fig. 13 f) resulted in another peak of the applied torque at ~916 ps.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

In the simulations presented here, SMD applications were extended from inducing an unbinding of ligands from proteins (Grubmüller et al., 1996; Izrailev et al., 1997; Isralewitz et al., 1997; Lüdemann et al., 1997; Hermans et al., 1998; Wriggers and Schulten, 1999; Kosztin et al., 1999) to inducing a domain motion in a very large aggregate of proteins, lipids, and water. Domain motion occurs during the functional cycle of many multidomain protein complexes, e.g., in ATPase (Abrahams et al., 1994; Elston et al., 1997; Wang and Oster, 1998). SMD opens the possibility of investigating such motions on time scales accessible to molecular dynamics simulations.

The simulations described in this paper met new challenges arising in the analysis of the SMD results due to the complexity and multidimensionality of the simulated system. Protein domain movement was induced by applying a torque to 126 atoms simultaneously. The behavior of the total torque applied to the protein revealed general features of the ISP rotation process; however, analysis of the contributions of individual residues to the total torque was required to discern the detailed mechanism of the rotation.

For each restrained residue, a force constant corresponding to a soft spring, K  70 pN/Å, was employed (Izrailev et al., 1997; Balsera et al., 1997). The applied torque exhibited profiles (see Fig. 12) typical for soft springs, but since application of n springs of stiffness K to a rigid system is equivalent to application of one spring of stiffness nK, the profile of the torque averaged over all restrained residues featured large fluctuations typical for stiff springs (Izrailev et al., 1997). To reduce the magnitude of the fluctuations of the total applied torque, and thus increase the level of detail obtained about the simulated process from the torque profile, it is desirable to use even smaller force constants in SMD simulations when external forces are exerted on many atoms. However, in a 100-ps test simulation of the ISP rotation that employed a smaller force constant K  5 pN/Å (results not shown), the restrained atoms of the ISP lagged significantly behind the restraint points, so that the forces applied to the atoms were not directed along the tangent of the assumed rotation path and did not produce an appropriate torque.

The results of the SMD simulations presented in the previous section reinforce the suggestion that a motion of the soluble head of ISP controls the functional cycle of the cytochrome bc₁ complex. It was shown that the motion can be described by a continuous rotational path connecting two conformations of ISP observed in the crystal structures of the bc₁ complex from chicken heart mitochondria. The fact that in order to follow the rotational path ISP had to push a loop of cytochrome b composed of residues 263-268 from its original position by ~2 Å indicated that the motion of ISP in natural conditions may not comprise a rotation about only one axis, but may be rather more complex. Since a complete detachment of ISP from cytochrome b appears to involve a potential energy barrier, reflected by the peak of the torque applied to the ISP at ~600 ps (see Fig. 10), the ISP is likely to remain in close contact with cytochrome b. Possible deviations from a purely rotational motion should not affect the qualitative picture of the contacts made by the ISP to the other subunits, as captured in the SMD simulations.

The presence of several peaks of the torque applied to the system reflects the possibility of several potential energy barriers that the ISP has to surmount during its motion. At present it is impossible, however, to determine quantitative features of these barriers from SMD data due to the short (1 ns) time scale of the simulation and the lack of appropriate analysis methods. The peak at 50 ps resulted mainly from the rupture of an inner ISP backbone hydrogen bond. This event may not happen on the time scale of the natural movement of the ISP, so this peak may to a certain degree be an artefact of the short time scale of the simulation compared to that of the natural functional cycle of cytochrome bc₁ complex.

The magnitude of the average torque applied to ISP residues during simulation is presented in Fig. 11. Multiplication of the average torque by the number of residues to which the torque was applied (123) yields a total torque on the order of 12,300 pN Å. This value of the torque is 30-50 times larger that the torques measured for the ATPase (230-450 pN Å) at the rotation rates ~17 revolutions per second (Noji et al., 1997). In our simulations, the 56° rotation of ISP was induced within 1 ns, which corresponds to a rotation rate of 155 revolutions per microsecond, i.e., seven orders of magnitude faster. To induce this extremely fast rotation, one needs to overcome the friction forces arising from the movement against the viscosity of the water. These friction forces account for the large applied torque.

An estimate of the energy required to accomplish the rotation can be obtained by the amount of mechanical work performed by the torque. Integration of the time-averaged torque presented in Fig. 10 over the rotation angle yields the work performed by this torque of ~98 kcal/mol. Due to the nonequilibrium nature of the simulation, most of this work is irreversible and is transformed into thermal energy. However, this energy is much smaller than the total thermal energy of the system (on the order of 42,000 kcal/mol at 310 K) and, therefore, has very little effect on the dynamics.

The contacts observed between ISP and the other subunits of the bc₁ complex along the rotational path identified a possible mechanism of guiding the motion of ISP from one stable conformation observed in the crystal structures to the other through a series of metastable conformations, transitions between which correspond to the observed peaks of the applied torque. Such a correspondence resembles the results of SMD simulations of protein-ligand complexes (Grubmüller et al., 1996; Izrailev et al., 1997; Isralewitz et al., 1997; Wriggers and Schulten, 1999; Kosztin et al., 1999) where an external harmonic force was applied to a ligand in a chosen direction in order to induce the unbinding of the ligand, and the maxima of the applied force reflected, e.g., the breaking of hydrogen bond networks formed between the protein and the ligand.

Formation of hydrogen bonds by Tyr-279 of cytochrome b first to Cys-160, and then to His-161 of ISP, as well as subsequent formation of a hydrogen bond between His-161 of ISP and Pro-267 of cytochrome b in the course of the simulation lead to the suggestion that Tyr-279 and Pro-267 of cytochrome b guide the ISP into its docking positions near cytochrome b and cytochrome c₁. Pro-267 is not conserved among the sequences of cytochrome b from different species, so the hydrogen bond formed with this residue by His-161 of ISP in the bc₁ complex from chicken heart mitochondria considered in this study could be replaced by a hydrogen bond with another residue of the loop ef of cytochrome b during the motion of ISP in bc₁ complexes from other species. The importance of this loop in the mechanism is highlighted by the fact that the mutations of Thr-265 to Met (Brandt, 1998) and Ile-269 to Phe (Crofts et al., 1995) severely interfere with quinol oxidation, although these residues do not impinge on the volume of the Q_o site.

The simulations also indicated that ISP makes hydrogen bond contacts with other subunits of the bc₁ complex, namely with Glu-76 and Lys-86 of cytochrome c₁ of monomer I, Arg-102 and Arg-144 of the cytochrome c₁ of monomer II, and Glu-59 and Glu-51 of subunit 8 of monomer II. Breaking of these contacts contributed to the observed peaks of the applied torque. The mentioned residues are, therefore, also likely to stabilize the mobile head of ISP in intermediate positions during its motion.

We would like to emphasize that a certain degree of caution should be used in the interpretation of the results of SMD simulations. While the simulations suggested possible structural changes of the bc₁ complex, some of these changes may reflect the limitations of the model employed.

The displacement of the ISP was completed in 1 ns, compared with a natural displacement that may take as long as 100 µs (Crofts and co-workers, submitted for publication), and therefore the system was not able to equilibrate at each step of the movement. In addition, the displacement was forced by application of a rotational torque, which constrained the movement within limits that would not apply in the natural system. These limitations introduced some artefacts into the simulations, e.g., some displacements of side chains and flexible loops are likely exaggerated compared to the natural system. In particular, the following events might be unnatural. 1) Some of the contacts of the ISP with cytochrome c₁ observed in the simulation resulted from the movement of flexible loops of cytochrome c₁ by ~3 Å. This movement was caused in part by the pressure from water molecules that did not have enough time to escape the space between cytochrome c₁ and the ISP during rotation. In a longer simulation, water molecules would diffuse away without affecting cytochrome c₁. 2) The collision of the ISP with the loop ef of cytochrome b (residues around Leu-263) would likely be less dramatic, and softened by lateral diffusion in an unconstrained motion. 3) The constraints introduced by our method of applying torque effectively restricted the motions within ISP to flexibility about the fixed configuration represented by the refined structure based on the ISP soluble fragment (Iwata et al., 1996). This might have prevented the simulation from exploring alternative configurations, such as the "intermediate" state seen in P6₅ crystals from the recent structure of the complete bc₁ complex from beef (Iwata et al., 1998). The constraints were applied to residues 73-196, and did not therefore apply to the hinge span (see below). 4) In the crystallographic structure of the native complex, the hinge region of the ISP shows a helical configuration, which is not found at the end of the simulation. We ascribe this to the short time scale of the simulation that was insufficient for the helix to form. 5) Several other differences between the native structure and the structure arising at the end of the SMD simulation also reflect the nonequilibrium nature of the simulations.

A second limitation relates to a recognition that our attempts to place the protein in a membrane must inevitably include some short-cuts. Thus, the lipid composition of the membrane used in the simulations (DLPE) is certainly different from that of a natural membrane. Unfortunately, it is not clear what the effects of the lipid composition are on the protein structure, since although it is known that the activity of the isolated complex can be modified by specific phospholipids, it is not known how. Some displacements from the crystallographic positions observed during the preliminary equilibration after removal of stigmatellin from the structure may reflect this inadequacy. In particular, the 3.0 Å displacement of the loop of cytochrome b containing residues 156-166 is not seen in any of the crystals, does involve interaction with the lipid membrane, and may be artificially induced by these interactions.

Other limitations relate to the adequacy of our treatment of the point charge distributions on the heme groups and the Fe₂S₂ cluster, which were assumed to maintain their oxidation states during the simulation. The heme groups are fairly well removed from the ISP, and we do not anticipate that the trajectory seen would be much affected by a more refined treatment. However, many of the events occurring during the simulated movement are the breaking and making of hydrogen bonds with residues coordinating the Fe₂S₂ cluster. The relative strength of these bonds will clearly be affected by the charge distribution on the cluster and on these residues. Since the simulation involved only the reduced form of the complex, we need not be concerned with changes in hydrogen bonding associated with pK changes accompanying the changes of the redox state of the Fe₂S₂ cluster. Future simulations in which the redox state of the cluster is different will need to address this issue more critically.

A concern may arise due to our choice of the crystal structure used as the initial model in the simulations. The atomic coordinates of the cytochrome bc₁ complex from chicken heart mitochondria reported by Zhang et al. (1998) were unrefined and the electron density of the ISP was originally fit using the coordinates of the ISP soluble fragment from beef enzyme (Iwata et al., 1996). However, the structure of the stigmatellin-containing crystals had the ISP at close to 100% occupancy and the data were therefore quite well-defined, and could be well-fitted by the model. In our simulations we used the coordinates from a later data set. The coordinates for the ISP were refined after the initial fitting, so that this data set corresponds to an independently refined model.

The stigmatellin-containing structure from chicken (Zhang et al., 1998) is the only one with the ISP in the proximal position that does not have crystallographic contacts involving the ISP mobile domain and provides a reasonable starting point for the SMD simulations presented in this paper. Although the sequence used for ISP in the chicken structure was that of the beef enzyme, the interfacial surface that interacts with the cytochrome b is composed of completely conserved residues when alignments are made including bacteria species as well as known mitochondrial sequences. It therefore seems highly unlikely that the structure of this interfacial region is different in chicken.

One could also argue that the simulation trajectory should explicitly take into account the "intermediate" position of ISP seen in the P6₅ structure from beef (Iwata et al., 1998) and consists of two segments; one rotating the ISP head from the proximal position to the "intermediate" position, and another from the "intermediate" position to the distal position. However, since the ISP head has been found in six different positions in different native crystals of mitochondrial complexes (including four different positions in beef), it is apparent that these positions are more an accident of crystallography than an indication of physiological significance. Analysis of the chicken mitochondrial complex has shown that a substantial fraction of the ISP head is unaccounted for by electron density in all native structures from chicken (Zhang et al., 1998), and the same is true in the analysis of the structures from beef (Kim et al., 1998). The structures are likely disordered because of weak occupancy of a number of other states. Therefore, there is no justification for supposing that the "intermediate" position is more physiologically relevant than any of the other positions observed. This view is also supported by the fact that the ISP head in the "intermediate" position is seen only in one monomer of the bc₁ complex and is constrained by crystal contacts, while the ISP head in the other monomer in the same crystals is in a disordered state. Because we do not intend to demonstrate that the pathway mapped in the simulations reported here is the only, or even a preferred trajectory, forcing the ISP head to visit the "intermediate" position, i.e., only one position among an infinity of other possible states, would not help to answer the