# A jump-diffusion stochastic formalism for muscle contraction models at multiple timescales

Louis-Pierre Chaintron, François Kimmig, Matthieu Caruel, Philippe Moireau

► To cite this version:

Louis-Pierre Chaintron, François Kimmig, Matthieu Caruel, Philippe Moireau. A jump-diffusion stochastic formalism for muscle contraction models at multiple timescales. Journal of Applied Physics, 2023, 10.1063/5.0158191 . hal-04264293

HAL Id: hal-04264293

https://hal.science/hal-04264293v1

Submitted on 30 Oct 2023

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 CC BY 4.0 - Attribution - International License
# A jump-diffusion stochastic formalism for muscle contraction models at multiple timescales

Louis-Pierre Chaintron<sup>1,2,3,4</sup>, François Kimmig<sup>2,3</sup>, Matthieu Caruel<sup>5</sup>, and Philippe Moireau<sup>2,3</sup>

<sup>1</sup>DMA, École normale supérieure, Université PSL, CNRS, Paris, France

<sup>2</sup>Inria, Team MΞDISIM, Inria Saclay, Palaiseau, France

<sup>3</sup>LMS, CNRS, École polytechnique, Institut Polytechnique de Paris, France

<sup>4</sup>CERMICS, École des ponts, 77420 Champs-sur-Marne, France

<sup>5</sup>Univ Paris Est Creteil, Univ Gustave Eiffel, CNRS, UMR 8208, MSME, 94010 Créteil

October 30, 2023

## Abstract

Muscle contraction at the macro level is a physiological process that is ultimately due to the interaction between myosin and actin proteins at the micro level. The actin-myosin interaction involves slow attachment and detachment responses and a rapid temporal change in protein conformation called *power-stroke*. Jump-diffusion models that combine jump processes between attachment and detachment with a mechanical description of the power-stroke have been proposed in the literature. However, the current formulations of these models are not fully compatible with the principles of thermodynamics.

To solve the problem of coupling continuous mechanisms with discrete chemical transitions, we rely on the mathematical formalism of Poisson random measures. First, we design an efficient stochastic formulation for existing muscle contraction PDE models. Then, we write a new jump-diffusion model for actin-myosin interaction. This new model describes both the behavior of muscle contraction on multiple time scales and its compatibility with thermodynamic principles. Finally, following a classical calibration procedure, we demonstrate the ability of the model to reproduce experimental data characterizing muscle behavior on fast and slow time scales.

## 1 Introduction

In muscle cells, force is generated by a highly ordered structure of contractile units called sarcomeres. The sarcomeres are arranged in series and form a dense network of fibrils. Each sarcomere is an array of interlocking myosin (thick) and actin (thin) filaments that slide past each other in opposite directions during contraction. The sliding is caused by the ATP-consuming cyclic interaction of the myosin heads, which protrude from the thick filament, with the actin filament. Under steady-state activation conditions, cardiac cells require about 100 ms to reach the peak force of a single twitch contraction at fixed sarcomere length (Janssen & Hunter, 1995). This duration reflects the characteristic timescale of the actin-myosin interaction cycle and, in particular, the typical rates (30-40 s<sup>-1</sup>) at which myosin heads attach to and detach from actin filaments. The process of attachment and detachment alone does not provide a complete explanation of how the actin-myosin interaction generates tension. It has indeed been shown that the force itself is generated by the power-stroke: a conformational change of the myosin executed while the myosin is attached to the actin. This fundamental event of the actin-myosin interaction can be interpreted as the swift relaxation towards thermal equilibrium of the internal degrees of freedom characterizing the conformation of the myosin protein (Huxley & Simmons, 1971). This rapid relaxation can be monitored in situ during the transient response of a muscle cell to rapid changes in loading
conditions (Ford et al., 1977; Caremani et al., 2016; Irving et al., 2000). The timescale associated with this response is a few milliseconds, supporting the idea that the power-stroke is a purely mechanical phenomenon operating independently of ATP hydrolysis (Marcucci & Truskinovsky, 2010; Caruel & Truskinovsky, 2018).

This multi-timescale nature of actin-myosin interaction raises the question of choosing an ap-
propriate modeling framework. The dynamics of protein-protein interactions and conformational
changes can be simulated using molecular dynamics techniques that provide access to the motion
of protein structure in a high-dimensional environment at short time scales (typically on the order
of picoseconds) (Tuckerman, 2023; Lelièvre et al., 2010). This type of framework is therefore not
adequate to model the actin-myosin attachment-detachment process, nor the power-stroke. Molec-
ular processes occurring at a millisecond timescale can be simulated with coarse grained molecular
dynamics, based on the definition of collective variables diffusing in a lower dimension energy land-
scape. These models often take the form of Langevin dynamics in non-convex multi-dimensional
potentials. In such potentials, transition between metastable states can be further simplified if
energy barriers are high enough for the so-called Kramers approximation to be valid (Schuss, 2010;
Van Kampen, 1992). In this case, inter-basin dynamics are adequately described by classical first-
order reaction kinetics, in which the remaining information about the underlying energy landscape
is embedded in detailed balance laws between forward and reverse reaction rates.

The latter formulation provides the framework for the vast majority of available muscle con-
traction models, in which the entire actin-myosin interaction cycle is represented by a sequence
of chemical-mechanical reactions (Månsson, 2010; Caremani et al., 2015; Månsson, 2016), each
reaction being equipped with transition rates that satisfy the detailed balance. This formulation
guarantees the compatibility of the model with the principles of thermodynamics, as shown in the
seminal work of T. L. Hill (Hill, 1974).

In the last decade, another approach was proposed to model the power- stroke (Marcucci &
Truskinovsky, 2010; Caruel & Truskinovsky, 2018), based on the definition of a onedimensional
(coarse) collective variable representing the conformation of the myosin motor. The dynamics of
this variable is a one-dimensional continuous stochastic process driven by the energy gradient flow
and thermal noise (Marcucci & Truskinovsky, 2010; Caruel & Truskinovsky, 2018). The associated
energy landscape has two wells representing the pre- and the post-power- stroke conformations.

More recently, this approach has been coupled with the more classical representation of the
(slower) attachment-detachment process, viewed as chemical reactions (Caruel et al., 2019; Kimmig
& Caruel, 2020). The resulting multi-timescale model of the actin-myosin interaction is param-
eterized by two continuous variables, and one discrete binary variable. The continuous variables
represent the position of the myosin head and its conformation, respectively. The discrete variable,
inherited from classical models, describes the binary state of attachment of the head (attached or
detached). Its dynamics is a jump process with transition rates parameterized by the continuous
variables and by the relative position of specific binding sites on the actin filament.

A salient feature of this class of jump-diffusion models is that when a binding event occurs, the position of the head instantaneously jumps to the position of the binding site and follows its deterministic motion until it unbinds. In this context, the binding event leads to the projection of a two-dimensional continuous dynamics in the detached state (position and conformation of the myosin head) onto a one-dimensional manifold in the attached state (conformation of the myosin head only, the position being deterministic). The reverse of this event is a jump from this one-dimensional space to the detached two-dimensional space. In particular, the position of the head immediately after a detachment event is probabilistic, which makes the formulation of the classical Hill-type detailed balance relation between forward and reverse rates nontrivial.

To circumvent this difficulty, the original paper (Caruel et al., 2019) makes the simplifying
assumption that the detachment occurs without changes in the internal state of the myosin head.
In particular, the position of the myosin head does not change; detachment occurs in place. This
assumption comes at the cost of not being able to demonstrate general agreement with thermody-
namic principles. In this paper, we propose a solution to this problem by proving that generalized
detailed balance relations can be used to ensure compatibility of our multi-timescale jump-diffusion model with the second principle of thermodynamics, thereby extending the approach initially proposed by T. L. Hill.

Our method is based on a reformulation of the model proposed by Caruel et al. (2019), using Poisson random measures, which provides an appropriate framework to describe the jump rates between multidimensional spaces.

We show that the resulting model is capable of quantitatively reproducing experimental benchmark observations that characterize both the timescale of the fast transient and the time scale of the entire actin-myosin interaction cycle. In addition, we also make predictions regarding the energy balances, and the thermodynamic efficiency associated with steady-state isotonic shortening.

That this new jump-diffusion model also provides an alternative to the fully mechanical description, including for the attachment-detachment process, proposed by Sheshka & Truskinovsky (2014). Our jump-diffusion model indeed shares with this latter model the desirable property of compatibility with thermodynamics but encompass the discrete description of the attachment-detachment transitions, which is widely used in the literature.

Recently, a new approach for modeling the actin-myosin interaction has been proposed (Chairtron et al., 2023). It relies on a new paradigm reformulating the Huxley-Hill framework (Huxley, 1957; Eisenberg & Hill, 1978) in an Eulerian setting. From a mathematical point of view, the Eulerian formulation provides relaxed conditions to establish the well-posedness of the equations. Moreover, for the same modeling complexity, the myosin tip can be located at an arbitrarily distant position from its anchor point. The model formulated in this framework, called the *h*-model, also ensures the compatibility with thermodynamics. It however differentiates itself from the present model on several points. On the one hand, the present model uses the classical Huxley-Hill framework, which uses the reference position of the myosin head as the model parametrization, while the *h*-model uses a Eulerian framework. On the other hand, the *h*-model does not include a description of the power-stroke and is thus unable to capture the fast timescale of muscle contraction.

The manuscript is organized as follows. We begin with a reformulation of the jump process between binding and unbinding using Poisson random measures, using the classical Huxley-57 model as a prototypical example of dynamics on the slow time scale, see Section 2. In the following section, we present the integration of the fast time scale by invoking the previously proposed mechanical model of power-stroke to formulate a complete jump-diffusion model that can be made compatible with thermodynamic principles using appropriate detailed balance relations, see Section 3. Emphasis is placed on the discretization of the model in Section 4, where we again begin with the prototypical example of Huxley-57 before turning to the full jump-diffusion model. Finally, in Section 5 we show how the simulation of the model compares with experimental benchmark results on fast and slow time scales.

# 2 Revisiting Huxley in a stochastic framework

## 2.1 The Huxley-Hill model

The seminal work of A.F. Huxley, begun in the 1950s and completed by T.L. Hill in the 1970s, led to the establishment of the so-called Huxley-Hill model (Huxley, 1957; Hill, 1974, 1976). This model describes the interaction between myosin filaments and actin filaments, both assumed to be rigid. The description focuses on the myosin head, and the population of myosin heads is assumed to be represented by a single individual interacting with a single actin filament.

This classical assumption is justified by the large number of myosin molecules within a half-sarcomere and by the probable absence of strong short-range elastic interaction between them along thick and thin filaments (Ford et al., 1981; Pertici et al., 2019). This implies that the motors effectively operate in parallel in physiological conditions. In such framework, we consider that the myosin heads are independent and that a large cluster is likely well represented by a representative element. This point is further discussed in Section 6. The binding actin sites are
<figure><img src="image_2.png" /><figcaption>Figure 1: (a) Huxley-Hill model, displayed in the attached state. The state of the myosin head is characterized by a single discrete stochastic variable $\alpha$ characterizing the state of attachment: $\alpha = 0$ for the detached state and $\alpha = 1$ for the attached state. The attachment and detachment rates are denoted by $f$ and $g$, respectively. The energy of the head is $w_0$ in the detached state and $w_1$ in the attached state. In the attached state, the head is located at the nearest binding side whose position is denoted by $s$, and its energy is usually defined as quadratic in $s$ (linear spring). (b) Proposed jump-diffusion model, displayed in the detached state. The state of the myosin head is characterized by three stochastic variables: $X$, $Y$ and $\alpha$, which represent the position of the head, the conformation and the state of attachment, respectively. The position of the nearest binding site is still denoted by $s$. The attachment-detachment process is characterized by four jumps whose rates are denoted $k_{0\to 1}$, $k_{0\to 1}^{\text{rev}}$, $k_{1\to 0}$ and $k_{1\to 0}^{\text{rev}}$. The energy landscape has two main ingredients: a quadratic term representing the elasticity of the head and a double-well potential associated with the conformational change.</figcaption></figure>

periodically located along the actin filament with a spatial period of length $d$. We define a possibly non-symmetric binding window [$s^-, s^+$] of length $d$ around the myosin head anchor point. Without loss of generality, in the theoretical representation of the model, we assume that the attachment window is equal to $[-d/2, d/2]$. We denote by $s$ the distance between the anchor point of the myosin head and the actin site located in the attachment window (see Figure 1(a)). In the context of the Huxley-Hill model, at any given time, the myosin head can interact with one actin site only: the actin site located in the attachment window. As a consequence, the myosin head cannot be attached to an actin site that leaves the attachment window. The actin site velocity is assumed to be given and is denoted by $\dot{x}_c(t)$ with the convention that $\dot{x}_c(t) < 0$ in the case of contraction.

The discrete variable $\alpha$ characterizes the binding state of the myosin head ($\alpha = 1$ when the myosin head is bound to the actin site and $\alpha = 0$ when it is detached). The dynamics of $\alpha$ are determined by the rate functions for attachment and detachment. Attachment (transition $\alpha = 0 \to \alpha = 1$) occurs at the rate $f(s)$ and detachment (transition $\alpha = 1 \to \alpha = 0$) occurs at the rate $g(s)$. The functions $f$ and $g$ are defined in the interval $(-d/2, +d/2)$ and are non-negative. The ratio of myosin heads whose actin site is at $s$ at time $t$ is denoted by $P_1(s, t)$. The dynamics of $P_1(s, t)$ is governed by

$$ \left\{ \begin{array}{ll} \partial_t P_1(t, s) + \dot{x}_c(t) \partial_s P_1(t, s) & t > 0, s \in (-d/2, d/2), & (1a) \\ P_1(t, -d/2) = P_1(t, d/2), & t \ge 0, & (1b) \\ P_1(0, s) = P_1^{\text{ini}}(s) \in [0, 1], & s \in [-d/2, d/2]. & (1c) \end{array} \right. $$

The ratio of myosin heads whose actin site is in $s$ at time $t$ and which are detached is simply given by $P_0(s, t) = 1 - P_1(s, t)$. We enforce periodic boundary conditions on the domain of definition of $P_1$ to ensure the well-posedness of the model (1) (Chaintron et al., 2023). However, describing the physiological system with the assumption that a myosin head can interact with a single actin requires not only that the values of the ratio of attached heads $P_1$ on both sides of the definition
interval are the same, but also that this value is zero. To obtain this property, the following condition is imposed on the attachment rates

$$ \lim_{s \to \pm d/2} f(s) = 0, \quad \lim_{s \to \pm d/2} g(s) = +\infty. \qquad (2) $$

## 2.2 Stochastic formalism using Poisson random measures

In this section, the Huxley-Hill model is described in terms of stochastic processes. This formulation will be useful to describe a more elaborated model in Section 3.2. We have chosen a formal approach here, but a more rigorous representation is possible by following the more mathematically sound presentation that has recently been proposed (Chaintron et al., 2023).

**The state variable as a stochastic process.** Our aim is to describe the myosin head at time $t$ by defining a random variable $\alpha_t \in \{0, 1\}$. The detached state corresponds to $\{\alpha_t = 0\}$, while $\{\alpha_t = 1\}$ stands for the attached state. We recall that $t \mapsto s(t)$ is the location in time of the closest actin site in the window $[-d/2, +d/2]$. This site slides with a given speed $\dot{x}_c(t)$, so that $\dot{s}(t) = \dot{x}_c(t)$. The following behavior is to be modeled:

* when $\alpha_t = 0$, $\alpha_t$ jumps to 1 within the infinitesimal time-interval $[t, t + \delta t]$ with probability $\delta t f(s(t))$ (attachment event).

* when $\alpha_t = 1$, $\alpha_t$ jumps to 0 within the infinitesimal time-interval $[t, t + \delta t]$ with probability $\delta t g(s(t))$ (detachment event).

We assume that the number of jumps within disjoint time intervals are independent, which means that when $\alpha_t = 0$, the waiting time $\Delta T_t$ for the next attachment event $0 \to 1$ follows an exponential law with rate $f(s(t))$. This jump rate is time inhomogeneous since it depends on the current time $t$ through $s(t)$. Equivalently, for any $t' \ge 0$,

$$ \mathbb{P}(\Delta T_t \le t') = 1 - \exp \left[ - \int_t^{t+t'} f(s(r)) dr \right]. \qquad (3) $$

The same is true for the detachment event $1 \to 0$. We now set up a stochastic equation for $\alpha_t$ that accounts for this behavior. This requires the mathematical notion of *Poisson random measure* which relies on an elaborated stochastic theory (Ikeda & Watanabe, 2014; Del Moral & Penev, 2017) that we will not detail here. Instead, we will rather focus on practical aspects. Formally, this theory allows us to write that

$$ \alpha_t = \int_0^t \mathbb{1}_{\alpha_{r-}=0} \int_{\mathbb{R}_+} \mathbb{1}_{z \le f(s(r))} N_{0 \to 1}(dr, dz) - \int_0^t \mathbb{1}_{\alpha_{r-}=1} \int_{\mathbb{R}_+} \mathbb{1}_{z \le g(s(r))} N_{1 \to 0}(dr, dz), \quad (4) $$

where $N_{0 \to 1}(dt, dz)$ and $N_{1 \to 0}(dt, dz)$ are two independent Poisson random measures on $\mathbb{R}_+ \times \mathbb{R}_+$ with intensity measure $dt \otimes dz$. The integral terms $\int_{\mathbb{R}_+} \mathbb{1}_{z \le f(s(r))} N_{0 \to 1}(dr, dz)$ and $\int_{\mathbb{R}_+} \mathbb{1}_{z \le g(s(r))} N_{1 \to 0}(dr, dz)$ can be seen as stochastic entities that select random jump times at the desired rate, as described below.

**Differential formulation using Poisson random measures.** The random measure $N_{0 \to 1}(dt, dz)$ can be defined from a sequence $(T_n^{\text{jump}}, Z_n)_{n \ge 1}$ of random variables. The $T_n^{\text{jump}}$ are candidate jump times, while $Z_n \in \mathbb{R}_+$ is a random variable uniformly chosen at each candidate jump time $T_n^{\text{jump}}$ (since $\mathbb{R}_+$ has infinite Lebesgue mass, the term “uniformly chosen” requires further explanation, for which we refer to (Del Moral & Penev, 2017)). The $T_n^{\text{jump}}$ are default jump times that are chosen with rate 1: the probability that a default jump occurs in the infinitesimal time-interval $[t, t + \delta t]$ is $\delta t$. If $Z_n \le f(s(T_n^{\text{jump}}))$, the candidate jump time $T_n^{\text{jump}}$ is accepted and a jump $0 \to 1$ occurs for $\alpha_t$. Otherwise, nothing happens for $\alpha_t$ until the next candidate jump
time. This way of selecting jump times for $\alpha_t$ guarantees that waiting times satisfy (3), since the probability that a jump $0 \to 1$ for $\alpha_t$ occurs within $[t, t + \delta t]$ is $\delta t \int_{\mathbb{R}_+} 1_{z \le f(s(t))} dz = \delta t f(s(t))$, as desired. The same is true for the detachment times $1 \to 0$ for $\alpha_t$. We can write down the Poisson random measure $N_{0 \to 1}(dt, dz)$ using Dirac masses as

$$N_{0 \to 1}(dt, dz) := \sum_{n \ge 1} \delta_{(T_n^{\text{jump}}, Z_n)}(dt, dz),$$

shedding lights on the integrals in (4). The same goes for the random measure $N_{1 \to 0}(dt, dz)$. This way of selecting jump times from a given random collection allows for a unified description of jump processes using equations like (4), the source of randomness being isolated in the default stochastic clocks $N_{0 \to 1}(dt, dz)$ and $N_{1 \to 0}(dt, dz)$. Further generalizations are possible, choosing $Z_n$ at random in a general set $E$ instead of $\mathbb{R}_+$. Moreover, the random variables $Z_n$ can be involved in the jump amplitudes, and not only in the selection of jump times. This allows us to define more general processes $X_t$ of the kind:

$$X_t := \int_0^t \int_E \Phi(r, z) N_{0 \to 1}(dr, dz) = \sum_{n \ge 1} 1_{T_n^{\text{jump}} \le t} \Phi(T_n^{\text{jump}}, Z_n), \quad (5)$$

for a function $\Phi : \mathbb{R}_+ \times E \to \mathbb{R}$. Equation (4) corresponds to the case $E = \mathbb{R}_+$ and $\Phi(r, z) = 1_{z \le f(s(r))}$. In Section 3.2, we will use this stochastic formalism with $E = (-d/2, +d/2) \times \mathbb{R}_+$. In particular, we will directly write $N_{0 \to 1}(dt, du, dv)$ rather than $N_{1 \to 0}(dt, dz)$ with $z = (u, v)$. To conclude this informal presentation, we point out that the stochastic theory rather writes

$$dX_t = \int_E \Phi(t, z) N_{0 \to 1}(dt, dz),$$

instead of (5), so that (4) reads in differential form:

$$d\alpha_t = 1_{\alpha_{t-}=0} \int_{\mathbb{R}_+} 1_{z \le f(s(t))} N_{0 \to 1}(dt, dz) - 1_{\alpha_{t-}=1} \int_{\mathbb{R}_+} 1_{z \le g(s(t))} N_{1 \to 0}(dt, dz). \quad (6)$$

This stochastic equation prescribes the (random) increment $d\alpha_t \in \{-1, 0, 1\}$ of $\alpha_t$ at time $t$, as a differential equation.

**Recovering the Huxley-Hill equations on the density.** From (6), we now show how to recover the Huxley-Hill system (1) on $P_\alpha(t, s)$; this will assess the consistency of our stochastic formulation with (1). For $t \ge 0$ and $\alpha$ in $\{0, 1\}$, $P_\alpha(t, s)$ will be seen as the probability that $\alpha_t = \alpha$ at $t$ while knowing that $s(t) = s$. This means that $\alpha \mapsto P_\alpha(t, s(t))$ is the density of the law of the random variable $\alpha_t$. As a consequence, for every test function $\varphi : \{0, 1\} \to \mathbb{R}$,

$$\mathbb{E}\varphi(\alpha_t) = \int_{\{0,1\}} \varphi(\alpha) P_\alpha(t, s(t)) d\alpha = \varphi(0) P_0(t, s(t)) + \varphi(1) P_1(t, s(t)).$$

Using that jump events $0 \to 1$ and $1 \to 0$ occur with respective rates $f(s(t))$ and $g(s(t))$, we get that

$$\mathbb{E}\varphi(\alpha_{t+\delta t}) - \mathbb{E}\varphi(\alpha_{t-}) = \delta t f(s(t)) \mathbb{E}1_{\alpha_{t-}=0}[\varphi(1) - \varphi(0)] + \delta t g(s(t)) \mathbb{E}1_{\alpha_{t-}=1}[\varphi(0) - \varphi(1)] + o(\delta t), \quad (7)$$

as $\delta t \to 0$. Let us define the operator $L_t$ which maps a function $\varphi$ to $L_t\varphi : \{0, 1\} \to \mathbb{R}$ such that

$$L_t\varphi(\alpha) := f(s(t))1_{\alpha=0}[\varphi(1) - \varphi(0)] - g(s(t))1_{\alpha=1}[\varphi(1) - \varphi(0)]. \quad (8)$$

In the stochastic literature this operator is called the *infinitesimal generator* (Del Moral & Penev, 2017) of the stochastic process $\alpha_t$. Equation (7) rewrites

$$\forall \varphi, \quad \frac{d}{dt} \int_{\{0,1\}} \varphi(\alpha) P_{\alpha}(t, s(t)) d\alpha = \int_{\{0,1\}} L_t \varphi(\alpha) P_{\alpha}(t, s(t)) d\alpha. \quad (9)$$
By duality, we now define the operator $L_t^*$ on $P_t: \alpha \mapsto P_\alpha(t, s(t))$ by

$$ \forall \varphi, \quad \int_{\{0,1\}} \varphi(\alpha)L_t^*P_t(\alpha)d\alpha = \int_{\{0,1\}} L_t\varphi(\alpha)P_t(\alpha)d\alpha. \qquad (10) $$

From (8)-(10), we get that

$$ L_t^*P_t(\alpha) = \mathbb{1}_{\alpha=0}[-f(s(t))P_t(0) + g(s(t))P_t(1)] + \mathbb{1}_{\alpha=1}[-g(s(t))P_t(1) + f(s(t))P_t(0)]. \quad (11) $$

Since (9) holds for every $\varphi$, we deduce from (10) that

$$ \partial_t P_t(\alpha) = L_t^* P_t(\alpha), $$

that is, using (11),

$$ \left\{ \begin{aligned} \frac{d}{dt} P_0(t, s(t)) &= -f(s(t))P_0(t, s(t)) + g(s(t))P_1(t, s(t)), \\ \frac{d}{dt} P_1(t, s(t)) &= -g(s(t))P_1(t, s(t)) + f(s(t))P_0(t, s(t)). \end{aligned} \right. $$

Since $\frac{d}{dt}P_\alpha(t, s(t)) = \partial_t P_\alpha(t, s(t)) + \dot{x}_c(t)\partial_s P_\alpha(t, s(t))$, this exactly recovers (1a). As usual, the boundary condition on $P_1(t, s(t))$ stems from the compulsory detachment condition $g(\pm d/2) = +\infty$ on the boundary. The above developments show the consistency of our stochastic formulation with the Huxley-Hill model (1).

## 2.3 Including thermodynamic constraints

The myosin motors convert chemical energy extracted from the hydrolysis of ATP molecules into mechanical work. Accurately representing these energetic aspects of the physiology is a challenge for biophysical models. The thermodynamic properties of the Huxley-Hill models, which are computed as averages over the whole population of heads, are well established in the literature (Hill, 1977; Eisenberg & Hill, 1978; Eisenberg et al., 1980; Kimmig et al., 2019). In this section, we show that these properties can be obtained from the stochastic reformulation of the model (6).

We denote by $w_\alpha(s)$, defined on $(-d/2, d/2)$, the energy of myosin head in the attached ($\alpha = 1$) and the detached ($\alpha = 0$) configurations, see Figure 1(a). The Huxley-Hill model does not describe the different configurations of myosin head in the detached state; the energy function $w_0$ is thus a constant (green line in Figure 1(a)). The energy function $w_1$ can be seen as the – possibly non-linear – elastic potential of the myosin head in the attached configuration (dashed red line in Figure 1(a)).

The concentration of ATP molecules being maintained constant in the cell, the input of energy by the hydrolysis of ATP molecules is represented by a constant chemical potential $\mu_T$.

As shown in the seminal works of T.L. Hill (Hill, 1977; Eisenberg & Hill, 1978; Eisenberg et al., 1980), the compatibility with the principles of thermodynamics can be obtained if the system allows, for any transition, that a reverse transition returning to the initial state can also take place. The hydrolysis of ATP molecules occurs during the detachment of the myosin head from the actin site, and the associated energy brought to the system modifies its energy landscape. Consequently, the apparent attachment – whose rate is given by $f$ – and the apparent detachment – whose rate is given by $g$ – cannot be considered as reverse from one another. Therefore, we define the following transitions:

* direct attachment: a jump $0 \to 1$ with rate $k_{0\to 1}(s)$;

* reverse attachment: a jump $1 \to 0$ with rate $k_{0\to 1}^{\text{rev}}(s)$;

* direct detachment: a jump $1 \to 0$ with rate $k_{1\to 0}(s)$;

* reverse detachment: a jump $0 \to 1$ with rate $k_{1\to 0}^{\text{rev}}(s)$.
The apparent transition rates f and g can be recovered as:

$$f(s) = k_{0 \to 1}(s) + k_{1 \to 0}^{\text{rev}}(s), \quad g(s) = k_{1 \to 0}(s) + k_{0 \to 1}^{\text{rev}}(s). \qquad (12)$$

The boundary condition (2) now becomes

$$\lim_{s \to \pm d/2} k_{0 \to 1}(s) = \lim_{s \to \pm d/2} k_{1 \to 0}^{\text{rev}}(s) = 0 \quad \text{and} \quad \lim_{s \to \pm d/2} k_{1 \to 0}(s) = \lim_{s \to \pm d/2} k_{0 \to 1}^{\text{rev}}(s) = +\infty.$$

**First principle.** Let $w_0$ denote the energy of the detached myosin head, and let $w_1(s)$ be the energy of the head when attached to an actin site located at $s$. The internal energy of the system is the average energy over the window:

$$\mathcal{U}(t) := \frac{1}{d} \int_{-d/2}^{+d/2} [w_0(s)P_0(t, s) + w_1(s)P_1(t, s)] ds.$$

The instantaneous power of external efforts is defined as

$$\dot{\mathcal{W}}(t) := \dot{x}_c(t)\tau_c(t),$$

the force generated by the attached motors being

$$\tau_c(t) := \frac{1}{d} \int_{-d/2}^{+d/2} \partial_s w_1(s) P_1(t, s) ds.$$

Each time a myosin head transitions from the attached state to the detached state, an ATP molecule is consumed and the system is provided with an energy supply $\mu_T$. The energy supply flux is thus given by

$$\dot{\mathcal{E}}(t) := \frac{1}{d} \int_{-d/2}^{+d/2} \mu_T [k_{1 \to 0}(s)P_1(t, s) - k_{1 \to 0}^{\text{rev}}(s)P_0(t, s)] ds,$$

where $[k_{1 \to 0}(s)P_1(t, s) - k_{1 \to 0}^{\text{rev}}(s)P_0(t, s)]$ is the net detachment flux. The heat dissipation eventu-
ally reads

$$\begin{aligned} \dot{\mathcal{Q}}(t) := & -\frac{1}{d} \int_{-d/2}^{+d/2} [w_0 - w_1(s)][k_{0 \to 1}(s)P_0(t, s) - k_{0 \to 1}^{\text{rev}}(s)P_1(t, s)]ds \\ & -\frac{1}{d} \int_{-d/2}^{+d/2} [w_1(s) - (w_0 - \mu_T)][k_{1 \to 0}(s)P_1(t, s) - k_{1 \to 0}^{\text{rev}}(s)P_0(t, s)]ds, \end{aligned} \quad (13)$$

so that the model satisfies the first principle of thermodynamics

$$\frac{d}{dt}\mathcal{U}(t) = \dot{\mathcal{W}}(t) + \dot{\mathcal{E}}(t) + \dot{\mathcal{Q}}(t), \qquad (14)$$

which appears as a direct consequence of (1a). Notice that the terms in the right hand side are positive when received by the system. In particular, mechanical power is produced when $\dot{\mathcal{W}} \le 0$, which is achieved when producing a positive force during shortening.

**Second principle.** The chemical potential of the state $\alpha$ at time $t$ is defined (Hill, 1977) as

$$\mu(t, \alpha, s) = w_{\alpha}(s) + k_{B}T \ln P_{\alpha}(t, s),$$

where $T$ is the absolute temperature and $k_B$ is the Boltzmann constant. The free energy of the system then reads

$$\mathcal{F}(t) := \frac{1}{d} \int_{-d/2}^{+d/2} [\mu(t, 0, s)P_0(t, s) + \mu(t, 1, s)P_1(t, s)] ds.$$
The system satisfies the second principle of thermodynamics if

$$ \frac{d}{dt} \mathcal{F}(t) = \dot{\mathcal{E}}(t) + \dot{\mathcal{W}}(t) - T \dot{\mathcal{S}}_{\text{prod}}(t), \qquad (15) $$

with $\dot{\mathcal{S}}_{\text{prod}}(t) \ge 0$. Using (15), it is possible to obtain $\dot{\mathcal{S}}_{\text{prod}}(t)$ from the above formulae for $\mathcal{F}(t)$, $\dot{\mathcal{E}}(t)$ and $\dot{\mathcal{W}}(t)$. This classical computation (Hill, 1977; Caruel et al., 2019; Chaintron et al., 2023) yields

$$ \dot{\mathcal{S}}_{\text{prod}}(t) = \frac{1}{Td} \int_{-d/2}^{+d/2} [\mu(t, 1, s) - \mu(t, 0, s)] [k_{0 \to 1}^{\text{rev}}(s) P_1(t, s) - k_{0 \to 1}(s) P_0(t, s)] ds \\ + \frac{1}{Td} \int_{-d/2}^{+d/2} [\mu(t, 1, s) - (\mu(t, 0, s) - \mu_T)] [k_{1 \to 0}^{\text{rev}}(s) P_1(t, s) - k_{1 \to 0}^{\text{rev}}(s) P_0(t, s)] ds. $$

As desired, the above quantity is always non-negative if forward and reverse rates satisfy the well-known detailed balance conditions (Hill, 1977)

$$ k_{0 \to 1}(s) = \exp \left[ - \frac{w_1(s) - w_0}{k_B T} \right] k_{0 \to 1}^{\text{rev}}(s), \qquad (16a) $$

$$ k_{1 \to 0}(s) = \exp \left[ - \frac{(w_0 - \mu_T) - w_1(s)}{k_B T} \right] k_{1 \to 0}^{\text{rev}}(s). \qquad (16b) $$

# 3 The power-stroke fast time-scale integration

## 3.1 Mechanical model of the power-stroke

The Huxley-Hill model presented in the previous section dates back to the original description of actin-myosin interaction by A.F. Huxley in 1957 (Huxley, 1957). This model was originally formulated to reproduce isotonic shortening experiments, in particular the force-velocity curve obtained experimentally by A.V. Hill (Hill, 1938; Huxley, 1957), and therefore does not describe the internal conformational change of the power-stroke. The model in its original form is thus unable to reproduce the fast transient response observed after rapid changes in loading conditions.

The Huxley-Hill model is usually supplemented by additional states that account for the conformational changes associated with the power-stroke. However, it has been shown that a simple model with only two attached states, pre- and post-power-stroke, cannot be calibrated with physiological parameters to reproduce the fast transients (Huxley & Tideswell, 1996). One way to circumvent this problem is to divide the power-stroke into several sub-steps, each with as many additional transition rates (Huxley & Simmons, 1971; Piazzesi & Lombardi, 1995; Caremani et al., 2015; Månsson, 2016; Smith et al., 2008).

An alternative approach was proposed in (Marcucci & Truskinovsky, 2010). It consists in representing the conformational change as a one-dimensional continuous stochastic process in a non-convex energy landscape that includes both the mechanics and the kinetics of the power-stroke.

Further elaboration have been proposed where this continuous representation of the power-stroke has been combined with the stochastic formulation of the classical discrete Huxley-Hill model presented in the previous section (Caruel et al., 2019; Kimmig & Caruel, 2020).

The model is parameterized by three stochastic variables, see Figure 1(b). The first variable, denoted $X_t$, is the position of the myosin head with respect to its anchor point on the thick filament. The second variable, denoted by $Y_t$, represents the conformation of the head (power-stroke). The third variable is $\alpha_t$ and defines the attachment state of the myosin head, as in the Huxley-Hill model. The energy of the head comprises an elastic term and a non-convex term whose form depends on the attachment state $\alpha$.
## 3.2 Complete jump-diffusion model

We now formulate our jump-diffusion model describing the dynamics of the variables $X_t$, $Y_t$ and $\alpha_t$. The position in time of the nearest actin site $t \mapsto s(t) \in [-d/2, +d/2]$ is a given function with speed $\dot{x}_c(t)$. The random variable $X_t$ is confined in the open interval $(-d/2, +d/2)$, and $X_t$ is discontinuous at attachment and detachment times. The internal variable $Y_t$ is continuous with values in $\mathbb{R}$. The full Markov process $(Z_t)_{t \ge 0}$ reads

$$Z_t = (\alpha_t, X_t, Y_t) \in \{0, 1\} \times (-d/2, +d/2) \times \mathbb{R}.$$

The internal variable $Y_t$ follows an overdamped Langevin dynamics with energy landscape $y \mapsto w_{\alpha_t}(X_t, y)$ and viscous damping coefficient $\eta_y$. When the head is detached ($\alpha_t = 0$), $X_t$ follows a continuous overdamped dynamics with energy landscape $x \mapsto w_0(x, Y_t)$ and damping coefficient $\eta_x$. When the head attaches, $X_t$ jumps to the location $s(t)$ of the nearest actin site, and then $X_t$ follows a deterministic dynamics with speed $\dot{x}_c(t)$ until the next detachment event, so that $X_t = s(t)$ while $\alpha_t = 1$. To ensure thermodynamic consistency, four jumps are involved: binding and unbinding reactions, together with their reverse reactions. The process $(Z_t)_{t \ge 0}$ is then defined through the following stochastic differential equation (SDE):

$$ \begin{cases} 
\begin{aligned}
d\alpha_t = {}& 1_{\alpha_{t-}=0} \left[ \int_{\mathbb{R}_+} 1_{z \le K_{0 \to 1}(X_{t-}, Y_t, s(t)) N_{0 \to 1}(dt, dz)} \right. \\
& \qquad \left. + \int_{\mathbb{R}_+} 1_{z \le K_{1 \to 0}^{\text{rev}}(X_{t-}, Y_t, s(t)) N_{1 \to 0}^{\text{rev}}(dt, dz)} \right] \\
& - 1_{\alpha_{t-}=1} \left[ \int_{-d/2}^{+d/2} \int_{\mathbb{R}_+} 1_{z \le K_{1 \to 0}(x, Y_t, s(t)) N_{1 \to 0}(dt, dx, dz)} \right. \\
& \qquad \left. + \int_{-d/2}^{+d/2} \int_{\mathbb{R}_+} 1_{z \le K_{0 \to 1}^{\text{rev}}(x, Y_t, s(t)) N_{1 \to 0}^{\text{rev}}(dt, dx, dz)} \right], \\
dX_t = {}& 1_{\alpha_{t-}=0} \left[ -\eta_x^{-1} \partial_x w_0(X_t, Y_t) dt + \sqrt{2\eta_x^{-1} k_B T} dB_t^x \right] + 1_{\alpha_{t-}=1} \dot{x}_c(t) dt \\
& + 1_{\alpha_{t-}=0} \left[ \int_{\mathbb{R}_+} (s(t) - X_{t-}) 1_{z \le K_{0 \to 1}(X_{t-}, Y_t, s(t)) N_{0 \to 1}(dt, dz)} \right. \\
& \qquad \left. + \int_{\mathbb{R}_+} (s(t) - X_{t-}) 1_{z \le K_{1 \to 0}^{\text{rev}}(X_{t-}, Y_t, s(t)) N_{1 \to 0}^{\text{rev}}(dt, dz)} \right] \\
& + 1_{\alpha_{t-}=1} \left[ \int_{-d/2}^{+d/2} \int_{\mathbb{R}_+} (x - s(t)) 1_{z \le K_{1 \to 0}(x, Y_t, s(t)) N_{1 \to 0}(dt, dx, dz)} \right. \\
& \qquad \left. + \int_{-d/2}^{+d/2} \int_{\mathbb{R}_+} (x - X_{t-}) 1_{z \le K_{0 \to 1}^{\text{rev}}(x, Y_t, s(t)) N_{0 \to 1}^{\text{rev}}(dt, dx, dz)} \right], \\
dY_t = {}& -\eta_y^{-1} \partial_y w_{\alpha_t}(X_t, Y_t) dt + \sqrt{2\eta_y^{-1} k_B T} dB_t^y.
\end{cases} \tag{17} $$

The values of diffusion coefficients for overdamped Langevin dynamics obey Einstein's law for diffusion. The stochastic processes $(B_t^x)_{x \ge 0}$ and $(B_t^y)_{y \ge 0}$ are independent Brownian motions which model thermal fluctuations (white noise), $T$ being the absolute temperature of the system. Four Poisson random measures are used, each of them modelling one chemical reaction corresponding to a jump for $(\alpha_t, X_t)$. The internal variable $Y_t$ is continuous, and is not affected by the jumps. As in (6), an extra variable $z$ is used to account for inhomogeneous jump rates. In the present situation, jump rates depend on the current configuration $(X_t, Y_t, s(t))$ of the entire system.

* $N_{0 \to 1}(dt, dz)$ is a Poisson random measure on $\mathbb{R}_+ \times \mathbb{R}_+$ with intensity measure $dt \otimes dz$. This measure models the direct attachment jump $0 \to 1$ for $\alpha_t$. At attachment events, the
myosin head binds at the nearest actin whose location is $s(t)$, so that the head position $X_t$ undergoes a jump of amplitude $s(t) - X_{t-}$. The rate of the jump $(0, x, y, s) \to (1, s, y, s)$ for $(\alpha_t, X_t, Y_t, s(t))$ is $K_{0\to1}(x, y, s)$.

* $N_{0\to1}^{\text{rev}}(dt, dx, dz)$ is a Poisson random measure on $\mathbb{R}_+ \times (-d/2, +d/2) \times \mathbb{R}_+$ with intensity measure $dt \otimes dx \otimes dz$. This measure models a jump $1 \to 0$ for $\alpha_t$ that corresponds to the reverse of the attachment reaction. In this reverse reaction, the myosin head detaches from the actin fiber and jumps to a new random location $x$. The jump amplitude is thus $x - s(t)$. This random location $x$ is chosen by the Poisson random measure, in such a way that the rate of the jump $(1, s, y, s) \to (0, x', y, s)$ for $(\alpha_t, X_t, Y_t, s(t))$ where $x'$ belongs to an infinitesimal cell of size $dx$ around $x$ is $K_{0\to1}^{\text{rev}}(x, y, s)dx$.

* $N_{1\to0}(dt, dx, dz)$ is a Poisson random measure on $\mathbb{R}_+ \times (-d/2, +d/2) \times \mathbb{R}_+$ with intensity measure $dt \otimes dx \otimes dz$. This measure models a jump $1 \to 0$ for $\alpha_t$ that corresponds to the direct detachment reaction. In this reaction, the myosin head detaches from the actin fiber and jumps to a new random location $x$. The jump amplitude is thus $x - s(t)$. This random location $x$ is chosen by the Poisson random measure, in such a way that the rate of the jump $(1, s, y, s) \to (0, x', y, s)$ for $(\alpha_t, X_t, Y_t, s(t))$ where $x'$ belongs to an infinitesimal cell of size $dx$ around $x$ is $K_{1\to0}(x, y, s)dx$.

* $N_{1\to0}^{\text{rev}}(dt, dz)$ is a Poisson random measure on $\mathbb{R}_+ \times \mathbb{R}_+$ with intensity measure $dt \otimes dz$. This measure models a jump $0 \to 1$ for $\alpha_t$ that corresponds to the reverse of the detachment reaction. In this reverse reaction, the myosin head attaches to the actin fiber at the location $s(t)$ of the nearest actin site, so that the jump amplitude is $s(t) - X_{t-}$. Similarly, the rate of this jump $(0, x, y, s) \to (1, s, y, s)$ for $(\alpha_t, X_t, Y_t, s(t))$ is $K_{1\to0}^{\text{rev}}(x, y, s)$.

**Detailed balance conditions.** Compared to previous works (Caruel et al., 2019), a major challenge here is to write balance conditions for continuous-to-discrete jumps and their reverse. A way to ensure that the model obeys the second principle of thermodynamics is to formulate detailed balance conditions between the jump rates defined above that are analogous to (16a)-(16b) used in Section 2.3.

The reverse attachment reaction involves a discrete-to-continuous jump and must select a new random location $x$ for the detached myosin head. The probability law that selects this new location should be compatible with balance conditions. To solve this issue, it is convenient to look at jumps at the level of the entire system $(\alpha_t, X_t, Y_t, s(t))$. The detailed balance conditions are then written on the total jumps $(0, x, y, s) \to (1, s, y, s)$ and $(1, s, y, s) \to (0, x, y, s)$:

$$K_{0\to1}(x, y, s) = h_x \exp \left[ - \frac{w_1(s, y) - w_0(x, y)}{k_B T} \right] K_{0\to1}^{\text{rev}}(x, y, s), \quad (18a)$$

$$K_{1\to0}(x, y, s) = h_x^{-1} \exp \left[ - \frac{(w_0(x, y) - \mu_T) - w_1(s, y)}{k_B T} \right] K_{1\to0}^{\text{rev}}(x, y, s), \quad (18b)$$

where $h_x$ is the characteristic length for $x$. The length $h_x$ is taken equal to the characteristic length of the power stroke stochastic dynamics. This characteristic length is needed to account for the discrete-to-continuous jumps and guarantees that (18a)-(18b) are homogeneous. It is sufficient to impose e.g. $h_x^{-1}K_{0\to1}(x, y, s)$ and $h_x^{-1}K_{1\to0}^{\text{rev}}(x, y, s)$ to prescribe all the rates.

**Global jump rates.** From the above rate $K_{1\to0}(x, y, s)$, it is possible to recover the global detachment rate - independently of the position resulting from the detachment jump - by integrating
(18b):

$$
\begin{align}
k_{1 \to 0}(y, s) &:= \int_{-d/2}^{+d/2} K_{1 \to 0}(x, y, s) \mathrm{d}x \nonumber \\
&= \int_{-d/2}^{+d/2} h_x^{-1} \exp \left[ - \frac{(w_0(x, y) - \mu_T) - w_1(s, y)}{k_B T} \right] K_{1 \to 0}^{\text{rev}}(x, y, s) \mathrm{d}x. \tag{19}
\end{align}
$$

The above rate is the analogous of $k_{1\to0}(s)$ that appears in Section 2.3. By integrating $K_{0\to1}^{\text{rev}}(x, y, s)$ over $x$, we similarly define the global rate $k_{0\to1}^{\text{rev}}(y, s)$.

**PDE formulation.** Let $p(t, s, \alpha, x, y)$ denote the probability density of $(\alpha_t, X_t, Y_t)$. When $\alpha = 1$, we define $p(t, s, 1, y) := p(t, s, 1, s, y)$. As in Section 2.2, we want to write a PDE system coupling $p(t, s, 0, x, y)$ and $p(t, s, 1, y)$. To do so let $\varphi : \{0, 1\} \times (-d/2, +d/2) \times \mathbb{R} \to \mathbb{R}$ denote any bounded $C^2$ test function whose derivatives are also bounded. By definition,

$$
\mathbb{E}[\varphi(\alpha_t, X_t, Y_t)] = \int_{-d/2}^{d/2} \int_{\mathbb{R}} [\varphi(0, x, y)p(t, s(t), 0, x, y)dx + \varphi(1, x, y)p(t, s(t), 1, y)\delta_{s(t)}(dx)] dy.
$$

As in Section 2.2, we look for a linear operator $L_t$ such that

$$
\forall \varphi, \quad \frac{d}{dt} \mathbb{E}[\varphi(\alpha_t, x_t, Y_t)] = \mathbb{E}[L_t \varphi(\alpha_t, x_t, Y_t)],
$$

in order to obtain that

$$
\partial_t p = L_t^* p, \tag{20}
$$

where $L_t^*$ is the dual operator of $L_t$. In the present situation, the infinitesimal generator $L_t$ can be decomposed as

$$
L_t = L_{t,\text{cont}} + L_{t,0\to1} + L_{t,0\to1}^{\text{rev}} + L_{t,1\to0} + L_{t,1\to0}^{\text{rev}},
$$

where the operators $L_{t,0\to1}$, $L_{t,0\to1}^{\text{rev}}$, $L_{t,1\to0}$ and $L_{t,1\to0}^{\text{rev}}$ correspond to the contributions of each jump. From the Ito formula, it is standard (Del Moral & Penev, 2017; Caruel et al., 2019) that the contribution related to the continuous part of the dynamics reads

$$
\begin{align*}
L_{t,\text{cont}}\varphi(\alpha, x, y) = {}& \mathbb{1}_{\alpha_t=0}[-\eta_x^{-1}\partial_x w_0(X_t, Y_t)\partial_x\varphi(0, x, y) + \eta_x^{-1}k_B T \partial_{xx}^2\varphi(0, x, y)] \\
& + \mathbb{1}_{\alpha_t=1}\dot{x}_c(t)\partial_x\varphi(1, x, y) - \eta_y^{-1}\partial_y w_\alpha(X_t, Y_t)\partial_y\varphi(\alpha, x, y) + \eta_y^{-1}k_B T \partial_{yy}^2\varphi(\alpha, x, y).
\end{align*}
$$

For each jump event, we then reason as in (7) to obtain the corresponding operator. The analogous
of (11) for the direct attachment jump is

$$
L_{t,0 \to 1} \varphi(\alpha, x, y) = \mathbb{1}_{\alpha=0} K_{0 \to 1}(x, y, s(t)) [\varphi(0, s(t), y) - \varphi(0, x, y)],
$$

and

$$
L_{t,1 \to 0} \varphi(\alpha, x, y) = \mathbb{1}_{\alpha=1} \int_{-d/2}^{+d/2} [\varphi(1, x', y) - \varphi(1, x, y)] K_{1 \to 0}(x', y, x) dx',
$$

for the direct detachment jump. We proceed similarly for $L_{t,0 \to 1}^{\text{rev}}$ and $L_{t,1 \to 0}^{\text{rev}}$. We then compute the dual operators $L_{t,\text{cont}}^*$, $L_{t,0 \to 1}^*$, $(L_{t,0 \to 1}^{\text{rev}})^*$, $L_{t,1 \to 0}^*$ and $(L_{t,1 \to 0}^{\text{rev}})^*$ for each component, and we sum them to obtain $L^*$. After straightforward computations, (20) eventually reads, for $(t, x, y)$ in
$$ \mathbb{R}_{+} \times (-d/2, +d/2) \times \mathbb{R}, $$

$$
\left\{
\begin{aligned}
\partial_t p(t, s, 0, x, y) &= -\dot{x}_c(t) \partial_s p(t, s, 0, x, y) \\
&\quad + \partial_x [\eta_x^{-1} \partial_x w_0(x, y) p(t, s, 0, x, y) + \eta_x^{-1} k_B T \partial_x p(t, s, 0, x, y)] \\
&\quad + \partial_y [\eta_y^{-1} \partial_y w_0(x, y) p(t, s, 0, x, y) + \eta_y^{-1} k_B T \partial_y p(t, s, 0, x, y)] \\
&\quad - [K_{0 \to 1}(x, y, s) + K_{1 \to 0}^{\text{rev}}(x, y, s)] p(t, s, 0, x, y) \\
&\quad + [K_{1 \to 0}(x, y, s) + K_{0 \to 1}^{\text{rev}}(x, y, s)] p(t, s, 1, y), \\
\partial_t p(t, s, 1, y) &= -\dot{x}_c(t) \partial_s p(t, s, 1, y) \\
&\quad + \partial_y [\eta_y^{-1} \partial_y w_1(s, y) p(t, s, 1, y) + \eta_y^{-1} k_B T \partial_y p(t, s, 1, y)] \\
&\quad - p(t, s, 1, y) \int_{-d/2}^{+d/2} [K_{1 \to 0}(x, y, s) + K_{0 \to 1}^{\text{rev}}(x, y, s)] dx \\
&\quad + \int_{-d/2}^{+d/2} [K_{0 \to 1}(x, y, s) + K_{1 \to 0}^{\text{rev}}(x, y, s)] p(t, s, 0, x, y) dx, \\
p(t, \pm d/2, 0, x, y) &= 0, \quad p(t, s, 0, \pm d/2, y) = 0, \quad p(t, \pm d/2, 1, y) = 0.
\end{aligned}
\right.
\tag{21}
$$

This provides a PDE formulation for our model that is equivalent to (17).

**Recovering the Huxley-Hill model back.** To obtain the Huxley-Hill PDE system (1) as a limit case of (21), let us assume that rates $K_{0 \to 1}$, $K_{0 \to 1}^{\text{rev}}$, $K_{1 \to 0}$ and $K_{1 \to 0}^{\text{rev}}$ do not depend on y anymore. We then define the averaged rates

$$ f(s) := \frac{1}{d} \int_{-d/2}^{+d/2} \int_{\mathbb{R}} [K_{0 \to 1}(x, y, s) + K_{1 \to 0}^{\text{rev}}(x, y, s)] dy dx, $$

$$ g(s) := \int_{-d/2}^{+d/2} \int_{\mathbb{R}} [K_{1 \to 0}(x, y, s) + K_{0 \to 1}^{\text{rev}}(x, y, s)] dy dx, $$

together with the integrated densities

$$ P_0(t, s) := \int_{-d/2}^{+d/2} \int_{\mathbb{R}} p(t, s, 0, x, y) dy dx \quad \text{and} \quad P_1(t, s) := \int_{\mathbb{R}} p(t, s, 1, y) dy. $$

Integrating (21) over x and y then exactly recovers the Huxley-Hill system (1) on $P_0$ and $P_1$. There exist other possible ways to recover (1) which are based on thermalisation assumptions on the internal variable y (Caruel et al., 2019; Kimmig & Caruel, 2020).

### 3.3 Thermodynamic balances

Following Section 2.3, it is now possible to compute the energy and entropy balances. The detail of the computations is provided in appendix.

**First principle.** The first principle of thermodynamics reads

$$ \frac{d}{dt} \mathcal{U}(t) = \dot{\mathcal{W}}(t) + \dot{\mathcal{E}}(t) + \dot{\mathcal{Q}}(t), \tag{22} $$

where the above quantities stand for:

* The internal energy:

$$ \mathcal{U}(t) := \int_{-d/2}^{+d/2} \int_{\mathbb{R}} [w_0(x, y) p(t, s(t), 0, x, y) dx + w_1(s(t), y) p(t, s(t), 1, y) \delta_{s(t)}(dx)] dy. $$
• The power of external efforts:

$$ \dot{W}(t) := \dot{x}_c(t)\tau_c(t) := \dot{x}_c(t) \int_{\mathbb{R}} \partial_s w_1(s(t), y)p(t, s(t), 1, y)dy. $$

• The flux term:

$$ \dot{\mathcal{E}}(t) = \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \mu_T[K_{1\to0}(x, y, s(t))p(t, s(t), 1, y) - K_{1\to0}^{\text{rev}}(x, y, s(t))p(t, s(t), 0, x, y)]dx dy. \quad (23) $$

• The heat dissipation:

$$ \begin{aligned} \dot{Q}(t) = & - \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \eta_x^{-1} (\partial_x w_0(x, y))^2 p(t, s(t), 0, x, y) dy dx \\ & - \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \eta_y^{-1} (\partial_y w_0(x, y))^2 p(t, s(t), 0, x, y) dy dx \\ & - \int_{\mathbb{R}} \eta_y^{-1} (\partial_y w_1(s(t), y))^2 p(t, s(t), 1, y) dy \\ & + \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \frac{k_B T}{\eta_x} \partial_{xx}^2 w_0(x, y) p(t, s(t), 0, x, y) dy dx \\ & + \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \frac{k_B T}{\eta_y} \partial_{yy}^2 w_0(x, y) p(t, s(t), 0, x, y) dy dx \\ & + \int_{\mathbb{R}} \frac{k_B T}{\eta_y} \partial_{yy}^2 w_1(s(t), y) p(t, s(t), 1, y) dy \\ & - \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \left\{ [w_0(x, y) - w_1(s(t), y)] \right. \\ & \qquad \left. \times [K_{0\to1}(x, y, s(t))p(t, s(t), 0, x, y) - K_{0\to1}^{\text{rev}}(x, y, s(t))p(t, s(t), 1, y)] \right\} dy dx \\ & - \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \left\{ [w_1(s(t), y) + \mu_T - w_0(x, y)] \right. \\ & \qquad \left. \times [K_{1\to0}(x, y, s(t))p(t, s(t), 1, y) - K_{1\to0}^{\text{rev}}(x, y, s(t))p(t, s(t), 0, x, y)] \right\} dy dx. \end{aligned} $$

The flux term accounts for the energy input due to ATP consumption. The heat dissipation involves the energy dissipated by the gradient system and thermal dissipation due to the noise, together with the heat productions which stem from jumps.

**Second principle.** Using (21), analogous developments give the free energy balance

$$ \frac{d}{dt} \mathcal{F}(t) = \dot{W}(t) + \dot{\mathcal{E}}(t) - T \dot{S}_{\text{prod}}(t) \quad (24) $$

The definition free energy here reads

$$ \mathcal{F}(t) := \int_{\mathbb{R}} \int_{-d/2}^{+d/2} [\mu(t, s(t), 0, x, y)p(t, s(t), 0, x, y)dx + \mu(t, s(t), 1, y)p(t, s(t), 1, y)\delta_{s(t)}(dx)]dy, $$

where the state chemical potentials at time t are

$$ \begin{aligned} \mu(t, s(t), 0, x, y) &:= w_0(x, y) + k_B T \ln[h_x h_y p(t, s(t), 0, x, y)], \\ \mu(t, s(t), 1, y) &:= w_1(s(t), y) + k_B T \ln[h_y p(t, s(t), 1, y)], \end{aligned} $$
where $h_x$ and $h_y$ respectively denote the characteristic lengths for $x$ and $y$. Both lengths are taken equal to the characteristic length of the power stroke stochastic dynamics. The scale $h_y$ plays no role in the following since it vanishes when computing differences. On the opposite, $h_x$ is needed to account for the discrete-to-continuous jumps. A direct computation gives the produced entropy

$$ 
\begin{aligned}
T \dot{S}_{\text{prod}}(t) = & \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \eta_x^{-1}(\partial_x \mu_0(x, y))^2 p(t, s(t), 0, x, y) dy dx \\
& + \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \eta_y^{-1}(\partial_y \mu_0(x, y))^2 p(t, s(t), 0, x, y) dy dx \\
& + \int_{\mathbb{R}} \eta_y^{-1}(\partial_y \mu_1(s(t), y))^2 p(t, s(t), 1, y) dy \\
& + \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \left\{ \begin{aligned}
& [\mu(t, s(t), 0, x, y) - \mu(t, s(t), 1, y)] \\
& \times [K_{0 \to 1}(x, y, s(t)) p(t, s(t), 0, x, y) - K_{0 \to 1}^{\text{rev}}(x, y, s(t)) p(t, s(t), 1, y)]
\end{aligned} \right\} dy dx \\
& + \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \left\{ \begin{aligned}
& [\mu(t, s(t), 1, y) + \mu_T - \mu(t, s(t), 0, x, y)] \\
& \times [K_{1 \to 0}(x, y, s(t)) p(t, s(t), 1, y) - K_{1 \to 0}^{\text{rev}}(x, y, s(t)) p(t, s(t), 0, x, y)]
\end{aligned} \right\} dy dx.
\end{aligned}
$$

Using the detailed balance conditions (18), this quantity is always non-negative. Therefore, the entropy production rate is always positive, as required by the second principle of thermodynamics.

The above developments show that the model can be made compatible with the first and the second principles of thermodynamics, if the condition of detailed balance (18) is satisfied.

To illustrate a consequence of this result, we consider a closed actin-myosin interaction cycle $C$, and integrate the free energy balance (24) over that cycle, we obtain

$$ \int_C \frac{d}{dt} \mathcal{F}(t) dt = \int_C \dot{\mathcal{W}}(t) dt + \int_C \dot{\mathcal{E}}(t) dt - T \int_C \dot{S}_{\text{prod}}(t) dt = 0. $$

In the absence of metabolic energy, $\mu_T = 0$ which implies $\dot{\mathcal{E}}(t) = 0$, this equality becomes

$$ \int_C \dot{\mathcal{W}}(t) dt = T \int_C \dot{S}_{\text{prod}}(t) dt \geq 0. $$

Hence, as expected, no net work can be produced in the absence of ATP. The result that the myosin heads produce work even though the rates locally verify the detailed balance condition comes from the fact that the energy level at the beginning and at the end of the interaction differ by $\mu_T$, hence breaking the overall balance of the cycle (Jülicher et al., 1997).

Finally, as the transition rates usually depend in a complex way on the state of the myosin head, the detailed balance constraint makes the calibration of the model more robust as transitions that are unfavorable from a thermodynamic point of view will be negligible “by design”, see Section 5.1.

# 4 Discretization aspects

In this section, we detail how the system (27) can be simulated numerically. The following presentation does not use a rigorous mathematical framework. Our objective is to allow the interested reader to implement our model and to give an additional point of view on the Poisson random measures. We first present how to deal with a single Poisson jump process and we then introduce a numerical scheme approximating the complete system coupling jump processes with over-damped Langevin dynamics.
## 4.1 Discretization of a single jump

We first present how to numerically simulate a stochastic process driven by a Poisson random measure. We consider a $\mathbb{R}$-valued stochastic process $(X_t)_{t \ge 0}$ whose dynamics is given by the equation

$$dX_t = \int_{\mathbb{R}} x \mathbf{1}_{z \le K(x,y(t))} N(dt, dz, dx), \qquad (25)$$

where $y$ is a continuous time-dependent parameter, and $N(dt, dz, dx)$ is a Poisson random measure on $\mathbb{R}_+ \times \mathbb{R}_+ \times \mathbb{R}$ with intensity measure $dt \otimes dz \otimes dx$. Each possible jump $X_t \to X_t + x$ has individual rate $K(x, y(t))$. The global jump rate for $X_t$ at time $t$ (according to the definition given in Section 3.2) is

$$k(y(t)) := \int_{\mathbb{R}} K(x, y(t)) dx.$$

We have seen in Section 2.2 that the waiting time $\Delta T_t$ for the next jump of $X_t$ is an exponential variable with (time-inhomogeneous) rate $k(y(t))$ (see Equation (3) for more details). Such a variable $\Delta T_t$ can be simulated by

$$\Delta T_t = \inf \left\{ t' > 0, \int_t^{t+t'} k(y(r)) dr \ge e_J \right\}, \qquad (26)$$

where $e_J \sim \mathcal{E}(1)$ is an exponential variable with parameter 1. At the jump time $t + \Delta T_t$, the jump amplitude is a random variable with law $k^{-1}(y(t + \Delta T_t))K(x, y(t + \Delta T_t))dx$.

We consider a numerical scheme with fixed time step $\delta t$ and we define the discretized times $t_n$ by $t_n := t_0 + n\delta t$, $\forall n \in \mathbb{N}$. We denote by $X_n$ the approximation of $X_{t_n}$. We assume that the time step is adapted to the rate of the considered random process, i.e. $\delta t \cdot \max_y [k(y)] \ll 1$.

To numerically approximate the waiting time $\Delta T_t$ between two jump events, we draw an exponential random variable $e_J$ with parameter 1 and we discretize the integral defining $\Delta T_t$ in (26) with the time step $\delta t$. The approximation of the integral is cumulatively updated along the simulated path $X_1, \dots, X_n$.

When the next jump time $t_n$ is reached, the increment $\Delta X_n$ for $X_n$ is drawn from the probability law $k^{-1}(y(t_n))K(x, y(t_n))dx$. Drawing a uniform variable $U \sim \mathcal{U}[0, 1]$, we simulate the increment by

$$\Delta X_n = \left[ x \mapsto \int_{-d/2}^{x} k^{-1}(y(t_n)) K(x', y(t_n)) dx' \right]^{-1} (U).$$

The integral can be computed using standard quadrature methods, and an approximated inversion can be performed using a dichotomy algorithm. The numerical scheme to simulate $X_n$ is summarized in Algorithm 1.

## 4.2 Complete model discretization

To simulate (17), we couple the numerical scheme for the jump process described in 4.1 with a classical Euler-Maruyama scheme for over-damped Langevin equations. The variables $\alpha$ and $X$ share the same jump clock. Therefore, only four jump events need to be considered. A counter $c_J$ and a threshold $e_J$ are associated with each of the four jumps. However, the four jumps cannot be simultaneously triggered at the same time. Indeed, in the detached state only attachment jumps are admissible, and conversely in the attached state. Therefore, the terms that update the counters are weighted by the activation factors $\alpha_n$ and $1 - \alpha_n$. Furthermore, in the case where two jumps are triggered in the same time interval $\delta t$, we select the jump for which the residual $c_J + k\delta t - e_J$ is the greatest. Once a jump occurs, the counters are reset to 0 and the thresholds are replaced by new unit exponential variables.
**Algorithm 1:** Numerical scheme of time step $\delta t$ to simulate an approximation $X_n$ of the variable $X_{t_n}$ which is driven by a Poisson random measure. The jump dynamics is defined by the rate $K$ with associated global rate $k$. The variable $c_J$ is a counter, and the unit exponential variable $e_J$ acts as a threshold. The final simulation time is given by $n_{\max}\delta t$.

Initial state: $n \leftarrow 0, c_J \leftarrow 0, e_J \sim \mathcal{E}(1)$

while $n < n_{\max}$ do

while $c_J < e_J$ do

$$c_J \leftarrow c_J + k(y(t_n))\delta t$$

$$X_{n+1} \leftarrow X_n$$

$$n \leftarrow n + 1$$

end

$$\Delta X_n \sim k^{-1}(y(t_n))K(x, y(t_n))dx$$

$$X_{n+1} \leftarrow X_n + \Delta X_n$$

$$c_J \leftarrow 0, e_J \sim \mathcal{E}(1)$$

$$n \leftarrow n + 1$$

end

The numerical scheme for the complete system (17) is given by

$$
\left\{
\begin{aligned}
X_{n+1} &= X_n + \delta t \left[ \alpha_n \dot{x}_c(t_n) + (1 - \alpha_n) \eta_x^{-1} \partial_x w_{\alpha_n}(X_n, Y_n) \right] \\
&\quad + (1 - \alpha_n) \sqrt{2\eta_x^{-1} k_B T \delta t} G_{n+1}^x \\
&\quad + (1 - \alpha_n) [s_n - X_n] \mathbb{1}_{c_{J0 \to 1} + (1 - \alpha_n) K_{0 \to 1} \delta t \ge e_{J0 \to 1}} \\
&\quad + (1 - \alpha_n) [s_n - X_n] \mathbb{1}_{c_{J1 \to 0}^{\text{rev}} + (1 - \alpha_n) K_{1 \to 0}^{\text{rev}} \delta t \ge e_{J1 \to 0}^{\text{rev}}}, \\
&\quad + \alpha_n [X_{n+1}^{+, \text{rev}} - X_n] \mathbb{1}_{c_{J0 \to 1}^{\text{rev}} + \alpha_n k_{0 \to 1}^{\text{rev}} \delta t \ge e_{J0 \to 1}^{\text{rev}}} \\
&\quad + \alpha_n [X_{n+1}^{-} - X_n] \mathbb{1}_{c_{J1 \to 0} + \alpha_n k_{1 \to 0} \delta t \ge e_{J1 \to 0}^{\text{rev}}}, \\
Y_{n+1} &= Y_n + \delta t \eta_y^{-1} \partial_y w_{\alpha_n}(X_n, Y_n) + \sqrt{2\eta_y^{-1} k_B T \delta t} G_{n+1}^y, \\
\alpha_{n+1} &= \alpha_n \\
&\quad + (1 - \alpha_n) [1 - 2\alpha_n] \mathbb{1}_{c_{J0 \to 1} + (1 - \alpha_n) K_{0 \to 1} \delta t \ge e_{J0 \to 1}} \\
&\quad + (1 - \alpha_n) [1 - 2\alpha_n] \mathbb{1}_{c_{J1 \to 0}^{\text{rev}} + (1 - \alpha_n) K_{1 \to 0}^{\text{rev}} \delta t \ge e_{J1 \to 0}^{\text{rev}}}, \\
&\quad + \alpha_n [1 - 2\alpha_n] \mathbb{1}_{c_{J0 \to 1}^{\text{rev}} + \alpha_n k_{0 \to 1}^{\text{rev}} \delta t \ge e_{J0 \to 1}^{\text{rev}}} \\
&\quad + \alpha_n [1 - 2\alpha_n] \mathbb{1}_{c_{J1 \to 0} + \alpha_n k_{1 \to 0} \delta t \ge e_{J1 \to 0}}, \\
s_{n+1} &= s_n + \delta t \dot{x}_c(t_n) \quad \text{in } \mathbb{T}_d^0,
\end{aligned}
\right.
\qquad (27)
$$

where

$$G_{n+1}^x \sim \mathcal{N}(0, 1),$$

$$G_{n+1}^y \sim \mathcal{N}(0, 1).$$

## 5 Numerical validations

In this section we show that the stochastic model is able to quantitatively reproduce the fast and slow timescale responses of a cardiac muscle fiber submitted to rapid change in loading conditions applied at the peak of a twitch contraction. More details about the experimental methodology can be found in Caremani et al. (2016).
<figure><img src="image_3.png" /><figcaption>Figure 2: Energy landscape of the myosin head. (a) Double quadratic functions representing the energy $u_α$, intrinsic to the power-stroke. (b) & (c) Energy of the myosin head in the detached and attached states, respectively. The horizontal dashed lines show the separation between the pre- and the post-power-stroke conformations. The oblique line in (c) of equation $y = -s$ shows the separation between the region where the tension generated by the cross-bridge is positive (right, traction) and the region where it is negative (left, compression). Detailed expressions and parameters values of the functions are available in Table 1.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Global jump rates of the four jumps. (a) & (b) Direct attachment and reverse detachment, respectively. (c) & (d) Reverse attachment and direct detachment, respectively. Notice that the rates are strictly positive, even in the white regions. The level sets of the detachment rates are intentionally cut-of at $4 \text{ µs}^{-1}$ for visibility, but the rate value are in fact above this threshold in the corresponding regions. Detailed expressions and parameters values of the functions are available in Table 1.</figcaption></figure>
## 5.1 Calibration principles

The calibration of the model follows the procedure detailed in Kimmig & Caruel (2020), of which we here recall the main principles. The detailed expressions of the different model ingredients can be found in Table 1.

Our model is based on the definition of a fixed binding window [$s^-, s^+$] of length $d$. We have chosen a non-symmetric interval with $s^- = -30$ nm and $s^+ = 10$ nm. This choice is made to allow the heads to remain attached over a long distance before reaching the limit of the interval during contraction. The distance $d = 40$ nm corresponds to the length of the actin regulatory units (Gordon et al., 2000).

The energy landscapes $w_\alpha$ are decomposed as follows

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

where $u$ is a double quadratic function (see Figure 2(a), showing the energy landscape along the line $x + y = 0$) and the second term accounts for the elasticity of the myosin head. The level set of the two landscapes $w_0(x, y)$ and $w_1(s, y)$ are shown in Figure 2(b) and (c), respectively.

The parameters of $w_1$ are adjusted so that the tension obtained from the thermal equilibrium distribution of $Y$ at a given $s$ matches the tension measured at the end of the fast transient phases (phase I and phase II), taking into account the distribution of the attached heads before the application of the fast load change, i.e., under isometric conditions.

The parameters of the detached potential $w_0$ are chosen such that the gradient flow brings the particles towards the pre-power-stroke configuration corresponding to low values of $y$, see Figure 2(a) and (b). Fine tuning can be made on this detached landscape parameters as there are no direct experimental constraint available.

The model is parameterized by four rates. According to the detailed balance conditions (18), only two of them need to be calibrated. We assume that the rate $K_{0\to 1}$ depends only on $y$ and $s$. We choose for $K_{0\to 1}$ the regularized piecewise constant function shown in Figure 3(a). Attachment is allowed with a nonzero value if the head is in the pre-power stroke (low values of $y$) and if the site is close to $s = 0$. The value of this increased rate and the corresponding $s$ interval are adjusted so that the fraction of attached myosin heads in the isometric state is equal to the experimentally measured value of 0.15 (Caremani et al., 2016). The reverse attachment rate $K_{0\to 1}^{\text{rev}}$ is then obtained from the detailed balance relation (18a).

The global rate $k_{1\to 0}(y, s)$ (see equation 19) is the direct detachment rate $1 \to 0$ knowing the value $(y, s)$ of $(Y_t, s(t))$. Writing that

$$k_{1\to 0}(y, s)J_{1\to 0}(y, s, dx) := K_{1\to 0}(x, y, s)dx,$$

defines the jump measure $J_{1\to 0}(y, s, dx)$. At detachment events of global rate $k_{1\to 0}(y, s)$, this jump measure is the probability law that selects the new value of $X_t$ after detachment. Other jump measures could be similarly defined for the three other jumps. In the case of attachment events, these measures degenerate into Dirac masses, since $x = s$ is imposed after attachment.

Following the classical approach (Eisenberg & Hill, 1978; Piazzesi & Lombardi, 1995), also used in a previous work (Kimmig & Caruel, 2020), we assume that the global detachment rate $k_{1\to 0}$ is a regularized piecewise constant function of $s$, shown in Figure 3(d). The jump measure is then defined as $J_{1\to 0}(y, s, dx) = (1/\lambda)\mathbb{1}_{|s-x|\le\lambda/2} dx$, which corresponds to a uniform distribution of $X_t$ after detachment on the interval $[s - \lambda/2, s + \lambda/2]$, where $\lambda$ is a parameter. Hence the position of the head after detachment remains close to its previous attached position. The detachment rate $k_{1\to 0}^{\text{rev}}$ is finally computed from the detailed balance relation (18b), and is represented in Figure 3(b). This rate is positive only around the line $s = x$, and in a $\sim 5$ nm interval around $s = -17$ nm. The location around the line $s = x$ is a direct consequence of the choice of the jump measure $J_{1\to 0}$, which defines the possible starting points for reverse detachment. The location around $s = -17$ nm is the consequence of the thermodynamic constraint imposed by the detailed balance conditions:
it is in this region where the energy levels $w_1(y, s)$ and $w_0(x, y, s) - \mu_T$ are close to each other, see Figure 2.

In the next sections, we compare the simulation of our model with experimental data from two benchmark experiments: the A.V. Hill force-velocity curve and the quick force recovery after a fast sarcomere length change. In the first case, we impose a constant negative sliding velocity ($\dot{s} = \dot{x}_c \le 0$), and in the second case, we maintain the system in isometric condition ($\dot{s} = 0$) until it reaches a steady state and then apply an instantaneous length change $\delta s$. In both cases, the average active tension $\tau_c$ is computed from $5 \times 10^5$ independent realizations of the stochastic process. The macroscopic stress is then given by $T_c = \rho_{\text{surf}}\tau_c$, where $\rho_{\text{surf}} = 1.25 \times 10^{17} \text{ m}^{-2}$ (Pinzauti et al., 2018) is the number of myosin heads in a longitudinal portion of a cell that has a thickness equal to the half-sarcomere length $\ell_{hs} \approx 1 \text{ µm per cross-sectional area, both quantities being defined in the reference configuration.}$

## 5.2 Slow time-scale validation

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

Figure 4: Slow timescale validation of the model. (a) Force-velocity relation showing the steady state force obtained while imposing the shortening velocity, averaged over $1 \times 10^5$ realizations. On this curve, three points, marked A, B and C have been selected to illustrate the internal steady state of the system in these shortening regimes. Experimental data are marked with open circles: circles (Caremani et al., 2016), squares (Van Heuningen et al., 1982), diamonds (de Tombe & ter Keurs, 1992), pentagons (de Tombe & ter Keurs, 1990) and triangles (Daniels et al., 1984). (b) Distribution of the system in the attached state. The level sets reproduce the energy landscape shown in figure 2(c), and the colored dots show three snapshots of the stochastic realizations corresponding to the points A, B and C shown in (a) for a subset of $5 \times 10^3$ randomly chosen realizations. (c) Distributions $P_1$ of attached heads at points A, B and C. (d) Distributions in the double well potential at point A, B and C, for the selected values of $s$ marked the star symbols beside the vertical dashed lines in (b).

The force-velocity relation produced by the calibrated model is in good agreement with available experimental data gathered from several sources Figure 4(a). The average tension $T_c$ is normalized by the stall force $T_0$ obtained in isometric condition (point A). In our simulation we obtained $T_0 = 116 \text{ kPa}$, which is in agreement with physiological measurements in situ (Caremani et al., 2016) (118 kPa). This value corresponds to a total fraction of attached heads of 0.146, which is also compatible with the value inferred from experimental data (Pinzauti et al., 2018) (0.15).
In Figure 4(b-d), we illustrate the microscopic configuration of the model in the steady state regimes corresponding to the points labeled A, B and C in panel (a). The simulations are performed using $1 \times 10^5$ stochastic realizations to produce the average results. Each colored dot in panel (b) represents one of the stochastic realizations for a subset of $5 \times 10^3$ randomly chosen realizations. The distributions of attached heads and the distribution of the variable $Y_t$ are represented in panels (c) and (d), respectively for each of the three sliding velocities. The distribution of $Y_t$ is shown for selected values of $s$: $s = 2$ nm for isometric condition (point A, blue) and $s = -3$ nm for intermediate and fast shortening velocities (point B, orange, and C, green).

In isometric condition (Point A, illustrated in blue), the distribution of attached heads reflects the form of the attachment rate, which is high only in the vicinity of the line $s = 0$. A few particles are also attached at $s \approx -12$ nm as a consequence of the reverse detachment jump, see Figure 3(b). The conformation of the attached heads is split between pre- and post-power-stroke, see panel(d) with a majority of heads in pre-power-stroke at $s = 2$ nm.

During steady state shortening (points B and C), the particles shift toward lower values of $s$ and higher values of $y$, following the minimum of the energy landscape (see panel (b)). As the shortening velocity increases, less particles are attached compared to the isometric condition, as shown by the lower density of dots in panel (b) and by the reduction in the integral of the distribution $P_1$ of attached heads in panel (c). For small to medium shortening velocities, the conformation distribution favors the post-power-stroke state (see the orange curve in panel (d)), but at larger velocities, for the same value of $s$ (here $-3$ nm), we observe that a significant fraction of attached heads remain in the pre-power-stroke conformation, see the green curve in (d). This phenomenon reflects the fact that for large velocities, the heads that are attached in per-power-stroke have less time to undergo their conformational change before being dragged by the sliding.

Because the model is compatible with thermodynamics, energetic quantities can be computed to better characterize the bioenergetics of muscle contraction and to further validate the capacity of the model to capture the physiology. We first compute the prediction of the power of external forces, whose generation is the main purpose of cardiac muscles. The model predictions show a good match with experimental data (see Figure 5(a)). Then, we turn our attention to the consumption of ATP, which constitutes the energy source of muscle cell for the production of mechanical work. The model predicts a steady increase of energy consumption with the shortening velocity. Here, only the data point in isometric conditions is available in the literature, and it is correctly predicted by the model. Finally, we compute the efficiency of the contraction, which is defined as $\eta = \dot{W}/\dot{E}$. We obtain a maximum efficiency of contraction around 40%, which is compatible with experimental estimates (Smith et al., 2005; Han et al., 2013).

The energy input rate shown in Figure 5(b) is by definition proportional to the net detachment flux, see (23). We show in Appendix6 that the detailed conditions imply that if ATP is withdrawn from the system by setting $\mu_T = 0$ without changing the attachment and reverse detachment rates – see Figure3 and Table1 – the net attachment flux becomes zero in steady state and so is the net detachment flux. Hence, the overall steady state cycling rate is also zero in the absence of ATP, as expected from thermodynamics.

## 5.3 Fast time-scale validation

The fast timescale response of the system after instant change of the sarcomere length is illustrated in Figure 6. We simulate the response of the system to a step applied on the variable $s$ with $1 \times 10^5$ stochastic realizations. Since the filaments elasticity is not taken into account explicitly in the model the results of our simulations must be post-processed to incorporate this effect and be compared to the experimental data, see Kimmig & Caruel (2020); Caruel et al. (2019) for more details. We use here a filament compliance of $17 \text{ nm MPa}^{-1}$ obtained from experimental data (Pinzauti et al., 2018) and the resulting step size is denoted by $\delta$hs and expressed in nanometer per half-sarcomere.

When the step is applied, the tension drops from the isometric value $T_0$ down to a minimum value to $T_1$ which depends linearly on the step size. The value of $T_1$ resulting from a step in $s$ of
<figure><img src="image_6.png" /></figure>

Figure 5: Slow timescale thermodynamic predictions of the model in the steady state at various imposed shortening velocities. The mechanical work, consumed ATP energy flux and the contraction efficiency are presented in panels (a), (b) and (c), respectively. They are computed according to their definitions given in Section 3.3 with averages computed over $1 \times 10^5$ thermodynamics realizations. Panel (a) presents the power of external efforts compared with experimental data. Experimental data are marked with open circles: circles (Caremani et al., 2016), squares (Van Heuningen et al., 1982), diamonds (de Tombe & ter Keurs, 1992), pentagons (de Tombe & ter Keurs, 1990) and triangles (Daniels et al., 1984). Panel (b) presents the consumed ATP energy flux. The associated experimental data in isometric conditions is computed from the measurements on skinned cardiac trabeculae (de Tombe & Stienen, 2007). The observed ATP consumption is 0.1085 ATP/s/kPa per myosin head. For an isometric stress of 118 kPa (Caremani et al., 2016) and an energy input of 100 zJ per ATP molecule (Barclay, 2015), this translates into an ATP energy consumption of 1.28 zJ ms⁻¹. Panel (c) presents the efficiency of contraction defines as the quotient of the power of external efforts by the consumed ATP energy flux.

$\delta_s = -4$ nm, which corresponds to $\delta h_s = -5.6$ nm per half-sarcomere when the filament elasticity is taken into account, is indicated by the labeled box in Figure 6(a).

Simultaneously with the step, the isometric internal distribution translates leftward, see the blue and violet dots and curves in panels (b-d)). In particular, the distribution of the conformation variable $Y_t$ does not change significantly.

After the quick tension recovery, the tension reaches a level $T_2$, see the labeled box in Figure 6(a). During this transition, the distribution $P_1$ of attached heads does not change, see (b and c), but the conformation relaxes towards the post-power-stroke energy well, see the upward displacement of the dots in (b) and the corresponding shift in the distribution in (d). This observation is fully compatible with the experimental observation that the fast transient response is essentially due to the internal relaxation of the power-stroke, without significant change of the number of attached myosin heads (Reconditi et al., 2004). It also validates the method used to calibrate the attached potential $w_1$.
<figure><img src="image_7.png" /><figcaption>Figure 6: Fast times-scale validation of the model against experimental data (Caremani et al., 2016). The simulations are performed using $1 \times 10^5$ stochastic realizations. (a) Tension obtained at the end of phase I ($T_1$, dark blue circles) and at the end of phase II ($T_2$, orange triangles) for different instantaneous half-sarcomere shortening ($\delta$hs) applied in isometric condition ($T_0$). The shortening values this take into account the filament compliance (see first paragraph of Section 5.3). The insert shows a typical evolution of the tension in response such shortening step, here $\delta$hs = 5.6 nm per half-sarcomere. (b) Scatter plot of the attached heads superimposed on the energy landscape $w_1$. Each dot represents one stochastic realization for a subset of $5 \times 10^3$ randomly chosen stochastic realizations. The three sets of dots represent the state of the system at $T_0$ (blue, isometric contraction), $T_1$ (dark blue, end of phase I) and $T_2$ (orange, end of phase II) for a step $\delta s = -4$ nm. The corresponding points on panel (a) are marked with the square boxes. (c) & (d) Distribution of attached heads and the distribution of $Y_t$ corresponding to the same points. For the distribution of $Y_t$ in (d), the value selected for $s$ are 2, for the isometric state $T_0$ and -2 for the two other stages $T_1$ and $T_2$. These values are indicated by the vertical dashed lines in (b).</figcaption></figure>

# 6 Conclusion and perspectives

In this paper, we propose a model of the actin-myosin interaction in the framework of muscle contraction. The model enhances a previously published model (Caruel et al., 2019) by improving it in a comprehensive formulation, which makes it consistent with thermodynamics. It describes the population of myosin heads inside the sarcomere by a representative head, which interacts with its nearest actin site. Our model couples the fast timescale internal mechanical behavior of the myosin to the slow timescale attachment-detachment dynamics of the myosin head. The internal mechanical state myosin head is characterized by two continuous variables, one being associated with the position of the myosin head tip and the other one with the power-stroke. The binding state (attached or detached) is accounted for by a discrete binary variable.

We first formulate the model in a stochastic framework. To do so, we rely on Poisson random measures. Using this mathematical object, the model equations can be unified in a jump-diffusion equations system. The associated Fokker-Planck equations are then established and the compatibility with the first and second principles of thermodynamics is demonstrated. This property was not guaranteed with the original model. The model is then calibrated using simple rules following strategies proposed in the literature (Caruel et al., 2019; Kimmig & Caruel, 2020). Using this calibration, we display its ability to reproduce experimental data characterizing the various timescales at which muscles operate. In addition, our formulation allows to predict energy balances over the whole range of the force-velocity curve, which could to be validated when the corresponding
experimental data become available.

Moreover, its thermodynamic properties make our model directly compatible with space multi-scale framework (Kimmig et al., 2019) paving the way for the coupling of the nano-scale power-stroke conformation change to the macroscopic organ deformation.

Several extensions of our newly proposed model may constitute valuable improvements.

In the literature, models sometimes rely on the hypothesis that the periodicity of the actin sites available for attachment along the actin filament is smaller than the periodicity of the actin filament helix (Piazzesi & Lombardi, 1995; de Tombe & Stienen, 2007; Pertici et al., 2018). In this case, the assumption that the representative myosin interacts with a single actin site is no longer valid and a multi-site framework must be considered. Our model could also be extended to a multi-site framework.

A natural perspective would also be to couple our current model with a model of actin activation by calcium, which is essential to simulate the contraction of cardiac cells. Such coupling has already been formalized in the framework of the Huxley'57 model (Kimmig et al., 2022) and could be adjusted for our jump-diffusion model.

In this paper, we chose to remain within the classical Huxley-Hill framework, where motor-motor interactions are neglected. As mentioned in our introduction, this assumption is justified by the large number of myosins within a sarcomere and by the relatively small linear compliance of the myofilaments, which is compatible with mechanical models where the motors operate in parallel, thus without nearest neighbor interactions (Linari et al., 1998; Powers et al., 2020). However, the parallel arrangement induces long-range (mean-field-type) interactions which generate cooperative effects, indeed. For instance a lump representation of the filament elasticity in series with a parallel cluster of myosin motors, allows the cross-bridges to synchronize their working stroke during short-timescale force recovery experiments, where the number of attached myosin heads is constant (Caruel & Truskinovsky, 2018; Caruel et al., 2013; Caruel & Truskinovsky, 2017). At longer timescales, these long-range interactions can generate dynamic instabilities and some interesting finite size effects that have been studied (Jülicher et al., 1997; Guérin et al., 2010, 2011; Vilfan & Duke, 2003; Erdmann & Schwarz, 2012; Wagoner & Dill, 2021) though without complete quantitative comparison with physiological data from muscle contraction. The consequences of motor-motor interactions are also studied in the context of cargo transport by processive molecular motors, where finite size effects are more prominent (Leighton & Sivak, 2022; Wijeratne et al., 2022; Geyer & Diez, 2023; Syed & Lee, 2019).

Furthermore, the model developed in this paper is used for computing the force in response to an imposed velocity. A natural perspective of our work is then to simulate the model in imposed force conditions or in imposed length conditions with explicit account of the filament elasticity, and study for instance whether dynamic instabilities are obtained with our physiological calibration.

These further developments provide the opportunity to take into account the long-range interactions between myosin heads working in parallel, and quantify their effect on the thermodynamic properties of the system in different loading conditions. Finally, using the tools of statistical physics, the average force generated by the population of myosin heads can then be related to the macroscopic deformation of the sarcomere. In this way, an additional meso-scale would be introduced into the spatial hierarchical modeling framework between the nano-scale and the macro-scale.

# Appendices

## Thermodynamic balances computation

We compute the thermodynamic balances along a trajectory $s(t)$. The sliding velocity is given by $\dot{x}_c = \dot{s}(t)$. Equivalent thermodynamic balances can also be obtained averaging over all possible values of the parameter $s$.
First principle

The internal energy is defined as

$$
\mathcal{U}(t) = \int_{\mathbb{R}} \int_{-d/2}^{+d/2} [w_0(x, y)p(t, s(t), 0, x, y)dx + w_1(s(t), y)p(t, s(t), 1, y)\delta_{s(t)}(dx)] dy.
$$

Differentiating with respect to time, we obtain

$$
\begin{align*}
\frac{d\mathcal{U}}{dt}(t) = & \int_{-d/2}^{+d/2} \int_{\mathbb{R}} w_0(x, y) [\partial_t p(t, s(t), 0, x, y) + \dot{x}_c \partial_s p(t, s(t), 0, x, y)] dy dx \\
& + \int_{\mathbb{R}} \dot{x}_c \partial_x w_1(s(t), y) p(t, s(t), 1, y) dy \\
& + \int_{\mathbb{R}} w_1(s(t), y) [\partial_t p(t, s(t), 1, y) + \dot{x}_c \partial_s p(t, s(t), 0, x, y)] dy.
\end{align*}
$$

Using (21), we obtain

$$
\begin{align*}
\frac{d\mathcal{U}}{dt}(t) = & \int_{-d/2}^{+d/2} \int_{\mathbb{R}} w_0(x, y) \Biggl[ \partial_x[\eta_x^{-1} \partial_x w_0(x, y) p(t, s(t), 0, x, y) + \eta_x^{-1} k_B T \partial_x p(t, s(t), 0, x, y)] \\
& + \partial_y[\eta_y^{-1} \partial_y w_0(x, y) p(t, s(t), 0, x, y) + \eta_y^{-1} k_B T \partial_y p(t, s(t), 0, x, y)] \\
& - [K_{0 \to 1}(x, y, s(t)) + K_{1 \to 0}^{\text{rev}}(x, y, s(t))] p(t, s(t), 0, x, y) \\
& + [K_{1 \to 0}(x, y, s(t)) + K_{0 \to 1}^{\text{rev}}(x, y, s(t))] p(t, s(t), 1, y) \Biggr] dy dx \\
& + \int_{\mathbb{R}} \dot{x}_c \partial_s w_1(s(t), y) p(t, s(t), 1, y) dy \\
& + \int_{\mathbb{R}} w_1(s(t), y) \Biggl[ \partial_y[\eta_y^{-1} \partial_y w_1(s(t), y) p(t, s(t), 1, y) + \eta_y^{-1} k_B T \partial_y p(t, s(t), 1, y)] \\
& - p(t, s(t), 1, y) \int_{-d/2}^{+d/2} [K_{1 \to 0}(x, y, s(t)) + K_{0 \to 1}^{\text{rev}}(x, y, s(t))] dx \\
& + \int_{-d/2}^{+d/2} [K_{0 \to 1}(x, y, s(t)) + K_{1 \to 0}^{\text{rev}}(x, y, s(t))] p(t, s(t), 0, x, y) dx \Biggr] dy.
\end{align*}
$$
Integrating by parts, we have

$$
\begin{aligned}
\frac{d\mathcal{U}}{dt}(t) = & \int_{\mathbb{R}} \dot{x}_c \partial_s w_1(s(t), y) p(t, s(t), 1, y) dy \\
& - \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \left[ \eta_x^{-1} (\partial_x w_0(x, y))^2 p(t, s(t), 0, x, y) + \eta_y^{-1} (\partial_y w_0(x, y))^2 p(t, s(t), 0, x, y) \right] dy dx \\
& - \int_{\mathbb{R}} \eta_y^{-1} (\partial_y w_1(s(t), y))^2 p(t, s(t), 1, y) dy \\
& + \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \left[ \frac{k_B T}{\eta_x} \partial_{xx}^2 w_0(x, y) p(t, s(t), 0, x, y) + \frac{k_B T}{\eta_y} \partial_{yy}^2 w_0(x, y) p(t, s(t), 0, x, y) \right] dy dx \\
& + \int_{\mathbb{R}} \frac{k_B T}{\eta_y} \partial_{yy}^2 w_1(s(t), y) p(t, s(t), 1, y) dy \\
& - \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \left\{ \begin{aligned}[t]
& [w_0(x, y) - w_1(s(t), y)] \\
& \times [K_{0 \to 1}(x, y, s(t)) p(t, s(t), 0, x, y) - K_{0 \to 1}^{\text{rev}}(x, y, s(t)) p(t, s(t), 1, y)]
\end{aligned} \right\} dy dx \\
& - \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \left\{ \begin{aligned}[t]
& [w_1(s(t), y) - w_0(x, y)] \\
& \times [K_{1 \to 0}(x, y, s(t)) p(t, s(t), 1, y) - K_{1 \to 0}^{\text{rev}}(x, y, s(t)) p(t, s(t), 0, x, y)]
\end{aligned} \right\} dy dx
\end{aligned}
$$

Introducing $\mu_T$ leads to

$$
\begin{aligned}
\frac{d\mathcal{U}}{dt}(t) = & \int_{\mathbb{R}} \dot{x}_c \partial_s w_1(s(t), y) p(t, s(t), 1, y) dy \\
& + \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \mu_T [K_{1 \to 0}(x, y, s(t)) p(t, s(t), 1, y) - K_{1 \to 0}^{\text{rev}}(x, y, s(t)) p(t, s(t), 0, x, y)] dy dx \\
& - \int_{-d/2}^{+d/2} \int_{\mathbb{R}} [\eta_x^{-1} (\partial_x w_0(x, y))^2 p(t, s(t), 0, x, y) + \eta_y^{-1} (\partial_y w_0(x, y))^2 p(t, s(t), 0, x, y)] dy dx \\
& - \int_{\mathbb{R}} \eta_y^{-1} (\partial_y w_1(s(t), y))^2 p(t, s(t), 1, y) dy \\
& + \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \left[ \frac{k_B T}{\eta_x} \partial_{xx}^2 w_0(x, y) p(t, s(t), 0, x, y) + \frac{k_B T}{\eta_y} \partial_{yy}^2 w_0(x, y) p(t, s(t), 0, x, y) \right] dy dx \\
& + \int_{\mathbb{R}} \frac{k_B T}{\eta_y} \partial_{yy}^2 w_1(s(t), y) p(t, s(t), 1, y) dy \\
& - \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \left\{ \begin{aligned}[t]
& [w_0(x, y) - w_1(s(t), y)] \\
& \times [K_{0 \to 1}(x, y, s(t)) p(t, s(t), 0, x, y) - K_{0 \to 1}^{\text{rev}}(x, y, s(t)) p(t, s(t), 1, y)]
\end{aligned} \right\} dy dx \\
& - \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \left\{ \begin{aligned}[t]
& [w_1(s(t), y) + \mu_T - w_0(x, y)] \\
& \times [K_{1 \to 0}(x, y, s(t)) p(t, s(t), 1, y) - K_{1 \to 0}^{\text{rev}}(x, y, s(t)) p(t, s(t), 0, x, y)]
\end{aligned} \right\} dy dx,
\end{aligned}
$$

and we can identify the first term with the power of external efforts $\dot{W}$, the second term with the ATP energy flux term $\dot{\mathcal{E}}$ and the last terms with the heat dissipation $\dot{Q}$ and we retrieve the first principle (22)

$$
\frac{d\mathcal{U}}{dt}(t) = \dot{W}(t) + \dot{\mathcal{E}}(t) + \dot{Q}(t).
$$
that are presented in Section 3.3.

## Second principle

The free energy is defined as

$$ \mathcal{F}(t) = \int_{\mathbb{R}} \int_{-d/2}^{+d/2} [\mu_0(t, s(t), x, y)p(t, s(t), 0, x, y)dx + \mu_1(t, s(t), y)p(t, s(t), 1, y)\delta_{s(t)}(dx)] dy. $$

Differentiating with respect to time, we obtain

$$ \begin{aligned} \frac{d\mathcal{F}}{dt}(t) = & \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \frac{d}{dt}\mu_0(t, s(t), x, y)p(t, s(t), 0, x, y)dydx \\ & + \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \mu_0(t, s(t), x, y)[\partial_t p(t, s(t), 0, x, y) + \dot{x}_c\partial_s p(t, s(t), 0, x, y)]dydx \\ & + \int_{\mathbb{R}} \frac{d}{dt}\mu_1(t, s(t), y)p(t, s(t), 1, y)dy \\ & + \int_{\mathbb{R}} \mu_1(t, s(t), y)[\partial_t p(t, s(t), 1, y) + \dot{x}_c\partial_s p(t, s(t), 1, y)]dy. \end{aligned} $$

Noticing that

$$ \left| \begin{aligned} \frac{d}{dt}\mu_0(t, s(t), x, y)p(t, s(t), 0, x, y) &= k_B T \frac{d}{dt}p(t, s(t), 0, x, y), \\ \frac{d}{dt}\mu_1(t, s(t), y)p(t, s(t), 1, y) &= \dot{x}_c\partial_s\omega_1(s(t), y)p(t, s(t), 1, y) + k_B T \frac{d}{dt}p(t, s(t), 1, y), \end{aligned} \right. $$

we obtain

$$ \begin{aligned} \frac{d\mathcal{F}}{dt}(t) = & \int_{-d/2}^{+d/2} \int_{\mathbb{R}} k_B T [\partial_t p(t, s(t), 0, x, y) + \dot{x}_c \partial_s p(t, s(t), 0, x, y)] dydx \\ & + \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \mu_0(t, s(t), x, y) [\partial_t p(t, s(t), 0, x, y) + \dot{x}_c \partial_s p(t, s(t), 0, x, y)] dydx \\ & + \int_{\mathbb{R}} \dot{x}_c \partial_s \omega_1(s(t), y) p(t, s(t), 1, y) dy \\ & + \int_{\mathbb{R}} k_B T [\partial_t p(t, s(t), 1, y) + \dot{x}_c \partial_s p(t, s(t), 1, y)] dy \\ & + \int_{\mathbb{R}} \mu_1(t, s(t), y) [\partial_t p(t, s(t), 1, y) + \dot{x}_c \partial_s p(t, s(t), 1, y)] dy. \end{aligned} $$
Using (21) the first and fourth terms cancel out. Using also (21) in the second and fifth terms, we have

$$
\begin{aligned}
\frac{d\mathcal{F}}{dt}(t) = & \int_{\mathbb{R}} \dot{x}_c \partial_s \omega_1(s(t), y) p(t, s(t), 1, y) dy \\
& + \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \mu_0(t, s(t), x, y) \left[ \begin{aligned}[t]
& \partial_x [\eta_x^{-1} \partial_x w_0(x, y) p(t, s(t), 0, x, y) + \eta_x^{-1} k_B T \partial_x p(t, s(t), 0, x, y)] \\
& + \partial_y [\eta_y^{-1} \partial_y w_0(x, y) p(t, s(t), 0, x, y) + \eta_y^{-1} k_B T \partial_y p(t, s(t), 0, x, y)] \\
& - [K_{0 \to 1}(x, y, s(t)) + K_{1 \to 0}^{\text{rev}}(x, y, s(t))] p(t, s(t), 0, x, y) \\
& + [K_{1 \to 0}(x, y, s(t)) + K_{0 \to 1}^{\text{rev}}(x, y, s(t))] p(t, s(t), 1, y)
\end{aligned} \right] dy dx \\
& + \int_{\mathbb{R}} \mu_1(t, s(t), y) \left[ \begin{aligned}[t]
& \eta_y^{-1} \partial_y w_1(s, y) p(t, s, 1, y) + \eta_y^{-1} k_B T \partial_y p(t, s, 1, y) \\
& - p(t, s, 1, y) \int_{-d/2}^{+d/2} [K_{1 \to 0}(x, y, s) + K_{0 \to 1}^{\text{rev}}(x, y, s)] dx \\
& + \int_{-d/2}^{+d/2} [K_{0 \to 1}(x, y, s) + K_{1 \to 0}^{\text{rev}}(x, y, s)] p(t, s, 0, x, y) dx
\end{aligned} \right] dy.
\end{aligned}
$$

Noticing that

$$
\left|
\begin{aligned}
\partial_x w_0(x, y) p(t, s(t), 0, x, y) + k_B T \partial_x p(t, s(t), 0, x, y) &= \partial_x \mu_0(t, s(t), x, y) p(t, s(t), 0, x, y), \\
\partial_y w_0(x, y) p(t, s(t), 0, x, y) + k_B T \partial_y p(t, s(t), 0, x, y) &= \partial_y \mu_0(t, s(t), x, y) p(t, s(t), 0, x, y), \\
\partial_y w_1(s(t), y) p(t, s(t), 1, y) + k_B T \partial_y p(t, s(t), 1, y) &= \partial_y \mu_1(t, s(t), y) p(t, s(t), 1, y),
\end{aligned}
\right.
$$

and integrating by parts, we have

$$
\begin{aligned}
\frac{d\mathcal{F}}{dt}(t) = & \int_{\mathbb{R}} \dot{x}_c \partial_s \omega_1(s(t), y) p(t, s(t), 1, y) dy \\
& - \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \eta_x^{-1} (\partial_x \mu_0(x, y))^2 p(t, s(t), 0, x, y) dy dx \\
& - \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \eta_y^{-1} (\partial_y \mu_0(x, y))^2 p(t, s(t), 0, x, y) dy dx \\
& - \int_{\mathbb{R}} \eta_y^{-1} (\partial_y \mu_1(s(t), y))^2 p(t, s(t), 1, y) dy \\
& - \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \left\{ \begin{aligned}[t]
& [\mu(t, s(t), 0, x, y) - \mu(t, s(t), 1, y)] \\
& \times [K_{0 \to 1}(x, y, s(t)) p(t, s(t), 0, x, y) - K_{0 \to 1}^{\text{rev}}(x, y, s(t)) p(t, s(t), 1, y)]
\end{aligned} \right\} dy dx \\
& - \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \left\{ \begin{aligned}[t]
& [\mu(t, s(t), 1, y) - \mu(t, s(t), 0, x, y)] \\
& \times [K_{1 \to 0}(x, y, s(t)) p(t, s(t), 1, y) - K_{1 \to 0}^{\text{rev}}(x, y, s(t)) p(t, s(t), 0, x, y)]
\end{aligned} \right\} dy dx.
\end{aligned}
$$
Introducing $\mu_T$, we have

$$
\begin{aligned}
\frac{d\mathcal{F}}{dt}(t) = & \int_{\mathbb{R}} \dot{x}_c \partial_s \omega_1(s(t), y) p(t, s(t), 1, y) dy \\
& + \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \mu_T [K_{1\to0}(x, y, s(t)) p(t, s(t), 1, y) - K_{1\to0}^{\text{rev}}(x, y, s(t)) p(t, s(t), 0, x, y)] dy dx \\
& - \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \eta_x^{-1} (\partial_x \mu_0(x, y))^2 p(t, s(t), 0, x, y) dy dx \\
& - \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \eta_y^{-1} (\partial_y \mu_0(x, y))^2 p(t, s(t), 0, x, y) dy dx \\
& - \int_{\mathbb{R}} \eta_y^{-1} (\partial_y \mu_1(s(t), y))^2 p(t, s(t), 1, y) dy \\
& - \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \left\{ \begin{aligned}[t]
& [\mu(t, s(t), 0, x, y) - \mu(t, s(t), 1, y)] \\
& \times [K_{0\to1}(x, y, s(t)) p(t, s(t), 0, x, y) - K_{0\to1}^{\text{rev}}(x, y, s(t)) p(t, s(t), 1, y)]
\end{aligned} \right\} dy dx \\
& - \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \left\{ \begin{aligned}[t]
& [\mu(t, s(t), 1, y) + \mu_T - \mu(t, s(t), 0, x, y)] \\
& \times [K_{1\to0}(x, y, s(t)) p(t, s(t), 1, y) - K_{1\to0}^{\text{rev}}(x, y, s(t)) p(t, s(t), 0, x, y)]
\end{aligned} \right\} dy dx.
\end{aligned}
$$

and we can identify the first term with the power of external efforts $\dot{W}$, the second term with the ATP energy flux term $\dot{\mathcal{E}}$ and the last terms with $T\dot{\mathcal{S}}_{\text{prod}}$ and we retrieve the second principle (24)

$$ \frac{d}{dt} \mathcal{F}(t) = \dot{W}(t) + \dot{\mathcal{E}}(t) - T \dot{\mathcal{S}}_{\text{prod}}(t). $$

The quantity $\dot{\mathcal{S}}_{\text{prod}}$ is defined such that

$$
\begin{aligned}
T \dot{\mathcal{S}}_{\text{prod}}(t) = & \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \eta_x^{-1} (\partial_x \mu_0(x, y))^2 p(t, s(t), 0, x, y) dy dx \\
& + \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \eta_y^{-1} (\partial_y \mu_0(x, y))^2 p(t, s(t), 0, x, y) dy dx \\
& + \int_{\mathbb{R}} \eta_y^{-1} (\partial_y \mu_1(s(t), y))^2 p(t, s(t), 1, y) dy \\
& + \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \left\{ \begin{aligned}[t]
& [\mu(t, s(t), 0, x, y) - \mu(t, s(t), 1, y)] \\
& \times [K_{0\to1}(x, y, s(t)) p(t, s(t), 0, x, y) - K_{0\to1}^{\text{rev}}(x, y, s(t)) p(t, s(t), 1, y)]
\end{aligned} \right\} dy dx \\
& + \int_{-d/2}^{+d/2} \int_{\mathbb{R}} \left\{ \begin{aligned}[t]
& [\mu(t, s(t), 1, y) + \mu_T - \mu(t, s(t), 0, x, y)] \\
& \times [K_{1\to0}(x, y, s(t)) p(t, s(t), 1, y) - K_{1\to0}^{\text{rev}}(x, y, s(t)) p(t, s(t), 0, x, y)]
\end{aligned} \right\} dy dx.
\end{aligned}
$$
<figure><img src="image_8.png" /><figcaption>Figure 7: Net fluxes and cycling rate obtained with $\mu_T = 100$ zJ (black) and $\mu_T = 0$ zJ (blue) during an isometric contraction ($\dot{x}_c = 0$). The simulations are performed with $1 \times 10^5$ stochastic realizations. We assume an initial distribution in which all myosin heads are detached, the $(x, y)$ variables are distributed according to the Boltzmann distribution associated with the detached energy potential $\omega_0(x, y)$, and the values of $s$ are uniformly distributed in $[s^-, s^+]$.</figcaption></figure>

The positivity of the first three terms is straightforward. For the last two terms, the positivity follows from the detailed balance conditions (18). Indeed, we have

$$ 
\begin{aligned}
& [\mu(t, s(t), 0, x, y) - \mu(t, s(t), 1, y)] \\
& \quad \times [K_{0 \to 1}(x, y, s(t))p(t, s(t), 0, x, y) - K_{0 \to 1}^{\text{rev}}(x, y, s(t))p(t, s(t), 1, y)] \\
& = [\omega_0(x, y) + k_B T \ln[h_x h_y p(t, s(t), 0, x, y)] - \omega_1(s(t), y) - k_B T \ln[h_y p(t, s(t), 1, y)]] \\
& \quad \times K_{0 \to 1}(x, y, s(t)) \left( p(t, s(t), 0, x, y) - h_x^{-1} \exp\left[\frac{w_1(s(t), y) - w_0(x, y)}{k_B T}\right] p(t, s(t), 1, y) \right), \\
& = k_B T K_{0 \to 1}(x, y, s(t)) \exp\left[\frac{-\omega_0(x, y)}{k_B T}\right] h_x^{-1} h_y^{-1} \\
& \quad \times \left[\ln\left(h_x h_y \exp\left[\frac{\omega_0(x, y)}{k_B T}\right] p(t, s(t), 0, x, y)\right) - \ln\left(h_y \exp\left[\frac{\omega_1(s(t), y)}{k_B T}\right] p(t, s(t), 1, y)\right)\right] \\
& \quad \times \left(h_x h_y \exp\left[\frac{\omega_0(x, y)}{k_B T}\right] p(t, s(t), 0, x, y) - h_y \exp\left[\frac{w_1(s(t), y)}{k_B T}\right] p(t, s(t), 1, y)\right).
\end{aligned}
$$

Since $[\log a - \log b][a - b] \ge 0$ for every $a, b > 0$, we have the inequality

$$ [\mu(t, s(t), 0, x, y) - \mu(t, s(t), 1, y)] \times [K_{0 \to 1}(x, y, s(t))p(t, s(t), 0, x, y) - K_{0 \to 1}^{\text{rev}}(x, y, s(t))p(t, s(t), 1, y)] \ge 0. $$

Similarly,

$$ [\mu(t, s(t), 1, y) + \mu_T - \mu(t, s(t), 0, x, y)] \times [K_{1 \to 0}(x, y, s(t))p(t, s(t), 1, y) - K_{1 \to 0}^{\text{rev}}(x, y, s(t))p(t, s(t), 0, x, y)] \ge 0, $$

and therefore $\dot{S}_{\text{prod}}(t) \ge 0$.

## Cycling rate and influence of the ATP chemical potential

We define the net attachment and detachment fluxes by

$$ k_{0 \to 1}^* = \int_{-d/2}^{+d/2} \int_{\mathbb{R}} [K_{0 \to 1}(x, y, s(t))p(t, s(t), 0, x, y) - K_{1 \to 0}^{\text{rev}}(x, y, s(t))p(t, s(t), 1, y)] dx dy $$

$$ k_{1 \to 0}^* = \int_{-d/2}^{+d/2} \int_{\mathbb{R}} [K_{1 \to 0}(x, y, s(t))p(t, s(t), 1, y) - K_{1 \to 0}^{\text{rev}}(x, y, s(t))p(t, s(t), 0, x, y)] dx dy, $$
respectively. The cycling rate *r* can then be computed by counting positively the net fluxes in the direction of the Lymn-Taylor cycle and negatively the net fluxes in opposite direction:

$$r = \min(k_{0 \to 1}^*, k_{1 \to 0}^*) \mathbb{1}_{k_{0 \to 1}^* > 0} \mathbb{1}_{k_{1 \to 0}^* > 0} - \min(-k_{0 \to 1}^*, -k_{1 \to 0}^*) \mathbb{1}_{k_{0 \to 1}^* < 0} \mathbb{1}_{k_{1 \to 0}^* < 0}.$$

The evolution of the attachment and detachment fluxes, and of the cycling rate is shown in Figure 7, starting from a configuration where all myosin heads are detached. As expected from thermodynamics, the cycling rate rapidly converges to zero in the absence of ATP ($\mu_T = 0$, blue lines) and reaches a positive value in the presence of ATP ($\mu_T = 100$ zJ, black lines).

## Calibration parameters

In Table 1, we gather all models parameters with associated selected values. The calibration principles are presented in Section 5.1. More details about the calibration method can be found in a previous paper (Kimmig & Caruel, 2020).

<table><caption>Table 1: Calibrated jump-diffusion model parameters.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td colspan="3">Characteristic dimensions</td></tr><tr><td>Drag coefficients</td><td>η<sub>x</sub>, η<sub>y</sub></td><td>0.0972 ms pN nm<sup>-1</sup></td></tr><tr><td>Temperature</td><td>T</td><td>298 K</td></tr><tr><td>Boltzmann constant</td><td>k<sub>B</sub></td><td>1.38 &times; 10<sup>-23</sup> J K<sup>-1</sup></td></tr><tr><td>Characteristic length in x-direction</td><td>h<sub>x</sub></td><td>11 nm</td></tr><tr><td>Characteristic length in y-direction</td><td>h<sub>y</sub></td><td>11 nm</td></tr><tr><td colspan="3">Power stroke potentials</td></tr><tr><td>Bistable element in attached state (α = 1)</td><td></td><td></td></tr><tr><td>u&#773;<sub>1</sub>(y) = {&kappa;<sub>1pre</sub>/2(y - y<sub>1pre</sub>)<sup>2</sup> + v<sub>1</sub> if y &lt; &ell;<sub>1</sub>,</td><td>&kappa;<sub>1pre</sub></td><td>5.60 pN nm<sup>-1</sup></td></tr><tr><td>&kappa;<sub>1post</sub>/2(y - y<sub>1post</sub>)<sup>2</sup> otherwise,</td><td>&kappa;<sub>1post</sub></td><td>1.33 pN nm<sup>-1</sup></td></tr><tr><td></td><td>&ell;<sub>1</sub></td><td>1.42 nm</td></tr><tr><td>v<sub>1</sub> = &kappa;<sub>1post</sub>/2(&ell;<sub>1</sub> - y<sub>1post</sub>)<sup>2</sup> - &kappa;<sub>1pre</sub>/2(&ell;<sub>1</sub> - y<sub>1pre</sub>)<sup>2</sup></td><td>y<sub>1post</sub></td><td>11 nm</td></tr><tr><td></td><td>y<sub>1pre</sub></td><td>0 nm</td></tr><tr><td>Bistable element in detached state (α = 0)</td><td></td><td></td></tr><tr><td>u&#773;<sub>0</sub>(y) = {&kappa;<sub>0pre</sub>/2(y - y<sub>0pre</sub>)<sup>2</sup> + v<sub>0</sub> + E if y &lt; &ell;<sub>0</sub>,</td><td>&kappa;<sub>0pre</sub></td><td>6.5 pN nm<sup>-1</sup></td></tr><tr><td>&kappa;<sub>0post</sub>/2(y - y<sub>0post</sub>)<sup>2</sup> + E otherwise,</td><td>&kappa;<sub>0post</sub></td><td>0.60 pN nm<sup>-1</sup></td></tr><tr><td></td><td>&ell;<sub>0</sub></td><td>1.42 nm</td></tr><tr><td>v<sub>0</sub> = &kappa;<sub>0post</sub>/2(&ell;<sub>0</sub> - y<sub>0post</sub>)<sup>2</sup> - &kappa;<sub>0pre</sub>/2(&ell;<sub>0</sub> - y<sub>1pre</sub>)<sup>2</sup></td><td>y<sub>0post</sub></td><td>6 nm</td></tr><tr><td></td><td>y<sub>0pre</sub></td><td>0</td></tr><tr><td></td><td>E</td><td>80 zJ</td></tr><tr><td colspan="3">Energy landscapes w<sub>α</sub>(y) = u<sub>α</sub>(y) + 1/2&kappa;(x + y)<sup>2</sup> (Figure 2)</td></tr><tr><td>u<sub>α</sub>(y) = u&#773;<sub>α</sub>(y + s&#773;<sub>α</sub>)</td><td>&kappa;</td><td>1.34 pN nm<sup>-1</sup></td></tr><tr><td></td><td>s&#773;<sub>0</sub></td><td>1.2 nm</td></tr><tr><td></td><td>s&#773;<sub>1</sub></td><td>1.2 nm</td></tr></tbody></table>
<table><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td colspan="3">Attachment rates (Figure 3)</td></tr><tr><td>K<sub>0→1</sub>(x, y, s) = k<sub>max</sub> [1 - tanh(α<sub>y</sub>(y - y<sub>0</sub>))]</td><td>k<sub>max</sub></td><td>1.21 ms<sup>-1</sup></td></tr><tr><td>α<sub>y</sub></td><td>α<sub>y</sub></td><td>8 nm<sup>-1</sup></td></tr><tr><td>· [1/2(1 + tanh [α<sub>s</sub>(s + s<sub>0</sub><sup>ℓ</sup>)])] 1<sub>(s&lt;0)</sub>(s)</td><td>α<sub>s</sub></td><td>8 nm<sup>-1</sup></td></tr><tr><td>+ 1/2(1 - tanh [α<sub>s</sub>(s - s<sub>0</sub><sup>r</sup>)]) 1<sub>(s≥0)</sub>(s)]</td><td>s<sub>0</sub><sup>ℓ</sup>→1<br>s<sub>0</sub><sup>r</sup>→1</td><td>3.82 nm<br>3.82 nm</td></tr><tr><td>K<sup>rev</sup><sub>1→0</sub>(x, y, s) = exp [w<sub>0</sub>(x, y) - μ<sup>*</sup><sub>T</sub> - w<sub>1</sub>(s, y)]</td><td>μ<sup>*</sup><sub>T</sub></td><td>100 zJ</td></tr><tr><td>· {k<sub>wℓ</sub> [1/2(1 + tanh [α<sub>x,wℓ</sub>(x - s + ℓ<sub>wℓ</sub>)] 1<sub>x-s&lt;0</sub>(x, s) + 1/2(1 - tanh [α<sub>x,wℓ</sub>(x - s - ℓ<sub>wℓ</sub>)] 1<sub>x-s&lt;0</sub>(x, s)]}</td><td>k<sub>wℓ</sub></td><td>4 ms<sup>-1</sup></td></tr><tr><td>α<sub>x,wℓ</sub></td><td>α<sub>x,wℓ</sub></td><td>5 nm<sup>-1</sup></td></tr><tr><td>ℓ<sub>wℓ</sub></td><td>ℓ<sub>wℓ</sub></td><td>10 nm</td></tr><tr><td>α<sub>s,wℓ</sub></td><td>α<sub>s,wℓ</sub></td><td>8 nm<sup>-1</sup></td></tr><tr><td>s<sub>wℓ</sub></td><td>s<sub>wℓ</sub></td><td>-19 nm</td></tr><tr><td>k<sub>wr</sub></td><td>k<sub>wr</sub></td><td>5 ms<sup>-1</sup></td></tr><tr><td>α<sub>x,wr</sub></td><td>α<sub>x,wr</sub></td><td>5 nm<sup>-1</sup></td></tr><tr><td>ℓ<sub>wr</sub></td><td>ℓ<sub>wr</sub></td><td>10 nm</td></tr><tr><td>α<sub>s,wr</sub></td><td>α<sub>s,wr</sub></td><td>4 nm<sup>-1</sup></td></tr><tr><td>s<sub>wr</sub></td><td>s<sub>wr</sub></td><td>9 nm</td></tr><tr><td>k<sub>p</sub></td><td>k<sub>p</sub></td><td>1 ms<sup>-1</sup></td></tr><tr><td>α<sub>x,p</sub></td><td>α<sub>x,p</sub></td><td>5 nm<sup>-1</sup></td></tr><tr><td>ℓ<sub>p</sub></td><td>ℓ<sub>p</sub></td><td>1 nm</td></tr><tr><td>α<sub>s,p</sub></td><td>α<sub>s,p</sub></td><td>10 nm<sup>-1</sup></td></tr><tr><td>s<sub>p</sub></td><td>s<sub>p</sub></td><td>-9 nm</td></tr><tr><td>k<sub>b</sub></td><td>k<sub>b</sub></td><td>1.89 ms<sup>-1</sup></td></tr><tr><td>α<sub>b</sub></td><td>α<sub>b</sub></td><td>5 nm<sup>-1</sup></td></tr><tr><td>ℓ<sub>b</sub></td><td>ℓ<sub>b</sub></td><td>0.2 nm</td></tr><tr><td>k<sub>b</sub> [1/2(1 + tanh [α<sub>b</sub>(x - s + ℓ<sub>b</sub>)] 1<sub>x-s&lt;0</sub>(x, s) + 1/2(1 - tanh [α<sub>b</sub>(x - s - ℓ<sub>b</sub>)] 1<sub>x-s&lt;0</sub>(x, s)])</td><td>x<sub>ℓ</sub><br>x<sub>h</sub></td><td>-12 nm<br>-8 nm</td></tr><tr><td>· (1 - 1<sub>x<sub>h</sub>&gt;x&gt;x<sub>ℓ</sub></sub>(x) · 1<sub>y-y<sub>ℓ</sub>&gt;0</sub>(y) · 1<sub>s<sub>h</sub>&gt;s&gt;s<sub>ℓ</sub></sub>(s)</td><td>y<sub>ℓ</sub><br>s<sub>ℓ</sub><br>s<sub>h</sub></td><td>12.5 nm<br>-24 nm<br>-16 nm</td></tr><tr><td colspan="3">Geometric parameter</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>

## References

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

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.
URL https://link.aps.org/doi/10.1103/PhysRevLett.110.248103

Caruel, M., Moireau, P., & Chapelle, D. (2019). Stochastic modeling of chemical-mechanical coupling in striated muscles. *Biomechanics and Modeling in Mechanobiology*, 18(3), 563–587.

Caruel, M., & Truskinovsky, L. (2017). Bi-stability resistant to fluctuations. *Journal of the Mechanics and Physics of Solids*, 109, 117–141.
URL https://linkinghub.elsevier.com/retrieve/pii/S0022509616309103

Caruel, M., & Truskinovsky, L. (2018). Physics of muscle contraction. *Reports on Progress in Physics*, 81(3), 036602.
Chaintron, Louis-Pierre, Caruel, Matthieu, & Kimmig, François (2023). Modeling acto-myosin interaction: beyond the huxley-hill framework. *MathematicS in Action, 12*, 191-226.

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., & Stienen, G.J.M. (2007). Impact of temperature on cross-bridge cycling kinetics in rat myocardium. *The Journal of Physiology, 584*(2), 591-600.
URL http://doi.wiley.com/10.1113/jphysiol.2007.138693

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.

Del Moral, P., & Penev, S. (2017). *Stochastic Processes: From Applications to Theory*. Chapman and Hall/CRC.

Eisenberg, E., & Hill, T. L. (1978). A cross-bridge model of muscle contraction. *Progress in Biophysics and Molecular Biology, 33*(1), 55-82.

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

Erdmann, Thorsten, & Schwarz, Ulrich S. (2012). Stochastic Force Generation by Small Ensembles of Myosin II Motors. *Physical Review Letters, 108*(18), 188101.

Ford, L. E., Huxley, A. F., & Simmons, R. M. (1977). Tension responses to sudden length change in stimulated frog muscle fibres near slack length. *The Journal of Physiology, 269*(2), 441-515.

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.
URL https://onlinelibrary.wiley.com/doi/10.1113/jphysiol.1981.sp013582

Geyer, Veikko F., & Diez, Stefan (2023). Horizontal Magnetic Tweezers to Directly Measure the Force-Velocity Relationship for Multiple Kinesin Motors. *Small, 19*(30), 2300558.

Gordon, A. M., Homsher, E., & Regnier, M. (2000). Regulation of contraction in striated muscle. *Physiological reviews, 80*(2), 853-924.
URL https://www.physiology.org/doi/10.1152/physrev.2000.80.2.853

Guérin, T., Prost, J., & Joanny, J.-F. (2010). Dynamic Instabilities in Assemblies of Molecular Motors with Finite Stiffness. *Physical Review Letters, 104*(24), 248102.

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.

Han, June-Chiew, Taberner, Andrew J., Nielsen, Poul M. F., & Loiselle, Denis S. (2013). Interventricular comparison of the energetics of contraction of trabeculae carneae isolated from the rat heart. *The Journal of Physiology, 591*(3), 701-717.

Hill, A. V. (1938). The heat of shortening and the dynamic constants of muscle. *Proceedings of the Royal Society of London. Series B - Biological Sciences, 126*(843), 136-195. Publisher: The Royal Society.

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.

Huxley, A. F. (1957). Muscle structure and theories of contraction. *Progress in Biophysics and Biophysical Chemistry*, **7**, 255–318.

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

Huxley, A. F., & Tideswell, S. (1996). Filament compliance and tension transients in muscle. *Journal of Muscle Research and Cell Motility*, **17**(4), 507–511.

Ikeda, N., & Watanabe, S. (2014). *Stochastic differential equations and diffusion processes*. Elsevier.

Irving, M., Piazzesi, G., Lucii, L., Sun, Y.-B., Harford, J.J., Dobbie, I.M., Ferenczi, M.A., Reconditi, M., & Lombardi, V. (2000). Conformation of the myosin motor during force generation in skeletal muscle. *Nature Structural Biology*, **7**(6), 482–485.

Janssen, P.M., & Hunter, W.C. (1995). Force, not sarcomere length, correlates with prolongation of isosarcometric contraction. *AJP: Heart and Circulatory Physiology*, **269**(2), H676–H685.

Jülicher, F., Ajdari, A., & Prost, J. (1997). Modeling molecular motors. *Reviews of Modern Physics*, **69**(4), 1269–1282.

URL https://link.aps.org/doi/10.1103/RevModPhys.69.1269

Kimmig, F., & Caruel, M. (2020). Hierarchical modeling of force generation in cardiac muscle. *Biomechanics and Modeling in Mechanobiology*, **19**(6), 2567–2601.

Kimmig, François, Caruel, Matthieu, & Chapelle, Dominique (2022). Varying thin filament activation in the framework of the huxley'57 model. *International Journal for Numerical Methods in Biomedical Engineering*, **38**(12), e3655.

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.

URL https://amses-journal.springeropen.com/articles/10.1186/s40323-019-0128-9

Leighton, Matthew P., & Sivak, David A. (2022). Dynamic and Thermodynamic Bounds for Collective Motor-Driven Transport. *Physical Review Letters*, **129**(11), 118102.

Lelièvre, Tony, Stoltz, Gabriel, & Rousset, Mathias (2010). *Free Energy Computations: A Mathematical Perspective*. Imperial College Press.

Linari, Marco, Dobbie, Ian, Reconditi, Massimo, Koubassova, Natalia, Irving, Malcolm, Piazzesi, Gabriella, & Lombardi, Vincenzo (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.

URL https://linkinghub.elsevier.com/retrieve/pii/S0006349598779548

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. (2016). Actomyosin based contraction: one mechanokinetic model from single molecules to muscle? *Journal of Muscle Research and Cell Motility*, **37**(6), 181–194.

Marcucci, L., & Truskinovsky, L. (2010). Muscle contraction: A mechanical perspective. *The European Physical Journal E*, **32**(4), 411–418.
Pertici, I., Bongini, L., Melli, L., Bianchi, G., Salvi, L., & Falorsi, G. (2018). A myosin ii nanomachine mimicking the striated muscle. *Nature communications, 9*(1), 3532.
URL http://www.nature.com/articles/s41467-018-06073-9

Pertici, Irene, Caremani, Marco, & Reconditi, Massimo (2019). A mechanical model of the half-sarcomere which includes the contribution of titin. *Journal of Muscle Research and Cell Motility, 40*(1), 29-41.
URL http://link.springer.com/10.1007/s10974-019-09508-y

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

Pinzauti, F., Pertici, I., Reconditi, M., Narayanan, T., Stienen, G. 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.
URL http://doi.wiley.com/10.1113/JP275579

Powers, Joseph D., Bianco, Pasquale, Pertici, Irene, Reconditi, Massimo, Lombardi, Vincenzo, & Piazzesi, Gabriella (2020). Contracting striated muscle has a dynamic I-band spring with an undamped stiffness 100 times larger than the passive stiffness. *The Journal of Physiology, 598*(2), 331-345.
URL https://onlinelibrary.wiley.com/doi/10.1113/JP278713

Reconditi, Massimo, Linari, Marco, Lucii, Leonardo, Stewart, Alex, Sun, Yin-Biao, Boesecke, Peter, Narayanan, Theyencheri, Fischetti, Robert F., Irving, Tom, Piazzesi, Gabriella, Irving, Malcolm, & Lombardi, Vincenzo (2004). The myosin motor in muscle generates a smaller and slower working stroke at higher load. *Nature, 428*(6982), 578-581.

Schuss, Zeev (2010). *Theory and Applications of Stochastic Processes*, vol. 170 of *Applied Mathematical Sciences*. Springer New York.
URL http://link.springer.com/10.1007/978-1-4419-1605-1

Sheshka, R., & Truskinovsky, L. (2014). Power-stroke-driven actomyosin contractility. *Physical review. E, Statistical, nonlinear, and soft matter physics, 89*(1), 012708.

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, Nicholas P., Barclay, Christopher J., & Loiselle, Denis S. (2005). The efficiency of muscle contraction. *Progress in Biophysics and Molecular Biology, 88*(1), 1-58.

Syed, Christopher K., & Lee, Robert H. (2019). Slow and fast grouping of cargo velocities in axonal transport due to single versus multi-motor transport. *Journal of Theoretical Biology, 480*, 65-70.

Tuckerman, M.E. (2023). *Statistical Mechanics: Theory and Molecular Simulation: Second Edition*. Oxford Graduate Texts. Oxford University Press, new edition, second edition, new edition, second edition ed.

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.

Van Kampen, Nicolaas Godfried (1992). *Stochastic processes in physics and chemistry*, vol. 1.
Elsevier.

Vilfan, A, & Duke, T (2003). Instabilities in the transient response of muscle. *Biophysical Journal*, 85(2), 818-827.
Wagoner, Jason A., & Dill, Ken A. (2021). Evolution of mechanical cooperativity among myosin II motors. *Proceedings of the National Academy of Sciences*, 118(20), e2101871118.

Wijeratne, Sithara S., Fiorenza, Shane A., Neary, Alex E., Subramanian, Radhika, & Betterton, Meredith D. (2022). Motor guidance by long-range communication on the microtubule highway. *Proceedings of the National Academy of Sciences*, 119(28), e2120193119.