Article Information

• PDF (840 kb)
• Full Text with Large Figures
• Supporting Material
• Cited by in Scopus (1)

PubMed

• Articles by Fumiko Takagi
• Articles by Macoto Kikuchi

Related Articles

• Correction

  Correction Biophysical Journal, Volume 87, Issue 2, 1 August 2004, Pages 1398 Full Text | PDF (43 kb)

• Toward a Molecular Understanding of the Anisotropic Response...

  Toward a Molecular Understanding of the Anisotropic Response of Proteins to External Forces: Insights from Elastic Network Models Biophysical Journal, Volume 94, Issue 9, 1 May 2008, Pages 3424-3435 Eran Eyal and Ivet Bahar Abstract With recent advances in single-molecule manipulation techniques, it is now possible to measure the mechanical resistance of proteins to external pulling forces applied at specific positions. Remarkably, such recent studies demonstrated that the pulling/stretching forces required to initiate unfolding vary considerably depending on the location of the application of the forces, unraveling residue/position-specific response of proteins to uniaxial tension. Here we show that coarse-grained elastic network models based on the topology of interresidue contacts in the native state can satisfactory explain the relative sizes of such stretching forces exerted on different residue pairs. Despite their simplicity, such models presumably capture a fundamental property that dominates the observed behavior: deformations that can be accommodated by the relatively lower frequency modes of motions intrinsically favored by the structure require weaker forces and vice versa. The mechanical response of proteins to external stress is therefore shown to correlate with the anisotropic fluctuation dynamics intrinsically accessible in the folded state. The dependence on the overall fold implies that evolutionarily related proteins sharing common structural features tend to possess similar mechanical properties. However, the theory cannot explain the differences observed in a number of structurally similar but sequentially distant domains, such as the fibronectin domains. Abstract | Full Text | PDF (1915 kb)

• Multiphysics Simulation of Left Ventricular Filling Dynamics...

  Multiphysics Simulation of Left Ventricular Filling Dynamics Using Fluid-Structure Interaction Finite Element Method Biophysical Journal, Volume 87, Issue 3, 1 September 2004, Pages 2074-2085 Hiroshi Watanabe, Seiryo Sugiura, Hidenobu Kafuku and Toshiaki Hisada Abstract To relate the subcellular molecular events to organ level physiology in heart, we have developed a three-dimensional finite-element-based simulation program incorporating the cellular mechanisms of excitation-contraction coupling and its propagation, and simulated the fluid-structure interaction involved in the contraction and relaxation of the human left ventricle. The FitzHugh-Nagumo model and four-state model representing the cross-bridge kinetics were adopted for cellular model. Both ventricular wall and blood in the cavity were modeled by finite element mesh. An arbitrary Lagrangian Eulerian finite element method with automatic mesh updating has been formulated for large domain changes, and a strong coupling strategy has been taken. Using electrical analog of pulmonary circulation and left atrium as a preload and the windkessel model as an afterload, dynamics of ventricular filling as well as ejection was simulated. We successfully reproduced the biphasic filling flow consisting of early rapid filling and atrial contraction similar to that reported in clinical observation. Furthermore, fluid-structure analysis enabled us to analyze the wave propagation velocity of filling flow. This simulator can be a powerful tool for establishing a link between molecular abnormality and the clinical disorder at the macroscopic level. Abstract | Full Text | PDF (688 kb)

• …more

Copyright © 2007 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 93, Issue 11, 3820-3827, 1 December 2007

doi:10.1529/biophysj.106.103796

Biophysical Theory and Modeling

Structural Change and Nucleotide Dissociation of Myosin Motor Domain: Dual Gō Model Simulation

Fumiko Takagi^*, †, ,   and Macoto Kikuchi^†, *

^* Formation of Soft Nanomachines, Core Research for Evolutional Science and Technology, Japan Science and Technology Agency, Osaka, Japan
^† Cybermedia Center Osaka University, Osaka, Japan

Address reprint requests to F. Takagi.

The mechanism of biomolecular motors is one of the major topics in biophysics. Among a number of such systems that have been found so far, the actomyosin motor is of particular interest, because it is responsible for muscle contraction and cellular movements in eukaryotic cells. Myosin moves unidirectionally along the actin filament using chemical energy released by ATP hydrolysis 1,2,3,4,5. In contrast to macroscopic artificial machine, biomolecular motors work under a noisy environment in the cell. In fact, the free energy released during each ATP hydrolysis is only ∼20k_BT; therefore the thermal fluctuation should be appreciable. Although recent progress in imaging and of nanomanipulation has enabled the observation of single molecules, the movement mechanism of the actomyosin motor is still not understood.

There has been a long-standing controversy between the tight-coupling (lever-arm) model and the loose-coupling model. X-ray crystallographic studies have revealed that the angle of the neck domain changes relative to the motor domain, depending on the nucleotide state. The lever-arm model was proposed based on these observations, in which the structural change of myosin is tightly coupled with the ATP hydrolysis cycle and directly causes a stepwise sliding motion. It was shown, however, that the sliding distance of the myosin along the actin filament per ATP at the muscle contraction can be much longer than the displacement predicted by the lever-arm-like structural change of a single myosin molecule 6. Moreover, it is questionable whether a material as soft as proteins can accurately switch its conformation in the same way as a macroscopic machine under thermal fluctuation.

In the loose-coupling model, in contrast, the structural change does not always correspond to a step in a one-to-one correspondence; the motion is intrinsically stochastic and thermal fluctuation is an essential ingredient for its mechanism 7. The simplest class of models that produce the loose-coupling mechanism is based on a thermal ratchet, in which a myosin molecule is treated as a Brownian particle that moves along a periodic and asymmetric potential under both thermal noise and nonthermal perturbations 8,9. Although ratchet systems can, in fact, exhibit unidirectional flow in a noisy environment, a high efficiency comparable to that of the actomyosin system is found difficult to achieve using only a simple ratchet system. Even if the ratchet models do express some essence of the mechanism of the biomolecular motor, they are too much simplified and we should say that the connection with the real actomyosin system is rather vague. In particular, since the myosin is expressed as a particle, the effect of its conformational change is, at best, only implicitly taken into account. A somewhat more realistic modeling is desirable, which can reflect the chain conformation.

Recently, it was revealed by single molecule experiments that the chemomechanical cycle of the myosin head is controlled by a load on the actomyosin cross bridge 10,11,12. The observation suggested that the rate of ADP release from the myosin head depends on the force acting on myosin; namely, the chemical reaction rate varies with the deformation of the myosin. If the myosin head indeed acts as such a strain sensor, this would be reminiscent of a classical model by A. Huxley 13; in this model, the myosin head is supposed to undergo Brownian motion and change into a tightly bound state to the actin filament triggered by a structure dependent chemical reaction.

The relationship between structure and function of proteins has long been investigated. Thus far, mainly the static aspects of proteins have been considered; for example, the classical lock-and-key model of an enzyme. Recently, the role of structural fluctuations, or that of a more drastic structural change including partial unfolding on protein functions, has become a subject of growing interest. Although there have been many experimental studies to clarify the dynamical processes of protein at a functional level, it is still difficult to observe the structural changes with high resolution in both space and time. Computer simulations serve as possible alternatives. Typical computational studies treated equilibrium fluctuations near a crystal structure using the all-atom model 14,15,16 or the elastic network model 17,18,19,20; although these types of simulations are not suitable to study large-scale structural change, the low-frequency fluctuation modes were found to be consistent with the direction of motion of the structural change associated with the function.

Some attempts have been made to simulate a larger structural change beyond the elastic regime using a class of models called the Gō-like model 21. According to the recently developed theory of spontaneous protein folding, the protein energy landscape has a funnel-like global shape toward the native structure. The Gō-like model is certainly the simplest class of model that can realize a funnel-like landscape 22,23 and has successfully described the folding process of small proteins. It is, however, not suitable for the study of a change between two conformations, because only the interaction between the pairs of residues that contact each other in the native state are taken into account; the conformation other than the native state becomes too unstable as a result. Here, we introduce a dual Gō model in which two conformations can be embedded as a variant of the coarse-grained Gō-like model.

In this article, we investigate the dynamics of the spontaneous myosin conformational change from the pre-power stroke state to the near-rigor state by molecular dynamics simulations of dual Gō model. This process is called power-stroke because the angle of the lever-arm changes remarkably, and it is considered in the lever-arm model that this structural change directly causes the force generation. The spontaneous dissociation process of the nucleotide that accompanies the conformational change is also involved in the simulation by introducing a coarse-grained nucleotide as a short chain. To our knowledge, the ligand at the binding site has not been considered explicitly in coarse-grained protein simulations, possibly because the primary role of ATP is considered to be the release of chemical energy through hydrolysis, and the excluded volume effect of the molecule has not been investigated. We, however, consider that the presence or absence of nucleotides in the binding site would profoundly affect the structural fluctuation, which may play a key role in the protein function. Therefore, it is important to perform simulations including the nucleotide.

We introduce the dual Gō model, a variant of the C_α Gō-like model 24,25. In this model, a protein chain consists of spherical beads that represent C_α atoms of amino-acid residues connected by virtual bonds. While only the native structure is taken as a reference structure for the potential energy function in the conventional Gō-like models, the dual Gō model takes two reference structures, structure 1 and structure 2, in the effective potential.

The effective protein energy U_p at a conformation Γ is given as

(1)

(2)

(3)

(4)

For potential functions, and V_nnc, we use the same functions as Clementi et al. 24,

(5)

(6)

(7)

(8)

(9)

Display large version of this figure

Figure 1

Gō-potential: the solid line is dual Gō-potential energy profile for the ij pair. C₁₂=U₂/U₁.

If a residue pair ij is a native-contact pair in structure 12 but not in structure 21, the interaction potential between the i^th and j^th residues is a Lennard-Jones potential, with a single minimum. Other relevant parameters are k_b=100.0, k_θ=20.0, k_ϕ=1.0, k_nc=k_nnc=0.25, and C=4.0. The cutoff length for calculating is taken to be If a residue is included in the gap region of one of the reference structures, a potential energy function that has only a single minimum corresponding to the other reference structure is assigned. That is, a potential function for such a residue is the same as that of the standard Gō-like model.

To understand the mechanism of the actomyosin motor, it is desirable to study the conformational change to the rigor state. However, no x-ray crystal structure is currently available for a true rigor complex with actin; therefore, we study the structural relaxation from the pre-power stroke state to the near-rigor state. Structures 1 and 2 thus correspond to the near-rigor and pre-power stroke structures, respectively. We use 1Q5G 26, which is the nucleotide-free structure of Dictyostelium discoideum myosin II, as the near-rigor state. Although a few nucleotide-free structures of myosin II have been determined, only 1Q5G is regarded to be the near-rigor state, because both switch I and switch II are in the open position. We choose 1VOM 27 for the pre-power stroke structure. It includes an ADP·P_i analog (ADP·VO₄) in the nucleotide binding site. While 1Q5G consists of residues 2–765 without a gap region, 1VOM includes only residues 2–747 and has gap regions where the structure has not been determined by x-ray crystallography. Therefore, we use the structures of residues 2–747 for simulations. The residues in the gap regions of 1VOM are set randomly in the initial conformation and with self-avoiding conditions being satisfied under the condition that the bond length is 3.8Å. (See Fig. 2.)

Display large version of this figure

Figure 2

We chose (a) near-rigor structure, 1Q5G, and (b) pre-power stroke structure, 1VOM for structure 1 and 2, respectively. 1VOM contains the ADP·P_i analog, ADP·VO₄ (yellow). Also shown are the N-terminal (green), 50-kDa subdomain (red) and the converter (cyan) included in C-terminal subdomain (blue) that is connected to the lever arm.

The nucleotide molecule is also included in the simulation explicitly as a coarse-grained chain (Fig. 3). The coarse-grained nucleotide is represented as a short linear chain of five beads, corresponding to a purine base, a sugar (ribose), two phosphates, and VO₄ (a phosphate analog), respectively. Explicit treatment of nucleotide molecule as a short chain makes it possible for us to investigate the relation between nucleotide dissociation and structural change of the myosin motor domain.

Display large version of this figure

Figure 3

(a) ATP and (b) coarse-grained ATP.

The intranucleotide interaction is defined as

(10)

(11)

We introduce an index to characterize the state of the nucleotide binding, Q_nucl(Γ); we count how many of the nucleotide contacts formed in structure 2 (the pre-power stroke) remain in a given conformation, Γ. Then, Q_nucl(Γ) is defined as this number divided by the number of nucleotide contacts in structure 2. Q_nucl ∼ 1 when a nucleotide is bound, and Q_nucl=0 if the nucleotide-binding site is empty.

The total energy of the system, U_tot is the sum of three terms:

(12)

(13)

For a given protein conformation, Γ, we note that the native contact of structure 1 (or 2) between i and j is formed if the C_α distance r_ij=|r_ij| satisfies

Simulations are started from the pre-power-stroke structure. The initial velocities of each bead is given to satisfy the Maxwell distribution. The temperature is set lower than the folding temperature for structure 1.

We also run an implicit nucleotide simulation, in which the nucleotide is not explicitly included. Instead, the relative distance of nucleotide-contact residues are constrained by virtual bonds in all-to-all connection to keep the pre-power-stroke form. The natural length of the virtual bonds i-j, is the C_α distance between the i^th and j^th residues in the pre-power stroke conformation (structure 2). The total effective potential energy, U_tot is U_p+U_con, and

(14)

The dRMSD from the near-rigor structure is defined as

(15)

Typical time courses of dRMSD are shown in Fig. 4 for (a) k_p-n=0.6 and (b) 0.7. At the initial conformation, the dRMSD is ∼3.9Å, and decreases rapidly to ∼3Å, and stays there for a while. After that, the conformation abruptly relaxes into the final state, dRMSD ∼1.5Å. This final state is actually the near-rigor state, judging from its average structure; in fact, the dRMSD of the average structure is ∼1Å. In short, the myosin motor domain in the pre-power stroke conformation relaxes at first into the intermediate state and then to the near-rigor conformation. The time courses of Q_nucl are also shown in Fig. 4. The structural relaxation occurs after, or at the earliest, at the same time, as the nucleotide dissociation. Furthermore, the relaxation tends to synchronize with the dissociation as k_p-n increases.

Display large version of this figure

Figure 4

Relaxation time courses of dRMSD (from 1Q5G) (solid line) and Q_nucl. (dotted line) for five trajectories are shown. Different colors distinguish different runs. (a) k_p-n=0.6; (b) k_p-n=0.7.

The largest difference between the intermediate state and the initial (pre-power stroke) conformation is the position of the converter domain with respect to the other subdomains; the relative positions among other subdomains (for example, the N-terminal and the 50-kDa subdomain) are similar to those in the pre-power stroke conformation (see Supplementary Material Figs. 1 and 2 ). The average dRMSD of the intermediate state varies slightly with the parameter k_p-n (Fig. 5), reflecting little difference in the position of the converter relative to the other subdomains.

Display large version of this figure

Figure 5

Time sequence of dRMSD. The red line is the trajectory of the no-constraint simulation, and the other lines are trajectories of simulations in which nucleotide-binding site are constrained.

To discuss cooperativity of the final relaxation, we defined the contact formation rate, P^ij(d_L, d_H), as a probability that the native-contact in structure 1 between i and j residues is formed for d_L<dRMSD<d_H. Fig. 6 shows location of the native pairs that make contacts at the final relaxation. The residues involved in pairs that the contact is not formed in the intermediate state, P^ij(2Å, 3Å)<0.3, and is formed in the final state, P^ij(0Å, 2Å)>0.7, at k_p-n=0.7, as indicated by beads. The figure shows that these relaxation-associated pairs are concentrated at the boundary between N-terminal and 50-kDa subdomain, that is, the region around the nucleotide binding site. Thus, the final relaxation process consists of a rearrangement of the N-terminal subdomain against the other part. For the relaxation-associated contacts, we calculated the dRMSD dependence of the contact formation rate at intervals of 0.1Å( P^ij(d, d+Δd) with 1<d<4(Å) and Δd=0.1(Å)). The result is shown in Fig. 7. In this resolution, virtually all these contacts form at the same dRMSD, and thus simultaneously in time. We found that the native-contact formation takes place cooperatively at the final relaxation irrespective of k_p-n (see also Supplementary Material Fig. 3 ).

Display large version of this figure

Figure 6

The relaxation-associated contacts satisfying the conditions P^ij(2Å, 3Å)<0.3 and P^ij(0Å, 2Å)>0.7 are shown. Residues included in the contacts are represented by spheres; residue numbers are indicated at the left. Green, red, and blue correspond to N-terminal, 50-kDa, and C-terminal subdomain, respectively.

Display large version of this figure

Figure 7

dRMSD dependence of the contact formation rate of the relaxation-associated contacts at k_p-n=0.7. Residue numbers of the contacts are indicated at the left.

To clarify the k_p-n dependency of the synchronization, we plot the histograms of τ_d and of Δτ=τ_d–τ_r from 500 independent runs for each value of k_p-n, where τ_d is the number of steps taken before the nucleotide dissociates, τ_r is the number of steps taken before the conformation relaxes to the near-rigor state, and Δτ is the delay in the relaxation after the dissociation takes place (Fig. 8). For small k_p-n, both the histograms of τ_d and of Δτ show exponential decay; thus, the nucleotide dissociation and the relaxation of the myosin conformation are considered to be decoupled. For larger k_p-n, on the other hand, the histogram of τ_d cannot be fitted to an exponential decay. In addition, the average of τ_d is shifted to the right and the delays Δτ become shorter; in other words, the nucleotide stays longer at the binding site and the conformational relaxation tends to occur immediately after the nucleotide dissociation. For k_p-n=0.7, dissociation and relaxation occur nearly simultaneously in >70% of 500 trajectories. Note that apparent Δτ<0 cases are caused simply from the numerical ambiguity of τ_d and τ_r and actually correspond to coincidental dissociation-relaxation.

Display large version of this figure

Figure 8

Histogram of (a) the number of steps before dissociation and (b) the delay of the relaxation after the dissociation from 500 independent runs for each k_p-n.

As already mentioned, the myosin motor domain relaxes to the near-rigor conformation only after the dissociation of the nucleotide and not before. Thus, it seems that the nucleotide must be unbound for the final relaxation to occur. This observation leads to a speculation that the nucleotide blocks the deformation of myosin around the nucleotide-binding site by its volume. To investigate the case of when the nucleotide cannot dissociate, we attempt an implicit nucleotide simulation. In this simulation, instead of treating the nucleotide molecule explicitly, we connect the residues that would contact with the nucleotide by virtual bonds, to force the nucleotide-binding site to keep the pre-power stroke form. Fig. 5 shows the time courses obtained by the simulations. We find that the conformation remains at the intermediate state and that the relaxation toward the near-rigor state is not completed.

We have shown that the myosin motor domain does not relax to the near-rigor conformation before the nucleotide dissociates. Ishijima et al. 29 showed by simultaneous observation of ADP release and mechanical events that force is generated at the same time as or several hundreds of milliseconds after the dissociation of ADP. Our results are consistent with their experimental findings if force generation is preceded by structural relaxation. Moreover, the results from the simulations in which the conformation of the nucleotide-binding site is constrained also indicate that the relaxation is indeed prevented if the nucleotide cannot dissociate.

Based on these observations, we now suggest that the primary role of the nucleotide in the power stroke process is to suppress relaxation through blocking deformation around the nucleotide-binding site by its volume. In this scenario, hydrolysis is required to alter the affinity of the nucleotide to the binding site.

Our simulations have also shown that the structural relaxation is synchronous with nucleotide dissociation when the nucleotide is tightly bound to the myosin head. In other words, the nucleotide dissociates cooperatively with the motion of the subdomain. This strong coupling of deformation and dissociation seems to be relevant to the function of the strain sensor, in which nucleotide dissociation is controlled by the strain induced by an external force. The correlation depends on the binding strength, k_p-n; the relaxation is only loosely coupled with the dissociation for weak binding conditions. A large kinetic diversity among myosins observed by experiment 4,30 may be attributed to this binding strength dependence.

The intermediate state observed in the relaxation process should also be discussed. Although several intermediate states have been revealed by kinetic experiments 4, their structural aspects, except for ADP·P_i state, are little known. Shih et al. 31 reported, from their FRET study, that there are two pre-power stroke conformations; while one conformation corresponds to the crystal structure of the complex with the ADP·P_i analog, the other conformation has not yet been observed using crystallography. We found that the average structure of the intermediate state observed in this study is consistent with the latter conformation.

Recently, another process in a reaction cycle, the recovery stroke, which is the process from ATP-bound structure to pre-power stroke structure, was studied using the conjugate peak refinement method of all-atom model 32. In this method, structural change is described as an adiabatic process. Although it is suitable to investigate the detail of the side-chain interaction, it cannot treat dynamics of spontaneous nucleotide dissociation (our main interest) or distribution of dwell time.

In this work, actin filament was not included in the simulation. It is in fact known that the rate of the nucleotide release and of the structural change are affected by actin. But it is also known that the myosin motor domain alone undergoes structural changes depending on status of the nucleotide binding site. Thus, we regard that the relation between the nucleotide release and the structural change can, at least, be captured qualitatively without actin. Detailed study on effect of actin will be left for a future study.

We did not treat the releases of ADP and P_i separately. Although it is believed that a powerstroke is associated with the P_i release, as far as we are aware there is no direct observation that the P_i release is tightly linked to the myosin structural change or the force-generation. Since the P_i release is known as the rate-limiting step for myosin II, and ADP dissociates immediately after P_i is released 4,33, we suppose that the releases of ADP and P_i can be treated as a single nucleotide-release event in our coarse-grained model.

Some attempts to extend the Gō-like model to treat the structural change have been made. In the model of Hyeon et al. 34, the energy potential is gradually switched from that of the initial state to that of the final state. During the switching process, the potential function has a minimum at an intermediate structure of the initial and the final structure. The approach of Zuckerman 35 is similar to our model in the sense that each interaction potential function contains information on two structures. He introduced a double-square-well potential for the native-contact pairs into a lattice protein model. Okazaki et al. 36 proposed another type of model in which total effective energy has a double-basin. These types of models, in which two conformations are embedded in an energy potential function, will be useful in understanding the dynamics of protein conformational change.