This Article Abstract Full Text (PDF) Supplemental Material Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Li, G. Articles by Cui, Q. Search for Related Content PubMed PubMed Citation Articles by Li, G. Articles by Cui, Q.

Analysis of Functional Motions in Brownian Molecular Machines with an Efficient Block Normal Mode Approach: Myosin-II and Ca²⁺-ATPase

Department of Chemistry and Theoretical Chemistry Institute, University of Wisconsin, Madison, Wisconsin

Correspondence: Address reprint requests to Qiang Cui, Dept. of Chemistry, University of Wisconsin at Madison, 1101 University Ave., Madison, WI 53706. Tel.: 608-262-9801; E-mail: cui{at}chem.wisc.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION COMPUTATIONAL METHODS STRUCTURAL FLEXIBILITY AND... DISCUSSION CONCLUSIONS SUPPORTING MATERIALS ACKNOWLEDGEMENTS REFERENCES

 
The structural flexibilities of two molecular machines, myosin and Ca²⁺-ATPase, have been analyzed with normal mode analysis and discussed in the context of their energy conversion functions. The normal mode analysis with physical intermolecular interactions was made possible by an improved implementation of the block normal mode (BNM) approach. The BNM results clearly illustrated that the large-scale conformational transitions implicated in the functional cycles of the two motor systems can be largely captured with a small number of low-frequency normal modes. Therefore, the results support the idea that structural flexibility is an essential part of the construction principle of molecular motors through evolution. Such a feature is expected to be more prevalent in motor proteins than in simpler systems (e.g., signal transduction proteins) because in the former, large-scale conformational transitions often have to occur before the chemical events (e.g., ATP hydrolysis in myosin and ATP binding/phosphorylation in Ca²⁺-ATPase). This highlights the importance of Brownian motions associated with the protein domains that are involved in the functional transitions; in this sense, Brownian molecular machines is an appropriate description of molecular motors, although the normal mode results do not address the origin of the ratchet effect. The results also suggest that it might be more appropriate to describe functional transitions in some molecular motors as intrinsic elastic motions modulating local structural changes in the active site, which in turn gets stabilized by the subsequent chemical events, in contrast with the conventional idea of local changes somehow getting amplified into larger-scale motions. In the case of myosin, for example, we favor the idea that Brownian motions associated with the flexible converter propagates to the Switch I/II region, where the salt-bridge formation gets stabilized by ATP hydrolysis, in contrast with the textbook notion that ATP hydrolysis drives the converter motion. Another useful aspect of the BNM results is that selected low-frequency normal modes have been identified to form a set of collective coordinates that can be used to characterize the progress of a significant fraction of large-scale conformational transitions. Therefore, the present normal mode analysis has provided a stepping-stone toward more elaborate microscopic simulations for addressing critical issues in free energy conversions in molecular machines, such as the coupling and the causal relationship between collective motions and essential local changes at the catalytic active site where ATP hydrolysis occurs.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION COMPUTATIONAL METHODS STRUCTURAL FLEXIBILITY AND... DISCUSSION CONCLUSIONS SUPPORTING MATERIALS ACKNOWLEDGEMENTS REFERENCES

 
Motions are essential to the function of macromolecules (Brooks III et al., 1988; Gerstein and Krebs, 1998; McCammon and Harvey, 1987). Large-scale motions are commonly found in many biological molecules such as multidomain enzymes (Bahar et al., 1999; Thomas et al., 1999), signal transduction proteins (Ma and Karplus, 1997), and the ribosome (Frank, 2003). The large structural rearrangements are essential for the functions of these systems: e.g., binding and dissociation of substrates (Joseph et al., 1990); allosteric responses to signaling events; and release of synthesized proteins (peptides). The most dramatic class of systems where large-scale motions are crucial, however, has to include molecular machines (Banting and Higgins, 2000; Blumenfeld and Tikhonov, 1994; Hill, 1977; Schliwa, 2003) that convert free energies among various forms (see below). Although the mechanism of energy conversion has been a topic of great interest some 20 years ago (Hill, 1977; Hill and Eisenberg, 1981; Jencks, 1980; Simmons and Hill, 1976), the subject has been revived recently (Blumenfeld and Tikhonov, 1994). The free energy transduction processes can now be discussed in more structural and kinetic detail (Bustamante et al., 2001; Rees and Howard, 1999; Vale and Milligan, 2000), due in part to the observation of motor functions at unprecedented resolution offered by rapidly developing single molecule techniques (Ishii et al., 2003; Jung et al., 2002; Xie, 2002). Nevertheless, due to the challenges in making concurrent measurements at both high spatial and temporal resolutions, there remain many open questions concerning the detailed mechanisms of molecular motors.

The textbook description (Alberts et al., 1994) of molecular motors is that they harness the free energy of adenosine triphosphate (ATP) hydrolysis or phosphate transfer (i.e., phosphorylation) to perform net work. For example, kinesin motors can drag vesicle cargoes along microtubules in an unidirectional manner with ATP binding and/or hydrolysis (Hirokawa and Takemura, 2003). In each functional cycle, typically one ATP molecule in solution binds to the motor protein and is hydrolyzed (or phosphate-transferred), followed by the release of products back into solution; the net free energy change associated with the entire functional cycle is equivalent to the free energy of ATP hydrolysis in solution. This is the free energy that the molecular motor can, from the thermodynamics point of view, harness to generate useful work, although the actual hydrolysis obviously occurs in the protein rather than in solution (Hill, 1977). The precise molecular mechanisms for such free energy harness or conversion process have not been fully clarified. It is clear that it is not the free energies of the chemical reactions that are directly utilized because the exothermicity of these reactions in the protein, even if significant, will most likely dissipate to the rest of the molecule in pico-/nanoseconds without strongly preferred directions (Sagnella and Straub, 2001; Wang et al., 1998), which is too short compared to the timescale of motor functions (µs-ms) to have any functional significance. Moreover, the hydrolysis reaction in the hydrolyzing state of the ATPases often has an equilibrium constant close to 1 (Boyer, 1993; Holmes and Geeves, 1999). The key to the free energy conversion, therefore, has to do with the coupling of the protein structure and the ATP molecule in its various chemical states (Yang et al., 2003); for motors with tracks, such as myosin, the protein/track interaction is also an important element. Different protein or protein-track conformations could be preferred with different ATP states (solvated ATP, bound ATP, hydrolysis transition state, and ADP+P_i), and free energy is transferred between ATP and protein/track through well-coordinated conformational changes of the protein and binding/chemical modification of ATP.

It is our long-term goal to employ various types of simulations to provide the missing link between static high-resolution crystal structures and other lower-resolution studies such as single molecule investigations and phenomenological models for molecular motors. This, in light of the above discussion, necessarily involves understanding the coupling between the conformational properties of molecular motors and the chemical events; e.g., what are the structural flexibilities of various conformational states, what are the conformational transition pathways among these states, how does ATP affect these conformational properties and, in turn, how do different conformations determine the chemical fate of ATP? In the current work, we take the initial step to discuss the structural flexibilities in two distinct types of molecular machines: myosin-II and Ca²⁺-ATPase. The former is a prototypical molecular motor that utilizes the free energy of ATP binding and/or hydrolysis to perform mechanical work (Bagshaw, 2000; Keike and Titus, 2003), and the latter is the first P-type ion pump that has been characterized at atomic resolution (Lee, 2002; Toyoshima et al., 2000; Toyoshima and Nomura, 2002); Ca²⁺-ATPase transports calcium ions across the cell membrane against a concentration gradient with the help of ATP. Both undergo substantial conformational transitions during their functional cycles (see below). We set out to explore the structural flexibility inherent in these two systems and the corresponding functional significance; i.e., we inquire whether the requirement of making specific conformational transitions poses any constraints on the construct of those systems. Another important motivation is to explore useful schemes to characterize the functional transitions in conformations and identify the key degrees of freedom (Bustamante et al., 2001) such that, ultimately, the corresponding energetics and kinetics can be probed. The particular methodology we employed is an approximate normal mode analysis approach, the block normal mode (BNM) method (Li and Cui, 2002), which is based on the rotational-translational block (RTB) method of Tama and co-workers (Durand et al., 1994; Tama et al., 2000). As discussed below, although these normal mode results are highly approximate, they provided useful insights into the intrinsic structural features of myosin and Ca²⁺-ATPase; the results also revealed a set of collective coordinates that can be used in more elaborate microscopic simulations that can potentially offer more complete descriptions of the conformational transitions that occur in the functional cycle of molecular motors. In a recent article (Li and Cui, 2002), the low-frequency modes of Ca²⁺-ATPase in the calcium-loaded state (Toyoshima et al., 2000) have been discussed with the BNM method; here we focus on the correlation between low-frequency modes and the conformational transitions involving both the calcium-loaded and empty states, which was made possible by the recently solved x-ray structure for Ca²⁺-ATPase in the empty state (Toyoshima and Nomura, 2002).

The article is organized as follows. We first briefly reiterate the general value and limitation of normal mode analysis in studying large-scale conformational transition problems in biomolecules. We will then concisely review the approximate normal mode approach used in the calculations. Next, we present the normal mode results of myosin-II and Ca²⁺-ATPase, and discuss the relation between low-frequency modes of these proteins and their functions. Finally, we conclude with a number of summarizing statements.

	COMPUTATIONAL METHODS

TOP ABSTRACT INTRODUCTION COMPUTATIONAL METHODS STRUCTURAL FLEXIBILITY AND... DISCUSSION CONCLUSIONS SUPPORTING MATERIALS ACKNOWLEDGEMENTS REFERENCES

 
Normal mode analysis of large-scale conformational transitions: values and limitations
Although normal mode analysis has been applied to a large number of macromolecules (Bahar et al., 1997a; Case, 1994; Kitao and Go, 1999; Ma and Karplus, 1998; Marques and Sanejouand, 1995), it is useful to reiterate the value and limitation of this simple approach for analyzing conformational transitions in biological systems. The major challenge that hampers us from directly simulating large-scale conformational transitions in macromolecules is that these processes are diffusive and occur at a timescale (microseconds to milliseconds) much longer than what can be achieved using regular atom-based molecular dynamics with standard computational facilities (Elcock, 2002; Karplus and McCammon, 2002). One solution is to employ equilibrium simulations to explore thermodynamical quantities such as potential of mean force, in which time is not explicitly involved and kinetic properties can be estimated based on various rate theories (Hanggi et al., 1990); such types of calculations proved to be valuable for exploring the folding mechanism of small peptides and proteins (Shea and Brooks III, 2001). Another possibility is to use nonequilibrium simulations in which the biomolecule is biased along certain degrees of freedom, such as the root-mean square (RMS) difference from a given conformation (Schlitter et al., 1993); this type of calculation has been used to explore qualitative features associated with conformational transitions in molecular motors (Bockmann and Grubmuller, 2002; Ma et al., 2002) or ligand dissociations/translocations (Grubmuller et al., 1996; Isralewitz et al., 2001). In both schemes, the key is to identify appropriate reaction coordinate or order-parameters that concisely characterize the slow process. A fundamental assumption in those simulations is that all orthogonal degrees of freedom are in equilibrium during the simulation, which is usually difficult to justify for macromolecules but has to be made for practical purposes.

Compared to molecular dynamics-based methods such as principal component analysis (Amadei et al., 1993; Hayward et al., 1993; Kitao et al., 1991), normal mode analysis (NMA) is a simple but powerful approach for exploring the intrinsic structural flexibility (motion) of macromolecules (Brooks et al., 1995; Case, 1994; Go et al., 1983). The most significant advantage of NMA is that it is not limited by the timescale of the process under study, which makes it a unique tool for analyzing large-scale conformational transitions (Gerstein and Krebs, 1998). For example, several previous studies (Ma and Karplus, 1998; Tama and Sanejouand, 2001; Thomas et al., 1996, 1999) have illustrated that large-scale conformational changes can often be represented concisely by a very small number of low-frequency normal modes, which makes NMA an efficient approach for defining collective coordinates (Basu et al., 1995; Hayward and Go, 1995; Kitao and Go, 1999) that can be used in either equilibrium or nonequilibrium MD simulations; this is one of the important motivations for us to carry out NMA for the molecular motors of interest before more elaborate MD studies.

The most serious limitation of NMA is the harmonic approximation around a single conformational minima (Brooks et al., 1995), which seems to be a tremendously oversimplified model for macromolecules considering the rugged energy landscape associated with the high-dimensionality of these systems (Frauenfelder et al., 1991); another practical constraint is that NMA is usually carried out without the explicit representation of solvent molecules, which are well known to play important roles in determining the equilibrium and dynamical properties of macromolecules. Application of normal mode analysis to a large number of medium-size proteins, however, have shown that the characters of collective motions that correspond to the very low-frequency modes (Kitao and Go, 1999) are rather robust and NMA often gives similar results compared to molecular dynamics simulations (e.g., quasiharmonic analysis, Doruker et al., 2000; or singular value decomposition, Kundu et al., 2002). Physically, these results indicate that the motions associated with the low-frequency modes of macromolecules are most sensitive to the packing density of atoms; indeed, it has been shown that fairly robust results can be obtained for the low-frequency modes and patterns in the thermal atomic fluctuations by considering simplified models that contain only the packing distributions of atoms (Atilgan et al., 2002; Bahar et al., 1997b; Halle, 2002; Tirion, 1996). We emphasize that what is conserved in NMA is the direction of the collective motions, but not their actual timescales or thermal amplitudes; after all, these motions were found to be typically overdamped by the solvent molecules, which are usually neglected in NMA calculations. Moreover, for an analysis of conformational transitions, the directions of low-frequency modes are more relevant than their thermal amplitude (see below). Another important effect of the solvent molecules is that their presence increases the number of local minima on the potential energy surface, thus there are frequent transitions among various minima in solution. However, at least for small systems such as melittin, calculations indicated that the envelope of the potential energy surface in solution can be predicted correctly by normal mode analysis in vacuum (Kitao et al., 1991). In other words, although there may not be 1:1 correspondence between normal modes and the more realistic principal components that include anharmonic motions of the system, the subspace that dominates the atomic fluctuations can be captured by a linear combination of the low-frequency, collective normal modes (Kitao and Go, 1999; Kitao et al., 1998).

Approximate normal mode analysis: the block normal mode (BNM) approach
Despite the observed success of simplified elastic models in describing low-frequency modes of macromolecules (Doruker et al., 2000; Kundu et al., 2002; Tama and Brooks III, 2002), it is still desirable to be able to explore these motions with atomic level details and a physical representation of intermolecular interactions to ensure robustness of the calculations. This is especially true when the molecular system involves heterogeneous distributions of different chemical groups, such as in ribosome and RNA polymerase, which involves both protein and nucleic acids. To make connections between the normal modes of biomolecules and their mechanical properties (e.g., elastic modulus), a model with physical interactions is also more useful.

The computational expense associated with normal mode analysis in an all-atom representation, unfortunately, scales rather steeply with the size of the system (N atom); nominally, the Hessian matrix scales as N², and the diagonalization of this Hessian scales as N³ with the most standard diagonalization procedures. To reduce the computational cost, we have recently made an efficient implementation of the BNM approach (Li and Cui, 2002), which was based on the RTB method originally proposed by Tama and co-workers (Durand et al., 1994; Tama et al., 2000). In RTB, the macromolecule is first partitioned into blocks, where each block may range from a single residue to a secondary structure element depending on the problem in hand. The full atomic Hessian matrix is then projected onto a subspace spanned by the translation and rotation vectors associated with these blocks, which results in a much smaller matrix to be diagonalized. BNM is more general than RTB in that other eigenvectors of the blocks can also be included, which is appropriate when blocks are large in size. The merit of the RTB/BNM approach is that the cost of calculation is dramatically reduced with the full account of intermolecular interactions; the coarse-graining procedure does not compromise the accuracy of results when only the low-frequency, collective modes of the macromolecule are of interest (Li and Cui, 2002; Tama et al., 2000). Test calculations using small protein systems (<200 residues) and one residue per block indicated that a very good correlation (R² 0.997) between the BNM frequencies and standard NMA frequencies was observed up to at least 40 cm^-1; although the overlap between the BNM eigenvectors and standard NMA eigenvectors become <0.5 after the frequency reaches 10 cm^-1, the subspace spanned by the NMA eigenvectors can be represented >90% by the BNM eigenvectors for frequencies up to 25 cm^-1.

Compared to the original work of Tama and co-workers (Durand et al., 1994; Tama et al., 2000), our BNM implementation has the essential advantage that the projected Hessian matrix is constructed in a direct fashion (i.e., without the construction of the entire atomic Hessian), which makes the application to fairly large systems (N < 1000 residues) very straightforward. Recently, we have substantially increased the applicability of BNM by taking advantage of a sparse matrix diagonalization procedure (Cullum and Willoughby, 1985; Lanczos, 1950) implemented in the PARPACK package (Maschhoff and Sorensen, 1996) and parallel computations. These refinements made it possible to analyze collective motions of supermolecular assemblies such as the ribosome at atomic resolution with modest computational cost; BNM calculations for a 50S ribosome (128,000 atoms with a polar hydrogen set of force field, Brooks et al., 1983) with one residue per block took <17 h on a single node of a 2.0-GHz Athlon. The details of the new implementation and application to large molecular assemblies will be given separately.

Computational models and protocols
The newly implemented BNM method in CHARMM (Brooks et al., 1983) was used to analyze the structural flexibilities in two molecular motors where large-scale motions are critical to their functions: Ca²⁺-ATPase (Lee, 2002) and myosin (Holmes and Geeves, 1999; Steffen et al., 2003). For myosin, calculations were done for the prehydrolyzing state (see Scheme 1) with bound Mg²⁺·ATP (PDB code 1FMW, Bauer et al., 2000); the structures for the motor domain without ATP or ADP were found to be very similar to 1FMW (Bauer et al., 2000; Gulick et al., 1997) and therefore have not been calculated. Substantial conformational transition occurs in the ATP hydrolyzing state (PDB code 1VOM, Smith and Rayment, 1996), and therefore the normal modes for the prehydrolysis state was analyzed in terms of the conformational difference with the hydrolyzing state. Due to the many missing residues in the 1VOM structure, no normal mode analysis has been done for that structure.

View larger version (14K): [in this window] [in a new window]   SCHEME 1  The functional cycle of myosin-II-actin (only a myosin-monomer is shown).

View larger version (23K): [in this window] [in a new window]   SCHEME 2  The functional cycle of Ca2+-ATPase in the E1-E2 model.

Starting from the x-ray structures (see Table 1 for a summary), the systems were first energy-minimized before normal mode analysis. The minimization protocol started with 500 steps of steepest descent followed by 5000 steps of adapted basis Newton-Rhapson (Brooks et al., 1983). To prevent undesired structural distortions, gradually decreasing harmonic constraints with respect to the x-ray structure as the reference was used in the minimization. Another 3000 adapted-basis Newton-Rhapson minimization steps were finally carried out without constraints, which produces the final structure with a typical RMS gradient on the order of 0.01 kcal/(mol x Å). Such a threshold is sufficient for the current purpose, as reflected by the fact that the normal mode analysis does not produce imaginary modes. Moreover, the translation and rotation modes have frequencies on the order of 0.05 to 0.10 cm^-1. The RMS deviations between x-ray and minimized structure are typically 1.5 Å (see Table 1).

The results of the BNM calculations were analyzed in terms of the root-mean square atomic fluctuations (RMSF, r²), participation ratios (Leitner, 2001; Sagnella and Straub, 1999), and visual inspections. The RMSF and participation ratios associated with the k^th normal mode are defined by

(1)

(2)

where M^res is the number of residues in the system, m_j is the mass of the j^th atom, and N_i is the number of atoms in the i^th residue; L_jk and _k are the j^th components and vibrational frequency associated with the k^th eigenvector, respectively. As discussed extensively in previous work (Leitner, 2001; Sagnella and Straub, 1999), the reciprocal of R_k corresponds approximately to the number of residues involved in the k^th normal mode. To further characterize the motions involved in the normal modes, the Dyndom package (Hayward et al., 1997) was used to divide the structure into dynamical domains (to be distinguished from the structural domain identified based on static geometry) and identify hinge-bending residues associated with the relative displacements of these dynamical domains; for each normal mode analyzed, two structures were generated by following the corresponding eigenvector with phase angles of 90° and 270°, respectively, which were then used as inputs for Dyndom.

For myosin and Ca²⁺-ATPase, where multiconformations are known from x-ray crystallography, the individual involvement coefficients were also computed to reveal the connection between the normal modes and the conformational transitions (Tama and Sanejouand, 2001),

(3)

where X₁–X₂ is the displacement vector between two conformations. A large involvement coefficient indicates that the motion along this normal mode is highly relevant to the conformational transition being examined. A related quantity that indicates the weight of a set of normal modes in a conformational transition is the cumulative involvement coefficients,

(4)

Since the normal mode eigenvectors of a macromolecule form a complete orthonormal set for the 3N-dimensional space, CI_3N is expected to be 1; in other words, measures the maximal percentile contribution of exciting motions along the k^th mode to the conformational difference between the two states.

For an analysis of conformational transitions, the involvement coefficients from normal mode analysis are of great interest because they measure the relative importance of motions along a given normal mode. If the involvement coefficients are dominated by low-frequency normal modes, it is reasonable to expect that Brownian motions along these directions can accomplish large-scale conformational transitions with significant probability/rate. We note that thermal amplitude of a normal mode, which is an equilibrium quantity defined in terms of the mean-square displacement of the oscillation, is less relevant (compared to the direction of the mode) in conformational transition problems; therefore, the fact that normal mode analysis with simplified solvent treatment does not give quantitative thermal amplitudes is of little concern for our purpose. Take a simple example of gauche-trans isomerization of butane. The isomerization is dominated by a single torsional mode, which implies an involvement coefficient close to 1 for that mode. Although a normal mode analysis cannot directly predict any rate constant, the fact that the isomerization motion is dominated by a very low-frequency torsional mode implies that the barrier is likely to be low, relative to other type of processes such as breaking of the C–C bond. The thermal amplitude of the torsional mode is small at low temperatures, but this is obviously not contradictory to the importance of the torsional mode in the isomerization.

	STRUCTURAL FLEXIBILITY AND FUNCTION OF BROWNIAN MOLECULAR MACHINES: MYOSIN AND CA2+-ATPASE

TOP ABSTRACT INTRODUCTION COMPUTATIONAL METHODS STRUCTURAL FLEXIBILITY AND... DISCUSSION CONCLUSIONS SUPPORTING MATERIALS ACKNOWLEDGEMENTS REFERENCES

 
In this section, the structural flexibilities of myosin and Ca²⁺-ATPase are first discussed in terms of their low-frequency modes from BNM calculations. The functional relevance of these motions is illustrated with involvement coefficients, which provide a semiquantitative scheme for characterizing large-scale conformation transitions that are valuable for more elaborate MD simulations in the future. In the last subsection, we put the results into a more general framework and discuss the value of low-frequency modes in understanding motor functions.

Low-frequency modes of the motor domain in myosin
As shown in Fig. 1 a, the majority of conformational transitions in the myosin motor domain occurs in the C-terminal region, which is known as the converter; during the functional cycle of myosin, it is the large-scale displacement of this converter region that guides the motion of the longer level-arm (absent from the structure used here). When the two states of myosin (ATP prehydrolysis state and ATP hydrolyzing state, which will be referred to by their PDB code names, 1FMW, Bauer et al., 2000; and 1VOM, Smith and Rayment, 1996, respectively) are best fit with the first 650 amino acids, the total root-mean square difference (RMSD) for these residues is only 2.4 Å; the rest of the residues in the C-terminal region after such alignment, by contrast, have a much larger RMSD of 20.8 Å (Fig. 1 a). The displacement of the C-terminal region can be approximately described as a rigid body motion because residues in this region have a small RMSD of 2.5 Å after they are superimposed in the two states. In addition to the C-terminal part, two nearby regions, the N-terminal region (Asp-23) and a long -helical region (Ser-465 to Glu-497, the "relay" helix) adjacent to the ATP binding site also undergo substantial changes (Fig. 2 a). Among the three structural motifs essential to ATP binding and hydrolysis (P-loop, Switch I, and Switch II), the Switch II region undergoes a substantially larger change compared to the other two (Fig. 1 b); it corresponds to the formation and break of a salt-bridge (Arg-238 and Glu-459), which was proposed to be important for ATP hydrolysis based on mutation (Onishi et al., 2002) and computational studies (G. Li and Q. Cui, unpublished). An essential issue for the mechanochemical coupling in myosin concerns the relation between the larger-scale C-terminal motion and the more localized changes in the ATP binding site, and to what extent these motions are reflected by the intrinsic structural flexibility of the motor domain (i.e., low-frequency modes).

View larger version (44K): [in this window] [in a new window]   FIGURE 1  The difference between the prehydrolysis (1FMW) and hydrolyzing (1VOM) state of the Dictyostelium discoideum myosin-II motor domain. (Left) The superposition of the motor domain in the two states based on the backbone atoms in residues 1–650; red and blue corresponds to 1FMW and green and yellow indicates the 1VOM state. As illustrated by the arrow, the majority of the conformational difference occurs in the C-terminal converter region. (Right) The structural difference between the active sites of the two states, with the same color coding; the ATP molecule is shown in van der Waals representation, and the three important motifs for ATP hydrolysis are shown in ribbon-type representation.

View larger version (36K): [in this window] [in a new window]   FIGURE 2  The correspondence between structural difference and low-frequency normal modes in myosin-II motor domain. (a) The displacements of C atoms between the prehydrolysis (1FMW) and hydrolyzing (1VOM) structures, with the essential active site motifs highlighted; (b) the root-mean square fluctuations (RMSF) of C atoms in the prehydrolysis state calculated using the 100 lowest-frequency modes at 300 K (Eq. 1); (c) the individual involvement coefficients (Eq. 3) for the first 100 modes; and (d) the cumulative involvement coefficients (Eq. 4). In c and d, the overall translational and rotational modes are shown with negative mode indices; they contribute because the structures are superimposed with residues 1–650 (see text).

The relevance of the low-frequency modes in myosin to the conformational transition is more quantitatively illustrated by the involvement coefficients (Fig. 2, c and d). The individual involvement coefficients (I_k) are dominated by a number of low-frequency modes; e.g., only two intramolecular modes have I_k > 0.20. Since the two states (1FMW and 1VOM) were superimposed based on the first 650 residues, the overall translation and rotation also have notable contributions, although superimposing the two structures based on all heavy atoms gave very similar results in terms of involvement coefficients of internal normal modes. For example, one rotational mode has a large involvement coefficient 0.30 (Fig. 2 c), which highlights the fact that the conformational transition substantially changes the center of mass and moments of inertia of the system (note that intramolecular normal modes conserve the center of mass properties). The cumulative involvement coefficients (Fig. 2 d) indicate that nearly 50% of the overall conformational changes can be described in terms of motions in the direction of the 10 lowest-frequency modes. After the 10th mode, the increase of the cumulative involvement coefficients becomes substantially slower (Fig. 2 d); including all 100 modes, nearly 60% of the conformational changes can be accounted for.

The lowest-frequency internal mode (₁ = 1.8 cm^-1) has the largest involvement coefficient (I₁ = 0.32), and it mainly involves the displacement of the C-terminal converter region and the upper 50-K domain (Fig. 3 a). As shown in Fig. 4 a, the converter region undergoes a rigid-body rotation, which is precisely the motion that corresponds to the essential conformational transition between 1VOM and 1FMW. Dynamical domain analysis (Hayward et al., 1997) using the two end-structures found the hinge axis, which is also shown in Fig. 4 a; the hinge regions (see Table 2) include residues Gln-112–Thr-117, Phe-482, Ile-514, and Gly-680, which are shown as the dark green region in the Dyndom plot (also see Supporting Materials for the positions of the hinge residues in the structure); Phe-482 and Gly-680 are in close contact in the crystal structure, with Phe-482 on the relay helix whereas Gly-680 is on the so-called SH1–SH2 helix (Fig. 1 a). To conserve the center of mass, a loop and the long -helix in the upper 50-K domain rotate in the opposite direction. Although the difference in the converter is substantially larger than that in the upper 50-K domain between the two conformational states (Fig. 1 a), the magnitude of their motion is comparable in the normal modes, which is imposed by the conservation of the center-of-mass position. The magnitude of rotation in the converter domain predicted from the normal mode analysis at 300 K is much smaller than the difference between the actual x-ray structures of the two states (Fig. 1), which reflects the limitation of the harmonic assumption inherent in the normal mode analysis. As mentioned above, the active site has relatively small degree of motion in the low-frequency modes. In the first mode, for example, the observed displacement in Switch II (Fig. 4 b) has little overlap with the functional transition between the prehydrolysis to the hydrolyzing states (Fig. 1 b).

View larger version (50K): [in this window] [in a new window]   FIGURE 3  The RMSF of C atoms in selected normal modes of the prehydrolysis state myosin-II motor domain (1FMW) that have large individual involvement coefficients. For each mode, the three numbers above each plot are the vibrational frequency (in cm-1), reciprocal participation ratio (Eq. 2), and the individual involvement coefficient (Eq. 3).

View larger version (111K): [in this window] [in a new window]   FIGURE 4  Superposition of end-structures (in cartoon format) in selected normal modes of the prehydrolysis state myosin-II motor domain (1FMW) that have large individual involvement coefficients. The end-structures are generated based on the corresponding equation of motion, and they correspond to a phase angle of 90° and 270°, respectively; to clearly illustrate the structural variations, a high temperature of 1500 K was used. In each mode, red and blue indicate one structure; green and yellow correspond to the other. The arrows indicate qualitative in-plane motion of the modes, whereas "" indicates out-of-plane rocking motion. Similar to Fig. 3, the three numbers in each plot are the vibrational frequency (in cm-1), reciprocal participation ratio (Eq. 2), and the individual involvement coefficient (Eq. 3). Also shown are the results (in Trace format) of dynamical domain analysis using the Dyndom software with these end-structures. Different domains are colored differently, with dark green indicating regions undergoing significant bending going from one end-structure to the other; for the hinge residues, see Table 2 and Supporting Materials. The arrows indicate the hinge-bending axis, where the colors of the axis and arrowhead indicate the two dynamical domains involved in the relative motion. In b, displacements in the ATP (in vander Waals representation) binding site are shown (see text).

Low-frequency modes of the Ca²⁺-ATPase in different functional states
As shown in Fig. 5, large changes are observed in the TMD and cytoplasmic domains between E1·Ca²⁺ and E2, and the magnitude of conformational change depends on the scheme used to superimpose the two structures. If best fitted based on TMD, the N, P, and A domains in the two states have very large RMSDs of 39.9, 14.3, and 20.2 Å, respectively; the TMDs have a RMSD of 5.8 Å. If best fitted based on all heavy atoms, the RMSDs for N, P, and A are 17.0, 9.4, and 19.3 Å, respectively, and that for TMDs is 10.3 Å. Therefore, aligning the structures based on TMD significantly amplifies the displacement of the N domain (Fig. 5). It is not immediately clear which is the better way to characterize the conformational transition: on one hand, the TMD is held by the lipid and thus it appears to be more physical to minimize the displacements of TMD during the best-fit procedure; on the other hand, normal mode analysis indicates that the low-frequency modes contain significant components in TMD that correspond to the transition between the two states when best fitted based on all heavy atoms. Qualitative trends for the correlation between low-frequency modes and conformational transition are the same using either scheme for alignment. Therefore, in the following discussion, we focus on results obtained by fitting the E1·Ca²⁺ and E2 structures based on all heavy atoms. It should be noted that the internal structures of the N, P, and A domains change little in the two states, especially for N and P; the internal RMSD in the three domains are 1.5, 1.2, and 3.6 Å, respectively. Therefore, the dominant components in the conformational transition between E1·Ca²⁺ and E2 involve rigid domain-motions of N, P, and A relative to TMD, and to a somewhat smaller extent, the internal structural changes in TMD.

View larger version (53K): [in this window] [in a new window]   FIGURE 5  The superposition of calcium-loaded and calcium-free Ca2+-ATPase structures when they are aligned based on (a) backbone atoms in the transmembrane domain (TMD) and (b) all backbone atoms. In c and d, the corresponding displacements in C are shown respectively. In a and b, brown-red-purple-dark blue corresponds to the E1·Ca2+ state, and gray-green-blue-yellow corresponds to E2. In c and d, pale blue indicates the A domain, yellow indicates the N domain, pink indicates the P domain, and the rest corresponds to the TMD.

View larger version (14K): [in this window] [in a new window]   FIGURE 6  Individual (a, b) and cumulative (c) involvement coefficients calculated when the calcium-loaded/-free forms are aligned based on all backbone atoms. In a and b, the eigenvectors associated with the calcium-loaded and calcium-free structures are used, respectively. The overall translational and rotational modes are shown with negative mode indices, although they do not contribute significantly (see text).

View larger version (67K): [in this window] [in a new window]   FIGURE 7  The RMSF of C atoms in selected normal modes of the calcium-loaded form (E1·Ca2+) of Ca2+-ATPase that have large individual involvement coefficients. For each mode, the three numbers above each plot are the vibrational frequency (in cm-1), reciprocal participation ratio (Eq. 2), and the individual involvement coefficient (Eq. 3).

View larger version (55K): [in this window] [in a new window]   FIGURE 11  Similar to Fig. 10, but for the normal modes of the calcium-free form (E2) of Ca2+-ATPase.

View larger version (49K): [in this window] [in a new window]   FIGURE 10  Superposition of end-structures in selected normal modes of the calcium-loaded form (E1·Ca2+) of Ca2+-ATPase that have large individual involvement coefficients. The end-structures are generated based on the corresponding equation of motion, and they correspond to a phase angle of 90° and 270°, respectively; to clearly illustrate the structural variations, a high temperature of 1500 K was used. The dotted circles indicate that the enclosed domains move almost like a rigid body. The brown-red-purple-dark blue corresponds to one end-structure, and gray-green-blue-yellow corresponds to the other. Also shown are the results (in Trace format) of dynamical domain analysis using the Dyndom software with these end-structures. Different domains are colored differently, with dark green indicating regions undergoing significant bending going from one end-structure to the other; for the hinge residues, see Table 3 and Supporting Materials. The arrows indicate the hinge-bending axis, where the colors of the axis and arrowhead indicate the two dynamical domains involved in the relative motion.

In E2, the modes with large involvement coefficients are the first, second, fourth, and eighth modes, and the 14th, 16th, 17th, and 18th modes (Fig. 6 b); the involvement coefficients are the largest for the fourth and 18th modes, in contrast to the situation in E1·Ca²⁺ (Fig. 6 a). Due to the more compact structures, especially in terms of the orientation of the N domain and its interaction with the A domain, the normal modes in E2 are more complex and usually involve multiple hinges. Moreover, the modes in E2 exhibit evident internal motions of the TMD, with many hinge residues in the middle of transmembrane helices (see Table 3 and Supporting Materials). In the first mode, the P and A domains move approximately as a group, and the motion involves the rocking of this group, the N domain, and the twisting of the TMD (Figs. 8 a and 11 a). The hinges for the P-A motions are at the interface between P and M5 (Ile-743–Ala-746), whereas the hinges for the TMD twist involves a number of residues on M4 (Ile 315) and M5 (Lys-758–Ile-761). The motion of the second mode is rather complex and involves a number of relative displacements of various groups (Figs. 8 b and 11 b); these include the rocking of the N and A domains in directions perpendicular to those in mode 1, and twists of the TMD (for the displacement axis, see Fig. 11 b). The hinge residues exist on the interface between N/P, A/M1, A/M3, and A/P; two residues on M5, Tyr-763, and Leu-764 were found as hinges for the internal motion of the TMD. In the fourth mode (Figs. 8 c and 11 c), the N, P, and A domains and TMD move almost like rigid bodies, and the hinge residues occur at the interfaces between them (Table 3). In the eighth mode (Fig. 8 d and 11 d), the motion is rather distinct: it involves the relative opening and closing of the N and A domains, and complex motion of the TMD; there are many hinge residues that occur in different parts of TMD, including M3, M4, M5, and M8 (Table 3). The few higher-frequency modes (3 cm^-1) that have large involvement coefficients (Fig. 9) are more complex in terms of atomic motions, and include various relative displacements of the N, P, and A domains and TMD; the TMD also undergo notable intradomain rearrangement (see the figure in Supporting Materials).

View larger version (80K): [in this window] [in a new window]   FIGURE 8  The RMSF of C atoms in selected normal modes of the calcium-free form (E2) of Ca2+-ATPase that have large individual involvement coefficients. For each mode, the three numbers above each plot are the vibrational frequency (in cm-1), reciprocal participation ratio (Eq. 2), and the individual involvement coefficient (Eq. 3).

View larger version (68K): [in this window] [in a new window]   FIGURE 9  Similar to Fig. 8, but for a few modes of substantially higher frequencies that have large involvement coefficients for E2.

	DISCUSSION

TOP ABSTRACT INTRODUCTION COMPUTATIONAL METHODS STRUCTURAL FLEXIBILITY AND... DISCUSSION CONCLUSIONS SUPPORTING MATERIALS ACKNOWLEDGEMENTS REFERENCES

For the two systems studied, which represent two distinct classes of molecular machines that convert free energy among various forms (i.e., chemical mechanical in myosin, chemical ion gradient in Ca²⁺-ATPase), it was found that the various conformational states indeed have the flexibilities that are highly correlated with the conformational transitions required by the function. The direct numerical evidence is that the individual involvement coefficients are significant for only a small number of low-frequency modes. In myosin, 45% of the conformational change that converts the system from the prehydrolysis to the hydrolyzing state can be described with motions in the direction of the first 10 normal modes. In the E1·Ca²⁺ state of Ca²⁺-ATPase, strikingly, 70% of the conformational difference from the E2 state can be described with only the 10 lowest-frequency modes. With the 10 lowest-frequency normal modes of the E2 state, however, <40% of the expected conformational transition can be accounted for, which implies that energy rises more rapidly around E2 compared to E1·Ca²⁺ during the isomerization (i.e., the initial transitions involves motions along modes of higher frequencies).

Intrinsic structural flexibilities in proteins that perform functions through conformational transitions have been probed with NMR for signal transduction proteins (Volkman et al., 2001). For several receiver domains involved in two component systems, for example, it was argued based on NMR relaxation and chemical shifts measurements (Volkman et al., 2001) that these proteins have the intrinsic structural flexibility in regions that undergo large-scale conformational transitions. Therefore, in contrast to the traditional mechanism in which conformational transition is triggered by chemical events (e.g., phosphorylation), it was argued that the protein has certain population in the active conformation even without the chemical activation; the role of the chemical event is simply to shift more populations toward the active conformation. In light of the normal mode results discussed here, we speculate that many molecular motors work through such a population-shift mechanism; i.e., the motor has intrinsic structural flexibility such that Brownian motions can lead to significant population of specific conformations where chemistry (e.g., ATP hydrolysis) occurs. For example, myosin has to switch from the open to the closed state to allow ATP hydrolysis, whereas Ca²⁺-ATPase has to swing the N (by >25 Å) and A domains to allow ATP binding and the subsequent phosphorylation. For myosin, this has been demonstrated with bulk tryptophan fluorescence studies for a single tryptophan-containing mutant of Dictyostelium myosin (Malnasi-Csizmadia et al., 2001a,b, 2000), which revealed that the open (1FMW)-closed (1VOM) transition does not require ATP hydrolysis per se and can be induced by a nonhydrolyzable ATP analog. The flexible constructs of myosin and Ca²⁺-ATPase is an important component of the Brownian ratchet model, which has been proposed to rationalize the unidirectional character of molecular motors (Astumian and Hanggi, 2002; Julicher et al., 1997). The ratchet effect, however, cannot be addressed directly with the normal mode results. In the literature (Astumian and Hanggi, 2002; Julicher et al., 1997), it is often loosely stated that the ratchet effect arises because ATP hydrolysis breaks detailed balance by driving the system locally out of equilibrium. We emphasize that the term out-of-equilibrium should not refer to the hydrolysis reaction per se, because the equilibrium constant associated with the ATP hydrolysis step is usually measured at close to 1 in motor systems (Boyer, 1993; Holmes and Geeves, 1999). In other words, ATP and ADP·Pi are in constant equilibrium, which has been illustrated by ¹⁸O isotope exchange results in systems such as F₁-ATPase (Bagshaw et al., 1975) and myosin (Dale and Hackney, 1987). As already mentioned in the Introduction, this makes sense, because otherwise the energy released will be dissipated to the environment rather rapidly, which compromises the efficiency of the motor. Therefore, it is likely that the ratchet effect is related to other features of the functional cycle; e.g., it might be that products of the chemical events stabilize the new conformation over the initial one, or further events can occur with the new conformation; in the case of myosin, these correspond to the possibilities that ADP·Pi stabilize the closed conformation over the open (G. Li and Q. Cui, unpublished), and that the ADP·Pi state favors rebinding of myosin to actin, which is practically irreversible.

For most important conformational transitions, one expects a combination of collective displacements and local rearrangements in the active site. In myosin, the collective displacement mainly concerns the converter region, whereas the local rearrangement involves the salt-bridge region between Switch I and Switch II (Fig. 1 b), which is essential for the chemistry in the mechanochemical cycle (Holmes and Geeves, 1999; G. Li and Q. Cui, unpublished; Rayment, 1996). In calcium pump, the collective motion mainly involves displacements of the N and A domains, as well as some adjustments in the TMD, whereas local changes involve the reorientation of the side chains in Glu-58, Glu-771, and Glu-908 (Lee, 2002; Toyoshima and Nomura, 2002