## Abstract

In an activated muscle, binding sites on the thin filament and myosin heads switch frequently between different states. Because the status of the binding sites influences the status of the heads, and vice versa, the binding sites and myosin heads are dynamically coupled. The functional consequences of this coupling were investigated using MyoSim, a new computer model of muscle. MyoSim extends existing models based on Huxley-type distribution techniques by incorporating Ca^{2+} activation and cooperative effects. It can also simulate arbitrary cross-bridge schemes set by the researcher. Initial calculations investigated the effects of altering the relative speeds of binding-site and cross-bridge kinetics, and of manipulating cooperative processes. Subsequent tests fitted simulated force records to experimental data recorded using permeabilized myocardial preparations. These calculations suggest that the rate of force development at maximum activation is limited by myosin cycling kinetics, whereas the rate at lower levels of activation is limited by how quickly binding sites become available. Additional tests investigated the behavior of transiently activated cells by driving simulations with experimentally recorded Ca^{2+} signals. The unloaded shortening profile of a twitching myocyte could be reproduced using a model with two myosin states, cooperative activation, and strain-dependent kinetics. Collectively, these results demonstrate that dynamic coupling of binding sites and myosin heads is important for contractile function.

## INTRODUCTION

When a muscle is activated, the increasing intracellular Ca^{2+} concentration raises the proportion of troponin C molecules that have Ca^{2+} ions bound to their N domain. The binding of Ca^{2+} alters the shape of the troponin complex and causes some of the tropomyosin molecules to alter their position on the thin filament. This movement exposes binding sites on the actin monomers. If myosin heads attach to these sites, the heads can undergo conformational changes and generate force (Gordon et al., 2000; Kobayashi and Solaro, 2005).

This basic mechanism is complicated by the fact that the status of the myosin heads influences the status of the binding sites, and vice versa. For example, a myosin head cannot attach to a binding site that is in the “off” or unavailable position. Similarly, a binding site is less likely to switch off if a myosin head is already bound to it. The binding sites and the myosin heads are therefore dynamically coupled.

The functional significance of this coupling depends on the relative speeds at which the binding sites and myosin heads switch between different states. For example, the activation of a muscle in which the myosin heads cycle much faster than the binding sites change states will be determined primarily by thin filament kinetics. This will no longer be the case if the cycle times of the binding sites and the cross-bridges become comparable. In real muscles, the kinetics of individual binding sites and myosin heads are influenced by the status of near neighbor molecules. This produces cooperative behavior, which adds additional complexity (Fitzsimons and Moss, 2007). Collectively, these effects can make it challenging to predict contractile function from data and/or assumptions about the behavior of individual proteins.

Mathematical models, and associated computer programs, can be used to predict the behavior of dynamically coupled systems and may therefore help to improve understanding of muscle contraction. An appropriate model of muscle might help researchers produce new insights into molecular mechanisms and predict the consequences of potential therapies. For example, researchers could use the model to predict how changes in myosin isoform expression will affect muscle stiffness, or how changes in Ca^{2+} sensitivity might impact twitch contractions. These calculations could then be used to develop new hypotheses and to refine future experiments.

This paper introduces MyoSim, an innovative new model of muscle that includes molecular-level effects and incorporates dynamic coupling of binding sites and myosin heads. MyoSim extends on previous models based on the distribution techniques originally developed by A.F. Huxley (1957) by including Ca^{2+} activation and cooperative effects. This is achieved using techniques similar to those originally described by K. Campbell of Washington State University (Campbell, 1997; no relation to this author). The cross-bridge distribution approach provides greater control of the strain dependence of myosin kinetics than alternative distortion-based methods (Thorson and White, 1983; Razumova et al., 1999; Rice et al., 2008), and means that MyoSim can be used to predict the contractile behavior of muscles under a wide range of experimental conditions. As examples, this paper includes simulations of muscles during isometric activations and release/restretch maneuvers, and also of measurements performed under tension control. All of the source code and the computer programs that have been developed for this work can be downloaded from the supplemental material, along with tutorials and examples that explain how to use the software. Future updates to the MyoSim software will be available at http://www.myosim.org (freely available under the GNU General Public License).

This paper presents simulations of idealized experiments investigating the effects of sudden changes in the intracellular Ca^{2+} concentration (Figs. 1 and 2), and then shows how MyoSim models can be fitted to real tension-recovery data obtained with chemically permeabilized myocardial preparations (Figs. 3 and 5) and unloaded shortening profiles measured using living cardiac myocytes (Fig. 6). The figures show how dynamic coupling influences contractile kinetics and provide additional insights into how twitch contractions are regulated. Other researchers may be able to use the software to explore these effects and to investigate the consequences of dynamic coupling and cooperative mechanisms in their own work.

## MATERIALS AND METHODS

### Overview

MyoSim, a new computational model of muscle, was used to simulate the mechanical properties (for example, force and length) of chemically permeabilized multicellular cardiac preparations and single enzymatically isolated cardiac myocytes. In this work, these mechanical properties were assumed to be equivalent to those of single half-sarcomeres. The force in the preparations was thus calculated as the sum of the active force produced by a single population of cycling myosin heads and a passive force representing the elastic behavior of titin molecules and, in the case of the multicellular cardiac preparations, extracellular collagen filaments.

Because the precise molecular mechanisms underlying muscle contraction are still unclear, MyoSim was specifically designed to be flexible. For example, the software can simulate the behavior of myosin heads cycling through schemes with two, three, or more different biochemical states. Similarly, the user can define how the rate constants governing transitions between these states vary with the length of the actin–myosin links. In many cases, different assumptions will lead to divergent mechanical behaviors, which means that MyoSim can be used to test hypotheses about specific molecular mechanisms. However, all current MyoSim simulations share the following features: (a) Binding sites on actin filaments switch between a state that myosin heads can attach to and a state where attachment is forbidden; (b) Myosin heads can attach to available binding sites and subsequently cycle through a kinetic scheme chosen by the researcher; (c) The software calculates population distributions *n _{i}*(

*x*) that define the number of myosin heads that are attached to actin in the ith attached state with a length x for each time step in the simulation; (d) Active force is calculated by summing the forces produced by the attached cross-brides, that is

*k*is the stiffness of an actin–myosin link, and Δ

_{cb}*x*is the length of the actin–myosin link when a myosin head in the ith-attached state is bound to a directly opposed binding site. For example, in the two-state scheme illustrated in Fig. 1 A, Δx

_{i}_{1}is equal to x

_{ps}.

### Dynamic coupling

Binding sites and myosin heads are dynamically coupled in muscle because the status of the sites influences how the heads will transition between states. Similarly, the status of the heads influences the kinetic behavior of the binding sites. These effects were modeled using coupled differential equations.

Because the MyoSim software allows the user to choose the kinetic scheme for myosin heads, the precise equations depend on the system that is being studied. For a muscle that is held isometric, the equations for the two-state myosin system illustrated in Fig. 1 A are as follows:*A*(*x*,*t*) is the number of myosin heads attached to the thin filament at time *t* with elastic springs of length *x*, *D*(*t*) is the number of myosin heads in the detached state, *k _{1}*(

*x*) and

*k*(

_{−1}*x*) are the strain-dependent rate constants for the attachment and detachment transitions, respectively, and

*N*(

*t*) is the number of myosin heads that can participate in cross-bridge cycling. Because this number is set by the status of the binding sites in MyoSim simulations,

*N*(

*t*) can also be defined as the number of sites that have heads bound plus the number of empty sites that are available for additional heads to bind to. This is the convention that will be used in the remainder of this paper. The system of equations defined in Eq. 2 obviously has to be modified for more complex myosin schemes (three or more cross-bridges states), but the underlying principles remain the same.

The two-state thin filament system assumed in MyoSim models is simpler than the three-state thin filament system proposed by McKillop and Geeves (1993). The simpler system was chosen to reduce complexity and to ensure that binding sites and arbitrary cross-bridge schemes can always be coupled by the following equation,*a _{on}* and

*a*are rate constants,

_{off}*N*is the maximum number of binding sites that heads could potentially interact with (defined by the arrangement of the thick and thin filaments within the half-sarcomere),

_{overlap}*N*is the number of sites that have a head bound, and

_{bound}*k*and

_{plus}*k*are constants defining the strengths of two potential cooperative effects. The first term on the right hand side of Eq. 3 dictates that binding sites activate at a rate that is proportional to the intracellular Ca

_{minus}^{2+}concentration. The second term implies that binding sites cannot switch to the unavailable state if they are already occupied by a myosin head. If

*k*is >0, the third term causes bound myosin heads to recruit additional binding sites to the active state. Evidence for this type of mechanism includes the increase in force at low Ca

_{plus}^{2+}concentrations induced by NEM-S1 (Swartz and Moss, 1992). Finally, the fourth term describes an effect whereby sites that are unavailable for myosin attachment suppress other sites from switching to the available state. This is a simple way of mimicking some of the functional effects of increasing the stiffness of tropomyosin molecules. According to the semi-flexible chain model of Smith et al. (2003), this modification will reduce the thermal fluctuations required for initial myosin binding (Smith and Geeves, 2003).

*N _{bound}* is calculated by integrating the population distributions for each attached myosin state. Thus,

*N*reduces to the total number of myosin heads in state A. For the more complex schemes shown in Figs. 4 and S1, the summation encompasses all of the attached states.

_{bound}### Interfilamentary movement

In addition to changing the force contributed by passive elastic elements, interfilamentary movement alters the active force produced by cross-bridges. As originally described by A.F. Huxley (1957), the attached cross-bridges are modeled as linear springs and change length when the filaments move. This effect was incorporated into the MyoSim calculations by using polynomial interpolation to displace the distributions that describe the number of myosin heads bound in each attached state with a given spring length. As noted previously (Campbell and Moss, 2000), separating the effects of interfilamentary movement and cross-bridge cycling in this manner produces a more stable computational algorithm than direct integration of the appropriate partial differential equations.

Filament compliance effects (Huxley et al., 1994; Wakabayashi et al., 1994) were incorporated by assuming that a half-sarcomere length change of Δx displaced each actin–myosin link by (1/2)Δx (Campbell, 2009). This simplifies the complex realignment of actin-binding sites and myosin heads that occurs in real muscles (Daniel et al., 1998) but is a useful simplification for Huxley-type distribution models (Getz et al., 1998).

### Control modes

MyoSim calculations can be switched between length and tension control modes by flags in the instruction file defining the simulation protocol. This mimics the situation in real biological experiments that control either the length of the preparation or the force that it generates.

In length-control mode, the MyoSim software updates the number of myosin heads and binding sites in each state by integrating the differential equations that describe their kinetics over a small time step Δt. The software then displaces the cross-bridge distributions to account for any interfilamentary movement that occurs during the time step and finally calculates the sum of the passive and active forces. Tension-control calculations proceed similarly except that the final stage of the calculation for each time step uses an iterative procedure to deduce the length of the half-sarcomere at which the sum of the active and passive forces equals the imposed load.

### Computational methods

MyoSim calculations were performed using computer code written in C++ and compiled as a multi-threaded console application for Windows-based machines (Visual Studio 2010; Microsoft). Routines from the GNU Scientific Library (GSL; Galassi et al., 2009) were used for nontrivial calculations. The differential equations describing the kinetics of the myosin heads and binding sites were integrated using a Runge–Kutta method with adaptive step-size control. Root finding (required for simulating tension-control conditions) was performed using Brent’s method.

Parameter optimization and sensitivity analyses were performed using a graphical user interface written in MATLAB (MathWorks). Optimizations were performed as described previously using both simplex search methods and particle swarm techniques (Campbell, 2006a, 2009; Campbell and Moss, 2002). Additional details about the sensitivity tests are described in the Results.

### Availability of software

All of the computer code is freely available under the GNU General Public License and can be downloaded from the supplemental material. The graphical user interface (which can also be downloaded) makes it easier for researchers without programming experience to use MyoSim. We also provide tutorials explaining how to reproduce the results presented in this paper and how to run new simulations. Future updates to the MyoSim software will be available at http://www.myosim.org.

### Experimental data

Figs. 3, 5, and 6 show simulated records that have been fitted to real experimental data. The data plotted in Figs. 3 and 5 were collected using chemically permeabilized myocardial preparations isolated from rats that had been previously injected with streptozotocin. As a result of this treatment, the myocardial samples contained only the β isoform of myosin heavy chain. Summary data describing the functional properties and the myosin heavy chain content of these preparations have been published (Chung et al., 2013).

The experimental data shown in Fig. 6 were collected at 37°C using a living myocyte isolated from a 3-mo-old Sprague Dawley rat. Full details of the experimental technique have been published (Campbell et al., 2013). The smooth Ca^{2+} transient used to drive the simulation was obtained from the experimental record by fitting the recorded values to a function with two exponentials, as described by Rice et al. (2008; Eqs. 54 and 55 in their publication).

The experimental protocols were approved by the Institutional Animal Care and Use Committee at the University of Kentucky. Animals were anesthetized with pentobarbital (50 mg kg^{−1}) before their hearts were removed.

### Online supplemental material

Tables showing parameter values and figures describing additional simulations are provided as online supplemental material. Fig. S1 shows a three-state myosin scheme. Fig. S2 shows the best fit obtained with this type of model to tension recovery data. Table S1 lists model parameters for the simulations in Figs. 1 and 2. Tables S2–S5 list model parameters for the simulations in Figs. 3, 5, 6, and S2, respectively. A computer file that will install the MyoSim model on Windows systems is also available online, along with a copy of the current documentation for the software. The online supplemental material is available at http://www.jgp.org/cgi/content/full/jgp.201311078/DC1.

## RESULTS

### Dynamic coupling

Fig. 1 illustrates dynamic coupling of regulated binding sites and cycling cross-bridges by showing how a very simple two-state cross-bridge system (Fig. 1 A) responds to step increases in the activating Ca^{2+} solution. Two different situations are shown.

The left-hand panels of Fig. 1 (C–E) show predictions for a simple system where cross-bridges cycle faster than binding sites switch between their available and unavailable states. In this situation, force rises at nearly the same rate as the number of available binding sites when the activating Ca^{2+} concentration is increased. This is because the cross-bridges cycle sufficiently quickly that the number of myosin heads that are bound always remains close to equilibrium with the number of available binding sites. As more binding sites are activated, they are rapidly occupied by newly attaching heads.

Similar to real muscles, this quickly cycling system exhibits activation-dependent contractile kinetics. Specifically, the rate at which force rises after the step change in Ca^{2+} concentration increases approximately threefold from 39 s^{−1} for the pCa 5.8 activation to 125 s^{−1} in the pCa 4.5 condition. The rate of force development is activation dependent in the simulations because the rate at which binding sites became available (*k*_{on} [Ca^{2+}]; Eq. 3) increases in proportion with the Ca^{2+} concentration.

The right-hand panels of Fig. 1 (C–E) show simulations for the contrasting situation in which myosin heads cycle more slowly than the binding sites switch between their available and unavailable states. The force responses predicted for this system differ in two important ways from those discussed above. First, force lags behind the number of available binding sites when the activating Ca^{2+} concentration is suddenly increased. Second, the rate at which force develops is not as sensitive to the activating Ca^{2+} concentration. It increases by only 35% (from 5.3 to 7.1 s^{−1}) as the pCa changes from 5.8 to 4.5. Both of these effects reflect the slower myosin cycling rate. The myosin attachment rate in the right-hand panels is simply too slow for the number of bound heads to remain in equilibrium with the number of available binding sites immediately after the step change in Ca^{2+} concentration. As a result, the rate of force development is limited by the rate (*k*_{1}(x); Fig. 1, A and B) at which individual myosin heads attach to binding sites. Because this rate is independent of the prevailing Ca^{2+} concentration, the time course of force development is not dramatically dependent on the level of Ca^{2+} activation.

Further analysis of these simulations shows that both systems develop the same final force levels. This happens because even though the myosin kinetics are different, the ratio of k_{1} to k_{−1} was chosen to be the same for both the fast and the slow cycling systems (see Fig. 1 B and Table S1). Myosin heads thus spend an identical proportion of their time bound to the thin filament even though the mean attached lifetimes in the two cases are very different. This simple demonstration reinforces the fact that systems that exhibit the same steady-state properties can produce very different responses to rapid perturbations.

In real muscles, the coupling between regulated binding sites and cycling myosin heads is further modulated by cooperative mechanisms. Fig. 2 summarizes some of the simplest implications of these effects. The left-hand panels show the force responses predicted for a system where bound heads increase the rate at which binding sites switch to the available state (that is, k_{plus} = 100 s^{−1} and k_{minus} = 0 s^{−1}; Eq. 3). This cooperative effect has the greatest impact on function when the Ca^{2+}-dependent on-rate (*k*_{on} [Ca^{2+}]; Eq. 3) is slow (that is, at a low Ca^{2+} concentration). In this situation, new binding sites are activated each time new myosin heads bind, and force thus creeps progressively upwards. As originally described by K. Campbell (1997), this slows the apparent rate of force development because the system takes longer to reach steady state. In the specific example shown in the left-hand panels of Fig. 2, the rate of force development for the pCa 5.8 condition decreased from 5.3 s^{−1} in the absence of cooperativity to 1.6 s^{−1} when k_{plus} was set to 100 s^{−1}. The cooperative effect also perturbs the steady-state force produced by the system resulting in a leftward shift of the pCa–tension relationship (an increase in Ca^{2+} sensitivity).

The right-hand panels in Fig. 2 show simulations for the contrasting system where the rate at which binding sites switch to the unavailable state increases in proportion with the number of unavailable sites (that is, k_{plus} = 0 s^{−1} and k_{minus} = 15 s^{−1}; Eq. 3). This mimics a biological situation in which binding sites are less likely to activate if their adjacent neighbors on the thin filament are in the unavailable state. The net effects are to depress steady-state force for a given activating Ca^{2+} concentration (reduced Ca^{2+} sensitivity) without substantially affecting the rate of force development.

Together, Figs. 1 and 2 illustrate some of the basic features of a simple myofibrillar system when it is held at a constant length. Even under these simplified conditions, the dynamic coupling of regulated binding sites and cycling myosin heads produces complex mechanical properties. Moreover, predicting how the myofibrils will respond to a biological intervention requires complete quantitative knowledge of the entire system. For example, increasing the attachment rate of single myosin heads can have a large or a negligible effect on the rate of force development depending on how quickly binding sites transition between their available and unavailable states.

### Length perturbations and tension recovery

The behavior of a myofibrillar system can become even more complicated if it changes length. This is because the rate constants that dictate how quickly myosin heads transition between different states must depend on the length of the connection between the thick filament and the binding site. That is, the rate constants are strain dependent. For example, in real muscles, myosin heads will not be able to reach distant binding sites and will therefore be more likely to bind to sites that they are close to. Similarly, a bound head that is stretched by an extremely large movement will detach more quickly than it will if the filaments remain stationary. Note that, for simplicity, this effect was omitted in the simulations of isometric activations shown in Figs. 1 and 2.

The functional consequences of strain-dependent kinetics depend on the speed of the interfilamentary movement. For example, if the filaments move slowly and the cross-bridges cycle quickly, the population distributions will not be markedly perturbed and the length change will not have a large effect on force production. Fast movements, on the other hand, will displace myosin heads from their normal isometric distributions. This can perturb the normal dynamic coupling and can produce additional complexities.

Fig. 3 investigates these issues by showing predicted force responses for a system that is activated in four different Ca^{2+} concentrations and then subjected to a rapid shortening/restretch perturbation. This protocol mimics that used in many laboratories to measure the rate of tension recovery (commonly called *k*_{tr}) (Brenner, 1988). Note that these simulations were performed under length control except during the rapid shortening/restretch maneuver where the force generated by the half-sarcomere was constrained to be greater than or equal to zero. This mimics the real experimental situation in which muscle preparations that are rapidly shortened below their slack length fall into a loop and shorten against zero load until tension redevelops or the preparations are restretched.

Because the behavior of the myofibrillar system is too complex to calculate easily by hand, a computer program was used to search for the set of parameter values (rate constants, cooperative strengths, etc.) that produced the optimal fit to real experimental data. Specifically, multidimensional minimization was used to fit the predicted force responses to tension-recovery records previously measured using myocardial preparations isolated from diabetic rats. These samples were useful for this analysis because, unlike most myocardial preparations from rodents, they only contained a single isoform of myosin heavy chain—in this case, the slower cycling β type (Chung et al., 2013).

The best fits obtained with this approach are compared with the experimental data in Fig. 3 D. The simulated force records lie mostly within one standard error of the mean data value (gray regions) and had an average root-mean-square deviation from the experimental data of 4.9%. Importantly, the simulated records also reproduced the experimentally recorded activation dependence so that the predicted *k*_{tr} values increased from 2.0 s^{−1} in pCa 5.8 solution to 5.5 s^{−1} in pCa 4.5 solution.

Analysis of the simulations shows that the activation dependence of the rate of tension recovery is produced by dynamic coupling of regulated binding sites and cycling myosin heads. During the period of unloaded shortening, cross-bridges are displaced by the interfilamentary movement. Because the rate constant for myosin detachment is strain dependent in this model (see Table S2), this perturbs the cycling kinetics and the number of attached myosin heads falls (Fig. 3 B, dashed lines). This in turn impacts the binding sites because the transition to the unavailable state is no longer blocked by bound myosin heads. As a result, the number of sites in the available state falls (Fig. 3 B, solid lines). Cooperative effects caused by the k_{plus} and k_{minus} terms (Eq. 3) augment this effect. However, because of the Ca^{2+} dependence of the on rate (*k*_{on} [Ca^{2+}] term; Eq. 3), the reduction in the number of available binding sites is proportionally larger and longer lasting at the lower levels of activation than it is in the saturating pCa 4.5 condition.

The net effect of these complex interactions is that in the pCa 4.5 activation, myosin heads can attach to nearly as many binding sites immediately after the restretch as they could before the perturbation. The rate of tension recovery under these conditions is thus limited by the intrinsic kinetics of individual myosin heads. This contrasts with the situation at lower levels of activation, where the muscle partially deactivates during the period of unloaded shortening. Immediately after the restretch, relatively few binding sites are available for myosin attachment, and the rate of tension recovery is therefore limited by the speed at which the thin filaments reactivate.

In summary, these simulations suggest that the rate of force recovery at maximum activation is limited by myosin cycling kinetics, whereas at lower levels of activation, the rate is limited by how quickly the binding sites switch to the available state.

### Limitations of the two-state myosin model

The simulations shown in Fig. 3 are significant because they describe a molecular mechanism that explains the activation dependence of the rate of tension recovery. However, the calculated records do not precisely match the experimental data. For example, the half-sarcomeres shorten too far (to 80% of their initial length) during the imposed perturbation. This is because the myosin attachment rate is sufficiently slow that the interfilamentary movement detaches all of the myosin heads when the muscle starts to shorten. Once this happens, the half-sarcomere is constrained only by its passive elastic properties and shortens by the maximum amount allowed by the protocol.

Similarly, the simulated force records start the tension recovery from a very low level and do not reproduce the large residual forces (∼50% of steady-state force) measured in the real experiments. This is evident in Fig. 3 D where the largest difference between the simulated forces and the experimental data occurs immediately after the restretch. This discrepancy has important implications for molecular mechanisms because the relative residual force has been shown to predict the rate of subsequent tension recovery in a wide variety of experimental conditions (Campbell, 2006b; Campbell and Holbrook, 2007). The experimental correlation implies that residual forces reflect the kinetic status of the binding sites and myosin heads and do not result solely from viscous drag of the myofilaments.

### Residual forces and more complex myosin schemes

To produce realistic residual forces, myosin heads must attach very quickly when the muscle re-lengthens. This means that unless attachment is assumed to be velocity dependent (for example, Campbell and Lakie, 1998; Månsson, 2010), two-state myosin schemes that produce realistic residual forces also predict rates of tension recovery that are too quick. Models based on more complex kinetic schemes can potentially overcome this problem by separating the processes of myosin attachment (which has to be quick to produce residual forces) and myosin force development (which has to be slow to match the rate of tension recovery).

The simplest myosin scheme that separates attachment and force development has three states (one detached and two attached; Fig. S1). However, Fig. S2 shows that the best fits obtained with this model were not satisfactory. One problem is that the muscle does not shorten quickly enough during the length perturbations. This is because the myosin heads that are producing force and driving the motion are not replaced quickly enough to sustain fast shortening.

Additional simulations were therefore performed using a six-state kinetic scheme (Fig. 4) similar to that described by Lombardi and Piazzesi (1990). This scheme assumes that myosin heads generate force in two steps (Fig. 4, k_{3} and k_{4}) and also that heads that are forcibly detached from the intermediate state reattach very quickly to the thin filament (Fig. 4, k_{6}). In biochemical terms, this might correspond to quick rebinding of heads from a detached M.ADP state. Mechanically, the rapid attachment produces a frictional resistance during rapid length changes.

The simulations with the six-state myosin scheme (Fig. 5) fit the experimental data better than the two-state scheme (the average root-mean-square deviation from the experimental traces was reduced from 4.9 to 2.9%) and reproduce both the residual forces and the activation dependence of the rate of tension recovery. Myosin heads can bind to the thin filament quickly enough (Fig. 4, rates k_{1} and k_{6}) to generate appropriate residual forces during the restretch and yet develop force slowly enough (Fig. 4, rates k_{2} and k_{3}) to match the experimental recordings. The scheme is also consistent with biochemical and mechanical data that suggest that force development is a two-step process (Sweeney and Houdusse, 2010).

The six-state simulations thus have significant advantages over those performed with the two-state model. The simulations fit the experimental data more closely and are based on a more realistic biochemical scheme. The main disadvantage of the six-state scheme is that it is more complicated. The calculations take longer to perform, and they require setting values for more model parameters. Fitting the model to experimental data thus becomes more difficult and requires making additional assumptions. In this example, there is therefore a tradeoff between a model that fits the data (particularly the residual forces) well and a model that is simple.

### Unloaded contractions

All of the figures to this point describe protocols in which the length of the muscle system is controlled. As described previously, the MyoSim framework can also simulate experiments performed under isotonic conditions. Fig. 6 shows an example where the external load is maintained at zero to mimic an unloaded contraction. In this case, the model is based on the two-state cross-bridge scheme shown in Fig. 1 A with added strain dependence (Table S4) and is activated by a Ca^{2+} transient (Fig. 6 A) calculated from experimental data measured from an enzymatically isolated rat cardiac myocyte. The calculations assume that the force developed by the attached cross-bridges is always balanced by the passive elastic force produced by titin molecules. As a result, the simulated muscle shortens as it activates and subsequently re-lengthens when it relaxes. Fig. 6 C shows that, after parameter optimization, the system reproduces the experimentally recorded length profile with a root-mean-square deviation between the simulated and recorded traces of 3.8%.

### Sensitivity analysis

Fitting MyoSim models to experimental data produces numerical estimates for the parameters (rate constants, cooperative strengths, etc.) that the calculations are based on. It is also possible to deduce confidence limits for these estimates by performing sensitivity analyses. These calculations can be very complicated if the model has many parameters. However, a practical approach is to assume that the parameters are mathematically independent and then adjust each of them in turn until the simulated record deviates from the experimental data by a predefined amount (Press et al., 2007).

Fig. 7 shows the results of this type of analysis for the simulations shown in Fig. 6. Each plot shows that the fit to the experimental data becomes worse as the parameter value is progressively increased or decreased from its optimal value. The values when the fit deviates from the optimum by 5% are used as confidence limits. This type of analysis is useful because it shows that the model is more sensitive to some parameters than others. For example, the confidence limit for a_{on} (the parameter that defines how quickly binding sites are activated by Ca^{2+}; Eq. 3) is 0.17%. In contrast, the values of k_{−1,a} and k_{−1,b} (parameters that define the strain dependence of myosin detachment) can deviate by 17 and 53%, respectively, without substantially affecting the fit. This suggests that subtle changes to calcium sensitivity will have a larger impact on unloaded shortening profiles than modifying the strain dependence of myosin kinetics.

## DISCUSSION

This paper shows that dynamic coupling of regulated binding sites and cycling myosin heads can modulate the mechanical properties of striated muscle, and that this effect may help to regulate the contractile function of striated muscles. The mechanical effects of dynamic coupling can be predicted using mathematical modeling, and the MyoSim framework described in this paper was designed to help with this process. Researchers may be able to use the model in their own work to gain insights into myofibrillar mechanisms, to develop new hypotheses, and to refine future experiments.

### Challenges to selecting a model

One of the challenges when developing simulations that include myofibrillar-level effects is choosing the right model. Many physiologists would agree with the frequently quoted statement, “The model should be as simple as possible, but no simpler.” In practice, this is a difficult requirement to meet. Missing data are a key problem. For example, it is not yet clear whether myosin heads in real muscles have a linear spring constant, as originally suggested by Huxley and Simmons (1971), or exhibit nonlinear behavior, as suggested by laser trap experiments (Kaya and Higuchi, 2010). Because models that make different assumptions about molecular properties will produce divergent predictions, it is difficult to identify the simplest model that is consistent with a given set of experimental data. Changing the assumptions might lead to a simpler solution.

A second challenge to choosing the simplest model is deciding exactly what the model should be able to simulate. For example, scientists who focus primarily on the rate at which muscles develop force might choose to start modeling their results with a very simple, nonstrain-dependent two-state kinetic scheme. This is the approach that was chosen by Brenner in his initial pioneering study (Brenner, 1988). However, models lacking strain-dependent kinetics will not produce realistic force transients during length perturbations because the kinetics of the myosin heads will not be influenced by the movement. This means that researchers who wish to simulate real mechanical experiments need to incorporate strain-dependent effects. Even here, the situation becomes complicated. Fig. 3 shows that adding strain-dependent kinetics and cooperative mechanisms to a two-state myosin system can explain the activation dependence of the rate of tension recovery. However, a six-state cross-bridge scheme is needed to explain the large residual forces measured immediately after the restretch (Fig. 5).

It is also possible that cross-bridge transition rates depend on factors other than strain. For example, Månsson (2010), Yadid and Landesberg (2010), and Yadid et al. (2011) have suggested that cross-bridge transition rates may depend on the velocity of interfilamentary movement. Similarly, studies using single molecule techniques, including those by Visscher et al. (1999), Fisher and Kolomeisky (2001), Vilfan (2005), and Veigel et al. (2003), suggest that some transition rates may be load dependent. These possibilities complicate the interpretation of cross-bridge models.

Deciding how many states to include in muscle models is also difficult. For simplicity, MyoSim assumes that binding sites on the thin filament can exist in one of two states (available or unavailable for myosin heads to attach to). However, users have complete control over the kinetic scheme for myosin molecules. A simple system with just two myosin states seems to be sufficient to predict some unloaded shortening profiles (Fig. 6). However, schemes with multiple attached states are probably required to investigate the effects on contraction of ADP and Pi. This is because these metabolites directly influence transitions between different bound states (Sweeney and Houdusse, 2010), and a model that did not incorporate these effects would not be sensitive to changes in the local ionic environment.

Even more complicated models are required to explain the molecular basis of cooperativity. MyoSim can incorporate cooperative effects using the k_{plus} and k_{minus} parameters (Eq. 3), but it does not explain the molecular mechanisms by which they arise. Spatially explicit models can contribute to this shortfall because Daniel et al. (1998) and several subsequent groups (Chase et al., 2004; Smith et al., 2008; Tanner et al., 2008, 2012) have shown that realignment of extensible thick and thin filaments can augment myosin binding and contribute to cooperative activation. As shown by Mijailovich et al. (2012), electrostatic interactions between actin and tropomyosin are also important for cooperativity and can be simulated using mathematical techniques developed for continuous flexible chains.

Spatially explicit models use stochastic numerical techniques and are at the leading edge of myofibrillar modeling. When implemented carefully, they can provide superb fits to experimental data and offer important insights into fundamental molecular mechanisms. However, they do have some disadvantages. Spatially explicit models require sophisticated algorithms, take longer to run on computers, and necessitate setting values for more numerical parameters. The MyoSim model is much simpler and was developed by extending A.F. Huxley’s original cross-bridge distribution technique (Huxley, 1957) by adding Ca^{2+} dependence and cooperative effects. In many cases, this may provide a suitable compromise between utility and complexity.

In short, selecting the simplest model in muscle physiology requires deciding what the goal of the work is and which assumptions can be made. Neither of these tasks is straightforward.

### Simulations of tension recovery measurements

In 1994, Landesberg and Sideman (1994) showed that a loose-coupling model can explain the activation dependence of the rate of force development. (In this context, loose coupling means that myosin heads can exist in force-generating states when the associated troponin molecule is not binding Ca^{2+}.) Subsequent work by K. Campbell from Washington State University also showed that cooperative effects can produce activation-dependent rates of force development (Campbell, 1997). The current work extends these studies by presenting simulations that integrate Ca^{2+}-dependent activation, cooperative mechanisms, and strain-dependent myosin kinetics. The new simulations are significant because they show that the time course of the force record after a rapid shortening/restretch maneuver is critically dependent on the behavior of the muscle during the mechanical perturbation. In the two-state model (Fig. 3), the number of attached cross-bridges declined markedly during the period of unloaded shortening because the strain dependence of the myosin kinetics (Table S2) dictated that heads were being pulled off the thin filaments faster than they could reattach. When the muscle was restretched to its original length, force development therefore started from a value near zero. The situation was different for the six-state model because myosin heads could reattach quickly through the *k*_{6} transition (Fig. 4). This produced a frictional resistance that allowed the half-sarcomere to shorten at a constant velocity during the latter stages of the perturbation. This is consistent with the behavior recorded in real experiments and differs from that predicted by the two-state model, where the half-sarcomere shortens unrealistically far once the myosin heads have detached and the system is constrained only by its passive elastic properties.

The simulations for the six-state cross-bridge scheme also suggest why previous experiments have shown that the rate of tension recovery correlates with the relative residual force, that is, the force measured immediately after the restretch divided by the isometric force for the condition (Campbell, 2006b; Campbell and Holbrook, 2007). In the model, the residual force results from cross-bridges that have been strained during the rapid restretch of the half-sarcomere. Because myosin heads can attach very quickly through the *k*_{6} transition, the number of attached cross-bridges is limited primarily by the number of binding sites in the available state. Thus, the magnitude of the residual force is an indirect measure of the number of activated binding sites. If the residual force is large, most of the binding sites are already available and, as argued previously, the rate of subsequent force development will be limited by the intrinsic kinetics of myosin heads. If, on the other hand, the residual force is low, many binding sites are unavailable for attachment, and the rate of force development will be limited by how quickly the thin filament reactivates.

### Time-varying Ca^{2+} concentrations

Another advantage of the MyoSim model is that it can be used to simulate the effects of time-varying Ca^{2+} concentrations. Figs. 1 and 2 show simple tests where the pCa of the activating solution undergoes a rapid step change near the beginning of the simulation. Fig. 6 shows a more complicated situation where the model is driven by a Ca^{2+} transient derived from experimental data. In this simulation, the half-sarcomere was mechanically unloaded so that it changed length during the experiment, shortening as the muscle activated, and re-lengthening as it relaxed. Fig. 6 D shows that after parameter optimization, the simulated half-sarcomere length record matched that recorded in the original experiments. If the length of the half-sarcomere had been constrained in the simulations (which is achieved by a simple setting in the computer code), the calculations would have predicted the force developed during an isometric twitch.

### Future directions

MyoSim was developed to try and help produce insights into muscle contraction at a myofibrillar level, and to provide a tool that researchers may be able to use in their own work. The computer code is freely available and, because of its modular nature, should be relatively easy to incorporate into larger frameworks. It will be relatively straightforward, for instance, to link MyoSim half-sarcomeres into parallel chains. This will facilitate the investigation of emergent properties arising from inter–half-sarcomere interactions (Telley et al., 2006; Campbell, 2009; Stoecker et al., 2009; Campbell et al., 2011; Rassier, 2012). It may also be possible to integrate MyoSim into finite element models of entire muscles. This might improve the fidelity of these systems by incorporating sarcomere-level effects such as short-range stiffness (Hu et al., 2011). Finally, sensitivity analyses can predict molecular-level effects that have particularly large impacts on function. This data might help scientists and researchers develop better therapies for muscle diseases.

## Acknowledgments

The author thanks Mihail I. Mitov and Charles S. Chung (both University of Kentucky) for performing the experiments with the chemically permeabilized and living myocyte preparations, respectively, and members of his laboratory for many helpful discussions.

This work was supported by the University of Kentucky Research Challenge Trust Fund and National Institutes of Health (NIH) grants HL090749 and UL1 TR000117. The content is solely the responsibility of the author and does not necessarily represent the official views of the NIH.

The author declares no competing financial interests.

Richard L. Moss served as editor.

- Submitted: 5 August 2013
- Accepted: 16 January 2014

This article is distributed under the terms of an Attribution–Noncommercial–Share Alike–No Mirror Sites license for the first six months after the publication date (see http://www.rupress.org/terms). After six months it is available under a Creative Commons License (Attribution–Noncommercial–Share Alike 3.0 Unported license, as described at http://creativecommons.org/licenses/by-nc-sa/3.0/).