Active behavior of skeletal muscles is associated with time scales of about 30 ms  [1,2]. At shorter times ($\sim$1 ms) muscles exhibit a nontrivial passive response: if a tetanized muscle is suddenly extended, it comes loose, and if it is shortened, it tightens up with apparently no involvement of Adenosine Triphosphate (ATP)  [3–7]. This unusual mechanical behavior, associated with the unfolding of the attached myosin cross-linkers (cross-bridges), qualifies muscles as metamaterials  [8]. As we argue below, an important factor in this behavior is the dominance of parallel connections with multiple shared links entailing cooperative effects; see Fig. 1. Similar mean-field coupling can be found in many hierarchical biological systems  [9–11]; in particular, it plays a crucial role in cell adhesion, where individual binding elements interact through a common elastic background  [12].

Interaction-induced synchronization during muscle contractions reveals itself through macroscopic fluctuations and spatial inhomogeneities  [13–17]. In ratchet-based and chemomechanical models such collective behavior is usually attributed to breaking of the detailed balance  [18–21], with long range interactions entering the problem implicitly as a force dependence of the chemical rates  [22–25]. However, the cooperative behavior of myosin cross-bridges can be detected during short time force recovery  [26–29], and therefore the origin of synchronization should be within reach of models disregarding disequilibrium and activity. In this Letter, we show that already equilibrium response of “muscle material” is associated with highly synchronized behavior at the microscale which explains its unusual passive response.

In particular, we show that an order-disorder phase transition is displayed by the celebrated Huxley-Simmons (HS) model  [3] if, instead of physiological isometric loading conditions (length clamp) also known as a hard device, one considers isotonic (load clamp) loading conditions or a soft device. While a considerable difference in behavior of a muscle loaded in these two different ways can be deduced from experimental data  [30–32], the origin of the disparity has been so far unexplained. We argue that behind it is a nonequivalence of equilibrium ensembles ubiquitous in systems with long range interactions  [33].

Most remarkably, we find that a careful parameter fitting places the actual skeletal muscle almost exactly into a ferromagnetic Curie point. This agrees with the observation  [3] that the effective stiffness of skeletal muscles associated with fast force recovery is close to zero in the state of isometric contractions and strongly suggests that muscles are finely tuned to perform near marginal stability. Other experimentally observed manifestations of criticality include kinetic slowing down and large scale macroscopic fluctuations near the stall force conditions  [26–29].

We demonstrate the robustness of our predictions by comparing the HS model, which we interpret as a hard spin description, with a regularized (RHS) model where filament elasticity is taken into consideration and conventional spins are replaced by elastic snap springs.

The HS model.— We consider a prototypical model of a half-sarcomere with $N$ attached cross-bridges arranged in parallel  [3]; at time scales of fast force recovery $N$ can be considered constant  [5,31]. Each cross-bridge is modeled as a series connection of a bistable spin unit and a linear (shear) spring; see Fig. 2. We use dimensionless variables with the power-stroke size $a$ as a unit of displacement and $\kappa a^{2}$ as the unit of energy, where $\kappa$ is the stiffness of the series spring. Then the spin variable takes values $x_{i}=0$ (pre-power-stroke state) and $x_{i}=-1$ (post-power-stroke state), and the total energy per particle in the hard device is $v(\boldsymbol{x},z) = (1/N) \sum [(1+x_{i})v_0 + (1/2)\left( z-x_{i} \right)^{2}]$, where $v_0$ is the energetic bias of the pre-power-stroke state. In this formulation, the HS model describes the simplest paramagnetic spin system.

At finite temperature $\theta$, the equilibrium behavior of this system is characterized by the free energy per particle $\hat{f}(z,\beta)=-[1/ (N \beta)]\ln \int \exp\left[-\beta Nv(\boldsymbol{x},z)\right]d\boldsymbol{x}$ where $\beta = \kappa a^{2}/(k_{b}\theta)$. At fixed $p=-(1/N) \sum x_{i}$, representing the fraction of cross-bridges in the post-power-stroke state, the macroscopic ( $N\to\infty$ ) free energy takes the form $$f = p\biggl[\frac{1}{2}( z+1)^{2}\biggr] + (1-p)\biggl[ \frac{1}{2}z^{2} + v_{0}\biggr] + \frac{1}{\beta}S(p)$$ where $S(p)= p\log(p) + (1-p)\log(1-p)$. The function $f(p)$ is always convex with a minimum at $\hat{p}(z,\beta) = 1/2 -(1/2)\tanh[(\beta/2)(z-v_{0}+1/2)]$. The equilibrium tension per cross-bridge is then $t = (z+\hat{p} )$, which is exactly the formula found by Huxley and Simmons. The equilibrium free energy is $$\hat{f} = - \frac{1}{\beta}\ln\!\biggl\{\exp\biggl[-\frac{\beta}{2} (z+1)^2\biggr]\,+\,\exp\biggl[-\beta\biggl(\frac{z ^{2}}{2} +v_0\biggr)\biggr]\biggr\}$$ and the susceptibility (stiffness) is $t'(z)=\hat{f}''(z)$. The function $\hat{f}(z)$ is convex for $\beta<4$ and is nonconvex for $\beta>4$ exhibiting a range with negative equilibrium stiffness (metamaterial behavior). If we take the values from  [3], $a=8$ nm and $\kappa a^{2}/2 = 2 k_{b}\theta$, we obtain $\beta = 4$, which corresponds to zero stiffness at the state of isometric contractions $z_0=v_{0}-1/2$.

In the soft device setting, not studied by Huxley and Simmons, the energy becomes $w(\boldsymbol{x},z,t) = v(\boldsymbol{x},z) - tz$, where $t$ is the applied force per particle. Now the variable $z$ plays the role of an internal parameter whose adiabatic elimination produces a Curie-Weiss mean-field potential depending on $\left(\sum x_{i}\right)^2$. The equilibrium Gibbs free energy is now $\hat{g}(t,\beta)=-[1/ (N \beta)]\ln \int \exp\left[-\beta Nw(\boldsymbol{x},z,t)\right]d\boldsymbol{x}dz.$ At fixed $p$ we obtain in the thermodynamic limit $$g = -\frac{1}{2} t^{2} + pt +(1-p) v_{0}+ \frac{1}{2}p(1-p) + \frac{1}{\beta}S(p),$$ where the “regular solution” term is responsible for cooperative (ferromagnetic) behavior.

In Fig. 3, we show the position of the minima of $g(p)$ when $t$ is chosen to ensure that in the paramagnetic phase $\hat{p}(t,\beta)=1/2$. In the disordered (high temperature) state all cross-bridges are in random conformations, while in the ordered (low temperature) state the system exhibits coherent fluctuations between post-power-stroke and pre-power-stroke configurations. These fluctuations describe temporal microstructures responsible for the plateau in the force-elongation relation $\hat{z} = t - \hat{p}$, where $\hat{p}$ is a solution of $t-\hat{p}+1/2-v_{0}+ (1/\beta)\ln[\hat{p}/(1-\hat{p})]=0$.

The equilibrium Gibbs energy is concave because $\hat{g}''(t) = -1-\beta N \langle\left(p-\hat{p}\right)^{2}\rangle\leq 0$, so in the soft device the stiffness is always positive. Since in the hard device the stiffness is sign indefinite, the two ensembles are not equivalent. This is expected for systems with strong long range interactions that are inherently nonadditive  [33–35]. Negative stiffness in the hard device HS model has been known for a long time  [3,23,36–39]; however, it was not previously associated with the particular internal architecture of muscle material.

As we have already mentioned, the original HS fit of experimental data  [3] places the system exactly into the critical state (Curie point). In this state the correlation length diverges and fluctuations become macroscopic, which is consistent with observations at stall force conditions  [14,40–42]. This suggests that skeletal muscles, as many other biological systems, may be tuned to criticality. The proximity to the critical point would then be the result of either evolutionary or functional self-organization. The marginal stability of the critical state allows the system to amplify interactions, ensure strong feedback, and achieve considerable robustness in front of random perturbations. In particular, it is a way to quickly switch back and forth between highly efficient synchronized stroke and stiff behavior in the desynchronized state.

The ensembles nonequivalence in the HS model has also a kinetic signature. Experiments on quick recovery reveal that muscle fibers react to load steps much slower than to length steps  [7,31,32]. This agrees with our model, where coherent response (in isotonic conditions) is always slower than disordered response (in isometric conditions).

Indeed, by using Kramers approximation Huxley and Simmons obtained in a hard device the kinetic equation $\dot{\langle p \rangle} = -k_{-}\langle p \rangle + k_{+}(1-\langle p \rangle)$, where $\langle p \rangle$ is the average over ensemble. The constants $k_{+}, k_{-}$ satisfy the detailed balance $k_{+}/k_{-}=\exp\left[-\beta(z-v_{0}+1/2)\right]$, and the recovery rate is $1/\tau= k_{-}\{1+\exp[-\beta(z-v_{0}+1/2)]\}$. In a soft device we may use the same model with $z = t-\langle p \rangle$, which accounts for force-dependent chemistry and introduces nonlinear feedback. The characteristic rate around a given state $\langle p \rangle$ is then $1/\tau = k_{-}\{1+[1-\beta(1-\langle p \rangle)]\exp[-\beta(t-\langle p \rangle-v_{0}+1/2)]\}$. When $\langle p \rangle$ is small, $t-\langle p \rangle> z$, and the relaxation in a soft device is slower than in a hard device. In Fig. 4, we show the rates obtained from the HS model; in the case of a soft device, the nonlinear kinetic equation was solved numerically for the duration 10 ms. We see that in a soft device the rates are indeed slower than in a hard device, however, the experimental measurements are not matched quantitatively.

The RHS model.— To test the robustness of the HS mechanism of synchronization and to achieve quantitative agreement with kinetic data, we now consider a natural regularization of the HS model. First, following  [38] we replace hard spins by soft spins described by a piecewise quadratic double well potential – see Fig. 2(a)–  $$u_{\text{\tiny RHS}}(x) = \begin{cases} \frac{1}{2}\lambda_0(x)^{2} + v_{0} & \text{if $ x> l $,}\\ \frac{1}{2}\lambda_1(x + 1)^{2} & \text{if $ x\leq l $,} \end{cases}$$ where $\lambda_{1} = \kappa_{1}/\kappa$, $\lambda_{0} = \kappa_{0}/\kappa$. Second, we introduce a mixed device (mimicking myofilament elasticity  [18,43,44]) by adding to our parallel bundle a series spring. The resulting energy per cross-bridge in a hard device is $$v(\boldsymbol{x},\!y;\!z)=\frac{1}{N} \sum_{i=1}^{N} \big[u_{\text{\tiny RHS}}(x_{i}) + \frac{1}{2}( y-x_{i})^{2}\big]+\frac{\lambda_{f}}{2}(z-y)^{2},$$ where $y$ is a new internal variable and $\lambda_{f} = \kappa_{f}/(N\kappa)$; see Fig. 2(b). It is clear that our lump description of filament elasticity misrepresents short range interactions  [45–47]; however, this should not affect our main results  [48,49].

To study the soft device case we must again consider the energy $w(\boldsymbol{x},y,z,t)=v(\boldsymbol{x},y,z)-tz$, where $t$ is the applied force per cross-bridge. A transition from hard to soft ensemble is made by taking the limit $\lambda_{f}\rightarrow 0$, $z \rightarrow \infty$ with $\lambda_{f}z \rightarrow t$. At finite $\lambda_{f}$ the RHS model can be viewed as a version of the mean-field $\varphi^4$ model studied in  [10,11,33].

The HS model is a limiting case of the RHS model with $\lambda_{1,0}\rightarrow\infty$ and $\lambda_{f}\rightarrow\infty$. The first of these limits allows one to replace continuous dynamics by jumps and use the language of chemical kinetics; however, it also erases information about the barriers; see  [38]. The second limit eliminates the Curie-Weiss (mean-field) interaction among individual cross-bridges at fixed $z$, and that is why the synchronized behavior was overlooked in  [3].

Equilibrium behavior in the RHS model can be again described analytically, because it is just a redressed HS model. In the limit $\lambda_{f}\to\infty$ the function $f (p)$ is convex as in HS model, while at finite values of $\lambda_{f}$ it is now nonconvex. This shows that in the RHS model the account of filament elasticity brings about phase transition (and bistability) also in the hard device.

The bistable nature of the macroscopic free energy in both soft and hard devices implies that the system can be in two coherent states, and therefore within a large set of half-sarcomeres one should expect observable spatial inhomogeneities. This prediction is in agreement with ubiquitous “off-center” displacements of M lines recorded during isometric contractions  [14].

The phase diagram showing the role of filament elasticity in hard device is shown in Fig. 5. The dependence of the critical temperature on $\lambda_{f}$ suggests that actomyosin systems can control the degree of cooperativity by tuning the internal stiffness; likewise, variable stiffness of the loading device may be used in experiments to either activate or deactivate the collective behavior. Notice that the realistic choice of parameters again selects a near critical state; the exact criticality is compromised since the symmetry between the pre- and post-power-stroke states is now broken (as $\lambda_{1} \neq \lambda_{0}$) and the phase transition becomes weakly first order; see Fig.6.

A behavior similar to our synchronization has been also observed in the models of passive adhesive clusters, where the elastic feedback appears as strain- or force-dependent chemistry  [12]. Given that the two systems exhibit almost identical cooperative behavior, we expect criticality to be also a factor in the operation of focal adhesions.

The two ensembles, soft and hard, remain inequivalent in the RHS model. Thus, in the soft device the equilibrium Gibbs free energy $\hat{g}$ is concave since $\hat{g}''(t) = -1/\lambda_{f} - \beta N \langle(y-\langle y\rangle)^{2}\rangle\leq 0$ which means that the stiffness is always positive. Instead in the hard device $\hat{f}''(z) = \lambda_{f}[1-\beta N \langle(y- \langle y\rangle)^{2}\rangle]$, and the stiffness can be both positive and negative. While negative stiffness should be a characteristic feature of realistic half-sarcomeres (see Fig. 7), it has not been observed in experiments on whole myofibrils. The reason may be that in myofibrils a single half-sarcomere is never loaded in a hard device. The effective dimensionless temperature may also be higher because of the quenched disorder, and the stiffness may be smaller due to nonlinear elasticity. One can also expect the unstable half-sarcomeres to be stabilized actively through processes involving ATP hydrolysis  [50].

To study kinetics in the RHS model we perform direct numerical simulations by using a Langevin thermostat. We assume that the macroscopic variables $y$ and $z$ are fast and are always mechanically equilibrated which is not an essential assumption. The response of the remaining variables $x_i$ is governed by the system $dx_{i} = b(x_{i})dt + \sqrt{2\beta^{-1}}dB_{i}$, where the drift is $b (\boldsymbol{x},z) = -u'_{\text{\tiny RHS}} (x_{i}) + (1+\lambda_{f})^{-1}( \lambda_{f}z + N^{-1}\sum x_{i}) - x_{i}$ in a hard device and $b (\boldsymbol{x},t) = -u'_{\text{\tiny RHS}} (x_{i}) + t + N^{-1}\sum x_{i} - x_i$ in a soft device. In both cases the diffusion term $dB_i$ represents a standard Wiener process.

In Fig. 7, we show the results of stochastic simulations imitating quick recovery experiments  [3]. The system, initially in thermal equilibrium at fixed $z_{0}$ (or $t_0$), was perturbed by applying fast ($\sim100\ \mu$s) length (or load) steps with various amplitudes. In a soft device the system was not able to reach equilibrium within the experimental time scale. Instead, it remained trapped in a quasistationary (glassy) state because of the high energy barrier associated with collective power stroke. Such kinetic trapping which fits the pattern of two-stage dynamics exhibited by systems with strong long range interactions  [33,51,52] may explain the failure to reach equilibrium in experiments reported in  [27–29]. In the hard device case, the cooperation among the cross-bridges is weaker and kinetics is much faster, allowing the system to reach equilibrium at the experimental time scale. A quantitative comparison of the rates obtained in our simulations with experimental values (see Fig. 4) shows that the RHS model reproduces the kinetic data in both hard and soft ensembles rather well.

In conclusion, we mention that the prototypical nature of our model implies that passive collective behavior should be a property common to general cross-linked actomyosin networks. We have shown that the degree of cooperativity in such networks can be strongly affected by elastic stiffness of the filaments. This suggests that a generic system of this type can be tuned to criticality by an actively generated prestress  [53].

References