A quantitative molecular model

The previous lecture in this topic (‘The molecular mechanism of muscle’) introduced the actomyosin molecular motor (Slide 1) and showed how it can be modelled using the formalism of diffusion and chemical reactions. This lecture starts with the quantitative model introduced there. This is solved analytically in the special case of constant velocity. In general, such models need to be solved numerically. The lecture concludes with a discussion of one numerical approach to modelling, using Monte Carlo simulations.

2.1 A quantitative model of acto-myosin

The previous lecture introduced a quantitative model of actomyosin. Slide 2 (which is a repeat of Slide 12 from the previous lecture) lists reaction rate constants and free energies assembled from experimental data, plus some educated guesswork. These rates are averages, measured under a variety of physiologically relevant conditions. The experiments best represent the equilibrium position (minimum) of each state.

States in the mechanochemical cycle are defined by whether actin is bound to myosin, and which if any of ATP and its hydrolysis products are bound, which in turn determines the most stable shape of the myosin molecule. Typical rate constants for transitions between the states are given, as well as free energies relative to the reference state labelled 2.

Note that binding is the slowest forward step (highlighted by the red circle): each head spends most of its time unbound – idling, while other heads pull. This motor is described as having a ‘low duty ratio’.

2.1.1 Energy profiles

Slide 3 shows the energies involved in the actomyosin cycle. Each chemical state has a distinct free energy profile as a function of the mechanical coordinate. States are labelled with the same letters and numbers as in Slide 2. The energies of unbound states, and energy minima of bound states are taken from the table in the previous slide. (We have included both a static slide and a video of an animated slide.)

The position-dependent energies for bound states are approximated by parabolae – that is, it is assumed that the protein acts as a linear spring. The power stroke is built in by incorporating a shift in the equilibrium position between states 4 and 5 (the transition that accompanies phosphate release, described in Lecture 1, slides 8, 9 and 11). The spring constants are unknown combinations of short-range interactions between atoms, long-range electrostatic forces and entropic forces. (Force = –∂G/∂x, an extension of the normal definition –∂U/∂x to include entropic forces.) The same spring constant is assumed for all states; it is adjusted such that the transition at the beginning of the power stroke (4→5) is close to neutral in free energy. If ΔG45 were more negative, the transition would waste energy as heat; if ΔG45 were positive, then the probability of phosphate rebinding before the power stroke would become significant.

N.B. states 2 and 2´ are identical, except that they are separated by one cycle and thus the hydrolysis of one ATP molecule.

2.1.2 Rates

Also in Slide 3, actin binding (3→4) and phosphate release (4→5), which is followed by the power stroke, are only allowed in a narrow range of positions near the relaxed position of the unbound motor (through the gate).

Release is only allowed after the power stroke is completed. (It could be achieved by controlling any of the cascade of transitions 5→6→1→2.)

Note that an external force could be included in the model by adding a work term Fx to the free energies of all states. This would raise the energy minima of bound states 4,5,6,1 and move them to the left, shortening the power stroke.

Detailed balance requires that all reactions are reversible. To make a reaction go mainly one way, a free energy drop is needed (the system must be out of equilibrium). This mechanochemical cycle is driven by the chemical free energy of ATP hydrolysis. Small free energy drops ‘waste’ free energy as heat but ensure that the cycle proceeds forwards.

Unbinding from the actin filament after the power stroke is accompanied by a significant free energy drop. This makes it unlikely that rebinding will occur at positions where the crossbridge would produce a negative force. As soon as the power stroke is over, rapid and irreversible release is desirable.

2.2 Reaction–diffusion equations

2.2.1 Setting up the equations

The next step is to find a formalism that will allow quantitative predictions to be drawn from the energy profiles of Slide 3. We can write down a set of reaction–diffusion equations (Slide 4) that describe the evolution of the probability density that describes the state of the system (cf the derivation of the Einstein relation in Lecture 2 ‘Brownian Motion’ in the Statistical Thermodynamics topic).

Let P(x,t)dx be the probability of finding a motor in the range x→x+dx at time t, and consider diffusion in the presence of an external force. We can write an expression for the probability flux J:

Equation 1
Equation 1


Equation 2
Equation 2

and D and γ: are diffusion and drag coefficients, respectively.

(Note that the Einstein relation can be derived very simply from Equation 1 by considering equilibrium, where P is proportional to exp(–U/ kBT) and J = 0. This derivation is set in one of the tutorial problems accompanying this lecture.)

The model is constrained by requiring that total probability density is conserved:

Equation 3
Equation 3

Hence we can write down the Fokker–Planck equation, which describes the evolution of the probability density for a particular state, without any chemical transitions:

Equation 4
Equation 4

Transitions are assumed instantaneous, so they are modelled by rate constants that are functions of position (as in slides 2 and 3). So, including the possibility of transitions between states, we have:

Equation 5
Equation 5

Equation 5 is a reaction–diffusion equation. The first two terms are ‘mechanical’ and describe motor movement while in state i, while the third term is ‘chemical’ and describes transitions in and out of state i at position x.

2.2.2 Solving the reaction–diffusion equations for actomyosin

To simplify the model of Slide 3 yet further, we can ignore all transitions except those that involve mechanical motions. The model has only two states, with myosin bound to actin or not (free). Similar to states 2 and 2´ in Slide 3, states free and free´ are identical except that they are separated by one cycle and thus the hydrolysis of one ATP molecule. Transitions retained in Slide 5 are:

  • binding to the actin filament/phosphate release;
  • power stroke;
  • unbinding.

  • Other elements of the model are:

  • binding is allowed only around x = –5 nm (5 nm is the length of the power stroke);
  • unbinding is allowed only after the power stroke is finished (x > 0).
  • These two constraints are realistic and are imposed by the structure of the myosin motor and the actin–myosin interface. Note that binding and unbinding are not simply the reverse of each other but follow different pathways and are coupled to different chemical reactions. The arrows in Slide 5 (right) illustrate this. The cycle hydrolyses one molecule of ATP as in the full model of Slide 3.

    Note also that the reverse reactions, represented by the dashed lines in Slide 5 (right), cannot have zero rate constants without violating the principle of detailed balance (see Lecture 1 of this topic). However, the free energy of ATP hydrolysis biases the cycle in the forward direction, allowing us to neglect the backward reactions.

    As shown in Slide 6, we can use this minimal model to obtain an analytical force–velocity relationship in the special case that the muscle is contracting at constant velocity. This is a very good approximation in a muscle, where many myosin heads are connected in parallel, averaging out the stochastic behaviour of individual heads.

    The motion of the particular myosin molecule that we are considering is imposed by the collective behaviour of all motors attached to the same thick filament – the appropriate diffusion coefficient is therefore that of the thick filament, which we will take to be effectively zero:

    To find the average force during steady muscle contraction at a constant speed, we need a time-independent (‘steady-state’) solution:

    Before we start, we can deduce that the muscle will generate less force at higher speeds for two reasons:

  • the probability of binding each time the head passes the ‘target site’ is proportional to the time that it spends within the target range, which is less at higher speed;
  • the drag exerted by the crossbridge after it has gone past the energy minimum but before it has released will be more at higher speed because the time to release is independent of speed and the drag force increases with distance.
  • In the special case of the simplified model of Slide 5, at constant speed, Equation 5 reduces to a single ordinary differential equation, with solution P(x), as follows:

    Constant velocity


    so that x and t are no longer independent variables, and we can write

    Equation 6
    Equation 6

    2 states

    We define

    Equation 7
    Equation 7

    Strictly, P(x) is a probability density, with dimensions of inverse length. However, constant velocity means that the motor spends equal time at all values of x within one period of the actin filament’s repeating structure, so that

    Equation 8
    Equation 8

    We are free to set this constant to 1 for convenience. Doing this makes P(x) equal to the probability that myosin is bound to actin at any given position x. An alternative convention that the integral of (Pfree(x)+ Pbound(x) )dx over the period, d, of the actin filament should equal 1 would require us to set our constant equal to 1/d rather than to 1. (Here the period, d, is the distance between adjacent binding sites.)

    2.2.3 Steady state

    For a steady-state solution to Equation 5 we set


    Equation 9
    Equation 9

    Solving for P(x) is part of the problem set that accompanies these lectures. The reaction–diffusion equation is solved piecewise (Slide 7): the probability of binding increases from zero as the motor crosses the target zone; unbinding occurs after the power stroke is complete.

    The three pieces of the solution are:

    Binding zone around x = x0 (= –5 nm).

    Assume initial condition P(x0) = 0 (binding zones are widely separated).

    Power stroke zone x0 <x < 0.

    Drag-stroke zone x > 0, rate constant koff.

    Each section is described by a simple first-order differential equation with exponential relaxation as its solution; the boundary condition is that the probability density is continuous. Even this extremely simple model gives a reasonably accurate prediction of the measured force–velocity relationship for muscle.

    If the stiffness of each crossbridge is κ, and d is the distance between binding sites, the average force exerted per crossbridge is:

    Equation 10
    Equation 10

    The graph in Slide 8 shows that the simple model for actomyosin gives a good fit to the data.

    2.2.4 Other molecular motors

    Reaction–diffusion equations for other molecular motors are solved using a variety of numerical methods. Typically, the spatial variables are also discretised, so the free energy landscape is reduced to a set of discrete transitions between states. The steady-state approximation is usually made, in which case the numerical solution consists of a set of probabilities for the motor to be found in a particular state at a particular mechanical position, averaged over long periods under unchanging conditions.

    Advantage: average properties obtained directly.

    Disadvantage: no microscopic information (e.g. distributions of step times).

    2.3 Monte Carlo simulation

    An alternative approach to modelling molecular motors is to use Monte Carlo simulations, in which representative motor trajectories (chemical state and position versus time) that have correct probability distributions for transitions are generated from sequences of random numbers.

    2.3.1 Chemical transitions

    Slide 9 defines parameters for Monte Carlo simulations of chemical transitions. At each time step (Slide 10), random numbers are used to determine the chemical state of the motor at the next step based on transition probabilities appropriate to the current state of the motor. Rate constants for chemical transitions are input parameters, as before, and must obey detailed balance (as before, though the minimal model was simplified by setting small reverse reaction rates to zero).

    The time step must be small compared with the reaction time so that it is reasonable to assume that the motor is in one particular chemical state during the time step.

    At time t, the motor is in state i at position x. State i has chemical transitions to states j, j´, etc. Rate constants for these transitions ensure detailed balance (cf Equation 2 in the previous lecture):

    Equation 11
    Equation 11

    The probability, P, of making a transition from state i to state j during a time step Δt is

    Equation 12
    Equation 12

    For sufficiently small Δt:

    Equation 13
    Equation 13

    A random number between 0 and 1 decides whether any transition should occur in time Δt and, if so, which.

    Slide 10 shows just one way to use a random number to model chemical transitions, with a constant time step. There are other methods, including some where the time step is variable.

    2.3.2 Motion

    Finally we can show how the Monte Carlo technique can simulate motion. The motor is described by an instantaneous position and velocity. Potentials U(x) and drag coefficient γ are the input parameters to the model (Slide 11). The time step must be small enough that the motor does not move a significant distance (compared with the length scales of the energy profiles and rate-constant functions), so that it is reasonable to consider the motor to be at one position during the step.

    The Langevin Equation is the normal equation of motion for a particle of mass m and drag coefficient γ, but with an added thermal noise force (the second term) Fi:

    Equation 14
    Equation 14

    Using the Einstein relation, it is possible to work out the power spectrum of the thermal force F(t). It is ‘white noise’ with a uniform power density of 4kBTγ (see Howard Chapter 4).

    The inertial term is negligible. (Reynold’s number, R = VL /ν, is a measure of the importance of inertial relative to viscous terms, where V and L are characteristic speed and length scales, and ν is the kinematic viscosity. For water, ν = 10–6 m2s–1; a typical sarcomere contraction speed is of the order of 1 μm s–1; the sizes of a myosin head and a sarcomere are of the order of 10 nm and 1 μm, respectively. The Reynolds number is therefore in the range
    10–6 ~10–8. Equation 14 thus reduces to

    Equation 15
    Equation 15

    In Slide 12, integrating over one time step gives two terms in the displacement: deterministic movement due to the motor potential; and random motion due to the thermal force. The second term is simply free diffusion:

    Equation 16
    Equation 16

    leading to

    Equation 17
    Equation 17

    where N(0, 2DΔt) is a random number drawn from a distribution with mean 0 and variance 2DΔt.

    In Equation 17 the first term describes deterministic motion due to a potential, while the second describes Brownian (stochastic) motion. The second term in Equation 17a is the distance the particle would have travelled by free diffusion in time Δt; the simulation substitutes a random distance drawn from a Gaussian distribution with mean square 2DΔt.

    The output of a Monte Carlo simulation is a trajectory of the position and chemical state of the motor as functions of time – sometimes called Brownian dynamics simulation.

    Advantage: individual runs reveal the microscopic behaviour of the motor: step sizes, distribution times for intervals between steps, etc.

    Disadvantage: many or long runs are needed to determine average properties, such as motor speed, ion flux, dependence upon chemical driving force or external load. This can be computer intensive (but still far less so than molecular dynamics simulations).