Muscle activation described with a differential equation model for large ensembles of locally coupled molecular motors

Molecular motors, by turning chemical energy into mechanical work, are responsible for active cellular processes. Often groups of these motors work together to perform their biological role. Motors in an ensemble are coupled and exhibit complex emergent behavior. Although large motor ensembles can be modeled with partial differential equations (PDEs) by assuming that molecules function independently of their neighbors, this assumption is violated when motors are coupled locally. It is therefore unclear how to describe the ensemble behavior of the locally coupled motors responsible for biological processes such as calcium-dependent skeletal muscle activation. Here we develop a theory to describe locally coupled motor ensembles and apply the theory to skeletal muscle activation. The central idea is that a muscle filament can be divided into two phases: an active and an inactive phase. Dynamic changes in the relative size of these phases are described by a set of linear ordinary differential equations (ODEs). As the dynamics of the active phase are described by PDEs, muscle activation is governed by a set of coupled ODEs and PDEs, building on previous PDE models. With comparison to Monte Carlo simulations, we demonstrate that the theory captures the behavior of locally coupled ensembles. The theory also plausibly describes and predicts muscle experiments from molecular to whole muscle scales, suggesting that a micro- to macroscale muscle model is within reach.

Figure 1

Cartoons of various aspects of activation. (a) Molecular events of vertebrate skeletal muscle activation. Tropomyosin can be found in three positions on actin: the blocked, closed, and open states. Calcium binding to troponin favors a shift from the blocked to the closed states. Assuming these states are in equilibrium, we model both with an equivalent blocked or closed state (dark). Myosin binding causes a further shift in tropomyosin position into the open state (light), a transition that occurs more easily at high calcium when the equilibrium position of tropomyosin favors the closed state. (b) Schematic of the continuous flexible chain model. Tropomyosin is treated as a slender infinite linearly elastic beam in a plane. The interaction of tropomyosin with actin and troponin is treated as a potential energy density $W\left(z\right)$. Myosin binding induces a displacement into the open state and the shape of tropomyosin minimizes the overall energy (which includes the energy of elastic deformation and the change in the actin-troponin-tropomyosin interaction energy). (c) Simplification of the continuous flexible chain model. When nearby myosins bind to actin (myosins 2 and 4), they induce a uniform displacement of tropomyosin into the open state. When distant myosins bind to actin (myosins 4 and 7), they induce independent deformations of tropomyosin. (d) Simplified continuous flexible chain model predicts binding energy $(\Delta V$ measured in units of $k_{B}T)$ as a function of the distance to a bound myosin $L$, which can be defined by two parameters $\mathcal{C}$ and $\varepsilon$ that then define local coupling. (e) At low levels of activation, myosin ensembles with strong local coupling show clusters of locally active molecules. The middle shows 1500 molecules that are either unbound (indicated by $u)$ or bound (indicated by $b)$. The bottom shows the attachment rate constant for every molecule. There are short regions where $k_a=40{\mathrm{s}}^{-1}$. Here tropomyosin is in the open state. The top shows a more detailed picture of a cluster of activated molecules that spans 24 myosins.

Figure 2

The theory of locally coupled ensembles is consistent with Monte Carlo simulations of a two-state model in steady state. The model is defined by two parameters that determine cooperativity $\mathcal{C}$ and $\varepsilon$ and three parameters that determine myosin's binding to actin $D,K$, and $E$. In all plots, $\mathcal{C}=11$ and $D=100$. (a) Steady-state binding probability $P_b$ as a function of the coupling parameter $\varepsilon$ in the absence of an external force. With $E=0.5$, three different values of $K$ are shown. Predictions of the theory in the sparse binding limit (solid line) and dense binding limit (dashed line) fit the simulations. (b) Steady-state binding probability $P_b$ as a function of $\varepsilon$ under isometric conditions. With $K=0.25$, three different values of $E$ are shown. Predictions of the theory in the sparse binding limit (solid line) and dense binding limit (dashed line) fit the simulations. (c) The isometric steady-state cluster distribution predicted by the theory (solid line) agrees with Monte Carlo simulations (dots are the mean of 20 simulations and error bars indicate standard error). The parameters are $ln\left(\varepsilon \right)=-6,E=1$, and $K=0.25$. (d) The isometric steady-state void distribution predicted by the theory (solid line) agrees with Monte Carlo simulations (dots are the mean of 20 simulations and error bars indicate standard error). The parameters are $ln\left(\varepsilon \right)=-5,E=1$, and $K=0.25$.

Figure 3

The theory of locally coupled ensembles is consistent with the time dependence of Monte Carlo simulations of a two-state model. In all plots, $\mathcal{C}=11,D=100,K=0.25$, and $E=1$ and time is measured in units of $1/k_a$. In (a) and (b), long simulations were run with ensembles of $M=1000$ molecules; short simulations (shown in insets) had $M=5000$. We did 20 repetitions of these simulations. (a) Under isometric conditions, the dense binding limit is consistent with simulations. Simulations started from steady state with tropomyosin in the open state and local coupling $\left[ln\left(\varepsilon \right)=-5\right]$ was initiated at $t=0$. The left plot shows the average binding probability for the theory (solid line) and simulations (the dark gray line is the mean and the light gray region the standard deviation). The inset shows a larger ensemble size. The right plot shows the void distribution at two different times (symbols are the mean and error bars the standard error; lines are the theory). The inset shows the time dependence of voids of size 1 (the dark gray line is the mean, the light gray region is the standard error, and the black line is the theory). (b) The sparse binding limit is consistent with isometric simulations starting with no bound molecules with $ln\left(\varepsilon \right)=-6$. The left plot shows the average binding probability for the theory (solid line) and simulations (the dark gray line is the mean and the light gray region the standard deviation). The inset shows a shorter set of simulations with a larger ensemble size. The right plot shows the cluster distribution at two times (symbols are the mean and error bars the standard error; lines are the theory). The inset shows the time dependence of clusters of size 1 as a function of time (the dark gray line is the mean and the light gray region the standard error; the black line is the theory). (c) The theory is consistent with simulations where force is varied dynamically. In the top two plots, the mean of ten simulations is the dark gray line and the standard deviation is the light gray region. The sparse binding limit is the black line. The bottom plot shows force (in units of $kd)$ as a function of time. The theory successfully predicts actin position (in units of $d$, middle plot) and binding probability (top plot) as a function of time. In these simulations, $ln\left(\varepsilon \right)=-6$.

Figure 4

Splicing together the sparse and dense binding limits allows the model to fit simulations that transition from one limit to the other. The time evolutions of (a) an ensemble with no initially bound molecules and (b) an ensemble starting from steady state with tropomyosin in the open state (bottom) are shown. In all plots, $K=0.25,E=1,D=100$, and $\mathcal{C}=11$. With the exception of the top curve in (b) (which has 20 simulations), simulations contain 40 simulations of ensemble size $M=1000$. The mean is the dark gray line and the standard deviation is the light gray region. The spliced theory is shown as a black line with splices (when present) indicated by an arrow; for reference, the sparse and dense binding limits are shown as dashed lines in (a) and (b), respectively.

Figure 5

The theory describes realistic model simulations and published experiments at various size scales. Here $p\mathrm{Ca}$ is the negative logarithm of the calcium concentration. (a) Simulations of the model with $\mathcal{C}=13$ and $\varepsilon =0.0012$ are consistent with rigor activation measurements, where actin speed as a function of ATP is measured at high and low calcium [40] (dots show the mean and error bars the standard deviation for ten simulations). (b) One-parameter fit (solid line, a gap separates the sparse and dense limits) to isometric force $p\mathrm{Ca}$ data from skinned muscle fibers [81]. To demonstrate that the theory reasonably reflects the model's behavior, the results of Monte Carlo simulations are also shown (closed circles and inset). (c) The predicted motility speed $p\mathrm{Ca}$ relationship, using parameters from the fits to isometric force $p\mathrm{Ca}$ and rigor activation data, is consistent with recent measurements [40, 82]. (d) Theory and simulation agree when $\varepsilon$ is varied dynamically. To show that the theory can describe the activation of intact fibers or whole muscle, where a train of electrical impulses of variable frequency activates a muscle, we simulated square waves in $\varepsilon$ at slow (S, bottom), medium (M, middle), and fast (F, top) frequency. In all cases, the theory (black line) reasonably reproduces the average of three simulations (gray line). A splice occurs in the theoretical curve for the fast frequency at $F_{\mathrm{iso}}/F_{\max }=0.5$. (e) The theory reproduces activation of intact fibers and whole muscle. The top left plot shows the theoretical response from a single electrical impulse at $t=0.1\mathrm{s}$, showing calcium concentration (bottom) and force (top) as a function of time. Shown on the right is the theoretical force response to a short series of electrical impulses (doublets and quintets) at variable frequency (cf. Fig. 3 in Ref. [83]). The bottom left plot shows that the predicted force response to a long series of electrical impulses of variable frequency is consistent with experiments [84, 85, 86] $(f_{50}$ is the frequency of half-maximal force).