An implementation of the model and input data is available on our group website [www.cemrg.co.uk](http://www.cemrg.co.uk)

Introduction {#sec001}
============

Tension generation in cardiac muscle is a highly cooperative process, with significant increases in tension caused by relatively small increases in the calcium concentration. The Hill coefficient (*n* ~*H*~) describing the degree of cooperativity of the force-calcium relationship is typically around *n* ~*H*~ = 3 in experiments on skinned muscle cells \[[@pcbi.1004376.ref001]\], and as high as *n* ~*H*~ = 10 in intact cells \[[@pcbi.1004376.ref002]\]. Our understanding of the molecular mechanisms giving rise to this cooperative activation and the precise regulation of tension generation required for effective cardiac pump function remains incomplete. However, there is a general agreement on the potential types of interactions involved in cooperative activation between regulatory units (RU) and crossbridges (XB) \[[@pcbi.1004376.ref003]--[@pcbi.1004376.ref005]\]. Each half-sarcomere in a myocyte contains 26 RU's, and each RU consists of 7 actin monomers, one long tropomyosin molecule spanning the actin monomers, and a complex of troponin (troponin I, troponin C and troponin T) which regulates local activation. Within an RU, calcium (Ca^2+^) bind to troponin C (TnC), causing a conformational change in tropomyosin, unblocking actin for myosin crossbridge (XB) binding \[[@pcbi.1004376.ref006]\].

Underlying cooperative activation, three types of interactions are proposed between regulatory units and crossbridges. Cooperative effects between RU's are known as 'RU-RU cooperativity', where unblocking of tropomyosin in one RU leads to an increased probability of unblocking in a nearby RU, due to overlap of tropomyosin molecules between neighbouring RU's. Evidence in support of these effects includes experimental data which shows a significant decrease in cooperativity when the overlap between neighbouring tropomyosin units is removed or reduced \[[@pcbi.1004376.ref007]--[@pcbi.1004376.ref009]\], dependence on nearest neighbour interactions in cardiac muscle \[[@pcbi.1004376.ref010]\], and modifications to long-range cooperativity by phosphorylation of tropomyosin \[[@pcbi.1004376.ref011]\]. In addition there are cooperative interactions in which the binding of crossbridges increases the rate at which further crossbridges bind, known as 'XB-XB cooperativity' \[[@pcbi.1004376.ref012], [@pcbi.1004376.ref013]\]. Although XB-XB interactions can increase the steady-state force per activated RU, more complex interactions with neighbouring RU's are involved in their effect on cooperativity. These more complex interactions by which crossbridges affect RU activation are known as 'XB-RU cooperativity' \[[@pcbi.1004376.ref013], [@pcbi.1004376.ref014]\]. Evidence for the importance of these effects on muscle activation can be seen from various experiments in which calcium sensitivity is affected by changes to crossbridge affinity using crossbridge inhibitors and enhancers \[[@pcbi.1004376.ref001], [@pcbi.1004376.ref015]--[@pcbi.1004376.ref017]\]. A potential factor in XB-RU cooperativity are the proposed effects of tension generation on the affinity of TnC for Ca^2+^ \[[@pcbi.1004376.ref018], [@pcbi.1004376.ref019]\]. The mechanisms and significance of this interaction remain controversial, with some researchers claiming this effect appears mainly from non-physiological rigor crossbridges \[[@pcbi.1004376.ref001]\], while others point to it as a key component of normal muscle function \[[@pcbi.1004376.ref018], [@pcbi.1004376.ref020], [@pcbi.1004376.ref021]\].

In addition to uncertainties in the biophysical basis for cooperativity, the exact link between calcium binding to TnC and the movement of tropomyosin has remained obscure, troponin I (TnI) is known to play a key role in transmitting this signal \[[@pcbi.1004376.ref022]--[@pcbi.1004376.ref024]\], and in recent years this link has been clarified with research on crystal structures of troponin \[[@pcbi.1004376.ref025], [@pcbi.1004376.ref026]\]. These studies showed that calcium binding to TnC opens up a hydrophobic patch on TnC which has a high affinity for the switch region of TnI \[[@pcbi.1004376.ref027]\]. The movement of the switch region also moves the nearby inhibitory ('C-terminal') region of TnI which is responsible for pinning tropomyosin in the blocking position on actin in resting muscle \[[@pcbi.1004376.ref028]\]. The competitive binding of these TnI regions to both TnC and actin results in the unblocking of actin at higher Ca^2+^ concentration, allowing myosin crossbridges to bind and generate force. Further support for the critical role of TnI is given by its numerous phosphorylation sites and role in regulating muscle function through *β*-adrenergic stimulation and the response to length-dependent activation \[[@pcbi.1004376.ref028]--[@pcbi.1004376.ref030]\]. [Fig 1](#pcbi.1004376.g001){ref-type="fig"} gives an overview of a regulatory unit (RU) and its states in this competitive binding framework.

![Overview of RU states.\
Each regulatory unit (RU) along the actin-tropomyosin thin filament contains troponin C, which binds calcium, and troponin I, which can bind either to actin or troponin C. These reactions, in addition to the binding of crossbridges, define each of the states. In the schematic, state names with + have calcium bound, state names which include an 'X' have crossbridge(s) bound, and B, U, S refer to the labels 'blocked', 'unblocked', 'stable unblocked' at the bottom of the schematic. 'Blocked' refers to troponin I bound to actin, which blocks myosin binding. 'Unblocked' refers to troponin I being not bound to actin, and 'stable unblocked' refers to troponin I being held in place by troponin C. Each state allows for crossbridge binding, although this is very improbable in the 'blocked' states, such that states ('BX' and 'BX^+^') rarely occur. Note that all transitions between neighbouring states exist, in addition to transitions between the top and bottom rows. The green arrows indicate the main pathway during activation, with Ca^2+^ binding to TnC, and TnI moving from actin to TnC⋅Ca^2+^ to allow crossbridge binding. Red arrows indicate the main deactivation pathway, TnI detaching from TnC⋅Ca^2+^, followd by Ca^2+^ detaching from TnC and TnI binding to actin to block crossbridge binding.](pcbi.1004376.g001){#pcbi.1004376.g001}

In the three-state framework proposed by McKillop and Geeves \[[@pcbi.1004376.ref031]\] RU's can be either in the 'blocked' state with TnI pinning tropomyosin to actin, in the neutral 'closed' state where myosin crossbridges are able to bind, or in the 'open' state with crossbridges having moved tropomyosin in the opposite direction compared to TnI binding. The continuous flexible chain models represent the spatial deformation of tropomyosin along the whole thin filament. At points along the chain with a TnI binding site or crossbridge, the chain is in a fully 'blocked' or 'open' position respectively. However, in the space between bound sites, the chain can occupy a continuum of intermediate states. Due to the single TnI binding site per RU, we can still unambiguously refer to an RU as blocked based on TnI binding. Describing an RU itself as 'open' becomes more problematic in this modelling framework, as there are 2--3 crossbridges per RU and any combination of these can be bound to actin at any one time. In the rest of this paper we refer to the state of RU's only as 'blocked' and 'unblocked' based on TnI-actin binding, regardless of the tropomyosin deformation induced or number of crossbridges bound near the RU.

There are several challenges in applying these advances in physiology to create a computational model of cardiac contraction that is both tractable for a wide range of simulation and analysis, and captures the critical physiological features of the underlying proteins. Firstly, computational models which include tension-dependent feedback mechanisms often suffer from non-physiological hysteresis, in which tension generation is higher for decreasing calcium compared to increasing calcium \[[@pcbi.1004376.ref032]\]. Secondly, in the absence of a clear mechanistic explanation for cooperativity, computational models based on ordinary differential equations (ODE) tend to use phenomenological representation of cooperativity to achieve adequate tension development \[[@pcbi.1004376.ref005], [@pcbi.1004376.ref033]--[@pcbi.1004376.ref036]\]. Some recent developments have begun to address these shortcomings, including detailed models of the thin filament based on spatial interaction of tropomyosin \[[@pcbi.1004376.ref037], [@pcbi.1004376.ref038]\]. Firstly, the work by Campbell et al. includes a model of tropomyosin interaction between neighbouring RU \[[@pcbi.1004376.ref038]\], and is based on ODEs. However, it is limited to approximately 9 RU's, and requires the assumption that calcium bound to TnC does not unbind in the tropomyosin 'closed' state. Extending the model beyond these assumptions quickly leads to an increase in the required number of states beyond what is computationally tractable. Nevertheless, the model is arguably the most biophysically detailed contraction model to have been applied in the context of a whole-organ cardiac mechanics \[[@pcbi.1004376.ref039]\]. Secondly, a more detailed underlying model of cooperativity is given by models of tropomyosin as a continuous flexible chain, based on the work by Smith et al. \[[@pcbi.1004376.ref040]--[@pcbi.1004376.ref043]\]. These approaches assume that tropomyosin, which consists of many molecules overlapping end-to-end to form a long filament, can be modelled as a homogeneous flexible chain. The deformation of the tropomyosin chain in these models is determined by a combination of weak electrostatic interactions with actin and elastic deformation of the chain. Although still a simplification that ignores potential inhomogeneities arising from end-to-end overlap, these models provide a more detailed description of tropomyosin kinetics which are able to describe longer range interactions in the thin filament, compared to models which assume only nearest-neighbour interactions. Solving these more detailed models remains computationally challenging, and results are typically given by developing approximations for the equations of the deformation of tropomyosin, or applying stochastic approaches to predict a steady-state force-calcium relationship.

Our goal in this article is to create an ODE-based model of cardiac contraction with a biophysically detailed representation of cooperativity based on the competitive binding model of troponin I and the continuous flexible chain model for tropomyosin. The formulation of an ODE-based model facilitates modelling of a wide range of simulations of dynamic function of muscle, and will allow us to link this model to whole organ mechanics in the future.

This paper is organized as follows: We start with a general theory on modelling tropomyosin as a continuous flexible chain and the use of Boltzmann's law. The section "Steady-state models" describe our model for the steady-state blocking and unblocking of RU's in the absence of myosin crossbridges. We extend this model to include myosin crossbridges, and develop techniques to make this approach computationally tractable. This extended steady-state model is then used to explain the sources of cooperativity, and the effects on myosin binding in producing XB-RU effects and changes in Ca^2+^-TnC affinity. The section "Dynamic models" develops the dynamic models of cardiac contraction, which we use to investigate the role of cooperative activation in isometric tension development across different species, as well as the influence of cooperative effects on the rate of tension redevelopment and its dependence on Ca^2+^.

Models {#sec002}
======

Our framework combines a model of the deformation of the tropomyosin filament with the more typical protein-protein interactions. This is accomplished using Boltzmann's law, which says that given a system of molecules with different states *S* ~1~, *S* ~2~, ... *S* ~*n*~ with corresponding energies *E* ~1~, *E* ~2~, ... *E* ~*n*~, the probability *P*(*S* ~*i*~) of being in a state that has energy *E* ~*i*~ when the system is in thermal equilibrium is $$P\left( S_{i} \right) \sim e^{- \frac{E_{i}}{k_{\! B}T}}$$ where *k* ~*B*~ is Boltzmann's constant and *T* the absolute temperature in Kelvin \[[@pcbi.1004376.ref044]\].

For our model, we consider interactions between TnI, TnC, Ca^2+^, actin, and myosin and the deformation of tropomyosin to be the significant interactions \[[@pcbi.1004376.ref043]\]. We introduce four constants to represent differences in free energy related to the different protein-protein interactions that regulate cooperative tension development: *E* ~*C*~ is the energy required for Ca^2+^ binding to TnC to form TnC⋅Ca^2+^, *E* ~*A*~ is the energy required for TnI binding to actin to form TnI⋅A, *E* ~*I*~ is the energy required for TnI binding to TnC⋅Ca^2+^ to form TnI⋅TnC⋅Ca^2+^, and *E* ~*M*~ is the energy required for myosin binding to actin. For example, *E* ~*A*~ is the difference in free energy between the state with TnI bound to actin (TnI⋅A) and unbound from actin (TnI+A). This constant can be linked to the ratio between occupation of the states and the dissociation constant *K* ~DA~ via the Boltzmann term: $$e^{- \frac{E_{A}}{k_{B}T}} = e^{- \frac{({E_{\text{TnI} \cdot \text{A}} - E_{\text{TnI} + \text{A}}})}{k_{B}T}} = \frac{P\left( \text{TnI} \cdot \text{A} \right)}{P\left( {\text{TnI} + \text{A}} \right)} = \frac{1}{K_{\text{DA}}}$$ Our model assumes that binding of TnI to TnC in the absence of Ca^2+^ is improbable enough to be negligible, as the hydrophobic patch on TnC is not opened. In addition this implies that the unbinding of Ca^2+^ from TnI⋅TnC is similarly negligible, because of thermodynamic consistency with the high energy required to form TnI⋅TnC in the absence of Ca^2+^. This assumption leads to the absence of the states 'S' and 'SX' in [Fig 1](#pcbi.1004376.g001){ref-type="fig"} and the corresponding transitions ('UX → SX' and 'SX^+^ → SX'). However, note that the unbinding of TnI from TnC and subsequent unbinding of Ca^2+^ is still possible in unblocked RU's, it simply leaves TnI bound to neither actin nor TnC as in the unblocked states (middle column) in [Fig 1](#pcbi.1004376.g001){ref-type="fig"}. Likewise, TnI can unbind from actin, leading to unblocking of RU's and even binding of myosin crossbridges in the absence of Ca^2+^, corresponding to the transitions B → U → UX in [Fig 1](#pcbi.1004376.g001){ref-type="fig"}. This is improbable under normal circumstances, but has been observed experimentally in conditions of low ATP \[[@pcbi.1004376.ref045]\].

Unlike similar stochastic frameworks in which each transition can be handled separately \[[@pcbi.1004376.ref043]\], for our ODE-based approach we combine all components to give the total free energy of a half-sarcomere. Combined with a model which gives the free energy associated with tropomyosin deformation *E* ~tm~, this total free energy is given by: $$E_{\text{tot}}\left( i,j,k,l \right) = E_{\text{tm}} + \left( n - i \right)E_{A} + jE_{M} + k\left( E_{I} + E_{C} \right) + lE_{C}$$ Where *n* is the number of RU's, and *i* = *N* ~*u*~(tm) the number of unblocked RU's for a tropomyosin state tm (and thus (*n* − *i*) the number of blocked RU's),*j* = *N* ~xb~(tm) the total number of crossbridges bound for a tropomyosin state tm,*k* is the number of RU's with both Ca^2+^ and TnI bound to TnC (states S^+^, SX^+^ in [Fig 1](#pcbi.1004376.g001){ref-type="fig"}),*l* is the number of RU's with Ca^2+^, but not TnI bound to TnC (states B^+^, BX^+^, U^+^, UX^+^). We represent the energy related to tropomyosin deformation (*E* ~tm~) using a continuous flexible chain model, which approximates all end-to-end connected tropomyosin molecules as a single chain. The displacement of tropomyosin is regulated by troponin complexes in 26 RU's per half-sarcomere spaced 38.5 nm apart \[[@pcbi.1004376.ref046]\]. With respect to tropomyosin deformation, regulatory units can be either in the 'blocked' state with TnI pinning tropomyosin to actin at -25°, or in an 'unblocked' state resulting in an angle determined by neighbouring units, and tending towards the neutral 0° position. In addition, myosin crossbridges displace tropomyosin in the opposite direction to TnI, at +10°. The chain will assume a minimal energy configuration constrained by these 'fixed points' introduced by TnI and myosin. More formally, the deformation of this chain is represented by the angle *ϕ*(*x*) by which tropomyosin is displaced from its neutral position in the helical groove of the actin filament, and we solve the energy minimization problem: $$\text{minimize} E_{\text{tm}}\left( \phi\left( x \right) \right)\text{with}\text{Dirichlet}\text{boundary}\text{conditions}\phi\left( x_{1} \right) = \phi_{1},\phi\left( x_{2} \right) = \phi_{2}\,\ldots$$ The boundary conditions represent stable points that force tropomyosin to have specific angles at specific locations along the actin filament introduced by troponin I (*ϕ*(*x* ~*i*~) = −25°) and myosin binding (*ϕ*(*x* ~*i*~) = +10°) \[[@pcbi.1004376.ref047], [@pcbi.1004376.ref048]\]. Solving this minimization problem results in the deformation *ϕ*(*x*) along with the free energy *E* ~tm~ for a certain tropomyosin 'state' dependent only on these stable points. The full description of the equations and solution using a finite element model is given in the Supporting Information ([S1 Text](#pcbi.1004376.s001){ref-type="supplementary-material"}, [S1 Fig](#pcbi.1004376.s002){ref-type="supplementary-material"}).

As we are interested in the global properties of the thin filament which determine RU unblocking and force generation, rather than the probability of calcium being bound to any specific RU, we sum over all possibilities of Ca^2+^ or TnI being bound in RU's, neither of which affect the deformation of tropomyosin. Specifically, for a tropomyosin state tm with *i* = *N* ~*u*~(tm) unblocked RU's, there are $\left( \frac{i}{k} \right)$ ways to have *k* TnI bound to TnC⋅Ca^2+^. In addition there are $\left( \frac{n - k}{l} \right)$ ways to have *l* Ca^2+^ bound to TnC without TnI being bound. This allows us to calculate the probability *P*(tm) of being in the thin filament state tm as: $$\begin{aligned}
{P\left( \text{tm} \right)} & {\sim \sum\limits_{k = 0}^{i}\sum\limits_{l = 0}^{n - k}\begin{pmatrix}
i \\
k \\
\end{pmatrix}\begin{pmatrix}
{n - k} \\
l \\
\end{pmatrix}e^{- \frac{E_{\text{tot}}{({i,j,k,l})}}{k_{B}T}}\mspace{180mu}\text{where}\mspace{180mu} i = N_{u}\left( \text{tm} \right),j = N_{\text{xb}}\left( \text{tm} \right)} \\
 & {= e^{- \frac{E_{\text{tm}}}{k_{B}T}}\frac{1}{K_{\text{DA}}^{n - i}}\frac{1}{K_{\text{DM}}^{j}}\left( {1 + \frac{\left\lbrack \text{Ca}^{2 +} \right\rbrack}{K_{\text{DC}}}} \right)^{n - i}\left( {1 + \frac{1}{K_{\text{DI}}}\frac{\left\lbrack \text{Ca}^{2 +} \right\rbrack}{K_{\text{DC}}} + \frac{\left\lbrack \text{Ca}^{2 +} \right\rbrack}{K_{\text{DC}}}} \right)^{i}} \\
\end{aligned}$$

Where the equality follows by applying the binomial theorem (for details see [S1 Text](#pcbi.1004376.s001){ref-type="supplementary-material"}) and the relation between the energies and dissociation constants (see [Table 1](#pcbi.1004376.t001){ref-type="table"}). This result can also be understood more intuitively by considering there are *n* − *i* RU's which have nothing bound, or Ca^2+^ bound to TnC (states 'B' and 'B^+^' in [Fig 1](#pcbi.1004376.g001){ref-type="fig"}), and *i* RU's where there is a TnI⋅TnC⋅Ca^2+^ state ('S^+^') in addition to the states with nothing or Ca^2+^ bound (states 'U' and 'U^+^'). As in [Fig 1](#pcbi.1004376.g001){ref-type="fig"}, crossbridges can be bound in each of these states, although the term *E* ~tm~ makes some of these combinations less probable (e.g. those corresponding to states 'BX' and 'BX^+^'). Note that although states with detached TnI exist, where TnI is bound to neither TnC nor actin ('U' states in [Fig 1](#pcbi.1004376.g001){ref-type="fig"}), they are predicted by our model to be transient and unpopulated (at ∼ 1%) as both *K* ~DI~ and *K* ~DA~ are small.

10.1371/journal.pcbi.1004376.t001

###### Index of model parameters.

![](pcbi.1004376.t001){#pcbi.1004376.t001g}

  Parameter                                                                         Description
  --------------------------------------------------------------------------------- ------------------------------------------------------------
  $K_{\text{DA}} = e^{- \frac{E_{A}}{k_{\! B}T}}$                                   Dissociation constant for troponin I binding to actin
  $K_{\text{DI}} = e^{- \frac{E_{I}}{k_{\! B}T}}$                                   Dissociation constant for troponin I binding to TnC⋅Ca^2+^
  $K_{\text{DC}} = e^{- \frac{{\lbrack\text{Ca}^{2 +}\rbrack}/E_{C}}{k_{\! B}T}}$   Dissociation constant for Ca^2+^ binding to TnC
  $K_{\text{DM}} = e^{- \frac{E_{M}}{k_{\! B}T}}$                                   Dissociation constant for myosin binding to actin
  *k* ~A+~, *k* ~A-~                                                                On- and off-rate for troponin I binding to actin
  *k* ~I+~, *k* ~I-~                                                                On- and off-rate for troponin I binding to TnC⋅Ca^2+^
  *k* ~C+~, *k* ~C-~                                                                On- and off-rate for Ca^2+^ binding to TnC
  *k* ~M+~, *k* ~M-~                                                                On- and off-rate for myosin binding to actin
  *γ*                                                                               Scaling parameter for tropomyosin properties

Steady-state models {#sec003}
-------------------

### Steady-state model of thin filament kinetics {#sec004}

We start by developing a model which ignores crossbridge binding, and only calculates the number of unblocked RU's. For this case *j* = 0 in [Eq 4](#pcbi.1004376.e007){ref-type="disp-formula"}, and the tropomyosin state 'tm' is defined only by the points where TnI is bound to actin and moves tropomyosin to the 'blocked' position in each of the *n* = 26 RU's, giving rise to 2^26^ ≈ 67 million states. We can further group the tropomyosin states from [Eq 4](#pcbi.1004376.e007){ref-type="disp-formula"} by their number of unblocked RU's, while setting the number of crossbridges bound to *j* = 0. As only *E* ~tm~ is dependent on the specific tropomyosin state, and the other terms are only dependent on the number of unblocked RU's, this results in: $$P\left( {N_{u}\left( \text{tm} \right) = i} \right) \sim \text{SE}_{i}\mspace{180mu}\frac{1}{K_{\text{DA}}^{n - i}}\ \left( {1 + \frac{\left\lbrack \text{Ca}^{2 +} \right\rbrack}{K_{\text{DC}}}} \right)^{n - i}\mspace{180mu}\left( {1 + \frac{\left\lbrack \text{Ca}^{2 +} \right\rbrack}{K_{\text{DC}}K_{\text{DI}}} + \frac{\left\lbrack \text{Ca}^{2 +} \right\rbrack}{K_{\text{DC}}}} \right)^{i}$$ $$\text{SE}_{i} = \sum\limits_{\{\text{tm} \mid N_{u}{(\text{tm})} = i\}}e^{- \frac{E_{\text{tm}}}{k_{\! B}T}}$$

Thus, despite the large number of thin filament states, the *n* + 1 constants (SE~0~ to SE~*n*~) are sufficient to determine the relation between calcium concentration and (un)blocking of RU's in the absence of crossbridges. Furthermore, these constants are sums of *E* ~tm~ terms that can be readily determined from the flexible chain model of the thin filament. Although this process requires calculating 2^26^ finite element solutions to *E* ~tm~, this is computationally tractable.

### Steady-state model including crossbridges {#sec005}

The model described by [Eq 5](#pcbi.1004376.e012){ref-type="disp-formula"} does not yet include the effect of myosin crossbridges, which also displace the tropomyosin filament. Including these effects will be key in predicting the effects of XB-RU and XB-XB interactions. Crossbridges binding and unbinding to actin-tropomyosin are represented similarly to TnI binding to actin-tropomyosin, using a dissociation constant and an effect on tropomyosin deformation through the term *E* ~tm~. For a full model which includes XB-XB and XB-RU interactions through the effects of both troponin I and crossbridges binding to actin on the deformation of the tropomyosin filament, we can use [Eq 4](#pcbi.1004376.e007){ref-type="disp-formula"} to determine the probability of being in a state with *i* unblocked RU's and *j* crossbridges bound as: $$P\left( {N_{u}\left( \text{tm} \right) = i \land N_{\text{xb}}\left( \text{tm} \right) = j} \right) \sim \ \text{SE}_{i,j}\ \frac{1}{K_{\text{DA}}^{n - i}K_{\text{DM}}^{j}}\ \left( {1 + \frac{\left\lbrack \text{Ca}^{2 +} \right\rbrack}{K_{\text{DC}}}} \right)^{n - i}\mspace{180mu}\left( {1 + \frac{\left\lbrack \text{Ca}^{2 +} \right\rbrack}{K_{\text{DC}}K_{\text{DI}}} + \frac{\left\lbrack \text{Ca}^{2 +} \right\rbrack}{K_{\text{DC}}}} \right)^{i}$$ $$\text{SE}_{i,j} = \sum\limits_{\{\text{tm}\mspace{600mu} \mid \mspace{600mu} N_{u}{(\text{tm})} = i,\mspace{600mu} N_{\text{xb}}{(\text{tm})} = j\}}e^{- \frac{E_{\text{tm}}}{k_{\! B}T}}$$

[Eq 7](#pcbi.1004376.e014){ref-type="disp-formula"} is similar to [Eq 5](#pcbi.1004376.e012){ref-type="disp-formula"}, but now depends on both the number of unblocked RU's and the number of crossbridges in a half-sarcomere. The SE~*i*,*j*~ values denote the sum of Boltzmann terms for all tropomyosin states with *i* RU's unblocked and *j* crossbridges bound. In this initial investigation we do not consider details of sarcomere geometry and filament overlap effects, but simply model myosin binding sites as evenly spaced every 14.5 nm along the 1001 nm long filament \[[@pcbi.1004376.ref043], [@pcbi.1004376.ref049]\], resulting in *m* = 69 potential crossbridges per half-sarcomere. Including the displacement of tropomyosin introduced by these cross-bridges results in 2^69^ myosin states for each of the 2^26^ configurations of the thin filament RU's. Unlike in the previous section, this large number of states makes the full computation of the state space computationally untraceable, and an approximation is required.

### Sampling crossbridge states {#sec006}

To solve for the thin filament activation kinetics in the presence of XB's requires evaluation of the 26 × 69 values SE~*i*,*j*~. Evaluating these values requires the solution of 2^26^ ⋅ 2^69^ tropomyosin states, which is not tractable. We combine two techniques for approximating these terms without requiring a brute-force calculation.

Firstly, in exploring smaller models we found that the crossbridge binding properties of the thin filament are dominated by the number and length of adjacent stretches of unblocked RU's. Thus, if 'B' and 'U' indicate blocked and unblocked RU's respectively, the crossbridge binding properties of the states 'UBBUUU' and 'BUUUBU' can be well approximated as identical. However, unlike some previous models \[[@pcbi.1004376.ref038]\], crossbridge binding properties can not be inferred from only the number of unblocked RU's, such that (e.g.) the states 'BUUUBU' and 'UUBBUU' have significantly different crossbridge binding properties even though they both have four unblocked RU's. [Fig 2](#pcbi.1004376.g002){ref-type="fig"} illustrates these example states. This reduction results in 3010 *classes* of thin filament states which are equivalent in terms of crossbridge binding properties. We designate the state with the Boltzmann term $e^{- \frac{E_{\text{tm}}}{k_{\! B}T}}$ closest to the mean of the class as the 'representative state' for that class. For example, the 'fully blocked' and 'fully unblocked' states have their own class, the class with one unblocked RU represents 26 states (one for each position), and the class with two disconnected unblocked RU's represents 300 states (each pair of disconnected positions). The standard deviation of the free energy within a class was on average 0.06% and in the worst case 0.11% of the free energy of the representative state of a class. This procedure allows us to perform subsequent calculations only on the 3010 representative states, instead of all 2^26^ thin filament states, but does not reduce the large number of crossbridge states, and still leaves us with 2^69^ cross-bridge states to be calculated for each of these representative classes. Additional results for the representative classes and states are shown in [S2 Fig](#pcbi.1004376.s003){ref-type="supplementary-material"}.

![Equivalent tropomyosin states.\
The top two states are considered equivalent, and part of the same class of states, as they both have three adjacent unblocked RU's and one isolated unblocked RU. The bottom state is part of a different class, even though it also has four unblocked RU's, as it has two stretches with two adjacent unblocked RU's each.](pcbi.1004376.g002){#pcbi.1004376.g002}

Secondly, for each of these representative classes, we approximate the sums of Boltzmann terms in [Eq 8](#pcbi.1004376.e015){ref-type="disp-formula"} by using a Monte Carlo approximation of the sum by random sampling. Technical details of this sampling procedure are described in [S1 Text](#pcbi.1004376.s001){ref-type="supplementary-material"}.

### Independent crossbridge approximation {#sec007}

An alternative strategy for reducing the number of crossbridge states assumes that crossbridge binding does not significantly affect RU unblocking or binding of further crossbridges. Specifically, the model includes RU-RU cooperativity, and each crossbridge binding is affected by the current state of RU's and corresponding tropomyosin deformation, but the resulting deformation of tropomyosin after crossbridge binding does not affect RU kinetics or other crossbridges. Creating a model which assumes XB-RU and XB-XB effects are negligible provides a baseline for comparing results of more detailed models and in determining the importance of XB-XB and XB-RU effects. Given this assumption, we can solve for the tropomyosin deformation without any crossbridges bound to get the free energy *E* ~tm~, and with a single crossbridge bound to get $E_{\text{tm}}^{\prime}$ and the difference in tropomyosin free energy due to crossbridge binding $\Delta E_{\text{tm}} = E_{\text{tm}}^{\prime} - E_{\text{tm}}$. Using Boltzmann's law, the *duty ratio* of a crossbridge, or the fraction of time it is expected to be bound to actin, is given by: $$\frac{P\left( \text{xb}\text{on} \right)}{P\left( \text{xb}\text{on} \right) + P\left( \text{xb}\text{off} \right)} = \frac{1}{1 + K_{\text{DM}}\mspace{600mu} e^{\frac{\Delta E_{\text{tm}}}{k_{\! B}T}}}$$ Where *P*(xb off), *P*(xb on) denote the probability of a specific crossbridge xb being off or on, and *K* ~DM~ the dissociation constant for myosin as used previously. Tension generation in this simplified model is proportional to sum of duty ratios of all potential crossbridges for a representative tropomyosin state.

### Parametrization of the steady-state model {#sec008}

The steady state model has five parameters: four dissociation constants and a parameter for the finite element model of the continuous flexible tropomyosin chain. In this section we determine these parameters for the full 26 RU model based on cooperative activation in intact muscle at body temperature, to reproduce general cooperative activation as observed across different species. Whenever possible we give parameters to a number of significant digits which reflects parameter sensitivity and experimental constraints.

We start by setting the dissociation constant of myosin *K* ~DM~ based on the average duty ratio of crossbridges in a fully activated thin filament (determined by the average fraction of crossbridges bound), which is largely independent of other dissociation constants.

Estimates for the duty ratio of myosin vary from 5--10% for myosin heads in an actin-activated myosin ATPase assay \[[@pcbi.1004376.ref050]\], 14% in-vivo based on power-stroke distance \[[@pcbi.1004376.ref044]\], and up to approximately 30% at high force in experiments on human muscle fibers \[[@pcbi.1004376.ref051]\]. We set *K* ~DM~ = 2, which results in a crossbridge duty ratio of 25% at full activation (1000 μM Ca^2+^), which is consistent with previous modelling work \[[@pcbi.1004376.ref038]\]. This value is at the higher end of experimental measurements but includes XB-XB interactions, and represents a myosin dimer rather than an isolated head.

Next, we determine the dissociation constants for competitive binding of TnI to TnC⋅Ca^2+^ (*K* ~DI~) and actin (*K* ~DA~) along with the scaling parameter *γ* of tropomyosin properties (bending stiffness and electrostatic interactions, see [S1 Text](#pcbi.1004376.s001){ref-type="supplementary-material"}). These parameters all influence cooperativity, while the remaining parameter *K* ~DC~ only affects calcium sensitivity. Varying *K* ~DI~ and *K* ~DA~ between 10^−4^ and 0.1 shows the maximum Hill coefficient is *n* ~*H*~ ≈ 4 when *γ* = 1, rising to *n* ~*H*~ ≈ 6 for *γ* = 2 and *n* ~*H*~ ≈ 8 for *γ* = 3. Based on measurements of the Hill coefficient below and above Ca~50~ (*n* ~2~, *n* ~1~ respectively), and the average Hill coefficient *n* ~*H*~ required to replicate dynamic muscle function in phenomenological models \[[@pcbi.1004376.ref035], [@pcbi.1004376.ref036]\], we use *γ* = 2.

[Fig 3](#pcbi.1004376.g003){ref-type="fig"} shows results for force-calcium relationships as a function of *K* ~DI~ and *K* ~DA~. This parameter sensitivity study shows that there is a relatively large triangular region in which cooperativity is high. We choose parameters sufficiently far away from regions where activation is impaired as indicated in [Fig 3B](#pcbi.1004376.g003){ref-type="fig"}, such that contractile function is maintained even if *K* ~DA~, *K* ~DI~ are varied (by e.g. phosphorylation of TnI). Additionally, given the lack of a clear lower bound, we consider very high affinities (i.e. very low *K* ~DA~, *K* ~DI~) to be physiologically less plausible due to larger differences in free energy between states. Within these constraints we choose *K* ~DI~ = 4 ⋅ 10^−3^ and *K* ~DA~ = 10^−3^, which results in high cooperativity (*n* ~2~ = 7.5, *n* ~1~ = 2.7) consistent with experimental data \[[@pcbi.1004376.ref052], [@pcbi.1004376.ref053]\],

![Influence of TnI affinity for actin and TnC on muscle cooperativity.\
Panel A shows cooperativity plotted as a function of the dissociation constant of TnI for actin (*K* ~DA~) and the dissociation constant of TnI for TnC⋅Ca^2+^ (*K* ~DI~). There is a relatively large triangular region in parameter space in which cooperativity is high, with a slight tendency for higher cooperativity at very low *K* ~DA~, *K* ~DI~ reflecting more extreme competitive binding of TnI. Panel B shows calcium sensitivity, which follows a smooth gradient The yellow 'X' indicates our choice of parameters, and the red contours indicate the regions within which *n* ~*H*~ ≥ 5. At high *K* ~DA~, affinity for actin is insufficient to block tropomyosin effectively, leading to a permanent high level of activation (indicated by the blue text and contour line for minimum force greater than 1% of maximum force at the top of the plot). When the affinity of TnI for actin is much lower than for TnC, muscle activation is decreased, (indicated by the green text and contour line for maximum force less than 95% of overall maximum force in the bottom right of the plot).](pcbi.1004376.g003){#pcbi.1004376.g003}

Finally, we set *K* ~DC~ = 5.9 μM based on a half-activation value Ca~50~ for the force-calcium relationship of approximately 0.5 μM, consistent with requirements of dynamic models with a peak calcium at the lower end of typical physiological range of 0.5--1 μM \[[@pcbi.1004376.ref054]\], experimental data (*K* ~DC~ ≈ 5 μM \[[@pcbi.1004376.ref019]\]) and previous estimates of *K* ~DC~ in models (between 1 and 10 μM \[[@pcbi.1004376.ref037]\]).

### Testing approximation strategies for the crossbridge model {#sec009}

We proposed two strategies in making the model with crossbridges computationally tractable: the reduction of tropomyosin states to representative states based on connected stretches of unblocked RU's (c.f. [Fig 2](#pcbi.1004376.g002){ref-type="fig"}), and subsequent Monte Carlo sampling with *n* ~*s*~ = 1000 samples per representative state. In addition we proposed a much simpler 'independent crossbridge model' which does not include XB-RU and XB-XB interactions. To determine which of these models provides the best compromise between accuracy and computational tractability, we compare them with a brute-force approach on a smaller filament with 7 RU's and 18 crossbridges. With 2^7+18^ ≈ 33 million calculations for the tropomyosin bending energy *E* ~tm~, this is the largest thin filament for which a brute force approach is currently tractable.

Results in [Fig 4](#pcbi.1004376.g004){ref-type="fig"} are based on four simulation results: explicit calculation of all tropomyosin states, approximation of RU (un)blocking with representative states with with brute-force crossbridge calculation, further approximation of crossbridge states with Monte Carlo sampling, and an approximation with the assumption of independent crossbridge binding. Results indicate the 'representative state' approximation (with brute-force crossbridge calculation) overlaps completely with an exhaustive brute-force approach, i.e. the state of tropomyosin and its ability to bind myosin is very well approximated from the number and length of connected unblocked regulatory units. The Monte Carlo sampling approach also approximates the accurate solution well, with maximal differences of 1.2% at higher force levels. By contrast, the independent crossbridge approximation is shown to significantly underestimate force, by approximately 50%.

![Comparison of different modelling approaches on a short filament.\
Shown are the different modelling approaches on a filament with *n* = 7 RU's and 18 crossbridges. Results for the 'brute-force solution', and the representative state approximation overlap and are indicated by a single line. The Monte Carlo approximation performs well, with only a ∼ 1% difference at higher force levels. Comparison with the independent crossbridge approximation shows the importance of including XB-RU cooperativity, which increases calcium sensitivity, as well as XB-XB cooperativity, which increases maximum force development as shown by the difference in the number of crossbridges bound per half-sarcomere at high calcium. Note that due to the lower number of 7 RU's, cooperativity is significantly lower than in realistic models with 26 RU's presented in other results, and the duty ratio is moderately reduced to approximately 2.7/18 = 15%.](pcbi.1004376.g004){#pcbi.1004376.g004}

The difference between the independent crossbridge model and the brute force model can be attributed to two effects. Firstly, as nearly all RU's are unblocked at the maximum calcium concentration, the difference in maximal force can be attributed to XB-XB cooperativity, i.e. the shifting of tropomyosin by a crossbridge makes it easier for neighbouring crossbridges to attach. Secondly, there is a significant shift towards lower calcium sensitivity in the independent crossbridge approximation, which is attributed to XB-RU cooperativity, where bound crossbridges inhibit the transition of tropomyosin to the 'blocked' state, effects which are also absent in the independent crossbridge approximation. All subsequent steady state analysis is performed on the representative state model with Monte Carlo crossbridge sampling.

### Results for cooperative activation and XB-RU interactions {#sec010}

In this section we use the full 26 RU model developed using the Monte Carlo sampling approach to investigate steady-state cooperative activation and XB-RU effects. [Fig 5A](#pcbi.1004376.g005){ref-type="fig"} shows cooperative tension development of the model, and compares it with the unblocking of RU's. These results show that, at lower Ca^2+^, unblocking of RU's is significantly less cooperative than force. This difference in cooperativity results in around 5% of RU's still being unblocked at points when force is reduced to nearly zero. This difference is explained by low affinity of myosin for 1--4 neighbouring unblocked RU's, shown in [Fig 5B](#pcbi.1004376.g005){ref-type="fig"}, and is a significant difference compared to most previous models which use a linearly increasing probability for crossbridge binding as a function of RU unblocking.

![Steady-state model behaviour.\
Panel A shows the Force-calcium curve relationship of the model alongside the RU unblocking as a function of calcium. This shows that the Force-calcium curve is significantly steeper than RU activation, as indicated by its higher Hill coefficient. Panel B explains this effect by looking at the expected number of crossbridges on tropomyosin chains with a fixed number of RU's unblocked, where the probability of tropomyosin sub-states is according to Boltzmann's law. This shows that there is significant inhibition of crossbridge binding at low numbers of RU's, compared to models where the 'linear response' of crossbridge binding being proportional to the number of unblocked RU's is assumed.](pcbi.1004376.g005){#pcbi.1004376.g005}

[Fig 4](#pcbi.1004376.g004){ref-type="fig"} presented in the previous section showed a shift in calcium sensitivity when XB-XB and XB-RU effects are included. There are several experiments which reveal important interactions between crossbridge affinity and calcium sensitivity, including experiments in which blebbistatin or sodium vanadate are used to decrease crossbridge affinity \[[@pcbi.1004376.ref001], [@pcbi.1004376.ref016], [@pcbi.1004376.ref054]\], and experiments with dATP where crossbridge affinity is increased \[[@pcbi.1004376.ref055]\]. In general, an increase in crossbridge affinity leads to higher calcium sensitivity, i.e. the muscle activating at lower \[Ca^2+^\]. We modeled the effects of blebbistatin by a decrease in myosin head affinity (3× higher *K* ~DM~), and the effects of dATP by a higher myosin affinity (25% decrease in *K* ~DM~). Increasing myosin affinity led to a leftward shift of the force-calcium curve, and vice versa, as shown in [Fig 6A](#pcbi.1004376.g006){ref-type="fig"}. The changes in crossbridge affinity for blebbistatin and dATP were fitted to the change in maximum force shown in experiments, resulting in good quantitative agreement for the predicted shift in calcium sensitivity (ΔpCa~50~), as indicated in [Table 2](#pcbi.1004376.t002){ref-type="table"}. In addition, we tested the models ability to activate due a high myosin affinity as observed in conditions of low ATP \[[@pcbi.1004376.ref045]\]. For the rigor test we vary *K* ~DM~ and record the force generation predicted by the model. The results ([Fig 6B](#pcbi.1004376.g006){ref-type="fig"}) are qualitatively similar to experimental data for pCa 4.5 \[[@pcbi.1004376.ref045]\], with a decreasing sigmoidal relationship. Thus, our model replicates the shifts in calcium sensitivity shown in experiments where crossbridge affinity is modified, and is also able to activate in the absence of significant calcium due to rigor crossbridges.

![Force-pCa curves for the model and effects of changes in crossbridge affinity.\
The effects of crossbridge inhibition by substances such as blebbistatin and sodium vanadate were simulated by changing the dissociation constant of myosin (*K* ~DM~), resulting in significant changes in calcium sensitivity. Shown in panel (A) are the default, highly cooperative, force-calcium relationship of the model (*n* ~*H*~ = 5.1), along with the following virtual experiments: (1) In red: a factor 3 decrease in crossbridge affinity. This reproduces data from experiments with the cross-bridge inhibitor blebbistatin \[[@pcbi.1004376.ref016]\], showing a decrease in calcium sensitivity and a mild decrease in cooperativity (*n* ~*H*~ = 4.2). (2) In blue: a 33% increase in crossbridge affinity. This reproduces data from experiments with the cross-bridge augmenter 2-deoxy-ATP (dATP) \[[@pcbi.1004376.ref055]\], showing an increase in calcium sensitivity and a small increase in cooperativity (*n* ~*H*~ = 5.3). Panel (B) shows the *K* ~DM~-dependence of force at zero Ca^2+^ and at pCa 4.5 \[[@pcbi.1004376.ref045]\], showing the model produces maximal force at both calcium levels for a sufficiently high myosin affinity, and a sigmoidal relationship between *K* ~DM~ and force. The dashed line indicates the value of *K* ~DM~ used in the model, which intersects the pCa 4.5 curve at approximately 0.25, the duty ratio of myosin used in the model. Thus, the maximal force generated in panel (B) for *K* ~DM~ → 0 is approximately 4× higher than the 'default model' curve in panel (A).](pcbi.1004376.g006){#pcbi.1004376.g006}

10.1371/journal.pcbi.1004376.t002

###### Model predictions for shifts in pCa~50~ after changes in crossbridge affinity.

![](pcbi.1004376.t002){#pcbi.1004376.t002g}

  Experiment                  Blebbistatin \[[@pcbi.1004376.ref016]\]            dATP \[[@pcbi.1004376.ref055]\]
  --------------------------- -------------------------------------------------- ---------------------------------
  Experimental conditions     Skinned, room temp., mouse                         Skinned, 15°C, rat
  Experimental max. force     -73%                                               +31%
  Experimental ΔpCa~50~       -0.34                                              +0.13
  Model conditions            Both models represent intact muscle, 37°C, mouse   
  Model max. force (fitted)   -73.6%                                             +30.9%
  Model ΔpCa~50~              -0.248                                             +0.117

Predictions of changes in calcium sensitivity (ΔpCa~50~) when fitting crossbridge affinity in the model to reproduce changes in maximal force shown experimentally. See [Fig 6](#pcbi.1004376.g006){ref-type="fig"} for full traces.

Dynamic model {#sec011}
-------------

The steady-state models developed in the previous sections replicate a range of experimental measurements related to steady-state cooperative activation. Cooperative effects also have an important impact on beat-to-beat dynamic tension generation, and may have different roles in different species due to differences in heart rate and calcium dynamics. To be able to investigate the role of cooperative activation in dynamic tension generation, in this section we extend our proposed framework to simulate dynamic changes in tension in response to transient changes in Ca^2+^.

For a dynamic model of *n* RU's and *m* crossbridges, we use: A regular grid of (*n* + 1) ⋅ (*m* + 1) state variables TmXB~i,j~ which represent the fraction of half sarcomeres with *i* RU's unblocked and *j* crossbridges bound.The state variable TnC~B~, the fraction of RU's that are blocked with Ca^2+^ bound to TnCThe state variable TnC~U~, the fraction of RU's that are unblocked with Ca^2+^ (but not TnI) bound to TnC,The state variable TnITnC, the fraction of RU's that are blocked with Ca^2+^ and TnI bound to TnC

Several dependent variables are useful in the formulation in the differential equations for these states: $$U = \sum\limits_{i = 0}^{n}\sum\limits_{j = 0}^{m}\mspace{180mu}\frac{i}{n}\text{TmXB}_{\text{i,j}}\mspace{126mu}(\text{fraction\ of\ RU}’\text{s\ in\ the\ unblocked\ state})$$ $$B = 1 - U\mspace{252mu}(\text{fraction\ of\ RU}’\text{s\ in\ the\ blocked\ state})$$ $$\text{TnI} = U - \text{TnITnC}\mspace{72mu}(\text{fraction\ of\ RU}’\text{s\ with\ TnI\ not\ bound\ to\ actin\ or\ TnC})$$

The kinetics of the tropomyosin and crossbridge states can now be defined using a standard Markov model approach, with transition rates defined using the ratio of Boltzmann terms: $$\begin{aligned}
\frac{\text{dTmXB}_{\text{i,j}}}{\text{d}t} & {= \overset{\text{if}{i \neq 0}}{\overset{︷}{k_{\text{i-1,j}}^{\text{b}\rightarrow\text{u}}\text{TmXB}_{\text{i-1,j}}-k_{\text{i,j}}^{\text{u}{\mspace{-200mu}\rightarrow\mspace{-270mu}}\text{b}}\text{TmXB}_{\text{i,j}}}}\mspace{600mu} + \mspace{600mu}\overset{\text{if}{i \neq n}}{\overset{︷}{k_{\text{i+1,j}}^{\text{u}{\mspace{-200mu}\rightarrow\mspace{-270mu}}\text{b}}\text{TmXB}_{\text{i+1,j}}-k_{\text{i,j}}^{\text{b}\rightarrow\text{u}}\text{TmXB}_{\text{i,j}}}}} \\
 & {+ \mspace{180mu}\overset{\text{if}{j \neq 0}}{\overset{︷}{k_{\text{i,j-1}}^{\text{xb+}}\text{TmXB}_{\text{i,j-1}}-k_{\text{i,j}}^{\text{xb-}}\text{TmXB}_{\text{i,j}}}}\mspace{600mu} + \mspace{600mu}\overset{\text{if}{j \neq m}}{\overset{︷}{k_{\text{i,j+1}}^{\text{xb-}}\text{TmXB}_{\text{i,j+1}}-k_{\text{i,j}}^{\text{xb+}}\text{TmXB}_{\text{i,j}}}}} \\
\end{aligned}$$ $$k_{\text{i,j}}^{\text{xb+}}/k_{\text{i,j+1}}^{\text{xb-}} = \frac{\text{SE}_{i,j + 1}}{\text{SE}_{i,j}} \cdot \frac{1}{K_{\text{DM}}}$$ $$k_{\text{i,j}}^{\text{b}\rightarrow\text{u}}/k_{\text{i+1,j}}^{\text{u}\rightarrow\text{b}} = \frac{\text{SE}_{i + 1,j}}{\text{SE}_{i,j}} \cdot K_{\text{DA}} \cdot \frac{U}{\text{TnI}}$$

Where $k_{\text{i,j}}^{\text{u}\rightarrow\text{b}},k_{\text{i,j}}^{\text{b}\rightarrow\text{u}}$ are the transition rates from unblocked to blocked, and blocked to unblocked, respectively, for the state with *i* unblocked RU's. The ratio $k_{\text{i,j}}^{\text{b}\rightarrow\text{u}}/k_{\text{i+1,j}}^{\text{u}\rightarrow\text{b}}$ follows naturally from the ratio of Boltzmann terms between the states, requiring only the addition of the probability that TnI is free to bind to actin *P*(TnI free∣RU unblocked) = TnI/*U*.

The framework based on Boltzmann's law can be used to determine the ratios of transition rates, but does not result in an absolute on- and off-rate. To determine how the energy difference influences the on- and off-rates, we use a similar approach as proposed by Campbell et al. \[[@pcbi.1004376.ref038]\]. Given two states *S* ~1~, *S* ~2~ with energies *E* ~1~, *E* ~2~ and Boltzmann terms $B_{1} = e^{- \frac{E_{1}}{k_{\! B}T}},B_{2} = e^{- \frac{E_{2}}{k_{\! B}T}}$, the transition rates between them are given by: $$k_{S_{1}\rightarrow S_{2}}/k_{S_{2}\rightarrow S_{1}} = B_{2}/B_{1}$$ $$k_{S_{1}\rightarrow S_{2}} = \left( B_{2}/B_{1} \right)^{r}$$ $$k_{S_{2}\rightarrow S_{1}} = \left( B_{2}/B_{1} \right)^{- (1 - r)}$$ The parameter *r* represents how strongly the on-rate and off-rate depend on the difference in energy between the states *E* ~1~/*E* ~2~, ranging from only the off-rate (*r* = 0) to only the on-rate (*r* = 1). We apply this model to the effect of tropomyosin deformation on the rates $k_{\text{i,j}}^{\text{b}\rightarrow\text{u}},k_{\text{i,j}}^{\text{u}\rightarrow\text{b}}$, with the assumption that both rates are equally affected (*r* = 0.5). The affinity of TnI for actin *K* ~DA~ is handled separately from the influence of tropomyosin, and is split into two rate constants *k* ~A-~, *k* ~A+~. In addition, as our state TmXB~i,j~ is a combination of several sub-states with *i* unblocked RU's, what remains is to take into account is the difference in the number of transitions in the average state, given by (*n* − *i*) potential RU's for a blocked-to-unblocked transition and *i* potential RU's for a unblocked-to-blocked transition, for any state. Combined, these considerations result in the following equations for the rate constants: $$k_{\text{i,j}}^{\text{b}\rightarrow\text{u}} = k_{\text{A-}} \cdot \left( {n - i} \right) \cdot \left( \frac{\text{SE}_{i + 1,j}/\left( {n - i} \right)}{\text{SE}_{i,j}/\left( {i + 1} \right)} \right)^{r}$$ $$k_{\text{i,j}}^{\text{u}\rightarrow\text{b}} = k_{\text{A+}} \cdot i \cdot \left( \frac{\text{SE}_{i,j}/\left( {n - i + 1} \right)}{\text{SE}_{i - 1,j}/i} \right)^{- {({1 - r})}} \cdot \frac{\text{TnI}}{U}$$ For crossbridge binding and unbinding rates we introduce a parameter *q* to represent the effect of tropomyosin deformation on the crossbridge binding and unbinding rate. Two different choices will be considered. Firstly, tropomyosin deformation affecting both on- and off-rate equally (*q* = 0.5) similar to RU unblocking. Secondly, the choice of a constant unbinding rate (*q* = 1), where only the on-rate is affected by tropomyosin deformation. The impact of this choice will be considered in the next section. Taking into account the number of potential crossbridges, equations for the transition rates are given by: $$k_{\text{i,j}}^{\text{xb+}} = k_{\text{M+}} \cdot \left( {m - j} \right) \cdot \left( \frac{\text{SE}_{i,j + 1}/\left( {m - j} \right)}{\text{SE}_{i,j}/\left( {j + 1} \right)} \right)^{q}$$ $$k_{\text{i,j}}^{\text{xb-}} = k_{\text{M-}} \cdot j \cdot \left( \frac{\text{SE}_{i,j}/\left( {m - j + 1} \right)}{\text{SE}_{j - 1,i}/j} \right)^{- {({1 - q})}}$$ We model TnC and TnI kinetics using simplified global state variables for the blocked and unblocked regulatory units. The equations for these kinetics are given by: $$\frac{\text{dTnITnC}}{\text{d}t} = k_{\text{I+}}\text{TnC}_{\text{U}} - k_{\text{I-}}\text{TnITnC}$$ $$\frac{\text{dTnC}_{\text{B}}}{\text{d}t} = k_{\text{C+}}\left\lbrack \text{Ca}^{2 +} \right\rbrack\mspace{180mu}\left( B - \text{TnC}_{\text{B}} \right) - k_{\text{C-}}\text{TnC}_{\text{B}} - J_{\text{bu}} + J_{\text{ub}}$$ $$\frac{\text{dTnC}_{\text{U}}}{\text{d}t} = k_{\text{C+}}\left\lbrack \text{Ca}^{2 +} \right\rbrack\mspace{180mu}\left( U - \text{TnC}_{\text{U}} - \text{TnITnC} \right) - k_{\text{C-}}\text{TnC}_{\text{U}} - \frac{\text{dTnITnC}}{\text{d}t} + J_{\text{bu}} - J_{\text{ub}}$$

These equations represent standard Michaelis-Menten kinetics, apart from the terms *J* ~bu~ and *J* ~ub~, which represent the 'flux' of Ca^2+^ bound to TnC between the global TnC buffers for blocked and unblocked RU's. Each time an RU blocks or unblocks, we need to consider the probability of a calcium ion moving between TnC~B~ and TnC~U~. These fluxes are given by: $$J_{\text{bu}} = \frac{1}{n}\frac{\text{TnC}_{\text{B}}}{B}\sum\limits_{j = 0}^{m}\sum\limits_{i = 0}^{n - 1}k_{\text{i,j}}^{\text{b}\rightarrow\text{u}}\text{TmXB}_{\text{i,j}}$$ $$J_{\text{ub}} = \frac{1}{n}\frac{\text{TnC}_{\text{U}}}{\text{TnI}}\sum\limits_{j = 0}^{m}\sum\limits_{i = 1}^{n}k_{\text{i,j}}^{\text{u}\rightarrow\text{b}}\text{TmXB}_{\text{i,j}}$$ The sums represent the total transition rates between blocked and unblocked RU's from [Eq 13](#pcbi.1004376.e023){ref-type="disp-formula"}, multiplied by 1/*n* to represent one RU out of *n* changing for each transition. This is multiplied by the probability of a calcium ion being present on a closing or blocking RU, which is $\frac{\text{TnC}_{\text{B}}}{B}$ for blocked units and $\frac{\text{TnC}_{\text{U}}}{\text{TnI}}$ for unblocked units as this latter probability needs to be considered over unblocked RU's which do not have TnI bound to TnC·Ca^2+^ (and $k_{\text{i,j}}^{\text{u}\rightarrow\text{b}}$ already contains a factor TnI/*U*).

For a numerical implementation, quantities such as TnI/*U* (in [Eq 20](#pcbi.1004376.e034){ref-type="disp-formula"}) should be calculated as TnI/ max(*U*, *ɛ*) for a small constant *ɛ* to avoid undefined 0/0 quantities. Secondly, many of the states TmXB~i,j~ are not populated in practice and can be removed from the formulation for improved efficiency and stability. Specifically, for every *i* we include states TmXB~i,j~ up to the value of *j* for which SE~*i*,*j*~ \< 10^−6^∑~*i*~ SE~*i*,*j*~, which reduces the number of states from 1794 to ∼ 750 without significantly affecting the solution. The model's initial condition should be determined by pacing for a specific calcium transient and parametrization, starting from the completely de-activated state (all state variables set to 0 except TmXB~0,0~ = 1).

Results {#sec012}
=======

Tension development in different species {#sec013}
----------------------------------------

In a physiological setting, cooperativity is a key component of normal activation and relaxation of the heart. Thus, an important test of the effectiveness of a model in reproducing physiological cooperativity is its ability to reproduce realistic tension based on experimentally measured calcium transients. In this section, we investigate if tension development across different species is consistent with identical cooperativity, despite significant differences in heart rate. Although the equations in the dynamic model appear to have a large number of free parameters, all transition rates follow from the 9 parameters listed in [Table 3](#pcbi.1004376.t003){ref-type="table"}. We first parametrize our model to reproduce twitch tension at 37°C in mouse, as we have found this the most challenging test case in practice, and start with the *q* = 0.5 case. Active tension was calculated by setting the maximal tension developed to 120 kPa as in previous work \[[@pcbi.1004376.ref036]\] (see [S1 Text](#pcbi.1004376.s001){ref-type="supplementary-material"} for details). We vary *k* ~A-~, *k* ~I-~, *k* ~M-~, *k* ~C+~ between 0.01/ms and 100/ms, and found that tension development and relaxation are all sensitive to the choice of these parameters. As *k* ~A-~, *k* ~I-~ are generally not thought to be rate-limiting \[[@pcbi.1004376.ref056]\], we set both of them to 10/ms. Parameters *k* ~M-~, *k* ~C+~ are then determined by requirements for time to peak tension and relaxation times according to the ranges of experimental measurements determined in previous work \[[@pcbi.1004376.ref036]\] and summarized in [Table 4](#pcbi.1004376.t004){ref-type="table"}, resulting in *k* ~C+~ = *k* ~M-~ = 0.5/ms.

10.1371/journal.pcbi.1004376.t003

###### List of model parameters.

![](pcbi.1004376.t003){#pcbi.1004376.t003g}

  Parameter   Value        Determined based on
  ----------- ------------ ---------------------------------------------------------------------------------------------------------
  *K* ~DA~    10^−3^       Cooperativity, effective competitive binding
  *K* ~DI~    4 ⋅ 10^−3^   Cooperativity, effective competitive binding
  *K* ~DC~    5.9 μM       Required Ca~50~, literature data, previous modelling \[[@pcbi.1004376.ref019], [@pcbi.1004376.ref037]\]
  *K* ~DM~    2            Duty ratio in fully activated muscle
  *k* ~A-~    10/ms        High enough to not be rate-limiting
  *k* ~I-~    10/ms        High enough to not be rate-limiting
  *k* ~C+~    0.5/ms       Fitting to TPT, RT, Maximum twitch force
  *k* ~M-~    0.5/ms       Fitting to TPT, RT, Maximum twitch force
  *γ*         2            Cooperativity *n* ~2~ ≈ 8

Shown here are all model parameters for the mouse model with the method used for parametrization.

10.1371/journal.pcbi.1004376.t004

###### Parametrization of the model for different species.

![](pcbi.1004376.t004){#pcbi.1004376.t004g}

                                                                 Mouse                                                Rat                          Human
  -------------------------------------------------------------- ---------------------------------------------------- ---------------------------- ----------------------------------------------------
  Parametrization with variable XB unbinding rates (*q* = 0.5)                                                                                     
  *K* ~DI~                                                       4 ⋅ 10^−3^                                           8.7 ⋅ 10^−3^                 2.6 ⋅ 10^−3^
  *k* ~M-~ (ms^−1^)                                              0.5                                                  0.09                         0.016
  Model TPT (ms)                                                 33                                                   41                           170
  Model RT~50~ (ms)                                              29                                                   30                           124
  Model RT~90~ (ms)                                              49                                                   51                           217
  Model RT~95~ (ms)                                              56                                                   57                           247
  Parametrization with constant XB unbinding rates (*q* = 1)                                                                                       
  *K* ~DI~                                                       4 ⋅ 10^−3^                                           8.7 ⋅ 10^−3^                 3.3 ⋅ 10^−3^
  *k* ~M-~ (ms^−1^)                                              2                                                    0.3                          0.049
  Model TPT (ms)                                                 33                                                   37                           157
  Model RT~50~ (ms)                                              28                                                   29                           125
  Model RT~90~ (ms)                                              49                                                   55                           250
  Model RT~95~ (ms)                                              57                                                   64                           298
  Experimental data                                                                                                                                
  TPT (ms)                                                       30--40                                               39 ± 6                       147--172
  RT~50~ (ms)                                                    22--35                                               31 ± 5                       109--125
  RT~90~ (ms)                                                    45--60                                                                            
  RT~95~ (ms)                                                                                                                                      334 ± 43
  References                                                     \[[@pcbi.1004376.ref057]--[@pcbi.1004376.ref061]\]   \[[@pcbi.1004376.ref062]\]   \[[@pcbi.1004376.ref063], [@pcbi.1004376.ref064]\]

Shown in this table are the changes to *K* ~DI~, *k* ~M-~ needed to replicate mouse, rat and human dynamic muscle function. Also shown are metrics of time to peak tension (TPT) and 50%/90%/95% relaxation times (RT~50~, RT~90~, RT~95~) of the dynamic model. Results for force transient are also shown in [Fig 7](#pcbi.1004376.g007){ref-type="fig"}. The bottom section shows experimental data ranges used for parametrization of the models. Two sets of results are shown, corresponding to a constant crossbridge unbinding rate (*q* = 1 in [Eq 22](#pcbi.1004376.e036){ref-type="disp-formula"}) and a crossbridge unbinding rate that is variable, modified by the average difference in free energy for the tropomyosin deformation (*q* = 0.5 in [Eq 22](#pcbi.1004376.e036){ref-type="disp-formula"})

Based on this initial parametrization, we investigate dynamic function of the model across species, applying it to both tension generation in isometric twitches and tension redevelopment. Firstly, we parameterize the model for three different species as well as both choices of crossbridge binding rates (*q* = 1 or *q* = 0.5). Troponin C is a highly conserved protein and its kinetics are not sensitive to temperature or (mammalian) species \[[@pcbi.1004376.ref056]\] while crossbridge cycling rates are highly dependent on species, temperature and myosin heavy chain isoform. However, despite the lack of variation in TnC properties, calcium transients vary between the different species while a similar peak isometric tension of approximately 40 kPa is required to be consistent with whole organ contraction across different species \[[@pcbi.1004376.ref036], [@pcbi.1004376.ref070]\]. Inspired by these observations we use our model to test the hypothesis that contraction across different mammalian species (mouse, rat and human) can be reproduced using only differences in TnI-TnC affinity as given by *K* ~DI~ and the rate of crossbridge kinetics as given by *k* ~M-~. We determined these two parameters by using a two-dimensional parameter sweep which shows the influence of different constraints on both parameters. Constraints used include species-specific constraints for time to peak tension and relaxation times, and constraints on minimum force \< 1 kPa and maximum force between 35--45 kPa in all species. [Table 4](#pcbi.1004376.t004){ref-type="table"} shows the parametrization for the three different species \[[@pcbi.1004376.ref057]--[@pcbi.1004376.ref064]\], and [Fig 7](#pcbi.1004376.g007){ref-type="fig"} shows the corresponding force transients. We were able to capture the different twitch kinetics between species with variations in only *K* ~DI~ and *k* ~M-~. The resulting parameters show that changes in crossbridge kinetics are consistent with differences in heart rate (mouse \> rat \> human), while changes to the parameter *K* ~DI~ for TnI affinity for TnC⋅Ca^2+^ correspond to differences in the calcium transients (c.f. [Fig 7B](#pcbi.1004376.g007){ref-type="fig"}). With respect to the choice of crossbridge unbinding rate on the thin filament state (*q* = 1 or *q* = 0.5), the parameter *K* ~DI~ could be kept the same for the different choices although it was not fixed a priori. However, *k* ~M-~ needs to be significantly higher when constant unbinding rates are used. Overall, model results suggest contractile function across these different species are consistent with a common mechanism and kinetics for thin filament cooperative activation.

![Isometric tension in mouse, rat, and human.\
Shown in panel A are the tension transients described in [Table 4](#pcbi.1004376.t004){ref-type="table"}. The last 500 ms of the human tension transients are not shown, as force is very low throughout. These results were driven by fixed calcium transients shown in panel B, based on recent data in mouse \[[@pcbi.1004376.ref066]\] and rat (see [S1 Text](#pcbi.1004376.s001){ref-type="supplementary-material"} \[[@pcbi.1004376.ref067], [@pcbi.1004376.ref068]\]) measured at 37°C and 6 Hz, and data from Coppini et al. \[[@pcbi.1004376.ref069]\] (see [S1 Text](#pcbi.1004376.s001){ref-type="supplementary-material"} \[[@pcbi.1004376.ref070]\]). Two sets of results are shown, corresponding to a constant crossbridge unbinding rate (*q* = 1 in [Eq 22](#pcbi.1004376.e036){ref-type="disp-formula"}) and a crossbridge unbinding rate that is variable, modified by the average difference in free energy for the tropomyosin deformation (*q* = 0.5 in [Eq 22](#pcbi.1004376.e036){ref-type="disp-formula"})](pcbi.1004376.g007){#pcbi.1004376.g007}

Simulating tension redevelopment {#sec014}
--------------------------------

The rate of tension redevelopment (*k* ~tr~) has been shown to vary significantly depending on \[Ca^2+^\], with experiments showing a range of 4--8× between the lowest and highest rates observed \[[@pcbi.1004376.ref055], [@pcbi.1004376.ref071]--[@pcbi.1004376.ref074]\]. These differences in tension redevelopment rates have been linked to the effect of thin filament activation kinetics \[[@pcbi.1004376.ref055]\], and is often not well reproduced by computational models. To test our model's ability to replicate and explain these more complex dynamic effects which involve multiple contractile proteins, and isolate the role of cooperativity in tension redevelopment, we have performed simulations of the calcium dependence of the redevelopment rate of tension.

As our model does not account for dynamic length changes, we apply the following procedure to simulate the *k* ~tr~ protocol. For each calcium level the model is run to steady-state, and 50% of crossbridges are instantly detached to simulate the state of the filament after a rapid shorten-relengthen protocol (c.f. trace in \[[@pcbi.1004376.ref071]\]). Crossbridges are detached by setting TmXB~i,j/2~ to TmXB~i,j~ for even *j*, and TmXB~i,⌈j/2⌉~ and TmXB~i,⌊j/2⌋~ to $\frac{\text{TmXB}_{\text{i,j}}}{2}$ for odd *j*, where ⌈⋅⌉, ⌊⋅⌋ denote rounding up and down, respectively. Subsequently the model is run with normal binding rates to simulate tension redevelopment, and mono-exponential curve is fitted to the resulting force: $$\frac{F}{F_{\text{max}}} = 1 - ae^{- k_{\text{tr}}t}$$ To investigate the importance of XB-RU cooperative activation, this procedure is repeated with disabled TnI-actin dynamics (*k* ~A-~ = *k* ~A+~ = 0) during tension redevelopment to prevent RU (un)blocking. Results in [Fig 8](#pcbi.1004376.g008){ref-type="fig"} show that a constant unbinding rate independent of tropomyosin deformation (*q* = 1) shows a larger difference between the minimal *k* ~tr~ and the value at high Ca^2+^, with around a 6.5× difference in mouse and rat. Interestingly, unlike experimental data, our results show a clear minimum rather than the typical monotonically increasing *k* ~tr~. We attribute this effect to the smaller range of \[Ca^2+^\] investigated in experiments and the tendency for noise to dominate measurements at low force levels. For a variable unbinding rate, this ratio is only around 3--4 in mouse and rat. In addition, *k* ~tr~ rises steeply at low \[Ca^2+^\], as the unbinding rate becomes very fast. In both cases, transient blocking of RU's is the main mechanism leading to the minimum in *k* ~tr~. At very low and high \[Ca^2+^\], the state of the thin filament is not much affected by crossbridges, being mostly blocked or unblocked regardless. This can be seen by the similar *k* ~tr~ for results with and without RU (un)blocking disabled. However, at intermediate \[Ca^2+^\], unblocking of crossbridges will cause RU's to start moving to the 'blocked' position, requiring additional time to be re-activated by XB-RU interaction resulting in the lower *k* ~tr~.

![Tension redevelopment rates in mouse, rat and human.\
Panel A shows results for *q* = 1 (constant crossbridge unbinding rate), and Panel B plot results for *q* = 0.5 (variable crossbridge unbinding rate). Each plot shows both normal kinetics, and with RU (un)blocking due to TnI-actin binding disabled for mouse, rat, and human parametrizations of the model. The difference in *k* ~tr~ when disabling RU (un)blocking reveals the influence of transient blocking of RU's on tension redevelopment kinetics. This transient blocking causes a significant difference between the minimum *k* ~tr~ and the *k* ~tr~ at high \[Ca^2+^\], the ratio of which is indicated with the results. Note the shift in the minimum *k* ~tr~ is not significant here as it is a direct result of the difference in peak calcium shown in [Fig 7A](#pcbi.1004376.g007){ref-type="fig"} and subsequent parametrization.](pcbi.1004376.g008){#pcbi.1004376.g008}

Discussion {#sec015}
==========

In this paper we developed a novel model of cardiac contraction and used it to investigate the effects of different cooperative mechanisms on both the steady-state and dynamic behaviour of cardiac muscle. We have been able to show the relative importance of the different potential mechanisms for generating the steeply cooperative force-calcium relationship in muscle. Firstly, RU-RU cooperativity is the clear dominant mechanism for cooperativity near and above Ca~50~, caused by progressively easier unblocking for any individual RU as more of their close neighbours are unblocked. This can be seen from the cooperativity of the RU unblocking in the absence of crossbridges ([Fig 4](#pcbi.1004376.g004){ref-type="fig"}) and RU activation ([Fig 5](#pcbi.1004376.g005){ref-type="fig"}), with moderate cooperativity which is not significantly biphasic. Although RU-RU cooperativity deriving from end-to-end interactions between tropomyosin is a major part of cooperative activation, it does not explain all of the cooperative activation seen.

Secondly, cooperativity between crossbridges, i.e. of the 'XB-XB' type, approximately doubles the number of crossbridges bound at maximal activation ([Fig 4](#pcbi.1004376.g004){ref-type="fig"}), and in doing so increases the effects of XB-RU interactions, but does not in itself significantly increase the steepness of the force-calcium relationship. Thirdly, we have shown that activation of regulatory units by crossbridges, i.e. 'XB-RU' cooperativity, is important for determining calcium sensitivity. Our model reproduces the effects of changes in the affinity of myosin XB on calcium sensitivity shown experimentally \[[@pcbi.1004376.ref001], [@pcbi.1004376.ref016], [@pcbi.1004376.ref054]\] ([Fig 6](#pcbi.1004376.g006){ref-type="fig"}). Despite differences in experimental conditions due to temperature and permeabilization of muscle, there is good quantitative agreement in the shift in calcium sensitivity (ΔpCa~50~, [Table 2](#pcbi.1004376.t002){ref-type="table"}). Our results show a mild decrease in cooperativity for simulations with decreased crossbridge affinity (*n* ~*H*~ = 5.1 → 4.2) which is almost entirely in the lower half (below Ca~50~) of the force-calcium relationship (*n* ~2~ = 7.5 → 6.1), with nearly identical values for the upper half (*n* ~1~ = 2.7 → 2.4). The modest change in Hill coefficient may explain some of the contradictory results in the literature so far, as Hill fits are inherently sensitive to the choice of fitting method and calcium concentrations used, in addition to experimental noise. The ratio of *n* ~2~/*n* ~1~ varies between 2.5--2.7 in these simulations, compared to approximately 2.0 in experiments \[[@pcbi.1004376.ref052], [@pcbi.1004376.ref053]\]. This ratio is very sensitive to the choice of fitting methods and \[Ca^2+^\] window used for fitting, which could explain these differences. Alternatively, these differences could be the result of the simplified sarcomere geometry in the current model. In addition, our model reproduces the activation of muscle by high-affinity 'rigor' crossbridges, such as in conditions of low ATP \[[@pcbi.1004376.ref045]\] or special experimental preparations with NEM-S1 myosin \[[@pcbi.1004376.ref075]\], even in the complete absence of Ca^2+^. Although results from these experiments are less important for contraction models to reproduce as they represent conditions far from physiological, they are a direct result of XB-RU cooperativity in our model and thus increases confidence in the choice of biophysical framework. In addition this is a novel feature for computational models, as in the vast majority of models in this area RU's are mathematically unable to activate in the absence of Ca^2+^ \[[@pcbi.1004376.ref034]--[@pcbi.1004376.ref036], [@pcbi.1004376.ref038]\]. Our model also explains effects of force on Ca^2+^-TnC binding. The proposed effects of crossbridges on Ca^2+^-TnC affinity appear in our model as an emergent property of the competitive binding of TnI and pinning of tropomyosin by both TnI and crossbridges. Specifically, crossbridge binding prevents tropomyosin to move to the blocked position, and this effectively prevents TnI from binding to tropomyosin-actin. This in turn increases the relative time TnI is bound to TnC⋅Ca^2+^, thus increasing the effective affinity of TnC for Ca^2+^, as Ca^2+^ is highly unlikely to dissociate from TnC when TnI is bound. As a result, our framework represents the observed tension-dependent feedback on TnC affinity without hysteresis in the steady-state force-pCa relationship.

Finally, we have identified an important effect by which approximately five neighbouring unblocked tropomyosin units are required to significantly bind myosin ([Fig 5](#pcbi.1004376.g005){ref-type="fig"}). These effects are also partly responsible for the very steep force-calcium relationship in the region below Ca~50~ and the strongly biphasic shape of the force-calcium curve. Specifically, it allows for near-zero force in a state with a significant fraction of unblocked RU's with Ca^2+^ bound to TnC, as there are not enough consecutive unblocked RU's to generate significant force. Achieving low force at physiological diastolic \[Ca^2+^\] of approximately 0.1--0.2 μM is particularly important for effective diastolic filling in the heart. In addition this result shows that the common modelling assumptions of crossbridge binding properties being linearly proportional to the number of unblocked RU's is potentially inaccurate.

Thus, cooperative activation of muscle derives in part from the end-to-end interactions of tropomyosin, with additional cooperativity mostly below Ca~50~ driven by XB-RU effects and nonlinear crossbridge binding properties. Building on the steady-state models, we developed different models for dynamic muscle function in mouse, rat and human, to investigate the ability of our model of cooperative activation in reproducing physiological activation and relaxation of muscle. We have parameterized the model to these three different species using only changes in the affinity of TnI for TnC (*K* ~DI~) and crossbridge cycling rates (*k* ~M-~). Our results in [Table 4](#pcbi.1004376.t004){ref-type="table"} and [Fig 7](#pcbi.1004376.g007){ref-type="fig"} show it is possible to reproduce the different kinetics and accommodate large changes in Ca^2+^ transients without changing the properties of TnC which are generally thought to be highly conserved. Our results are consistent with a highly conserved underlying mechanism for cooperative activation, while the relative order of crossbridge cycling rates is highly species dependent as expected from the differences in heart rates, with mouse faster than rat, and rat faster than human.

The model also naturally reproduces the Ca^2+^-dependence of the rate of force redevelopment *k* ~tr~, and is able to explain the strong Ca^2+^ dependence of force redevelopment as a result of transient blocking of RU's ([Fig 8](#pcbi.1004376.g008){ref-type="fig"}). Specifically, detachment of crossbridges causes a steady-state activation that is lower (similar to the blebbistatin experiments in the steady-state models), and results in blocking of RU's due to a decrease in XB-RU effects. At high and low Ca^2+^, the filament state is mostly blocked or unblocked, and is not much affected by the number of crossbridges. However, near Ca~50~, this state is particularly sensitive to crossbridges and XB-RU effects. As crossbridges are detached by the fast length change, some RU's move to a blocked position, and require additional time to become unblocked again through XB-RU cooperative interactions. This causes a higher *k* ~tr~ compared to that seen at high Ca^2+^ concentrations. The choice of constant or tropomyosin-deformation dependent crossbridge unbinding rates significantly affects our results for *k* ~tr~. Firstly, the Ca^2+^-dependence is stronger for the case with constant unbinding rate, and better represents experimental observations of 4--8 fold changes, which suggests that the choice of a constant crossbridge unbinding rate may be the more physiological one. However, reproducing physiological relaxation rates requires particularly high unbinding rates (≈ 2/ms in mouse). A potential explanation for this is that velocity dependent effects are important even in isometric tension relaxation due to internal sarcomere shortening \[[@pcbi.1004376.ref076]\]. Secondly, the behaviour at low Ca^2+^ is different, with the variable unbinding rate showing increasing *k* ~tr~. This can be explained by considering the high energy barrier for a single crossbridge binding to a fully blocked thin filament. Experimental results rarely show this region, most likely due to the experimental noise dominating the near-zero force. Regardless, our results and interpretation suggest a Hill fit as applied previously (e.g. \[[@pcbi.1004376.ref055]\]) may not be a suitable choice for these results, as the Ca^2+^-*k* ~tr~ curves show a clear theoretical minimum.

The current model is designed to investigate the protein-protein interactions responsible for cooperativity. As a computational model can not capture every detail, we have made a range of assumptions to make the model computationally tractable, such as the regular spacing of fixed myosin crossbridges, homogeneity of tropomyosin properties, the choice of two global TnC buffers, and several techniques for reducing the number of states. For a general model of contraction, the most important limitations in the current formulation of the model is the absence of velocity- and length-dependent effects. The current model uses a very simplified sarcomere geometry with regularly spaced crossbridges which are all capable of binding to actin, and does not include the effects of filament overlap which prevents crossbridges from binding to specific regions on the thin filament, even at resting length. Extending the model to represent length-dependence of tension would require an accurate model of filament overlap and other details of sarcomere geometry. This in turn would require changes to the current Monte Carlo sampling strategies and representative state approach, as currently equivalent states would have different crossbridge binding properties depending on the location of unblocked units on the thin filament. However, with respect to the effects of length-dependent activation, the explicit representation of troponin I in our model makes it especially suitable for testing current hypotheses relating to length-dependent modulation of troponin I \[[@pcbi.1004376.ref077]--[@pcbi.1004376.ref079]\]. The addition of velocity dependence would allow us to investigate the interaction between velocity-dependent unbinding and XB-RU effects in ventricular relaxation. Although spatially detailed velocity dependence would likely be too computationally demanding for a practical model, the current dynamic model framework is suitable for extensions using averaged strain-distortion approaches \[[@pcbi.1004376.ref080], [@pcbi.1004376.ref081]\].

To limit model complexity, the current model uses simple two-state crossbridge kinetics in which the transitions which affect tropomyosin deformation are considered to be most relevant. Another possible extension of the model is to incorporate more detailed crossbridge dynamics, using a 3-state (weak, strong, unbound) \[[@pcbi.1004376.ref035]\] or more complex model which incorporates ATP/ADP/P~i~ kinetics \[[@pcbi.1004376.ref082]\]. Such extensions are not expected to affect the steady-state results, as our current model can be interpreted as a steady-state approximation of a more complex modelling framework. These extensions could significantly affect dynamic model results and would also allow for more detailed velocity-dependent binding rates between different crossbridge states. Extensions with ATP kinetics would also allow for more quantitative comparison with experimental data showing ATP-dependence of force, as shown in the rigor crossbridge tests in [Fig 6B](#pcbi.1004376.g006){ref-type="fig"}. However, spatially detailed effects would have to be simplified for an extended dynamic model to keep such an extended model tractable, using techniques such as those applied for blocked and unblocked TnC states (TnC~U~, TnC~B~) in the current model.

Finally, the model currently consists of 750 ODEs, which is a relatively high number compared to many phenomenological contraction models currently available. Models which explicitly represent interactions between RU's range in complexity from a few ODEs to several thousands. For example, the work by Razumova et al. uses the assumptions of a uniform spatial distribution of RU states, which allows them to reproduce cooperative effects with a lower number of ODEs \[[@pcbi.1004376.ref083]\]. However, cooperative activation of RU's breaks this assumption as unblocked units tend to cluster together. Furthermore, states with identical number of unblocked RU's can have very different crossbridge binding properties (see [Fig 2](#pcbi.1004376.g002){ref-type="fig"}). More recent work from this group no longer assumes a uniform spatial distribution, but requires several thousand ODEs to represent 9 RU's \[[@pcbi.1004376.ref038]\], as does the model by Dobrunz et al. which takes a similar approach \[[@pcbi.1004376.ref084]\]. Although a high number of ODEs has been used before in whole organ models \[[@pcbi.1004376.ref039]\], and is thus not an immediate computational problem, it limits their portability and ease of use. Thus, model reduction strategies will form an important part of future work. A potential advantage of our modelling approach is its basis in an underlying biophysical model of the tropomyosin chain. This approach allows for longer-range interactions than these previous ODE models including only nearest-neighbour interaction, and results in fewer free parameters compared to approaches where RU-RU, XB-RU and XB-XB effects are defined by independent parameters. A potential disadvantage of our model is the requirement for a pre-calculation using Monte-Carlo sampling, which is computationally costly. To increase the models usability, an implementation including the results of the Monte Carlo sampling procedures are made available online at [cemrg.co.uk](http://cemrg.co.uk). Future work in both model reduction and deformation-dependence will further increase the usability of this biophysically detailed model in a whole-organ setting, and improve the predictive power of multi-scale cardiac models in biological and clinical applications.

Supporting Information {#sec016}
======================

###### Supporting information with model development details.

The supporting information describes: (1) Details for the model for tropomyosin as a flexible chain and the finite element solution of these equations. (2) Details for the derivation of [Eq 5](#pcbi.1004376.e012){ref-type="disp-formula"} for the tropomyosin state energy. (3) The Monte Carlo sampling algorithm used to approximate the result of [Eq 8](#pcbi.1004376.e015){ref-type="disp-formula"}. (4) Details of force calculation and Hill curve fitting. (5) A description of the experimental data used for rat and human calcium transients.

(PDF)

###### 

Click here for additional data file.

###### Effect of spatial arrangement of the thin filament.

This figure demonstrates the importance of the spatial arrangement of the thin filament in the continuous flexible chain model. The left column represents the normal 38.5 nm spacing on which all our models are based. The middle column uses half spacing (19.25 nm), effectively halving the width of an RU, but shows more RU's unblocked to result in an identical length of the tropomyosin chain being unblocked. The right column also uses half spacing and shows an identical number of 5 RU's unblocked as the left column. All differences in free energy Δ*E* are relative to the fully blocked state shown in the top plot in each column. There is a minor difference in the energy freed in unblocking between the left and middle column (A to D vs B to E) due to effect of the spacing of the blocked RU's flanking the unblocked region. However, this effect is significantly smaller than the effect of differences in RU size when an identical number of smaller size RU's are unblocked (e.g. compare D and F). Crossbridge binding is even more similar in the left and middle column (D to G vs E to H) with both changes the free energy by around 1.6*k* ~*B*~ *T*, while in the right column (F to I) the difference is 5*k* ~*B*~ *T*. Thus, the spatial arrangement of the thin filament in the flexible chain model is more influential than the particular number of RU's, because flexible chain models of tropomyosin are not limited to nearest-neighbour interactions. Note that the closer spacing does not represent any particular physiological condition, as troponin spacing is highly conserved across species.

(TIF)

###### 

Click here for additional data file.

###### Description of representative classes and states.

The left plot shows the variation in free energy within each representative class, which is at most 0.11%. The right plot shows the number of states within each class, varying from 1 to 1.2 million. There is a clear tendency for classes with lower energy to have fewer states in each class, as these correspond to longer connected unblocked regions, leading to fewer possible locations within 26 RU. However, due to the symmetries in certain states, classes with a low number of states appear along the entire energy range. Notably, other than the 'fully blocked' and 'fully unblocked' states, which each have their own class, so do the two additional stats `UUUUUUUUBUUUUUUUUBUUUUUUUU`, and `UUBUUBUUBUUBUUBUUBUUBUUBUU`, due to this kind of symmetry.

(TIF)

###### 

Click here for additional data file.

[^1]: The authors have declared that no competing interests exist.

[^2]: Conceived and designed the experiments: SL SAN. Performed the experiments: SL. Analyzed the data: SL. Wrote the paper: SL SAN.
