The effect of molecule tethering in three-dimensional (3-D) space on bimolecular binding kinetics is rarely addressed and only occasionally incorporated into models of cell motility. The simplest system that can quantitatively determine this effect is the 3-D sarcomere lattice of the striated muscle, where tethered myosin in thick filaments can only bind to a relatively small number of available sites on the actin filament, positioned within a limited range of thermal movement of the myosin head. Here we implement spatially explicit actomyosin interactions into the multiscale Monte Carlo platform MUSICO, specifically defining how geometrical constraints on tethered myosins can modulate state transition rates in the actomyosin cycle. The simulations provide the distribution of myosin bound to sites on actin, ensure conservation of the number of interacting myosins and actin monomers, and most importantly, the departure in behavior of tethered myosin molecules from unconstrained myosin interactions with actin. In addition, MUSICO determines the number of cross-bridges in each actomyosin cycle state, the force and number of attached cross-bridges per myosin filament, the range of cross-bridge forces and accounts for energy consumption. At the macroscopic scale, MUSICO simulations show large differences in predicted force-velocity curves and in the response during early force recovery phase after a step change in length comparing to the two simplest mass action kinetic models. The origin of these differences is rooted in the different fluxes of myosin binding and corresponding instantaneous cross-bridge distributions and quantitatively reflects a major flaw of the mathematical description in all mass action kinetic models. Consequently, this new approach shows that accurate recapitulation of experimental data requires significantly different binding rates, number of actomyosin states, and cross-bridge elasticity than typically used in mass action kinetic models to correctly describe the biochemical reactions of tethered molecules and their interaction energetics.

## Introduction

The effect of concentration ratio of tethered molecules and their ligands, constrained by their spatial arrangements and occupancy, is rarely addressed and only sporadically incorporated in models of cell contraction and adhesion. Almost all models include concepts originally introduced by Kramers (Kramers, 1940; Evans and Ritchie, 1997) where the strain dependence of a model state transition is taken in only one dimension. The spatial occupancy of interacting species is only partially considered in mass action kinetic models or by invoking simplified spatial rules that do not reflect accurately the local spatial arrangement and occupancy of the species. The best defined system to quantitatively determine this effect is in striated muscle because of its well-ordered hexagonal lattice of interacting actin and myosin filaments. Because the myosin molecules are tethered to myosin-containing thick filaments and actin monomers are arranged in filaments with the myosin-binding sites associated with each monomer, the flux of the myosin binding is not directly associated with the concentrations of myosin and actin species, but is instead strongly affected by the number of myosin-binding sites on actin filaments that can be reached by the tethered myosin molecules and whether or not any of these sites is already occupied by bound myosin. This behavior can be quantitatively addressed by a new generation of multiscale models that bridge scales from atoms to cells and to muscle tissues and provide detailed information about how structure relates to function.

Attempts to computationally model muscle contraction started soon after the discovery of sliding filaments (Huxley and Niedergerke, 1954; Huxley and Hanson, 1954), but most of these attempts are so-called mass action kinetic models of muscle contraction (e.g., Huxley [1957]) which do not take into account the discrete geometrical relationship between individual myosin heads and sites on actin in the 3-D sarcomere lattice. Instead, most of the models used state probability density distribution functions along the reaction coordinate, usually denoted as bond strain (e.g., Hill [1974] and Smith and Geeves [1995a]). In the first and best known mass action (two-state) model, the strain is expressed as the relative displacement, *x*, between the positions of a myosin-binding site on an actin filament (for simplicity a site on actin) and of unstrained myosin head (Huxley, 1957). This strain can be related via bond stiffness to bond force and stored elastic energy. In this and later mass action models, strain-dependent transitions between the states are governed by rate transition functions, i.e., prescribed functions of the stretch of the tethered myosin molecule to reach a binding site or being bound to it. The derivation of these strain-dependent rate functions from the reaction energy landscapes was elegantly formalized by Hill (1974) and later used in most subsequent mass action models.

The main problem for all previous mass action models is in the calculation of the transition flux between the unbound and bound myosin states. In the original Huxley 1957 model (Huxley, 1957) and in most of the later mass action models, it is assumed that the sum of state probability density distribution functions at each strain, *x*, is equal to one. Consequently, the state transition fluxes at the coordinate *x* are proportional to the sum of products of state transition rates at *x* and the corresponding state probability density distribution function at the same location (Huxley, 1957; Hill, 1974). Thus, if a subset of the bound myosins moves out of the binding region, then the sum of the state probability density distribution functions within the binding region, necessary for accounting of the net flux between states, cannot add up to one unless an equivalent number of myosin heads are drawn into the binding region. However, strictly speaking, this approach is ill posed because the pool of unbound myosin heads is finite and also modulated by the number of already attached heads. To get around this problem, most authors, including A.F. Huxley, focused exclusively on the state probability density distribution functions of myosin bound to actin that have finite state probabilities. This approach accounts only for the fraction of bound heads to the total number of heads, i.e., the probability density of myosin head to be in a bound state. Another difficulty with these approaches is to ensure that the total number of myosin heads is preserved, i.e., that the sum of the probabilities of myosin being in any state of actomyosin cycle is equal to one. During shortening, lengthening, or sudden change of length, the number of myosin-binding sites on actin filaments coming into the range of myosin binding typically does not match the number of detached myosins; thus, the sum of bound and detached myosins is not preserved. Piazzesi and Lombardi (1995) were the first to attempt to resolve this problem by proposing a periodicity rule to reduce the fraction of available sites on actin for the fraction of bound myosin outside the myosin-binding region. Using a similar approach, Mijailovich et al. (2000) correctly accounted for conservation of myosin species in a four-state model of smooth muscle. These approaches enable the conservation of the myosin species and allow for a more precise accounting of the fluxes between the states. However, they may not reflect correctly the populations of bound myosin to sites on actin, of free myosin-binding sites on actin, and of the unbound myosin species in the sarcomere lattice. Also, even with these approaches, mass action models do not account for the number of available sites on actin filaments that are several times larger than the number of active myosin heads.

In the hexagonally packed, myofilament array in the sarcomere, myosins are tethered to the thick filament backbone and can only bind to a few, geometrically accessible, actin-binding sites, positioned within the limited range of movement of myosin heads. These geometrical constraints strongly distort the classical view of biochemical reactions. There have been several models proposed to prescribe the rules determining the interactions between pairs of myosin molecules and multiple myosin-binding sites on actin filament, so called target zones. Daniel et al. (1998) formulated a spatially explicit model of interaction for only two filaments, defining that binding is only possible at collinear sites. This approach was extended by Chase et al. (2004) to a 3-D lattice of three myosin filaments interacting with 11 actin filaments and with binding restricted to collinear sites. Taking into account only collinear sites provides for an unrealistically low number of bound cross-bridges. To overcome this problem, we proposed a more comprehensive approach (Smith et al., 2008) in a spatially explicit model of muscle contraction that defines the rules of actin–myosin interactions in the filament lattice of the half sarcomere. In this approach, we introduced longitudinal and azimuthal (Steffen et al., 2001; Smith et al., 2008) selection rules in a 3-D sarcomere lattice for mapping heads to target zones of adjacent sites on actin. This methodology yielded a more realistic number of bound cross-bridges and magnitude of force per myosin filament (Smith and Mijailovich, 2008), but the overall approach was limited by the use of semi-probabilistic binding kinetics that did not account for the stochastic nature of myosin binding to actin and, therefore, could not precisely account for the numbers of myosin and actin species entering the reaction.

To overcome the limitations of these earlier models, a better method to account for realistic matching of available myosin heads to actin-binding sites is needed. To this end, we have developed the computational platform MUSICO (Muscle Simulation Code) as a tool that can be used by the muscle community to interpret experimental data, identify gaps in our current understanding, and provide information that can be used in the development of novel hypotheses and the design of critical new experiments to test them. The structure of MUSICO allows quantitative assessment of the effect of variations in sarcomere geometry, incorporation of new models with kinetic and structural details of the actomyosin cycle, provide more detailed descriptions of thin filament regulation and provide new knowledge on the role of auxiliary muscle proteins in muscle functions. The platform explicitly incorporates the 3-D sarcomere structure with extensible actin and myosin filaments, along with various actomyosin cycles, thin filament regulation models, and protocols of biochemical and mechanical loading (Mijailovich, S.M., et al. 2007. Biophysical Society 51st Annual Meeting. Abstr. 155a; Mijailovich, S.M., et al. 2008. Biophysical Society 52nd Annual Meeting. Abstr. 404a; Mijailovich, S.M., et al. 2009. Biophysical Society 53rd Annual Meeting. Abstr. 201a). Fig. 1 B shows the modular structure of the MUSICO platform, which permits various combinations of modules. The MUSICO simulations can predict the spatial positions and connectivity of sarcomeric proteins at the current sarcomere configuration, local and global mechanical forces, and energetics at any instant of time, taking into account the exact number of unbound myosin molecules interacting only with unoccupied binding sites on actin filaments in their proximity. The modular structure of the platform permits incorporation of the newest findings in a particular module without need for any change of other components and, therefore, provides a quick way to implement and test the implications of new hypotheses or findings. These may include, for example, the packing of myosin heads onto the thick filament backbone or the role of MyBP-C.

In this first paper, we focus on one very important aspect of this modeling approach: we quantitatively assess how the concentration ratio of tethered molecules and their ligands affect predictions of classical muscle contraction experiments. We compare the predictions of a spatially explicit stochastic model of muscle contraction (MUSICO) that takes into account the interaction of each myosin head with a few nearest binding sites on actin in the 3-D sarcomere lattice and those of mass action models that take into account only axial strain dependence and neglect the discrete positions of myosin heads and the sites on actin. To illustrate the differences between these models, we compared the two simplest models proposed for cross-bridge kinetics: (1) the two-state A.F. Huxley 1957 model (Huxley, 1957) and (2) a three-state model including a power stroke (Duke, 1999). The same state transition rates are used for both the mass action models and in the new stochastic formulation in the 3-D sarcomere lattice. This comparison allows determination of the effects of the distortion of interacting molecules on biochemical reactions, i.e., modulation of actomyosin cycle state transition rates by the imposed geometrical constraints of tethered myosin molecules and the accessibility of binding sites on actin filament.

This is the first step toward using MUSICO to build a comprehensive multiscale model of a contracting multi-sarcomere structure or a myofibril. However, for simplicity and without loss of generality, here we limit our current simulations to contractions of a half sarcomere. Even with this restriction, this typically corresponds to 500 myosin and 1,000 actin filaments, which is approximately the number of filaments in a cross section of a myofibril, and creates realistic fluctuations in force that we expect to find in muscle fibers. We demonstrate that this 3-D lattice model is essential to predict the correct number of attached cross-bridges and hence the force per cross-bridge, both necessary first steps before going on to develop the next stage of the model.

## Materials and methods

### Actin–myosin interactions in the filament lattice

Actin–myosin interactions in a sarcomere lattice incorporate the discrete lattice structure of interdigitated actin and myosin filaments, and the lattice structure varies in different kinds of muscles (Luther and Squire, 1980). In most muscle types, myosin and actin filaments are packed into a regular 3-D lattice with three myosin filaments around each actin filament and six actin filaments around a myosin filament (Luther and Squire, 1980). Each myosin filament is decorated with crowns of myosin dimers, spaced by ∼14.3 nm along the filament (Fig. 2). In vertebrate muscle, each crown consists of three myosin dimers with transverse orientations spaced by 120°, and successive crowns along the filament are rotated by 40° looking toward the Z-line. Thus, crown orientations are repeated every 42.9 nm; i.e., every fourth crown has the same orientation as the initial.

The actin monomers in a thin filament form a double-stranded helical structure associated with the regulatory proteins tropomyosin and troponin with binding sites apart by ∼5.5 nm on each strand, with a half-period of ∼35.75 nm in the relaxed state (Bordas et al., 1999). The difference in periodicities between actin-binding sites (∼35.75 nm) and myosin crowns (42.9 nm) creates a vernier of longitudinal spacing between myosin heads and actin-binding sites. The 3-D sarcomere geometry with the extensible filaments and myosin-binding sites on actin filaments require both longitudinal position matching (Fig. 2 B) and angular matching in the azimuthal plane (Fig. 2, C and D). A myosin head and the closest binding site on actin form the most probable pair of these molecules that can create a cross-bridge interconnecting actin and myosin filaments. The concept of matching and realignment of myosin heads and its binding sites on actin in a compliant filament sarcomere lattice was first introduced in a Monte Carlo model by Daniel et al. (1998) and represents the basis for the more comprehensive approach introduced below.

#### A spatially discrete 3-D model

The model of 3-D sarcomere lattice (Fig. 2 A) is composed of thick and thin filaments arranged in an overlapping hexagonal lattice with a geometry that is consistent with the mean spacing measured in vertebrate striated muscle (Higuchi et al., 1995; Millman, 1998). The minimal contractile unit can be considered as a myofibril. For a typical myofibril of 1.2 µm in diameter, the number of myosin filaments is ∼500–650, depending on interfilament distance *d*_{10} (Millman, 1998). A myofibril runs the entire length of the muscle fiber and consists of a large number of sarcomeres in series. In a half sarcomere, each myosin filament with ∼150 myosin molecules faces six actin filaments, where each actin filament has ∼360 myosin-binding sites (∼260 active at maximum overlap) arranged in a double helix (Fig. 2 B). Each myosin molecule has two myosin heads. The kinetics of the two heads of the dimer can be modeled by assuming that they compete for sites from the same target zone (Smith and Mijailovich, 2008; Smith et al., 2008), but this process is complex and beyond the scope of this paper. For simplicity, we consider here that only one myosin head is active at any one time and is denoted as a cross-bridge.

Even with this simplification, this model fully represents the contractile behavior of the 3-D filament lattice unlike previously studied one-dimensional models (Huxley, 1957; Wood and Mann, 1981; Pate and Cooke, 1989; Piazzesi and Lombardi, 1995; Smith and Geeves, 1995a,b) where the discrete position of myosin heads and their binding sites on actin were ignored. The 3-D sarcomere structure is viewed as an array of thin and thick filaments connected by cross-bridges and other elastic elements (e.g., titin and MyBP-C) in a lattice network and that all of these elements are represented as linear springs (Daniel et al., 1998) or beams if radial forces are taken into account.

For balancing only axial forces, the spring system is sufficient for obtaining the displacement field in any instant of time, assuming that the movement of the filaments is restricted in both radial and azimuthal directions. However, radial constraints imposed on the lattice of intact muscle fibers through shrinking, osmotic swelling, or by changes in muscle length may induce significant radial forces on cross-bridges (Millman, 1998; Williams et al., 2012, 2013). In this case, the 3-D sarcomere structure should be represented by interconnected beams in order to form a mechanically stable system. The resulting radial forces may also affect the strain dependence of cross-bridge kinetics and, therefore, overall muscle dynamics. In either case, the instantaneous equilibrium between actin and myosin filaments interconnected by cross-bridges is described as

where **K**(*t*) is a stiffness matrix that includes the elasticity of thick and thin filaments and attached cross-bridges, **u**(*t*) is the displacement vector of all myosin-binding sites along actin filaments and all myosin links to the backbone of the myosin filament, and **F**(*t*) is a vector of all external forces (load) and internal forces generated by the action of cross-bridges (e.g., power stroke). The number of these linear equations and, therefore, the unknown displacements is equal to the total number of all actin and myosin sites. If the actin and myosin filaments are assumed to be homogeneous along the filament length, the number of linear equations can be significantly reduced by fusing all finite elements between two attached cross-bridges into a single element, i.e., taking in account only interconnected nodes on actin and myosin filaments by the cross-bridges. The stiffness matrix, **K**(*t*), is constantly changing as actomyosin bonds are created or disrupted. In addition, **F**(*t*) also changes during conformational changes in attached myosins. Thus, the matrices **K**(*t*) and **F**(*t*) must be updated after any chemo-mechanical transition occurs in order to accommodate structural system changes caused by these transitions. The formulation and solution of these equations is obtained by standard finite element procedures for linear systems (Eq. 1) or an incremental procedure for nonlinear systems (Bathe and Mijailovich, 1988; Mijailovich et al., 1993; Bathe, 1996; Kojic, 1996; Kojić et al., 1996; Mijailovich, S.M., et al. 1998. Biophysical Society 42nd Annual Meeting. Abstr. A156).

#### Calculation of cross-bridge strain and the strain dependent kinetics of actomyosin cycle

The elasticity of myosin molecules and thermal agitation define how far a given myosin head can move to reach actin-binding sites. The 3-D nature of myosin binding in the sarcomere lattice requires that myosin heads move not only axially but also azimuthally to reach actin monomers with the correct orientation (Squire, 1992). The discrete nature of the myosin-binding sites in two helically arranged strands in each actin filament and the regular arrangement of myosin crowns along the myosin filaments, each crown having three myosin molecules, form a large number of possible combinations of relative spatial positions between myosin heads and binding sites on the actin filament. These relative positions can be calculated from the coordinates of deformed actin and myosin filaments, where the deformation of the elastic filaments is obtained from the solution of Eq. 1 for appropriate boundary conditions.

In the 3-D sarcomere lattice, the relative distance between a myosin head and the adjacent binding site on actin is defined by four factors (Figs. 3 and 4): (1) the axial displacements along myosin and actin filaments, (2) the transverse distance between myosin and actin filaments, $d M-A =2 3 d 10 , $ (3) the angle α defining relative position of myosin to actin filament, and (4) the angle β defining how much a myosin head needs to turn to reach an actin monomer in the correct orientation (Figs. 3 and 4). These spatial distances and angles are essential information for formulation of the strain dependence of the state transition rates, reflecting the discrete geometric relationships between myosin heads and actin-binding sites.

In all sliding filament models hitherto, the strain dependence of the actomyosin cycle has been prescribed only in the axial direction, except in Williams et al. (2012) where radial forces are also taken into account. We follow this approach for convenience of comparing sarcomeric contractions predicted by MUSICO with mass action models, but in addition, we consider three other geometric factors, expressed as weight functions associated with myosin binding to actin in 3-D sarcomere lattice.

Let us start with the myosin binding step because myosin heads can interact with only one myosin-binding site on actin out of several reachable adjacent sites (target zone). From the solution of Eq. 1, at any instant of time, the displacements **u**(*t*) and positions of each myosin molecule **X*** _{M}*(

*t*) and binding site on the actin filament

**X**

*(*

_{A}*t*) in the deformed sarcomere configuration are known. In order for a myosin head to reach the adjacent binding sites on actin filament, the elastic component of a cross-bridge needs to be strained depending on relative position of a myosin head (cross-bridge) $XMm(t)$ and adjacent binding sites on actin $XAn+1(t),$ where superscripts

*m*and

*n*+ 1 denote a myosin molecule and the adjacent sites on actin, where

*l*= 0, 1, 2, …, ($\mathcal{L}a$ − 1) is an index of accessible sites on actin in neighborhood of

*n*, respectively. The strain vectors for $\mathcal{L}a$ adjacent binding sites on actin reachable by myosin

*m*are $xml(t)=XAn+1(t)\u2212Xm(t).$ Here, $xml(t)$ has four components: the axial strain $xml,$ the radial spacing between centers of actin and myosin filament $dml$ and the relative angles $\alpha ml$ and $\beta ml.$ Note that for each myosin molecule

*m*there would be $\mathcal{L}a$ binding sites on actin in the neighborhood of

*m*. $\mathcal{L}a$ associated sites on actin can be located on one or two actin filaments, depending on the angle $\beta ml;$ thus, $dml$ can include multiple

*d*

_{M-A}s, one for each associated actin filament. In mass action models, the strain-dependent rate constants between actomyosin states, denoted as (

*k*), exclusively depend on the strain component

_{ij}*x*, whereas all other components are ignored. In the stochastic 3-D sarcomere used here, we adapt the same approach, but for binding we calculate $\mathcal{L}a$ rate constants $k(xml).$ These constants are further modulated by weight factors that take into account the lattice spacing between the filaments,

*d*

_{M-A}, and the azimuthal angles α and β. The resulting $\mathcal{L}a$ binding rate constants per myosin molecule are then used for construction of binding probabilities for the stochastic process as explained below. To match the overall binding flux in mass action models, the binding rate distribution is scaled down by the factor

*f*, which takes into account the mean number of sites on actin reachable by each cross-bridge for the prescribed strain-dependent binding rate function.

_{sk}### Cross-bridge rate kinetics

The modular structure of the computational platform MUSICO can accommodate many different actomyosin cycles. To quantitatively assess the effect of the constrained interaction of myosin and actin in the sarcomere lattice, we consider, without loss of generality, the two simplest cross-bridge models: the two-state A.F. Huxley 1957 model (Huxley, 1957) and the minimal three-state model that can predict fast tension transients (Huxley and Simmons, 1971; Ford et al., 1977). The same kinetic models of myosin binding and the state transition rate constants are used in both mass action kinetics and stochastic simulations of the 3-D spatially explicit model.

#### The two-state A.F. Huxley 1957 model

The two-state A.F. Huxley 1957 model has only two states: detached when the heads can freely move and attached when a myosin head is attached to actin and can only move with actin. The state transition rate constant between the detached and the attached sate is defined as *k _{f}*(

*x*) =

*k*

_{f}_{1}

*x*/

*h*in range 0 ≤

*x*≤

*h*, where

*k*(

_{g}*x*) =

*k*

_{g}_{1}

*x*/

*h*for

*h*≥

*x*≥ 0,

*k*(

_{g}*x*) =

*k*

_{g}_{2}for

*x*< 0, and

*k*(

_{g}*x*) =

*f*·

_{Zah}*k*

_{g}_{2}for

*x*≥

*h*, where

_{Zah}*x*is the axial component of cross-bridge strain,

*h*is the range for positive probability of attachment (Huxley, 1957), and

*f*is the Zahalak factor (Zahalak, 1986), which accounts for higher detachment rates during lengthening for cross-bridges stretched more than

_{Zah}*h*.

_{Zah}#### Three-state model with power stroke

The minimal model capable of predicting transients during isotonic shortening (Lombardi et al., 2004) and the *T*_{1} − *T*_{2} transition (Huxley and Simmons, 1971; Ford et al., 1977) must have at least three states, including a swinging lever arm step or power stroke. For simplicity, we use here a three-state actomyosin cycle (Fig. 5) similar to that proposed by Duke (1999) and Daniel et al. (1998). These two models differ only in the prescribed strain dependence of the state transition rates that are in a large degree arbitrary and designed to fit the experimental data. In Fig. 5 A, we outline a more complete six-state model that includes all essential steps: binding to actin, the power stroke, *d*, associated with Pi release, the ADP release stroke, δ, ATP binding and cross-bridge detachment, and M.ATP hydrolysis, but for simplicity, we can reduced it to the three-state model by combining the states A.M.ADP, A.M, and A.MT together as strongly attached post-power stroke states and M.ATP and M.ADP.Pi together as detached states. The six-state scheme is illustrative in constructing effective state transition rate after combining the low populated states. Because the forward rates are very fast, the equivalent forward rate is set to be limited by ADP release and the reverse rates are negligible because they are very slow. In addition, from the solution data, it can be expected that the binding rate could be limited by rate hydrolysis of M.ATP, but in practice, this step is irrelevant because the second myosin head is most probably already in the M.ADP.Pi state and available for binding after the first head of the cross-bridge detaches.

The physics behind models of strain-dependent thermally activated chemical reaction rates originates from the transition-state theory (Eyring, 1935; Glasstone et al., 1941), as updated by Kramers (1940). To logarithmic accuracy, the rate of a reaction is limited by the energy barrier presented by the point of highest Gibbs energy on the reaction path. The general form of the state transition rates in models with a power stroke is formulated by Hill (1974) as the ratio of forward to backward rates that must satisfy Gibbs’ thermodynamic identity

in terms of the Gibbs energies of the initial and final states, including the elastic strain energy derived from the cross-bridge tension (Hill, 1974; Wood and Mann, 1981). Here, *K _{ij}*(

*x*) is the equilibrium rate constant between states

*i*and

*j*, and each forward or backward rate constant

*K*(

_{ij}*x*) is composed of a strain-independent rate $k ij o , $ which, in principle, is the rate observed in a solution-kinetic experiment under the same conditions, and an

*x*-dependent, i.e., axial strain-dependent, function that is equal to 1 when the molecules are not tethered.

The functions for the free energy of the three states (Fig. 5 B) in the simplest form, after setting the unbound state energy 0, are defined as (Duke, 1999):

where −Δ*G _{bind}* is reduction in free energy caused by myosin binding to actin, −Δ

*G*is a large (negative) change in chemical free energy associated with Pi release, κ is the elasticity of the cross-bridges in pN/nm,

_{stroke}*k*is the Boltzmann constant,

_{B}*T*is absolute temperature in °K, and

*d*is displacement of the lever arm after carrying out the power stroke, i.e., the length of the power stroke in nm.

From Eqs. 2 and 3, the forward and backward state transitions rate constants are defined as follows: For the reaction M.ADP.Pi $\u21c4$ A.M.ADP.Pi, the activation energy is supplied by thermal (Brownian) fluctuations. Following the approach of Kramers (Kramers, 1940; Papoulis, 1991; Hunt et al., 1994; Daniel et al., 1998), the binding rate in quadratic form is derived from a Langevin type of equation balancing thermal fluctuations, elastic restoring, inertial, and viscous drag forces. The elastic element of the detached myosin molecule fluctuates thermally, and the strain-dependent binding rate varies with the axial cross-bridge strain, *x*:

where *k _{bind}* is the overall binding rate of M.ADP.Pi to actin, κ is cross-bridge stiffness,

*T*is absolute temperature, and

*k*is Boltzmann’s constant. The reverse reaction occurs at a constant rate:

_{B}The two attached states A.M.ADP.Pi and A.M.ADP can be assumed to be in dynamic equilibrium if [Pi] > 0 because the forward/backward reaction between them is rapid. Also Pi release is accompanied by a large (negative) change in chemical free energy Δ*G _{stroke}*. Thus, to define the population of these states in some cases, it suffices to specify the ratio of forward and reverse rates (i.e., ratio of the relative probabilities of these post- and prepower stroke states) that also includes the change in elastic energy of the power stroke:

The exponential forms of the forward and backward rate constants can reach huge values and cause numerical problems in calculating the transition probabilities. In these cases, it is convenient to cap the maximums of forward and backward constants and still preserve the equilibrium constants (Smith and Mijailovich, 2008). For example, for a prescribed $k 23 cap $ and $k 32 cap , $ the cross-bridge strain where the capping of large exponentially growing rate constants is $xo=\u22120.5d\u2212[\Delta Gstroke+ln(k23cap/k32cap)]kBT/(\kappa d),$ and the forward rate is

The backward rate, *k*_{32}(*x*), can be simply calculated as the ratio of *K*_{23}(*x*) (Eq. 4b) and *k*_{23}(*x*).

For strain-dependent ADP release (i.e., A.M.ADP → A.M), the difference between the initial and final strain energies gives the forward rate as

where δ is the displacement that the lever arm must move to open the nucleotide pocket and allow ADP release, and $k ADP o $ is the rate of ADP release when the elastic element is relaxed. ATP binding is much faster than the reverse rate *k*_{13}, and thus any strain dependence caused by the conformational change has almost no effect. It is therefore reasonable to take *k*_{13} to be small and constant.

### Monte Carlo simulations of rate dependent stochastic processes

In the stochastic model, we used the standard Metropolis algorithm where a kinetic transition in time step Δ*t* is implemented when a random number in (0, 1) lies in the range (0, *k*Δ*t*), where *k* is the first-order transition rate constant. This algorithm generates a Markov process if *k*Δ*t* ≪ 1, so that at most one transition occurs per Monte Carlo time step in a single subsystem, here considered as one myosin filament and the associated actin filaments. Thus, Δ*t* must be less than the inverse of the fastest rate constant of the system, *k _{max}*, and in practice

*k*Δ

_{max}*t*< 0.001 was required to achieve satisfactory statistics. The transition rates between actin and myosin states are applied to

*N*myosin heads per myosin filament. For all myosin heads,

_{mh}*N*, at any instant, is assigned its state in the actin–myosin cycle, whether it is free or attached to actin. For attached myosin heads, it is also assigned the index of the myosin-binding site on actin and also which strand on the actin filament (1 or 2) it is bound to. In general, both myosin heads of a dimer can be active, but in skeletal muscle, it is occasionally found that both heads attached to different actin filaments. For simplicity, we consider here that only one head of dimeric myosin is active and interacts with actin at a time. Further, we also assume that the active myosin head, usually denoted as a cross-bridge, can bind to any one of the reachable discrete binding sites on actin filament (∼5.5 nm apart, following the unloaded actin spacing). The reachable sites can be on one or two actin filaments depending on angle α and one or both strands on each filament depending on angle β (see Figs. 3 and 4).

_{mh}Each active myosin head (cross-bridge) can be at time *t* in only one state. During a small time step, each myosin head can change its state to neighboring states defined by the actin–myosin cycle or stay in its current state. The Monte Carlo drawing defines whether transitions between actomyosin states will occur or not. In the computational algorithm, for each myosin (cross-bridge), we calculate probabilities for each possible change of its current state to neighboring states as *P* = *k*Δ*t* for each transition, and complement probability up to 1 is attributed to no change of the state. Because each myosin attaches to one of the $\mathcal{L}a$ sites on actin in its proximity, the probability of myosin *m* to attach to *l*th site of $\mathcal{L}a$ available is

where $P j/1 m $ is probability of binding cross-bridge *m* to a single site when it is the only possibility.

In the Huxley 1957 model (Huxley, 1957) there are only two states, attached and detached. In this case, only one probability for the change of state is assigned. For the myosin heads in attached state, the probability of a cross-bridge to detach during time step Δ*t* is equal to *P _{det}* =

*k*(

_{a}*x*)Δ

*t*. For myosin heads in the detached state, the attachment probability is shared between all reachable states defined as

where $Patt(xml)=kf*(xml)\u22c5C\alpha (\alpha ml)\u22c5C\beta (\beta ml)\Delta t$ and $k f * =k f /f sk . $ Because *P _{att}* is the sum of probabilities of attaching myosin heads to each of reachable sites

*l*= 1 to

*l*, the equivalent axial strain–dependent binding rate $k f * $ is set to provide the same flux as the mass action binding rate

_{max}*k*. The scaling factor,

_{f}*f*, decreases the magnitude of

_{sk}*k*at each

_{f}*x*by a mean number of reachable binding sites to myosin heads. The weight factors

*C*

_{α}and

*C*

_{β}are associated with the azimuthal position of actin filaments in the sarcomere lattice relative to the myosin head, angle α, and azimuthal angle of myosin-binding site on actin β (Figs. 3 and 4).

Similarly, in the three-state model, the probabilities of changing state are constructed in the same way, but each state can transition to two neighboring states; therefore, the transitions from the attached states (2 and 3) are defined by probabilities to *P*_{21} = *k*_{21}(*x*)Δ*t* and to *P*_{23} = *k*_{23}(*x*)Δ*t* for the pre-stroke state 2 (A.M.ADP.Pi) and to *P*_{32} = *k*_{32}(*x*)Δ*t* and to *P*_{31} = *k*_{31}(*x*)Δ*t* for the post-stroke state 3 (A.M.ADP). The transition states from the detached state M.ADP.Pi include two attachment probabilities *P*_{12} and *P*_{13} associated with axial strain–dependent rates with multiple binding sites on actin, $k 12 * ( x m l )=k 12 ( x m l )/f sk $ and $k 13 * ( x m l )=k 13 ( x m l )/f sk . $ The scaling factor *f _{sk}* and weight factors

*C*

_{α}and

*C*

_{β}are defined in the same way for all models.

For each cross-bridge, we use one Monte Carlo drawing to define whether the cross-bridge remains in its current state or it will change its state into one of the possible states within current time step Δ*t*. For each cross-bridge, the probability in the range from 0 to 1 is divided into probability bins, *P _{ij}*, in a specified order, including the set of probability bins associated with a cross-bridge changing state and a bin associated with the probability of remaining in the current state. Depending on which bin the drawn random number falls in, the fate of a particular cross-bridge is defined and set for the following time step. The time step, Δ

*t*, is set to be sufficiently small such that at most only one transition can occur per a subsystem so interference between multiple transitions within a single subsystem is avoided and between the systems is negligibly small.

Once the Monte Carlo drawing is completed over all cross-bridges, the connectivity matrix, stiffness matrix, internal (active cross-bridge) forces, external (boundary) forces, and constraints are updated, the iterative solution of the equilibrium equation provides the configuration at time *t* + Δ*t*. The stepwise process is repeated until the maximum number of prescribed time steps is reached.

For simplicity and without loss of generality, all simulations presented here are performed at full calcium activation ([Ca^{2+}] ≅ 2.5 × 10^{−5}), where the thin filament proteins troponin and tropomyosin are assumed to have a negligible effect on the cross-bridge cycle (Mijailovich, S.M., et al. 2015. Biophysical Society 59th Annual Meeting. Abstr. 337a). For these simulations, therefore, we used a version of MUSICO without the calcium regulation module.

### A probabilistic formulation of cross-bridge kinetics in mass action models

The mathematical formulation of sliding filament system at the level of a single myosin head is given in terms of the state probability density function, *p _{i}*(

*x*,

*t*), i.e., the probability of the head being attached to the thin filament (actin) at a binding site displaced by distance

*x*at time

*t*(Eulerian formulation). The subscript

*i*identifies the particular actomyosin state and runs from 1 to the total number of states

*n*used in the cross-bridge cycle. The myosin may be detached or bound to actin. For simplicity, in the mass action simulations presented here, we simplified the spatially distributed elasticity of the filaments into a (discrete) series elastic component and consider myosin binding to actin as in the inextensible filament model. Note that in the MUSICO simulations, myofilament, and cross-bridge compliances are treated explicitly as a function of position along the filaments. The strain-dependent state transitions between myosin states are governed by conservation laws expressed as field equations. In vector form, this system of partial differential equations becomes

_{n}where **M**(x) is the state transition rate matrix. The matrix size, *n* × *n*, is defined by the number of states in the actomyosin cycle, and, in general, this matrix in 3-D depends on the relative spatial position of unstrained myosin heads and their binding sites on actin, **x**∈(*x*, *r*, α, β). We consider here only the axial strain dependence (exclusively of *x*). In this case **M**(**x**) ≡ **M _{x}**(

*x*), where the components of

**M**(

_{x}*x*) are built from rate constants

*k*(

_{ij}*x*) for the transition from state

*i*to state

*j*. The diagonal elements of matrix

**M**(

_{x}*x*) are defined as

where δ* _{il}* = 1 if

*i*≠

*l*and 0 otherwise; and off diagonal terms as $Mijx$ =

*k*δ

_{ji}*, where δ*

_{ij}*= 1 if*

_{ij}*i*≠

*j*and 0 otherwise. The rate constants

*k*(

_{ij}*x*) have finite positive values (>0) only if transition between states

*i*and

*j*are possible, otherwise they are equal to 0. Each component of the vector

**p**(

*x*,

*t*) is the probability of finding the myosin head in one of its

*n*states at time

*t*and at strain

*x*(Huxley, 1957). The operator

*D*/

*Dt*is the material derivative $\u2202/\u2202t\u2212v( t )\u2202/\u2202x, $ where

*v*(

*t*) is the shortening velocity of the actin filament relative to myosin filament. Note that for simplicity, here we assume uniform shortening velocity along the filaments, i.e., $\u2202v/\u2202x=0. $ This is achieved, without loss of generality, by including a series elastic component equivalent to the elasticity of the filaments and apparently rigid filaments in the overlap region where myosin interacts with actin. The normalized shortening velocity,

*V*(

*t*), denotes

*v*(

*t*) in half-sarcomere lengths per second.

The state transition matrix **M _{x}**(

*x*), formed exclusively from rate constants

*k*(

_{ij}*x*), is singular because the sum of the elements in each column is equal to zero, i.e., determinant $|Mx(x)|=0.$ Thus, to find a unique solution of Eq. 5, it is necessary to replace one row in

**M**(

_{x}*x*) by a constraint equation. In the original Huxley 57 model (Huxley, 1957), the constraint is that at any instant of time

*t*, the sum of all state probability densities at particular

*x*is equal to 1, i.e.,

This constraint was later used in almost all mass action models (Hill, 1974; Pate and Cooke, 1989; Smith and Geeves, 1995b,a; Mijailovich et al., 1996).

For Huxley’s two-state formulation, after setting up *p*_{1}(*x*,*t*) as the attached state and *p*_{2}(*x*,*t*) as the detached state and including the above constraint equation, the state transition matrix becomes:

and for the three-state model, after setting *p*_{1}(*x*,*t*) as the detached state (M.ADP.Pi), *p*_{2}(*x*,*t*) as the attached pre-stroke state (A.M.ADP.Pi), and *p*_{3}(*x*,*t*) as the post-stroke state (A.M.ADP), it becomes

### Conservation of the number of cross-bridges (myosin heads)

The assumption in mass action kinetic models that the sum of all state probability densities, *p _{i}*(

*x*,

*t*), at any

*x*is equal to 1 does not guarantee that the sum of probabilities of all states of myosin within actin–myosin cycle or sum of all actin states add to 1. Thus, if the sum of all state probabilities is not equal to 1, one cannot preserve the numbers of species entering or leaving the reaction. Piazzesi and Lombardi (1995) were the first to attempt to partially address this problem, and their approach is successfully implemented in the Hai and Murphy four-state model of smooth muscle contraction (Mijailovich et al., 2000) as an essential condition to preserve the number of myosin species (i.e., cross-bridges).

In mass action kinetic models, we do not track each cross-bridge, but rather the binding of the population of cross-bridges that can reach actin-biding sites within a region $R$. If there is no relative movement between actin and myosin filaments, i.e., in truly isometric conditions, all bound myosins are within the region $R$, forming probability density distributions, *p _{i,att}*(

*x*,

*t*), and the detached probabilities complement these probability density distributions to 1 at any

*x*according to the original Huxley assumption (Huxley, 1957). Because the sum of

*p*(

_{i}*x*,

*t*) at any

*x*is equal to 1, then the sum of all probabilities of the myosin states,

is also equal to 1. Here $lR$ is equal to the length of the binding region $R$. In the original Huxley 1957 formulation, the region $R$ is defined as 0 ≤ *x* ≤ *h*, and the probability of myosin to be attached in fully activated muscle at steady-state is equal to *P _{att}* =

*k*/(

_{f}*k*+

_{f}*k*). The probability of a cross-bridge to be in the detached state is equal to 1 −

_{g}*P*. Piazzesi and Lombardi (1995) in their multistate model defined the biding region $R$ to be within binding range dictated by binding and detachment strain–dependent rate functions. Their approach is equally applicable under isometric conditions to other multistate models, but the proportion of bound and unbound myosins strongly depend on the width of the myosin-binding region $lR$. Thus, defining the region $lR$ needs to be done carefully.

_{att}The above formulations are correct for the isometric case, but when shortening or lengthening is allowed, they are not because some myosin heads are drawn away from the attachment region $R$ while still attached, and at the same time, some unoccupied binding sites on actin are moved within $R$ (Mijailovich et al., 2000). Consequently, the number of available myosins for binding within region $R$ is less than a sum of all attached states,

for the number of the attached myosin heads that were drawn outside of region $R$. Here *n _{s}* is number of states in the actomyosin cycle. Following the approach of Piazzesi–Lombardi (Piazzesi and Lombardi, 1995; Mijailovich et al., 2000), we define ξ as a local coordinate of the binding site on actin available in $R$ (thus 0 ≤ ξ ≤ $lR$). We also assume that after detachment a cross-bridge rapidly (of the order of a microsecond) regains its original configuration; thus, the myosin head can reattach to one site on actin in $R$, and therefore all detached states,

*p*

_{i}_{,}

*(*

_{det}*x*,

*t*), are also only within $R$. Thus, all of the heads detaching at a given

*x*beyond region $R$, regain the position, within $R$, by shifting the change of Δ

*p*

_{i}_{,}

*(*

_{att}*x*,

*t*) to Δ

*p*

_{i}_{,}

*(ξ + $nlR,t$), where subscripts*

_{det}*att*and

*det*denote the attached and detached states, where $lR$ is the length of binding region $R$ and also taken as a spatial period and $nR$ is the maximal number of times the attached cross-bridge has exceeded one boundary of $R$ in the same direction, taking values negative for shortening and positive for lengthening. Thus, the requirement that the sum of the probability densities over all possible states includes translation of detached states, for all ξ, according to the condition

This condition ensures that the number of cross-bridges (myosins) is conserved, but the number of detached states does not precisely reflect the number of the unoccupied myosin-binding sites on actin. The correct number of unoccupied sites on actin can only be achieved by spatially explicit models such as MUSICO.

### Model parameters

#### Sarcomere geometry and myofilament elasticity

In the frog (sartorius) muscle 3-D sarcomere lattice, actin filaments are ∼1 µm long, having 182 monomers in a strand or 364 in total (Burgoyne et al., 2008). The actin monomer spacing is 2.735 nm and a half-period of one strand is 35.55 nm (Huxley et al., 1994; Wakabayashi et al., 1994; Prodanovic et al., 2016). The length of a myosin filament is ∼1.58 µm, having 50 crowns, i.e., 150 myosin molecules per half-thick filament, with crown spacing of 14.3 nm (Luther et al., 2008). The sarcomere slack length is set at ∼2.175 µm, with ∼0.7 µm actin–myosin overlap. The lattice interfilament spacing is *d*_{10} = 37 nm, recalculated to the above sarcomere length from Matsubara and Elliott (1972) and Millman (1998), the actin radius is *r _{a}* = 3.5 nm and myosin radius is

*r*= 7.8 nm (see Table S1 for complete set of model parameters). For simplicity, we limited all stochastic simulations to a half sarcomere with 500 myosin and 1,000 actin filaments. This number of filaments is comparable with the number of filaments in a cross-section of a typical myofibril and provides sufficient statistical averaging without running the simulation multiple times.

_{m}Actin and myosin filaments are extensible with filament moduli (elastic modulus times cross-section area) derived from x-ray diffraction or direct measurement: for actin, K* _{a}* = 0.65 × l0

^{5}pN; and for myosin K

*= 1.32 × 10*

_{m}^{5}pN (Huxley et al., 1994; Kojima et al., 1994). The elasticity of the filaments in mass action models is, for simplicity, lumped into a series elastic component normalized to a subsystem containing one myosin and two actin filaments. The stiffness of the series elastic component in the Huxley and Huxley PL models, equivalent to the actin and myosin filament elasticity in MUSICO, is $K SE Hux $ = 144 pN/nm; and for Duke and Duke PL models is $K SE Duke $ = 198 pN/nm (Table S1).

#### Cross-bridge model parameters

To contrast the inclusion of the geometric and steric constraints of the actomyosin cycle in the explicit 3-D sarcomere lattice versus mass action kinetic models, we compared the simulations of classical mechanical experiments in muscle using the two simplest cross-bridge models: (1) the two-state A.F. Huxley 1957 model (Huxley, 1957) and (2) a three-state, including a power stroke, model similar to that proposed by Duke (1999) (see also Table S2).

#### Cross-bridge model parameters: A.F. Huxley 1957 model parameters

For simplicity, we used here the original parameters of A.F. Huxley’s model: the attachment and detachment kinetic rate constants are *f*_{1}, *g*_{1}, and *g*_{2} equal to 43.3, 10, and 209 s^{−1}, respectively. The myosin distortion displacement scale, *h*, was taken to be 15.6 nm, as reported in Huxley’s original paper (Huxley, 1957). For reaching experimentally observed lengthening maximum velocity, we included a Zahalak factor *f _{Zah}* = 1.8 for

*x*>

*h*= 15.6 nm (Zahalak, 1986). The Piazzesi–Lombardi spatial period is $lR$ =

_{Zah}*h*. In all simulations, cross-bridge stiffness is taken to be κ = 0.58 pN/nm, in order to match maximum isometric tension per myosin filament of ∼530 pN corresponding to the tension observed in frog skeletal muscle of 270 kPa.

#### Cross-bridge model parameters: Three-state model parameters

According to the approach of Duke (1999), the state transition rate constants are as follows: for binding, the equilibrium constant is $Kbind=kbind/kunbind=e\u2212\Delta Gbind/kBT\u224520,$ where Δ*G _{bind}* = −3

*k*and forward rate constant at zero cross-bridge strain is

_{B}T*k*= 170 s

_{bind}^{−1}; for power stroke, equilibrium constant (Eq. 6b) is defined by Δ

*G*= −15

_{stroke}*k*and power stroke

_{B}T*d*= 10.6 nm; and for ADP release/detachment, $k ADP o $ = 56 s

^{−1}and second power stroke δ = 0.9 nm. Because of the exponential forms in Eqs. 4b and 4c, the state transition rates can become very large and can generate numerical problems—they are capped to $k 23 cap $ = 1,000 s

^{−1}, $k 32 cap $ = 100 s

^{−1}, and $k 31 cap $ = 10

^{4}s

^{−1}. These values are chosen as optimal values to satisfy Monte Carlo statistics for time steps of the order of 1 µs. When the cap value is reached, the reverse rates are changed to decay exponentially to satisfy the equilibrium constant,

*K*(

_{ij}*x*). In all simulations, cross-bridge stiffness is taken to be κ = 1.3 pN/nm (as used by Duke [1999]) and the value for

*k*= 4.14 pN·nm at

_{B}T*T*= 300°K. The Piazzesi–Lombardi spatial period is $lR$ ∼ 9.4 nm.

#### Normalization of isometric force

The spatially explicit model in MUSICO calculates the force in each actin and myosin filament, and the force is directly related to stochastic kinetics of the actomyosin cycle and the elasticity of cross-bridges and myofilaments. The simulations are done by setting up the initial conditions and letting the system evolve over time (for example, 1 s with 10-µs steps). Because of the stochastic process of myosin interactions with actin, the forces in the myofilaments fluctuate in time, and each filament experiences somewhat different force. Overall, each simulation with 500 myosin filaments per a half sarcomere provides not only sufficient averaging of the muscle fiber tension but also fluctuations in force in each myosin and actin filament. It is important to note that these forces at any instant of time are different and continuously change in time although the fully developed muscle fiber force under isometric conditions fluctuates very little.

The mean force per myosin filament, *F*, is shown in all graphs where one axis is force. In contrast, the force from mass action kinetic models cannot be directly obtained from kinetic model parameters including the elasticity of cross-bridges and myofilaments. It is typically normalized by a factor that relates the integral of the equilibrium isometric cross-bridge force distribution and maximum isometric tension, *T _{o}* (Huxley, 1957). In our plots, we show both the myofilament force and the isometric tension. The scales are related by a factor that takes into account how many myosin filaments there are per unit of the fiber cross-sectional area (Linari et al., 1998). In intact fibers from frog skeletal muscle, the

*d*

_{10}lattice spacing is in the range of 36.0 to 38.5 nm at a sarcomere length of 2.3 µm (Matsubara and Elliott, 1972; Millman, 1998), and the fraction of cross-section occupied by myofibrils as 0.83 (Mobley and Eisenberg, 1975). For spacing at sarcomere optimal length of 2.175 µm used in our calculations, the number of myosin filaments is ∼505 per µm

^{2}of the muscle fiber cross section, and the mean force per myosin filament,

*F*, is between 530 and 600 pN for corresponding muscle isometric tensions,

_{o}*T*, between 270 and 300 kPa, respectively. For simplicity, all force data are plotted in terms of force per myofilament,

_{o}*F*, and the scaling relationship between the maximum isometric force,

*F*, and maximum isometric tension,

_{o}*T*, is displayed in the figure legends.

_{o}#### Normalization of the cross-bridge distributions

Understanding of the complex muscle fiber response to activation and mechanical load challenges is rooted in the distribution of cross-bridges versus the strain of each cross-bridge and its actomyosin state. The differences in responses between different models for the same activation and mechanical protocols can also be explained by comparing the cross-bridge distributions. The distribution in a spatially explicit model such as MUSICO can be represented as the frequency of cross-bridges in each actomyosin state within a bin of the prescribed width, Δ*x _{b}*, at the mean bin strain

*x*. Thus, these frequencies can be, after summing, related to the mean number of myosin heads in a particular state per myosin filament. The state probability density functions,

_{b}*p*(

_{i}*x*,

*t*), in mass action kinetic models are frequently denoted as distributions of cross-bridge fractions in each state over a range of strains

*x*at time

*t*. An integral of an attached state probability density function represents the probability of a cross-bridge to be in that state. Thus, the problem is in defining the number of cross-bridges in detached states because the probability density functions of these states cannot be integrated over the specified domain of

*x*and therefore cannot be related to the fraction of the bridges in any particular state.

To compare the distributions between mass action kinetic models and MUSICO, we need to limit the range of detached states. Integrating the distributions within these limits can provide the fraction of attached versus detached cross-bridges comparable with the ratio of number of attached versus number of detached cross-bridges calculated from MUSICO. Setting up these limits is not simple and should be defined for each type of actomyosin cycle. Fully developed isometric tension is probably the best state for establishing the normalization factors.

Because in different MUSICO simulations we can use different numbers of myosin filaments, we adjusted the bin width so that the height of the bin represents the number of the cross-bridges in each attached state per myosin filament. The comparative plots of the attached cross-bridge distributions from MUSICO and state probability distribution functions from the mass action kinetic models require the scaling factors for each actomyosin model used.

#### Normalization of the cross-bridge distributions: Huxley 1957 model

The distribution of attached cross-bridges at fully developed force calculated by the Huxley mass action kinetic model is defined as $p bound iso =f 1 /( f 1 +g 1 ) $ over the binding range 0 < *x* < *h* and as $p unbound iso =1\u2212p bound iso . $ Thus, limiting the detached states to the range 0 < *x* < *h* provides the fraction of attached cross-bridges that is proportional to the number of attached cross-bridges calculated from MUSICO. This proportionality factor is used for calculating normalizing factor relating the frequency of cross-bridges per bin to the probability density function in the same attached state.

#### Normalization of the cross-bridge distributions: Three-state model

The probability density distribution of two attached cross-bridge states at fully developed force is more complex than from the Huxley 1957 mass action kinetic model, and defining the equivalent detached state distribution is also much more complex. The problem with this actomyosin cycle model is that at fully developed isometric force, the distributions of bound states significantly differ in shape between the mass action kinetic model and MUSICO predictions. This problem is rooted in loose definitions of the detached state probability density function from mass action models. However, it can be resolved in part by defining Piazzesi–Lombardi binding region length, $lR$, using, for example, the data from MUSICO simulations.

In addition, normalization of the distributions of the attached states requires special procedures, including determination of the proportion factor relating bin heights from MUSICO simulations to the probability density distributions from mass action models. The scaling of frequencies from MUSICO simulations and probability density functions of each attached state is obtained from simulations of isometric force development at 5 and 20 ms where the shape of pre- and post-power stroke attached states distributions almost coincide. The distribution of the sum of bin heights of the two states is excellently fit by a modified Gaussian function (see Fig. S1). This same function also excellently fits the sum of fractions of attached cross-bridge states. After normalization of the peaks, these two fitted curves almost coincide but only at early times during force development. The ratio of these peaks is used to scale the axis for the frequency distributions from MUSICO simulations and the fraction of attached cross-bridges from simulations from Duke and Duke PL models.

### MUSICO software environment and simulation details

The MUSICO software has been developed as C++ object-oriented application that includes LAPACK linear algebra package and deal.II finite element library. Typical run times for these simulations depend on the number of actin and myosin filaments. For the simulation of 500 myosin filaments over 1 s with a time step of 10 µs is ∼20 h on the AEGIS04-KG grid site, consisting of 6 nodes, each equipped with 2 AMD Opteron 6276 16-core processors and 96 GB RAM, totaling 192 processors.

### Online supplemental material

The supplemental text shows normalization of the probability distributions of the attached states. Fig. S1 matches the scales of the number of attached cross-bridges per bin from MUSICO simulations to the fraction of attached cross-bridges from the Duke mass action model prediction at 5 ms after the onset of isometric force development. Fig. S2 shows transient velocities and cross-bridge distributions after quick release to isotonic force at *F*/*F _{o}* =

*T*/

*T*= 0.6 and 0.4. Fig. S3 shows three-state model predictions of transient velocities and cross-bridge distributions after a quick release to isotonic force at

_{o}*F*/

*F*=

_{o}*T*/

*T*= 0.6 and 0.4. Fig. S4 shows two-state model predictions of velocities and distributions for shortening velocities of 63.5 and 215.5 nm/s (i.e., 0.04 and 0.133 of

_{o}*v*) that correspond to MUSICO velocities for

_{max}*F*/

*F*=

_{o}*T*/

*T*= 0.6 and 0.4, respectively. Fig. S5 shows three-state model predictions of velocities and distributions for shortening velocity of 259 and 445 nm/s (i.e., 0.17 and 0.29 of

_{o}*v*) that correspond to MUSICO velocities for

_{max}*F*/

*F*=

_{o}*T*/

*T*= 0.6 and 0.4, respectively. Fig. S6 shows the cross-bridge distributions

_{o}*t*= 0

_{QRel}_{+}, 5, 10, 20, 80, and 200 ms after a quick release to zero isotonic load for model, where black bars represent weakly bound cross-bridges and red bars the post-power stroke bound myosins (MUSICO). Fig. S7 shows three-state model (Duke) predictions of the evolution of the cross-bridge distributions during

*T*

_{1}−

*T*

_{2}transitions after a quick decrease in length of 7.35 nm. Table S1 lists values of key sarcomere lattice parameters. Table S2 lists values of key cross-bridge kinetic parameters and constraints.

## Results

The spatially explicit model implemented in the MUSICO platform directly addresses the binding of myosin molecules to and detachment from actin in the sarcomere lattice. The simulations at any instant of time record all connections between myosin and actin and explicitly take into account how many myosin heads (cross-bridges) are available to interact with accessible sites on actin filaments. Using the MUSICO platform simulations, we have assessed the differences between our spatially explicit model and two versions of mass action kinetic models, the original Huxley 1957 model (Huxley, 1957) and a more realistic three-state model similar to that of Duke (1999). Both of these mass action models have an original version and with the constraint introduced by Piazzesi and Lombardi (Piazzesi and Lombardi, 1995; Mijailovich et al., 2000) that preserves the number of myosin heads, denoted as Huxley PL and Duke PL.

### Two-state (Huxley) model

#### Isometric force development

The Huxley 1957 model is the simplest and is used here to illustrate key differences among spatially explicit and mass action kinetic models. Fig. 6 A shows force development and isotonic shortening after a quick release to zero force (*F* = 0). For comparison, the tension in mass action models is scaled to force per myosin filament (Fig. 6 A); thus, the rescaled maximum isometric tension, *T _{o}*, from mass action models is the same as the maximum isometric force in the myosin filament,

*F*, predicted by MUSICO. During force development, the predicted tensions by Huxley and Huxley PL almost coincide, and only small differences are caused by a small amount of shortening because of the extensibility of actin and myosin filaments (Fig. 6 A). The MUSICO predicted force per myosin filament also matches for the first ∼50 ms, when the population of bound myosins is still low, but when the number of bound cross-bridges increases (

_{o}*t*> 60 ms), the MUSICO predicted force is slightly lower and followed by a long slow rising phase. At

*t*> 500 ms, the MUSICO prediction also plateaus and coincides with the other two models.

The insets in Fig. 6 A display comparisons among cross-bridge distributions from MUSICO and the mass action kinetic models. At 5 and 20 ms, all three models show about the same distributions reflecting negligible differences in force in the early period of force development. The fluctuations in the distributions calculated by MUSICO are caused by favorable or unfavorable spatial matching of discrete positions of myosin heads and binding sites on the actin filament in the 3-D sarcomere lattice. Because of the extensibility of the filaments during force development, a small relative movement between actin and myosin filaments (shortening) is observed even under isometric conditions (Mijailovich et al., 1996). This shortening shifts the distributions leftward and allows additional flux of myosin binding that differs between MUSICO and mass action kinetic models. At 80 ms, a significant number of myosin heads is bound, and some sites on actin are already occupied; thus, the rate of increase in force is a little slower than mass action models. Consequently, the MUSICO predicted force in this region is ∼5% smaller. Giving sufficient time, the slow rising phase of force becomes small and all the net fluxes become smaller. At 700 ms, the distributions of cross-bridges at almost all cross-bridge strains, *x*, become the same.

#### Isotonic shortening

After a quick release from maximal isometric force, *F _{o}*, to isotonic force,

*F*= 0, the shortening response of the half sarcomere is the same in all three models because the series elastic components in the Huxley and Huxley PL mass action kinetic model match the overall elasticity of the extensible filaments used in MUSICO simulations. After a quick release, the distributions shift for about the same amount (Fig. 6 A, inset), exposing parts of the distribution to high binding rates (see inset in Fig. 6 B at a time immediately after quick release,

*t*= 0

_{QRel}_{+}). Although the displacements are the same at

*t*= 0

_{QRel}_{+}, the velocity traces are different (Fig. 6 B). The largest velocity is predicted by Huxley, whereas the MUSICO and Huxley PL have similar values. The reason for these much larger velocities predicted by the Huxley model is the large binding flux in the region of high binding rates where there are no bound myosins. In contrast, the binding flux predicted by MUSICO and Huxley PL is significantly smaller because of the significantly reduced number of myosin heads available for binding (Fig. 6 B, vertical arrows in inset at 10 ms). The resulting difference is best shown at

*t*= 10 ms after the quick release (Fig. 6 B, inset), showing that the cross-bridge distributions for the fraction of newly bound cross-bridges in the Huxley model, displayed at the right hand side of the distributions, is much higher than in MUSICO or Huxley PL distributions and significantly more shifted leftward reflecting higher shortening velocities after the quick release (Fig. 6 B, inset, open horizontal arrows). At later times, when the number of attached cross-bridges is significantly reduced, the velocities predicted by all three models become closer together, having only a minimal difference in steady-state at

_{QRel}*t*= 200 ms after the release.

_{QRel}The predicted force-velocity curves by MUSICO and mass action kinetic models for the same sets of parameters (Huxley, 1957) showed quite different values of isotonic force at low velocity of shortening (Fig. 7 A). For example, for isotonic force *F*/*F _{o}* =

*T*/

*T*= 0.6 MUSICO velocity is significantly smaller than predicted by Huxley and Huxley PL mass action kinetic models (63.5 nm/s, 181nm/s, and 147.5 nm/s, respectively). At

_{o}*F*/

*F*= 0.4, Huxley PL velocity is midway between Huxley and MUSICO. At higher values of isotonic force, Huxley PL velocities are close to Huxley’s, but as the isotonic force drops to low values,

_{o}*F*/

*F*=

_{o}*T*/

*T*< 0.3, the Huxley PL velocities approach the MUSICO predictions, reaching about the same value at

_{o}*F*/

*F*=

_{o}*T*/

*T*≤ 0.2. Interestingly, although the predictions of the original Huxley model (Huxley, 1957) agree very well with A.V. Hill’s force-velocity relationship (Hill, 1938), MUSICO predications and, to a lesser degree, those of Huxley PL underestimate the shortening velocities at the same isotonic force (Fig. 6 A).

_{o}The differences in calculated force-velocity curves originate from neglecting conservation of the numbers of myosin heads and taking into account the occupancy of binding sites on actin in mass action kinetic models. Huxley PL takes into account conservation of the number of myosin heads but still neglects the occupancies of sites on actin. This effect can be seen in the bound cross-bridge distributions at different times during transients after quick release at, for example, *F*/*F _{o}* =

*T*/

*T*= 0.4 (Fig. 8 A). In both cases, immediately after a quick release (at

_{o}*t*= 0

_{QRel}_{+}ms), the bound cross-bridge distributions shift for the same amount in all models, following the same instantaneous (elastic) shortening after release at

*F*/

*F*=

_{o}*T*/

*T*= 0.4 (Fig. 8 A, inset). However, the velocities after release are not the same because the effective flux of myosin binding differs significantly among the models. The cumulative effect is shown in distributions at

_{o}*t*= 20 ms after quick release (Fig. 8 C, left). The fastest growth of the number of newly bound myosin heads, on the right hand side of the distributions, is predicted by the Huxley model, much slower growth is predicted by Huxley PL, and the slowest is predicted by MUSICO. The Huxley binding flux is highest (Fig. 8 C, green arrows) because it does not take into account the reduction in the number of myosin heads already bound and drawn out from the region of binding (0 <

_{QRel}*x*<

*h*) and, therefore, shows the highest shortening velocity (Fig. 8 C, horizontal open green arrows). The Huxley PL model takes into account the fraction of bound heads outside the region 0 <

*x*<

*h*and reduces the effective flux of myosin binding and therefore predicts lower shortening velocity (Fig. 8 C, open red arrows) than Huxley at the same isotonic force. But Huxley PL does not take in account the reduced number of the occupied sites on actin and predicts higher velocity than MUSICO. At

*t*= 200 ms after the release, the distributions reached steady-state, showing about the same pattern for all three models, but the differences in velocities remain because the fine balance of modulated binding flux by attached cross-bridges (Fig. 8 C, arrows up) and the detachment flux (Fig. 8 C, arrows down). The distinct differences in distributions are visible at the right hand side as a reflection of the difference in shortening velocities (for details see Fig. S2). At low isotonic force, the difference in velocities becomes smaller and smaller, reflecting a low number of bound myosin heads and therefore its influence on binding flux and predicted velocities. Consequently, when

_{QRel}*F*/

*F*≪ 0.1, the steady-state distributions almost coincide (at

_{o}*t*≥ 80 ms, Fig. 6 B, inset).

_{QRel}If we now look at the same velocities of shortening as in MUSICO at 211.5 nm/s (*F*/*F _{o}* =

*T*/

*T*= 0.4), Huxley and Huxley PL mass action kinetic models achieved the same velocities at different isotonic forces (Fig. 7 A, dashed pink lines). Consequently, the decrease in a half-sarcomere length is different and reflects the drop in force, whereas the slopes of the displacements are about the same (at t ≥ 0.715 s, Fig. 8 B, insets). Interestingly, Huxley and Huxley PL achieved about the same steady-state velocities, and overall velocities differ slightly in early transient phase (Fig. 8 B), although cross-bridge distributions after release to different isotonic forces displayed different initial shifts (Fig. S4 C at

_{o}*t*= 0

_{QRel}_{+}ms). These different shifts compensate for the higher effective binding rates of Huxley 1957 and Huxley PL models, resulting in about the same binding and detachment fluxes and therefore similar velocities. The cumulative effect of this behavior is best shown at

*t*= 20 ms (Fig. 8 C, right panel). At steady-state (

_{QRel}*t*= 200 ms), the cross-bridge distributions in detachment regions (

_{QRel}*x*< 0) are almost the same in all three cases as a reflection of all having about the same shortening velocity. For details, see Fig. S4.

#### Isotonic lengthening

With the same model parameters, MUSICO predicted more realistic lengthening velocities than Huxley 57 (Fig. 7 B). The Huxley 1957 model predicted a maximum lengthening velocity about three times larger than observed of 1.8 *T*/*T _{o}* (Katz, 1939; McMahon, 1984), whereas MUSICO overshoots experimental observations by only 50%, and Huxley PL slightly more (Fig. 7 B, solid lines and closed symbols). Inclusion of a Zahalak factor of 1.8 brings the predicted lengthening velocities by MUSICO and Huxley PL close to the experimental value of 1.8

*F*/

*F*and 1.8, respectively (Fig. 7 B, dashed lines and open symbols). To achieve the same values, Huxley’s model requires a Zahalak factor of >4 (Fig. 7 B, dashed green line). These data demonstrate that for lengthening it is important to take into account both the conservation of myosin heads and availability of actin-binding sites. However, on closer look, the Huxley 1957 predictions are close to the Katz data (Katz, 1939) at low lengthening velocities (Fig. 7 B), but at lengthening velocities >60 nm/s, the Huxley predictions continue to grow at the same pace, whereas the experimental data level off at ∼1.8

_{o}*T*/

*T*. The MUSICO predictions are even better at the low lengthening velocities (Fig. 7 B, inset) but also show departures to higher

_{o}*T*/

*T*than observed. MUSICO predictions with

_{o}*f*= 1.8 reaches the correct maximum force, but the initial slope of force-velocity curve during lengthening and at low velocities is significantly lower than observed (Fig. 7 B, inset). Overall, these data show that inclusion of the conservation of interacting species and the geometrical constraints significantly improve predictions of lengthening force-velocity relationship, but it may require inclusion of more precise strain dependence of the actomyosin cycle and more than two states to achieve more realistic predictions.

_{Zah}### Three-state (Duke) model

#### Isometric force development

A three-state model including a power stroke is the simplest realistic model that can explain *T*_{1} − *T*_{2} transitions (Huxley and Simmons, 1971). Fig. 9 A shows force development and isotonic shortening after a quick release to *F* = 0. During force development, the predicted isometric forces by MUSICO and mass action three-state models show significant differences (Fig. 9 A). The MUSICO prediction shows a much faster rise of force than the mass action models and has an overshoot peak at ∼80 ms. The force in the Duke PL model also rises faster than Duke’s model and shows consistently higher values. At longer times (*t* > 0.5 s), the predicted force by all three models merges at the same constant value. The MUSICO predicted force decreases after the peak, and Duke PL approaches to a constant value the fastest, whereas the force predicted by the Duke model has a long slow rising phase after initial fast growth.

The differences in the dynamics of force development can be explained by comparing cross-bridge distributions at different instants of time (Fig. 9 A, insets). The largest differences are at ∼80 ms when the MUSICO force reaches a peak. The reason for the fastest growth of force and the peak is caused by faster binding of cross-bridges during early times of development that increases the number of bound cross-bridges and intrinsic shortening of isometric muscle caused by extensibility of actin and myosin filaments. This intrinsic shortening brings a significant number of attached pre-stroke bridges in the region of power stoke firing and therefore contributes to a further increase of force. Because the rise of force is rapid, there is no sufficient time during that period for detachment of the cross-bridges; thus, the force overshoots its steady-state value. However, when a sufficient number of cross-bridges is reached, the effective binding flux decreases, intrinsic shortening stops, and the detachment of cross-bridges becomes faster, whereas the attachment of new cross-bridges becomes slower; thus, force decreases and intrinsic shortening reverses to intrinsic lengthening.

At *t* > 0.5 s, the MUSICO predicted force also plateaus and coincides with the other two models. However, the cross-bridge distributions are quite different. The MUSICO distributions are much narrower than Duke’s. The reason is that Duke’s model does not have the restriction of available myosin heads and binding sites on actin; thus, it continues to bind at higher positive and negative strains and shows much wider distribution of bound cross-bridges. Consequently, Duke predicted force after initial fast growth continues with a long slow growth. Duke PL does not show this behavior because the binding was permitted only in the region $\xb1lR/2$ = ±4.7 nm, reducing the artificial effect of unrestricted binding of Duke’s model that results in a wide equilibrium cross-bridge steady-state distribution. Consequently, Duke PL reaches the plateau fastest and does not have an overshoot, but if not normalized, the Duke PL model has smaller fully developed tension than the Duke model because cross-bridges that can bind are restricted within the region $lR$.

#### Isotonic shortening

After a quick release to unloaded shortening at *t* > 0.7 s (Fig. 9 A, inset), the shortening response of the half sarcomere for all three models is also the same because the series elastic components of Duke and Duke PL models match the overall effect of myosin and actin elasticity used in MUSICO simulations. All three models predicted transient phases after a filament force or tension drop from maximum isometric *F _{o}* or

*T*to the isotonic

_{o}*F*or

*T*: phase 1 reflects the undamped elastic shortening during an instantaneous drop of force, phase 2 shows rapid shortening, which is followed by phase 3 displaying a period of shortening at reduced speed, and finally phase 4, where shortening velocity evolves into steady velocity

*v*(Piazzesi et al., 2002). During the transient phase 3, the duration of the reduced velocity region is shortest in the Duke model, significantly longer in Duke PL, and the longest in MUSICO, causing a delay in the steady-state decrease rate in displacement traces (Fig. 9 A, inset).

_{max}Although the initial shortenings are the same, the velocity traces at *t _{QRel}* = 0

_{+}ms after the quick release are different (Fig. 9 B). The fast transient phase of velocities of each model reaches high values of ∼5 µm/s (Fig. 9 B, inset) and then decrease quickly to low values >1 µm/s for Duke and to significantly lower values (to ∼0.2 µm/s) for Duke PL and MUSICO. The reason for very high velocities just after release at

*t*= 0

_{QRel}_{+}is the rapid transition of a large number of cross-bridges from the pre- to post-power stoke state increasing the force driving the shortening and also decreasing the resistance to shortening (Fig. 9, A and B). At

*t*> 5 ms after the release, a large number of the cross-bridges in the post-power stroke state detaches quickly, significantly decreasing shortening velocity. MUSICO and Duke PL have a longer period at low velocities (phase 3) than does the Duke model. Duke’s cross-bridge distributions at

_{QRel}*t*∼ 10 ms show a large shift to the left and a different profile compared with the MUSICO and Duke PL distributions (Fig. 9 C and Fig. S6). Thus, at this time, the shortening velocity predicted by the Duke model rapidly increases, whereas the shortening velocities of MUSICO and Duke PL still have low velocities for an additional 8–10 ms when they also begin to increase (Fig. 9 B). After a few damped oscillations between fast and slow shortening, all three models reach steady-state velocities (Fig. 9 B).

_{QRel}Similarly to the predictions from models with the Huxley actomyosin cycle, the force-velocity curves for three-state actomyosin cycle showed large differences between MUSICO and the mass action kinetic models (Fig. 10). At higher values of isotonic forces after release, the predicted velocities by the Duke model diverge, whereas MUSICO and Duke PL have the similar values, consistent with observation of Edman (1988). The reason for a good agreement between MUSICO and Duke PL is the carefully chosen width of the myosin-binding region $lR$ = 9.4 nm. However, if the width of the region $lR$ increases, the predicted velocities by Duke PL at higher isotonic forces become larger, approaching the divergent Duke model predictions. The reason for this divergent behavior could be in using state transition rate constants, which are quite different than originally used by Duke (1999).

The predicted velocities using Duke’s original constants show no divergence of the Duke model, the velocities were close to those predicted by Duke PL, and overall, all three models show force-velocity relations similar to the observations of Edman (1988), but only after renormalization (Fig. 10, inset). Also, the differences in model predictions have the same trend as in Fig. 10 but to a lesser degree. MUSICO again predicted lower velocities than Duke or Duke PL for the same isotonic forces. Although the original constants used by Duke provided more reasonable force-velocity curves, they provided only about the half of the expected force per myosin filament (Fig. 10, inset) and a large overshoot during *T*_{1} − *T*_{2} transitions (Fig. 12 A). Thus, we will focus here only on the set of constants that provides very good predictions by MUSICO for all experiments. It is also important to emphasize that at higher isotonic tensions (*T*/*T _{o}* > 0.7) the shortening velocities show instabilities, and only an approximate value of the shortening can be deduced from the stepwise displacements traces.

In the middle range of isotonic forces for the new set of parameters (Fig. 10), for example at *F*/*F _{o}* =

*T*/

*T*= 0.4, the MUSICO predicted velocity is slightly smaller than that predicted by Duke PL but significantly smaller than predicted by Duke (445 nm/s, 543 nm/s, and 763 nm/s, respectively). Similar behavior is observed at lower or slightly larger isotonic force than

_{o}*T*/

*T*= 0.4, but for

_{o}*T*/

*T*> 0.7, Duke velocities diverge more and more from the MUSICO predicted velocities (Fig. 10). Velocity transients after a quick release show similar behavior as for isotonic unloaded shortening, displaying a very fast increase in velocities because of rapid power stroke transitions followed by a large drop and period of low velocities, damped oscillations, and finally reaching steady-state (Fig. 11, A and B). At low levels of isotonic force, MUSICO predicts ∼9% smaller velocities than Duke and 8% smaller than Duke PL (Fig. 10). These differences are much larger than observed in the comparative simulations with the Huxley’s actomyosin cycle (Figs. 6 B and 7 A). The reason for these larger differences is that the number of attached bridges at low isotonic forces (at

_{o}*F*/

*F*=

_{o}*T*/

*T*≈ 0) is much larger in simulations with the three-state model compared with Huxley’s two-state model. This result demonstrates the importance of conserving the number of myosin heads (Mijailovich et al., 2000) and the reduction of the number of available sites on actin for myosin binding.

_{o}As shown for the Huxley 1957 actomyosin binding model, neglecting conservation of myosin heads and occupancy of binding sites on actin in mass action kinetic models could explain the differences in the force-velocity curves in Fig. 11. The effect of the simplifications inherent in these mass action kinetic models can be seen in the bound cross-bridge distributions at different times during transients after quick release at *F*/*F _{o}* =

*T*/

*T*= 0.4 (and at 0.6, shown in Figs. S3 and S5). In both cases, at

_{o}*t*= 0

_{QRel}_{+}ms after a quick release, the bond distributions shift for the same amount in all models, following the same shortening after release at

*F*/

*F*=

_{o}*T*/

*T*= 0.4. However, the velocities after release are not the same (Fig. 11 A, inset) because the effective flux of the transitions between power stroke states, detachment of bound myosin, and the flux of binding myosins significantly differs among the models. The cumulative effect of the history of velocities on the cross-bridge distributions at

_{o}*t*= 20 ms (after the release) shows approximately the same distributions predicted either by MUSICO and Duke PL, and at that time the velocities are in phase 3, whereas the Duke distributions are significantly different, reflecting largely increased shortening velocities that are already in phase 4. Consequently, the distribution is shifted to the left, and the magnitude of post-power stroke state distributions is decreased because of a large flux of detaching post-stroke cross-bridges. At

_{QRel}*t*= 200 ms after the release, the distributions reached steady-state, showing about the same pattern for all three models, but the key difference between the Duke and the other two models is a larger flux of attaching cross-bridges and a larger fraction of both pre- and post-power stroke states (see Fig. S3). Note that both pre- and most post-power stroke cross-bridges contribute to contraction force and therefore higher steady-state velocities. At higher and middle range of isotonic forces, the cross-bridge distributions of MUSICO and Duke PL are similar, and their velocities are also similar. At low isotonic force, the difference in velocities becomes smaller, but the Duke PL distributions become closer to Duke’s distributions and slowly move away from the MUSICO distribution. The velocities follow the same trend.

_{QRel}Similar differences can be observed at the same velocities of shortening as in MUSICO at, for example, 445 nm/s (i.e., 0.29 of *V _{max}* at

*F*/

*F*= 0.4). To achieve the same velocities, each model requires release to different isotonic force levels (Fig. 10). For example, Duke PL achieved the same steady-state velocity at slightly higher, but Duke at much higher, isotonic force (Fig. 10). Consequently, cross-bridge distributions after the release to different isotonic forces display different initial shifts in the cross-bridge distributions. These different shifts create large differences in the flux of myosins binding to actin, causing large differences in shortening velocities at early times after release (

_{o}*t*up to 40 ms; Fig. 11 B, insets). The cumulative effect of the history of these different shortening velocities is best represented at the cross-bridge distributions at

_{QRel}*t*= 20 ms (Fig. 11 C, right), showing approximately the same distributions predicted by MUSICO and Duke PL (phase 3) but significantly different Duke distributions, reflecting much higher binding flux and higher shortening velocities typical for phase 4. When all three models reach phase 4, both the distribution and shortening velocities become similar (see Fig. S5). It is interesting to note that at steady-state the same velocities are achieved at different isotonic forces, i.e., by quantitatively different underlying binding processes.

_{QRel}#### Fast force transients (*T*_{1} − *T*_{2})

The three-state model predictions of the tension recovery after rapid small length shortening per half sarcomere, *y _{o}*, to zero force or tension (

*T*

_{1}∼ 0) also showed significant differences between the models (Fig. 12 A). MUSICO and Duke PL showed rapid force or tension recovery, after initial drop in force and early plateau (

*T*

_{2}), followed by a slow force redevelopment phase to isometric force or tension. The response of all three models almost coincides during ∼1.2 ms after the release but shows large differences during the plateau phase (from 2 to 80 ms after release). During the plateau phase, the balance between fast detachment and reattachment fluxes of cross-bridges defines differences in the time courses of force recovery. The key difference in the response originates in the reattachment flux, which is much larger in Duke’s model. The cumulative effect of this process is best shown in Fig. 12 B at $t T 1 T 2 $ = 5 ms after rapid shortening, where the distribution of the fraction of reattached cross-bridges by Duke model (green lines) far exceeds that predicted by Duke PL (pink lines) and MUSICO (cross-bridge frequencies, black bars), generating much larger force,

*T*

_{2}, overshooting the maximum isometric force or tension,

*T*(Vilfan and Duke, 2003). In contrast, the steep initial fast rise of tension up to 1.2 ms after the release showed almost no difference among the models (Fig. 12 A, inset). This rapid tension recovery immediately after release is predominately caused by dominant transition of the attached cross-bridges from pre- to post-power stokes states and a minor reverse transition flux. These power stroke transitions increase the force with only negligible attachment or detachment of cross-bridges, minimizing the effect of the number of available cross-bridges or occupied sites on actin and therefore in attachment and detachment fluxes. The main difference between the models is in the late phase of rapid force redevelopment and the plateau phase. The difference in the time course of the responses originates from the differences in binding fluxes (Fig. 12 B, arrows). Unrestricted rebinding in the Duke mass action model generates a much larger binding flux (Fig. 12 B, green arrows) than Duke PL that accounts for reduced number of available myosin heads for binding (Fig. 12 B, red arrows) and an even smaller net flux in the MUSICO prediction, which in addition accounts for the occupied sites on the actin filament (Fig. 12 B, black arrows). Consequently, the three models show quite different

_{o}*T*

_{2}values. The much larger binding flux shown at 5 ms (Fig. 12 B, green arrows) causes a much higher rate for the late phase of rapid force redevelopment and an erroneous overshoot of force above the isometric value and an unrealistically high

*T*

_{2}. In the final phase, reattachment and detachment of cross-bridges continues creating distributions of cross-bridges toward the steady-state distribution at isometric force or tension,

*T*, observed during isometric force development. See Fig. S7 for detailed cross-bridge distributions during the transitions.

_{o}## Discussion

The large amount of available structural, biochemical, and biophysical data on muscle contraction provides an extraordinary environment to allow the development of a more comprehensive model that can translate the findings of simple and well controlled experiments into physiological insights to the working muscle in health and disease. The development of the MUSICO platform and other 3-D explicit sarcomeric models (Daniel et al., 1998; Chase et al., 2004; Tanner et al., 2007, 2008, 2012; Williams et al., 2012) provides a significant step toward integration of the broad sweep of current experimental observations. Incorporation of new findings into the model provides an opportunity to unveil important features overlooked in simplified models. Here, we define the effect of the concentration ratio of tethered molecules and their ligands, their geometric constraints, and the occupancy of species on binding kinetics in contracting sarcomeres.

Mass action kinetic models provided important insights for understanding complex muscle behavior. However, the success of these models was limited to explaining typically one or a very few experiments and provided only the apparent relationship between molecular and fiber data. Although such models have been useful in defining important parameters of contraction, the key missing characteristics are (a) number of bound cross-bridges, (b) the relationship between cross-bridge forces and muscle tension, (c) conservation of the number of cross-bridges and actin-binding sites, (d) nonlinear cross-bridge stiffness, (e) local interactions between myosin, actin, and regulatory proteins, and (f) role of sarcomeric accessory proteins, for example titin, nebulin, and MyBP-C. These must be included to achieve a more comprehensive model of contraction and its regulation, which is the long-term goal of this work. This paper represents a first step toward this goal.

In the sarcomere lattice, binding of tethered myosin heads to the available sites on actin is constrained by local geometry and, as such, largely differs from myosin binding in solution or in mass action kinetic models. To quantify these differences, we compared the predictions of the 3-D spatially explicit MUSICO model with the simplest two mass action models for three classical experiments in muscle: force development, isotonic shortening, and fast transients after a sudden change in muscle length. Without loss of generality, these two models provide the most transparent display of quantifiable differences between mass action models and MUSICO simulations.

The key features of our explicit 3-D model of sarcomere contraction include myosin binding to several available sites on actin and explicit accounting of myosin heads in each state of the ATPase cycle. The largest differences between the 3-D explicit model predictions (MUSICO) and mass action kinetic models are in cases when sarcomere length or force suddenly change from the state where a large fraction of cross-bridges are attached or during shortening at medium to high isotonic tension after release. In all cases, the number of myosins available for binding is significantly reduced and large number of sites on actin is already occupied, thus the effective rate of binding, i.e., the binding flux, is also significantly reduced in comparison with the mass action models. The differences are a direct reflection of an incorrect formulation of all mass action models derived from the original A.F. Huxley 1957 model (Huxley, 1957), namely that although the sum of probability densities functions at a given *x* is 1, the global sum of probabilities is not 1; thus, the governing equations are not properly constrained. This results in unrealistic fluxes of myosin rebinding to actin during transient muscle responses or even in steady-state responses when a large number of myosin heads are drawn outside the binding region. The reformulation of the original Huxley approach by Piazzesi and Lombardi (1995) reduces, in part, the effect of the incorrect formulation of the models, but it is still incomplete, and as a consequence, differences still remain. In addition, the availability of multiple binding sites to each myosin head reduces the myosin binding rate in MUSICO by a factor *f _{sk}* = 2.5–3 times from the binding rates in mass action models in order to achieve equivalent binding fluxes to that of the mass action models. This scaling factor,

*f*can be calculated and takes into account the span of actin-binding sites that myosin can reach for a prescribed binding rate function. Notably, this step is essential in translating rates obtained in solution into binding kinetics in muscle fibers.

_{sk}### Isotonic shortening

During isotonic shortening, large differences between models are observed immediately after a quick release at all isotonic forces because during a quick release a significant number of cross-bridges shift out of the binding region, and the effective binding flux strongly depends on the number of available myosin heads and the unoccupied actin-binding sites. MUSICO takes into account these steric restrictions that are not present in original mass action kinetic models. Quantitative differences in the unloaded shortening velocities are displayed in Fig. 6 B and 9 B where the differences are largest immediately after release, i.e., when the number of the bound cross-bridges is the highest and is significantly reduced at unloaded steady-state velocities where the number of attached cross-bridges is significantly reduced. Similar behavior is apparent in force-velocity relationships (Figs. 7 and 10) where the differences in steady-state velocities are higher at higher isotonic tensions and much smaller at low isotonic tensions. Thus, both the transient response and the steady-state force-velocity relation expose the larger differences between the MUSICO and the mass action models in cases when the number of attached cross-bridges is large. Consequently, at larger isotonic forces, both mass action kinetic models achieve the same shortening velocities at much higher isotonic tensions than MUSICO (Figs. 7 and 10). For example, at isotonic forces *F* > 0.3*F _{o}*, the difference is between 20 and 30% of

*F*. Inclusion of the Piazzesi–Lombardi condition somewhat reduces these differences, but at higher levels of isotonic forces, the differences are still large. These results demonstrate that in order to fit the experimental data, the state transition rates and the rates of ATPase cycle should be reevaluated, and new quantitative relationships for the energetics need to be derived. Notably, this may not be possible with simple models like Huxley 1957 because in order to match Hill’s experimental data (Fig. 7), MUSICO may require more than two actomyosin states and more realistic strain dependence of state transition rates than Huxley originally proposed (Huxley, 1957). This is in part achieved with MUSICO predictions with Duke’s three-state cycle (Fig. 10), but for better matching the experimental data, it may be necessary to use a more realistic multistate actomyosin cycle and updated elastic and geometric parameters.

_{o}The oscillations of velocities are observed in both stochastic (MUSICO) and mass action three-state models (Figs. 9 B and 11, A and B) and are caused by instabilities during the transient change of tension from isometric to isotonic force. These fluctuations are also observed in traces of isotonic shortening (Piazzesi et al., 2002). The additional oscillations in MUSICO predictions are caused by combination of the stochastic process of binding, causing the fluctuation in force over the time, and variations in the degree of favorable matching of myosin binding in the 3-D sarcomere lattice, i.e., the effect of periodic structures in myofilament lattice. These effects are amplified in the simplified half-sarcomere system and would become less prominent if more sarcomeric structures in series are included.

Muscle power output, *vT*(*v*), per myofilament (Fig. 13) reflects the differences in velocities between MUSICO and mass action kinetic model predictions. For both actomyosin models, the MUSICO prediction of power output shows lower peaks than the mass action models, and the peaks occur at different velocities. The Huxley PL and Duke PL models are closer to MUSICO prediction, signifying the importance of taking into account the correct kinetic relationship between tethered myosins and binding sites on actin in a 3-D sarcomere lattice. Because the rate constants of the original Huxley model (Huxley, 1957) are derived from fits of A.V. Hill force velocity curves (Hill, 1938), the peak power of the Huxley mass action model coincides with the observed peak at *v* ≈ *v _{max}*/3 (Fig. 13 A, blue dashed line), whereas the MUSICO predicted peak is shifted to the right (Fig. 13 A, black arrow) because of much smaller isotonic velocities and higher values of isotonic forces (Fig. 7 A).

For the three-state model with commonly used rate constants, MUSICO predicted the peak power at the correct velocity (∼*v _{max}*/3) but with peak power slightly lower than the 140 aW value reported in Edman’s data (Edman, 1988; Smith and Mijailovich, 2008). The Duke PL prediction is also in the same range, whereas the Duke predicted peak was shifted to the right and the peak power was greatly overestimated.

Overall, the large number of myosin and actin filaments and the extensibility of the filaments are not sufficient to erase the fluctuations in distributions, caused by spatial mismatching of discrete positions of myosin heads and the binding sites, in a half-sarcomere system. However, in a multiple-sarcomere system, the fluctuations in distributions would be expected to be significantly reduced by inter-sarcomere fluctuations in length caused by variations in force induced by the stochastic binding process.

*T*_{1} − *T*_{2} transients

Similarly, the *T*_{2} response also shows large differences among the models (Fig. 12). The fast recovery phase and the “shoulder” differ not only in magnitude but also in shape. For example, Duke’s original model predicts an overshoot (Vilfan and Duke, 2003), whereas MUSICO and Duke PL showed more realistic behavior. In this case, the different responses are more complex and consist of two processes. The first process, during fast recovery, represents an interplay between cross-bridge and filament elasticity and generates significantly different force responses in MUSICO than in mass action kinetic models because each cross-bridge experiences a different change in strain, Δ*x*, after a rapid shortening of *y _{o}* per half sarcomere, dependent on the position of the cross-bridge along the filaments. This contrasts with mass action models where the change in strain is uniform for all cross-bridges after correction for the effect of series elasticity, $\Delta x=yo*=yo\u2212yoSE.$ Therefore, in MUSICO simulations, different Δ

*x*values after rapid shortening impose variable power stroke transition rates along the myofilaments, and each cross-bridge contributes in a different manner to fast force recovery. During the later stages of force recovery, where the flux of the reattachment of cross-bridges is modulated by the number of attached cross-bridges and available binding sites on actin, MUSICO simulations do not show an overshoot and look more realistic.

In the analysis of the differences in responses between the models, we compare predictions assuming the same (or equivalent) input parameters. We made no attempt to fit any particular set of data; however, the mass action model predictions show the features of the original published modeling results that were fit to experimental data. Each of the presented predictions, however, illustrates significant differences among the models. The Duke mass action model predictions that include an overshoot are reported in Vilfan and Duke (2003). Though, the model that includes the Piazzesi–Lombardi condition (Piazzesi and Lombardi, 1995) predicted only modest overshoot and a more realistic time course than Duke model. The MUSICO predictions show similar trends, but differences between the models are significant. Regarding the spatially explicit models, Daniel et al. (1998) showed similar predictions as MUSICO but reported only one noisy time course, which was difficult to compare.

### Isometric force development

Differences between MUSICO and mass action kinetic models are also observable during isometric force development but to a much smaller degree. Duke’s model shows large differences during the early phase of the force development (Fig. 10) that are primarily caused by poorly defined detached states within the binding region; thus, the distribution of available binding sites on actin does not match the fraction of the detached state. Fixing the binding region to $lR,$ i.e., an equivalent range to Huxley’s *h* provides slightly better match of Huxley PL to MUSICO predictions because it, in certain ways, preserves the number of myosin heads, but the number of available sites on actin remains undefined and still causes large differences. Even at steady-state (*t* ≥ 0.7 s), the differences in distribution of both bound states are quite different: the MUSICO frequency distributions of bound cross-bridges show much sharper peaks than Duke or Duke PL. These differences in distributions are caused by different sampling of unoccupied sites on actin by the available cross-bridges: in MUSICO, each cross-bridge samples all available sites on actin within its reach, whereas mass action models are limited to the available sites on actin at a particular *x*. The excess number of binding sites on actin to the number of the unattached cross-bridges plays a pivotal role in binding showing sharper and narrower peaks in cross-bridge (frequency) distributions than in the state probability density functions in the mass action models.

It is interesting to compare the time courses of the force development. With the three-state model (Duke), MUSICO predicted force reaches a peak during the early phase of the force development (Fig. 9 A). This is caused by the combination of a fast increase in force and instabilities associated with local shortening between the actin and myosin filaments caused by the extensibility of the filaments. In contrast, during force development with the Huxley model, there is almost no difference between mass action and MUSICO simulations (Fig. 6 A). The only small differences in forces are seen in the early phase of development caused primarily by lumped versus distributed filament compliance in a 3-D sarcomere lattice. Looking over the very small differences during the development of bound cross-bridge distributions shows that the original Huxley model conserves the number of myosin heads and available actin-binding sites within the binding range, 0 → *h*. However, the small differences during a fast change in force are caused by differences in binding fluxes resulting from the small amount of shortening during isometric contractions caused by filament extensibility (Mijailovich et al., 1996).

### Number of attached cross-bridges and mean cross-bridge force

One of the shortcomings of mass action kinetic models is their inability to predict the number of attached cross-bridges and estimate mean cross-bridge force. In contrast, the MUSICO simulations continuously trace the number of cross-bridges on each myosin filament. For example, at fully developed muscle tension, MUSICO with Huxley’s two-state model predicts a mean number of attached bridges of 110 per half-myosin filament and a mean cross-bridge force of 4.8 pN. Under the same conditions, Huxley 1957 predicts a fraction of attached cross-bridges of 0.812 that translates to ∼122 attached cross-bridges per half myosin filament and a mean cross-bridge force of 4.35 pN. Both models likely overestimate the number of attached cross-bridges. This overestimation is caused by the specified magnitude and relatively broad strain dependence of Huxley’s original rate constants (Huxley, 1957). MUSICO with the three-state model predicted a lower number of attached cross-bridges, 97 (or ∼65%), and a larger mean cross-bridge force of 6.2 pN. However, the Duke mass action kinetic model is not normalizable; thus, the number of attached cross-bridges is not accessible for comparison. With Duke PL, the number of attached cross-bridges can be estimated by a similar approach to that used for the Huxley 1957 model, but the estimated number of cross-bridges strongly depends on the binding region length, $lR.$

### Comparison with other models

There are only a few models that estimated the number of attached cross-bridges and force per myosin filament. For example, Smith and Mijailovich (2008) predicted ∼110 bound cross-bridges per half myosin filament and a mean cross-bridge force of 5.9 pN. The other estimates vary from 60 to 100 bound cross-bridges and mean cross-bridge forces in range from 4.5 to 7 pN depending on type of muscle, experimental conditions, and temperature (Linari et al., 1998, 2007; Decostre et al., 2005). The general view is that 33% of myosin heads (or 66% of cross-bridges or ∼100 cross-bridges per a half of myosin filament) are attached in fully developed tension. Thus, our predictions using the three-state model are within the expected range.

The elasticity of filaments plays an important role in coupled strain dependence among the attached cross-bridges affecting transition rates of actomyosin cycle and overall contribution of cross-bridge stiffness to the muscle stiffness (Mijailovich et al., 1996; Daniel et al., 1998). During a fast change of force or muscle length, the change of cross-bridge strain Δ*x* is nonuniform because of the coupling between nonuniform strains in actin and myosin filaments and the elasticity of cross-bridges. Thus, the nonuniform cross-bridge strains along the overlap region affect the strain rates along the overlap. In mass action models, for simplicity, the filament extensibility is replaced with a series elastic component, and the change of the cross-bridge strains Δ*x* along the filaments is not taken into account. The effect of nonuniform change of strain Δ*x* along the overlap, however, is small compared with the effect of modulated myosin binding to actin. We performed simulations with rigid filaments in MUSICO versus mass action models without a series elastic component and showed similar differences between MUSICO and mass action model predictions as reported here. Thus, we confirmed that the responses of mass action models were not significantly affected by use of series elastic element in the mass action models instead of explicit filament elasticity.

The stochastic approach of Daniel et al. (1998) could potentially account for the effect of constrained binding, but limiting the system to only two filaments and binding to collinear sites involves a low number of attached cross-bridges and will not show differences as large as reported here. We have not included a large variation in elastic properties of myosin filament as explored by Daniel et al. (1998), but rather we used the elasticity of the filaments taken from the experiment of Kojima et al. (1994) for actin and assessed from x-ray diffraction (Huxley et al., 1994; Wakabayashi et al., 1994) for both actin and myosin filaments. The validity of these input data for the compliances of filaments used in MUSICO were tested by showing that MUSICO predicted x-ray patterns consistent with x-ray data, and these findings are reported in our publications (Prodanovic, M., et al. 2014. Biophysical Society 58th Annual Meeting. Abstr. 768a; Prodanovic, M., et al. 2014. 40th Annual Northeast Bioengineering Conference [NEBEC]. Abstr. 6972910; Prodanovic, M., et al. 2015. Biophysical Society 59th Annual Meeting. Abstr. 442a; Prodanovic et al., 2016).

The current development of MUSICO is not designed for fitting experimental data per se, but rather to incorporate sufficiently accurate representations of muscle structural and kinetic data to predict muscle function without empirical parameters. Here, we presented only a small part of projected powerful possibilities that are currently under development. Because of its modular structure, the model can be made more and more realistic as better and better information is incorporated. As the model becomes more realistic, it can generate emergent properties that can motivate new and more insightful experiments. In particular, the basic modular form of MUSICO provides an opportunity to develop comprehensive multiscale models of various myopathies. The advantage of this approach is that any new application only requires development of the specific parts associated with specificities of the myofilament system of interest. Extending MUSICO to simulate pathological states can be done relatively simply by incorporating known effects of, for example, pathological mutations in myosin heavy chain, titin, or MyBP-C into the appropriate module.

### Conclusions

In summary, 3-D sarcomere model predictions using MUSICO show large differences compared with the predictions of the two simplest actomyosin cycles in force-velocity curves, isotonic and isometric transients, including the velocity transients after quick release and the tension recovery after rapid small length shortening. Including in mass action kinetic models the conservation of myosin heads proposed by Piazzesi–Lombardi (Huxley PL and Duke PL models) partially reduced the differences in overall predicted responses, but at the molecular scale, the predictions are still disconnected from realistic actomyosin interactions in the sarcomere lattice. The origin of the differences in predicted muscle response is rooted in the flawed mathematical description, which does not ensure conservation of species, originally introduced by Huxley (1957) and used later in almost all mass action models. The consequences of this flaw in the model description is best visible in the differences between the myosin binding fluxes predicted by MUSICO versus the mass action models and corresponding differences in instantaneous cross-bridge distributions. Consequently, the kinetic rate parameters from the experiments using the MUSICO versus mass action model predictions show significantly different values. Furthermore, the explicit 3-D sarcomere model simulations provides not only quantitative differences in muscle response but a plethora of additional essential information including (a) the number of attached cross-bridges per myosin filament; (b) the force per myosin filament that matches the value approximately estimated from the fiber tension; (c) the range of cross-bridge forces within reasonable values; (d) the number of cross-bridges in each actomyosin cycle state; and (e) a precise accounting of energy consumption. Thus, inclusion of tethered molecule kinetics, geometrical constraints, and explicit accounting for the occupancies of interacting species is essential for interpretations in terms of interaction energetics.

## Acknowledgments

This project was supported by National Institutes of Health grant R01 AR048776 (S.M. Mijailovich) and grant P41 GM103622 (to T.C. Irving) and British Heart Foundation grant 30200 (M.A. Geeves). We also would like to thank the graceful support of the Mijailovich family, especially Dragica Mijailovic, Esq. (LLM).

The authors declare no competing financial interests.

Richard L. Moss served as editor.