Passive viscoelastic response of striated muscles

TL;DR

This paper develops a novel second-order rheological model for the passive viscoelastic response of striated muscles, incorporating long-range interactions and validated with frog muscle data.

Contribution

It introduces a new second-order continuum model for muscle mechanics derived from a stochastic micromodel, improving upon previous first-order approaches.

Findings

Model shows excellent agreement with experimental data

Second-order equation captures complex muscle behavior

Microscopic parameters accurately predict macroscopic response

Abstract

Muscle cells with sarcomeric structure exhibit highly nontrivial passive mechanical response. The difficulty of its continuum modeling is due to the presence of long-range interactions transmitted by extended protein skeleton. To build a rheological model for muscle 'material' we use a stochastic micromodel and derive a linear response theory for a half-sarcomere. Instead of the first order rheological equation, anticipated by A.V. Hill on the phenomenological grounds, we obtain a novel second order equation. We use the values of the microscopic parameters for frog muscles to show that the proposed rheological model is in excellent quantitative agreement with physiological experiments.

Figures

Equations

$$H=\sum _{i=1}^N\left[V(x_i)+\frac{(y-x_i{)}^2}{2}+\frac{{\lambda }_f}{2}{\left(z-y\right)}^2\right]-zf,$$

$$\frac{\nu }{N}\dot{z}={\lambda }_f\frac{⟨x⟩-z}{1+{\lambda }_f}+\frac{f}{N}.$$

$$\frac{\partial p}{\partial t}=\frac{\partial }{\partial x}\left[\left(\frac{\partial V}{\partial x}+x-\frac{⟨x⟩+{\lambda }_{f}z}{1+{\lambda }_f}+\frac{1}{\beta }\frac{\partial }{\partial x}\right)p\right]\text{mbox},$$

$$U(x,z)=V(x)+\frac{1}{2}{\left(x-\frac{⟨x{⟩}_s(z)+{\lambda }_{f}z}{1+{\lambda }_f}\right)}^2$$

$$Lp=\frac{\partial }{\partial x}\left[\left(\frac{\partial V}{\partial x}+x-\frac{⟨x{⟩}_s+z{\lambda }_f}{1+{\lambda }_f}+\frac{1}{\beta }\frac{\partial }{\partial x}\right)p\right]\text{mbox}.$$

$$\frac{\partial \delta p}{\partial t}=L\delta p-\frac{1}{1+{\lambda }_f}\frac{\partial p_s}{\partial x}\left(\int \text{mbox}dxx\delta p+{\lambda }_f\delta z\right).$$

$$⟨\delta x⟩(t)={\int }_{-\infty }^{\infty }\frac{⟨\delta x(t^{\prime })⟩+{\lambda }_f\delta z(t^{\prime })}{1+{\lambda }_f}{\chi }_{xx}(t-t^{\prime })\text{mbox}d{t}^{\prime },$$

$${\chi }_{xx}(t)=-\Theta (t)\int \text{mbox}dxx{e}^{Lt}\frac{\partial }{\partial x}{p}_s(x)\text{mbox},$$

$$S_{xx}(t)\simeq c{e}^{-\frac{t}{{\tau }_x}},$$

$$\begin{gathered} \left[-{\omega }^2{\tau }_x{\tau }_z+i\omega \left({\tau }_z+{\tau }_x-\frac{\beta c{\tau }_z}{1+{\lambda }_f}\right)+1-\beta c\right]\delta z(\omega )= \\ \left[i\omega {\tau }_x\frac{1+{\lambda }_f}{{\lambda }_{f}N}+\frac{1+{\lambda }_f-\beta c}{{\lambda }_{f}N}\right]\delta f(\omega )\text{mbox}. \end{gathered}$$

$$\theta \nu \ddot{\delta z}+(\theta E+\nu )\dot{\delta z}+C\delta z=\delta f+\theta \dot{\delta f}$$

$$\begin{gathered} \delta z(t)=\delta f_0[-\frac{{\tau }_{+}{\tau }_{-}}{\nu ({\tau }_{+}-{\tau }_{-})}\left(1-\frac{{\tau }_{+}}{\theta }\right)e^{-\frac{t}{{\tau }_{+}}} \\ +\frac{{\tau }_{+}{\tau }_{-}}{\nu ({\tau }_{+}-{\tau }_{-})}\left(1-\frac{{\tau }_{-}}{\theta }\right)e^{-\frac{t}{{\tau }_{-}}}+\frac{1}{C}]\text{mbox}, \end{gathered}$$

$${\tau }_{\pm }=\frac{{\eta }_1+\nu }{2E_1}+\frac{\nu }{2E}\pm \sqrt{{\left(\frac{{\eta }_1+\nu }{2E_1}-\frac{\nu }{2E}\right)}^2+\frac{{\nu }^2}{E_{1}E}}\text{mbox}.$$

Full text

One of the simplest biological systems, that still defies the attempts to reproduce it artificially as a macroscopic material, is a striated muscle Kidambi et al. (2018). Its mechanical complexity is due to the presence of a large number of nonlinear, hierarchically organized microscopic sub-systems that are strongly coupled through long-range interactions Caruel and Truskinovsky (2018). This makes the task of reconstructing the macroscopic constitutive relations describing even its passive mechanical response rather challenging Caruel et al. (2019); Regazzoni et al. (2018).

A broadly used phenomenological theory of the passive visco-elastic response of striated muscles, proposed by A.V.Hill Hill (1938, 1949); Glantz (1974); Gerazov and Garner (2016), does not rely on coarse graining techniques DeVault and McLennan (1965); Bavaud (1987) and therefore does not offer a link between macro and micro parameters. Since Hill’s rheological relation involves a single characteristic time scale, it also does not capture the difference in the passive response exhibited by striated muscles abruptly loaded in soft (isotonic) and hard (isometric) loading devices Huxley and Simmons (1971); Reconditi et al. (2004).

A microscopically guided stochastic approach to muscle visco-elasticity was proposed by Huxley and Simmons Huxley and Simmons (1971); Eisenberg et al. (1980) who assumed that the individual force producing units (myosin cross-bridges) are stochastically independent. A mean-field interaction between the cross-bridges was incorporated in a closely related model by Shimizu Shimizu and Yamada (1972); Kometani and Shimizu (1975); Bonilla (1987); Frank (2004). The two approaches have been recently unified Marcucci and Truskinovsky (2010a, b); Caruel et al. (2013); Caruel and Truskinovsky (2018). In the present article we use this framework to rigorously derive from a micro-model a linear rheological response theory for a muscle half-sarcomere. Our analysis builds on the work of Shiino Shiino (1987) who obtained a similar linear response theory for the related model of Desai and Zwanzig Desai and Zwanzig (1978); Bonilla (1987); Frank (2004); other relevant out-of-equilibrium systems were studied in Patelli et al. (2012); Patelli and Ruffo (2014).

Our main result is the linear spring-dashpot scheme which reproduces adequately the mechanical behavior of a muscle fiber subjected to a time dependent perturbation. In contrast to the classical model of Hill Hill (1938, 1949), the proposed rheological equation contains not only the first but also the *second * time derivatives of the macroscopic displacement. We use the values of the microscopic parameters for frog muscles to show that our macroscopic model, which does not rely on any fitting parameters, is in excellent quantitative agreement with physiological experiment.

The model contains two time scales that can be associated with the transient stages of muscle response known as phases 1 and 2, see Fig. 1 and Huxley and Simmons (1971); Piazzesi et al. (2003); Decostre et al. (2005); Pia (2007), with both of them captured adequately by our rheological model. The presence of the two time scales reflects involvement of the two parallel passive processes: the (microscopic) conformational change in $N$ myosin heads and the macroscopic relaxation of the myofilaments in a viscous environment. The obtained linear rheological equation does not capture the collective barrier crossing which can be associated with the slower phase 3 explained in Caruel and Truskinovsky (2016, 2017, 2018); it also does not address the active phase 4, see Fig. 1.

We recall that striated muscle is a hierarchical chemo-mechanical system with the smallest scale represented by force generating half-sarcomeres Alberts et al. (2013); Caruel and Truskinovsky (2018). To reproduce the passive response of this molecular machine we use the simplest microscopic model developed in Marcucci and Truskinovsky (2010a, b); Caruel et al. (2013); Caruel and Truskinovsky (2016, 2017). We assume that behind the muscle power stroke, responsible for fast force recovery, is a double well potential of Landau type, $V(x)$, describing two conformational states of a cross-bridge Marcucci and Truskinovsky (2010a); it is known that myosin heads can be in other configurations Huxley (1974); Decostre et al. (2005); Caremani et al. (2013), so here we simplify the real physical picture Molloy et al. (1995); Reconditi et al. (2004). We also assume that the potential $V(x)$ is asymmetric, which allows the system to generate (stall) force in the physiological regime of isometric contractions. While this asymmetry is maintained actively Sheshka et al. (2016), we can still interpret the short time response in such system as passive. The active response involving detachment of the cross-bridges becomes dominant at time-scales of the order of $40$ ms Howard (2001).

To define the dimensionless units we rescale the length by the size of the maximum working stroke $a$, the energy by $\kappa\,a^{2}$, where $\kappa$ is the cross-bridge stiffness, and the time by the ratio of the cross-bridge drag coefficient $\gamma_{x}$ and $\kappa$. To model a bundle of thick and thin filaments linked by $N$ cross-bridges and loaded with the force $f$ we use the dimensionless energy Caruel and Truskinovsky (2018)

[TABLE]

where the variables $x_{i}$ represent the configuration of individual cross bridges. The latter are elastically coupled through the cross-bridge stiffness to a collective variable $y$ representing an actin filament; the associated quadratic term in (1) represents mean field type interaction. The combined elasticity of actin and myosin filaments is described by the quadratic term coupling the variable $y$ and $z$ with the dimensionless coefficient $\lambda_{f}$ representing the overall filamental stiffness, see Caruel and Truskinovsky (2018) for more details.

We assume that the meso-scopic collective variable $y(t)$ relaxes instantaneously Huxley et al. (1994); Wakabayashi et al. (1994) and we can therefore eliminate such (Weiss-type) variable adiabatically. We obtain $y=(\langle x\rangle+\lambda_{f}\,z)/(1+\lambda_{f}),$ where we introduced notation $\langle x\rangle=N^{-1}\,\sum_{i=1}^{N}x_{i}.$ Substituting the obtained expression for $y$ into (1) we obtain the redressed energy $\tilde{H}(\textbf{x},z)=H(\textbf{x},y(\textbf{x},z),z)$ whose main feature vis-à-vis the model of Huxley and Simmons is the presence of all-to-all interaction between the cross-bridges expressed through the term $\langle x\rangle^{2}$, see also Desai and Zwanzig (1978); Shiino (1987).

Next we assume that the macro-variable $z(t)$ is also deterministic, however now its dynamics is governed by the relaxation equation: $\nu\dot{z}=-\partial\tilde{H}/\partial z,$ where the superimposed dot denotes time derivative and $\nu$ is the dimensionless macroscopic friction coefficient obtained by normalizing the filamental drag coefficient $\gamma_{z}$ by the myosin head drag coefficient $\gamma_{x}$. We can rewrite this equation as

[TABLE]

The time scale $\tau_{z}=\nu(1+\lambda_{f})/(\lambda_{f}N)$ will then characterize the evolution of the macro-variable $z$.

In the evolution of the micro-variables $x_{i}$ we need to account for the thermal noise: $\dot{x}_{i}=-\partial\,\tilde{H}/\partial\,x_{i}+\xi_{i},$ where $\langle\xi_{i}\rangle=0$ and $\langle\xi_{i}(t)\,\xi_{j}(t^{\prime})\rangle=(2/\beta)\,\delta(t-t^{\prime})\,\delta_{ij}$. Here we introduced the inverse dimensionless temperature $\beta$. Under the assumption that $N$ is large, the single particle probability density $p(x,t)$ can be found from the nonlinear Fokker-Planck equation Desai and Zwanzig (1978); Shiino (1987); Frank (2005)

[TABLE]

where we use approximation $\langle x\rangle\approx\int\mbox{d}x\,p\,x$. The stationary solution of (3) is $p_{s}(x)=Z^{-1}\,\exp{\left[-\beta U\right]}$, where

[TABLE]

is an effective double well potential and $Z$ is a normalization constant Caruel and Truskinovsky (2018). In (4) we must use the self-consistence condition $\langle x\rangle_{s}=\int\mbox{d}x\,p_{s}(x)\,x$; such closure can be justified rigorously in the thermodynamic limit Dauxois et al. (2003).

To develop the linear response theory we use as a starting point the approach developed in Shiino (1987). First, we linearize (3) to obtain the propagator

[TABLE]

A perturbation $\delta p(x,t)=p(x,t)-p_{s}(x)$ associated with a small change of the macroscopic strain variable $\delta z(t)$ will then satisfy a linear equation

[TABLE]

In particular, the collective variable $\langle\delta x\rangle(t)=\int\mbox{d}x\,x\,\delta p(x,t)$, will evolve according to

[TABLE]

which in Fourier space reads $\langle\delta x\rangle(\omega)=\lambda_{f}\,\chi_{xx}(\omega)\delta z(\omega)/(1+\lambda_{f}-\chi_{xx}(\omega)).$ The susceptibility in (7) is defined by the relation

[TABLE]

where $\Theta(t)$ is the Heaviside function.

Next we write the conventional fluctuation-dissipation type identity $\chi_{xx}=-\beta\,\Theta(t)\,dS_{xx}/dt$, where $S_{xx}(t)=\langle x(t)\,x(0)\rangle-\langle x\rangle_{s}^{2}$ is the auto-correlation function of a single element $x_{i}$ which is assumed to be evolving in the effective potential (4). Given that in linear approximation each cross-bridge can be viewed as conducting independently a simple Brownian motion in a double well potential we can use the Kramers approximation to obtain Skinner and Wolynes (1978); Pratolongo et al. (1995):

[TABLE]

where $\tau_{x}\simeq\frac{2\pi}{\sqrt{|U^{\prime\prime}(x_{M})|}}\frac{e^{\beta U(x_{M})}}{\sqrt{U^{\prime\prime}(x_{0})}\,e^{\beta U(x_{0})}+\sqrt{U^{\prime\prime}(x_{1})}\,e^{\beta U(x_{1})}}\,\mbox{,}$ with $x_{0,1}$, being the minima of the potential (4) and $x_{M}$, the local maximum between them. The pre-factor can be also computed analytically $c=(x_{0}-x_{1})^{2}\,q/(1+q)^{2}$ with $q=\sqrt{(U^{\prime\prime}(x_{0})/U^{\prime\prime}(x_{1}))}\,\exp[\beta(U(x_{0})-U(x_{1}))]$. In this computation we neglected the effect of the relaxation within a single well because it is much faster in physiological conditions than the barrier crossing.

If we rewrite (9) in the Fourier space $\chi_{xx}(\omega)=c\,\beta/(1+i\,\omega\,\tau_{x}),$ and use (2), we obtain the desired linear response relation between the macro-variables $\delta z(t)$ and $\delta f(t)$:

[TABLE]

The corresponding spring-dash-pot model is shown in Fig. 2. At the inner level we have a parallel bundle of a spring with stiffness $E_{1}=N\,(1-\beta\,c)/(\beta\,c)$ and a dash-pot with viscosity $\eta_{1}=N\,\tau_{x}/(\beta\,c)$. This subsystem accounts for the cross-bridge dynamics. The outer level contains a spring with stiffness $E=N\,\lambda_{f}/(1+\lambda_{f})$ and the dash-pot with viscosity $\eta=\nu$; the latter represents viscous response of the effective backbone.

In the real space the rheological relation (10) takes the form

[TABLE]

where $C=E_{1}\,E/(E_{1}+E)$ and $\theta=\eta_{1}/(E_{1}+E)$. If $\tau_{x}=0$ we obtain the Kelvin-Voigt model, and if $\tau_{z}=0$ Eq. 11 reproduces the rheological structure of the (passive) Hill’s model.

To make a more detailed comparison in a hard device we can define the storage modulus $G^{\prime}$ and the loss modulus $G^{\prime\prime}$ as real and imaginary parts of the ratio $\delta f(\omega)/\delta z(\omega)$; in a soft device we must consider instead the ratio $\delta z(\omega)/\delta f(\omega)$. The frequency dependence of these parameters is illustrated in Fig. 3. Note the divergence of $G^{\prime\prime}_{H}$ at large $\omega$ in qualitative difference with the Hill’s model where this parameter tends to zero. Similarly, in the Hill’s model $G^{\prime}_{S}$ has a finite limit at large $\omega$ while in our model it tends to zero.

Consider now the response of the system (11) to canonical step-like perturbations Mainardi and Spada (2011) as in the typical muscle experiments. In the hard device, the response to the input $\delta z(t)=\delta z_{0}\,\Theta(t)$, is described by the equation $\delta f(t)=\delta z_{0}\left[\,\left(E-C\right)\,e^{-\frac{t}{\theta}}+C+\nu\,\delta(t)\right]\mbox{,}$ were $\delta(t)$ is the Dirac delta function, whose effect cannot be detected in experiments unless the perturbation is strictly instantaneous. Note the first jump in tension $\lim_{t\rightarrow 0}\delta f(t)=E\,\delta z_{0}$, taking place simultaneously with the applied length step (phase 1 in Fig. 1b and 4a), which is the signature of a purely elastic behavior Reconditi et al. (2004). The elastic phase is followed by an exponential relaxation (phase 2 in Fig. 1b and 4a) with the time scale $\theta$. The condition $\theta>0$, or equivalently, $1+\lambda_{f}-\beta\,c>0$, then serves as a bound (since we are using an approximate expression for the correlation function $S_{xx}$) for the mechanical stability threshold of the equilibrium system in the hard device.

In the soft-device the response to a small step like perturbation $\delta f(t)=\delta f_{0}\,\Theta(t)$ is described by the equation

[TABLE]

where we introduced two new effective time scales

[TABLE]

Observe that according to (12), $\lim_{t\rightarrow 0}\delta z(t)=0$, indicating that there is no synchronous response to a step-like perturbation. The fact that in the soft device the relaxation starts when the perturbation is still being delivered Decostre et al. (2005), explains the difficulty in separating stages 1 and 2 (shown in Fig. 1a and 4b) of the transient response in isotonic conditions: since $\tau_{-}\ll\tau_{+}$ the first relaxation process with the time scale $\tau_{-}$ was sometimes interpreted as purely elastic Huxley and Simmons (1971); Piazzesi et al. (2003); Decostre et al. (2005); Pia (2007). The approximate stability condition is now $1-\beta\,c>0$. Note the difference between the stability thresholds in soft and hard device reflecting the ensemble inequivalence in this mean field system Barré et al. (2001); Caruel et al. (2013).

Using the basic physiological constraints on the parameters, $\beta\geq 0$, $\nu\geq 0$, $\lambda_{f}\geq 0$, one can show that in stable regimes $\tau_{+}>\theta$, which is in agreement with the fact that the relaxation in the soft the device is slower than in the hard device. In unstable regimes the expression under the square root in (13) may become negative which would indicate the possibility of oscillatory relaxation (in the soft device). Damped oscillations have been indeed observed in some mechanical experiments Edman and Curtin (2001); Sugi and Tsuchiya (1981), however, they appeared at larger timescales where the implied mechanical instability might have been suppressed actively Sheshka et al. (2016).

To make quantitative predictions and compare our results with physiological measurements we need to calibrate the model using the physical values of parameters.

For the cross-bridge stiffness we take the value $\kappa=3.29$ pN/nm Ford et al. (1977), while the combined stiffness of actin and myosin filaments can be estimated at the value $\kappa_{f}\simeq 153$ pN/nm Huxley et al. (1994); Wakabayashi et al. (1994). Since $N\simeq 100$ we obtain $\lambda_{f}=\kappa_{f}/(\kappa\,N)\simeq 0.46$. Given that the maximum working stroke is $a=11$ nm, the energy scale is $\kappa\,a^{2}=398\,\mbox{zJ}$, and we can conclude that at $T=273$ K the value of the dimensionless inverse temperature is $\beta\simeq 105$.

To estimate parameter $\gamma_{x}$ we assume that a typical myosin head has the diameter $h\simeq 6$ nm and that the cytoplasm has the effective dynamical viscosity $\eta\simeq 2.3\times 10^{-3}\,\mbox{Pa}\,\mbox{s}$ Kushmerick and Podolsky (1969); Arrio-Dupont et al. (1996). Viewing it as a sphere we obtain $\gamma_{x}\simeq 3\,\pi\,\eta\,h\simeq 1,30\times 10^{-4}$ ms pN/nm, so that the corresponding time unit is $\gamma_{x}/\kappa=3.44\times 10^{-5}$ ms. The dimensionless viscosity of a thick filament $\nu=\gamma_{z}/\gamma_{x}$ can be now calculated using the estimate for the drag coefficient $\gamma_{z}=\xi/\rho_{f}\simeq 0.480\,\mbox{ms pN}/\mbox{nm}$, where $\xi\simeq 3\times 10^{8}\,\mbox{N s}/\mbox{m}^{3}$ is a viscosity coefficient obtained in Ford et al. (1977) for Rana temporaria. We use here the expression $\rho_{f}=2/(\sqrt{3}\,d^{2})$ for the density of thick filaments in the section of a sarcomere where the thick filaments form a triangular pattern and lie at a distance $d=43$ nm Matsubara and Elliott (1972). The drag coefficient for the bundle of thick filaments constituting the half-sarcomere can be taken in the form $\gamma_{z}^{HS}\simeq\xi\,S$, where $S\simeq 8\,\mu\mbox{m}^{2}$ is the cross section of the sarcomere. Substituting this expression into the estimate for $\gamma_{z}$ we can extend the predictions of the model to the case of the whole half-sarcomere. The number of crossbridges $N$ has to be then multiplied by number of thick filaments in a half-sarcomere with the parallel rescaling $\delta f\rightarrow\delta f/(S\,\rho_{f})$ .

To model the double well potential $V(x)$ we use the simplest quartic function with one minimum $V(0)=0$ (pre-power-stroke state) and another minimum $V(-0.25)=-0.03$ (post-power-stroke state). The barrier is chosen at the level $V(x_{*})=0.11$ because the activation energy for the muscle power stroke was previously estimated to be either $55.3$ zJ or $95.1$ zJ Elangovan et al. (2012); Anson (1992) and we have chosen the smaller value as more relevant for the the fast time transient response.

The rheological equation (11) obtained for a single half-sarcomere must be now renormalized to the scale of a muscle fiber. To this end we assume that the response is affine, at least when perturbations are sufficiently small. We can view a myofibril as a chain of $L\sim 10^{4}$ half-sarcomeres connected in series, and represent a muscle fiber by a parallel arrangement of $M\sim 200-2000$ such myofibrils. The renormalization will then reduce to the substitutions $\delta z\rightarrow\delta z/L$ and $\delta f\rightarrow\delta f/M$.

Using these values of parameters we compute $\theta\simeq 10200$ ($\simeq 0.35$ ms) which is compatible with the relaxation time measured in Piazzesi et al. (2003). Next we estimate the elastic modulus for the instantaneous response of a half-sarcomere in the hard device. Given that $l_{0}=95.5$ ($\simeq 1.05\,\mu\mbox{m}$) is the length of a half-sarcomere and $r\simeq 0.83$ is the ratio of the cross-section of the muscle fiber occupied by sarcomeres Mobley and Eisenberg (1975), we obtain $E_{Y}=\lim_{t\rightarrow 0}\delta f(t)\,\rho_{f}\,r\,l_{0}/\delta z_{0}\simeq N\,\rho_{f}\,r\,\lambda_{f}\,l_{0}/(1+\lambda_{f})=190$ ($\simeq 57\,\mbox{MPa}$), close to the value measured in Piazzesi et al. (2003). Finally we compute $\tau_{+}\simeq 12700$ ($\simeq 0.44$ ms) and $\tau_{-}\simeq 116$ ($\simeq 0.004$ ms), which is also in good agreement with experimental observations Piazzesi et al. (2002); Decostre et al. (2005).

We are now in the position to evaluate to what extent the microscopic stochastic model and the macroscopic deterministic model can reproduce the outcomes of the realistic experiments. Numerical simulations of the microscopic model were conducted with a second order stochastic Runge-Kutta algorithm. We simulated the response of the microscopic model to a step-like perturbation $\delta z$ in a hard device, and $\delta f$ in a soft device, computing the corresponding responses $\delta f(t)$ and $\delta z(t)$.

The results are summarized in Fig. 4 where we used the physiological values of parameters except that in the simulation of the stochastic micro-model of a half sarcomere, instead of the realistic value $N\sim 500000$, we used the computationally reachable value $N=32768$, while appropriately rescaling the parameter $\nu$ to ensure that the time scale $\tau_{z}$ remains at its realistic value. Numerical experiments aimed at a single bundle of thick and thin filaments with $N=128$ and using the physiological value of $\nu$ give practically the same results. In both cases the agreement between the stochastic model and the rheological equation  (11) is excellent for both hard and soft devices.

To conclude, we have shown that the time dependent passive viscoelastic behavior of striated muscles can be understood at the quantitative level starting from the microscopic structure of a half sarcomere. From a stochastic microscale model we derived a deterministic rheological model which describes linear response of muscle fibers under various loading conditions and found explicit relations between the macroscopic and the microscopic parameters. The fact that the derived rheological relation involves not only first but also second derivatives allowed us to explain the qualitative differences in the mechanical response of isometrically and isotonically loaded muscles. The model was calibrated based on independent data and excellent quantitative agreement was reached with physiological observations. It would be of interest to check experimentally the predicted link between the the mechanical response to fast perturbationsPiazzesi et al. (2002, 2003) and the power spectrum of mechanical fluctuations Malmerberg et al. (2015).

Acknowledgements. The authors thank M. Caruel, R. Garcia Garcia and S. Ruffo for helpful discussions. FS was supported by a postdoctoral fellowship from Ecole Polytechnique. The work of LT was supported by the grant ANR-10-IDEX-0001-02 PSL.