From Single Molecule Fluctuations to Muscle Contraction: A Brownian Model of A.F. Huxley's Hypotheses

Muscular force generation in response to external stimuli is the result of thermally fluctuating, cyclical interactions between myosin and actin, which together form the actomyosin complex. Normally, these fluctuations are modelled using transition rate functions that are based on muscle fiber behaviour, in a phenomenological fashion. However, such a basis reduces the predictive power of these models. As an alternative, we propose a model which uses direct single molecule observations of actomyosin fluctuations reported in the literature. We precisely estimate the actomyosin potential bias and use diffusion theory to obtain a Brownian ratchet model that reproduces the complete cross-bridge cycle. The model is validated by simulating several macroscopic experimental conditions, while its interpretation is compatible with two different force-generating scenarios.


Introduction
Forty years ago, A.F. Huxley and R.M. Simmons [1] reported the behaviour of skeletal muscle fiber by observing near instantaneous force changes following fast and small length perturbations. Accordingly, Huxley and Simmons proposed that the force is generated from an actomyosin complex that exists in (at least) two stable but different force-generating conformations, or states, that oscillate between one another in response to external perturbations. This theory has since been analysed by several studies examining muscle structure [2][3][4], macroscopic muscle fiber behaviour [5][6][7][8][9], and single molecule experiments (SME) [10][11][12][13] including recent SME studies that have directly observed the oscillations [14,15]. In their same seminal paper, Huxley and Simmons proposed a mathematical model that described an explicit relationship between the probability of the actomyosin state and the external perturbation applied. That model has since been refined by several works and extended to the entire crossbridge cycle [5,[16][17][18][19][20].
As informative as these contemporary models are, they are limited due to the lack of direct information on the state transition dynamics. While they introduce an explicit definition of the energy in each state, they define the transition between states as rate functions with ad-hoc dependencies on the chemical reaction coordinate or base the transitions on hidden assumptions of the actomyosin potential energy. In other words, actomyosin properties are deduced by fitting the macroscopic behaviour of the muscle fiber even though the goal of the model is to interpret the muscle behaviour from the actomyosin properties. We believe that this approach reduces the predictive power of such models because the oscillatory behaviour of the actomyosin complex is a fundamental property of muscle function such that its description cannot be restricted to phenomenological rate functions. Therefore, we here describe an alternate model based on diffusion that considers the entire Lymn-Taylor cycle (Fig. 1). Although a diffusional approach for the attachment-detachment process is not feasible because of insufficient experimental data for skeletal myosin, configuration changes in the actomyosin complex can be described by thermal diffusion if the periodic potential is properly designed for there to be a clear relationship between the chemical configurations and the external force applied. Thus, we approach the problem by defining a diffusional model using SME data [14,15] and validate our hypotheses by simulating experimental macroscopic muscle fiber data. This new approach (see also [21][22][23]), and its future application to the attachment-detachment process, aims to reduce the degree of freedom in the parameters to more accurately explain muscle behaviour by using SME data.

Cross-bridge Cycle and Model Description
In its most basic form, the cross-bridge cycle can be described as a four-state model with two detached and two attached states (upper and lower parts of Fig. 1, respectively). The left side of figure 1 shows a detached myosin head that initially has a low affinity for the actin filament (M-ATP). The myosin next hydrolyses ATP to generate a high energetic state (M-ADP-Pi). Brownian fluctuations then allow the myosin to bind to the preferred actin-binding position to form the actomyosin complex AM-ADP-Pi. This attachment is probably driven by the ''Brownian search-and-catch mechanism'', which assumes that the affinity for actin increases when the myosin elastic element is stretched forward [24,25]. Force generation stretches the elastic element even more so that the complex takes the AM-ADP state. Detachment occurs when a new ATP molecule substitutes for the exhausted ADP molecule (M-ATP) to relax the elastic element and begin a new cycle.
Oscillating sub-steps during the force generating (AM-ADP) state have been observed experimentally by attaching a single myosin II molecule to a large micro-needle and associating it with an actin bundle [14,15]. This results in several 5:5 nm steps per ATP cycle that are biased in one direction, as predicted in [1]. The steps are conventionally explained by the lever arm theory: the long portion of the myosin head rotates to slide the actin and myosin filaments past each other while the globular portion remains firmly attached to the actin monomer (Scenario 1 in Fig. 1) [1,17,26,27]. Despite that, the constant size and frequency of the steps observed experimentally may be generated by a different mechanism where the helical shape of the actin filament drives the sliding of the myosin head on different actin monomers to pull the thick filament (Scenario 2 in Fig. 1) [14,15] and the rotation of the lever arm takes place only at the end of sliding during Pi release. This scenario has been recently supported by MD simulations of skeletal actomyosin in the rigor configuration [28]. In that work, the preferential direction of thermal fluctuations in the actomyosin complex was shown to arise without ATP consumption.
In both scenarios, the steps are generated by stable actomyosin states which mathematically correspond to several minima in the potential energy, E a (x), which describes the actomyosin interaction ( Fig. 1, right). Thus, our model does not discriminate between the two scenarios and defines the energy potential using only SME data. Force is generated when the actomyosin complex fluctuates between stable states in a manner that allows the complex to adapt efficiently to different external stimuli. The key feature of the model is that mathematically these fluctuations are described as diffusional motions inside E a (x). Therefore, we do not need AM-ADP-Pi to AM-ADP transition rate functions, nor do we need to consider the dependence of these functions on the head position when fitting the macroscopic behaviour of a contracting muscle. On the contrary, having E a (x) defined by SME data alone means we can use macroscopic behaviour to test the goodness of the shape of the potential.
To model the macroscopic response (see Text S1 for method), we assume N myosin heads are in parallel in a one-half sarcomere, and that several sarcomeres are in series and behave uniformly. The myosin heads are attached through an asymmetric elastic element [29] to a rigid backbone (thick filament, right side of Fig. 1) at uniformly spaced locations due to the incommensurability of the periodicity of the two filaments [30]. In order to maintain a clear chemo-mechanical relationship, we define a minimalistic model where the myosin can exist in only two different states, one force generating and one non-force generating (zero on average), which roughly relates to the attached and detached states, respectively. Transitions between these states are regulated by the stochastic variable v (pure number), which fluctuates between the values 0 and 1.
In the detached state (v~0), the myosin head is subjected to unbiased thermal fluctuations and to the force generated by the elastic element. The transition from M-ATP to M-ADP-Pi is a diffusion driven process that does not require a rate function. When myosin interacts with the actin filament (v~1), however, it is also subjected to the force generated by E a (x). Each myosin molecule is considered an over-damped particle whose geometry is completely defined by its drag coefficient, g x~7 0 pNns=nm [31]. The environment in which the myosin moves is defined by the Boltzmann constant, k B , and the absolute temperature, h~297K. The motion of the i-th myosin head leads to a flashing ratchet system with its dynamics described by the Langevin equation: In the force generating step, the actomyosin complex oscillates between stable states either by rotating the lever arm domain (scenario 1) or by sliding the myosin head along the actin filament (scenario 2). Right side: Mechanical representation of the half-sarcomere as many parallel myosin heads. In the detached state, thermal fluctuations by a myosin head are constrained only by an asymmetric elastic element [29]; the absence of other interactions is represented by the flat energy landscape. In the attached state, the head also experiences an actomyosin complex energy landscape, E a (x), which has a periodicity of 36 nm and four minima or stable states (see text where E e~1 2 K(x i ,l i )(x i {l i ) 2 is the asymmetric [29] elastic energy of the myosin globular portion, lever domain, and long tail in series. Stiffness is defined as [17] and l i is the reference position on the backbone of the i-th head. By setting l i (t)~z(t)ziL=N, where z(t) is the position of the backbone, continuous behaviour of the myofibril is ensured by having l uniformly distributed over ½0{L (L~36 nm). vC(t)w is a random term constrained by vC(t)w~0 and vC(t 1 )C(t 2 )w~2d(t 1 {t 2 ) (white noise). Total force generated by the myosin heads is given by modeled as a piecewise linear, multi-stable potential with minima equally spaced by a distance that is compatible with the actin monomer diameter, d, and a directional bias defined by DG. Although the potential is locally biased, it is flat on average, repeating itself every L, with each period containing six actin monomers (d~6 nm).
Experimentally, the number of steps ranged from one to five at low micro-needle stiffness (K n~0 :01{0:03 pN=nm [14]) and one to four at high micro-needle stiffness (K n~0 :1{0:6 pN=nm [15]). We impose that the myosin can interact with only four actin monomers, and attachment cannot occur in the region nLz½0{3d|½5:5d{6d, where n is an integer, a constraint that mimics steric effects in the actin double helix.
Since the potential is globally flat, net movement in the backbone will only occur by breaking the global equilibrium [23]. We introduce a jump process that constrains v i (x) to follow the rate functions k z (x) and k { (x) (described in Fig. 1) out of balanced equilibrium by setting k { =k z exp½E a (x)=k B h. Referring to the Brownian search-and-catch mechanism originally proposed by A.F. Huxley in 1957 [25] and seen in Myosin VI [24], we introduce an attachment process where a detached myosin head searches for its preferred actin-binding site in the forward direction by thermal fluctuations. The detachment rate function k { (x) we use is a slightly modified version from the one proposed in [25] in order to consider the more detailed geometry of the energy potential.

SME Simulations
To quantitatively define the energy bias for muscle in vivo, we simulated in vitro data from [14] and [15] (Fig. S1). Previously, a bias of DG~2{3 k B h is obtained in [15] assuming the microneedle applies a constant force, F , to the molecular motor during a jump and N f =N b~e xp½(DG{Fd)=k B h (N f and N b are the number of forward and backward jumps, respectively). Nevertheless, the drag coefficient, g n , of the micro-needle was three orders of magnitude higher than that of the molecule, meaning the external system applies a force that varies with the instantaneous position of the myosin during a jump. Therefore, the previous estimation may underestimate the real DG. Consequently, we numerically explore the effect of DG on N f =N b by simulating the micro-needle displacement for several DG and by comparing the resulting N f =N b with the experimental data (Fig. 2). Because in those experiments only the S-1 portion of the myosin head was used, we imposed a symmetric molecular stiffness K~2 pN=nm [17]. The micro-needle is also represented as an over-damped particle whose drag coefficient is g y~1 0 3 g x [32] and is maintained near its original position by stiffness K n [15].
The set-up is described by the following system of stochastic differential equations: where E a is the actomyosin potential shown in Fig. 1. Since the ratio of forward and backward jumps were analysed in [15], we impose the mean value of needle stiffness used in that experiment K n~0 :3 pN=nm. Three parameters are needed to define the potential: the energetic barrier, DG B , between two minima; the asymmetry of the potential, l, which is the ratio of the distance between an adjacent maximum and minimum and d; and the bias or energy difference, DG, between two minima. We impose an energy barrier of DG B~1 6k B h, which is compatible with myosin II conformational changes that occur in less than one millisecond in the unloaded condition [33], as obtained by the exact solution of the first passage time, t f , of a particle against a constant force F a~( DG B {DG)=ld [31]. During a change of state in a multi-stable energy potential, the particle spends most of its time surmounting the energy barrier [34]. We suppose that the maximum of the barrier between two minima is displaced toward the forward minimum (lv0:5), as in [1] where this distance (d 2 in Fig. S1) is zero. Numerically we cannot impose a zero distance, thus we do the following reasoning. The distance between the steps observed in [14,15] was 5:5 nm. This should equal the distance between two stable states, or minima in the energy potential. Assuming d~6 nm, we set l~0:1 so that the distance between a minimum and the subsequent Figure 2. Bias of actomyosin energy. The diffusion of a particle linked to an external micro-needle in a periodic potential tilted by a constant bias, DG=d, is simulated, and the ratio of forward, N f , and backward, N b , jumps for different DG is reproduced and compared to the experimental data. Red squares, DG~3kh; green triangles, DG~5kh; and blue rumble, DG~7kh. Black dots, experimental data from [15]. Continuous lines are the corresponding exponential fits. Simulations show that a bias higher than DG~3kh corresponds to the observed ratio. A value of DG~5kh is chosen as most probable due to the boundary conditions (see text). Insert: Experimental oscillating substeps as expected from [1] during the attached state (see fig. S2). Figure adapted from [15]. doi:10.1371/journal.pone.0040042.g002 maximum (d 1 in Fig. S1), the portion of the potential that influences the step, is (1{l) Ã d~5:4 nm. Then we simulated eq. (2) for different external forces based on the micro-needle position and different values of DG. A typical simulated trace is seen in Fig.  S2. This is comparable with the experimental trace shown in the inset of Fig. 2. The analysis shows that DG~3 k B h (red line) results in a N f =N b ratio that always underestimates the experimental value (grey line). On the other hand, if one assumes DG~5 k B h and F~0:1 pN, then the relation used in [15] should lead to N f =N b~e xp½(DG{Fd)=k B h~128, whereas our numerical results give a ratio one order of magnitude smaller, showing that the previous relation strongly underestimates the effect of the micro-needle size.
In order to minimize the number of parameters, we imposed some simplifications on our SME simulations. In particular, the myosin head is considered to be attached in the active state such that we ignore the rigor state. In other words, myosin always oscillates between attached stable states. This situation creates a boundary effect such that when the head reaches the last minimum it can only jump backward. Due to the bias of the potential, this means the longer the simulation, the lower the N f =N b ratio. To overcome this problem, we simulated the trajectory of each particle until the lowest minimum is reached and then interrupted the computation. In the experimental protocol of [15], jumps were measured until the rigor state occurred, which likely resulted in underestimating the ratio with respect to our method. We therefore concluded DG~5 k B h as the best estimate for the data and use this value to define the bias for the remainder of our analysis.
Summarizing the known and unknown parameters, the myosin stiffness, K, the drag coefficient, g x , and the energetic barrier, DG B , are known or accepted values from the literature, while the temperature, which defines kh, the needle stiffness, K n , and the drag coefficient, g y , are known experimental values. The energy bias, DG, is the variable under analysis and is varied over a range of realistic values. Then, the unknown or free parameters that remain are the potential asymmetry l and the distance between minima d. The product of these two values is chosen to satisfy the observed step length in the computational limit of our method (d=0).

Model Simulations
One limit of the diffusional approach is the time required to perform a real-time (up to hundreds of milliseconds for the force velocity curve) simulation of the muscle system. The reason is that for equation (1) to have any statistical meaning, the time step of the simulation must be at least two orders of magnitude less than the characteristic time scale of the system. We therefore introduced the following constraints. We set l~0:5 to increase the time step and DG B~1 1 k B h to reduce the energetic barrier maximum in the periodic potential. To evaluate the effects of these constraints, we compared the dwell times of a free particle to the two DG B values, 11 k B h and 16 k B h, finding the lower energy barrier made the dwell time approximatively 150 times shorter. Consequently, we also increased the attachment and detachment rates in a way that maintains the relative velocity to approximate those seen in experiments.
Numerically, we imposed a z~3 000 s {1 nm {1 , a {~az =12, x z~1 0 nm, x {~{ 3 nm, and k { T~1 5000 s {1 (see Fig. 1). Thus, when considering a distance ld, the maximum attachment rate is about 150 times faster than the realistic value of 67 s {1 [20]. These values result in a simulated behaviour that is faster than that seen experimentally. For example, the simulated maximum velocity (v s max ) is 80 mm=s, which is 30 to 50 times faster than the experimental value (v e max ) for frog muscle [33,35]. Thus, we must slow the simulated processes by v s max =v e max to obtain a corresponding quantity in real-time. To prove the feasibility of this approach, we simulated the rising phase curve assuming that all heads are initially detached in the isometric (z(t)~0) state, slowed the simulated curve by a factor of 40 (an average of the two v s max =v e max values), and compared it to that seen experimentally (Fig. S3), finding the two superimposed almost perfectly. Therefore, the attachment and detachment rates in our model are proportional to the velocity of contraction in a way compatible with the experimental results.
Summarizing, the protein drag coefficient, g x , its stiffness in the stretched configuration, K, and the temperature (or kh) are known or accepted values. The potential periodicity, L, the uniform random distribution of proteins, the number of monomers (six, d~6 nm), and the constraints for the four attachment sites, which relate to the helical shape of the actin filament, are either known values in the second scenario or are new experimental constraints for the first scenario. A potential bias, DG, is also imposed by the experimental data (see previous sub-section). The number of proteins, N, and the myosin filament drag coefficient, g z , are scale parameters that do not affect the normalized results below. Protein stiffness in the compressed configuration, even if qualitatively imposed in accordance with recent experimental observations, is used to fit the data and must be considered a free parameter. The potential asymmetry, l, and maxima DG B in the model simulation are also free parameters used to modify the time scale to make the calculation possible. Finally, the parameters needed to describe the attachment and detachment rates (k z and k { , in s {1 ), a z , a { , x z , x { , and k { T , are free parameters whose values are chosen in accordance with l and DG B to satisfy a coherent rising phase and maximum velocity (see above).
A comparison with other macroscopic models is difficult to make due to differences in fundamental conditions and to the wide range of simulated data. Despite that, the number of free parameters is relatively low and mostly attributed to the attachment-detachment process.

Results
We quantitatively tested our model for two classical skeletal muscle experiments: fast tension recovery after a small increment in isometric length (length clamp), and velocity of contraction against a constant load (load clamp). For the first case, a fast (microsecond) and small (few nanometers per half sarcomere) change in muscle length, d, during isometric contraction resulted in typical tension transient behaviour (Fig. 3, inset). Initially, an almost instantaneous change from the isometric tension T 0 to tension T 1 (d) is seen. This is followed by a slower (milliseconds) tension recovery that plateaus at T 2 (d). Over a longer time scale, a new population of myosins attach such that tension recovers to T 0 . In the second case, a muscle fiber bears a constant load which generates a tension TvT 0 and contracts at a constant velocity that depends on the load in an hyperbolic manner (Fig. 4) [36].
The model can simulate the length clamp behaviour by using equation (1) and assuming: which experimentally describes a length change at time t j . The force clamp behaviour is simulated by applying equation (1) and assuming where the drag coefficient of the thick filament, g z~N g x , considers the total number of motors and T~aT 0 is the external tension, with a ranging from 0 to 2.5. Simulations for both T 1 (d) and T 2 (d) curves (Fig. 3) and the V =V max vs. T=T 0 curve (Fig. 4) show very good fitting. The latter also shows qualitative fitting for the TwT 0 region.
To further investigate the model predictions, we simulated the number of attached myosins during isokinetic contraction under different external loads and compared the results with experimental data obtained from [8], where the molecular basis of the force velocity relationship was determined by X-ray interference and mechanical measurements on intact single cells. Assuming asymmetric stiffness, that study could only determine the number of stretched myosins, not all attached myosins. We simulated the average number of stretched myosins during the steady phase of an isotonic contraction at various T, finding the data could be fit excellently without any supplemental hypothesis (Fig. 5, lower). Also, the predicted mean strain per motor is in very good agreement with the experimental data at high and medium forces (Fig. 5, upper). The different behaviour at low loads can be attributed to the low number of attached motors, which may have compromised the experimental analysis.
Finally, we can conclude our model results are not due to inefficient cycling (only a small fraction of ATP energy being converted to work). Noting that each detachment process consumes one ATP molecule, we measured the efficiency of the flashing Brownian ratchet as the product of the force and the velocity of contraction divided by the total energy consumed [37].
Defining r as the number of detaching cross-bridges per unit time, we can write g eff~{ f ext v=rDG ATP . In Fig. 6, blue triangles show efficiency versus the velocity of contraction assuming    [8]. Mean strain is well fitted at high T, but overestimated at low T, where the low number of attached motors may have compromised the experimental analysis, as observed in [8]. Lower: Relative number of attached motors at different loads during the steady state phase of isotonic contraction. Simulation, triangles; experimental data, circles. Mean values+standard deviation over 11 trials are shown. Where error bars are not visible, errors are smaller than the symbol width. doi:10.1371/journal.pone.0040042.g005 DG ATP~2 0kh. The maximum efficiency was more than 30%, which is comparable with experimental values [38]. We observe that the model prediction is rather good even though it somewhat underestimates the efficiency. This is likely because we only consider two states to simplify the model. In other words, we ignore the condition where myosin heads attach and detach to the actin filament without generating force or consuming energy (weakly attached state). Incorporating this state should increase the efficiency. To estimate how much of an increase, we computed the physical upper bound of the efficiency produced by the ratchet according to the formula proposed in [39]: E t~vi E a zE e is the total energy computed before (t { i ) and after (t z i ) the i-th jump, j i . The summation is done for every jump that occurs per unit time in the region of steady shortening. The efficiency limit is thus given by g eff~{ f ext v=r s shown in Fig. 6 (squares). The experimental efficiency is between the two model efficiencies, indicating a more detailed description of the states can improve the fitting.

Discussion
In this work, we have used diffusion theory to model muscle mechanics, defining some parameters using SME data and quantitatively simulating several macroscopic behaviours. The model satisfies two different hypotheses regarding the origin of the actomyosin stable states, which are often viewed as stable states of the lever arm (scenario 1). This interpretation gives a clear picture of the role of the lever arm structure in muscle contraction and is compatible with X-ray reflections during isotonic contractions that describe the regular arrangement of myosin motors on the backbone [7,8]. However, if we assume scenario 2, where each minimum of the actomyosin potential energy corresponds to a given actin monomer on the thin filament, the role of the lever arm structure is more equivocal. That the actin diameter is about 5:5 nm means the vertical distance between the backbone and different attachment sites in one half-pitch varies considerably due to the three-dimensional structure of the actin filament. Thus, we can assume that the orientation of the lever arm plays an important role for myosin reaching the binding sites. In our model, we assumed one-dimensional geometry for simplification. Nevertheless, we can still infer some aspects of the lever arm orientation at different minima. In the insets of fig. S4, we describe the distribution of the attached heads in the four subsequent minima in one half-pitch shown in Fig. 1 from the highest (d) to lowest (4d) minimum at different times of an unloaded force-clamp simulation. We can see that at time t 1 during an isometric contraction, all four minima are populated, but the central ones are more populated. As expected, at the end of phase 2 in the load-clamp simulation (time t 2 ), the lowest minimum (4d) is the most populated, while the first two minima populations are near zero. By associating each minimum a different lever arm orientation, the periodic distribution of the myosin head mass along the thin filament becomes very different such that the distribution can justify the different experimental X-ray reflection observations produced at times t 1 and t 2 [7,8]. Also, during isokinetic contraction (time t 3 ), the population re-equilibrates toward a more homogeneous distribution, which is consistent with experimental observations.
It is worth noting that scenario 1 allows for more flexibility than scenario 2 when defining parameters such as the distance between minima. Nevertheless, the good fitting obtained shows that either scenario is consistent with the experimental data. Moreover, high speed AFM showed that a myosin V head undergoes brief translocation along actin by repeatedly detaching and attaching after Pi release, i.e., at the force generating state [40], giving creedance to scenario 2. To determine which of the two scenarios is the true source of force generation, more experimental evidence is needed.
One limit of our model is that the contraction velocity is slightly underestimated for mean values of applied force. Consequently there is an enhanced underestimation of the maximum power output (T|v)=(T 0 |v max ) (Fig. S5). This limit is common to other macroscopic models [20,35], and can be due to a low number of simulated attached heads in the central part of the force-velocity curve. The number of attached cross bridges can be increased by a higher attachment rate, but that would create an unrealistically fast rising phase during isometric contraction [35]. A way to resolve the problem is to hypothesize an attachment rate that depends on the velocity of the contraction in addition to the dependency on the position of the head [20]. Despite using a different physical description of force generation in the attached state (diffusion), our model shares with the above models the same chemical description of the attachment-detachment process (rate functions) because of a lack of experimental data. Thus, it is not surprising that our model suffers from underestimating the maximum power output. Resolving this problem would benefit from a description of the attachment-detachment process by incorporating diffusion into the energy landscape.

Conclusions
Although the mechanics of muscle contraction have been well studied, a comprehensive model that reproduces at the same time all macroscopic behaviour is still lacking. Furthermore, one weakness in muscle modelling is that the same behaviour can be reproduced starting from different working models. We believe that this is in part because the description of the oscillatory behaviour of the actomyosin complex, which is generally accepted Figure 6. Predicted efficiency at different velocities Simulation of the efficiency of contraction, as predicted by the model compared to experimental data (circles from [38], green line). Efficiency is computed as tension times velocity divided by the chemical energy consumed [37] (triangles) and following the method described in [39] (squares). Mean values+standard deviation over 11 trials are shown. Where error bars are not visible, errors are smaller than the symbol width. doi:10.1371/journal.pone.0040042.g006 as a fundamental property of muscle function [1] and whose properties are becoming more accessible from SME, is often restricted to rate functions based on the macroscopic behaviour of muscle fibers.
To emphasize the importance of this oscillatory behaviour, we have applied a new approach to muscle modelling, where actomyosin complex dynamics are described by the diffusion of particles constrained within a defined energetic landscape that is described by SME data. Having inferred the bias of actomyosin energy, we showed the feasibility of this approach for quantitatively describing several macroscopic behaviours. Furthermore, the new approach reduces the number of free parameters and makes them constant rather than empirical functions that are based on chemical reaction coordinates. This model can potentially be a basis for future developments where parameters that describe the behaviour of single molecules can be built upon to understand all contraction mechanism details, such as how molecular cooperation leads to macroscopic behavior.

Supporting Information
Figure S1 SME set-up from [32]. An external system (microneedle) with a very high drag coefficient is linked to a multistable particle (attached skeletal myosin). The motion of the particle is analysed through the position of the microneedle. A linear spring with symmetric stiffness links the two bodies. (EPS) Figure S2 Example of a simulated microneedle trace based on Figure S1. The myosin head is always in the attached configuration. The trace is comparable with the experimental results obtained in [15]. (EPS) Figure S3 Simulation of the isometric contraction. The simulated rising phase is slowed by a factor of v s r =v e r~4 0 (blue trace). Simulations for v s r =v e r~3 0 (green trace) and v s r =v e r~5 0 (yellow trace) are also shown. The slowed rising phase is in agreement with the experimental data (red trace from [42], the initial difference is due to the latency of force development and stretching of the tendons, which are not considered in the model). (EPS) Figure S4 Distribution of the attached myosin heads (insets) in the four minima of actomyosin energy E a during different phases of the unloaded force-clamp simulation (continuous line). Scale in the insets is normalized to the number of attached heads. (EPS) Figure S5 Power output predicted by the model (blue triangles) and comparison with experimental data (black circles). The lower velocity of contraction at intermediate tensions leads to an underestimation of the maximum power output. This drawback is common with other models and is related to the attachment-detachment process (see main text). Mean values+standard deviation over 11 trials are shown. Where error bars are not visible, errors are smaller than the symbol widths.