# Hierarchical modeling of force generation in cardiac muscle

François Kimmig, Matthieu Caruel

► To cite this version:

François Kimmig, Matthieu Caruel. Hierarchical modeling of force generation in cardiac muscle. Biomechanics and Modeling in Mechanobiology, 2020, 19, pp.2567-2601. 10.1007/s10237-020-01357-w. hal-02570784

HAL Id: hal-02570784

https://hal.science/hal-02570784

Submitted on 12 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Hierarchical modeling of force generation in cardiac muscle

François Kimmig¹,² and Matthieu Caruel³

¹LMS, CNRS, École polytechnique, Institut Polytechnique de Paris, France

²Inria, France

³MSME, CNRS, Université Paris-Est, France

## Abstract

Performing physiologically relevant simulations of the beating heart in clinical context requires to develop detailed models of the microscale force generation process. These models however may reveal difficult to implement in practice due to their high computational costs and complex calibration. We propose a hierarchy of three interconnected muscle contraction models – from the more refined to the more simplified – that are rigorously and systematically related with each other, offering a way to select, for a specific application, the model that yields a good trade-off between physiological fidelity, computational cost and calibration complexity. The three models families are compared to the same set of experimental data to systematically assess what physiological indicators can be reproduced or not and how these indicators constrain the model parameters. Finally, we discuss the applicability of these models for heart simulation.

**Keywords** – muscle modeling; model reduction; sliding filaments; cross-bridges; sarcomere

## 1 Introduction

Mechanical modeling of the micro-scale muscle contraction mechanisms is an essential component for patient-specific physiologically relevant in silico heart simulations, with the aim of providing effective diagnoses and treatment planning tools. The effectiveness of the developed models relies primarily on their ability to reproduce the biological processes at the origin of the muscle contraction with a level of detail adapted to the investigated clinical questions.

One of the most essential process to be modeled is the conversion of the metabolic energy extracted from ATP turnover into mechanical work by Myosin II molecular motors in the presence of actin [Hill, 1938; Alberts, 2015; Barclay et al., 2010; Barclay, 2015].

In the muscle tissue, this myosin-actin interaction occurs inside elementary contractile units called sarcomeres, where the myosin motors bundle into thick filaments. The heads of the myosin motors, protruding from the thick filaments, can attach to specific actin binding sites on the neighboring parallel actin (thin) filaments. This attachment is possible only after the actin filaments have been activated by calcium ions whose release in the cell cytosol triggers the contraction.

While attached to actin, a myosin head produces force via a large conformational change – the power stroke – which induces a relative displacement between the myosin filaments and the actin filaments (~ 10 nm at zero load) Lynn & Taylor [1971]; Huxley & Simmons [1971]; Rayment et al. [1993b,a]. The energy necessary to recharge the power stroke mechanism is provided by ATP hydrolysis inside the catalytic domain of the myosin protein, while the motor is detached from actin.

The concerted action of antagonistically oriented groups of molecular motors, pulling the thin filaments in opposite directions, allows for the shortening of the sarcomeres. These force generating units are arranged in series along the longitudinal direction of the muscle fibers which allows to generate and transmit force throughout the whole tissue.

The active force production is subjected to two types of regulation mechanisms that vary the level of the developed force. The first regulation mechanism is extrinsic, via the neuroendocrine system which affects various levels of the activation-contraction coupling. It results in particular in variations of the calcium supply inside the cell [Silverthorn et al., 2009] and of the thin filament responsiveness to
calcium [Solaro & Rarick, 1998], and it may also modify the actin-myosin interaction [de Tombe, 2003]. The second regulation mechanism is intrinsic, via the modulation of the actin filament responsiveness to calcium, ultimately affecting the level of thin filament activation, and the availability of the myosin heads as a function of the degree of sarcomere stretch [Ter Keurs et al., 1980; Dobesh et al., 2002; de Tombe et al., 2010].

Both regulation pathways are essential for the physiological heart contraction and the intrinsic regulation is, in particular, at the core of the so-called Frank-Starling effect at the macroscopic scale [Allen & Kentish, 1985; de Tombe et al., 2010].

The extrinsic regulation acts by means of complex inter-protein interactions, whose study is beyond the scope of this work, and we thus focus on the basal behavior of muscle cells, i.e. when isolated from the neuroendocrine system, as is the case in ex-vivo preparations. Furthermore, recent experimental studies performed on ex-vivo preparations of different types have demonstrated that the intrinsic regulation mechanisms do not alter the cycling kinetics [Amiad Pavlov & Landesberg, 2016], the response to fast load changes and the steady-state force-velocity relation [Caremani et al., 2016] nor the force per attached head [Pinzauti et al., 2018]. All these indicators being directly linked to the actin-myosin interaction, we conclude that the intrinsic regulation does not affect the internal mechanical properties of the actomyosin system, and only modulates the number of motors participating in the contraction. Hence, the basal force generation mechanism can be studied, both experimentally and theoretically, at any nominal sarcomere length and any level of calcium activation, by appropriately normalizing the measured force.

The most widely used experiment designed to characterize the mechanical output of the contractile system consists in measuring, at the sarcomere level, the transient isotonic shortening in response to a sudden force step applied within ~200 µs from an isometric state [Caremani et al., 2016]. Remarkably, the reaction to this perturbation allows to identify, from a single experiment, three essential physical characteristics of the actomyosin system – namely elasticity, power stroke and attachment and detachment processes – by using the fact that they operate at three different timescales: ~200 µs, ~1 ms and ~30 ms, respectively [Caremani et al., 2016].

The data obtained from these experiments have served as a benchmark for the theoretical modeling of the actomyosin system since the seminal work of A.F. Huxley [Huxley, 1957], later thermodynamically formalized by T.L. Hill [Hill, 1977]. In this framework, myosin heads are modeled as a spring that can reversibly attach to specific binding sites on the surrounding actin filaments [Eisenberg et al., 1980; Piazzesi & Lombardi, 1995; Smith et al., 2008; Smith & Mijailovich, 2008; Månsson, 2010; Caremani et al., 2015]. The spring itself has internal – discrete or continuous – degrees of freedom that represent the conformational – and ligand-binding – state of the protein. The actin and myosin filaments are then usually modeled as an independent spring in series with the molecular motors [Linari et al., 1998], which then effectively act on two rigid backbones [Ford et al., 1981].

A large variety of models of this type, though with different degrees of complexity, can be derived to account for a large body of experimental data. The choice of a particular theoretical model of the muscle tissue usually results from a compromise between the degree of model refinement – to reproduce the micro-scale physiology –, the possibility to calibrate the model with the available data, and the observability of the modeling assumptions effects on the macro-scale observables. Moreover, the use of refined representations of the actomyosin interaction in realistic models of the heart contraction comes at a high computational cost, which may limit the applicability of numerical simulation in a clinical context where real-time simulations are often required, indeed.

In this paper, we propose a hierarchical approach to this problem by formulating a series of interrelated models, from the more refined, capturing the finest effects observed in single cell experiments, to the more coarse grained aimed at efficient organ simulations. Each coarse graining step is obtained from simplifications of the more refined models, such that the loss of information is well controlled. For each model in this hierarchy, we discuss its ability to reproduce the experimental data and its relevance for heart simulations as a guidance for appropriate model choosing depending on the application sought.

Our starting point is the stochastic model formulated in [Caruel et al., 2019], which combines a standard Hill-type approach to describe the attachment and detachment processes with a more recently developed continuous stochastic dynamics for the internal variables describing the power stroke conformational change [Marcucci & Truskinovsky, 2010; Caruel et al., 2013; Marcucci et al., 2016].

The first level of simplification is obtained by adiabatically eliminating these internal variables based on the observation that the ~1 ms timescale of power stroke can be considered to be low compared to the ~30 ms timescale associated with the completion of the ATPase cycle [Caremani et al., 2015, 2016].
Hence, at each time one can make the approximation that internal degrees of freedom are in thermal
equilibrium. This procedure, already used in [Zahalak, 2000; Hill, 1977; Caruel et al., 2019], results in
a power stroke equilibrated model (PSE model), which can be viewed as an instance of the classical
Huxley'57 two-state (attached and detached) model.

Both the stochastic and the PSE models describe the evolution of the population of the available myosin motors – either from stochastic differential equations (SDE) or deterministic partial differential equations (PDE) – from which the mechanical macroscopic outputs such as force, stiffness or thermodynamic yield, can be derived by direct statistical averaging.

The second level of simplification is based on the classical expansion of the PDE solution of the PSE model into an infinite series of macroscopic moments, whose dynamics are obtained by solving an infinite set of coupled Ordinary Differential Equations (ODEs) [Chapelle et al., 2012; Bestel et al., 2001; Zahalak, 1981]. Simple closure relations can be formulated in order to limit the system at an arbitrary – ideally low – number of equations. In this case, the macroscopic output of the models appears as linear functions of the macroscopic moments.

Our first result is that, by restricting the space of functions describing the energy landscapes of the conformational variables and the strain dependence of the transition rates, the constraints imposed by the available experimental data define the set of parameters almost uniquely for both the stochastic and the PSE models [Caremani et al., 2016]. We show that, once calibrated, the stochastic model accurately reproduces the mechanical response to fast load clamp experiments. Interestingly, while the derived PSE model fails, by construction, to reproduce the physiological processes associated with the shortest timescales, it is still able to reproduce the most relevant markers of the mechanical performance of cardiac fibers [Caremani et al., 2016; Pertici et al., 2018; de Tombe & Stienen, 2007]. We emphasize that with both models the ATP consumption is not calibrated but correctly predicted, showing that the calibration procedure is robust.

The macroscopic models are calibrated by further restraining the PSE model space of functions for the attachment-detachment process and for the internal equilibrium energy landscape. The second result of our work is that the obtained macroscopic models remain in good agreement with the PSE model and experimental data over the range of loading conditions that is relevant for the heart functioning during a typical pressure-volume loop, the differences becoming significant only in regimes that are non physiological.

In conclusion, our stochastic model appears as a good candidate for representing the refined behavior of the actomyosin system with the aim at understanding, for instance, the impact of a mutation or a drug on its mechanical properties. Alteration of this type are typically investigated with tools such as structural crystallography [Robert-Paganin et al., 2019; Kuhlman & Bradley, 2019; Robert-Paganin et al., 2018], or single molecule experiments [Woody et al., 2018], whose results could indeed be reproduced using our stochastic model. However, the computational costs associated with the stochastic model makes it likely inadequate for organ scale finite elements simulation. The PSE model represents a significant gain in terms of computational complexity while still being able to capture the main physiological characteristics of the actomyosin system, therefore providing a good tradeoff between physiological relevance and computational efficacy. The fastest simulations can be performed with the ODE-based macroscopic models at the expanse of only a partial adequacy with the physiological functioning. Hence, despite their simplicity and the fact that they are not able to fully reproduce the behavior of the contractile system observed in vitro, the macroscopic models appear to be particularly fit for 3D heart simulation in clinical context [Sermesant et al., 2012; Chabiniok et al., 2011]. The originality of our approach is that the well established connexion between the different levels of our hierarchy allows to consistently enrich the simpler models if need be.

The paper is organized as follows. In Section 2 we derive our hierarchy of models, starting with the population models. Section 3 presents the calibration and the comparison of the models using the isotonic shortening experimental data. Finally, Sections 4 and 5 present a discussion of the results and our conclusions, respectively.
<figure><img src="image_2.png" /><figcaption>Figure 1: Summary of the stochastic model (a) and its reduction to a two state model PSE model (b). In (a), the states are indexed as follows: 1 detached pre-power stroke; 2 attached pre-power stroke, 3 attached post-power stroke, 4 detached post-power stroke. The PSE model is obtained by the adiabatic elimination of the stochastic model power stroke degrees of freedom X and Y.</figcaption></figure>

# 2 Hierarchy of models

## 2.1 Population models

### 2.1.1 Stochastic model

The model on which the proposed hierarchy is built on was formulated in [Caruel et al., 2019], and calibrated using data from frog skeletal muscle experiments. The main purpose of this model is to capture all the response timescales of a fiber submitted to rapid load changes. The framework was formulated by Huxley [1957] and Hill [1977, 1974, 1976], based on the hypothesis that the myosin heads and actin sites are regularly distributed along their respective filaments, and that the distance to their neighbors is constant. The compliance of the myofilaments is assumed to be lumped into an effective linear elastic spring in series with a set of myosin motors interacting with two rigid backbones [Ford et al., 1981]. We consider a single actin site model, in which the myosin head can bind only to the nearest actin site, located at a distance s from the anchor point of the head in the myosin filament, see Figure 1(a).

The state of the myosin head is parametrized by three stochastic variables :

* $X^t$, a continuous variable representing the displacement between the position of the head tip in the current state and at rest;

* $Y^t$, a continuous internal representing the internal conformation of the head (power stroke)

* $\alpha^t$, a discrete variable accounting for the attachment state of the head, taking the value $\alpha^t = 1$ when the myosin head is attached and the value $\alpha^t = 0$ when it is detached.

Following [Marcucci & Truskinovsky, 2010], the myosin head is modeled as a bistable snap-spring comprising an elastic component, in series with a bi-stable element accounting for the power stroke conformational change, see Figure 1(a). The internal energy of the head in state $\alpha$ is then defined as

$$w(x, y, \alpha) = w_{\alpha}(x, y) = \frac{\kappa}{2}(x + y)^2 + u_{\alpha}(y), \quad (1)$$

where $\kappa$ is the stiffness of the elastic component and $u_\alpha$ is a double well potential associated with the bistable element. The force developed by the attached myosin head on the thin filament is thus $\tau_c(x, y) = \partial_x w_1(x, y) = \kappa(x + y)$. Hence, an increase of y - representing the power stroke - effectively stretches the elastic element and increases the developed force, see Figure 1(a). When the myosin is detached, no force is developed on average.
In the overdamped regime, the dynamics of an individual myosin head is governed by the following system of stochastic equations [Caruel et al., 2019]

$$ \begin{cases} dX^t = [\alpha^t \dot{x}_c - (1 - \alpha^t)\eta^{-1} \partial_x w_\alpha(X^t, Y^t)] dt + \delta(t - t_s)[s(t) - X^t] dt + (1 - \alpha^t)\sqrt{2D} dB_x^t, & (2a) \\ dY^t = -\eta^{-1} \partial_y w_\alpha(X^t, Y^t) dt + \sqrt{2D} dB_y^t, & (2b) \end{cases} $$

$$ \mathbb{P}[\alpha^{t+dt} = 1 | \alpha^t = 0] = k_+(Y^t, s)dt, \quad (2c) $$

$$ \mathbb{P}[\alpha^{t+dt} = 0 | \alpha^t = 1] = k_-(Y^t, s)dt, \quad (2d) $$

where $\eta$ is a drag coefficient, $\dot{x}_c$ is the sliding velocity (being positive for an extending sarcomere). The parameter $D = k_B T / \eta$ where $T$ and $k_B$ denote the temperature and the Boltzmann constant, respectively, is a diffusion coefficient scaling the variance of the two standard brownian motion increments $dB_{x,y}^t$. In (2a), $t_s$ denotes the time when a transition from the detached to the attached state occurs and the term $\delta(t - t_s)[s(t) - X^t]dt$ ensures that between $t_s$ and $t_s + dt$, the variable $X^t$ jumps from its current value $X^{t_s}$ to the value $s(t_s + dt)$. Equations (2c) and (2d) represent the probabilities to switch between the attached and the detached states in the time interval $[t, t + dt]$. The associated attachment and detachment transition rates are denoted by $k_+$ and $k_-$, respectively.

We now consider a population of heads in a segment of a cardiomyocyte that has the thickness of a half-sarcomere. Since the density of myosin heads in such a segment is $1.25 \times 10^{17} \text{ m}^{-2}$ [Pinzauti et al., 2018], a statistical description of the system can be used. We denote by $p(x, y, \alpha; s, t)$ the probability distribution of the myosin head state for a subpopulation of the heads located at distance $s$ of the nearest actin site at time $t$. This probability distribution mixing continuous and discrete variables is normalized in the following way:

$$ \sum_{\alpha=\{0,1\}} \iint p(x, y, \alpha; s, t) dx dy = 1, \quad \forall t \ \forall s. $$

Since, in the attached state, the tip of the myosin head can only be located at the position of the actin site ($X^t = s(t)$), the probability distribution $p(x, y, \alpha = 1; s, t)$ is degenerated and becomes

$$ p(x, y, \alpha = 1; s, t) = \delta(x - s) \bar{p}(y; s, t). $$

The Fokker-Planck equations associated to the system of SDEs (2) are (see [Caruel et al., 2019])

$$ \left\{ \begin{aligned} \partial_t p(x, y, 0; s, t) + \dot{x}_c \partial_s p(x, y, 0; s, t) &= \partial_y \left[ \eta^{-1} \partial_y w_0(x, y) p(x, y, 0; s, t) + D \partial_y p(x, y, 0; s, t) \right] \\ &\quad + \partial_x \left[ \eta^{-1} \partial_x w_0(x, y) p(x, y, 0; s, t) + D \partial_x p(x, y, 0; s, t) \right] \\ &\quad + k_-(y, s) \delta(x - s) \bar{p}(y; s, t) - k_+(y, s) p(x, y, 0; s, t), \quad (3) \\ \\ \partial_t \bar{p}(y; s, t) + \dot{x}_c \partial_s \bar{p}(y; s, t) &= \partial_y \left[ \eta^{-1} \partial_y w_1(x, y) \bar{p}(y; s, t) + D \partial_y \bar{p}(y; s, t) \right] \\ &\quad + \int k_+(y, s) p(x, y, 0; s, t) dx - k_-(y, s) \bar{p}(y; s, t), \end{aligned} \right. $$

and their solution allows to compute the active force as the average force generated by the attached head

$$ \tau_c(t) = \frac{1}{d_a} \int_{s^-}^{s^+} \int_{-\infty}^{+\infty} \partial_s w_1(s, y) \bar{p}(y; s, t) dy ds, \quad (4) $$

where $d_a$ is the distance between two consecutive actin sites, and $[s^-, s^+]$ is the interval of reachable actin sites with $s^+ - s^- = d_a$. The cross-bridges being arranged in parallel in the proposed half-sarcomere model, the macroscopic tension developed by a cardiomyocyte is obtained by rescaling the active force (4):

$$ T_c(t) = \rho_{\text{surf}} \tau_c(t), \quad (5) $$
where $\rho_{surf}$ is the number of myosin heads in a longitudinal portion of a cardiomyocyte that has a thickness equal to the half-sarcomere length $\ell_{hs}$ per cross-section area, both quantities being defined in the reference configuration.

As noted by T.L. Hill [Hill, 1977], the assumption of a single actin site available at a given time for each myosin head imposes that the probability that a myosin head is attached on the boundary of the interval $[s^{-}, s^{+}]$ vanishes. We should thus have

$$\iint p(x, y, \alpha; s^{-}, t) dxdy = \iint p(x, y, \alpha; s^{+}, t) dxdy = 0, \text{ for } \alpha \in \{0, 1\}. \quad (6)$$

This property can be ensured by an adequate choice of the attachment and detachment rates $k_{+}$ and $k_{-}$.

We remind that this model does not take into account the calcium-induced thin filament activation [Kobayashi et al., 2008]. The recently observed OFF-state in which the myosin head is folded on the back bone of the thick filament and do not undergo the Lymn-Taylor cycle [Linari et al., 2015], [Reconditi et al., 2017] is also not considered. Hence, all myosin heads are considered available for attachment and all actin sites are considered activated.

### 2.1.2 Adiabatic elimination of the power stroke: the PSE model

The first level of simplification of the stochastic model is obtained by considering the experimental evidence of a separation between the timescale of the power stroke transition (~1 ms) and the characteristic timescale for the completion of the Lymn-Taylor cycle (~30 ms), observed in the muscle response to a rapid change in loading conditions (see Section 3.1.1).

Building on this separation of timescales, one can adiabatically eliminate the internal degrees of freedom $X$ and $Y$ – characterizing the power stroke conformational change – by assuming that their dynamics is infinitely fast compared to the timescale associated with the other variable dynamics ($\alpha^{t}$). Here we show that this elimination results in a two state Huxley’57 type model. In the limit of an infinitely fast relaxation, the distribution of the random variables $X$ and $Y$ follows the classical Boltzmann equilibrium distribution

$$p_{0}^{th}(x, y) = \frac{\exp(-w_{0}(x, y)/(k_{B}T))}{\iint \exp(-w_{0}(x, y)/(k_{B}T)) dy dx}, \quad (7a)$$

$$p_{1}^{th}(y; s) = \frac{\exp(-w_{1}(s, y)/(k_{B}T))}{\int \exp(-w_{1}(s, y)/(k_{B}T)) dy}. \quad (7b)$$

The probability distributions $p(x, y, \alpha; s, t)$ can then be decomposed as follows

$$\begin{cases} p(x, y, 0; s, t) = P_{0}(s, t)p_{0}^{th}(x, y), \\ p(x, y, 1; s, t) = P_{1}(s, t)\delta(s - x)p_{1}^{th}(y; s), \end{cases}$$

where

$$P_{1}(s, t) = \iint p(x, y, 1; s, t) dxdy \quad \text{and} \quad P_{0}(s, t) = \iint p(x, y, 0; s, t) dxdy = 1 - P_{1}(s, t)$$

define the probabilities that a given myosin head is attached and detached, respectively, or equivalently the ratio of attached and detached heads among the population. They are named hereafter population probability of being attached and detached. Remark that the ratio $P_{1}$ is denoted by $n$ in [Huxley, 1957]. The ensuing model has only two states corresponding to $\alpha = \{1, 0\}$, and characterized by the equilibrium free energies

$$\mathcal{F}_{0}^{th} = \iint \left[ w_{0}(x, y)p_{0}^{th}(x, y) + k_{B}T p_{0}^{th}(x, y) \ln(a^{2} p_{0}^{th}(x, y)) \right] dxdy, \quad (8a)$$

$$\mathcal{F}_{1}^{th}(s) = \int \left[ w_{1}(s, y)p_{1}^{th}(y; s) + k_{B}T p_{1}^{th}(y; s) \ln(a p_{1}^{th}(y; s)) \right] dy, \quad (8b)$$
where the parameter a denotes the characteristic size of the power stroke, used as the reference length in our system. The free energies $\mathcal{F}_{0,1}^{th}$ appear naturally as the averages of the equilibrium chemical potentials (internal energy plus entropic terms) defined by the more refined stochastic model, while they are postulated in the classical Huxley'57 type models. The two functions are represented schematically in Fig. 1(b) and shown for the calibrated model in Fig. 15

To derive the simplified dynamics, we integrate the Fokker-Planck equation (3) with respect to $x$ and $y$. Then, using the explicit definition of the internal energy levels (1) and the equilibrium distributions (7) and obtain

$$ \partial_t P_1(s,t) + \dot{x}_c \partial_s P_1(s,t) = f^{th}(s)(1 - P_1(s,t)) - g^{th}(s)P_1(s,t), \quad (9) $$

$$ P_0(s,t) = 1 - P_1(s,t) \quad (10) $$

where (9) corresponds to the classical Huxley'57 model conservation equation. Importantly again, as opposed to the Huxley'57-type model, [Huxley, 1957] the averaged transition rates

$$ f^{th}(s) = \iint k_+(y, s) p_0^{th}(x, y) dy dx, \quad (11a) $$

$$ g^{th}(s) = \int k_-(y, s) p_1^{th}(y; s) dy. \quad (11b) $$

are not defined per se but are instead derived from the rates $k_+$ and $k_-$ of the stochastic model.

The average force per head generated by a population of myosin heads (4) now simplifies into

$$ \tau_c^{th}(t) = \frac{1}{d_a} \int_{s^{-}}^{s^{+}} P_1(s,t) \tau^{th}(s) ds, $$

where the term $\tau_c^{th}(s)$ is the average equilibrium tension exerted by an attached head located at a distance $s$ from its binding site, which derives from the thermal equilibrium free energy in the attached state $\mathcal{F}_1^{th}$:

$$ \tau^{th}(s) = \frac{d\mathcal{F}_1^{th}}{ds}(s). \quad (12) $$

With the assumption of a linear neck elasticity used in our model (see (1)), we simply have

$$ \tau^{th}(s) = \int \kappa(y + s) p_1^{th}(y; s) dy. \quad (13) $$

The equivalence between (12) and (13) is given in Appendix (A.1), where the average stiffness per head is also computed.

The boundary conditions requirement that no head can remain attached outside the $[s^{-}, s^{+}]$ interval (see (6)) becomes here

$$ P_1(s^{-}, t) = P(s^{+}, t) = 0, \quad (14) $$

which is ensured by appropriate choices of the transition rates $f^{th}$ and $g^{th}$ or their antecedent $k_+$ and $k_-$. A more detailed assessment of the mathematical properties of the Huxley'57 equation solutions is presented in [Kimmig et al., 2019].

Finally, similarly to (5), the tension developed by a muscle fiber is

$$ T_c(t) = \rho_{\text{surf}} \tau_c^{th}(t). \quad (15) $$

To conclude this section we mention that the PSE model can be seen as an instance of a larger family of models within the Huxley'57 framework (see Appendix A.2 for the details). These models are characterized by only two states (detached and attached corresponding respectively to $\alpha = 0, 1$) associated with the free energies $\mathcal{F}_{0,1}$ and the transition rates $f$ and $g$, which are the direct inputs of the dynamical system (9) and the mechanical force (12). Furthermore the approach developed here to eliminate the fast equilibrating continuous internal degrees of freedom $X$ and $Y$ can be directly transposed to the more conventional modeling approach using a discrete internal degree of freedom to characterize the various states of the cross-bridge Huxley & Simmons [1971]; Eisenberg et al. [1980].
## 2.2 Macroscopic models

A further reduction of the model from Huxley's 57 equations has been proposed in [Bestel et al., 2001; Chapelle et al., 2012] and reused in [Caruel et al., 2019]. It aims at establishing the dynamics equations associated directly with relevant macroscopic quantities such as the tension developed by the population of myosin heads in a half-sarcomere. We apply this simplification here and extend it to non-linear constitutive behaviors. For a choice of the attached state free energy as a polynomial function, the active tension is defined as a combination of the population probability moments formally defined by $\int s^p P_1(s) ds$ (see (29)), which motivates the decomposition of the solution on these quantities.

We consider in all generality a multi-site model, i.e. the myosin head can bind to any actin sites, which are regularly located along the thin filament at distances $s + jd_a \forall j \in \mathbb{Z}$ (see [Kimmig et al., 2019] for more detail). The moment of order p of the population probability $P_1$ is defined by

$$M_p(t) = \frac{1}{d_a} \int_{s-}^{s^+} \sum_{j \in \mathbb{Z}} (s + jd_a)^p P_1(s + jd_a, t) ds.$$

Integrating (9) and performing an integration by parts – with the boundary conditions (14) –, we obtain the following unclosed system of ODEs, for $p \ge 1$

$$\dot{M}_p(t) = p \dot{x}_c M_{p-1}(t) + f_p - \frac{1}{d_a} \int_{s-}^{s^+} \sum_{j \in \mathbb{Z}} (f(s + jd_a) + g(s + jd_a))(s + jd_a)^p P_1(s + jd_a, t) ds, \quad (16)$$

with

$$f_p = \frac{1}{d_a} \int_{s-}^{s^+} \sum_{j \in \mathbb{Z}} (s + jd_a)^p f(s + jd_a) ds.$$

Solving the infinite system of coupled ODEs (16) is equivalent to solving the PDE (9). However, a significant gain in computational efficiency can obtained by reducing the number of relevant ODEs to a finite set through appropriate closure relations. For instance, one can assume that the sum of the transition rates $f(s) + g(s)$ does not depend on the space variable $s$ [Bestel et al., 2001]. In that case (16) simplifies into

$$\left\{ \begin{array}{l} \dot{M}_0(t) = -(f+g)M_0(t) + f_0, \\ \dot{M}_p(t) = p\dot{x}_c M_{p-1}(t) - (f+g)M_p(t) + f_p \quad \text{for } p > 0, \end{array} \right. \qquad (17)$$

where the dynamics of each moment only depends on the moments of inferior order, restricting the computation of any moment to the resolution of a finite number of ODEs.

### 2.2.1 Linear macroscopic models

Assuming, as in [Bestel et al., 2001] that the cross-bridges has a quadratic energy of the form

$$w_1^L(s) = \frac{\kappa_{\mathrm{xb}}}{2} (s + s_0)^2, \qquad (18)$$

where $s_0$ is the length of the unloaded spring, the macroscopic first Piola-Kirchhoff stress developed by the muscle fiber is then defined by

$$T_c(t) = \frac{\rho_{\mathrm{surf}}}{d_a} \int_{s-}^{s^+} \sum_{j \in \mathbb{Z}} \frac{d}{ds} w_1^L(s + jd_a) P_1(s + jd_a, t) ds \qquad (19)$$

which, given (18), gives $T_c(t) = \rho_{\mathrm{surf}}\kappa_{\mathrm{xb}}(s_0M_0(t) + M_1(t))$. Similarly, the macroscopic stiffness per unit surface is similarly given by

$$K_c(t) = \frac{\rho_{\mathrm{surf}}}{d_a} \int_{s-}^{s^+} \sum_{j \in \mathbb{Z}} \kappa_{\mathrm{xb}} P_1(s + jd_a, t) ds = \rho_{\mathrm{surf}}\kappa_{\mathrm{xb}} M_0(t).$$
The moment dynamics (17) then reduces to two ODEs describing the dynamics of $T_c$ and $K_c$, namely

$$ \begin{cases} \dot{K}_c(t) = -(f+g)K_c(t) + f_0 K_\infty, \\ \dot{T}_c(t) = -(f+g)T_c(t) + \dot{x}_c K_c(t) + (s_0 f_0 + f_1) K_\infty, \end{cases} \qquad (20) $$

with $K_\infty = \rho_{surf}\kappa_{xb}$ the stiffness in the rigor state where all cross-bridges are attached.

The relation between the force and the sliding velocity in steady-state shortening is classically measured by experimentalists and is usually termed *force-velocity relation*. With (20), the closed form expression of the steady state force-velocity relation is linear

$$ T_c^\infty = \frac{f_0 K_\infty}{(f+g)^2} \dot{x}_c + \frac{s_0 f_0 + f_1}{f+g} K_\infty. \qquad (21) $$

and reveals the close interplay between the cross-bridges stiffness (through the parameter $K_\infty$) and the myosin cycling rate (through $f$ and $g$) in the mechanical behavior. While the relation (21) is clearly incompatible with the observed hyperbolic shape of the force-velocity relation, we show below that it may still prove useful in some physiological conditions. Note that a non-linear force-velocity relation can still be obtained with a linear elasticity assumption [Huxley, 1957; Duke, 1999]. The linearity of the force-velocity relation predicted by this model results from the combined effect of the linear elasticity assumption and the specific choice of the transition rates.

### 2.2.2 Nonlinear macroscopic model

If, instead of the linear spring energy (18), we consider a nonlinear spring associated with a polynomial internal energy of the form

$$ w_1^{NL}(s) = \sum_{i=0}^{n} \frac{c_i}{\max(i, 1)} s^i, $$

then, adapting (19), the force can be written

$$ T_c(t) = \sum_{i=0}^{n-1} C_{i+1} M_i(t) \qquad (22) $$

where $C_i = \rho_{surf}c_i$. We show in appendix A.3 that considering the moment dynamics (17) in the permanent regime, (22) leads to the steady-state tension

$$ T_c^\infty = \frac{1}{f+g} \left[ \dot{x}_c \sum_{i=1}^{n-1} i C_{i+1} M_{i-1}^\infty(\dot{x}_c) + \sum_{i=0}^{n-1} C_{i+1} f_i \right], \qquad (23) $$

where $M_{i-1}^\infty$ is explicitly given by

$$ M_i^\infty(\dot{x}_c) = \sum_{j=0}^{i} \frac{\frac{i!}{j!} f_i}{(f+g)^{i+1-j}} \dot{x}_c^{i-j}. $$

The obtained force-velocity relation is then a polynomial function of the shortening velocity of order $n-1$ when the internal energy is of order $n$.

## 3 Results

In the previous section, we have derived two classes of models: the population models containing our most refined stochastic model (see Section 2.1.1) and the PSE model (see Section 2.1.2), and the macro-models (see Section 2.2). In what follows, we will first present the calibration and the results for the population models and compare the stochastic model to the PSE model; then we will present the calibration of the macro-models and study how they compare to the PSE model, the latter being considered at the pivot of the whole model hierarchy.
<table><caption>Table 1: Parameters whose value is not adjusted in the calibration process but directly taken from the literature. We present here typical values for our modeling purposes but it should be noted that uncertainties remain in the literature.</caption><thead><tr><th>Non calibrated parameters</th><th>symbol</th><th>Value</th><th>Reference</th></tr></thead><tbody><tr><td>Power stroke characteristic length</td><td>a</td><td>11 nm</td><td>[Holmes & Geeves, 2000]</td></tr><tr><td>Distance between two actin sites</td><td>d<sub>a</sub></td><td>40 nm</td><td>[Craig & Padrón, 2004]</td></tr><tr><td>Density of thick filament in the cross-section of thickness l<sub>hs</sub></td><td>ρ<sub>surf</sub></td><td>1.25 × 10<sup>17</sup> m<sup>-2</sup></td><td>[Pinzauti et al., 2018]</td></tr><tr><td>Myofilament compliance</td><td>C<sub>m</sub></td><td>17 nm MPa<sup>-1</sup></td><td>[Pinzauti et al., 2018]</td></tr></tbody></table>

## 3.1 Population models

### 3.1.1 Calibration

The properties of the system depend on three categories of model parameters. The first category characterizes the energy landscapes ($w_α$ for the stochastic model, see (1), and $\mathcal{F}_α^{th}$ for the PSE model, see (8)) and the viscosities associated with the internal degrees of freedom $X$ and $Y$. The second category determines the kinetics of the attachment and detachment interactions between the myosin heads and the actin filaments ($k_+$ and $k_-$ for the stochastic model or $f^{th}$ and $g^{th}$ for the PSE model). The last category concerns the energetic properties of the system related the metabolic reactions among which the most important one is the ATP turnover.

To limit the number of free parameters in the model and make the calibration process more robust, we impose the following restraints on the space of parameters. First, we assume a bi-quadratic form for the energy landscapes $u_α$ which, in addition to being analytically simple, have proven efficient in reproducing the fast transient response of skeletal muscle fibers, see [Marcucci & Truskinovsky, 2010; Caruel et al., 2013]. Second, the attachment and detachment transitions kinetics are described by (minimally) regularized piece-wise constant functions over the interior of the $[s^-, s^+]$ interval. Near the boundaries, the transitions rates either tend to zero or grow to infinity, exponentially.

The models calibration procedure uses the hierarchical link between the stochastic and the PSE model, which are thus calibrated in a interdependent and simultaneous manner. The parameters resulting from this calibration strategy are listed in Table 5.

**Fixed parameters** Before turning to the calibration procedure per se, we mention that some of the model parameters are not chosen through a calibration process but are directly taken from published experimental data. These parameters are listed in Table 1.

**Isometric indicators** We here refer to isometric conditions when the developed force reaches its peak in a twitch contraction in which the length of the fiber is maintained constant [Ter Keurs et al., 1980; Van Heuningen et al., 1982; Caremani et al., 2016]. It is not, as for skeletal muscles, a tetanised state.

The main indicator of the isometric state is the maximal isometric tension (or stress) $T_0$ generated by the fraction $\tilde{n}_{att}$ of the myosin motors bound to actin. In maximal activation conditions, this tension can reach the value of 118 kPa corresponding to a maximal ratio of attached heads $\tilde{n}_{att} = 0.15$ and to a tension per head of 6.14 pN [Caremani et al., 2016; Pinzauti et al., 2018].

For the PSE model, the isometric population probability $\tilde{P}_1$ is the steady state solution of (9) with $\dot{x}_c = 0$,

$$ \tilde{P}_1(s) = \frac{f^{th}(s)}{f^{th}(s) + g^{th}(s)} = \frac{f^{th}(s)/g^{th}(s)}{1 + f^{th}(s)/g^{th}(s)}, $$

which depends only on the ratio $f^{th}(s)/g^{th}(s)$, with $f^{th}$ and $g^{th}$ derived from $k_+$ and $k_-$, respectively, through (11).

To limit the number of parameters that characterize this ratio, we assume, first that the attachment rate $k_+$ is non-zero only in an interval $S_f = [-\ell_+, \ell_+]$ around the nearest binding site; second, that $k_+$ is almost constant in this interval; and third that the detachment rate $k_-$ is strictly constant on this interval, see Figure 14. Using these hypotheses in the definitions of $f^{th}$ and $g^{th}$ (see Equation (11))
<figure><img src="image_3.png" /><figcaption>Figure 2: Illustration of the calibration ground principles for the PSE model. Note that the model elements represented here are chosen for illustration purposes and do not reflect the actual calibration. (a) Isometric conditions. The ratio of attached heads $\tilde{n}_{att}$ depends on the width and the height of the isometric population probability $\tilde{P}_1$ (orange area). The isometric force per attached head $\tilde{\tau}_c^{th}/\tilde{n}_{att}$ is a function of the isometric population probability $\tilde{P}_1$ and the equilibrium tension $\tau^{th}$. (b) Illustration of the situation at the end of phase II in a fast length step experiment. The isometric distribution is shifted due to the applied step $\delta s$. The legend is the same as in panel (a). (c) Computation of the stiffness in isometric conditions. The stiffness depends on the width of the isometric population probability $\tilde{P}_1$ and the asymptotic stiffnesses of the two attached potential wells $\kappa_{1pre}$ and $\kappa_{1post}$. (d) When the filaments slide past each other, the population probability $P_1$ is spread towards negative values of $s$ (solid red curve). The force is then given by the integral of $P_1(s)$ against the thermal equilibrium tension $\tau^{th}(s)$ (blue area).</figcaption></figure>
together with the (regularized) piecewise constant shapes proposed in Table 5, leads to the following approximations: for the attachment rate

$$f^{th}(s) \approx \begin{cases} \bar{f}^{th} = k_{\max} \alpha_{0pre} & \text{if } s \in [-\ell_{+}, \ell_{+}], \\ 0 & \text{otherwise,} \end{cases}$$

where $\alpha_{0pre}$ denotes the fraction of detached myosin heads in the pre-power stroke conformation, which is independent of $s$ but depends on the energy landscape $w_0$. For the detachment rate we have $g^{th}(s) \approx \bar{g}^{th} = k_{\min}$ for $s \in [-\ell_{+}, \ell_{+}]$. In summary, the isometric population probability $\tilde{P}_1$, depends only on the ratio $\bar{f}^{th}/\bar{g}^{th}$ and on the support width $d_f = 2\ell_+$ (see illustration in Figure 2(a)). If we carry on with these assumptions, we obtain the following approximations: (i) for the isometric ratio of attached heads

$$\tilde{n}_{att} = \frac{1}{d_a} \int_{s^{-}}^{s^{+}} \tilde{P}_1(s) ds \approx \frac{d_f}{d_a} \frac{\bar{f}^{th}/\bar{g}^{th}}{1 + \bar{f}^{th}/\bar{g}^{th}} \quad (\text{see the orange area in Figure 2(a)}) \text{ and (ii) for the isometric}$$

force per attached head $\frac{\tilde{\tau}_c^{th}}{\tilde{n}_{att}} = \frac{1}{\tilde{n}_{att}} \frac{1}{d_a} \int_{s^{-}}^{s^{+}} \tau^{th}(s) \tilde{P}_1(s) ds \approx \frac{1}{d_f} \int_{-\ell_+}^{\ell_+} \tau^{th}(s) ds$.

Given the bi-quadratic form postulated for the energy landscapes, the equilibrium tension $\tau^{th}$ depends on three parameters (see Figure 2(a)): the asymptotic stiffnesses $\kappa_{1pre}$ and $\kappa_{1post}$ and the position of the potential barrier $\ell_1$ of the attached energy landscape $w_1$, see Table 5. The stiffness $\kappa_{1post}$ will be univocally determined from the fast isotonic transients (see next paragraph), which leaves only 3 parameters ($\ell_+$, $\kappa_{1pre}$ and $\ell_1$) to determine the tension per attached heads. Once these parameters are fixed, the ratio $\bar{f}^{th}/\bar{g}^{th}$ is calibrated to match the measured fraction of attached heads, see Table 2.

<figure><img src="image_4.png" /><figcaption>Figure 3: Experimental data from fast transient experiments in force control conditions. (a) Imposed tension on the preparation. (b) Transient response following the tension drop. The response display three different phases indexed by 1, 2 and 3 (see text for more explanations). (c) Signature of phase I and phase II. The amount of shortening (per half-sarcomere) at the end of phase I and phase II, denoted respectively by $L_1$ and $L_2$ is reported as a function of the relative tension drop.</figcaption></figure>

**Fast isotonic transients** A typical shortening response of a cardiomyocyte to a sudden force drop applied within ~200 µs at the peak force $T_0$ of an isometric twitch contraction is illustrated in Figure 3(a) & (b).

The response in length displays three phases. First (phase I) an instantaneous shortening – $L_1$ – is observed alongside the load step itself. Second (phase II), after the force has stabilized to its prescribed value, a second shortening – up to $L_2$ – happens within the next ~2 ms. After this second phase, the system enters a steady state characterized by a constant load-dependent shortening velocity (phase III). To reproduce the signature of phase I and phase II, namely the curves $L_1(T_c)$ and $L_2(T_c)$ shown in Figure 3(c), we compute the response of our system to a rapid change in the filament relative position (per half-sarcomere) $\delta s$ and compute $T_1(\delta s)$ (phase I) and $T_2(\delta s)$ (phase II). This approach is justified by the experimental observation that the relations $T_1(\delta \ell_{hs})$ (phase I) and $T_2(\delta \ell_{hs})$ (phase II), where $\delta \ell_{hs}$ is the applied step in length (per half-sarcomere), are identical to $L_1(T_c)$ and $L_2(T_c)$, respectively [Piazzesi et al., 2002].

Due to the filament compliance, the displacement seen by the cross-bridges $\delta s$ is different from the displacement $\delta \ell_{hs}$ measured or imposed in experiments. To account for this difference, we map the
displacement $\delta s$ applied in our simulations for phase I and phase II to its $\delta \ell_{hs}$ counterpart, which corresponds to the measured $L_1$ and $L_2$, using $L_{1,2}(\delta s) = \delta s + C_m[T_{1,2}(\delta s) - T_0]$ [Caruel et al., 2019], where the filament compliance $C_m$ is given in Table 1.

We assume that, during the first phase (from $\delta \ell_{hs} = 0$ to $\delta \ell_{hs} = L_1$), all internal conformational degrees of freedom can be considered “frozen”, and the fraction of attached heads constant. Phase I is then a purely elastic response characterized by the relation (see Appendix B.2)

$$T_1(\delta s) \sim T_0 + \kappa \rho_{\text{surf}} \tilde{n}_{\text{att}} \delta s, \qquad (24)$$

which allows to give a first estimate of the parameter $\kappa$, found to be 1.07 pN nm$^{-1}$ in [Pinzauti et al., 2018]. When computing the response of the system to an idealized instantaneous length step, this choice of $\kappa$ allows to match the experimental data. However, in the length step experiments, the step is not instantaneous (its duration being about 100 µs in [Pinzauti et al., 2018]). In this case, viscous effects come into play and $\kappa$ must be calibrated alongside the viscosity $\eta$.

The rapid shortening in phase II is the specific signature of the relaxation of the internal degrees of freedom parametrizing the power-stroke conformational change, not involving attachments nor detachments Huxley & Simmons [1971]. Consequently, the $T_2(\delta s)$ relation predicted by the PSE model can be computed as (see Figure 2(b))

$$T_2(\delta s) \approx \frac{\rho_{\text{surf}}}{d_a} \int_{s-}^{s^+} \tilde{P}_1(s - \delta s) \tau^{th}(s) ds = \frac{\rho_{\text{surf}}}{d_a} \int_{s-}^{s^+} \tilde{P}_1(s) \tau^{th}(s + \delta s) ds. \qquad (25)$$

For large shortening, one can further assume that the whole population of attached cross-bridges are in the post-power stroke conformation. The fact that the observed relation $L_2(T_c)$ is linear for large shortening (see Figure 3(c)) supports the choice of the bi-quadratic form of $u_1$, see (1) and Table 5. We can then derive the following large shortening approximation (see Appendix B.3)

$$T_2(\delta s) = \frac{\rho_{\text{surf}}}{d_a} \int_{s-}^{s^+} \tilde{P}_1(s) \tau^{th}(s + \delta s) ds \approx \rho_{\text{surf}} \tilde{n}_{\text{att}} \frac{\kappa \kappa_{1\text{post}}}{\kappa + \kappa_{1\text{post}}} (\delta s + \tilde{s}_1 + a), \qquad (26)$$

which allows to univocally calibrate both $\kappa_{1\text{post}}$ and $\tilde{s}_1$, independently of the form of the population probability $\tilde{P}_1$, indeed. Equation (26) is an illustration of the explicit link between the stochastic model—characterized by the stiffnesses defining the attached energy potential $w_1$—and the PSE model—characterized by the tension $\tau^{th}$.

The slope of the $T_2(\delta s)$ curve at the origin, which is given by $\frac{\rho_{\text{surf}}}{d_a} \int_{s-}^{s^+} \tilde{P}_1(y) \partial_s \tau^{th}(s) ds = \rho_{\text{surf}} \tilde{\kappa}_c^{th}$, aggregates the contributions of the isometric population probability $\tilde{P}_1$ and the attached energy landscape through the equilibrium average tension $\tau^{th}$. Its value results from the balance between heads contributing with a negative stiffness and heads contributing with a positive stiffness, see Figure 2(c). Increasing the width $d_f$ of $\tilde{P}_1$ leads to an increase in the stiffness and, conversely, decreasing $d_f$ decreases the stiffness, possibly to negative values. Therefore, the measured slope of the $T_2(\delta \ell_{hs})$ curve at the origin imposes a limitation not only on the support of the attachment rate $d_f$ but also indirectly on the ratio $\bar{f}^{th}/\bar{g}^{th}$ through the constraint on the isometric ratio of attached heads $\tilde{n}_{\text{att}}$, see above.

Altogether, the data obtained from the isometric state and the fast transient experiments, allow to fully calibrate the isometric population probability $\tilde{P}_1$ (characterized by $\ell_+$ and the ratio $\bar{f}^{th}/\bar{g}^{th}$) and the attached energy landscape (parametrized by $\kappa$, $\kappa_{1\text{pre}}$, $\kappa_{1\text{post}}$, $\tilde{s}$, $\ell_1$) along with the viscosity $\eta$.

The ratio $\bar{f}^{th}/\bar{g}^{th}$ is fixed but a particularization of the attachment rates $\bar{f}^{th}$ and $\bar{g}^{th}$ is still to be done. The shape of the detachment rate $g^{th}$ for value of $s$ non close to the origin $s = 0$, along with the detached energy landscape $w_0(x, y)$ also remain to the calibrated.

The parameters of $w_0(x, y)$ are less constrained than the other parameters of the model. Indeed, only few data characterizing the behavior of the detached myosin heads are available, most of them being obtained with solutions of myosin heads and actin sites in vitro, whose applicability in the context of modeling the actin-myosin interaction in the sarcomere may not be straightforward. We choose to calibrate the detached double well potential so that the pre-power stroke well position $y_{0\text{pre}}$ and the energy barrier position $\ell_0$ correspond to that of the attached double well potential $y_{1\text{pre}}$ and $\ell_0$, respectively. The stiffness of the two wells $\kappa_{0\text{post}}$ and $\kappa_{0\text{pre}}$ and the position of the post-power stroke well $y_{0\text{post}}$ are set such that the jumps occurring between the attached and detached states are associated with a
"maximal" energy loss given the constraint that $w_1$ is now fixed and that the energy brought by ATP in the detachment process $\mu_T$ is equal to 100 zJ [Barclay, 2015], so that the reverse rates, which have been neglected, are actually small compared to the forward rates.

<figure><img src="image_5.png" /><figcaption>Figure 4: Summary of steady-state shortening experimental data measured on rat cardiac muscles tra-
beculae with data points obtained in various sarcomere stretch and thin filament activation conditions.
These data have been registered for decades with a great consistency in the results. (a) Force-velocity
curve. (b) Produced work. (Black) Data from [Van Heuningen et al., 1982] at 25°C. (Yellow) Data from
[Daniels et al., 1984] at 25°C. (Lilac) Data from [de Tombe & ter Keurs, 1990] at 25°C. (Orange) Data
from [de Tombe & ter Keurs, 1992] at 25°C. (Blue) Data from [Caremani et al., 2016] at 27°C.</figcaption></figure>

**Global curvature of the force-velocity relation** The force-velocity curve is the last element that will allow to close the calibration. We have gathered experimental data from several works published over the past four decades, which show a remarkable consistency, see Figure 4. The typical shape of the curve is parametrized by its curvature, the maximum shortening velocity and the near isometric behavior.

The global curvature of the force-velocity curve results from the balance between two effects: (i) a “stiffness effect” as the shortening of the sarcomere reduces the force generated by attached myosin heads, and (ii) a “cycling effect” as the myosin heads can detach from a position where they exert low or even negative force and reattach somewhere else on the thin filament where they generate a positive force. Both effects are illustrated in Figure 5.

The cycling effect can be illustrated with the PSE model by varying the attachment and detachment
rates $k_{max}$ and $k_{min}$ (thus augmenting the cycling rate) while maintaining the ratio $k_{max}/k_{min}$ constant
(see Figure 5(a)). These variations also modify the thermal equilibrium transition rates $f^{th}$ and $g^{th}$,
with a ratio $\bar{f}^{th}/\bar{g}^{th}$ that remains constant. Since the heads always bind in a region where they exert
a high force (see Figure 2(a)), increasing the cycling rate increases the steady-state force produced at a
given velocity.

In the PSE model, the equilibrium state in the energy potentials $w_\alpha$ is assumed to be always reached.
Therefore the apparent force-shortening relation is always given by the $T_2(\delta s)$ curve whose slope is an
apparent stiffness, see (25) and Figure 2(a). To illustrate the “stiffness effect”, we vary this apparent
stiffness by changing the value of $\kappa$, see Figure 5(c). Since the value $\delta s = -(a + \tilde{s}_1)$ corresponding
to $T_2(\delta s) = 0$ is independent from $\kappa$ (see Equation (26)), increasing $\kappa$ increases the tension for low
shortening ($\delta s > -(a + \tilde{s}_1)$), and decreases the tension for large shortening ($\delta s < -(a + \tilde{s}_1)$).

At low shortening velocity, the heads experience a moderate displacement before detachment and therefore exert a higher force at higher apparent stiffnesses. Conversely for high shortening velocity, the attached heads can reach a region where they start to be under compression (negative forces) before they detach, and a higher apparent stiffness then lowers the average force. These tendencies are illustrated in Figure 5(d).

At this stage of the calibration process, the stiffness parameters are already fixed, the global shape
<figure><img src="image_6.png" /><figcaption>Figure 5: Parametric study on the force-velocity curve predicted by the PSE model. (a) & (b) Illustration of the “cycling effect”. (a) The attachment rate $k_{max}$ is varied while the ratio $k_{max}/k_{min}$ is maintained constant. The condition $k_{max}/k_{min}$ constant also implies that $\bar{f}^{th}/\bar{g}^{th}$ is constant. (b) Variation of the position $\ell_{-,l}$, which varies in the same way as the position of $g^{th}$ pseudo-asymptote $\bar{\ell}_l$ (see Figure 2(d)). (c) & (d) Illustration of the “stiffness effect”. Variation of the stiffness $\kappa$ and subsequent effects on the $T_2$-curve (c) and the force-velocity curve (d). The slope of the $T_2$-curve constitutes the apparent stiffness for slow time scale responses (as it is the case in steady-state shortening).</figcaption></figure>
of the experimental force-velocity curve is used to particularize the parameters $\bar{f}^{th}$ and $\bar{g}^{th}$.

**Large shortening velocity behavior** In steady state isotonic shortening, while the myosin and actin filaments slide pass each other, the population probability $P_1(s)$ is spread towards negative values of $s$. A zero force is obtained when the contributions of heads exerting a negative force and that of heads exerting a positive force balance, see Figure 2(d).

In the PSE model, the unloaded shortening sliding velocity is determined mainly by the position $\bar{l}_l$ of the vertical pseudo-asymptote of the detachment rate $g^{th}$ for negative values of $s$, see Figure 2(d). The value of $\bar{l}_l$ itself depends on the position of the detachment rate $k_-$ pseudo-asymptote – denoted $\ell_{-,l}$ – and on the reciprocal spatial characteristic length $\lambda_3$ for the stochastic model (see Table 5). Hence, different couples $(\ell_{-,l}, \lambda_3)$ may lead to the same force-velocity relation, and we present here only the effect of $\ell_{-,l}$, see Figure 5(b). When the vertical pseudo-asymptote of the detachment rate is located at a higher position $\ell_{-,l}$, the detachment rate is increased for all value of $s$ at the vicinity of $s^-$. This results in an increase in the cycling rate, which, as presented in Figure 5(b), leads to an increase of the developed force. Hence the maximum shortening velocity is increased by increasing $\bar{l}_l$ through an increase of $\ell_{-,l}$.

Similarly, the position $\bar{l}_r$ of the vertical pseudo-asymptote of the detachment rate $g^{th}$ for positive values of $s$ – parametrized by $\ell_{-,r}$ and $\lambda_3$ – is chosen such that we can numerically ensure that $P_1(s = s^+) \approx 0$ when the sliding velocity is positive and corresponds to the maximal physiological value observed in the filling phase of the heart or in the relaxation phase of a fiber twitch contraction ($\dot{x}_c \sim 2 \mu m s^{-1}$).

**Near isometric behavior** The shape of the force-velocity curve in near isometric conditions – i.e. for shortening with characteristic time $\ell_{hs}/\dot{x}_c$ that is small with respect to the transition rates – is determined by a balance between two effects: a change in the number of attached heads and a change in the averaged force per attached head. This balance can lead to an increase or a decrease of the force at slow sliding velocities with respect to the isometric force and therefore to dynamic instabilities, generally referred to as “anomalous force-velocity relation”, which have been studied both experimentally Edman [1988]; Edman et al. [1997]; Edman & Curtin [2001] and theoretically Jülicher & Prost [1995, 1997]; Vilfan et al. [1999]; Duke [1999]; Guérin et al. [2011]; Månsson [2010, 2014, 2016].

In our reference model, the detachment function $g^{th} = k_{min}$ is constant over the whole interval $[s^+, s^+]$ except near the boundary where it diverges. In this case, we observe that the tension actually increases at slow shortening velocities compared to the isometric configuration, see Figure 5 and Figure 16(b, black line) in Appendix B.4.

We also show in Appendix B.4 that this effect can be modulated by introducing a plateau in the detachment rate $g^{th}$ with only a marginal influence on the rest of the force-velocity curve and, importantly, without having to modify the other model parameters.

In summary, the isometric indicators and the fast transients allow to unambiguously determine the isometric population probability $\dot{P}_1$ and the attached energy landscape $w_1$. Adding information from the force-velocity curve leads to the unique definition of the transition rates $f^{th}$ and $g^{th}$. Only the detached potential $w_0$ remains subjected to variability in the calibration. More experimental data targeting the detached state would be needed to reduce this variability.

### 3.1.2 PSE model vs Stochastic model

We recall that the PSE model is derived as the asymptotic limit of the stochastic model by performing the adiabatic elimination of the fastest time scales. We analyze here in detail the differences between the two population models induced by this assumption.

**Isometric indicators** We show in Table 2 that, in isometric condition, the PSE model and the stochastic model are both in good agreement with the experimental macroscopic indicators.

Once calibrated, our models can be used to predict the energetic properties of the contractile system in isometric conditions. We choose as an indicator of the energetic performance the “ATP tension cost”, which we define as the ATP consumption rate per myosin head in the sarcomere per unit of produced force. Combining the experimentally measured ATP consumption per unit volume (obtained with rat cardiac cells at 25 °C in skinned conditions ensuring that ATP is solely consumed by the cycling myosin heads, see [de Tombe & Stienen, 2007]) with the density of myosin heads per unit volume of
<table><caption>Table 2: Isometric physiological indicators for cardiac isometric experiments and models predictions. The experimental maximal stress and force per attached head are taken from [Caremani et al., 2016] and [Pinzauti et al., 2018], respectively. The maximal stress is obtained with intact cells with an extra-cellular calcium concentration higher than in physiological conditions, and we assume that these conditions correspond to a full activation of the thin filament. The experimental ratio of attached heads value is inferred from data obtained in [Pinzauti et al., 2018] in sub-maximal activation conditions knowing that the relationship between the developed force and the ratio of attached heads is linear. The ATP tension cost is then derived from data obtained in [de Tombe & Stienen, 2007] (see text for more explanation).</caption><thead><tr><th>Isometric indicators</th><th>symbol</th><th>Experiments</th><th>Stochastic model</th><th>PSE model</th></tr></thead><tbody><tr><td>Ratio of attached heads</td><td>n&#x0303;<sub>att</sub></td><td>0.15</td><td>0.150</td><td>0.154</td></tr><tr><td>Force per attached head</td><td>&tau;&#x0303;<sup>th</sup><sub>c</sub> / n&#x0303;<sub>att</sub></td><td>6.14 pN</td><td>6.15 pN</td><td>6.17 pN</td></tr><tr><td>Maximal total stress</td><td>T<sub>0</sub></td><td>118 kPa</td><td>115 kPa</td><td>119 kPa</td></tr><tr><td>ATP tension cost</td><td>A<sub>T</sub></td><td>0.109 /s/head/kPa</td><td>0.0875 /s/head/kPa</td><td>0.0878 /s/head/kPa</td></tr><tr><td>ATP consumption<br>at T<sub>0</sub> (/s/head)</td><td>J<sub>ATP</sub>/&mu;<sub>T</sub></td><td>12.8</td><td>10.1</td><td>10.4</td></tr></tbody></table>

1.25 × 10<sup>23</sup> heads/m<sup>3</sup> [Pinzauti et al., 2018] decreased by 21% due to the skinning process [Konhilas et al., 2002], we obtain an estimate of 0.109 /s/head/kPa for the ATP tension cost.

In the Lymn-Taylor cycle, ATP binding promotes the detachment of myosin from actin. The ATP consumption rate per myosin head is therefore equal to the flux of detachment per myosin head. For the PSE model, the tension cost is therefore given by

$$A_T = \frac{\text{detachment flux per myosin head}}{\text{produced force}} = \frac{\frac{1}{d_a} \int_{s-}^{s^+} g^{th}(s) \tilde{P}_1(s) \, ds}{\frac{\rho_{\text{surf}}}{d_a} \int_{s-}^{s^+} \tau^{th}(s) \tilde{P}_1(s) \, ds},$$

and we obtain a value of 0.0878 /s/head/kPa in relatively good agreement with the experimental data and with the prediction of the stochastic model, see Table 2. The fact that the models reproduce this indicator without further adjustments confirms that our parametrization is well constrained.

**Fast isotonic transients** We continue by comparing the calibrated models predictions of the fast transient response with experimental data. The calibration results are presented in Figure 6.

While the relation $T_2(\delta s)$ is explicitly defined in the PSE model (see Equation (25)), it has to be inferred from the average trajectory of the internal variables $X^t$, $Y^t$ and $\alpha^t$ for the stochastic model. We used the method proposed in [Ford et al., 1981; Caruel et al., 2019] to reconstruct the $T_2(\delta s)$ relation.

The fundamental assumption underlying the derivation of the PSE model is precisely to neglect the time required for the equilibration of the internal mechanical degrees of freedom corresponding to the power stroke. There is thus a theoretical equivalence between the end of phase II for the stochastic model and the PSE model. This translates into the numerical simulations (see Figure 6) showing that the $(\delta s, T_2)$-curves obtained with the stochastic and the PSE models are effectively the same. As a result, the elastic response of the PSE model aggregates the dynamics of these two – fast – power stroke elements into a single – instantaneous – nonlinear elastic dynamics.

The accordance between the predictions of the calibrated PSE model (blue line in Figure 6(a)) and the predictions of the stochastic model (green triangles in Figure 6(a)) validates the hypotheses of the PSE model, that the force length relation characterizing the end of phase II corresponds to the thermal relaxation of the internal degrees of freedom.

The response in phase I, which corresponds to the fastest time scale, is already averaged into the equilibrium nonlinear elastic element of the PSE model (blue line in Figure 6(a)), and therefore it can only be captured by the stochastic model.

In our simulations the $T_1(\delta s)$ relation is obtained simply by reporting the value of the tension at the end of the applied length step, as in experiments. Since the length step is not applied instantaneously, the Phase I response is viscoelastic, and therefore the $T_1(\delta s)$ relation depends on both the stiffness $\kappa$ and the
<figure><img src="image_7.png" /><figcaption>Figure 6: Comparison of the calibrated model prediction and experimental data for the fast transient response. The simulation data where obtained by applying a shift in the binding site position *s* within 100 µs (marks in green) starting from the isometric steady state. (a) Tension-length relationship characterizing the end of phase I and II. The results are displayed as function of the half sarcomere length change, taking into account the filaments deformation corresponding to the applied step. (b) Influence of viscosity on the force *T*₁- and *T*₂-curves predicted by the stochastic model.</figcaption></figure>
drag coefficient η, which are then calibrated together. The result of this calibration is shown in Figure 6(a) (green symbols) where we have set κ = 1.34 pN nm⁻¹ (slightly larger than the value κ = 1.07 pN nm⁻¹ obtained by assuming a purely elastic response, see (24) and Figure 3) and η<sub>ref</sub> = 0.0972 ms pN nm⁻¹, see Table 5. The effect of changing the value of the drag coefficient is illustrated in Figure 6 (b). Remarkably none of the values tested affects the T₂(δs) response. Furthermore, one can recover an almost purely elastic response by doubling the drag coefficient (red symbols). The choice of the value of η has also a strong influence on the rate of the tension recovery in phase II, we checked that our choice allows to keep this rate in the physiological range, see section 4.1.1.

<figure><img src="image_8.png" /><figcaption>Figure 7: Steady state shortening response. (a) Force-velocity curve for the stochastic model and the PSE model. (b) Produced work rate and energy input rate per myosin head in steady state shortening conditions. (c) Efficiency of the transduction of ATP energy into work for the PSE model in steady state shortening conditions. (d) Influence of the viscosity on the force-velocity curve computed with the stochastic model.</figcaption></figure>

**Isotonic shortening velocity** The force-velocity relations obtained with both the stochastic model and the PSE model are shown in Figure 7(a). The predictions of the two models are relatively similar between $\dot{x}_c = 0$ and $-\dot{x}_c = 3 \mu\text{m s}^{-1}$. Beyond this value, the tension computed from the stochastic model is systematically lower than the one computed from the PSE model.

In general, the discrepancy between the two models lies in the fact that the fundamental assumption
that the distribution of the myosin heads in the attached and detached energy potentials are distributed
according to the thermal equilibrium distribution is not fully satisfied, implying that the equilibration
of the myosin heads in the attached and detached potential is not always very fast compared to other
dynamical processes. As an illustration of this statement, one can note that the discrepancy between the
force-velocity relations predicted by the two models can be reduced by a decrease in the drag coefficient
η (that is by speeding up the equilibration process) or conversely increased by increasing this coefficient.

To understand the origin of the difference between the response of the two models, we compare in
Figure 8 the ratio of attached heads in steady-state conditions (panels (a), (b) & (c)), which we denote by $P_1^\infty$, along with the associated tension (panels (d), (e) & (f)) in isometric conditions and at two different sliding velocities. In isometric conditions (see Figure 8(a)), the ratio $P_1^\infty$ for the two models coincide for $s > 3$ nm while below this value, the PSE model predicts a higher fraction of attached bridges. This observation was already made in [Caruel et al., 2019] and explained again by the fact that the characteristic relaxation time in the detached energy potential is not much faster than the inverse of the attachment rate (see Figure 21 in appendix). Cycling heads thus spend more time in the detached state resulting in a lower ratio of attached heads. For larger values of $s$ this limitation is not involved since the post-power-stroke configuration is not populated, which prevents the heads from completing their cycle. Altogether the ratio of attached heads is higher for the PSE model than for the stochastic mode (see also Table 2) with heads populating configurations that produce on average a positive force (see Figure 15(b) in appendix). The global developed tension is thus higher for the PSE model as can be seen in Figure 8(d).

For nonzero sliding velocities (see Figure 8(b and c)), the shapes of the ratios $P_1^\infty$ predicted by the two models are similar (with the stochastic model predicting slightly lower values than that from the PSE model) except for low values of $s$ where the fraction of heads vanishes in the PSE model while ~10% of cross-bridges remain attached within the stochastic model. This result is explained by the fact that, if the equilibration in the attached potential is very fast, which corresponds to the conditions of the PSE model, heads are present in configuration where the detachment rate $k_-$ that diverges, resulting in an averaged detachment rate $g^{th}$ for the PSE model that also diverges when $s \approx -20$ nm. In the stochastic model, attached heads can still encounter a detachment rate as low as $k_{min} + k_* = 228$ s$^{-1}$ for large values of $Y$, where the tension is negative (see Figure 13 in the appendix). Therefore, in the stochastic model for $s < -20$ nm some heads remain attached and produce on average a negative force, while they are already detached in the PSE model, leading to a lower global force, see Figure 8 (e) & (f) for values of $s$ around -20 nm.

An additional effect also contributes to the higher developed force in the PSE model. In Appendix C (Figure 17), we show that the distributions of the conformational variable $Y$ predicted by the two models start to differ over a limited range of $s$ as velocity increases, which shows that the sliding kinetics start to compete with the equilibration process kinetics. This behavior is naturally expected since increasing the filament sliding velocity impairs the validity of the PSE model fundamental assumption that the equilibration processes in the energy potentials is much faster than any other process in the system. At a given value of $s$, this competition leads for the attached heads to an increased likelihood to remain in the pre-power stroke configuration, in which less force is developed than in the post-power stroke configuration. This behavior is illustrated in Figure 8(e) & (f) for values of $s$ in the range [-4 nm, 0 nm] and see also Figure 17 for the same values of $s$ at non-zero sliding velocities.

However, we show that beyond these discrepancies, the steady state normalized probability distributions $p^\infty(x, y, 0; s) / \iint p^\infty(x, y, 0; s)dydx$ and $\bar{p}^\infty(y; s) / \int \bar{p}^\infty(y; s)dy$ obtained from the stochastic model fit very well (with the above mentioned exceptions) with the thermal equilibrium distributions $p_0^{th}(x, y)$ and $p_1^{th}(y; s)$ used in the PSE model, respectively. Note that in the region where the distributions display a large discrepancy between the two models, the real quantity of interest, which is the force, displays a less pronounced difference (see Figure 8 (d), (e) & (f)). Overall this match justifies the validity of the assumption made to derive the PSE model from the stochastic model. Besides, the force-velocity differ between by less than 30% in the physiological regions in terms of mechanical loading $T_c \sim 0.5T_0$ (see Section 3.2.2 for more explanations on the definition of the physiological region).

Finally, as in isometric contraction, we are able to make some prediction regarding the contraction energetics, see Figure 7(b). The work production rate and the associated ATP energy consumption rate per myosin head can straightforwardly be computed from the steady state isotonic shortening response. They are defined by

$$ \begin{cases} \dot{W} = \dot{x}_c \tau_c, \\ \dot{\varepsilon} = \frac{\mu_T}{d_a} \int_{s^-}^{s^+} g^{th}(s) \tilde{P}_1(s) ds, \end{cases} $$

respectively, where $\mu_T$ is the chemical potential brought by ATP. From the latter the yield $\eta_T = \dot{W}/\dot{\varepsilon}$ can be calculated (see Figure 7(c)).

Again without any particular adjustment of the parameter and using the validated ATP turnover
<figure><img src="image_9.png" /><figcaption>Figure 8: (Upper row) Fraction of attached cross-bridges predicted by the PSE model (in blue) and by the stochastic model (in black) at different shortening velocities, $-\dot{x}_c = 0$ (a), $3 \mu m s^{-1}$ (b) and $6 \mu m s^{-1}$ (c). (Lower row) Force developed by cross-bridges predicted by the PSE model (in blue) and by the stochastic model (in black) at different shortening velocities, $-\dot{x}_c = 0$ (d), $3 \mu m s^{-1}$ (e) and $6 \mu m s^{-1}$ (f).</figcaption></figure>
energy $\mu_T = 100$ zJ, our model recovers the experimental ATP consumption rate in isometric conditions obtained from the tension cost and the maximal tension (see Table 2). Our prediction of the consumption rate and the yield during shortening remains to be validated experimentally.

## 3.2 Macro models

### 3.2.1 Calibration

We recall that the main assumption made to obtain the simplified macro-models was to consider that the sum $f(s)+g(s)$ was constant over the whole $[s^{-}, s^{+}]$ interval, see Section 2.2. To compare the PSE model with the macro-models we thus need to define “equivalent” transition rates for the different versions of the macroscopic models.

<figure><img src="image_10.png" /><figcaption>Figure 9: (a) & (b) Ingredients of the macroscopic models and the PSE model: transition rates (a) and free energy of the cross-bridge (b). Note that the energy levels are defined up to a constant. (c) & (d) Models predictions and corresponding experimental data: elastic response (c) and force-velocity relation (d). Experimental data are as in Figure 4.</figcaption></figure>

We chose the attachment rates of the macro-models such that the maximal value $\bar{f}$ of the attachment rate is equal to $\bar{f}^{th}$ for the PSE model and we adjusted the width of region where the attachment function $f(s)$ does not vanish, such that the isometric ratio of attached heads $\tilde{n}_{att}$ matches the data while verifying the condition that $f + g$ is constant (see Figure 9(a)). The minimal value of the detachment rate $g$ is denoted by $\bar{g}$.

The first macro-model (linear macro-model) assumes a linear elastic cross-bridges characterized by two parameters: a stiffness $\kappa_{xb}$ and a reference length $s_0$. We consider two cases for the definition of these constants. In the first case ($T_1$-linear macro-model), we chose the value of the stiffness measured experimentally $\kappa_{xb} = \kappa = 1.34$ pN nm$^{-1}$ and the pre-strain $\tilde{s} = 4.59$ nm that leads to an isometric
<table><caption>Table 3: Target isometric physiological indicators for the calibration on cardiac data of the PSE model and the macro-models.</caption><thead><tr><th>Isometric indicators</th><th>Symbol</th><th>PSE model</th><th>T<sub>1</sub>-linear<br>macro-model</th><th>T<sub>2</sub>-linear<br>macro-model</th><th>nonlinear<br>macro-model</th></tr></thead><tbody><tr><td>Ratio of<br>attached head</td><td>n&#770;<sub>att</sub></td><td>0.154</td><td>0.153</td><td>0.153</td><td>0.153</td></tr><tr><td>Force per<br>attached head</td><td>&tau;&#770;<sup>th</sup><sub>c</sub> / n&#770;<sub>att</sub></td><td>6.17 pN</td><td>6.14 pN</td><td>6.32 pN</td><td>6.14 pN</td></tr><tr><td>Maximal total<br>stress</td><td>T<sub>0</sub></td><td>119 kPa</td><td>117 kPa</td><td>121 kPa</td><td>117 kPa</td></tr><tr><td>ATP tension cost<br>(/s/head/kPa)</td><td>A<sub>T</sub></td><td>0.0878</td><td>0.107</td><td>0.104</td><td>0.107</td></tr><tr><td>ATP consumption<br>at T<sub>0</sub> (/s/head)</td><td>J<sub>ATP</sub>/&mu;<sub>T</sub></td><td>10.4</td><td>12.6</td><td>12.6</td><td>12.6</td></tr></tbody></table>

force per attached head of 6.14 pN (see free energy profile in Figure 9(b)). In the second case ($T_2$-linear model), we aim at approximating the response at the end of phase II ($\delta s, T_2$) with a linear elastic response. The cross-bridge stiffness is then taken equal to $\kappa_{xb} = 0.52 \text{ pN nm}^{-1}$ and the pre-strain $\tilde{s}$ is then set to 12.15 nm to obtain an isometric force per attached head of 6.31 pN (see free energy profile in Figure 9(b)).

The second macroscopic model (nonlinear macro-model) assumes a nonlinear elastic energy for the cross-bridge. We choose to define the free energy of the attached state with a polynomial of order 6 ($n = 6$) and calibrate the model so that the elastic response follows the ($\delta s, T_2$) response (the model parameters are given in Table 6), see Figure 9((b) and (c)).

### 3.2.2 PSE model vs macro-models

**Isometric indicators** The isometric indicators are presented in Table 3. With the proposed calibrations, all models are able to reproduce the key isometric indicators. The good agreement between the four models is due to the fact that the population probabilities $P_1$ are similar (see rates in Figure 9 (a)) and that the forces developed in the interval $S_f$ are similar (similar slope of the four free energy potentials in Figure 9 (b)).

**Effective elastic response** By construction, the PSE model and the macro-models cannot reproduce the fast transients. The parameters of the stochastic model has been chosen such that the effective elastic response of the PSE model reproduces the ($\delta l_{hs}, T_2$) relation obtained from the isotonic fast transients.

With the calibration choices made in the previous section, the elastic response of the $T_1$-linear macro-model corresponds to the ($\delta s, T_1$)-curve, while the elastic response of the $T_2$-linear macro-model corresponds to a linear approximation of the ($\delta s, T_2$)-curve (see Figure 9(c)). As expected the nonlinear macro-model best reproduces the results of the PSE model as regards to the effective elasticity.

**Steady-state isotonic shortening** We compare the steady-state isotonic behavior of the macro-models in Figure 9(d). The force-velocity relation corresponding to both linear models are straight lines given by (21), see Figure 9(d, dotted purple and green dashed lines). Note that a curved force-velocity relation can be obtained with a macroscopic model using a linear elastic cross-bridge with the addition of a velocity dependent term in the definition of the detachment rate [Chapelle et al., 2012; Månsson, 2010]. However, in the hierarchical framework, these terms do not appear and are thus not considered here.

The force of the $T_2$-linear macro-model (green dashed line) is larger at any fixed shortening velocity than that of the $T_1$-linear macro-model, which has a larger stiffness, in accordance with the prediction of (21). This illustrates the “stiffness effect” discussed in Section 3.1.1.
For the nonlinear macro-model, the force-velocity curve is given by (23). It remains close to the PSE model at low shortening velocity as expected from the higher force at small displacements. However the polynomial approximation leads to a dramatic force drop for $s < 20$ nm therefore, as the velocity increases, the nonlinear macro-model rapidly loses force and thus has a low maximum shortening velocity, again because of the “stiffness effect”.

**Macro-models as effective heart contraction models** To put the development of the macro-models into perspective, we compare the implication of the behavior differences between the two proposed versions (linear elastic and nonlinear elastic) in the context of modeling the ejection phase of the heart.

Considering experimental pressure-volume loops and end-systolic pressure-volume relations obtained on rats, the range of physiological loads can be estimated to be between 45% and 65% of the maximal load [Sato et al., 1998; Pacher et al., 2004]. In this region the shortening velocity measured experimentally has a value of about $2 \mu m s^{-1}$ (the physiological region is indicated in blue in Figure 9(d)). The same value is predicted by the PSE model and the $T_2$-linear macro-model, whereas the $T_1$-linear macro-model predicts a shortening velocity of $0.6 \mu m s^{-1}$.

Denoting the duration of ejection by $\tau_e$ and assuming that the heart cavity is spherical and that the thickness is negligible, we obtain a end-systolic volume $V_{ES}$ of

$$V_{ES} = V_{ED} \left( 1 + \dot{x}_c \frac{\tau_e}{\ell_{hs}} \right)^3$$

The ejection fraction being defined by $(V_{ED} - V_{ES})/V_{ED}$ and taking a duration of ejection equal to 45 ms and end-diastolic volume of 200 mL, the $T_1$-linear macro-model would predict an ejection fraction that is 66% lower than that predicted by the PSE model and the $T_2$-linear macro-model. Note that the ejection fractions estimated with the PSE model and the $T_2$-linear macro-model are outside of the physiological range observed (50-65%) at 37 °C [Pacher et al., 2004]. This is due to the fact that the models are calibrated on data obtained at 25 °C. Considering the data measured by de Tombe & ter Keurs [1990] at 30 °C, we see that the shortening speed in the physiological range is approximately $5 \mu m s^{-1}$, which leads to an ejection fraction of 53%.

The $T_1$-linear macro-model is thus not applicable in the context of heart modeling. On the contrary, the $T_2$-linear macro-model has good properties to be used in heart simulations. Its force-velocity curve presents a good match with the data over the whole physiological range and the elastic response is in good agreement with the experimental $(\delta s, T_2)$-curve.

The PSE model has the capability of reproducing the end of phase II response and the whole force-velocity curve. The computational complexity is further increased because the dynamics is governed by the partial differential equation (9).

The use of one or the other of the models will naturally depend on the goal of the simulation. For a heart simulation in which a default physiological active behavior is sufficient, the $T_2$-linear macro-model is probably a good trade-off between the physiological consistency and the computational cost. However, if the heart simulation is targeting the active behavior of the cardiac tissue, the PSE model should be recommended. Indeed, with its tighter link to the physiology, the PSE model main advantage is to be able to give a better physiological interpretation of the simulated behavior. Moreover, it can consistently transfer properties across time and space scales, for instance the impact of a cardiomyopathy that affect the actin-myosin interaction can then be incorporated in the model by an adjustment of the calibration and then brought to the macroscopic space scale.

# 4 Discussion

In this section we put our results in perspective with other theoretical and experimental works.

## 4.1 Limitations of the models

Despite its ability to reproduce the most fundamental physiological indicators of muscle contraction with a rather limited number of parameters, we list here the limitations our model hierarchy.
### 4.1.1 Viscoelasticity

Internal viscous damping controlling the kinetics of the relaxation of the conformational variables $Y$ is the essential link between the stochastic model and the PSE model. The lower the internal viscosity parameter $\eta$, the more valid the PSE approximation.

For the stochastic model, the parameter $\eta$ was calibrated together with the head stiffness $\kappa$ to match the observed $T_1$. The viscosity also controls the kinetics of the Phase II, but to our best knowledge data are only available for load step experiment, and it has been shown that the results may differ in length controlled experiments. Furthermore, in the length controled experiment the kinetics of the force recovery strongly depends on the filament elastic properties which are not taken into account in this paper Caruel et al. [2013]; Piazzesi et al. [2014].

Therefore we cannot conclude at this stage on the ability of the stochastic model to reproduce the rate of phase II. Nevertheless, in Figure 10 we show predictions of the model obtained with three different values of the viscosity. While data are lacking for a direct assessment of our viscosity, we note that the predicted range of values remains physiological Piazzesi & Lombardi [1995]; Piazzesi et al. [2002].

<figure><img src="image_11.png" /><figcaption>Figure 10: Influence of the internal viscosity parameter $\eta$ on the recovery rate associated with phase II.</figcaption></figure>

### 4.1.2 Rate of tension redevelopment

In the calibration process, we adjust the transition rates $k_+$ and $k_-$ to fit the experimental force-velocity curve. In particular, the parameters $k_{max}$ and $k_{min}$, or equivalently the value of the equilibrium transition rate $\bar{f}^{th}$ and $\bar{g}^{th}$, influence the global shape of the force-velocity curve (see Section 3.1.1).

de Tombe & Stienen [2007] propose an alternative way to calibrate these model parameters. They measure the tension redevelopment characteristic time in isometric conditions after applying a fast release-restretch maneuver to force the detachment of some myosin heads. With this experimental setup, the dynamics of the tension is supposedly not affected by the thin filament activation process. For the family of models derived from the Huxley'57 model, the tension redevelopment characteristic time $\tau_{tr}$ is given by $\tau_{tr} = \frac{1}{\bar{f} + \bar{g}}$. Coupling these data with the ATP consumption measurement, which is linked to the detachment rate $g$, they obtain an estimation of the transition rates $f$ and $g$. For rat cardiac muscle at 25°C, they obtain $\bar{f} = 35.0 \text{ s}^{-1}$ and $\bar{g} = 14.4 \text{ s}^{-1}$. These values are different from the ones obtained after our calibration. We have $\bar{f}^{th} = 205 \text{ s}^{-1}$ and $\bar{g}^{th} = 66 \text{ s}^{-1}$, which means that our model would fail to reproduce the kinetics of force redevelopment reported in de Tombe & Stienen [2007]. Nevertheless, if we use the value proposed by [de Tombe & Stienen, 2007] – all other things remaining equal –, we obtain the force-velocity curve presented in Figure 11, which does not match the experimental data. This comparison illustrate one of current challenges in the understanding and modeling of muscle contraction: how a given model can match both the observed power-output of the fiber – which necessitates a high cycling rate – and the rather low rate of force redevelopment – which necessitates low cycling rate. For a discussion of this conundrum, we refer to Månsson et al. [2015].

To resolve this issue Pertici et al. [2018] proposed a model (simplified version of Caremani et al. [2015]) that is able to capture the physiological skeletal muscle force-velocity curve without involving
<figure><img src="image_12.png" /><figcaption>Figure 11: Force-velocity curve using the transition rate values given in [de Tombe & Stienen, 2007], all other model parameters remain as given in Table 5.</figcaption></figure>

a cycling rate comparatively as high as in our calibration. To do so, they introduce the possibility for the attached myosin heads to slide to a neighbor actin site without detaching, thus restoring a part of their force generation capability without consuming ATP. Another attempt, consisting in prescribing a phenomenological dependence of the attachment rate on the shortening velocity, is considered in Månsson [2010]. While both attempts successfully reconcile the two apparently contradicting observations within a single model, to our best knowledge, further studies are required to assess the validity of their assumptions.

### 4.1.3 Compatibility with X-ray diffraction

As we discussed in Section 3.1.1, the width $d_f$ of the equilibrium attachment rate support $S_f$ (see Figure 2) is adjusted to match the indicators of the isometric contraction. We have $d_f \approx 8$ nm, which is larger than values reported for instance in [Piazzesi & Lombardi, 1995; Caremani et al., 2015; Månsson et al., 2015] but comparable to the values reported in [Smith & Mijailovich, 2008]. According to [Reconditi, 2006], the dispersion of the projected mass of the attached myosin heads on the actin filament in isometric contraction, which corresponds to the dispersion of the attachment rate, is about 5 nm and is thus incompatible with a wider value of $d_f$. To reconcile our model with these observations, one can consider a chain of half-sarcomeres as the elementary unit of the contraction instead of a single representative motor as it is done most of the time. While the tension-elongation relation $\tau_c^{th}(s)$ of a single motor shows a region of negative slope, Caruel & Truskinovsky [2018] have shown that, with the same parameters, the tension-elongation relation of a series arrangement of half-sarcomere is necessarily larger or equal to 0. If we take such a relation as a “constitutive behavior” of the PSE model, we would need a lower value of $d_f$ to match the isometric contraction indicators, indeed.

### 4.1.4 Comparison with other models

Our model can be compared to two other Huxley’57-type models, the first one developed in [Pertici et al., 2018] and the other in [de Tombe & Stienen, 2007], see Table 4.

One of the major differences is that our model considers the actin periodicity to be $d_a = 40$ nm, which corresponds to the periodicity of the double helix, whereas the other models consider the distance between monomers: $d_a = 5.5$ nm. To match the experimental indicators of the isometric contraction both Pertici et al. [2018] and de Tombe & Stienen [2007] define a positive attachment rate over the interval $d_f = d_a$ whereas our model has $d_f < d_a$. The consequence of this choice is that the fraction of the heads that are attached in isometric contraction $\tilde{n}_{att}$ can, in the PSE model, be decorrelated from the duty ratio $r$, i.e. the fraction of the cycle that a motor spends bound to actin. Indeed, the fraction
<table><caption>Table 4: Coarse model comparison in isometric conditions. Note that we consider in this table the model proposed by Pertici et al. [2018] without the possibility of actin shifting so that we can compare the calibration of the transition rates.</caption><thead><tr><th></th><th>Pertici et al. [2018]</th><th>PSE model</th><th>de Tombe & Stienen [2007]</th></tr></thead><tbody><tr><td>d<sub>a</sub></td><td>5.5 nm</td><td>40 nm</td><td>5.5 nm</td></tr><tr><td>d<sub>f</sub></td><td>5.5 nm</td><td>8 nm</td><td>5.5 nm</td></tr><tr><td>f&#773;</td><td>15 s<sup>-1</sup></td><td>205 s<sup>-1</sup></td><td>35 s<sup>-1</sup></td></tr><tr><td>g&#773;</td><td>40 s<sup>-1</sup></td><td>66 s<sup>-1</sup></td><td>14 s<sup>-1</sup></td></tr><tr><td>n&#770;<sub>att</sub></td><td>0.25</td><td>0.15</td><td>0.7</td></tr><tr><td>r</td><td>0.25</td><td>0.75</td><td>0.7</td></tr><tr><td>J<sub>ATP</sub>/μ<sub>T</sub></td><td>11 s<sup>-1</sup></td><td>10 s<sup>-1</sup></td><td>10 s<sup>-1</sup></td></tr><tr><td>τ<sub>tr</sub></td><td>2 × 10<sup>-2</sup> s ✔</td><td>3.5 × 10<sup>-3</sup> s ✗</td><td>2 × 10<sup>-2</sup> s ✔</td></tr><tr><td>F-V curve</td><td>✗</td><td>✔</td><td>✗</td></tr></tbody></table>

of attached heads and the duty ratio are given by

$$ 
\begin{aligned}
\tilde{n}_{\text{att}} &= \frac{1}{d_a} \int_{s^{-}}^{s^{+}} \tilde{P}_1(s) \, ds \approx \frac{d_f}{d_a} \frac{\bar{f}}{\bar{f} + \bar{g}}, \\
r &= \frac{\frac{1}{d_a} \int_{s^{-}}^{s^{+}} \frac{1}{g(s)} \tilde{P}_1(s) \, ds}{\frac{1}{d_a} \int_{s^{-}}^{s^{+}} \left( \frac{1}{f(s)} + \frac{1}{g(s)} \right) \tilde{P}_1(s) \, ds} \approx \frac{\bar{f}}{\bar{f} + \bar{g}},
\end{aligned}
$$

so that having $d_f = d_a$ necessarily implies $r = \tilde{n}_{att}$. Another consequence of having $d_f \neq d_a$ is that, a significant fraction of heads cannot attach and therefore do not consume energy. The consumption of ATP is defined by

$$ \frac{J_{\text{ATP}}}{\mu_T} = \frac{1}{d_a} \int_{s^{-}}^{s^{+}} g(s) \tilde{P}_1(s) \, ds \approx \bar{g} \tilde{n}_{\text{att}} = \bar{g} \frac{d_f}{d_a} \frac{\bar{f}}{\bar{f} + \bar{g}} $$

which shows that, even with high attachment and detachment rates, the energy consumption of the PSE model can be kept low, again by tuning the ratio $d_f/d_a$.

As explained above the two-state models presented in Table 4 cannot reproduce both the rate of force recovery $\tau_{tr}$ and the force velocity relation at the same time. In the case of Pertici et al. [2018], the two state model as to be supplemented by the addition of the possibility for actin shifting without detachment in order to reproduce the power output. In this paper, we chose to concentrate on reproducing the force-velocity relation assuming that other physiological mechanisms may play a role in the specific kinetics of the force redevelopment.

## 4.2 Limitations of our calibration

The experiments are performed at 25–27 °C and not at body temperature for technical reasons: at higher temperatures the first phase of the response become too fast to allow a precise measurement of $T_1$. Moreover, temperature is known to affect the shape of the force-velocity curve, the value of the unloaded shortening velocity $V_0$ and the value of the force developed in isometric conditions [de Tombe & ter Keurs, 1990], but also the apparent cycling rates of the cross-bridges [de Tombe & Stienen, 2007]. However, the only range of temperature, at which all types of quantitative characterization of the cardiac muscle behavior are available, is around 25 °C. We have therefore performed the calibration of our models with data measured at this temperature.
Furthermore, the actin-myosin interaction quantitative properties vary between species. Experiments performed on rat trabeculae may therefore not give precise information and the behavior of human cardiac muscles.

To give an idea of the potential variability induced by this parameters, we can compare the data in isotonic shortening protocols measured on rat cardiomyocytes at 25 °C [Caremani et al., 2016] and on skeletal muscles at 4 °C performed on frog [Piazzesi et al., 2002]. These two experimental conditions variations combined results in a rate of phase II that is about four times larger on cardiac sample than on skeletal samples, and the transition rates needed with our model to match the measured power-output of a fiber are about six times larger on cardiac sample than on skeletal samples, see Caruel et al. [2019].

# 5 Conclusion

In this paper we have developed a hierarchy of interconnected models of the actomyosin system, from the more refined, accounting for a large number of physiological indicators, to the more coarse-grained, allowing for fast simulation at the organ scale.

We have shown that our starting stochastic mechanical model can be calibrated to reproduce essential mechanical indicators of the cardiac contractile unit that can be obtained experimentally. The calibration procedure leads to a robust definition of the parameters value, which allow, in particular, for specifically identifying the effect of a parameter value on the output of the model.

Our second model (PSE model) is based on the assumption that the internal stochastic variables parametrizing the power stroke conformational change quickly relax towards their equilibrium distribution. It can be viewed as an instance landmark Huxley'57 family of models, whose dynamics takes the form of a local PDE in a finite element simulation environment. Unlike existing model of this type, its parameters are directly inherited from the stochastic model and not postulated a priori. Importantly this allows to rigorously calibrate the PSE model with experimental data obtained at fine time scales. Moreover, in the region of in vivo physiological loadings in vivo ranges, the force-velocity curve of the stochastic model and the PSE model are similar, see Fig. 7.

By imposing specific constrains on the parameters of the PSE model, we derived the system of ordinary differential equations, accounting for the moments dynamics characterizing our third model. This type of moment-based approach can be combined with standard geometrical simplification of the heart to construct reduced dimensional models of the organ that can be simulated in real time and therefore used as a primary estimation tool or for rapid pre-calibration purposes Caruel et al. [2014] or in clinical contexts [Le Gall et al., 2020].

The advantage of our approach is that the hierarchical relation between the models allow to predict the consequences of microscopic actin-myosin constitutive behavior variations, for instance as a result of a disease or a treatment affecting the molecular motor [Robert-Paganin et al., 2019; Woody et al., 2018], at the larger time and spaces scales.

In addition, the coarse graining approach is well controlled i.e. the loss of information inherent to the simplifications is quantified and the associated hypothesis can always be tested. Therefore, a more refined model can be called upon if needed for a specific application, without having to start over a tedious calibration procedure. Moreover, the method is not specific to our starting stochastic model: it can be applied to the widely used chemical-mechanical modeling framework [Eisenberg et al., 1980; Caremani et al., 2015], which, through elimination of the fastest chemical-like reactions, can be also reduced to a simple population model and further to a moment-based model [Zahalak, 1981].

Finally, we recall that only the basal mechanical behavior of the actomyosin system have been addressed in this work. In particular, all the regulation aspects, essential for a relevant simulation of the organ behavior, remain out of the scope of our study. Detailed activation models also involve the simulation of the dynamics of populations of active and inactive agents – actin sites or myosin heads – using similar chemical analogy. Therefore, building a similar model hierarchy, might also be a solution for an organ-scale simulation of the microscopic activation and regulation physiology.

# 6 Acknowledgements

The authors are grateful to Marco Linari and his collaborators (University of Florence), whose insightful comments helped to considerably improve this manuscript. This work also greatly benefited from regular
interactions with Dominique Chapelle and Lev Truskinovsky. F.K. thanks the AMX Ph.D. fellowship from École polytechnique for financial support.

# 7 Compliance with ethical standards

**Conflict of interest** The authors declare that they have no conflict of interest.

# Appendices

## A Model computations

### A.1 Computation of the PSE tension

To show the equivalence of the two expression of the equilibrium tension $\tau^{th}(s)$, let us compute the derivation of the equilibrium attached free energy $\mathcal{F}_1^{th}(s)$.

$$ \frac{d\mathcal{F}_1^{th}}{ds}(s) = \int \partial_s w_1(s, y) p_1^{th}(y; s) dy + \underbrace{\int \partial_s p_1^{th}(y; s) [w_1(s, y) + k_B T(\ln(a p_1^{th}(y; s)) + 1)] dy}_{\equiv A(s)}. \quad (27) $$

Defining $Z_1(s) = \int \exp(-w_1(s, y)/(k_B T)) dy$, we have

$$ \begin{aligned} A(s) &= \int \left[ \left[ -\frac{1}{k_B T} \partial_s w_1(s, y) \frac{e^{-w_1(s, y)/(k_B T)}}{Z_1(s)} + \frac{1}{k_B T} \frac{e^{-w_1(s, y)/(k_B T)}}{Z_1(s)} \int \partial_s w_1(s, y') \frac{e^{-w_1(s, y')/(k_B T)}}{Z_1(s)} dy' \right] \right. \\ & \qquad \left. \cdot \left[ w_1(s, y) + k_B T \left( -\frac{w_1(s, y)}{k_B T} - \ln(Z_1(s)/a) + 1 \right) \right] \right] dy \\ &= k_B T \left( 1 - \ln(Z_1(s)/a) \right) \left[ - \int \partial_s w_1(s, y) \frac{e^{-w_1(s, y)/(k_B T)}}{Z_1(s)} dy \right. \\ & \qquad \left. + \int \partial_s w_1(s, y') \frac{e^{-w_1(s, y')/(k_B T)}}{Z_1(s)} dy' \underbrace{\int p_1^{th}(y, s) dy}_{=1} \right] \\ &= 0. \end{aligned} $$

From the expression (1) of $w_1$, we thus finally obtain

$$ \tau^{th}(s) = \int \kappa(y + s) p_1^{th}(y; s) dy = \int \partial_s w_1(s, y) p_1^{th}(y; s) dy = \frac{d\mathcal{F}_1^{th}}{ds}(s). $$

Naturally, a calculation with the free energy $\mathcal{F}(s) = -k_B T \ln(Z(s))$, the partition function defined as $Z(s) = \int e^{-w_1(s, y)/(k_B T)} dy$ and a definition of $\tau^{th}$ as $\tau^{th}(s) = -\partial_s \mathcal{F}(s)$ leads to the same result.

Note that similarly the equilibrium average stiffness per head can be computed:

$$ \kappa_c^{th}(t) = \frac{1}{d_a} \int_{s-}^{s^+} P_1(y, t) \partial_s \tau^{th}(s) ds. $$

### A.2 Generalized Huxley’s 57 model

The PSE model can be seen as an instance of a larger family of models within the Huxley’s 57 framework. These models are characterized by only two states ($\alpha = 0, 1$) associated with the free energies $\mathcal{F}_{0,1}$ and the transition rates $f$ and $g$, which are the inputs of the model. The probability of being attached is ruled by the PDE

$$ \partial_t P_1(s, t) + \dot{x}_c \partial_s P_1(s, t) = f(s)(1 - P_1(s, t)) - g(s)P_1(s, t), \quad (28) $$
and the active tension is given by

$$ \tau_c(t) = \frac{1}{d_a} \int \frac{d\mathcal{F}_1}{ds}(s) P_1(s,t) ds. \qquad (29) $$

We see that our PSE model pertains to this larger family with the inputs derived from a more refined model and not prescribed directly as parameters.

Finally, we mention that a more widely used modeling approach can be retrieved by replacing the continuous internal degree of freedom y by a discrete variable i that is defined only in the attached state, the detached states being as before characterized by constant (space independent) energy levels [Hill, 1977; Eisenberg et al., 1980]. In that case, the internal energy of the attached state i is usually written as

$$ w_1(s, i) = \frac{\kappa_{\text{xb}}}{2} (s + s_i)^2 + u_1(i). $$

The original Huxley'57 model has only one of such states.

In this framework, the adiabatic elimination of the fast variable i leads to the definition of the equilibrium probabilities

$$ a^{th}(i; s) = \frac{\exp[-w_1(s, i)/(k_B T)]}{\sum_i \exp[-w_1(s, i)/(k_B T)]}, $$

and the resulting two-state model is then characterized by the free energy

$$ \mathcal{F}_1^{th}(s) = \sum_i \left[ w_1(s, i) a^{th}(y; s) + k_B T a^{th}(i; s) \ln(a \cdot a^{th}(i; s)) \right], $$

as an analog to (8b) and the tension

$$ \tau^{th}(t) = \sum_i \kappa_{\text{xb}}(s + s_i) a^{th}(i; s) $$

as an analog to (13).

## A.3 Steady-state tension in the non-linear macro-model

We recall that in the non-linear macro model

$$ T_c(t) = \sum_{i=0}^{n-1} C_{i+1} M_i(t) \qquad (30) $$

where $C_i = \rho_{\text{surf}}c_i$. The dynamics of this force then follows

$$ \dot{T}_c(t) = \sum_{i=0}^{n-1} C_{i+1} \dot{M}_i(t). $$

In the permanent regime, and using the moment dynamics (17), we obtain the relation

$$ 0 = C_1 \left[ -(f+g)M_0^\infty + f_0 \right] + \sum_{i=1}^{n-1} C_{i+1} \left[ i\dot{x}_c M_{i-1}^\infty - (f+g)M_i^\infty + f_i \right], $$

where $M^\infty$ is the steady-state value of the moments given by the recursive relation deduced from (17). They are given by

$$ \left\{ \begin{array}{l} M_0^\infty = \frac{f_0}{f+g}, \\ M_i^\infty(\dot{x}_c) = \frac{i}{f+g} \left[ \dot{x}_c M_{i-1}^\infty(\dot{x}_c) + f_i \right], \text{ for } i \ge 1. \end{array} \right. \qquad (31) $$

Applying recursively (31), we obtain a general form of the steady-state moment value

$$ M_i^\infty(\dot{x}_c) = \sum_{j=0}^{i} \frac{\frac{i!}{j!} f_i}{(f+g)^{i+1-j}} \dot{x}_c^{i-j}. $$

Noting finally that the steady-state tension is defined from (22) by $T_c^\infty = \sum_{i=0}^{n-1} C_{i+1} M_i^\infty$, we obtain equation (23).
# B Calibration

## B.1 Reference calibration

The double quadratic well potentials $u_\alpha$ are presented in Figure 12. Combining them with the myosin neck potential energy, we obtain the energy landscapes $w_\alpha$ of the attached and detached states, which are depicted in Figure 13. The transition rates $k_+$ and $k_-$ are presented in Figure 14.

From these calibration choices, we derive the constitutive elements of the PSE model. The energy landscapes $u_\alpha$ allow to compute the equilibrium probability densities $p_0^{th}$ and $p_1^{th}$. Then, the integration against the transition rates $k_+$ and $k_-$ leads to the derivation of the thermal equilibrium transition rates $f^{th}$ and $g^{th}$. The elements of the PSE model are presented in Figure 15.

<figure><img src="image_13.png" /><figcaption>Figure 12: Double well potential in both attached ($u_1$) and detached state ($u_0$). The value of ATP chemical potential $\mu_T$ is 100 zJ</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13: Contour lines of the myosin head energy landscape. The thin dashed line represents the separation between the pre-power stroke and post-power stroke conformations. (Left) Detached potential $w_0(x, y)$. (Right) Attached potential $w_1(s, y)$. The solid black line of equation $s + y = 0$ separates the regions where attached heads develop positive (traction) and negative (compression) forces.</figcaption></figure>
<table><caption>Table 5: Calibration of the model for cardiac data.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td colspan="3">Power stroke potentials</td></tr><tr><td>Power stroke characteristic length</td><td>a</td><td>11 nm</td></tr><tr><td>Bistable element in attached state (α = 1)</td><td></td><td></td></tr><tr><td rowspan="3">u<sub>1</sub>(y) = { κ<sub>1pre</sub>/2(y - y<sub>1pre</sub>)<sup>2</sup> + v<sub>1</sub> if y &lt; ℓ<sub>1</sub>,<br/>κ<sub>1post</sub>/2(y - y<sub>1post</sub>)<sup>2</sup> otherwise,</td><td>κ<sub>1pre</sub></td><td>5.60 pN nm<sup>-1</sup></td></tr><tr><td>κ<sub>1post</sub></td><td>1.40 pN nm<sup>-1</sup></td></tr><tr><td>ℓ<sub>1</sub></td><td>1.42 nm</td></tr><tr><td rowspan="2">v<sub>1</sub> = κ<sub>1post</sub>/2(ℓ<sub>1</sub> - y<sub>1post</sub>)<sup>2</sup> - κ<sub>1pre</sub>/2(ℓ<sub>1</sub> - y<sub>1pre</sub>)<sup>2</sup></td><td>y<sub>1post</sub></td><td>a</td></tr><tr><td>y<sub>1pre</sub></td><td>0</td></tr><tr><td>Bistable element in detached state (α = 0)</td><td></td><td></td></tr><tr><td rowspan="3">u<sub>0</sub>(y) = { κ<sub>0pre</sub>/2(y - y<sub>0pre</sub>)<sup>2</sup> + v<sub>0</sub> + E if y &lt; ℓ<sub>0</sub>,<br/>κ<sub>0post</sub>/2(y - y<sub>0post</sub>)<sup>2</sup> + E otherwise,</td><td>κ<sub>0pre</sub></td><td>11.5 pN nm<sup>-1</sup></td></tr><tr><td>κ<sub>0post</sub></td><td>0.45 pN nm<sup>-1</sup></td></tr><tr><td>ℓ<sub>0</sub></td><td>1.42 nm</td></tr><tr><td rowspan="3">v<sub>0</sub> = κ<sub>0post</sub>/2(ℓ<sub>0</sub> - y<sub>0post</sub>)<sup>2</sup> - κ<sub>0pre</sub>/2(ℓ<sub>0</sub> - y<sub>1pre</sub>)<sup>2</sup></td><td>y<sub>0post</sub></td><td>6 nm</td></tr><tr><td>y<sub>0pre</sub></td><td>0</td></tr><tr><td>E</td><td>80 zJ</td></tr><tr><td colspan="3">Energy landscapes w<sub>α</sub>(y) = u<sub>α</sub>(y) + <sup>1</sup>&frasl;<sub>2</sub>κ(x + y)<sup>2</sup> (Figure 13)</td></tr><tr><td rowspan="4">u<sub>α</sub>(y) = u&#773;<sub>α</sub>(y + s&#773;<sub>α</sub>)</td><td>κ</td><td>1.34 pN nm<sup>-1</sup></td></tr><tr><td>s&#773;<sub>0</sub></td><td>1.2 nm</td></tr><tr><td>s&#773;<sub>1</sub></td><td>1.2 nm</td></tr><tr><td>μ<sub>T</sub></td><td>100 zJ</td></tr><tr><td colspan="3">Stochastic dynamics</td></tr><tr><td>Drag coefficient</td><td>η</td><td>0.0972 ms pN nm<sup>-1</sup></td></tr><tr><td>Microscopic timescale</td><td>γ = η/κ</td><td>0.0725 ms</td></tr><tr><td>Temperature</td><td>T</td><td>298 K</td></tr><tr><td colspan="3">Attachment / detachment rates (Figure 14)</td></tr><tr><td>ψ<sub>0α</sub>(y) = 1/2 {1 + tanh[λ<sub>1</sub>(ℓ<sub>α</sub> - y)]}, ψ<sub>1α</sub> = 1 - ψ<sub>0α</sub></td><td>λ<sub>1</sub></td><td>7.27 nm<sup>-1</sup></td></tr><tr><td>k<sub>+</sub>(s, y) = ψ<sub>00</sub>(y) k&#773;<sub>+</sub>(s)</td><td>k<sub>max</sub></td><td>390 s<sup>-1</sup></td></tr><tr><td>k<sub>-</sub>(s, y) = ψ<sub>01</sub>(y) k&#773;<sub>pre</sub>(s, y) + ψ<sub>11</sub>(y) k&#773;<sub>post</sub>(s, y) + k&#773;<sub>*</sub>(s, y)</td><td>λ<sub>2</sub></td><td>1.60 nm<sup>-1</sup></td></tr><tr><td rowspan="5">k&#773;<sub>+</sub>(s) = <sup>k<sub>max</sub></sup>&frasl;<sub>2</sub> [tanh(λ<sub>2</sub>(s + ℓ<sub>+</sub>)) 1<sub>1</sub><sub>-∞,0</sub>(s) + tanh(λ<sub>2</sub>(s - ℓ<sub>+</sub>)) 1<sub>1</sub><sub>0,∞</sub>(s)]</td><td>ℓ<sub>+</sub></td><td>3.5 nm</td></tr><tr><td>k<sub>0</sub></td><td>1400 s<sup>-1</sup></td></tr><tr><td>λ<sub>3</sub></td><td>5 nm<sup>-1</sup></td></tr><tr><td>λ<sub>4</sub></td><td>5 nm<sup>-1</sup></td></tr><tr><td>ℓ<sub>-,l</sub></td><td>-6.5 nm</td></tr><tr><td>k&#773;<sub>post</sub>(s, y) = k<sub>0</sub> exp(-λ<sub>3</sub>(s + y - ℓ<sub>-,l</sub>))</td><td>ℓ<sub>-,r</sub></td><td>9 nm</td></tr><tr><td>k&#773;<sub>pre</sub>(s, y) = k<sub>0</sub> exp(λ<sub>4</sub>(s + y - ℓ<sub>-,r</sub>))</td><td>k<sub>min</sub></td><td>66 s<sup>-1</sup></td></tr><tr><td rowspan="3">k&#773;<sub>*</sub>(s, y) = k<sub>min</sub> + <sup>k<sub>*</sub></sup>&frasl;<sub>2</sub> [2 + tanh(λ<sub>*</sub>(s - ℓ<sub>*</sub>)) - tanh(λ<sub>*</sub>(s + ℓ<sub>*</sub>))]</td><td>k<sub>*</sub></td><td>162 s<sup>-1</sup></td></tr><tr><td>λ<sub>*</sub></td><td>72.7 nm<sup>-1</sup></td></tr><tr><td>ℓ<sub>*</sub></td><td>4.5 nm<sup>-1</sup></td></tr><tr><td colspan="3">Geometric parameter</td></tr><tr><td>Reference length of a half sarcomere</td><td>ℓ<sub>hs</sub></td><td>0.925 µm</td></tr><tr><td>Lower bound of the reachable actin sites interval</td><td>s<sup>-</sup></td><td>-30 nm</td></tr><tr><td>Upper bound of the reachable actin sites interval</td><td>s<sup>+</sup></td><td>10 nm</td></tr></tbody></table>
<figure><img src="image_15.png" /><figcaption>Figure 14: Transition rates $k_+(s, y)$ and $k_-(s, y)$. The thin dashed line represents the separation between the pre-power stroke and post-power stroke conformations. (Left) Attachment rate $k_+$ (color map) and detached energy landscape (contour line). (Right) Detachment rate $k_-$ (color map) and attached energy landscape (contour line). The solid black line of equation $s + y = 0$ separates the regions where attached heads develop positive (traction) and negative (compression) forces.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 15: Constitutive elements of the PSE model. (a) Equilibrium energy levels. (b) Equilibrium average tension. (c) Equilibrium transition rates.</figcaption></figure>

<table><caption>Table 6: Model parameters of the macro-models.</caption><tbody><tr><th colspan="2">T<sub>1</sub>-linear macro-model</th></tr><tr><td>κ<sub>xb</sub></td><td>1.34 pN nm<sup>-1</sup></td></tr><tr><td>s<sub>0</sub></td><td>4.59 nm</td></tr><tr><th colspan="2">T<sub>2</sub>-linear macro-model</th></tr><tr><td>κ<sub>xb</sub></td><td>0.52 pN nm<sup>-1</sup></td></tr><tr><td>s<sub>0</sub></td><td>12.15 nm</td></tr><tr><th colspan="2">Nonlinear macro-model</th></tr><tr><td>c<sub>0</sub></td><td>0.349 κa<sup>2</sup></td></tr><tr><td>c<sub>1</sub></td><td>0.534 κa<sup>2</sup>/a</td></tr><tr><td>c<sub>2</sub></td><td>0.155 κa<sup>2</sup>/a<sup>2</sup></td></tr><tr><td>c<sub>3</sub></td><td>-0.286 κa<sup>2</sup>/a<sup>3</sup></td></tr><tr><td>c<sub>4</sub></td><td>0.213 κa<sup>2</sup>/a<sup>4</sup></td></tr><tr><td>c<sub>5</sub></td><td>0.234 κa<sup>2</sup>/a<sup>5</sup></td></tr><tr><td>c<sub>6</sub></td><td>0.0511 κa<sup>2</sup>/a<sup>6</sup></td></tr></tbody></table>
## B.2 Asymptotic calculation of the $T_1$-curve

At the time scale of the phase I (assuming that this phase is instantaneous), the force predicted by the stochastic model (4) becomes

$$ T_1(\delta s) = \frac{\rho_{\text{surf}}}{d_a} \int_{s^{-}-\infty}^{s^{+}+\infty} \int_{s^{-}-\infty}^{s^{+}+\infty} \kappa(s + \delta s + y) \tilde{p}(y; s) \, dy ds = T_0 + \kappa \rho_{\text{surf}} \delta s \frac{1}{d_a} \int_{s^{-}-\infty}^{s^{+}+\infty} \int_{s^{-}-\infty}^{s^{+}+\infty} \tilde{p}(y; s) \, dy ds, \quad (32) $$

where $\tilde{p}(y; s)$ is the isometric probability density of being attached and is the steady-state solution of (3) with $\dot{x}_c = 0$. By definition, we have

$$ \frac{1}{d_a} \int_{s^{-}-\infty}^{s^{+}+\infty} \int_{s^{-}-\infty}^{s^{+}+\infty} \tilde{p}(y; s) \, dy ds = \tilde{n}_{\text{att}}, $$

and thus the tension $T_1$ is given by

$$ T_1(\delta s) = T_0 + \kappa \rho_{\text{surf}} \tilde{n}_{\text{att}} \delta s. $$

## B.3 Asymptotic calculation of the $T_2$-curve

We want to establish the equation of the asymptotic branch of $\tau^{th}$ for large and negative $s$. We consider the energy landscape $w_1$ defined with the quadratic double well potential $u_1$ by $w_1(s, y) = \kappa/2(s+y)^2 + u_1(y)$. We have

$$ \begin{aligned} \tau^{th}(s) & \approx \int_{-\infty}^{\infty} \kappa(s+y) p_1^{th}(y; s) \, dy = \int_{-\infty}^{\infty} \kappa(s+y) \frac{e^{-w_1(s,y)/k_B T}}{\int_{-\infty}^{\infty} e^{-w_1(s,y')/k_B T} \, dy'} \, dy \\ &= \kappa s + \frac{\kappa}{\int_{-\infty}^{\infty} e^{-w_1(y',s)/k_B T} \, dy'} \int_{-\infty}^{\infty} ye^{-w_1(y,s)/k_B T} \, dy. \end{aligned} \quad (33) $$

The integral over the internal variable $y$ can be split into two parts, one for each well of the attached potential. For large and negative $s$, all the heads are concentrated in the post-power stroke well of the energy landscape (see Figure 13). Therefore, we can neglect the contribution of the pre-power stroke well to the integral and we approximate the whole bistable potential by the post power stroke well, and obtain

$$ \begin{aligned} \int_{-\infty}^{\infty} e^{-w_1(s,y')/k_B T} \, dy' & \approx \int_{\ell_1}^{\infty} e^{-w_1(s,y')/k_B T} \, dy' \\ & \approx \int_{-\infty}^{\infty} \exp \left[ - \frac{1}{k_B T} \left[ \frac{\kappa}{2}(s+y')^2 + \frac{\kappa_{1\text{post}}}{2}(y'-y_{1\text{post}})^2 \right] \right] \, dy'. \end{aligned} $$

Writing the quadratic potential in the canonical form leads to

$$ \begin{aligned} \int_{-\infty}^{\infty} e^{-w_1(s,y')/k_B T} \, dy' & \approx \int_{-\infty}^{\infty} \exp \left[ - \frac{\kappa + \kappa_{1\text{post}}}{2k_B T} (y' + \frac{\kappa s - \kappa_{1\text{post}}y_{1\text{post}}}{\kappa + \kappa_{1\text{post}}})^2 - \frac{1}{2k_B T} \frac{\kappa_{1\text{post}}\kappa}{\kappa_{1\text{post}} + \kappa} (s + y_{1\text{post}})^2 \right] \, dy' \\ &= \exp \left[ - \frac{1}{2k_B T} \frac{\kappa_{1\text{post}}\kappa}{\kappa_{1\text{post}} + \kappa} (s + y_{1\text{post}})^2 \right] \sqrt{\frac{2\pi k_B T}{\kappa + \kappa_{1\text{post}}}}. \end{aligned} $$
We define the stiffness $\bar{\kappa}$ and the two lengths $y'_0(s)$ and $L_{\text{ref}}$ by

$$ \bar{\kappa} = \frac{\kappa \kappa_{1\text{post}}}{\kappa + \kappa_{1\text{post}}}, \quad y'_0(s) = \frac{\kappa s - \kappa_{1\text{post}} y_{1\text{post}}}{\kappa + \kappa_{1\text{post}}}, \quad \text{and} \quad L_{\text{ref}} = \sqrt{\frac{k_B T}{\kappa + \kappa_{1\text{post}}}}. $$

Using the previous calculation and its result, (33) becomes

$$ \tau^{th}(s) \approx \kappa s + \frac{\kappa}{\exp\left[-\frac{1}{2k_B T} \bar{\kappa}(s + y_{1\text{post}})^2\right] \sqrt{2\pi L_{\text{ref}}}} \\ \cdot \exp\left[-\frac{1}{2k_B T} \bar{\kappa}(s + y_{1\text{post}})^2\right] \int_{-\infty}^{\infty} y \exp\left[-\frac{1}{2L_{\text{ref}}^2} (y + y'_0(s))^2\right] dy. $$

With the change of variable $u = \frac{1}{L_{\text{ref}}} (y + y'_0(s))$, we obtain

$$ \tau^{th}(s) \approx \kappa s + \frac{\kappa}{\sqrt{2\pi} L_{\text{ref}}} \int_{-\infty}^{\infty} (L_{\text{ref}}u - y'_{0}(s)) \exp\left[-\frac{u^2}{2}\right] L_{\text{ref}} du. $$

Noting that

$$ \int_{-\infty}^{\infty} u \exp\left[-\frac{u^2}{2}\right] du = 0, $$

and expanding $y'_0(s)$, we obtain

$$ \tau^{th}(s) \approx \kappa s - \frac{\kappa}{\sqrt{2\pi}} \sqrt{2\pi} \frac{\kappa s - \kappa_{1\text{post}} y_{1\text{post}}}{\kappa + \kappa_{1\text{post}}}. $$

Finally, the asymptotic expression of $\tau^{th}(s)$ for large and negative values of $s$ is

$$ \tau^{th}(s) \underset{s \to -\infty}{\sim} \frac{\kappa \kappa_{1\text{post}}}{\kappa + \kappa_{1\text{post}}} (s + y_{1\text{post}}). $$

Similarly, we obtain the asymptotic branch for large and positive values of $s$, neglecting the contribution of myosin heads in the post-power stroke conformation, as

$$ \tau^{th}(s) \underset{s \to +\infty}{\sim} \frac{\kappa \kappa_{1\text{pre}}}{\kappa + \kappa_{1\text{pre}}} (s + y_{1\text{pre}}). $$

In conclusion, the $T_2$-curve being given by

$$ T_2(\delta s) = \frac{\rho_{\text{surf}}}{d_a} \int_{s^{-}}^{s^{+}} \check{P}_1(s) \tau^{th}(s + \delta s) ds, $$

the its asymptotic slopes are given by

$$ \left| \frac{\partial T_2}{\partial \delta s} \right|_{\delta s \to -\infty} = \rho_{\text{surf}} \check{n}_{\text{att}} \frac{\kappa \kappa_{1\text{post}}}{\kappa + \kappa_{1\text{post}}}, \\ \left| \frac{\partial T_2}{\partial \delta s} \right|_{\delta s \to +\infty} = \rho_{\text{surf}} \check{n}_{\text{att}} \frac{\kappa \kappa_{1\text{pre}}}{\kappa + \kappa_{1\text{pre}}}. $$

## B.4 Near isometric behavior

The shape of the force-velocity curve in near isometric conditions – ie. for shortening with characteristic time $\ell_{hs}/\dot{x}_c$ that is small with respect to the transition rates – is determined by a balance between two effects: a change in the number of attached heads and a change in the averaged force per attached head. This balance can lead to an increase or a decrease of the force at slow sliding velocities with respect to
<figure><img src="image_17.png" /><figcaption>Figure 16: Investigation of the force-velocity curve in near isometric conditions obtained with the PSE model. We compare three values of $k_*$ – black line: $k_* = 0$, orange dashed line: $k_* = 0.25(k_{max} - k_{min})$ and blue dashed line: $k_* = 0.5(k_{max} - k_{min})$. (a) Transition rates for the considered value of $k_*$. (b) Global force-velocity curve with focus on the near isometric region of the force-velocity curve. (c) Ratio of attached heads as a function of the shortening velocity. (d) Force per attached head as a function of the shortening velocity.</figcaption></figure>
the isometric force. In our calibration, once all other parameter are fixed, the near isometric properties of the force-velocity curve are adjusted with the function $\bar{k}_*$ (see Table 5).

In our reference model, $\bar{k}_* = 0$ and therefore the detachment function $g^{th} = k_{min}$ is constant over the whole $[s^+, s^+]$ except near the boundary where it diverges. With this detachment function, we observe that the tension actually increases at slow shortening velocities compared to the isometric configuration, see Figure 5 and Figure 16(b, black line).

This effect may be linked to the observed change of curvature of the force-velocity relation at high load Edman [1988], as well as the oscillatory behavior following load small load perturbations in near isometric condition Edman & Curtin [2001]. We shown in Figure 16 that this effect can be modulated by the function $\bar{k}_*$ with only a marginal influence on the FV-curve.

The function $\bar{k}_*$ depends on three parameters: the position of the transition between the two detachment rate regions $\ell_*$, the reciprocal characteristic length $\lambda_*$ and the detachment rate value $k_*$.

We here analyze only the influence of $\bar{k}_*$, see Figure 16. Increasing the detachment rate reduces the fraction of attached heads (see Figure 16(c)) but it also increases the tension per attached heads (see Figure 16(d)). The first effect has to be stronger than the second to eliminate the increase of the overall tension at low shortening velocity. Furthermore, this increase of the detachment rate only marginally affects the rest of the force-velocity curve, as shown in Figure 16(b). Therefore this introduction of a s-dependence in the detachment rate can be adjusted only at the end of the calibration procedure. Note that similar effect can be obtained by adjusting the parameters $\ell_*$ and $\lambda_*$, also without affecting the global shape of the force-velocity curve.

The deviation of the force-velocity relation from the canonical hyperbola has been extensively studied, both experimentally Edman [1988]; Edman et al. [1997]; Edman & Curtin [2001] and theoretically Jülicher & Prost [1995, 1997]; Vilfan et al. [1999]; Duke [1999]; Guérin et al. [2011]; Månsson [2010, 2014, 2016]. In our work, we show that once all other parameter have been calibrated, the near isometric behavior can be adjusted by modifying the detachment rate only. In [Vilfan et al., 1999], the authors also show that the shape of the force-velocity curve close to the stall force is very sensitive to the detachment rate. Other work have shown that this part of the force-velocity curve also depends on the cross-bridge stiffness (our parameter $\kappa$) [Månsson, 2014] or on the free energy difference between the two attached states [Månsson, 2010]. The first method does not seem to strongly affect the rest of the force velocity curve, as in our case, while the second method has a strong impact on the whole curve. We mention that another effect that might play a role is the viscous damping by the fluid surrounding the contractile apparatus Duke [1999]. This viscous term is not modeled in this work but is usually included in the behavior law of the tissue, and accounts for its passive visco-elastic response, see for instance [Chapelle et al., 2012; Caruel et al., 2014].

## B.5 Effect of the viscous drag on the fast transient response

In Figure 6(b), we show different $T_1(\delta\ell_{hs})$ obtained with the stochastic model for different values of the drag coefficient $\eta$.

As $\eta$ increases, the deviation from the purely elastic response (see black curve in Figure 6(b) and (24)) is reduced and $r_2$, the rate of phase II, decreases (see Figure 6(c)). Our reference value $\eta = 0.0972 \text{ ms pN nm}^{-1}$ was chosen to match the measured $T_1$ response while maintaining the order of magnitude of $r_2$ in accordance with experimental data [Caremani et al., 2016]. It must be noted however that our estimate of the rate of phase II cannot be directly compared with experimental data as it does not consider the effect of filament compliance, which has a strong impact on the rate of force recovery [Piazzesi et al., 2014].

## C Model comparison

In this section, we present a more detailed analysis of the validity of the fundamental assumption used to derive the PSE model from the stochastic model. This assumption supposes that the internal variables X and Y of the stochastic model are distributed according to the thermal equilibrium distribution for both the attached and the detached states. Equivalently this assumption can be characterized by the fact that the equilibration of the internal variables in their respective energy potential occurs within a time scale that is much faster than any other time scale in the system.
To evaluate validity of this assumption, we thus need to compare the distribution $p(x, y, 0; s, t)/P_0(s, t)$ and $\bar{p}(y; s, t)/P_1(s, t)$ with their respective counterparts for the PSE model $p_0^{th}(x, y)$ and $p_1^{th}(y; s)$. We recall that we have the following definitions

$$P_0(s, t) = \iint p(x, y, 0; s, t)dydx,$$

$$P_1(s, t) = \int \bar{p}(y; s, t)dy.$$

We compare the distributions of the two models in steady-state in isometric conditions and at two different shortening velocities $\dot{x}_c$. The results are presented in Figure 17 for attached state and in Figures 18, 19 and 20 for the detached state.

We can first notice that the two sets of distributions display a very good match in isometric conditions for both the attached and the detached state. Moreover, the distributions also show a good consistency in the detached state at all considered sliding velocities. This justifies the validity of our approach to derive the PSE model. However, for the attached state, a few notable differences appear when the sliding velocity increases.

First, for values of the parameter $s$ approximately in the range $[-2\,\text{nm}, 2\,\text{nm}]$, myosin heads are more likely to be in the pre-power stroke configuration (and thus less likely to be in the post-power stroke configuration) in the stochastic model compared to the PSE model. This is the result of the competition between the time scale of the transition from the pre-power stroke configuration to the post-power stroke configuration, which implies to overcome the energy barrier, and the time scale of the filament sliding. With increasing velocity, heads are more quickly moved to lower values of $s$ and thus have less time to accomplish the power stroke. Therefore they are more likely to remain in the pre-power stroke configuration. For values of $s$ outside of the range $[-2\,\text{nm}, 2\,\text{nm}]$, the energy barrier is lower and therefore this effect is less present.

Second, the distributions also slightly differs for low values of the parameter $s$ ($s < -20\,\text{nm}$) with a distribution for the variable $y$ that can reach lower value in the PSE model (we present the results for the value $s = -24\,\text{nm}$ in Figure 17 as an example). This difference is due to the fact that in the stochastic model, heads with these low values of $y$ are subjected to a very high detachment rate so their probability of having this configuration is reduced. This is not the case with the PSE model since the equilibration process is assumed to be much faster than any attachment-detachment process.

We will now see how the comparison of the distributions can be used to explain the discrepancy observed between the force-velocity relations predicted by the two models. Let us first analyze the attached distributions (see Figure 17). In the first distribution discrepancy between the stochastic and the PSE model, the configuration of the myosin heads in the PSE model is shifted from the pre-power stroke configuration to the post-power stroke configuration, which is associated with a higher force. Therefore, this discrepancy contributes to the fact that the developed force is higher in the PSE model than in the stochastic model.

In the second distribution discrepancy, lower values of the variable $y$ are more likely for the PSE model than for the stochastic model. However, since the detachment rate $k_-$ diverges in this region of $y$, the small difference in the distribution has a very large impact on the averaged transition rate $g^{th}$ and explains the differences in the shape of the ratio of attached heads (see Figure 8) and ultimately on the force-velocity relation at large sliding velocities.

We now analyze the detached distributions. At the first glance no noticeable differences appear between the stochastic and the PSE models. However, we have seen in Figure 8 that the shape of isometric ratio of attached heads differ slightly between the two models, which we interpret as a result of the competition between the attachment rate and the rate for returning to the pre-power stroke configuration in the detached state for cycling heads. To analyze this effect in more detail, we estimate the average transition time from the post-power stroke configuration to the pre-power stroke configuration in the detached state. The results are displayed in Figure 21(a). The average transition time is 1.7 ms which is not much faster than the inverse of the attachment rate $1/k_{max} = 2.6$ ms explaining the discrepancy in the ratio of attached heads $P_1^{\infty}$ but still shorter which explains that differences do not appear clearly in the distributions. Note that we also estimated the transition time from the pre-power stroke to the post-power stroke configuration in the attached state (see Figure 21(b)) and we obtain an average transition time of 0.68 ms.
<figure><img src="image_18.png" /><figcaption>Figure 17: Comparison of the distribution of the variable *y* in the attached state between the stochastic and the PSE models in steady-state conditions at various sliding velocities $\dot{x}_c$. The estimations of the distribution for the stochastic model are computed from a simulation using a population of 2 × 10⁶ myosin heads. For each sliding velocity, several values of the parameter *s* are displayed. Note that here we present the distribution of the variable *y* for heads that are attached, the actual ratio of attached heads can be read in Figure 8. Note also that for ($\dot{x}_c$ = 0 µm s⁻¹, *s* = −25 nm), the stochastic distribution cannot be evaluated because not heads are attached for this value of the parameter *s*.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 18: Comparison of the distribution of the variables *x*, *y* in the detached state between the stochastic and the PSE models in steady-state conditions at the sliding velocities $\dot{x}_c = 0 \mu m s^{-1}$. The estimations of the distribution for the stochastic model are computed from a simulation using a population of 2 × 10⁶ myosin heads. For each sliding velocity, several values of the parameter *s* are displayed. Note that here we present the distribution of the variable *y* for heads that are attached, the actual ratio of attached heads can be read as $P_0^\infty(s) = 1 - P_1^\infty(s)$ in Figure 8.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 19: Comparison of the distribution of the variables *x*, *y* in the detached state between the stochastic and the PSE models in steady-state conditions at the sliding velocities $\dot{x}_c = -3 \mu m s^{-1}$. The estimations of the distribution for the stochastic model are computed from a simulation using a population of 2 × 10⁶ myosin heads. For each sliding velocity, several values of the parameter *s* are displayed. Note that here we present the distribution of the variable *y* for heads that are attached, the actual ratio of attached heads can be read as $P_0^\infty(s) = 1 - P_1^\infty(s)$ in Figure 8.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 20: Comparison of the distribution of the variables *x*, *y* in the detached state between the stochastic and the PSE models in steady-state conditions at the sliding velocities $\dot{x}_c = -6 \mu m s^{-1}$. The estimations of the distribution for the stochastic model are computed from a simulation using a population of 2 × 10⁶ myosin heads. For each sliding velocity, several values of the parameter *s* are displayed. Note that here we present the distribution of the variable *y* for heads that are attached, the actual ratio of attached heads can be read as $P_0^\infty(s) = 1 - P_1^\infty(s)$ in Figure 8.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure 21: Estimation of the configuration transition time in the stochastic model for both the attached and detached state. (a) Detached state. Myosin heads are initially in the post-power stroke configuration and we register over time the ratio of the population that has transition to the post power stroke configuration. The average transition time is 1.7 ms (b) Attached state. Myosin heads are initially in the pre-power stroke configuration and we register over time, for various values of the parameter *s*, the ratio of the population that has transition to the post power stroke configuration. The average transition time is 0.68 ms</figcaption></figure>

D Summary of the main symbols used in the paper
<table><thead><tr><th>Symbol</th><th>Definition</th><th>First occurrence</th></tr></thead><tbody><tr><td>X<sup>t</sup>, x</td><td>location of the head tip</td><td>p.4</td></tr><tr><td>Y<sup>t</sup>, y</td><td>internal dof of the head (power stroke)</td><td>p.4</td></tr><tr><td>s</td><td>position of the nearest actin binding site</td><td>p.4</td></tr><tr><td>α<sup>t</sup></td><td>attachment state of the head (0 or 1)</td><td>p.4</td></tr><tr><td>w<sub>α</sub>, w<sub>α</sub><sup>L</sup>, w<sub>α</sub><sup>NL</sup></td><td>energy of the attached state</td><td>(1), p.8, p.9</td></tr><tr><td>F</td><td>free energy level</td><td>p.29</td></tr><tr><td>F<sup>th</sup></td><td>free energy level for the PSE model</td><td>(8)</td></tr><tr><td>u</td><td>double well potential</td><td>(1)</td></tr><tr><td>p(x, y, α; s, t)</td><td>probability function for (x, y, α) (densities in (x, y), discrete in α)</td><td>p.5</td></tr><tr><td>p̄</td><td>effective probability in the attached state</td><td>p.5</td></tr><tr><td>p<sub>α</sub><sup>th</sup></td><td>equilibrium probability distribution of the fast internal variables (x and y)</td><td>(7)</td></tr><tr><td>d<sub>a</sub></td><td>distance between two consecutive actin sites</td><td>(4)</td></tr><tr><td>P<sub>α</sub></td><td>population probability</td><td>p.6</td></tr><tr><td>dot_c</td><td>filaments relative sliding velocity</td><td>(2)</td></tr><tr><td>η</td><td>drag coefficient</td><td>(2)</td></tr><tr><td>k<sub>+</sub></td><td>attachment rate of the stochastic model</td><td>(2)</td></tr><tr><td>k<sub>-</sub></td><td>detachment rate of the stochastic model</td><td>(2)</td></tr><tr><td>f<sup>th</sup></td><td>attachment rate of the PSE model</td><td>(9)</td></tr><tr><td>g<sup>th</sup></td><td>detachment rate of the PSE model</td><td>(9)</td></tr><tr><td>f</td><td>attachment rate of the Huxley'57 like models</td><td>(28)</td></tr><tr><td>g</td><td>detachment rate of the Huxley'57 like models</td><td>(28)</td></tr><tr><td>S<sub>f</sub></td><td>support of the attachment rate f<sup>th</sup> or f</td><td>p.10</td></tr><tr><td>d<sub>f</sub></td><td>width of S<sub>f</sub></td><td>p.12</td></tr><tr><td>f̄<sup>th</sup></td><td>value of f<sup>th</sup> on S<sub>f</sub></td><td>p.12</td></tr><tr><td>ḡ<sup>th</sup></td><td>value of g<sup>th</sup> on S<sub>f</sub></td><td>p.12</td></tr><tr><td>f̄</td><td>value of f on S<sub>f</sub></td><td>p.22</td></tr><tr><td>ḡ</td><td>value of g on S<sub>f</sub></td><td>p.22</td></tr><tr><td>n̄<sub>att</sub></td><td>ration of attached heads in isometric conditions</td><td>p.12</td></tr><tr><td>l<sub>hs</sub></td><td>length of the half-sarcomeres in the reference configuration</td><td>p.6</td></tr><tr><td>ρ<sub>surf</sub></td><td>surface density of myosin heads in a longitudinal portion of a cardiomyocyte of thickness l<sub>hs</sub></td><td>(5)</td></tr><tr><td>τ<sub>c</sub></td><td>force per myosin heads</td><td>p.4</td></tr><tr><td>τ<sup>th</sup></td><td>average equilibrium tension</td><td>(12)</td></tr><tr><td>τ<sub>c</sub><sup>th</sup></td><td>force per myosin heads for the PSE model</td><td>p.7</td></tr><tr><td>τ̄<sub>c</sub><sup>th</sup></td><td>force per myosin heads for the PSE model in isometric conditions</td><td>p.12</td></tr><tr><td>κ<sub>c</sub><sup>th</sup></td><td>stiffness per myosin heads for the PSE model</td><td>p.29</td></tr><tr><td>κ̄<sub>c</sub><sup>th</sup></td><td>stiffness per myosin heads for the PSE model in isometric conditions</td><td>p.13</td></tr><tr><td>T<sub>c</sub></td><td>macroscopic tension</td><td>(5)</td></tr><tr><td>T<sub>0</sub></td><td>macroscopic tension in isometric conditions</td><td>p.10</td></tr><tr><td>T<sub>c</sub><sup>∞</sup></td><td>macroscopic tension in steady-state shortening</td><td>(21)</td></tr><tr><td>T<sub>1</sub></td><td>macroscopic tension at the end of phase I in fast transient experiments</td><td>p.12</td></tr><tr><td>T<sub>2</sub></td><td>macroscopic tension at the end of phase II in fast transient experiments</td><td>p.12</td></tr><tr><td>K<sub>c</sub></td><td>macroscopic stiffness</td><td>p.9</td></tr><tr><td>K<sub>∞</sub></td><td>macroscopic stiffness in the rigor state</td><td>(20)</td></tr><tr><td>L<sub>1</sub></td><td>shortening per half-sarcomere at the end of phase I in fast transient experiments</td><td>p.12</td></tr><tr><td>L<sub>2</sub></td><td>shortening per half-sarcomere at the end of phase II in fast transient experiments</td><td>p.12</td></tr><tr><td>μ<sub>T</sub></td><td>ATP chemical potential</td><td>p.14</td></tr><tr><td>J<sub>ATP</sub></td><td>ATP consumption flux per myosin head</td><td>p.27</td></tr><tr><td>ε̄</td><td>flux of internal energy per myosin head</td><td>p.20</td></tr><tr><td>W̄</td><td>flux of work produced per myosin head</td><td>p.20</td></tr><tr><td>M<sub>p</sub></td><td>moment of the population probability P<sub>1</sub></td><td>p.8</td></tr><tr><td>M<sup>∞</sup></td><td>moment of the population probability P<sub>1</sub> in steady-state</td><td>p.30</td></tr></tbody></table>
References

Alberts, B. (2015). *Molecular Biology of the Cell*. Garland Science, Taylor and Francis Group, 6th ed.

Allen, D.G., & Kentish, J.C. (1985). The cellular basis of the length-tension relation in cardiac muscle. *Journal of molecular and cellular cardiology, 17*(9), 821-840.

Amiad Pavlov, D., & Landesberg, A. (2016). The cross-bridge dynamics is determined by two length-independent kinetics: Implications on muscle economy and Frank-Starling Law. *Journal of molecular and cellular cardiology, 90*, 94-101.

Barclay, C.J. (2015). Energetics of contraction. *Comprehensive Physiology, 5*(2), 961-995.

Barclay, C.J., Woledge, R.C., & Curtin, N.A. (2010). Inferring crossbridge properties from skeletal muscle energetics. *Prog. Biophys. Molec. Biol., 102*(1), 53-71.

Bestel, J., Clément, F., & Sorine, M. (2001). A Biomechanical Model of Muscle Contraction. In *Medical Image Computing and Computer-Assisted Intervention – MICCAI 2001*, (pp. 1159–1161). Berlin, Heidelberg: Springer, Berlin, Heidelberg.

Caremani, M., Melli, L., Dolfi, M., Lombardi, V., & Linari, M. (2015). Force and number of myosin motors during muscle shortening and the coupling with the release of the ATP hydrolysis products. *The Journal of Physiology, 593*(15), 3313-3332.

Caremani, M., Pinzauti, F., Reconditi, M., Piazzesi, G., Stienen, G.J.M., Lombardi, V., & Linari, M. (2016). Size and speed of the working stroke of cardiac myosin in situ. *Proceedings of the National Academy of Sciences, 113*(13), 3675-3680.

Caruel, M., Allain, J.-M., & Truskinovsky, L. (2013). Muscle as a metamaterial operating near a critical point. *Physical Review Letters, 110*(24), 248103.

Caruel, M., Chabiniok, R., Moireau, P., Lecarpentier, Y., & Chapelle, D. (2014). Dimensional reductions of a cardiac model for effective validation and calibration. *Biomech. Model Mechanobiol., 13*(4), 897-914.

Caruel, M., Moireau, P., & Chapelle, D. (2019). Stochastic modeling of chemical-mechanical coupling in striated muscles. *Biomechanics and modeling in mechanobiology*, (pp. 1-25).

Caruel, M., & Truskinovsky, L. (2018). Physics of muscle contraction. *Reports on Progress in Physics, 81*(3), 036602.

Chabiniok, R., Moireau, P., Lesault, P.-F., Rahmouni, A., Deux, J.-F., & Chapelle, D. (2011). Estimation of tissue contractility from cardiac cine-MRI using a biomechanical heart model. *Biomechanics and modeling in mechanobiology, 11*(5), 609-630.

Chapelle, D., Le Tallec, P., Moireau, P., & Sorine, M. (2012). Energy-preserving muscle tissue model: formulation and compatible discretizations. *Journal for Multiscale Computational Engineering*.

Craig, R., & Padrón, R. (2004). Molecular structure of the sarcomere. *Myology, 3*, 129-144.

Daniels, M., Noble, M.I., ter Keurs, H.E.D.J., & Wohlfart, B. (1984). Velocity of sarcomere shortening in rat cardiac muscle: relationship to force, sarcomere length, calcium and time. *The Journal of Physiology, 355*, 367-381.

de Tombe, P.P. (2003). Cardiac myofilaments: mechanics and regulation. *Journal of Biomechanics, 36*(5), 721-730.

de Tombe, P.P., Mateja, R.D., Tachampa, K., Mou, Y.A., Farman, G.P., & Irving, T.C. (2010). Journal of Molecular and Cellular Cardiology. *Journal of molecular and cellular cardiology, 48*(5), 851-858.

de Tombe, P.P., & Stienen, G.J.M. (2007). Impact of temperature on cross-bridge cycling kinetics in rat myocardium. *The Journal of Physiology, 584*(2), 591-600.
de Tombe, P.P., & ter Keurs, H.E.D.J. (1990). Force and velocity of sarcomere shortening in trabeculae from rat heart. Effects of temperature. *Circulation Research, 66*(5), 1239-1254.

de Tombe, P.P., & ter Keurs, H.E.D.J. (1992). An internal viscous element limits unloaded velocity of sarcomere shortening in rat myocardium. *The Journal of Physiology, 454*(1), 619-642.

Dobesh, D.P., Konhilas, J.P., & de Tombe, P.P. (2002). Cooperative activation in cardiac muscle: impact of sarcomere length. *American Journal of Physiology-Heart and Circulatory Physiology, 282*(3), H1055-H1062.

Duke, T.A.J. (1999). Molecular model of muscle contraction. *Proceedings of the National Academy of Sciences, 96*(6), 2770-2775.

Edman, K.A.P. (1988). Double hyperbolic force velocity relation in frog muscle fibers. *J. Physiol., 404*, 301-321.

Edman, K.A.P., & Curtin, N.A. (2001). Synchronous oscillations of length and stiffness during loaded shortening of frog muscle fibres. *J. Physiol.*, **534**(265), 553-563.

Edman, K.A.P., Månsson, A., & Caputo, C. (1997). The biphasic force-velocity relationship in frog muscle fibres and its evaluation in terms of cross-bridge function. *The Journal of Physiology, 503*(1), 141-156.

Eisenberg, E., Hill, T.L., & Chen, Y. (1980). Cross-bridge model of muscle contraction. Quantitative analysis. *Biophysical Journal*, **29**(2), 195-227.

Ford, L.E., Huxley, A.F., & Simmons, R.M. (1981). The relation between stiffness and filament overlap in stimulated frog muscle fibres. *The Journal of physiology, 311*(1), 219-249.

Guérin, T., Prost, J., & Joanny, J. F. (2011). Dynamical behavior of molecular motor assemblies in the rigid and crossbridge models. *The European Physical Journal E, 34*(6), 60.

Hill, A.V. (1938). The heat of shortening and the dynamic constants of muscle. *Proc. R. Soc. Lond. B*, **126**(843), 136-195.

Hill, T.L. (1974). Theoretical formalism for the sliding filament model of contraction of striated muscle Part I. *Progress in biophysics and molecular biology*, **28**, 267-340.

Hill, T.L. (1976). Theoretical formalism for the sliding filament model of contraction of striated muscle part II. *Progress in biophysics and molecular biology*, **29**, 105-159.

Hill, T.L. (1977). *Free Energy Transduction in Biology*. Academic Press.

Holmes, K.C., & Geeves, M.A. (2000). The structural basis of muscle contraction. *Philosophical transactions of the Royal Society of London. Series B, Biological sciences*, **355**(1396), 419-431.

Huxley, A.F. (1957). Muscle structures and theories of contraction. *Progress in Biophysics and Biophysical Chemistry*.

Huxley, A.F., & Simmons, R.M. (1971). Proposed mechanism of force generation in striated muscle. *Nature*.

Jülicher, F., & Prost, J. (1995). Cooperative molecular motors. *Physical review letters*, **75**(13), 2618-2621.

Jülicher, F., & Prost, J. (1997). Spontaneous oscillations of collective molecular motors. *Physical review letters*, **78**(23), 4510-4513.

Kimmig, F., Chapelle, D., & Moireau, P. (2019). Thermodynamic properties of muscle contraction models and associated discrete-time principles. *Advanced Modeling and Simulation in Engineering Sciences*, **6**(1), 6.
Kobayashi, T., Jin, L., & de Tombe, P.P. (2008). Cardiac thin filament regulation. *Pflügers Archiv - European Journal of Physiology*, **457**(1), 37–46.

Konhilas, J.P., Irving, T.C., & deTombe, P.P. (2002). Myofilament calcium sensitivity in skinned rat cardiac trabeculae: role of interfilament spacing. *Circulation Research*, **90**(1), 59–65.

Kuhlman, B., & Bradley, P. (2019). Advances in protein structure prediction and design. *Nature Reviews Molecular Cell Biology*, (pp. 1–17).

Le Gall, A., Vallée, F., Mebazaa, A., Chapelle, D., Gayat, E., & Chabiniok, R. (2020). Monitoring of cardiovascular physiology augmented by a patient-specific biomechanical model during general anaesthesia. a proof of concept study. *PloS ONE*, *Public Library of Science, in press*.

Linari, M., Brunello, E., Reconditi, M., Fusi, L., Caremani, M., Narayanan, T., Piazzesi, G., Lombardi, V., & Irving, M. (2015). Force generation by skeletal muscle is controlled by mechanosensing in myosin filaments. *Nature*, **528**(7581), 276–279.

Linari, M., Dobbie, I., Reconditi, M., Koubassova, N., Irving, M., Piazzesi, G., & Lombardi, V. (1998). The stiffness of skeletal muscle in isometric contraction and rigor: The fraction of myosin heads bound to actin. *Biophysical Journal*, **74**(5), 2459–2473.

Lymn, R.W., & Taylor, E.W. (1971). Mechanism of adenosine triphosphate hydrolysis by actomyosin. *Biochemistry*, **10**(25), 4617–4624.

Månsson, A. (2010). Actomyosin-ADP states, interhead cooperativity, and the force-velocity relation of skeletal muscle. *Biophysical Journal*, **98**(7), 1237–1246.

Månsson, A. (2014). Hypothesis and theory: mechanical instabilities and non-uniformities in hereditary sarcomere myopathies. *Frontiers in physiology*, **5**, 350.

Månsson, A. (2016). Actomyosin based contraction: one mechanokinetic model from single molecules to muscle? *Journal of Muscle Research and Cell Motility*, **37**(6), 181–194.

Månsson, A., Rassier, D., & Tsiavaliaris, G. (2015). Poorly Understood Aspects of Striated Muscle Contraction. *BioMed Research International*, **2015**(1), 1–28.

Marcucci, L., & Truskinovsky, L. (2010). Mechanics of the power stroke in myosin II. *Physical Review E*, **81**(5), 051915.

Marcucci, L., Washio, T., & Yanagida, T. (2016). Including Thermal Fluctuations in Actomyosin Stable States Increases the Predicted Force per Motor and Macroscopic Efficiency in Muscle Modelling. *PLOS Computational Biology*, **12**(9), e1005083–20.

Pacher, P., Mabley, J.G., Liaudet, L., Evgenov, O.V., Marton, A., Haskó, G., Kollai, M., & Szabó, C. (2004). Left ventricular pressure-volume relationship in a rat model of advanced aging-associated heart failure. *AJP: Heart and Circulatory Physiology*, **287**(5), H2132–7.

Pertici, I., Bongini, L., Melli, L., Bianchi, G., Salvi, L., Falorsi, G., Squarci, C., Bozó, T., Cojoc, D., Kellermayer, M.S.Z., Lombardi, V., & Bianco, P. (2018). A myosin II nanomachine mimicking the striated muscle. *Nature Communications*, **9**(1), 3532.

Piazzesi, G., Dolfi, M., Brunello, E., Fusi, L., Reconditi, M., Bianco, P., Linari, M., & Lombardi, V. (2014). The myofilament elasticity and its effect on kinetics of force generation by the myosin motor. *Arch. Biochem. Biophys.*, **552-553**, 108–116.

Piazzesi, G., & Lombardi, V. (1995). A cross-bridge model that is able to explain mechanical and energetic properties of shortening muscle. *Biophysical Journal*, **68**, 1966–1979.

Piazzesi, G., Lucii, L., & Lombardi, V. (2002). The size and the speed of the working stroke of muscle myosin and its dependence on the force. *The Journal of Physiology*, **545**(1), 145–151.
Pinzauti, F., Pertici, I., Reconditi, M., Narayanan, T., S., Ger J.M., Piazzesi, G., Lombardi, V., Linari, M., & Caremani, M. (2018). The force and stiffness of myosin motors in the isometric twitch of a cardiac trabecula and the effect of the extracellular calcium concentration. *The Journal of Physiology, 596* (13), 2581-2596.

Rayment, I., Holden, H.M., Whittaker, M., Yohn, C.B., Lorenz, M., Holmes, K.C., & Milligan, R.A. (1993a). Structure of the actin-myosin complex and its implications for muscle contraction. *Science, 261*(5117), 58-65.

Rayment, I., Rypniewski, W.R., Schmidt-Base, K., Smith, R., Tomchick, D.R., Benning, M.M., Winkelmann, D.A., Wesenberg, G., & Holden, H.M. (1993b). Three-dimensional structure of myosin subfragment-1: a molecular motor. *Science, 261*(5117), 50-58.

Reconditi, M. (2006). Recent improvements in small angle x-ray diffraction for the study of muscle physiology. *Reports on Progress in Physics, 69*(23), 2709-2759.

Reconditi, M., Caremani, M., Pinzauti, F., Powers, J.D., Narayanan, T., Stienen, G.J.M., Linari, M., Lombardi, V., & Piazzesi, G. (2017). Myosin filament activation in the heart is tuned to the mechanical task. *Proceedings of the National Academy of Sciences*, (pp. 201619484-6).

Robert-Paganin, J., Auguin, D., & Houdusse, A. (2018). Hypertrophic cardiomyopathy disease results from disparate impairments of cardiac myosin function and auto-inhibition. *Nature communications*, (pp. 1 – 13).

Robert-Paganin, J., Pylypenko, O., Kikuti, C., Sweeney, H.L., & Houdusse, A. (2019). Force generation by myosin motors: A structural perspective. *Chemical Reviews*.

Sato, T., Shishido, T., Kawada, T., Miyano, H., Miyashita, H., Inagaki, M., Sugimachi, M., & Sunagawa, K. (1998). ESPVR of in situ rat left ventricle shows contractility-dependent curvilinearity. *The American journal of physiology, 274*(5), H1429-34.

Sermesant, M., Chabiniok, R., Chinchapatnam, P., Mansi, T., Billet, F., Moireau, P., Peyrat, J.-M., Wong, K., Relan, J., Rhode, K., Ginks, M., Lambiase, P., Delingette, H., Sorine, M., Rinaldi, C.A., Chapelle, D., Razavi, R., & Ayache, N. (2012). Patient-specific electromechanical models of the heart for the prediction of pacing acute effects in CRT: A preliminary clinical validation. *Medical Image Analysis, 16*(1), 201-215.

Silverthorn, D.U., Ober, W.C., Garrison, C.W., & Silverthorn, A.C. (2009). *Human physiology: an integrated approach*. Pearson.

Smith, D.A., Geeves, M.A., Sleep, J., & Mijailovich, S.M. (2008). Towards a unified theory of muscle contraction. I: foundations. *Annals of Biomedical Engineering, 36*(10), 1624-1640.

Smith, D.A., & Mijailovich, S.M. (2008). Towards a Unified Theory of Muscle Contraction. 2: Predictions with the Mean-Field Approximation. *Annals of Biomedical Engineering, 36*, 1353-1371.

Solaro, R.J., & Rarick, H. M. (1998). Troponin and tropomyosin: proteins that switch on and tune in the activity of cardiac myofilaments. *Circulation research, 83*(5), 471-480.

Ter Keurs, H.E.D.J., Rijnsburger, W.H., Van Heuningen, R., & Nagelsmit, M.J. (1980). Tension development and sarcomere length in rat cardiac trabeculae. evidence of length-dependent activation. *Circulation Research, 46*(5), 703-714.

Van Heuningen, R., Rijnsburger, W.H., & ter Keurs, H.E.D.J. (1982). Sarcomere length control in striated muscle. *The American journal of physiology, 242*(3), H411-20.

Vilfan, A., Frey, E., & Schwabl, F. (1999). Force-velocity relations of a two-state crossbridge model for molecular motors. *Europhysics Letters, 45*(3), 283-289.

Woody, M.S., Greenberg, M.J., Barua, B., Winkelmann, D.A., Goldman, Y.E., & Ostap, E.M. (2018). Positive cardiac inotrope omecamtiv mecarbil activates muscle despite suppressing the myosin working stroke. *Nature communications, 9*(1), 3838.
Zahalak, G.I. (1981). A distribution-moment approximation for kinetic theories of muscular contraction.
*Mathematical biosciences*, **55**(1-2), 89-114.

Zahalak, G.I. (2000). The two-state cross-bridge model of muscle is an asymptotic limit of multi-state models. *Journal of theoretical biology*, **204**(1), 67-82.