# Thermodynamic properties of muscle contraction models and associated discrete-time principles

François Kimmig, Dominique Chapelle, Philippe Moireau

► To cite this version:

François Kimmig, Dominique Chapelle, Philippe Moireau. Thermodynamic properties of muscle con-traction models and associated discrete-time principles. Advanced Modeling and Simulation in Engi-neering Sciences, 2019, 6, 10.1186/s40323-019-0128-9. [hal-02065957v3](https://hal-02065957v3)

HAL Id: hal-02065957

https://inria.hal.science/hal-02065957v3

Submitted on 23 Apr 2019

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Thermodynamic properties of muscle contraction models and associated discrete-time principles

François Kimmig, Dominique Chapelle,

Philippe Moireau

Received: date / Accepted: date

**Abstract** Considering a large class of muscle contraction models accounting for actin-myosin interaction, we present a mathematical setting in which solution properties can be established, including fundamental thermodynamic balances. Moreover, we propose a complete discretization strategy for which we are also able to obtain discrete versions of the thermodynamic balances and other properties. Our major objective is to show how the thermodynamics of such models can be tracked after discretization, including when they are coupled to a macroscopic muscle formulation in the realm of continuum mechanics. Our approach allows to carefully identify the sources of energy and entropy in the system, and to follow them up to the numerical applications.

**Keywords** muscle contraction · sliding filaments · thermodynamically consistent time-discretization · Clausius-Duhem inequality

**Mathematics Subject Classification (2010)** 74F25 · 74H15 · 65M12 · 35Q79 and 92C45

## Introduction

The modeling of the active mechanical behavior of muscles has been the object of intense research since the seminal work of A.F. Huxley [12] modeling the attachment-detachment process in the actin-myosin interaction responsible for sarcomere contraction. Then, numerous extensions – mostly based on refinements of the chemical process introduced by Huxley – of the previous model have been proposed in order to take into account different time scales of the actin-myosin interaction. In particular several models have been developed to account for the power stroke phenomenon [4, 5, 13, 20]. In parallel, the question of the thermodynamic balances associated with the chemical machinery was intensively studied, notably with the fundamental contributions of T.L. Hill [10, 11]. Note that these models are specific cases of molecular motors models without
the natural diffusion introduced by the Fokker-Plank equation [14, 3, 18, 2]. In this pa-
per, our objective is to develop a formalism allowing to derive these thermodynamic
balances for Huxley’s model and its extensions with an additional tracking of these
balances at the discrete level after time-discretizing the model dynamics. Moreover, we
present how these microscopic models can be incorporated into a macroscopic model of
muscle fibers in the spirit of [2] with the aim of following these thermodynamic balances
at the macroscopic level for the continuous-time dynamics but also after adequate time
discretization. This last part is general with respect to the chemical microscopic model
of interest and could also be extended to similar types of models [19, 3], or those mixing
mechanical and chemical modeling elements, for instance [22, 17, 1].

The outline of the paper is as follows. The first section presents the modeling ingre-
dients of the microscopic models of actin-myosin interaction and we derive in a second
section the fundamental properties of these models with the associated thermodynamic
balances, up to the coupling with the macroscopic mechanical formulation. The third
section then describes the discretization scheme and justifies its thermodynamic com-
patibility. Finally, the last section illustrates our results with numerical investigations.

# 1 Modeling of muscle contraction

1.1 Physiology of muscle contraction

Muscles are multi-scale structures in which motion is initiated at the cellular level by
the relative sliding between two types of filaments: actin filaments and myosin filaments.
At the surface of the myosin filament, myosin heads can bind to the actin filament.
The actin filament has a periodic structure with regularly spaced attachment sites.
The interaction between myosin heads and actin sites occurs in a cyclic manner [16],
see Figure 1. The cycle includes attachment and detachment of the myosin head to
and from an actin site and a conformation change of the attached myosin head called
the power stroke. The detachment stage requires an energy input obtained from ATP
molecules buffered inside the cell.

<figure><img src="image_2.png" /><figcaption>Fig. 1: Lymn-Taylor cycle representation. Each stage of the cycle can be seen as a
change of chemical state.</figcaption></figure>

Different levels of description of the actin-myosin interaction can be considered
benefiting from the fact that the power stroke occurs much faster than the attachment
and detachment processes.
## 1.2 Huxley'57 model

In his seminal work [12], Huxley describes the myosin head with two chemical states representing the attached and detached configurations. Each myosin can interact with its closest actin site only. The transition rates – for attachment and detachment – depend only on the distance from the myosin head rest position to its nearest attachment site denoted by $s$. We denote by $d_a$ the distance between two consecutive attachment sites. The distance $s$ thus lies in an interval of width $d_a$, not necessarily symmetric but containing 0, that we denote by $[s^-, s^+]$ – see Figure 2.

Considering, in a population of myosin heads, the subset of heads with rest position located at distance $s$ from their nearest attachment site, we define by $a(t, s)$ the ratio of actually attached heads at time $t$. Equivalently, the ratio of detached heads is denoted by $d(s, t) = 1 - a(s, t)$, due to the assumption that both types of filaments are rigid. The sliding velocity $v_c$ between the filaments is a macroscopic variable, hence independent of $s$ and often quasi-static with respect to the microscopic time scales. We refer to Section 2.5 for an illustrating example of coupling between a macroscopic formulation and such a microscopic model.

The detached state is associated with a constant energy level $w_0$ and the attached state with an energy $w_1$ that depends on the distance $s$ – the myosin head bound to actin is modeled as an elastic spring. This is where mechanics enters the model, and we point out that we extend here the original Huxley'57 model by allowing the spring to have a non-linear behavior. The myosin head is brought back to the initial detached energy level by the ATP energy input $\mu_T$.

Transition rates between the states satisfy the detailed balance, i.e. for a transition $i$ ($i = 1, 2$, see Figure 2) from a state of energy $w_j(s)$ to a state of energy $w_k(s)$, the forward and reverse rates – respectively denoted by $k_i$ and $k_{-i}$ – must satisfy the relation

$$ \frac{k_i(s)}{k_{-i}(s)} = \exp \left( \frac{w_k(s) - w_j(s)}{k_B T} \right), \qquad (1) $$

where $k_B$ is the Boltzmann constant and $T$ is the absolute temperature. A schematic of the model is presented in Figure 2.

The conservation of matter, assuming that there is no coupling between the myosin heads, leads to the following dynamical system, for all $t > 0$ and all $s \in [s^-, s^+]$

$$ \left\{ \begin{aligned} \frac{\partial a(s, t)}{\partial t} &= (k_1(s) + k_{-2}(s))d(s, t) - (k_2(s) + k_{-1}(s))a(s, t) - v_c \frac{\partial a(s, t)}{\partial s} \\ \frac{\partial d(s, t)}{\partial t} &= (k_2(s) + k_{-1}(s))a(s, t) - (k_1(s) + k_{-2}(s))d(s, t) - v_c \frac{\partial d(s, t)}{\partial s} \end{aligned} \right. \qquad (2) $$

from adequate initial conditions $a(s, 0) = a^0(s)$ and $d(s, 0) = d^0(s)$, to be specified later. The assumption that the myosin head can only interact with its nearest actin site imposes that the probability of being attached on the boundaries of the interval $[s^-, s^+]$ must be zero. Physically, the property $a(s^-, t) = a(s^+, t) = 0$ appears when the attachment rates $k_1(s)$ and $k_{-2}(s)$ vanish, while the detachment rates $k_{-1}(s)$ and $k_2(s)$ go to infinity on the boundaries of the interval $[s^-, s^+]$. Note that the energy levels and the transition rates are linked by the detailed balance (1), which implies that the energy of the attached level goes to infinity on the boundaries of the interval
<figure><img src="image_3.png" /><figcaption>Fig. 2: Huxley’57 model representation. Top: definition of the transitions between the attached state (A) and the detached state (D). Bottom left: model parametrization, representation in the detached state. Bottom right: model parametrization, representation in the attached state. The position of the actin site corresponding to s = 0 is represented by a thick dashed line.</figcaption></figure>

$[s^-, s^+]$. In a nutshell, the parameter functions must satisfy

$$ \left| \begin{array}{l} \lim_{s \to s^{-}} k_{1}(s) = \lim_{s \to s^{+}} k_{1}(s) = \lim_{s \to s^{-}} k_{-2}(s) = \lim_{s \to s^{+}} k_{-2}(s) = 0, \\ \\ \lim_{s \to s^{-}} \int_{s}^{0} (k_{-1} + k_{2})ds = \lim_{s \to s^{+}} \int_{0}^{s} (k_{-1} + k_{2})ds = +\infty, \end{array} \right. \qquad (3) $$

the second line enforcing that all heads are detached at the boundaries, see Section 2 (Equations (9) and (10)). This implies energetically that

$$ \lim_{s \to s^{-}} w_{1}(s) = \lim_{s \to s^{+}} w_{1}(s) = +\infty. $$

Actin sites and myosin heads are located at discrete locations separated by regular intervals along their respective filaments. The spatial periodicities are, however, different on each filament. Therefore, for a large population of heads, the distribution of their distance to the nearest actin site can be assumed to be uniform in the interval $[s^-, s^+]$, and the average tension developed per myosin head is given by

$$ \tau_c(t) = \frac{1}{d_a} \int_{s^{-}}^{s^{+}} a(s,t) \frac{\partial w_1}{\partial s}(s) ds. \qquad (4) $$

This force can then typically lead to a macroscopic active stress tensor and link to macro-scale models of muscle tissue as presented in [2] or in Section 2.5.
## 1.3 Extension of Huxley'57 model

To obtain a behavior closer to physiology, and in particular to capture the power stroke, various extensions of Huxley'57 model have been proposed [4, 5, 13, 20]. These extensions can use more than two states to describe the myosin head and allow interactions with an arbitrary number of attachment sites. In this section, our objective is to present these models in a general form, albeit close to the initial 2-state Huxley's model, in particular concerning their general mathematical and mechanical properties.

**Multi-state models** – A general formulation of these models considers $N_s$ chemical states $\{X_p\}_{1\le p\le N_s}$, that we can separate into two categories: attached states and detached states. The states are involved in $N_r$ reactions between the states in the form

$$X_{p_i} \underset{k_{-i}}{\overset{k_i}{\rightleftharpoons}} X_{q_i} \quad (5)$$

The collection of reactions between the states can be represented by a complete directed graph $G$. A complete directed graph is a set of vertices connected by edges, in which: edges have a direction; for each edge of the graph, the edge connecting the same vertices in the inverse direction also belongs to the graph; no vertex is connected to itself. We respectively denote by $V$ and $E$ the sets of vertices and edges of the complete directed graph, and we write $G = (V, E)$. This graph is made of $N_s$ vertices and $2N_r$ edges. A transition $X_{p_i} \to X_{q_i}$ is associated with the edge $E_{p_i q_i}$, whereas the reverse transition $X_{q_i} \to X_{p_i}$ is associated with the edge $E_{q_i p_i}$. The reaction presented in (5), which is bidirectional, is associated with the edges $E_{p_i q_i}$ and $E_{q_i p_i}$. The subsets of vertices corresponding to attached and detached states are respectively denoted by $V_a$ and $V_d$.

The ratio of heads in state $X_q$ located at $s$ at time $t$ is denoted by $x_q(s, t)$. We define the chemical flux between states $X_{p_i}$ and $X_{q_i}$ through transition $i$ by

$$J_{p_i q_i}(s, t) = k_i(s)x_{p_i}(s, t) - k_{-i}(s)x_{q_i}(s, t).$$

Note that we have

$$J_{q_i p_i}(s, t) = -J_{p_i q_i}(s, t) = k_{-i}(s)x_{q_i}(s, t) - k_i(s)x_{p_i}(s, t).$$

The system dynamics is then governed by

$$\frac{\partial x_q}{\partial t}(s, t) = \sum_{p | E_{pq} \in G} J_{pq}(s, t) - v_c \frac{\partial x_q}{\partial s}(s, t), \quad \forall q \in [1, N_s].$$

The Huxley'57 model presented in Section 1.2 can naturally be seen as a particular case of multi-state models with only one attached state and one detached state. The graph $G$ associated with the Huxley'57 model is given by

$$\left\{
\begin{array}{l}
V = \{A, D\}, \\
V_a = \{A\}, \\
V_d = \{D\}, \\
E = \{E_{AD}^1, E_{AD}^2, E_{DA}^1, E_{DA}^2\}.
\end{array}
\right.$$

Here, we use superscripts in the edges definition to denote that there are two reactions between the same vertices.
**Multi-site models** – In this further generalization, it is assumed that a myosin head can interact not only with its nearest actin site – located by definition at distance *s* – but also with all other actin sites located at distance {$s + jd_a$}$_{j \in \mathbb{Z}^*}$. A myosin head can thus be detached in state *q* – with a probability $x_q(s)$ – or attached in state *q* at a distance *s* + *jd*<sub>a</sub> of its rest position – with a probability $x_q(s + jd_a)$. We extend the definition of the ratio of heads in detached states by periodicity, i.e. if state *q* is a detached state $x_q(s + jd_a) = x_q(s)$, $\forall j$. We also refine the description of the graph defined for multi-site models, by splitting the set of edges *E* between the edges linking two detached states $\hat{E}$ and the remaining edges $\bar{E}$. We define the associated complete directed graphs $G_d = (V, \hat{E})$ and $G_a = (V, \bar{E})$. The system dynamics is governed by

$$ egin{cases} 
\frac{\partial x_q}{\partial t}(s+jd_a, t) = \sum_{p|E_{pq} \in G} J_{pq}(s+jd_a, t) - v_c \frac{\partial x_q}{\partial s}(s+jd_a, t), \forall V_q \in V_a, j \in \mathbb{Z}, \\ 
\frac{\partial x_q}{\partial t}(s,t) = \sum_{p|E_{pq} \in G_a} \sum_{j \in \mathbb{Z}} J_{pq}(s+jd_a, t) + \sum_{p|E_{pq} \in G_d} J_{pq}(s,t) - v_c \frac{\partial x_q}{\partial s}(s,t), \forall V_q \in V_d. 
\end{cases} 
\quad (6) $$

The ratios of attached head $x_q$ for $V_q \in V_a$ are defined on $\mathbb{R}$ and must vanish at infinity. By contrast, the ratios of detached head $x_q$ for $V_q \in V_d$ are defined on $[s^-, s^+]$ with periodic boundary conditions.

**Example of a multi-state, multi-site model: Piazzesi-Lombardi'95** – A specific representative – denoted PL95 – of this family of models has been derived in [20] with the aim of accounting for the energetics of muscle contraction. It describes the myosin head with five states arranged in two cycles of chemical reactions, see Figure 3. The five states are composed of three attached states $A_1$, $A_2$ and $A_3$ and two detached states $D_1$ and $D_2$. A first long cycle (cycle a) is meant to represent a complete power stroke, while a short cycle (cycle b) allows the myosin head to cycle at small or zero sliding velocity with incomplete power stroke.

<figure><img src="image_4.png" /><figcaption>Fig. 3: Left: graph associated with the model Piazzesi-Lombardi'95. The vertex indices are given in red. The transition indices are given in blue. Right: Piazzesi-Lombardi'95 model parametrization. The position of the actin site corresponding to $s = 0$ is represented by a thick dashed line.</figcaption></figure>

An energy $\mu_T$ is brought to the myosin head by ATP in the transitions $2 \to 5$ and $3 \to 4$.
The graph G associated with this model is given by

$$
\begin{array}{l}
V = \{V_1, V_2, V_3, V_4, V_5\}, \\
V_a = \{V_1, V_2, V_3\}, \\
V_d = \{V_4, V_5\}, \\
E = \bar{E} = \{E_{12}, E_{21}, E_{23}, E_{32}, E_{34}, E_{43}, E_{41}, E_{14}, E_{25}, E_{52}, E_{51}, E_{15}\}.
\end{array}
$$

Moreover, it is assumed in this model that the myosin can attach to an arbitrary number of actin sites, hence it is also multi-site.

We denote by $w_q$ the energy associated with the state of vertex $q$.

## 2 Model properties based on thermodynamics principles

2.1 From conservation of matter to boundary conditions and monotonicity properties

Let us consider the Huxley'57 model and derive its fundamental properties. System (2) was derived from the conservation of matter, hence we directly verify that for all $(s, t)$, $a(s, t) + d(s, t) = 1$ as soon as we choose our initial condition $\forall s \in [s^{-}, s^{+}]$, $a^{0}(s) \in [0, 1]$ and $d^{0}(s) = 1 - a^{0}(s)$, since

$$
\forall s \in [s^{-}, s^{+}], \quad \frac{d}{dt}(a(s,t) + d(s,t)) = 0,
$$

where we defined the total derivative by $d/dt(\bullet) = \partial/\partial t(\bullet) + v_c\partial/\partial s(\bullet)$. Therefore, we can rewrite the system (2) in the form of a single equation

$$
\frac{\partial a}{\partial t}(s,t) = k_{+}(s)(1 - a(s,t)) - k_{-}(s)a(s,t) - v_{c}\frac{\partial a}{\partial s}(s,t), \quad (7)
$$

where we denote the aggregated transition rates $k_+(s) = k_1(s) + k_{-2}(s)$ and $k_-(s) = k_2(s) + k_{-1}(s)$.

**Boundary values** – As explained in our model presentation, we expect the myosin head to only interact with the nearest actin site, which imposes that the probability of being attached must vanish on the boundaries of the interval $[s^-, s^+]$. However, the dynamics (2) is a first-order transport equation associated with only one boundary condition. Therefore, we can either consider one single Dirichlet boundary condition at one end of the interval – i.e., in $s^-$ if $v_c > 0$ and $s^+$ if $v_c < 0$ – and then rely on the conditions (3) to obtain the proper value of the solution at the other end – as a property – or alternatively consider periodic boundary conditions. As the first option yields a periodic solution, it is clear that the two options are equivalent. However, they differ at the discrete level, in which case we will have to make a choice, see Section 3.1.

In fact, closed-form expressions can be obtained for the solution. To fix the ideas in this derivation, we consider the case $v_c \ge 0$, although the same result can be obtained similarly for $v_c < 0$. As $v_c$ is assumed to be constant, the method of characteristic lines gives regular $C^1$ solutions from regular enough initial condition $a^0$. Considering $a(s, t)$ solution of (7), we define the function $\tilde{a}$ by

$$
\tilde{a}(s,t) = a(s,t) \exp \left( \frac{1}{v_c} \int_{0}^{s} (k_{+}(\xi) + k_{-}(\xi)) d\xi \right),
$$
which satisfies the equation

$$ \frac{\partial \tilde{a}}{\partial t}(s,t) + v_c \frac{\partial \tilde{a}}{\partial s}(s,t) = k_+(s) e^{1/v_c h(s)}, \qquad (8) $$

where we define $h(s) = \int_0^s (k_+(\xi) + k_-(\xi)) d\xi$. Solving (8) along a characteristic line and pulling back the result to $a(s,t)$ we obtain, defining $t_0 = \max(0, t - \frac{s-s_-}{v_c})$,

$$ a(s,t) = a^0(s - v_c t) \exp\left(\frac{-1}{v_c} \int_{s-v_c t}^{s} (k_+(\xi) + k_-(\xi)) d\xi\right) \\ + \int_{t_0}^{t} k_+(s - v_c t + v_c \tau) \exp\left(\frac{-1}{v_c} \int_{s-v_c t + v_c \tau}^{s} (k_+(\xi) + k_-(\xi)) d\xi\right) d\tau. \qquad (9) $$

We know that the aggregated attachment rate $k_+$ is a continuous function on $[s^-, s^+]$ and goes to zero on the boundaries of $[s^-, s^+]$. Therefore, under the condition on the aggregated detachment rate $\int_0^{s^+} k_-(\xi) d\xi = \infty$, we deduce that

$$ \lim_{s \to s^+} a(s,t) = 0 \qquad (10) $$

using the dominated convergence theorem for the second term of (9). Likewise, the property $\lim_{s \to s^-} a(s,t) = 0$ is obtained for $v_c < 0$, and a similar result can be obtained in a similar manner with periodic boundary conditions.

**Positivity and boundedness properties** – We want to check that the solution has values consistent with ratio quantities. More specifically, we want that, with an initial condition $a^0(s) \in [0, 1]$, the property $a(s,t) \in [0, 1]$ holds. Again, we rely on the solution obtained by the method of characteristic lines (9). As the transition rates and the initial condition are positive, we find that $a(s,t) \ge 0$. Then, noting that $1 - a(s,t)$ is governed by an equation of the same form as (7) with the initial condition $1 - a^0(s) \ge 0$, we similarly deduce that $a(s,t) \le 1$.

## 2.2 First principle

We now want to establish a first thermodynamic property of the Huxley'57 system (2), namely, a first principle, and in this respect we follow the approach proposed by [10]. We consider a system made of a population of myosin heads and define the average energy per myosin head, namely

$$ \mathcal{U}(t) = \frac{1}{da} \int_{s^{-}}^{s^{+}} \left[ w_1(s)a(s,t) + w_0d(s,t) \right] ds. \qquad (11) $$

We then define the chemical fluxes

$$ \begin{cases} J_1(s,t) = k_1(s)d(s,t) - k_{-1}(s)a(s,t), \\ J_2(s,t) = k_2(s)a(s,t) - k_{-2}(s)d(s,t). \end{cases} $$

We will henceforth make the natural assumption that the reaction rates are chosen in order for $w_1a$, $k_2a$ and $k_{-1}a$ to tend to zero when $s$ tends to $s^-$ and $s^+$, with the
physical interpretation that no finite energy ($w_1a$) is stored and no detachment flux ($k_2a$ and $k_{-1}a$) occurs at the ends of the interval. We will see in Section 4 that this assumption is easily satisfied in practice when (3) holds. Then, computing the time derivative, we obtain

$$ 
\begin{aligned}
\frac{d}{dt}\mathcal{U}(t) &= \frac{1}{d_a}\int_{s^{-}}^{s^{+}} \left[ w_1(s)\frac{\partial a(s,t)}{\partial t} + w_0\frac{\partial d(s,t)}{\partial t} \right] ds \\
&= \frac{1}{d_a}\int_{s^{-}}^{s^{+}} \left[ (J_1(s,t) - J_2(s,t))w_1(s) - w_1(s)v_c\frac{\partial a(s,t)}{\partial s} \right. \\
&\qquad \left. + (J_2(s,t) - J_1(s,t))w_0 - w_0v_c\frac{\partial d(s,t)}{\partial s} \right] ds.
\end{aligned}
$$

Using integrations by parts for the transport terms, the boundary properties of the solution – $w_1(s)a(t, s^{-}) = w_1(s)a(t, s^{+}) = 0$ and $d(t, s^{-}) = d(t, s^{+})$ – and considering that detachment is associated with the consumption of one ATP, we obtain

$$ 
\begin{aligned}
\frac{d}{dt}\mathcal{U}(t) = \frac{1}{d_a}\int_{s^{-}}^{s^{+}} \left[ (w_1(s) - w_0)J_1(s,t) + (w_0 - (w_1(s) + \mu_T))J_2(s,t) \right. \\
\left. \qquad + v_c a(s,t)\frac{\partial w_1(s)}{\partial s} \right] ds + \mu_T \overline{J_2}(t),
\end{aligned}
$$

where we denoted by

$$ \overline{J_2}(t) = \frac{1}{d_a}\int_{s^{-}}^{s^{+}} J_2(s,t) ds $$

the mean net influx of ATP. Finally, we derive the following formulation of the first principle

$$ \dot{\mathcal{U}}(t) = \dot{W}(t) + \dot{\mathcal{E}}(t) + \dot{\mathcal{Q}}(t), $$

where

$$ \dot{W}(t) = v_c \tau_c(t), \qquad (12a) $$

$$ \dot{\mathcal{E}}(t) = \mu_T \overline{J_2}(t), \qquad (12b) $$

$$ \dot{\mathcal{Q}}(t) = \frac{1}{d_a}\int_{s^{-}}^{s^{+}} \left[ (w_1(s) - w_0)J_1(s,t) + (w_0 - (w_1(s) + \mu_T))J_2(s,t) \right] ds, \qquad (12c) $$

with the active force $\tau_c$ defined in (4). The quantity $\dot{W}$ is the rate of work given to the system and $\dot{\mathcal{E}}(t) = \mu_T \overline{J_2}(t)$ corresponds to the input flux in energy brought by ATP hydrolysis. The remaining term $\dot{\mathcal{Q}}$ can be identified with a heat flux. In steady-state shortening ($\tau_c > 0$ and $v_c < 0$), the work is negative and in physiological conditions, we expect the energy input term to be positive, and the heat transfer to be negative (see numerical illustrations in Section 4). The energy balance can be interpreted as follows: the energy brought by ATP is for one part converted into work, the other part being dissipated as heat production.
## 2.3 Second principle

Let us now derive a second principle thermodynamic balance. We introduce the system entropy as

$$S(t) = -\frac{k_B}{d_a} \int_{s^-}^{s^+} \left[ a(s,t) \ln(a(s,t)) + d(s,t) \ln(d(s,t)) \right] ds,$$

where $k_B$ is the Bolztmann constant. The system remains at a constant temperature, the outside environment playing the role of a thermostat. We introduce the Helmholtz free energy

$$F(t) = U(t) - TS(t),$$

which therefore corresponds to

$$F(t) = \frac{1}{d_a} \int_{s^-}^{s^+} \left[ \mu_1(s,t)a(s,t) + \mu_0(s,t)d(s,t) \right] ds. \quad (13)$$

where

$$\left| \begin{array}{l} \mu_1(s,t) = w_1(s) + k_B T \ln(a(s,t)), \\ \mu_0(s,t) = w_0 + k_B T \ln(d(s,t)), \end{array} \right.$$

are the chemical potentials. Computing the time derivative, we get

$$\frac{d}{dt} F(t) = \frac{1}{d_a} \int_{s^-}^{s^+} \left[ \left( J_1(s,t) - J_2(s,t) - v_c \frac{\partial}{\partial s} a(s,t) \right) \cdot \left[ w_1(s) + k_B T \ln(a(s,t)) + k_B T \right] \right. \\ \left. + \left( J_2(s,t) - J_1(s,t) - v_c \frac{\partial}{\partial s} d(s,t) \right) \cdot \left[ w_0 + k_B T \ln(d(s,t)) + k_B T \right] \right] ds. \quad (14)$$

Using again integrations by part for the transport terms and the boundary properties of the solution, we obtain

$$\frac{d}{dt} F(t) = \tau_c(t)v_c + \mu_T \overline{J_2}(t) + \frac{1}{d_a} \int_{s^-}^{s^+} \left[ J_2(s,t) \left( \mu_0(s,t) - (\mu_1(s,t) + \mu_T) \right) \right. \\ \left. + J_1(s,t) \left( \mu_1(s,t) - \mu_0(s,t) \right) \right] ds. \quad (15)$$

Independently, we also have

$$\frac{d}{dt} F(t) = \frac{d}{dt} U(t) - T \frac{d}{dt} S(t). \quad (16)$$

Combining equations (15) and (16) and the first principle, we can write

$$\frac{d}{dt} S(t) = \frac{\dot{Q}(t)}{T} - \frac{1}{T} \frac{1}{d_a} \int_{s^-}^{s^+} \left[ J_2(s,t) \left( \mu_0(s,t) - (\mu_1(s,t) + \mu_T) \right) \right. \\ \left. + J_1(s,t) \left( \mu_1(s,t) - \mu_0(s,t) \right) \right] ds. \quad (17)$$
The second principle then reads

$$ \frac{d}{dt} S(t) = \frac{\dot{Q}(t)}{T} + \dot{S}_{\text{prod}}(t), \qquad (18) $$

where we naturally associate the second term of (17) with the entropy production $\dot{S}_{\text{prod}}(t)$.

The model will be compatible with the second principle if this entropy production is always positive. Using the relation (1) deduced from the detailed balance, we recall that

$$ \frac{k_1(s, t)}{k_{-1}(s, t)} = \exp\left(\frac{w_0 - w_1(s, t)}{k_B T}\right) \quad \frac{k_2(s, t)}{k_{-2}(s, t)} = \exp\left(\frac{w_1(s, t) + \mu_T - w_0}{k_B T}\right). \quad (19) $$

Thus, when introducing the ratio of the one-way fluxes for transition 1, $J_{1+}$ and $J_{1-}$, defined by $J_{1+}(s, t) = k_1(s, t)d(s, t)$ and $J_{1-}(s, t) = k_{-1}(s, t)a(s, t)$, we find

$$ \begin{aligned} \frac{J_{1+}}{J_{1-}} &= \frac{k_1(s)d(s, t)}{k_{-1}(s)a(s, t)} \\ &= \exp\left(\frac{w_0 - w_1(s, t)}{k_B T}\right) \cdot \exp\left(k_B T \left[\frac{\ln(d(s, t)) - \ln(a(s, t))}{k_B T}\right]\right) \quad (20) \\ &= \exp\left(\frac{\mu_0(s, t) - \mu_1(s, t)}{k_B T}\right). \end{aligned} $$

As a consequence, we have two cases. If $\mu_1(s, t) \ge \mu_0(s, t)$, we find that

$$ \frac{J_{1+}}{J_{1-}}(s, t) \le 1 \Rightarrow J_{1+}(s, t) = J_{1+}(s, t) - J_{1-}(s, t) \le 0. $$

Conversely, $\mu_1(s, t) \le \mu_0(s, t)$ implies that $J_{1+}(s, t) \ge 0$. Proceeding in the same way for the second reaction, we finally have

$$ \left| \begin{array}{l} J_{1}(s, t)(\mu_{1}(s, t) - \mu_{0}(s, t)) \le 0, \\ J_{2}(s, t)(\mu_{0}(s, t) - \mu_{1}(s, t) - \mu_{T}) \le 0. \end{array} \right. \qquad (21) $$

We thus obtain the conclusion that the entropy production term

$$ \dot{S}_{\text{prod}}(t) = -\frac{1}{T d_a} \int_{s^{-}}^{s^{+}} J_{2}(s, t) (\mu_{0}(s, t) - (\mu_{1}(s, t) + \mu_{T})) \\ + J_{1}(s, t) (\mu_{1}(s, t) - \mu_{0}(s, t)) \, ds \ge 0, \quad (22) $$

hence that the model is compatible with the second principle. We can summarize this property using (15) by the free energy balance

$$ \frac{d}{dt} F(t) = \tau_c(t)v_c + \mu_T \overline{J}_2(t) - T \dot{S}_{\text{prod}}(t) \le \tau_c(t)v_c + \mu_T \overline{J}_2(t). \quad (23) $$
## 2.4 Extension to multi-state, multi-site models

Let us now consider the Piazzesi-Lombardi'95 model. We want to establish the thermodynamic balances associated with this model.

Fundamental properties of the solution – in particular the monotonicity properties of the results established for the Huxley'57 model, see Section 2.1 – can be extended to this model. In particular, conservation of matter here reads

$$ \forall s \in [s^{-}, s^{+}] \quad \frac{d}{dt} \left[ \sum_{V_q \in V_a} \sum_{j \in \mathbb{Z}} x_q(s + jd_a, t) + \sum_{V_q \in V_d} x_q(s, t) \right] = 0. \quad (24) $$

*First principle* – We define the energy as

$$ U(t) = \frac{1}{d_a} \int_{s^{-}}^{s^{+}} \left[ \sum_{V_q \in V_a} \sum_{j \in \mathbb{Z}} x_q(s + jd_a, t) w_q(s + jd_a) + \sum_{V_q \in V_d} x_q(s, t) w_q \right] ds. $$

The time derivative reads, after integrating by parts the transport term

$$ \frac{d}{dt} U(t) = \frac{1}{d_a} \int_{s^{-}}^{s^{+}} \sum_{E_{pq} \in G} \sum_{j \in \mathbb{Z}} J_{pq}(s + jd_a, t) W_{pq}(s + jd_a) ds + v_c \tau_c(t), $$

with $W_{pq} = w_q - w_p$ and the active force defined as

$$ \tau_c(t) = \frac{1}{d_a} \int_{s^{-}}^{s^{+}} \sum_{V_q \in V_a} \sum_{j \in \mathbb{Z}} x_q(s + jd_a, t) \frac{\partial w_q}{\partial s}(s + jd_a) ds. $$

We then obtain

$$ \frac{d}{dt} U(t) = \mu_T (\overline{J_{25}}(t) + \overline{J_{34}}(t)) + v_c \tau_c(t) \\ + \frac{1}{d_a} \int_{s^{-}}^{s^{+}} \sum_{E_{pq} \in G} \sum_{j \in \mathbb{Z}} J_{pq}(s + jd_a, t) \tilde{W}_{pq}(s + jd_a) ds, $$

where the mean fluxes are given by

$$ \overline{J_{pq}}(t) = \frac{1}{d_a} \int_{s^{-}}^{s^{+}} \sum_{j \in \mathbb{Z}} J_{pq}(s + jd_a, t) ds, $$

and $\tilde{W}_{pq} = W_{pq} - \mu_T \delta_{pq} = \{25,34\}$ as an energy $\mu_T$ is brought to the myosin head by ATP during the transitions $2 \to 5$ and $3 \to 4$. Note that the introduction of the modified energy increments $\tilde{W}$ brings out the input energy fluxes $\mu_T(\overline{J_{25}}(t) + \overline{J_{34}}(t))$. The first principle then naturally reads

$$ \dot{U}(t) = \dot{W}(t) + \dot{\mathcal{E}}(t) + \dot{\mathcal{Q}}(t), \quad (25) $$
with

$$
\left|
\begin{aligned}
\dot{W}(t) &= v_c \tau_c(t), \\
\dot{\mathcal{E}}(t) &= \mu_T (\overline{J_{25}}(t) + \overline{J_{34}}(t)), \\
\dot{Q}(t) &= \frac{1}{d_a} \int_{s^{-}}^{s^{+}} \sum_{E_{pq} \in G} \sum_{j \in \mathbb{Z}} J_{pq}(s + j d_a, t) \tilde{W}_{pq}(s + j d_a) \, ds.
\end{aligned}
\right.
$$

Second principle – Following the work done for the Huxley’s7 model, we here define
the entropy of the system as

$$
\mathcal{S}(t) = -\frac{k_B}{d_a} \int_{s^{-}}^{s^{+}} \sum_{V_q \in V_a} \sum_{j \in \mathbb{Z}} x_q(s + jd_a, t) \ln(x_q(s + jd_a, t)) \\
\qquad + \sum_{V_q \in V_d} x_q(s, t) \ln(x_q(s, t)) \, ds.
$$

Then, we define the Helmholtz free energy as $\mathcal{F}(t) = \mathcal{U}(t) - TS(t)$, which can be
written as

$$
\mathcal{F}(t) = \frac{1}{d_a} \int_{s^{-}}^{s^{+}} \sum_{V_q \in V_a} \sum_{j \in \mathbb{Z}} \mu_q(s + jd_a, t) x_q(s + jd_a, t) + \sum_{V_q \in V_d} \mu_q(s, t) x_q(s, t) \, ds,
$$

with the definitions

$$
\left|
\begin{array}{l}
\mu_q(s + jd_a, t) = w_q(s + jd_a) + k_B T \ln(x_q(s + jd_a, t)), \quad \forall V_q \in V_a, j \in \mathbb{Z}, \\
\mu_q(s, t) = w_q(s) + k_B T \ln(x_q(s, t)), \quad \forall V_q \in V_d.
\end{array}
\right.
$$

We here define the entropy production as

$$
\dot{\mathcal{S}}_{\text{prod}}(t) = -\frac{1}{Td_a} \int_{s^{-}}^{s^{+}} \sum_{E_{pq} \in G} \sum_{j \in \mathbb{Z}} J_{pq}(s + jd_a, t) \tilde{M}_{pq}(s + jd_a, t) \, ds,
$$

where $\tilde{M}_{pq} = \mu_q - \mu_p - \mu_T \delta_{pq} = \{34,25\}$. Then, combining the first principle (25) with
the identity $\dot{F} = \dot{U} - T\dot{S}$, we finally obtain the second principle

$$
\frac{d}{dt} S(t) = \frac{\dot{Q}(t)}{T} + \dot{S}_{\text{prod}}(t).
$$

Similarly as in the Huxley’s7 model, the detailed balance ensures that the entropy
production is always positive and the model is thus thermodynamically compatible.
We have, as in (23), the free energy balance

$$
\begin{align}
\frac{d}{dt} \mathcal{F}(t) ={}& \tau_c(t)v_c + \mu_T(\overline{J_{25}}(t) + \overline{J_{34}}(t)) - T\dot{\mathcal{S}}_{\text{prod}}(t) \nonumber \\
& \leq \tau_c(t)v_c + \mu_T(\overline{J_{25}}(t) + \overline{J_{34}}(t)). \tag{26}
\end{align}
$$
<figure><img src="image_5.png" /><figcaption>Fig. 4: Muscle fiber configuration.</figcaption></figure>

## 2.5 Coupling with a macroscopic model of muscle fiber

The thermodynamic properties of these classes of models are very useful when coupling them with a macroscopic model, typically to represent a muscle fiber, as it will ensure a global consistent thermodynamic balance between macroscopic and microscopic contributions. Let us consider, indeed, a macroscopic model of muscle fiber modeled in the realm of non-linear continuum mechanics, as large deformations frequently occur in muscle fibers. The material points coordinates are denoted by $\underline{x} \in \Omega^0$ in the reference configuration. The displacement field associated with the deformation map is denoted by $\underline{y}$. We denote by $\underline{e}$ the Green-Lagrange strain tensor, i.e.

$$ \underline{e} = \frac{1}{2} (\underline{\nabla} \underline{y} + (\underline{\nabla} \underline{y})^{\mathrm{T}} + (\underline{\nabla} \underline{y})^{\mathrm{T}} \cdot \underline{\nabla} \underline{y}), $$

and the second Piola-Kirchhoff stress tensor is denoted by $\underline{\underline{\Sigma}}$. The fiber as shown in Figure 4 is subjected to a boundary force $\underline{t}_N$ on a boundary $\underline{\Gamma}_N^0$. The principle of virtual work (PVW) then reads: for any admissible virtual displacement field $\underline{w} \in \mathcal{V}_{ad}$,

$$ \int_{\Omega^0} \rho_0 \ddot{\underline{y}} \cdot \underline{w} \, d\underline{\Omega} + \int_{\Omega^0} \underline{\underline{\Sigma}} : \underline{d}_y \underline{e} \cdot \underline{w} \, d\underline{\Omega} = \int_{\underline{\Gamma}_N^0} \underline{t}_N \cdot \underline{w} \, d\underline{\Gamma}, $$

where the differential of the Green-Lagrange strain tensor with respect to the displacement field is given by

$$ \underline{d}_y \underline{e} \cdot \underline{w} = \frac{1}{2} (\underline{\nabla} \underline{w} + (\underline{\nabla} \underline{w})^{\mathrm{T}} + (\underline{\nabla} \underline{y})^{\mathrm{T}} \cdot \underline{\nabla} \underline{w} + (\underline{\nabla} \underline{w})^{\mathrm{T}} \cdot \underline{\nabla} \underline{y}). $$

In this formulation, we want to associate with each material point an active microscopic model based on the Huxley'57 model or its extensions. Typically, we want to incorporate the microscopic model into a 3D visco-hyperelastic constitutive behavior of hyperelastic potential $\Psi$ and viscous pseudo-potential $\Psi^v$ – taken here as $\Psi^v(\dot{\underline{e}}) = \frac{\eta}{2} \mathrm{tr}(\dot{\underline{e}}^2)$ – where $\eta$ denotes a viscosity modulus – to simplify the presentation. Following [2], which extends the classical Hill-Maxwell scheme [9] to nonlinear behavior, we gather all the constitutive ingredients by defining an adequate rheological scheme – presented in Figure 5 – valid for large deformations. The upper branch represents the sarcomere – including the above active behavior visualized by the collection of myosin heads in the figure – namely, constituents acting in the muscle fibre direction $\underline{\tau}$. The lower branch represents a 3D passive matrix, associated with the cellular envelope and the extracellular matrix. Each branch contains elastic and viscous constituents, respectively visualized by springs and dashpots, with specific constitutive equations given below.
<figure><img src="image_6.png" /><figcaption>Fig. 5: Fiber rheology combining a 1D active element (upper branch) and a 3D visco-hyperelastic element (lower branch).</figcaption></figure>

We consider the following natural rheological rule for the parallel branch

$$ \text{3D parallel law : } \underline{e} = \underline{e}_p = \underline{e}_a, \quad \underline{\Sigma} = \underline{\Sigma}_p + \underline{\Sigma}_a, \qquad (27) $$

where $\underline{e}, \underline{e}_a, \underline{e}_p$ denote Green-Lagrange tensors (global, active and passive) and $\underline{\Sigma}, \underline{\Sigma}_a, \underline{\Sigma}_p$ second Piola-Kirchhoff stress tensors (global, active and passive). However, we will depart from [2] for the series branch.

In fact, the natural view of muscle fibers made of a succession of active and passive segments points to a one-dimensional homogenization type of rheological interpretation. Let us denote by $e_{\text{fib}}$ the total (local) extension of a fiber, i.e the ratio of length change over initial length

$$ e_{\text{fib}} = \frac{\delta \ell_{\text{fib}}}{\ell_{\text{fib}}}, $$

where we will take for $\ell_{\text{fib}}$ the length of a half-sarcomere at rest and $\delta\ell_{\text{fib}}$ the variation thereof. Then, the length change of the half-sarcomere can be decomposed into

$$ \delta\ell_{\text{fib}} = \delta\ell_c + \delta\ell_s, $$

where $\delta\ell_c$ and $\delta\ell_s$ respectively denote the contributions of the active (rigid filaments) and passive parts in this length change. Note that $\delta\ell_c$ then represents the relative displacement of the actin and myosin filaments considered above, indeed. We also introduce the corresponding dimensionless extension quantities $e_c = \delta\ell_c/\ell_{\text{fib}}$ and $e_s = \delta\ell_c/\ell_{\text{fib}}$, so that $e_{\text{fib}} = e_c + e_s$. Of course, the two components carry the same tension, which we denote by $T_{\text{fib}}$ and define as the force – in the fiber direction – per unit area of transverse cross-section of tissue considered in the reference configuration. Therefore, we can summarize as

$$ \text{1D series law : } e_{\text{fib}} = e_c + e_s, \quad T_{\text{fib}} = T_c = T_s. \qquad (28) $$
Note that the reason why such a simple additive rule holds for $e_{\text{fib}}$ in this nonlinear framework is that we are considering extension quantities – scaled in an ad hoc manner – and not Green-Lagrange strains. Moreover, (28) must be complemented by relationships between 3D and 1D quantities. Considering the component of the Green-Lagrange strain tensor in the fiber direction, we directly have

$$1 + e_{\text{fib}} = (1 + 2\underline{\underline{\tau}} \cdot \underline{\underline{e}} \cdot \underline{\underline{\tau}})^{\frac{1}{2}}. \qquad (29)$$

Then, the tension $T_{\text{fib}}$ corresponds to a contribution in the first Piola-Kirchhoff stress tensor given by

$$\underline{T} = T_{\text{fib}} \frac{\underline{\underline{F}} \cdot \underline{\underline{\tau}}}{\|\underline{\underline{F}} \cdot \underline{\underline{\tau}}\|} \otimes \underline{\underline{\tau}}, \qquad (30)$$

where $\underline{\underline{F}} = \underline{\underline{1}} + \underline{\underline{\nabla}} y$ is the classical deformation gradient tensor, as can be easily verified by computing the resulting traction in the fiber direction $\underline{\underline{T}} \cdot \underline{\underline{\tau}}$. Hence, the associated contribution in the second Piola-Kirchhoff stress tensor reads

$$\underline{\underline{\Sigma}}_a = \underline{\underline{F}}^{-1} \cdot \underline{\underline{T}}_a = \frac{T_{\text{fib}}}{1 + e_{\text{fib}}} \underline{\underline{\tau}} \otimes \underline{\underline{\tau}}, \qquad (31)$$

since $\|\underline{\underline{F}} \cdot \underline{\underline{\tau}}\| = 1 + e_{\text{fib}}$.

Finally, the constitutive equations considered are

$$\underline{\underline{\Sigma}}_p = \frac{\partial \underline{\underline{\Psi}}}{\partial \underline{\underline{e}}} + \frac{\partial \underline{\underline{\Psi}}^v}{\partial \underline{\underline{\dot{e}}}}, \quad T_s = E_s e_s, \quad T_c = \nu \dot{e}_c + \check{T}_c, \qquad (32)$$

where $\check{T}_c$ represents the aggregation of forces contributed by actin-myosin cross-bridges as described above, i.e. $\check{T}_c = \rho_{\text{surf}} \tau_c$ with $\rho_{\text{surf}}$ the number of myosin heads in a layer of thickness $\ell_{\text{fib}}$ per unit of cross-section area. The series elastic element is here assumed to have a linear constitutive equation of elasticity modulus $E_s$. Note that a nonlinear hyperelastic behavior could be considered, at the price of having to deal with the dimensionless extension $e_s$ as an additional internal variable. Nevertheless, in physiological conditions the extension $e_s$ remains small and a linear behavior is adequate [1]. As regards viscosity, we here incorporate a simple component of viscous modulus $\nu$ in parallel with the active part in the sarcomere branch, and we recall that viscosity is also present in the parallel branch as provided by the term $\frac{\partial \underline{\underline{\Psi}}^v}{\partial \underline{\underline{\dot{e}}}} = \eta \dot{\underline{\underline{e}}}$. We can now summarize the 3D equations as

$$\int_{\Omega^0} \rho_0 \ddot{\underline{\underline{y}}} \cdot \underline{\underline{w}} \, d\Omega + \int_{\Omega^0} \underline{\underline{\Sigma}} : \underline{\underline{d}}_{\underline{\underline{y}}} \underline{\underline{e}} \cdot \underline{\underline{w}} \, d\Omega = \int_{\Gamma_N^0} \underline{\underline{t}}_N \cdot \underline{\underline{w}} \, d\Gamma, \quad \forall \underline{\underline{w}} \in \mathcal{V}_{\text{ad}} \qquad (33a)$$

$$\text{with } \underline{\underline{\Sigma}} = \frac{\partial \underline{\underline{\Psi}}}{\partial \underline{\underline{e}}} + \eta \dot{\underline{\underline{e}}} + \frac{T_{\text{fib}}}{(1 + 2\underline{\underline{\tau}} \cdot \underline{\underline{e}} \cdot \underline{\underline{\tau}})^{\frac{1}{2}}} \underline{\underline{\tau}} \otimes \underline{\underline{\tau}}$$

$$T_{\text{fib}} = \nu \dot{e}_c + \check{T}_c = E_s e_s \qquad (33b)$$

$$\text{with } \check{T}_c(\underline{x}, t) = \frac{\rho_{\text{surf}}}{d_a} \int_{s^{-}}^{s^{+}} a(\underline{x}, s, t) \frac{\partial w_1}{\partial s}(s, t) ds$$

$$\frac{\partial a}{\partial t}(\underline{x}, s, t) = k_{+}(s)(1 - a(\underline{x}, s, t)) - k_{-}(s)a(\underline{x}, s, t) - \ell_{\text{fib}} \dot{e}_c(\underline{x}, t) \frac{\partial a}{\partial s}(\underline{x}, s, t) \qquad (33c)$$

where (33c) is based on the Huxley'57 model with sliding velocity $v_c = \ell_{\text{fib}} \dot{e}_c$. This velocity is independent of the microscopic variable $s$, which justifies our above study.
Note, however, the dependency of $a(\underline{x}, s, t)$ on $\underline{x}$, which means that the microscopic model must be solved everywhere in the domain, i.e. at all numerical quadrature points in numerical simulations.

In order to establish a macroscopic energy balance for the system (33), we take the velocity field $\dot{\underline{y}}$ as an admissible virtual displacement field in (33a). We thus get

$$ \frac{d\mathcal{K}}{dt} + \int_{\Omega^0} \underline{\underline{\Sigma}} : \dot{\underline{\underline{\epsilon}}} d\Omega = \mathcal{P}_{\text{ext}}, $$

where $\mathcal{K} = \frac{1}{2} \int_{\Omega^0} \rho_0 |\dot{\underline{y}}|^2 d\Omega$ stands for the kinetic energy and $\mathcal{P}_{\text{ext}} = \int_{\Gamma_N^0} t_N \cdot \dot{\underline{y}} d\Gamma$ is the power of external forces. Then, we decompose

$$ \int_{\Omega^0} \underline{\underline{\Sigma}} : \dot{\underline{\underline{\epsilon}}} d\Omega = \frac{d}{dt} \left[ \int_{\Omega^0} \Psi d\Omega \right] + \int_{\Omega^0} \eta |\dot{\underline{\underline{\epsilon}}}|^2 d\Omega + \int_{\Omega^0} T_{\text{fib}} \frac{\underline{\underline{\tau}} \cdot \dot{\underline{\underline{\epsilon}}} \cdot \underline{\underline{\tau}}}{(1 + 2\underline{\underline{\tau}} \cdot \underline{\underline{\epsilon}} \cdot \underline{\underline{\tau}})^{\frac{1}{2}}} d\Omega, $$

where the first term is associated with the stored hyperelastic energy, the second term is a macroscopic viscous dissipation, and the last term – denoted $\mathcal{P}_{\text{fib}}$ – is the power of internal forces in the sarcomere. Then, using the rheological rules we find

$$ \begin{aligned} \mathcal{P}_{\text{fib}} &= \int_{\Omega^0} T_{\text{fib}} \frac{\underline{\underline{\tau}} \cdot \dot{\underline{\underline{\epsilon}}} \cdot \underline{\underline{\tau}}}{(1 + 2\underline{\underline{\tau}} \cdot \underline{\underline{\epsilon}} \cdot \underline{\underline{\tau}})^{\frac{1}{2}}} d\Omega \\ &= \int_{\Omega^0} T_{\text{fib}} \dot{\underline{\underline{\epsilon}}}_{\text{fib}} d\Omega = \int_{\Omega^0} [T_s \dot{\underline{\underline{\epsilon}}}_s + T_c \dot{\underline{\underline{\epsilon}}}_c] d\Omega = \int_{\Omega^0} [E_s e_s \dot{\underline{\underline{\epsilon}}}_s + \nu \dot{\underline{\underline{\epsilon}}}_c^2 + \tilde{T}_c \dot{\underline{\underline{\epsilon}}}_c] d\Omega \\ &= \frac{d}{dt} \left[ \underbrace{\int_{\Omega^0} \frac{E_s}{2} e_s^2 d\Omega}_{(1)} \right] + \underbrace{\int_{\Omega^0} \nu \dot{\underline{\underline{\epsilon}}}_c^2 d\Omega}_{(2)} + \underbrace{\int_{\Omega^0} \tilde{T}_c \dot{\underline{\underline{\epsilon}}}_c d\Omega}_{(3)}, \end{aligned} $$

where we recognize (1) an elastic energy stored in the series element of the sarcomere, (2) a viscous dissipation term in the sarcomere, and (3) the mechanical work of the actin-myosin bridges. Therefore, combining this energy balance with the free energy balance (23) computed from the Huxley'57 model – or identically from (26) for the extensions – we finally obtain a form of macroscopic Clausius-Duhem relation

$$ \begin{aligned} & \frac{d}{dt} \left[ \mathcal{K} + \int_{\Omega^0} \Psi + \frac{E_s}{2} e_s^2 + \rho_v \mathcal{F} d\Omega \right] \\ &= \mathcal{P}_{\text{ext}} - \int_{\Omega^0} \left[ \eta |\dot{\underline{\underline{\epsilon}}}|^2 + \nu \dot{\underline{\underline{\epsilon}}}_c^2 \right] d\Omega - \int_{\Omega^0} \rho_v T \dot{\mathcal{S}}_{\text{prod}} d\Omega + \int_{\Omega^0} \rho_v \dot{\mathcal{E}} d\Omega, \end{aligned} \quad (34) $$

where $\rho_v = \rho_{\text{surf}}/\ell_{\text{fib}}$ is the density of myosin head per unit volume in the reference configuration, and where we recall that $\mathcal{F}$ is the internal free energy of the bridges introduced in (13), $\dot{\mathcal{S}}_{\text{prod}}$ is the entropy production term defined in (22) corresponding to energy dissipation associated with chemical transitions, and $\dot{\mathcal{E}} = \mu_T \bar{J}_2$ as defined in (12b) corresponds to the input flux in energy provided by ATP hydrolysis.
### 3 Discretization and thermodynamic principles at discrete level

We now present the proposed discretization scheme for the muscle contraction models. Classical schemes are sufficient for our purposes, and the main originality of this work is to show their compatibility with discrete versions of the thermodynamical principles. Nevertheless, for the sake of completeness, some basic properties of the schemes are quickly re-established before focusing on thermodynamics.

#### 3.1 A numerical scheme for Huxley'57 model

To discretize the dynamics (2), we consider a regular grid for the simulation range $[s^-, s^+]$ of discretization length $\delta s$ and with the convention $s_0 = s^-$ and $s_\ell = s^+$. We then choose an upwind implicit scheme, that for the sake of simplicity we only present for a positive sliding velocity $v_c$, with a natural extension to negative sliding velocities $v_c$ by inverting the shift in space for the transport term to keep an upwind scheme. We initiate the discretization from an initial condition such that $a_0^0 = a_\ell^0 = 0$ and with the natural condition $a_i^0 \in [0, 1]$. The discretization scheme then reads

$$ \left\{ \begin{array}{ll} \frac{a_i^{n+1} - a_i^n}{\delta t} = k_{+,i} d_i^{n+1} - k_{-,i} a_i^{n+1} - v_c \frac{a_i^{n+1} - a_{i-1}^{n+1}}{\delta s}, & \forall i \in [1, \ell] \\ d_i^{n+1} = 1 - a_i^{n+1}, & \forall i \in [1, \ell] \end{array} \right. \quad (35) $$

with the definition

$$ |k_{+,i} = k_+(s^- + i\delta s), \quad i \in [1, \ell], \qquad (36a) $$

$$ |k_{-,i} = k_-(s^- + i\delta s), \quad i \in [1, \ell - 1], \qquad (36b) $$

$$ |k_{-,\ell} = 2k_{-,\ell-1}. \qquad (36c) $$

Note that the exact aggregated detachment rate goes to infinity on the boundary of the interval $[s^-, s^+]$. Numerically, we use a finite value defined as given in (36c), and we prove in the following section that this choice does not affect the convergence of the scheme.

For the numerical scheme (35), we also need to prescribe adequate boundary conditions. As the analytical solution of (7) vanishes on the boundaries of the interval $[s^-, s^+]$, we here again can choose: either a Dirichlet condition on one side, and check the consistency on the other side, or choose periodic boundary conditions and again ensure the consistency on the boundary of the interval $[s^-, s^+]$. From a numerical point of view, it is in fact more convenient for energy estimates to choose periodic boundary conditions for $a$, i.e. $a_0^n = a_\ell^n$. Note that, with this choice, we do not strictly have $a_0^n = a_\ell^n = 0$. This property is only satisfied approximately, or asymptotically when the spatial discretization length goes to zero.

Defining $\alpha = v_c\delta t/\delta s$ and $k_i = k_{+,i} + k_{-,i}$, the scheme can be written in a matrix form on the state vector $\mathbf{a}^n = [a_1^n \dots a_\ell^n]^T$

$$ \mathbb{D}\mathbf{a}^{n+1} = \mathbf{a}^n + \delta t \mathbf{k}_+, $$
where

$$ \mathbb{D} = \begin{bmatrix} 1 + \delta tk_1 + \alpha & & & -\alpha \\ -\alpha & 1 + \delta tk_2 + \alpha & & \\ & \ddots & \ddots & \\ & & -\alpha & 1 + \delta tk_\ell + \alpha \end{bmatrix}, \quad \mathbf{k}_+ = \begin{bmatrix} k_{+,1} \\ k_{+,2} \\ \vdots \\ k_{+,\ell} \end{bmatrix}. $$

## 3.2 Some fundamentals properties

We first present the basic – but essential – properties of the proposed scheme. This is done using classical strategies for the analysis of transport equations schemes (see for instance [21]).

**Uniform positivity and boundedness** – One first important property that must be satisfied by the discretization is that the natural bounds for ratio quantities be preserved at the discrete level, namely $\forall n, \forall i \in [1, \ell]$, $a_i^n \in [0, 1]$. To obtain this property, we need $\mathbb{D}$ to preserve the positivity – i.e. for $\mathbf{a} \in \mathbb{R}^\ell$, $\mathbb{D}\mathbf{a} \ge 0 \Rightarrow \mathbf{a} \ge 0$ (where we use the convention that a vector is positive if all its coefficients $a_i$ are positive). Let us then take $\mathbf{a} \in \mathbb{R}^\ell$ such that $\mathbb{D}\mathbf{a} \ge 0$. We have $\forall i \in [1, \ell]$

$$ (1 + \delta t k_i + \alpha) a_i - \alpha a_{i-1} \ge 0, \quad (37) $$

with the boundary condition $a_0 = a_\ell$. Multiplying (37) by $\alpha^{l-1}$ for $i = 1$ and by $\alpha^{l-1} \prod_{j=1}^{i-1} (1 + k_j + \alpha)$ for $i \in [2, \ell]$, and summing, we obtain

$$ \begin{aligned} & \sum_{i=2}^{\ell} \left[ \alpha^{\ell-i} \prod_{j=1}^{i} (1+k_j+\alpha)a_i - \alpha^{\ell-i+1} \prod_{j=1}^{i-1} (1+k_j+\alpha)a_{i-1} \right] \\ & \qquad + \alpha^{l-1}(1+k_1+\alpha)a_1 - \alpha^l a_l \ge 0 \\ \Leftrightarrow & \quad \prod_{j=1}^{\ell} (1+k_j+\alpha)a_\ell + \sum_{i=3}^{\ell} \alpha^{\ell-i+1} \left[ \prod_{j=1}^{i-1} (1+k_j+\alpha)a_{i-1} \right] \\ & \quad - \sum_{i=2}^{\ell} \alpha^{\ell-i+1} \left[ \prod_{j=1}^{i-1} (1+k_j+\alpha)a_{i-1} \right] + \alpha^{l-1}(1+k_1+\alpha)a_1 - \alpha^l a_l \ge 0. \end{aligned} $$

Noting that the middle term is a telescoping series, we obtain as expected

$$ \underbrace{\left[ \prod_{j=1}^{\ell} (1 + k_j + \alpha) - \alpha^{\ell} \right]}_{\ge 0} a_{\ell} \ge 0. $$

Then, recursively from (37), we get $\forall i \in [1, \ell]$, $a_i \ge 0$, which shows that the matrix operator $\mathbb{D}$ preserves the positivity. Knowing that the initial condition and the transition rates are positive, we obtain $\forall n \ge 0, \forall i \in [1, \ell]$, $a_i^n \ge 0$. Writing the numerical scheme for the variable $1 - a_i^n$ from (35), we similarly obtain that $\forall n \ge 0, \forall i \in [1, \ell]$, $a_i^n \le 1$. Therefore, the proposed numerical scheme preserves the adequate positivity and boundedness.
Consistency – Let us here denote $\bar{a}$ a sufficiently regular solution of (7). Note that this solution satisfies a Dirichlet boundary condition on one side of the simulation interval (i.e. $\bar{a}(s^{-}, t) = 0$) and that we showed in (10) that $\lim_{s \to s^{+}} \bar{a}(s, t) = 0$. We denote by $\bar{a}$ the vector of the values of $\bar{a}$ at the spatial discretization points at time $n \delta t$, $\bar{a}^{n} = [\bar{a}(\delta s, n \delta t) \dots \bar{a}(\ell \delta s, n \delta t)]^{\mathrm{T}}$. We define, as usual, the convergence error by

$$e^n = \mathbf{a}^n - \bar{\mathbf{a}}^n.$$

We have $\mathbb{D}e^{n+1} = e^n - \delta t \eta^n$, where the consistency error $\eta^n$ is given by

$$\eta^n = \frac{1}{\delta t} [\mathbb{D}\bar{\mathbf{a}}^{n+1} - (\bar{\mathbf{a}}^n + \delta t \mathbf{k}_+)].$$

Evaluating the continuous equation (7) at $s = s^{-} + i\delta s$ and $t = (n + 1)\delta t$, we obtain

$$\frac{\partial \bar{a}}{\partial t}(s^{-} + i\delta s, (n+1)\delta t) = k_{+,i} - k_{i}\bar{a}_{i}^{n+1} - v_{c}\frac{\partial \bar{a}}{\partial s}(s^{-} + i\delta s, (n+1)\delta t),$$

and this also holds for $i = \ell$ with the finite numerical value chosen for $k_{\ell}$, due to the fact that $\bar{a}(s^{+}) = 0$ and $k_{-}(s)\bar{a}(s)$ tends to zero in $s^{+}$. We directly infer, $\forall i \in [1, \ell]$, $\eta_{i}^{n} = \frac{\bar{a}(s^{-} + i\delta s, (n + 1)\delta t) - \bar{a}(s^{-} + i\delta s, n\delta t)}{\delta t} - \frac{\partial \bar{a}}{\partial t}(s^{-} + i\delta s, (n + 1)\delta t) + v_{c}\left\{\frac{\bar{a}(s^{-} + i\delta s, (n + 1)\delta t) - \bar{a}(s^{-} + (i - 1)\delta s, (n + 1)\delta t)}{\delta s} - \frac{\partial \bar{a}}{\partial s}(s^{-} + i\delta s, (n + 1)\delta t)\right\} = O(\delta s + \delta t)$.

$$\eta_i^n = \frac{\bar{a}(s^{-} + i\delta s, (n + 1)\delta t) - \bar{a}(s^{-} + i\delta s, n\delta t)}{\delta t} - \frac{\partial \bar{a}}{\partial t}(s^{-} + i\delta s, (n + 1)\delta t) + v_c \left\{ \frac{\bar{a}(s^{-} + i\delta s, (n + 1)\delta t) - \bar{a}(s^{-} + (i - 1)\delta s, (n + 1)\delta t)}{\delta s} - \frac{\partial \bar{a}}{\partial s}(s^{-} + i\delta s, (n + 1)\delta t) \right\} = O(\delta s + \delta t).$$

The numerical scheme is thus consistent at the first order with the continuous equation (7), although different – but compatible – boundary conditions are used. Note that a well-known result by Godunov states that we cannot have more than first-order convergence in time with a discrete scheme that satisfies the positivity and boundedness property [6].

$L^2$-Stability – We now analyze the $L^2$-stability by justifying the $\ell^2$-stability of the operator $\mathbb{D}^{-1}$. Multiplying (35) without the source term by $a_i^{n+1}$, we have

$$\frac{a_i^{n+1} - a_i^n}{\delta t} a_i^{n+1} + v_c \frac{a_i^{n+1} - a_{i-1}^{n+1}}{\delta s} a_i^{n+1} = -k_i (a_i^{n+1})^2.$$

Using the identity $-ab = \frac{1}{2}(a-b)^2 - \frac{1}{2}a^2 - \frac{1}{2}b^2$, we obtain

$$\frac{1}{2\delta t} \left( (a_i^{n+1})^2 - (a_i^n)^2 \right) + \frac{1}{2\delta t} \left( a_i^{n+1} - a_i^n \right)^2 = -\frac{v_c}{2\delta s} \left( (a_i^{n+1})^2 - (a_{i-1}^{n+1})^2 \right) - \frac{v_c}{2\delta s} \left( a_i^{n+1} - a_{i-1}^{n+1} \right)^2 - k_i (a_i^{n+1})^2.$$

Summing over $i$, and using the periodic boundary conditions $(a_{\ell}^{n+1})^2 = (a_0^{n+1})^2$, we find

$$\| \mathbb{D}^{-1} \mathbf{a}^n \|_{\ell^2} \le \| \mathbf{a}^n \|_{\ell^2}.$$

Note that this stability property is just a mathematical property.
Convergence – The stability analysis coupled to the consistency analysis gives directly the convergence error. Indeed, we find

$$e^n = (\mathbb{D}^{-1})^n e_0 - \delta t \sum_{k=0}^{n-1} (\mathbb{D}^{-1})^{n-k-1} \eta^k,$$

so that there exist $C, T > 0$ such that

$$\|e^n\|_{\ell^2} \le C \delta t \sum_{k=0}^{n-1} (\delta s + \delta t) = CT(\delta s + \delta t).$$

### 3.3 First principle

Our objective is more ambitious than numerical convergence, as we want *in fine* to establish thermodynamic balances at the discrete level. In this respect, let us first consider the energy balance. We recall that the average energy of a myosin head is given at the continuous level by (11). Similarly to (36), we assign a finite value to the energy of the attached state on the boundary of the interval $[s^-, s^+]$. With the notation

$$\begin{cases} w_{1,i} = w_1(s^- + i\delta s), \forall i \in [1, l-1], \\ w_{1,l} = 2w_{1,l-1}, \end{cases}$$

the energy is discretized as

$$\mathcal{U}^n = \frac{\delta s}{d_a} \sum_{i=1}^{\ell} w_{1,i} a_i^n + w_0 d_i^n.$$

Then, defining the fluxes as

$$\begin{cases} J_{1,i}^n = k_{1,i} d_i^n - k_{-1,i} a_i^n, \\ J_{2,i}^n = k_{2,i} a_i^n - k_{-2,i} d_i^n, \end{cases}$$

we have

$$\begin{aligned} \frac{\mathcal{U}^{n+1} - \mathcal{U}^n}{\delta t} &= \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \left[ w_{1,i} \left( \frac{a_i^{n+1} - a_i^n}{\delta t} \right) + w_0 \left( \frac{d_i^{n+1} - d_i^n}{\delta t} \right) \right] \\ &= \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \left[ w_{1,i} \left( J_{1,i}^{n+1} - J_{2,i}^{n+1} - v_c \frac{a_i^{n+1} - a_{i-1}^{n+1}}{\delta s} \right) \right. \\ & \qquad \left. + w_0 \left( J_{2,i}^{n+1} - J_{1,i}^{n+1} + v_c \frac{a_i^{n+1} - a_{i-1}^{n+1}}{\delta s} \right) \right] \\ &= \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \left[ J_{1,i}^{n+1} (w_{1,i} - w_0) + J_{2,i}^{n+1} (w_0 - w_{1,i} - \mu_T) \right. \\ & \qquad \left. + \mu_T J_{2,i}^{n+1} \right] - \frac{v_c \delta s}{d_a} \sum_{i=1}^{\ell-1} \frac{w_{1,i} - w_{1,i+1}}{\delta s} a_i^{n+1} + \frac{\delta s v_c}{d_a \delta s} a_0^{n+1} w_{1,1} \\ & \qquad - \frac{\delta s v_c}{d_a \delta s} a_{\ell}^{n+1} w_{1,\ell} + \frac{\delta s v_c}{d_a \delta s} w_0 (a_{\ell}^{n+1} - a_0^{n+1}). \end{aligned} \quad (38)$$
With the definition of the discrete force

$$ \tau_c^{n+1} = \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \frac{w_{1,i+1} - w_{1,i}}{\delta s} a_i^{n+1}, $$

where we define $w_{1,\ell+1} = w_{1,1}$, and using the periodicity of the solution, (38) becomes

$$ \frac{\mathcal{U}^{n+1} - \mathcal{U}^n}{\delta t} = v_c \tau_c^{n+1} + \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \left[ J_{1,i}^{n+1} (w_{1,i} - w_0) + J_{2,i}^{n+1} (w_0 - w_{1,i} - \mu_T) + \mu_T J_{2,i}^{n+1} \right]. $$

We thus obtain the discretized version of the first principle, namely

$$ \frac{\mathcal{U}^{n+1} - \mathcal{U}^n}{\delta t} = \frac{\mathcal{W}^{n+1} - \mathcal{W}^n}{\delta t} + \frac{\mathcal{Q}^{n+1} - \mathcal{Q}^n}{\delta t} + \frac{\mathcal{E}^{n+1} - \mathcal{E}^n}{\delta t}, \quad (39) $$

with

$$ \left| \begin{array}{l} \frac{\mathcal{W}^{n+1} - \mathcal{W}^n}{\delta t} = v_c \tau_c^{n+1}, \\ \\ \frac{\mathcal{Q}^{n+1} - \mathcal{Q}^n}{\delta t} = \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \left[ J_{1,i}^{n+1} (w_{1,i} - w_0) + J_{2,i}^{n+1} (w_0 - w_{1,i} - \mu_T) \right], \\ \\ \frac{\mathcal{E}^{n+1} - \mathcal{E}^n}{\delta t} = \mu_T \frac{\delta s}{d_a} \sum_{i=1}^{\ell} J_{2,i}^{n+1}. \end{array} \right. $$

## 3.4 Second principle

Let us now establish a discrete entropy balance. In this respect, we introduce the discrete entropy

$$ S^n = -k_B \frac{\delta s}{d_a} \sum_{i=1}^{\ell} (a_i^n \ln a_i^n + d_i^n \ln d_i^n), $$

and the free energy $\mathcal{F}^n = \mathcal{U}^n - TS^n$, which can be rewritten as

$$ \mathcal{F}^n = \frac{\delta s}{d_a} \sum_{i=1}^{\ell} (\mu_{1,i}^n a_i^n + \mu_{0,i}^n d_i^n), $$

by introducing the discrete chemical potentials

$$ \left| \begin{array}{l} \mu_{1,i}^n = w_{1,i} + k_B T \ln a_i^n, \\ \mu_{0,i}^n = w_0 + k_B T \ln d_i^n. \end{array} \right. $$
We then rewrite the previous calculation in a manner that closely follows the calculation in the continuous case. We have

$$ \frac{\mathcal{F}^{n+1} - \mathcal{F}^n}{\delta t} = \frac{1}{\delta t} \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \left[ \mu_{1,i}^{n+1} a_i^{n+1} - \mu_{1,i}^n a_i^n + \mu_{0,i}^{n+1} d_i^{n+1} - \mu_{0,i}^n d_i^n \right] \\ = \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \left[ \mu_{1,i}^{n+1} \left( \frac{a_i^{n+1} - a_i^n}{\delta t} \right) + a_i^n \left( \frac{\mu_{1,i}^{n+1} - \mu_{1,i}^n}{\delta t} \right) + \mu_{0,i}^{n+1} \left( \frac{d_i^{n+1} - d_i^n}{\delta t} \right) + d_i^n \left( \frac{\mu_{0,i}^{n+1} - \mu_{0,i}^n}{\delta t} \right) \right]. $$

Hence,

$$ \frac{\mathcal{F}^{n+1} - \mathcal{F}^n}{\delta t} = \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \left[ \mu_{1,i}^{n+1} \left( J_{1,i}^{n+1} - J_{2,i}^{n+1} - v_c \frac{a_i^{n+1} - a_{i-1}^{n+1}}{\delta s} \right) \right. \\ + a_i^n \left( \frac{w_{1,i} + k_B T \ln(a_i^{n+1}) - (w_{1,i} + k_B T \ln(a_i^n))}{\delta t} \right) \\ + \mu_{0,i}^{n+1} \left( J_{2,i}^{n+1} - J_{1,i}^{n+1} + v_c \frac{a_i^{n+1} - a_{i-1}^{n+1}}{\delta s} \right) \\ \left. + d_i^n \left( \frac{w_0 + k_B T \ln(d_i^{n+1}) - (w_0 + k_B T \ln(d_i^n))}{\delta t} \right) \right]. $$

Developing the expression of the chemical potentials, performing Abel transformations and using the periodic boundary conditions, we obtain

$$ \frac{\mathcal{F}^{n+1} - \mathcal{F}^n}{\delta t} = v_c \tau_c^{n+1} \\ + \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \left[ J_{1,i}^{n+1} (\mu_{1,i}^{n+1} - \mu_{0,i}^{n+1}) + J_{2,i}^{n+1} (\mu_{0,i}^{n+1} - \mu_{1,i}^{n+1} - \mu_T) + \mu_T J_{2,i}^{n+1} \right] \\ + \frac{v_c}{d_a} k_B T \sum_{i=1}^{\ell} \left[ a_i^{n+1} \left( \ln(a_{i+1}^{n+1}) - \ln(a_i^{n+1}) \right) \right] \\ + \frac{v_c}{d_a} k_B T \sum_{i=1}^{\ell} \left[ d_i^{n+1} \left( \ln(d_{i+1}^{n+1}) - \ln(d_i^{n+1}) \right) \right] \\ + \frac{k_B T \delta s}{d_a} \sum_{i=1}^{\ell} \left[ a_i^n \left( \frac{\ln(a_i^{n+1}) - \ln(a_i^n)}{\delta t} \right) \right] + \frac{k_B T \delta s}{d_a} \sum_{i=1}^{\ell} \left[ d_i^n \left( \frac{\ln(d_i^{n+1}) - \ln(d_i^n)}{\delta t} \right) \right]. $$
Since $x \mapsto \ln x$ is a concave function, we have $\ln a_j - \ln a_i \le \ln'(a_i)(a_j - a_i)$ so that

$$ \begin{aligned} \frac{\mathcal{F}^{n+1} - \mathcal{F}^n}{\delta t} &\le v_c \tau_c^{n+1} \\ &+ \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \left[ J_{1,i}^{n+1} (\mu_{1,i}^{n+1} - \mu_{0,i}^{n+1}) + J_{2,i}^{n+1} (\mu_{0,i}^{n+1} - \mu_{1,i}^{n+1} - \mu_T) + \mu_T J_{2,i}^{n+1} \right] \\ &+ \frac{v_c}{d_a} k_B T \sum_{i=1}^{\ell} [a_{i+1}^{n+1} - a_i^{n+1}] + \frac{v_c}{d_a} k_B T \sum_{i=1}^{\ell} [d_{i+1}^{n+1} - d_i^{n+1}] \\ &+ \frac{k_B T \delta s}{d_a \delta t} \sum_{i=1}^{\ell} [a_i^{n+1} - a_i^n] + \frac{k_B T \delta s}{d_a \delta t} \sum_{i=1}^{\ell} [d_i^{n+1} - d_i^n]. \end{aligned} $$

Using the fact that the scheme imposes $\forall n, \forall i \in [1, \ell]$, $a_i^n + d_i^n = 1$, the sums vanish two by two. We finally find

$$ \frac{\mathcal{F}^{n+1} - \mathcal{F}^n}{\delta t} \le v_c \tau_c^{n+1} \\ \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \left[ J_{1,i}^{n+1} (\mu_{1,i}^{n+1} - \mu_{0,i}^{n+1}) + J_{2,i}^{n+1} (\mu_{0,i}^{n+1} - \mu_{1,i}^{n+1} - \mu_T) + \mu_T J_{2,i}^{n+1} \right]. $$

As a point-wise evaluation of the continuous expression (21), we have $\forall i \in [1, \ell]$

$$ \left| \begin{aligned} J_{1,i}^{n+1} (\mu_{1,i}^{n+1} - \mu_{0,i}^{n+1}) &\le 0, \\ J_{2,i}^{n+1} (\mu_{0,i}^{n+1} - \mu_{1,i}^{n+1} - \mu_T) &\le 0. \end{aligned} \right. $$

Hence, in our case where $v_c > 0$, we finally obtain

$$ \frac{\mathcal{F}^{n+1} - \mathcal{F}^n}{\delta t} - v_c \tau_c^{n+1} - \mu_T \frac{\delta s}{d_a} \sum_{i=1}^{\ell} J_{2,i}^{n+1} \le 0. \quad (40) $$

To relate the decrease of the free energy to the creation of entropy, we first note

$$ \begin{aligned} \frac{\mathcal{F}^{n+1} - \mathcal{F}^n}{\delta t} - v_c \tau_c^{n+1} - \mu_T \frac{\delta s}{d_a} \sum_{i=1}^{\ell} J_{2,i}^{n+1} &= \frac{\mathcal{F}^{n+1} - \mathcal{F}^n}{\delta t} \\ &\quad - \left[ \frac{U^{n+1} - U^n}{\delta t} - \frac{Q^{n+1} - Q^n}{\delta t} \right] \\ &= -T \frac{S^{n+1} - S^n}{\delta t} + \frac{Q^{n+1} - Q^n}{\delta t}. \end{aligned} \quad (41) $$

Comparing (41) with the formal expression of the second principle (18), we define the discrete entropy creation by

$$ \frac{S_{\text{prod}}^{n+1} - S_{\text{prod}}^n}{\delta t} = \frac{S^{n+1} - S^n}{\delta t} - \frac{1}{T} \frac{Q^{n+1} - Q^n}{\delta t} \ge 0. \quad (42) $$
Note that the entropy creation is formally given by

$$ \frac{S_{\text{prod}}^{n+1} - S_{\text{prod}}^{n}}{\delta t} = -\frac{1}{T} \left( \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \left[ J_{1,i}^{n+1} (\mu_{1,i}^{n+1} - \mu_{0,i}^{n+1}) + J_{2,i}^{n+1} (\mu_{0,i}^{n+1} - \mu_{1,i}^{n+1} - \mu_T) \right] \right. \\ \left. + \frac{v_c k_B T}{d_a} \sum_{i=1}^{\ell} \left[ \left( \ln(d_i^{n+1}) - \ln(a_i^{n+1}) \right) \left[ a_i^{n+1} - a_{i-1}^{n+1} \right] \right] \right. \\ \left. + \frac{k_B T \delta s}{d_a} \sum_{i=1}^{\ell} \left[ a_i^n \left( \frac{\ln(a_i^{n+1}) - \ln(a_i^n)}{\delta t} \right) \right] \right. \\ \left. + \frac{k_B T \delta s}{d_a} \sum_{i=1}^{\ell} \left[ d_i^n \left( \frac{\ln(d_i^{n+1}) - \ln(d_i^n)}{\delta t} \right) \right] \right). $$

and finally combining (40), (41) and (42) we have the time-discrete counterpart of the free energy balance (23)

$$ \frac{\mathcal{F}^{n+1} - \mathcal{F}^n}{\delta t} = v_c \tau_c^{n+1} + \mu_T \frac{\delta s}{d_a} \sum_{i=1}^{\ell} J_{2,i}^{n+1} - T \frac{S_{\text{prod}}^{n+1} - S_{\text{prod}}^n}{\delta t} \quad (43) \\ \le v_c \tau_c^{n+1} + \mu_T \frac{\delta s}{d_a} \sum_{i=1}^{\ell} J_{2,i}^{n+1}. $$

## 3.5 Extension to multi-state, multi-site models

The case of multiple states and sites derives from the same principles, hence justifying that we developed precisely the computations for the Huxley'57 model. The developments are, however, not straightforward because the multi-site assumption implies an infinite number of attachment and detachment fluxes, which has to be properly integrated into the discrete thermodynamical balances. Indeed, in the case of a positive sliding velocity $v_c$, we discretize the system (6) with the following implicit upwind numerical scheme

$$ \left\{ \begin{array}{l} \frac{x_{q,i+j\ell}^{n+1} - x_{q,i+j\ell}^{n}}{\delta t} = \sum_{p|E_{pq} \in G} J_{pq,i+j\ell}^{n+1} - v_c \frac{x_{q,i+j\ell}^{n+1} - x_{q,i+j\ell-1}^{n+1}}{\delta s}, \forall V_q \in V_a, \\ \\ \frac{x_{q,i}^{n+1} - x_{q,i}^{n}}{\delta t} = \sum_{p|E_{pq} \in G} \sum_{j \in \mathbb{Z}} J_{pq,i+j\ell}^{n+1} - v_c \frac{x_{q,i}^{n+1} - x_{q,i-1}^{n+1}}{\delta s}, \forall V_q \in V_d, i \in [1, \ell] \end{array} \right. \quad (44) $$

with the notation $x_{q,i+j\ell}^n = x_q(s^- + (i+j\ell)\delta s, n\delta t)$ for $V_q \in V_a$, $x_{q,i}^n = x_q(s^- + i\delta s, n\delta t)$ for $V_q \in V_d$ and $J_{pq,i+j\ell}^n = J_{pq}(s^- + (i+j\ell)\delta s, n\delta t)$. Note that $x_q$ for $V_q \in V_d$ has periodic boundary conditions and thus $x_{q,i+j\ell} = x_{q,i}, \forall j \in \mathbb{Z}$, and we keep $j \in \mathbb{Z}$, albeit in practice we will bound the attachment zone, introducing a boundary consistency error.
Mass conservation – Defining the total quantity of matter

$$m_i^n = \sum_{V_q \in V_a} \sum_{j \in \mathbb{Z}} x_{q,i+j\ell}^n + \sum_{V_q \in V_d} x_{q,i}^n,$$

– which has the periodicity $m_0^n = m_\ell^n$ – and using the scheme (44), we find the classical implicit transport equation

$$\frac{m_i^{n+1} - m_i^n}{\delta t} = -v_c \frac{m_i^{n+1} - m_{i-1}^{n+1}}{\delta s},$$

hence, we retrieve, as in the Huxley’s7 model, the conservation of matter

$$\sum_{V_q \in V_a} \sum_{j \in \mathbb{Z}} x_{q,i+j\ell}^n + \sum_{V_q \in V_d} x_{q,i}^n = 1 \quad \forall n, \forall i \in [1, \ell]. \quad (45)$$

First principle – Concerning the energy balance, the internal energy is now defined as

$$\mathcal{U}^n = \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \left[ \sum_{V_q \in V_a} \sum_{j \in \mathbb{Z}} w_{q,i+j\ell} x_{q,i+j\ell}^n + \sum_{V_q \in V_d} w_q x_{q,i}^n \right].$$

and we find this time, with the notation $\tilde{W}_{pq,i+j\ell} = \tilde{W}_{pq}(s^{-} + (i + j\ell)\delta s)$, that

$$\frac{\mathcal{U}^{n+1} - \mathcal{U}^n}{\delta t} = \frac{\delta s}{d_a} \left[ \sum_{i=1}^{\ell} \sum_{E_{pq} \in G} \sum_{j \in \mathbb{Z}} J_{pq,i+j\ell}^{n+1} \tilde{W}_{pq,i+j\ell} \right. \\
+ \mu_T \sum_{i=1}^{\ell} \sum_{j \in \mathbb{Z}} \left( J_{25,i+j\ell} + J_{34,i+j\ell} \right) \\
- \frac{v_c \delta s}{d_a} \sum_{i=1}^{\ell} \sum_{V_q \in V_a} \sum_{j \in \mathbb{Z}} w_{q,i+j\ell} \frac{x_{q,i+j\ell}^{n+1} - x_{q,i+j\ell-1}^{n+1}}{\delta s} \\
\left. - \frac{v_c \delta s}{d_a} \sum_{i=1}^{\ell} \sum_{V_q \in V_d} w_q \frac{x_{q,i}^{n+1} - x_{q,i-1}^{n+1}}{\delta s} \right].$$

The last term vanishes with the periodic boundary conditions. Performing an Abel transformation on the penultimate term and defining the discrete force as

$$\tau_c^{n+1} = \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \sum_{V_q \in V_a} \sum_{j \in \mathbb{Z}} x_{q,i+j\ell} \frac{w_{q,i+j\ell+1}^{n+1} - w_{q,i+j\ell}^{n+1}}{\delta s},$$

we obtain the discrete first principle

$$\frac{\mathcal{U}^{n+1} - \mathcal{U}^n}{\delta t} = \frac{\mathcal{W}^{n+1} - \mathcal{W}^n}{\delta t} + \frac{\mathcal{Q}^{n+1} - \mathcal{Q}^n}{\delta t} + \frac{\mathcal{E}^{n+1} - \mathcal{E}^n}{\delta t},$$

with

$$\left| \frac{\mathcal{W}^{n+1} - \mathcal{W}^n}{\delta t} \right. = v_c \tau_c^{n+1}, \\
\frac{\mathcal{Q}^{n+1} - \mathcal{Q}^n}{\delta t} = \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \sum_{E_{pq} \in G} \sum_{j \in \mathbb{Z}} J_{pq,i+j\ell}^{n+1} \tilde{W}_{pq,i+j\ell}, \\
\left. \frac{\mathcal{E}^{n+1} - \mathcal{E}^n}{\delta t} \right. = \mu_T \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \sum_{j \in \mathbb{Z}} \left( J_{25,i+j\ell} + J_{34,i+j\ell} \right).$$
Second principle – We now define the discrete entropy as

$$S^n = -k_B \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \left[ \sum_{V_q \in V_a} \sum_{j \in \mathbb{Z}} x_{q,i+j\ell}^n \ln x_{q,i+j\ell}^n + \sum_{V_q \in V_d} x_{q,i}^n \ln x_{q,i}^n \right],$$

and the free energy $\mathcal{F}^n = \mathcal{U}^n - TS^n$, that we rewrite as

$$\mathcal{F}^n = \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \left[ \sum_{V_q \in V_a} \sum_{j \in \mathbb{Z}} \mu_{q,i+j\ell}^n x_{q,i+j\ell}^n + \sum_{V_q \in V_d} \mu_{q,i}^n x_{q,i}^n \right],$$

using discrete chemical potentials

$$\left| \begin{array}{l} \mu_{q,i+j\ell}^n = w_{q,i+j\ell} + k_B T \ln x_{q,i+j\ell}^n, \quad V_q \in V_a, \\ \mu_{p,i}^n = w_q + k_B T \ln x_{q,i}^n, \quad V_q \in V_d. \end{array} \right.$$

The discrete time derivative of the free energy is

$$\frac{\mathcal{F}^{n+1} - \mathcal{F}^n}{\delta t} = \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \sum_{V_q \in V_a} \sum_{j \in \mathbb{Z}} \left[ \mu_{q,i+j\ell}^{n+1} \left( \frac{x_{q,i+j\ell}^{n+1} - x_{q,i+j\ell}^n}{\delta t} \right) \right. \\ \left. + x_{q,i+j\ell}^n \left( \frac{\mu_{q,i+j\ell}^{n+1} - \mu_{q,i+j\ell}^n}{\delta t} \right) \right] + \sum_{i=1}^{\ell} \sum_{V_q \in V_d} \left[ \mu_{q,i}^{n+1} \frac{x_{q,i}^{n+1} - x_{q,i}^n}{\delta t} + x_{q,i}^n \frac{\mu_{q,i}^{n+1} - \mu_{q,i}^n}{\delta t} \right].$$

Using the notation $\tilde{M}_{pq,i+j\ell}^n = \mu_{q,i+j\ell}^n - \mu_{p,i+j\ell}^n - \mu_T \delta_{pq} = \{25,34\}$ and the results of the first principle, we get

$$\frac{\mathcal{F}^{n+1} - \mathcal{F}^n}{\delta t} = \frac{\delta s}{d_a} \left[ \sum_{i=1}^{\ell} \sum_{E_{pq} \in G} \sum_{j \in \mathbb{Z}} J_{pq,i+j\ell}^{n+1} \tilde{M}_{pq,i+j\ell}^{n+1} \right. \\ + \mu_T \sum_{i=1}^{\ell} \sum_{j \in \mathbb{Z}} \left( J_{25,i+j\ell} + J_{34,i+j\ell} \right) \\ - v_c \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \sum_{V_q \in V_a} \sum_{j \in \mathbb{Z}} k_B T \ln (x_{q,i+j\ell}^{n+1}) \frac{x_{q,i+j\ell}^{n+1} - x_{q,i+j\ell-1}^{n+1}}{\delta s} \\ - v_c \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \sum_{V_q \in V_d} k_B T \ln (x_{q,i}^{n+1}) \frac{x_{q,i}^{n+1} - x_{q,i-1}^{n+1}}{\delta s} \\ + \frac{k_B T \delta s}{d_a} \sum_{i=1}^{\ell} \sum_{V_q \in V_a} \sum_{j \in \mathbb{Z}} \left[ x_{q,i+j\ell}^n \left( \frac{\ln(x_{q,i+j\ell}^{n+1}) - \ln(x_{q,i+j\ell}^n)}{\delta t} \right) \right] \\ \left. + \frac{k_B T \delta s}{d_a} \sum_{i=1}^{\ell} \sum_{V_q \in V_d} \left[ x_{q,i}^n \left( \frac{\ln(x_{q,i}^{n+1}) - \ln(x_{q,i}^n)}{\delta t} \right) \right] + v_c \tau_c^{n+1}. \right]$$
Performing an Abel transformation on the transport terms, we obtain

$$
\begin{aligned}
\frac{\mathcal{F}^{n+1} - \mathcal{F}^n}{\delta t} ={}& \frac{\delta s}{d_a} \left[ \sum_{i=1}^{\ell} \sum_{E_{pq} \in G} \sum_{j \in \mathbb{Z}} J_{pq, i+j\ell}^{n+1} \tilde{M}_{pq, i+j\ell}^{n+1} \right. \\
& \left. \qquad + \mu_T \sum_{i=1}^{\ell} \sum_{j \in \mathbb{Z}} (J_{25, i+j\ell} + J_{34, i+j\ell}) \right] \\
& + \frac{v_c}{d_a} k_B T \sum_{i=1}^{\ell} \sum_{V_q \in V_a} \sum_{j \in \mathbb{Z}} \left[ x_{q, i+j\ell}^{n+1} (\ln(x_{q, i+j\ell+1}^{n+1}) - \ln(x_{q, i+j\ell}^{n+1})) \right] \\
& + \frac{v_c}{d_a} k_B T \sum_{i=1}^{\ell} \sum_{V_q \in V_d} \left[ x_{q, i}^{n+1} (\ln(x_{q, i+1}^{n+1}) - \ln(x_{q, i}^{n+1})) \right] \\
& + \frac{k_B T \delta s}{d_a} \sum_{i=1}^{\ell} \sum_{V_q \in V_a} \sum_{j \in \mathbb{Z}} \left[ x_{q, i+j\ell}^{n} \left( \frac{\ln(x_{q, i+j\ell}^{n+1}) - \ln(x_{q, i+j\ell}^{n})}{\delta t} \right) \right] \\
& + \frac{k_B T \delta s}{d_a} \sum_{i=1}^{\ell} \sum_{V_q \in V_d} \left[ x_{q, i}^{n} \left( \frac{\ln(x_{q, i}^{n+1}) - \ln(x_{q, i}^{n})}{\delta t} \right) \right] + v_c \tau_c^{n+1}.
\end{aligned}
$$

Using again that $\ln x_j - \ln x_i \le \ln'(x_i)(x_j - x_i)$, the detailed balance and the conservation of matter (45), we finally obtain

$$ \frac{\mathcal{F}^{n+1} - \mathcal{F}^n}{\delta t} - v_c \tau_c^{n+1} - \mu_T \sum_{i=1}^{\ell} \sum_{j \in \mathbb{Z}} (J_{25, i+j\ell} + J_{34, i+j\ell}) \le 0. \quad (46) $$

As in Section (3.4) this property is equivalent to

$$ \frac{S_{\text{prod}}^{n+1} - S_{\text{prod}}^{n}}{\delta t} = \frac{S^{n+1} - S^{n}}{\delta t} - \frac{1}{T} \frac{Q^{n+1} - Q^{n}}{\delta t} \geq 0. \quad (47) $$

Moreover, the entropy production is in fact given by

$$
\begin{aligned}
\frac{S_{\text{prod}}^{n+1} - S_{\text{prod}}^{n}}{\delta t} ={}& -\frac{1}{T} \left( \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \sum_{E_{pq} \in G} \sum_{j \in \mathbb{Z}} J_{pq, i+j\ell}^{n+1} \tilde{M}_{pq, i+j\ell}^{n+1} \right. \\
& \left. \qquad - \frac{v_c k_B T}{d_a} \sum_{i=1}^{\ell} \sum_{V_q \in V_a} \sum_{j \in \mathbb{Z}} \left[ \ln(x_{q, i+j\ell}^{n+1}) [x_{q, i+j\ell}^{n+1} - x_{q, i+j\ell-1}^{n+1}] \right] \right. \\
& \left. \qquad - \frac{v_c k_B T}{d_a} \sum_{i=1}^{\ell} \sum_{V_q \in V_d} \left[ \ln(x_{q, i}^{n+1}) [x_{q, i}^{n+1} - x_{q, i-1}^{n+1}] \right] \right. \\
& \left. \qquad + \frac{k_B T \delta s}{d_a} \sum_{i=1}^{\ell} \sum_{V_q \in V_a} \sum_{j \in \mathbb{Z}} \left[ x_{q, i+j\ell}^{n} \left( \frac{\ln(x_{q, i+j\ell}^{n+1}) - \ln(x_{q, i+j\ell}^{n})}{\delta t} \right) \right] \right. \\
& \left. \qquad + \frac{k_B T \delta s}{d_a} \sum_{i=1}^{\ell} \sum_{V_q \in V_d} \left[ x_{q, i}^{n} \left( \frac{\ln(x_{q, i}^{n+1}) - \ln(x_{q, i}^{n})}{\delta t} \right) \right] \right).
\end{aligned}
$$
and, by recombining (46) and (47), we finally get

$$ \frac{\mathcal{F}^{n+1} - \mathcal{F}^n}{\delta t} - v_c \tau_c^{n+1} - \mu_T \sum_{i=1}^{\ell} \sum_{j \in \mathcal{Z}} (J_{25,i+j\ell} + J_{34,i+j\ell}) + T \frac{S_{\text{prod}}^{n+1} - S_{\text{prod}}^n}{\delta t} = 0. \quad (48) $$

## 3.6 Discretization of the macroscopic model coupling

We can derive a full discretized version of the macroscopic model presented in Figure 4 and modeled by the dynamics (33). Here, we will rely – as in [2] – on mid-point rules for the discretization of the PVW, with additional corrections in order to guarantee the energy balance. Therefore, we will typically denote $\underline{v}^{n+\frac{1}{2}} = \frac{v^{n+1}+v^n}{2}$ for any variable $v$ and the notation $n + \frac{1}{2}\#^n$ will allow to indicate when we depart from this classical rule. First, as recommended in [7, 15, 8, 2], we will consider the following non-standard – albeit classical – mid-point quantities

$$ \underline{e}^{n+\frac{1}{2}\#^n} = \underline{e}(\underline{y}^{n+\frac{1}{2}}), \quad \dot{\underline{e}}^{n+\frac{1}{2}\#^n} = \frac{\underline{e}^{n+1} - \underline{e}^n}{\delta t}, \\ d_y \underline{e}^{n+\frac{1}{2}\#^n} \cdot \underline{w} = \frac{1}{2}(\underline{\nabla} \underline{w} + (\underline{\nabla} \underline{w})^T + (\underline{\nabla} \underline{y}^{n+\frac{1}{2}})^T \cdot \underline{\nabla} \underline{w} + (\underline{\nabla} \underline{w})^T \cdot \underline{\nabla} \underline{y}^{n+\frac{1}{2}}), $$

and a passive hyperelastic stress law discretization that includes an energy correction term, namely

$$ \left. \frac{\partial \Psi}{\partial \underline{e}} \right|^{n+\frac{1}{2}\#} = \frac{\partial \Psi}{\partial \underline{e}}(\underline{e}^{n+\frac{1}{2}\#}) + \left( \frac{\Psi(\underline{e}^{n+1}) - \Psi(\underline{e}^n)}{\delta t} - \frac{\partial \Psi}{\partial \underline{e}}(\underline{e}^{n+\frac{1}{2}\#}) : \dot{\underline{e}}^{n+\frac{1}{2}\#} \right) \frac{\dot{\underline{e}}^{n+\frac{1}{2}\#}}{\dot{\underline{e}}^{n+\frac{1}{2}\#} : \dot{\underline{e}}^{n+\frac{1}{2}\#}}. \quad (49) $$

Then, we propose the following discretization of (33)

$$ \frac{\underline{y}^{n+1} - \underline{y}^n}{\delta t} = \underline{v}^{n+\frac{1}{2}} = \frac{\underline{v}^{n+1} + \underline{v}^n}{2} \qquad (50a) $$

$$ \int_{\Omega^0} \rho_0 \frac{\underline{v}^{n+1} - \underline{v}^n}{\delta t} \cdot \underline{w} \, d\Omega + \int_{\Omega^0} \underline{\Sigma}^{n+\frac{1}{2}\#} : d_y \underline{e}^{n+\frac{1}{2}\#} \cdot \underline{w} \, d\Omega = \mathcal{P}_{\text{ext}}^{n+\frac{1}{2}}(\underline{w}), \forall \underline{w} \in \mathcal{V}_{\text{ad}} \quad (50b) $$

$$ \text{with } \underline{\Sigma}^{n+\frac{1}{2}\#} = \left. \frac{\partial \Psi}{\partial \underline{e}} \right|^{n+\frac{1}{2}\#} + \eta \frac{\underline{e}^{n+1} - \underline{e}^n}{\delta t} + \frac{T_{\text{fib}}^{n+\frac{1}{2}\#}}{(1 + 2\underline{\tau} \cdot \underline{e}^n \cdot \underline{\tau})^{\frac{1}{2}}} \underline{\tau} \otimes \underline{\tau} $$

$$ T_{\text{fib}}^{n+\frac{1}{2}\#} = \nu \frac{e_c^{n+1} - e_c^n}{\delta t} + \tilde{T}_c^{n+1} = E_s e_s^{n+\frac{1}{2}} \quad (50c) $$

$$ \text{with } \tilde{T}_c^{n+1} = \rho_{\text{surf}} \frac{\delta s}{d_a} \sum_{i=1}^{\ell} \frac{w_{1,i+1} - w_{1,i}}{\delta s} a_i^{n+1} $$

$$ \frac{a_i^{n+1} - a_i^n}{\delta t} = k_{+,i}(1 - a_i^{n+1}) - k_{-,i}a_i^{n+1} - \ell_{\text{fib}} \frac{e_c^{n+1} - e_c^n}{\delta t} \frac{a_i^{n+1} - a_{i-1}^{n+1}}{\delta s} \quad (50d) $$

Note here that (50c) and (50d) are defined at each quadrature point $\underline{x}_m$, albeit we omit – for the sake of brevity – this explicit dependence in the equations. If the 1D
elastic element is chosen nonlinear hyperelastic, the corresponding term in (50c) has to be treated as proposed for the 3D elastic element in (49).

To obtain a complete energy balance, we now proceed as in [2] by considering the mid-point velocity $\underline{v}^{n+\frac{1}{2}}$ as an admissible displacement field and recalling that $d_y \underline{e}^{n+\frac{1}{2}\#} \cdot \underline{v}^{n+\frac{1}{2}} = \frac{\underline{e}^{n+1}-\underline{e}^n}{\delta t}$. The balance associated with the hyperelastic contribution is handled by our choice in (49), and the viscous part directly gives a negative contribution, so that we have

$$ \frac{\mathcal{K}^{n+1} - \mathcal{K}^n}{\delta t} + \int_{\Omega^0} \frac{\Psi^{n+1} - \Psi^n}{\delta t} d\Omega \\ + \frac{1}{\delta t} \int_{\Omega^0} T_{\text{fib}}^{n+\frac{1}{2}\#} \frac{\underline{\tau} \cdot (\underline{\underline{e}}^{n+1} - \underline{\underline{e}}^n) \cdot \underline{\tau}}{(1 + 2\underline{\tau} \cdot \underline{\underline{e}}^n \cdot \underline{\tau})^{\frac{1}{2}}} d\Omega = \mathcal{P}_{\text{ext}}^{n+\frac{1}{2}} - \int_{\Omega^0} \eta \frac{|\underline{\underline{e}}^{n+1} - \underline{\underline{e}}^n|^2}{\delta t^2} d\Omega. $$

The function $x \mapsto (1+2x)^{\frac{1}{2}}$ is concave of derivative $(1+2x)^{-\frac{1}{2}}$, hence

$$ \frac{\underline{\tau} \cdot \underline{\underline{e}}^{n+1} \cdot \underline{\tau} - \underline{\tau} \cdot \underline{\underline{e}}^n \cdot \underline{\tau}}{(1 + 2\underline{\tau} \cdot \underline{\underline{e}}^n \cdot \underline{\tau})^{\frac{1}{2}}} \geq (1 + 2\underline{\tau} \cdot \underline{\underline{e}}^{n+1} \cdot \underline{\tau})^{\frac{1}{2}} - (1 + 2\underline{\tau} \cdot \underline{\underline{e}}^n \cdot \underline{\tau})^{\frac{1}{2}} = e_{\text{fib}}^{n+1} - e_{\text{fib}}^n. $$

Therefore we have, recalling that we have defined $e_s$ such that $e_{\text{fib}} = e_s + e_c$,

$$ \begin{aligned} \mathcal{P}_{\text{fib}} &= \frac{1}{\delta t} \int_{\Omega^0} T_{\text{fib}}^{n+\frac{1}{2}\#} \frac{\underline{\tau} \cdot (\underline{\underline{e}}^{n+1} - \underline{\underline{e}}^n) \cdot \underline{\tau}}{(1 + 2\underline{\tau} \cdot \underline{\underline{e}}^{n+1} \cdot \underline{\tau})^{\frac{1}{2}}} d\Omega \\ &\geq \int_{\Omega^0} T_{\text{fib}}^{n+\frac{1}{2}\#} \frac{e_{\text{fib}}^{n+1} - e_{\text{fib}}^n}{\delta t} d\Omega = \int_{\Omega^0} T_{\text{fib}}^{n+\frac{1}{2}\#} \left[ \frac{e_s^{n+1} - e_s^n}{\delta t} + \frac{e_c^{n+1} - e_c^n}{\delta t} \right] d\Omega \\ &\geq \int_{\Omega^0} E_s e_s^{n+\frac{1}{2}} \frac{e_s^{n+1} - e_s^n}{\delta t} d\Omega + \int_{\Omega^0} \left[ \nu \frac{(e_c^{n+1} - e_c^n)^2}{\delta t^2} + \tilde{T}_c^{n+1} \frac{e_c^{n+1} - e_c^n}{\delta t} \right] d\Omega. \end{aligned} $$

We now incorporate the Huxley’s75 discrete-time free energy balance (43) – and we could proceed identically with the other models using (48) – to get

$$ \begin{aligned} & \frac{\mathcal{K}^{n+1} - \mathcal{K}^n}{\delta t} + \int_{\Omega^0} \left[ \frac{\Psi^{n+1} - \Psi^n}{\delta t} + E_s \frac{|e_s^{n+1}|^2 - |e_s^n|^2}{2\delta t} + \rho_v \frac{\mathcal{F}^{n+1} - \mathcal{F}^n}{\delta t} \right] d\Omega \\ & \leq \mathcal{P}_{\text{ext}}^{n+\frac{1}{2}} - \int_{\Omega^0} \left[ \eta \frac{|\underline{\underline{e}}^{n+1} - \underline{\underline{e}}^n|^2}{\delta t^2} + \nu \frac{(e_c^{n+1} - e_c^n)^2}{\delta t^2} + T \rho_v \frac{S_{\text{prod}}^{n+1} - S_{\text{prod}}^n}{\delta t} \right] d\Omega \\ & \quad + \int_{\Omega^0} \rho_v \mu_T \frac{\delta s}{d_a} \sum_{i=1}^{\ell} J_{2,i}^{n+1} d\Omega, \quad (51) \end{aligned} $$

which is the discrete-time counterpart of the Clausius-Duhem relation (34), with here an inequality only due to numerical dissipation and consistent in $\delta t$. Note finally that in the case of models capturing the power stroke dynamics such as Piazzesi-Lombardi’95, time sub-iterations may be required. In this case it can be shown that the energy balance is preserved provided the active tension $\tilde{T}_c^{n+1}$ is redefined by weighing the intermediate states over all sub-iterations.
# 4 Numerical results and discussion

In this section, our goal is to illustrate the analysis of the discrete system presented in the previous section for the Huxley'57 model and the Piazzesi-Lombardi'95 model, which we chose as a representative of the multi-site, multi-state models. These illustrations serve several purposes. We first want to demonstrate that the thermodynamics identities established at the discrete level are satisfied in the numerical simulations. Then, we want to show that the ability to compute the thermodynamical balances numerically allows to gain additional insight into the physiology of muscle contraction. Additionally, for the Piazzesi-Lombardi'95 model, we compare our simulation results with that obtained in the original paper as a further validation of our approach.

## 4.1 Huxley'57 model

The choice of model parameters must satisfy the conditions (3) and the assumption that $w_1a$ and $k_-a$ tend to zero when $s$ tends to $s^-$ and $s^+$. We choose the energy levels and transition rates as follows

$$ \left\{ \begin{aligned} w_1(s) &= \kappa_w (s - s^*)^2 + \frac{\alpha_w}{s^+ - s} + \frac{\alpha_w}{s - s^-}, \\ k_1(s) &= k_{\max} \exp \left( -\lambda_1 \left[ (s - \bar{s})^8 + \frac{\alpha_{k_1}}{(s^+ - s)^2} + \frac{\alpha_{k_1}}{(s - s^-)^2} \right] \right), \\ k_2(s) &= k_{\text{mid}} - (k_{\text{mid}} - k_{\min}) \exp \left( -\lambda_2 (s - \bar{s})^8 \right) \\ &\quad + \frac{\alpha_{k_2}}{(s^+ - s)^2} + \frac{\alpha_{k_2}}{(s - s^-)^2}. \end{aligned} \right. \quad (52) $$

We choose here to prescribe $k_1$ and $k_2$ in addition to the energies $w_1$ and $w_0$. The reverse rates $k_{-1}$ and $k_{-2}$ are then derived from the detailed balance (1). The energy levels parametrization is shown in Figure 6. The transition rates are depicted in Figure 7. The models parameters used in the following simulations are presented in Table 1.

<figure><img src="image_7.png" /><figcaption>Fig. 6: Energy levels parametrization for the Huxley'57 model.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7: Transition rates parametrization for the Huxley'57 model.</figcaption></figure>

The asymptotic properties of the chosen transition rates and of the associated solutions can be easily obtained using the analytical solution (9) and the theorem of dominated convergence.
<table><caption>Table 1: Model parameters used in the simulations with the Huxley'57 model.</caption><thead><tr><th colspan="2">Model parameters</th></tr></thead><tbody><tr><td>μ<sub>T</sub></td><td>100 zJ</td></tr><tr><td>s<sup>+</sup></td><td>20 nm</td></tr><tr><td>s<sup>-</sup></td><td>-20 nm</td></tr><tr><td>s<sup>*</sup></td><td>0</td></tr><tr><td>tilde{s}</td><td>9 nm</td></tr><tr><td>tilde{s}</td><td>5 nm</td></tr><tr><td>κ<sub>w</sub></td><td>1.1 pN nm<sup>-1</sup></td></tr><tr><td>k<sub>max</sub></td><td>41.3 × 10<sup>-3</sup> ms<sup>-1</sup></td></tr><tr><td>k<sub>min</sub></td><td>10 × 10<sup>-3</sup> ms<sup>-1</sup></td></tr><tr><td>k<sub>mid</sub></td><td>30 k<sub>max</sub></td></tr><tr><td>λ<sub>1</sub></td><td>6.21 × 10<sup>-5</sup> nm<sup>-8</sup></td></tr><tr><td>λ<sub>2</sub></td><td>3 λ<sub>1</sub></td></tr><tr><td>α<sub>w</sub></td><td><sup>1</sup>&frasl;<sub>2</sub>κ<sub>w</sub>(tilde{s} - s)<sup>2</sup></td></tr><tr><td>α<sub>k<sub>1</sub></sub></td><td>k<sub>mid</sub>/λ<sub>1</sub></td></tr><tr><td>α<sub>k<sub>2</sub></sub></td><td>k<sub>mid</sub></td></tr></tbody></table>

We consider two simulation cases. First, we simulate the tension rise in isometric conditions ($v_c = 0$). Then, we compute the muscle response in contraction at constant shortening velocity ($v_c < 0$) starting from the isometric steady-state solution.

**Validation of the thermodynamical identities at discrete level** – We first want to verify that the discrete versions of the first principle (39) and the second principles (43) are satisfied numerically. To do so, we compute respectively the expressions

$$ \frac{U^{n+1} - U^n}{\delta t} - \left[ \frac{W^{n+1} - W^n}{\delta t} + \frac{Q^{n+1} - Q^n}{\delta t} + \frac{\varepsilon^{n+1} - \varepsilon^n}{\delta t} \right], $$

and

$$ \frac{\mathcal{F}^{n+1} - \mathcal{F}^n}{\delta t} - \frac{W^{n+1} - W^n}{\delta t} - \frac{\varepsilon^{n+1} - \varepsilon^n}{\delta t}. $$

The results for both simulation cases is presented in Figure 8. We notice that the first expression is ten orders of magnitude smaller than the individual terms that compose it (see Figures 9 and 10), showing that the first principle is satisfied at discrete level. The second expression is always negative showing the validity of the discrete second principle.

**Tension rise** – In our first illustration of the results obtained for the Huxley'57 model, we simulate the tension rise in isometric conditions ($v_c = 0$). We initialize all heads in the detached state and let the myosin heads evolve. Along the usual tension evolution, our scheme allows us to compute the thermodynamic fluxes associated with muscle contraction – see Figure 9. In the steady-state regime, the energy input remains positive and heat is dissipated. The force is sustained through the continuous cycling of the
<figure><img src="image_9.png" /><figcaption>Fig. 8: Validation of the discrete thermodynamical balances in two test cases for the Huxley'57 model.</figcaption></figure>

myosin heads in interaction with the actin filament. This process is fueled by the energy brought by ATP. We see here the active nature of muscle contraction. Force is produced when the muscle is supplied with energy. Naturally, as the velocity is zero no work is produced.

*Constant velocity contraction* – We now show a second illustrative example with a contraction at constant shortening velocity ($v_c < 0$) starting from the isometric steady-state solution. The simulation results are presented in Figure 10. After a transient phase, the system reaches a permanent regime in which the classical force-velocity curve is measured [9] (note that in the original experimental protocol force and not length is controlled). In this regime, we observe the energy mechano-transduction performed by the molecular motors: the energy input brought by ATP is for one part converted into work produced by the system ($W < 0$), and for the other part dissipated by entropy production.

## 4.2 Piazzesi-Lombardi'95 model

The Piazzesi-Lombardi'95 model reproduces the physiology of muscle contraction more precisely. In particular, it is able to capture the power stroke fast dynamics observed in length step experiments.

We simulate such an experiment starting from the isometric steady state with a length step of 8 nm. As in the experimental conditions, the length step is made by a ramp of duration 100 µs. Note that, here, the compliance of the myosin and actin filaments is neglected as in the original paper. We choose the energy levels as defined in [20]. We use modified transition rates to ensure that detachment rates diverge at
<figure><img src="image_10.png" /><figcaption>Fig. 9: Average tension and thermodynamic fluxes per myosin head in an isometric tension rise simulation with the Huxley'57 model. The isometric tension is denoted by $\tau_0$.</figcaption></figure>

infinity. The energy brought by ATP is set to 50 zJ following the model assumption that an ATP molecule can be used for the detachment of several myosin heads.

The results are presented in Figure 11. They match the results presented in the original paper [20], which shows the consistency of our approach with the original model, hence completed with thermodynamic balances.

## 5 Concluding remarks

Considering a large class of muscle contraction models based on actin-myosin interaction – i.e. the Huxley'57 model and various extensions thereof, including the Piazzesi-Lombardi'95 model – we have presented a mathematical setting in which solution properties can be established, including fundamental thermodynamic balances. Moreover, we have proposed a complete discretization strategy for which we were also able to obtain discrete versions of the thermodynamic balances and other properties. In addition, we have also shown how these models can be coupled with a macroscopic continuum mechanics formulation in such a way that these balances carry over to the macroscopic level, including for the discrete versions of the models. As muscle energetics are of major relevance in physiology, this is an important achievement, both from a fundamental and numerical point of view. This paves the way, indeed, for detailed
<figure><img src="image_11.png" /><figcaption>Fig. 10: Average tension and thermodynamic fluxes per myosin head for constant shortening velocity simulation with the Huxley'57 model. The isometric tension is denoted by $τ_0$.</figcaption></figure>

numerical studies of energy exchanges in various applications, such as with a complete realistic heart model.

## List of abbreviations

ATP: Adenosine triphosphate; PL95: Piazzesi-Lombardi'95; PVW: Principle of virtual work.

## Declarations

*Availability of data and material.* Nothing applicable.

*Competing interests.* The authors declare that they have no competing interests.

*Funding.* Nothing applicable.
<figure><img src="image_12.png" /><figcaption>Fig. 11: Tension evolution and thermodynamic fluxes per myosin head in a length step experiment simulation with Piazzesi-Lombardi'95 model. The isometric tension is denoted by $\tau_0$.</figcaption></figure>

*Authors' contributions.* FK was the main investigator in this work; he also performed the numerical implementation and simulations, and initiated the writing. DC was in a position of secondary supervision, and is responsible for some detailed technical aspects and writing. PM was in a position of primary supervision, and proposed initial ideas on topic and discretization strategy; he's also responsible for some detailed technical aspects and writing. In addition, DC and PM devised the proposed integration of muscle models in continuum mechanics and the associated discretization ingredients. All authors have read and approved the final manuscript.

*Acknowledgements.* Nothing applicable.

## References

1. M. Caruel, P. Moireau, and D. Chapelle. Stochastic modeling of chemical-mechanical coupling in striated muscles. *Biomechanics and Modeling in Mechanobiology*, pages 1–25, 2018.

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

3. M. Chipot, S. Hastings, and D. Kinderlehrer. Transport in a molecular motor system. *ESAIM: M2AN*, 38(6):1011–1034, December 2004.

4. E. Eisenberg and T.L. Hill. A cross-bridge model of muscle contraction. *Progress in Biophysics and Molecular Biology*, 33(1):55–82, 1978.
5. E. Eisenberg, T.L. Hill, and Y. Chen. Cross-bridge model of muscle contraction. Quanti-tative analysis. *Biophysical Journal*, 29(2):195-227, February 1980.

6. S.K. Godunov. A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics. *Matematicheskii Sbornik*, 89(3):271-306, 1959.

7. O. Gonzales. Exact energy and momentum conserving algorithm for general models in nonlinear elasticity. *Computer Methods in Applied Mechanics and Engineering*, 190(13-14):1763-1783, 2000.

8. P. Hauret and P. Le Tallec. Energy controlling time integration methods for nonlinear elastodynamics and low velocity impact. *Computer Methods in Applied Mechanics and Engineering*, 195:4890-4916, 2006.

9. A.V. Hill. The heat of shortening and the dynamic constants of muscle. *Proceedings of the Royal Society of London*, 1938.

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

11. T.L. Hill. *Free Energy Transduction And Biochemical Cycle Kinetics*. Dover, 2004.

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

13. A.F. Huxley and R.M. Simmons. Proposed mechanism of force generation in striated muscle. *Nature*, 1971.

14. F. Julicher, A. Ajdari, and J. Prost. Modeling molecular motors. *Reviews of Modern Physics*, 69(4):1269-1281, 1997.

15. P. Le Tallec and P. Hauret. Energy conservation in fluid structure interactions. *Numerical Methods for Scientific Computing - Variational Problems and Applications*, 2003.

16. R.W. Lynn and E.W. Taylor. Mechanism of adenosine triphosphate hydrolysis by acto-myosin. *Biochemistry*, 10(25):4617-4624, 1971.

17. L. Marcucci, T. Washio, and T. Yanagida. Including thermal fluctuations in actomyosin stable states increases the predicted force per motor and macroscopic efficiency in muscle modelling. *PLOS Computational Biology*, 12(9), 2016.

18. S. Mirrahimi and P.E. Souganidis. A homogenization approach for the motion of motor proteins. *Nonlinear Differential Equations and Applications*, 20(1):129-147, February 2012.

19. C.S. Peskin. *Mathematical Aspects of Heart Physiology*. Courant Institute of Mathematical Sciences - NYU, 1975.

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

21. R.D. Richtmyer and K.W. Morton. *Difference methods for initial-value problems*. Inter-science tracts in pure and applied mathematics. Interscience Publishers, 1967.

22. R. Sheshka and L. Truskinovsky. Power-stroke-driven actomyosin contractility. *Physical Review E*, 89(1), 2014.