Including Thermal Fluctuations in Actomyosin Stable States Increases the Predicted Force per Motor and Macroscopic Efficiency in Muscle Modelling

Muscle contractions are generated by cyclical interactions of myosin heads with actin filaments to form the actomyosin complex. To simulate actomyosin complex stable states, mathematical models usually define an energy landscape with a corresponding number of wells. The jumps between these wells are defined through rate constants. Almost all previous models assign these wells an infinite sharpness by imposing a relatively simple expression for the detailed balance, i.e., the ratio of the rate constants depends exponentially on the sole myosin elastic energy. Physically, this assumption corresponds to neglecting thermal fluctuations in the actomyosin complex stable states. By comparing three mathematical models, we examine the extent to which this hypothesis affects muscle model predictions at the single cross-bridge, single fiber, and organ levels in a ceteris paribus analysis. We show that including fluctuations in stable states allows the lever arm of the myosin to easily and dynamically explore all possible minima in the energy landscape, generating several backward and forward jumps between states during the lifetime of the actomyosin complex, whereas the infinitely sharp minima case is characterized by fewer jumps between states. Moreover, the analysis predicts that thermal fluctuations enable a more efficient contraction mechanism, in which a higher force is sustained by fewer attached cross-bridges.


S1.1 Model and numerical simulations
The myosin motor can exist in three states: detached (D), weakly attached (W), or strongly attached (S). The D and W states switch at constant transition rates k DW and k W D , respectively; the W-S (k W S ) and S-D (k SD ) transitions are functions of the myosin head position, as implemented in the Huxley 1957 model [1].
x < 0 α a x, 0 < x < l t 0, x > 0 The S state is described as a continuous energy landscape with three minima, corresponding to a pre-power stroke (S0) and a power stroke with two stages (S1 and S2). In SL, the central region of the biochemical energy landscape is defined as in Equation 5. ATP energy release favors the power stroke by biasing the biochemical energy such that E S0 c > E S1 c > E S2 c . To simulate this bias, we add a linear drop F atp to the flat sinusoidal part of E c . α d is a constant angle that adjusts for the constant ATP drop by shifting the first minimum to zero. The convex part is expressed by a polynomial, ensuring a continuous first derivative. The energy barrier in the sinusoidal function is chosen as H = 6κ b T 0 (in the two flat landscapes, the maximum is less important because it is included in the unloaded rate constants; see main text). We can then approximate the minimum in E c as a parabola with stiffness H(2π/d) 2 , which corresponds to k c = 48 pN/nm. Therefore, the energy bias due to the ATP is F atp = 8κ b T 0 /d (where T 0 = 37 • C is the reference temperature and d is the distance between two consecutive minima). Note that the total energy drop due to the ATP in the biochemical component is 16κ b T 0 , below the total amount of energy released by ATP hydrolysis (approximately 18κ b T 0 ). The ATP energy difference is here attributed to Huxley's 1957 [1] hypothesis, whereby preferential attachment at higher stretching levels requires energy to break the detailed balance. According to H57, the pre-power stroke state also generates a force, here justified by the above-mentioned ATP energy difference.

S1.2 Number of myosin molecules, TT units, Cooperativity and Geometry
The lengths of the half-actin and myosin filaments are 1224 nm and 925 nm, respectively. The myosin filament includes 50 nm of bare zone with no myosin heads, implying 60 myosin molecules per thick filament. Each thick filament is surrounded by six double-helix thin filaments, each of which can be reached by the myosin heads on three thick filaments; therefore, each myosin filament contains approximately 30 double-headed myosin molecules, which can interact with each actin filament. Because troponin-tropomyosin (TT) units are related to each single strand in the helix, this number is strictly halved, but is restored to 30 to accommodate two myosin heads per molecule. To allow for possible interaction from other myosin heads [3], we increase the number of myosin heads interacting with each actin filament to 38, as in [2]. The total force is then multiplied by two to recover the double-stranded helix.
Muscle contraction and release is governed by the amount of Ca in the troponin binding sites, which is able to shift the tropomyosin filament exposing F-actin to myosin motors. The occupancy of these binding sites is governed by the myoplasmic free [Ca 2+ ], which diffuses, following an action potential on the cell membrane, from the sarcoplasmatic reticulum (SR) release channels to the myofibrillar space, and is subsequently taken up by the SR Ca pumps. As a consequence, the myoplasmic free [Ca 2+ ] depends on the particular experimental or physiological situation. The actin filament houses several TT units, each with two Ca affinity sites, spaced at 36nm intervals. The rate of Ca ion occupancy depends on the instantaneous Ca concentration. When both Ca sites are empty (Ca-off in S1 Fig) or only one site is filled (Ca-on*), the TT unit is inactive and the corresponding portion of the actin filament has lower affinity for the myosin head (modeled through the variable Q Ca ). The affinity increases when both sites are filled (Ca-on). Moreover, if at least one myosin is already attached to an actin monomer corresponding to a TT unit, the deformation of the tropomyosin filament reduces the detachment rate of the captured Ca ions (by factors of 4 and 2 in Ca-on → Ca-on* and Ca-on* → Ca-off transitions, respectively).
In single-fiber simulations, we impose a tetanic state such that the myoplasmic free [Ca 2+ ] is constant and sufficiently high to saturate the Ca binding sites in the troponin. In the whole-heart simulations, the free [Ca 2+ ] transient is imposed as in S2 Fig In the whole-heart simulations, the free Ca at each side of the ventricle is computed by conducting an excitation propagation simulation [4], whereby the human cardiomyocyte model [5] is applied for the computation of the ion channel currents. S2  In the single-fiber simulations, the calcium concentration is fixed at its saturated value (tetanic situation), whereas in the whole-heart simulation, the Ca concentration follows a realistic time-transient trend. Moreover, the transient is shifted in each sarcomere, approximating the real spreading of an electrical signal.
The cooperativity of myosin attachment affects the transition rates k DW and k W D through a parameter γ n , where n is 0 if there are no nearest neighbors and 1 or 2 if one or two myosin heads are already weakly or strongly attached, respectively. Finally, to account for the overlapping of filaments in over-compressed or over-stretched situations, we multiply k W S by a sigmoidal factor χ: where x s is the overlapping or overstretched distance of the given myosin head and a s is a constant. In the absence of overlapping or overstretching, x s < 0. These three effects are modeled as described in [2], which gives: At each time step, the algorithm updates the position x i of the i − th myosin head with respect to its anchor on the thick filament. In this way, the total tension generated is given by: where N XB is the total number of myosin heads, N f il is the number of filaments in the sarcomere, and k(x i ) explicitly expresses the dependence of the asymmetric stiffness on the myosin position. The factor of two accounts for the second strand, as mentioned above. In SRI and SRII, the position of the myosin head is first updated based on the position of the actin filament, z, in the previous time step: Here, z a i and x a i are the positions of the actin filament and the myosin pre-stretch, respectively, at the time of attachment of the i − th myosin head. n ps equals 0, 1, or 2 when the myosin head resides in the S0, S1, or S2 minimum, respectively. Having updated the position of each myosin head, the algorithm updates the position of the actin filament using Newton-Raphson method. In each iteration, the total force is computed by Equation S1.
In the implicit scheme, an equation of the form is approximately solved so that the following two equations are simultaneously satisfied at every time step from t to t + ∆t.
where q is the internal force vector given as a function ofu and u, and f corresponds to the external force vector. In the single-fiber simulations, u is an (N x b * N f il + 1) vector of the positions of the myosin motors and actin filaments, whereasu represents their velocities. q is the force generated by myosin elastic elements and their drag force, as well as the drag force on the actin filament. f is the external force applied to the actin filament (when present). In SL, this also contains the random forces acting on each myosin motor. Finally,γ is a parameter of between 0 and 1. In our case,γ = 0.6. In the heart simulator, we use a similar but more complicated procedure, because we apply a number of Monte Carlo time steps (∆t = 1 µs) within each Finite Element Model time step (∆t = 1.25 ms). Refer to [2] for a detailed explanation. In Newton-Raphson method, Equations S4 and S5 are iteratively solved by linearizing the internal force term, as described in the following procedure. First, the initial guess is given by: Then, the first iteration is conducted as follows: Residual computation: PLOS Solution for the velocity update: Updating: t+∆t u (1) = t+∆t u (0) +γ∆t∆u (1) + ∆t tu (S10) t+∆tu(1) = t+∆tu(0) + ∆u (1) (S11) Here, C (0) = dq du ( t+∆tu(0) , t+∆t u (0) ), and K (0) = dq du ( t+∆tu(0) , t+∆t u (0) ).
From the second iteration (k > 0), the following procedure is conducted until the residual norm becomes sufficiently small.
Residual computation: Solution for the velocity update: (C(k) +γ∆tK(k))∆u (k+1) = r (k+1) (S13) Updating: t+∆t u (k+1) = t+∆t u (k) +γ∆u (k+1) (S14) t+∆tu(k+1) = t+∆tu(k) + ∆u (k+1) (S15) Here, C (k) = dq du ( t+∆tu(k) , t+∆t u (k) ), and K (k) = dq du ( t+∆tu(k) , t+∆t u (k) ). The difference between the first and subsequent iterations in the update of u is necessary to satisfy Equation S5. This generates the difference in the right-hand side of the solution for the velocity update. In SL, Equation S2 is not used; rather, the actomyosin complex is treated as a diffusing, over-damped particle with drag coefficient η. The complex is subjected to thermal, elastic, and chemical energies, the latter describing the actomyosin interaction and the ATP effect. In the equation in S1 Fig, the biochemical and elastic (mechanical) components are separated as E c and E e , respectively, and the third term describes the temperature effects. ω = 0 when the myosin head is detached and 1 when it is attached. Γ(t) is a random term satisfying < Γ(t) >= 0 and < Γ(t 1 ), Γ(t 2 ) >= 2δ(t 1 − t 2 ) (white noise). The total tension is computed as: where A f is the cross-sectional area of one filament. The work done by each ventricle is computed as the ejection flux multiplied by the pressure at the valves minus the work done by injection. The efficiency of the whole heart is computed as the ratio of total work to total ATP consumption.
At each time step ∆t and for each possible transition, a random number, uniformly distributed between zero and one, is generated. The transition occurs if this value is less than or equal to ∆tk, where k is the defined rate constant. ∆t = 1 µs in SRI, SRII, and the KS approximation of SL, and ∆t = 1 ns in SL. The value of k is modified accordingly. At the single cross-bridge and single sarcomere/fiber levels, the numerical simulations of scenarios SRI and SRII and scenario SL differ as follows.
The algorithm for SRI and SRII is:

PLOS
3. Define the stable state of each myosin head (D, W, S0, S1, or S2) using the rate constants k DW , k W D , k W S , k SD or the rate constants produced in step 1 (k 12 , k 21 , k 01 , k 10 ), depending on the stretching of the elastic component of the myosin head.
4. Update the positions of each myosin head and actin filament, depending on the setup and external conditions, and calculate the total force.

Increment the time step and return to
Step 2 until the total time is equal to the prescribed simulation time.
The algorithm for SL is slightly different: 1. Update the states of the TT units depending on the Ca concentration.
2. Define the stable state of each myosin head (D, W, S) using the rate constants k DW , k W D , k W S , and k SD .
3. Using an implicit method, determine the position of each myosin head from Equation 6.
4. Update the external conditions depending on the setup and calculate the total force in the current time step.
5. Increment the time step and return to Step 1. Repeat until the total time is equal to the prescribed simulation time.
To account for the different temperatures in the experimental observations of single fibers and cardiac contraction, we include a factor Q 10 in the definitions of k DW , k W D , k W S , and k SD . We assume that the biochemical and mechanical energy components are independent of temperature. The temperature effects on the oscillations are introduced through κ b T in the rate constants of SRI and SRII and in the Langevin equation for SL. The temperature effect on the TT units is neglected. The Langevin approach is unsuitable for multiscale analysis because of its intrinsic time step limitations. Instead, we use the KS approximation of the SL. The algorithm is the same as that used for SRII, but the rate constant matrix is generated by the KS approximation of the Langevin equation (see Equation 6). To account for the oscillations within the minima (and in the detached and weakly attached states), the position of each myosin head at each time step is computed by randomly selecting a position from the corresponding stationary probability distribution in the current state: where N is the normalization constant. This action modifies the last term in Equation S2 to account for the shifting of the total energy minima under the wide minimum hypothesis.

S1.3 Comparison between the Langevin and Kramers-Smoluchowski (KS) approximations
The numerical approximation to the system of Langevin equations describing the motion of each myosin molecule has several advantages, as it reproduces the detailed myosin dynamics, including the real effect of the energy barrier shape on the dwell time of the protein and other phenomena that cannot be simulated by different approaches. However, its applicability is limited by the long computation time, because the time step for the analysis must be smaller than the typical timescale of the process. In the PLOS myosin analysis, the time step is a few nanoseconds, which is feasible for simulations of single sarcomeres, but not for whole-heart simulations. For this reason, we adopt a rate constant approach based on the KS approximation of the energy landscape used in SL at the single sarcomere level. Here, we show that the approximation is sufficient for the purpose of our study by comparing the results of the KS approximation with the force clamp and the length step protocols in both SL and SRII. As shown in S3 Fig and S4  Fig, the

S2 Tension recovery after a small step
When a small, fast length change is imposed on an isometrically contracting muscle, the tension changes almost instantaneously from T 0 to T 1 , then exponentially recovers toward T 0 with a rate r, settling at a new value T 2 . The initial change is related to the mechanical relaxation of the elastic element, and the recovery is related to the equilibrium change caused by the power stroke [6]. r, T 1 , and T 2 depend on the imposed step δ. As experimentally observed, the T 2 value is somehow hidden by the detachment of old myosin heads and the attachment of new motors. To overcome this problem, we consider the population of cross-bridges in S2 to reach a maximum when the old cross-bridges have completed the power stroke. For this reason, we define the values of the variable T 2 (δ) and r(δ) at the time of the maximum in S2 (see S5 Fig).
The simulated data for T 1 and T 2 (see S6 Fig) are very similar in the three scenarios, and are comparable to the experimental behavior.