François Kimmig

Dominique Chapelle

Philippe Moireau

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

Keywords: muscle contraction, sliding filaments, thermodynamically consistent time-discretization, Clausius-Duhem inequality Mathematics Subject Classification (2010) 74F25, 74H15, 65M12, 35Q79 and 92C45

Considering a large class of muscle contraction models accounting for actinmyosin 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.

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 [START_REF] Huxley | Muscle structures and theories of contraction[END_REF] modeling the attachmentdetachment 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 [START_REF] Eisenberg | A cross-bridge model of muscle contraction[END_REF][START_REF] Eisenberg | Cross-bridge model of muscle contraction. Quantitative analysis[END_REF][START_REF] Huxley | Proposed mechanism of force generation in striated muscle[END_REF][START_REF] Piazzesi | A cross-bridge model that is able to explain mechanical and energetic properties of shortening muscle[END_REF]. 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,[START_REF] Hill | Free Energy Transduction And Biochemical Cycle Kinetics[END_REF]. Note that these models are specific cases of molecular motors models without the natural diffusion introduced by the Fokker-Plank equation [START_REF] Julicher | Modeling molecular motors[END_REF][START_REF] Chipot | Transport in a molecular motor system[END_REF][START_REF] Mirrahimi | A homogenization approach for the motion of motor proteins[END_REF][START_REF] Chapelle | Energy-preserving muscle tissue model: formulation and compatible discretizations[END_REF]. In this paper, 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 [START_REF] Chapelle | Energy-preserving muscle tissue model: formulation and compatible discretizations[END_REF] 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 [START_REF] Peskin | Mathematical Aspects of Heart Physiology[END_REF][START_REF] Chipot | Transport in a molecular motor system[END_REF], or those mixing mechanical and chemical modeling elements, for instance [START_REF] Sheshka | Power-stroke-driven actomyosin contractility[END_REF][START_REF] Marcucci | Including thermal fluctuations in actomyosin stable states increases the predicted force per motor and macroscopic efficiency in muscle modelling[END_REF][START_REF] Caruel | Stochastic modeling of chemical-mechanical coupling in striated muscles[END_REF].

The outline of the paper is as follows. The first section presents the modeling ingredients 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 compatibility. Finally, the last section illustrates our results with numerical investigations.

1 Modeling of muscle contraction

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 [START_REF] Lymn | Mechanism of adenosine triphosphate hydrolysis by actomyosin[END_REF], 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. 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.

power stroke conformation recovery attachment detachment

Huxley'57 model

In his seminal work [START_REF] Huxley | Muscle structures and theories of contraction[END_REF], 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 da the distance between two consecutive attachment sites. The distance s thus lies in an interval of width da, 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 vc 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 µ 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

k i (s) k -i (s) = exp w k (s) -w j (s) k B T , (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 ∈ [s -, s + ]

     ∂a(s, t) ∂t = k 1 (s) + k -2 (s) d(s, t) -k 2 (s) + k -1 (s) a(s, t) -vc ∂a(s, t) ∂s ∂d(s, t) ∂t = k 2 (s) + k -1 (s) a(s, t) -k 1 (s) + k -2 (s) d(s, t) -vc ∂d(s, t) ∂s (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 [START_REF] Caruel | Stochastic modeling of chemical-mechanical coupling in striated muscles[END_REF], which implies that the energy of the attached level goes to infinity on the boundaries of the interval [s -, s + ]. In a nutshell, the parameter functions must satisfy

D A 1 2 v c (s = 0) s s + s - v c (s = 0) s s + s -
lim s→s - k 1 (s) = lim s→s + k 1 (s) = lim s→s - k -2 (s) = lim s→s + k -2 (s) = 0, lim s→s - 0 s k -1 + k 2 ds = lim s→s + s 0 k -1 + k 2 ds = +∞, (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→s - w 1 (s) = lim s→s + w 1 (s) = +∞.
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

τc(t) = 1 da s + s - a(s, t) ∂w 1 ∂s (s)ds. (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 [START_REF] Chapelle | Energy-preserving muscle tissue model: formulation and compatible discretizations[END_REF] or in Section 2.5.

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 [START_REF] Eisenberg | A cross-bridge model of muscle contraction[END_REF][START_REF] Eisenberg | Cross-bridge model of muscle contraction. Quantitative analysis[END_REF][START_REF] Huxley | Proposed mechanism of force generation in striated muscle[END_REF][START_REF] Piazzesi | A cross-bridge model that is able to explain mechanical and energetic properties of shortening muscle[END_REF]. 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 Ns chemical states {Xp} 1≤p≤Ns , that we can separate into two categories: attached states and detached states. The states are involved in Nr reactions between the states in the form

Xp i ki ------ k-i Xq i . (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 Ns vertices and 2Nr edges. A transition Xp i → Xq i is associated with the edge Ep iqi , whereas the reverse transition Xq i → Xp i is associated with the edge Eq ipi . The reaction presented in [START_REF] Eisenberg | Cross-bridge model of muscle contraction. Quantitative analysis[END_REF], which is bidirectional, is associated with the edges Ep iqi and Eq ipi . The subsets of vertices corresponding to attached and detached states are respectively denoted by Va and V d .

The ratio of heads in state Xq located at s at time t is denoted by xq(s, t). We define the chemical flux between states Xp i and Xq i through transition i by

Jp iqi (s, t) = k i (s)xp i (s, t) -k -i (s)xq i (s, t).

Note that we have

Jq ipi (s, t) = -Jp iqi (s, t) = k -i (s)xq i (s, t) -k i (s)xp i (s, t).
The system dynamics is then governed by

∂xq ∂t (s, t) = p|Epq∈G Jpq(s, t) -vc ∂xq ∂s (s, t), ∀q ∈ [[1, Ns]].
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

V = {A, D}, Va = {A}, V d = {D}, E = {E 1 AD , E 2 AD , E 1 DA , E 2 DA }.
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 + jda} j∈Z * . A myosin head can thus be detached in state q -with a probability xq(s) -or attached in state q at a distance s + jda of its rest position -with a probability xq(s + jda). We extend the definition of the ratio of heads in detached states by periodicity, i.e. if state q is a detached state xq(s + jda) = xq(s), ∀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 Ê and the remaining edges E. We define the associated complete directed graphs G d = (V, Ê) and Ga = (V, E). The system dynamics is governed by

         
(s, t), ∀Vq ∈ V d . (6)
The ratios of attached head xq for Vq ∈ Va are defined on R and must vanish at infinity. By contrast, the ratios of detached head xq for Vq ∈ 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 [START_REF] Piazzesi | A cross-bridge model that is able to explain mechanical and energetic properties of shortening muscle[END_REF] 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. An energy µ T is brought to the myosin head by ATP in the transitions 2 → 5 and 3 → 4.

A 1 V1 A 2 V2 A 3 V3 D 1 V4 D 2 V5 2
The graph G associated with this model is given by

V = {V 1 , V 2 , V 3 , V 4 , V 5 }, Va = {V 1 , V 2 , V 3 }, V d = {V 4 , V 5 }, E = 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 }.
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 wq the energy associated with the state of vertex q.

2 Model properties based on thermodynamics principles

(s, t) = k + (s) 1 -a(s, t) -k -(s)a(s, t) -vc ∂a ∂s (s, t), (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 vc > 0 and s + if vc < 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 vc ≥ 0, although the same result can be obtained similarly for vc < 0. As vc is assumed to be constant, the method of characteristic lines gives regular C

where we define h(s) = s 0 (k + (ξ) + k -(ξ)) dξ. Solving (8) along a characteristic line and pulling back the result to a(s, t) we obtain, defining t 0 = max(0, t -s-s - vc ),

a(s, t) = a 0 (s -vct) exp -1 vc s s-vct k + (ξ) + k -(ξ) dξ + t t0 k + (s -vct + vcτ ) exp -1 vc s s-vct+vcτ k + (ξ) + k -(ξ) dξ dτ. (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

s + 0 k -(ξ) dξ = ∞, we deduce that lim s→s + a(s, t) = 0 (10)
using the dominated convergence theorem for the second term of [START_REF] Hill | The heat of shortening and the dynamic constants of muscle[END_REF]. Likewise, the property lim s→s -a(s, t) = 0 is obtained for vc < 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) ∈ [0, 1], the property a(s, t) ∈ [0, 1] holds. Again, we rely on the solution obtained by the method of characteristic lines [START_REF] Hill | The heat of shortening and the dynamic constants of muscle[END_REF]. As the transition rates and the initial condition are positive, we find that a(s, t) ≥ 0. Then, noting that 1-a(s, t) is governed by an equation of the same form as [START_REF] Gonzales | Exact energy and momentum conserving algorithm for general models in nonlinear elasticity[END_REF] with the initial condition 1 -a 0 (s) ≥ 0, we similarly deduce that a(s, t) ≤ 1.

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

U(t) = 1 da s + s - w 1 (s)a(s, t) + w 0 d(s, t) ds. (11)
We then define the chemical fluxes

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).
We will henceforth make the natural assumption that the reaction rates are chosen in order for w 1 a, k 2 a and k -1 a to tend to zero when s tends to s -and s + , with the physical interpretation that no finite energy (w 1 a) is stored and no detachment flux (k 2 a 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

d dt U(t) = 1 da s + s - w 1 (s) ∂a(s, t) ∂t + w 0 ∂d(s, t) ∂t ds = 1 da s + s - (J 1 (s, t) -J 2 (s, t)) w 1 (s) -w 1 (s) vc ∂a(s, t) ∂s + (J 2 (s, t) -J 1 (s, t)) w 0 -w 0 vc ∂d(s, t) ∂s
ds.

Using integrations by parts for the transport terms, the boundary properties of the so-

lution -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

d dt U(t) = 1 da s + s - (w 1 (s) -w 0 ) J 1 (s, t) + (w 0 -(w 1 (s) + µ T )) J 2 (s, t) + vc a(s, t) ∂w 1 (s) ∂s ds + µ T J 2 (t),
where we denoted by

J 2 (t) = 1 da s + s - J 2 (s, t) ds
the mean net influx of ATP. Finally, we derive the following formulation of the first principle

U(t) = Ẇ(t) + Ė(t) + Q(t),
where

Ẇ(t) = vcτc(t), Ė(t) = µ T J 2 (t), Q(t) = 1 da s + s - (w 1 (s) -w 0 ) J 1 (s, t) + (w 0 -(w 1 (s) + µ T )) J 2 (s, t) ds, (12a) (12b) (12c)
with the active force τc defined in (4). The quantity Ẇ is the rate of work given to the system and Ė(t) = µ T J 2 (t) corresponds to the input flux in energy brought by ATP hydrolysis. The remaining term Q can be identified with a heat flux. In steady-state shortening (τc > 0 and vc < 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.

Second principle

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

S(t) = - k B da s + s - a(s, t) ln(a(s, t)) + d(s, t) ln(d(s, t)) 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) -T S(t),
which therefore corresponds to

F(t) = 1 da s + s - µ 1 (s, t)a(s, t) + µ 0 (s, t)d(s, t) ds. ( 13
)
where

µ 1 (s, t) = w 1 (s) + k B T ln(a(s, t)), µ 0 (s, t) = w 0 + k B T ln(d(s, t)),
are the chemical potentials. Computing the time derivative, we get

d dt F(t) = 1 da s + s - J 1 (s, t) -J 2 (s, t) -vc ∂ ∂s a(s, t) • w 1 (s) + k B T ln(a(s, t)) + k B T + J 2 (s, t) -J 1 (s, t) -vc ∂ ∂s d(s, t) • w 0 + k B T ln(d(s, t)) + k B T ds. ( 14
)
Using again integrations by part for the transport terms and the boundary properties of the solution, we obtain

d dt F(t) = τc(t)vc + µ T J 2 (t) + 1 da s + s - J 2 (s, t) µ 0 (s, t) -µ 1 (s, t) + µ T + J 1 (s, t) µ 1 (s, t) -µ 0 (s, t) ds. (15)
Independently, we also have

d dt F(t) = d dt U(t) -T d dt S(t). (16)
Combining equations ( 15) and ( 16) and the first principle, we can write

d dt S(t) = Q(t) T - 1 T 1 da s + s - J 2 (s, t) µ 0 (s, t) -µ 1 (s, t) + µ T + J 1 (s, t) µ 1 (s, t) -µ 0 (s, t) ds. ( 17
)
The second principle then reads

d dt S(t) = Q(t) T + Ṡprod (t), (18)
where we naturally associate the second term of ( 17) with the entropy production

Ṡ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

k 1 (s, t) k -1 (s, t) = exp w 0 -w 1 (s, t) k B T k 2 (s, t) k -2 (s, t) = exp w 1 (s, t) + µ T -w 0 k B T . (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 J 1+ J 1- = k 1 (s)d(s, t) k -1 (s)a(s, t) = exp w 0 -w 1 (s, t) k B T • exp k B T ln(d(s, t)) -ln(a(s, t)) k B T = exp µ 0 (s, t) -µ 1 (s, t) k B T . ( 20
)
As a consequence, we have two cases. If µ 1 (s, t) ≥ µ 0 (s, t), we find that

J 1+ J 1- (s, t) ≤ 1 ⇒ J 1 (s, t) = J 1+ (s, t) -J 1-(s, t) ≤ 0.
Conversely, µ 1 (s, t) ≤ µ 0 (s, t) implies that J 1 (s, t) ≥ 0. Proceeding in the same way for the second reaction, we finally have

J 1 (s, t) µ 1 (s, t) -µ 0 (s, t) ≤ 0, J 2 (s, t) µ 0 (s, t) -µ 1 (s, t) -µ T ≤ 0. ( 21
)
We thus obtain the conclusion that the entropy production term

Ṡprod (t) = - 1 T da s + s - J 2 (s, t) µ 0 (s, t) -µ 1 (s, t) + µ T + J 1 (s, t) µ 1 (s, t) -µ 0 (s, t) ds ≥ 0, ( 22
)
hence that the model is compatible with the second principle. We can summarize this property using (15) by the free energy balance where the mean fluxes are given by

d dt F(t) = τc(t)vc + µ T J 2 (t) -T Ṡprod (t) ≤ τc(t)vc + µ T J 2 (t). (23
Jpq(t) = 1 da s + s -j∈Z
Jpq(s + jda, t) ds,

and Wpq = Wpq -µ T δ pq={25,34} as an energy µ T is brought to the myosin head by ATP during the transitions 2 → 5 and 3 → 4. Note that the introduction of the modified energy increments W brings out the input energy fluxes µ T J 25 (t) + J 34 (t) .

The first principle then naturally reads

U (t) = Ẇ(t) + Ė(t) + Q(t), (25)
with

Ẇ(t) = vcτc(t), Ė(t) = µ T J 25 (t) + J 34 (t) , Q(t) = 1 da s + s -Epq∈G j∈Z
Jpq(s + jda, t) Wpq(s + jda) ds.

Second principle -Following the work done for the Huxley'57 model, we here define the entropy of the system as We here define the entropy production as

S(t) = - k B da s + s -Vq∈Va
Ṡprod (t) = - 1 T da s + s -Epq∈G j∈Z
Jpq(s + jda, t) Mpq(s + jda, t) ds,

where Mpq = µq -µp -µ T δ pq={34,25} . Then, combining the first principle (25) with the identity Ḟ = U -T Ṡ, we finally obtain the second principle

d dt S(t) = Q(t) T + Ṡprod (t).
Similarly as in the Huxley'57 model, the detailled 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

d dt F(t) = τc(t)vc + µ T J 25 (t) + J 34 (t) -T Ṡprod (t) ≤ τc(t)vc + µ T J 25 (t) + J 34 (t) . (26) R Ω 0 τ L Γ 0 N t N x
Fig. 4: Muscle fiber configuration.

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 x ∈ Ω 0 in the reference configuration. The displacement field associated with the deformation map is denoted by y. We denote by e the Green-Lagrange strain tensor, i.e.

e = 1 2 ∇ y + (∇ y) T + (∇ y) T • ∇ y ,
and the second Piola-Kirchhoff stress tensor is denoted by Σ. The fiber as shown in Figure 4 is subjected to a boundary force t N on a boundary Γ 0 N . The principle of virtual work (PVW) then reads: for any admissible virtual displacement field w ∈ V ad ,

Ω 0 ρ 0 ÿ • w dΩ + Ω 0 Σ : dye • w dΩ = Γ 0 N t N • w dΓ,
where the differential of the Green-Lagrange strain tensor with respect to the displacement field is given by dye

• w = 1 2 ∇ w + (∇ w) T + (∇ y) T • ∇ w + (∇ w) T • ∇ 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 Ψ and viscous pseudo-potential Ψ v -taken here as

Ψ v (ė) = η
2 tr(ė 2 ) -where η denotes a viscosity modulus -to simplify the presentation.

Following [START_REF] Chapelle | Energy-preserving muscle tissue model: formulation and compatible discretizations[END_REF], which extends the classical Hill-Maxwell scheme [START_REF] Hill | The heat of shortening and the dynamic constants of muscle[END_REF] to nonlinear behavior, we gather all the constitutive ingredients by defining an adequate rheological scheme -presented in Figure 5 We consider the following natural rheological rule for the parallel branch

3D parallel law : e = e p = e a , Σ = Σ p + Σ a , (27)
where e, e a , e p denote Green-Lagrange tensors (global, active and passive) and Σ, Σ a , Σ p second Piola-Kirchhoff stress tensors (global, active and passive). However, we will depart from [START_REF] Chapelle | Energy-preserving muscle tissue model: formulation and compatible discretizations[END_REF] 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 fib the total (local) extension of a fiber, i.e the ratio of length change over initial length

e fib = δ fib fib ,
where we will take for fib the length of a half-sarcomere at rest and δ fib the variation thereof. Then, the length change of the half-sarcomere can be decomposed into

δ fib = δ c + δ s,
where δ c and δ s respectively denote the contributions of the active (rigid filaments) and passive parts in this length change. Note that δ c then represents the relative displacement of the actin and myosin filaments considered above, indeed. We also introduce the corresponding dimensionless extension quantities ec = δ c/ fib and es = δ c/ fib , so that e fib = ec + es. Of course, the two components carry the same tension, which we denote by T 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 1D series law : e fib = ec + es, T fib = Tc = Ts.

(28) Note that the reason why such a simple additive rule holds for e 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 fib = 1 + 2τ • e • τ 1 2 . ( 29
)
Then, the tension T fib corresponds to a contribution in the first Piola-Kirchhoff stress tensor given by

T a = T fib F • τ F • τ ⊗ τ , (30)
where F = 1 + ∇ y is the classical deformation gradient tensor, as can be easily verified by computing the resulting traction in the fiber direction T a • τ . Hence, the associated contribution in the second Piola-Kirchhoff stress tensor reads

Σ a = F -1 • T a = T fib 1 + e fib τ ⊗ τ , (31)
since

F • τ = 1 + e fib .
Finally, the constitutive equations considered are

Σ p = ∂Ψ ∂e + ∂Ψ v ∂ ė , Ts = Eses, Tc = ν ėc + Ťc, (32)
where Ťc represents the aggregation of forces contributed by actin-myosin cross-bridges as described above, i.e. Ťc = ρ surf τc with ρ surf the number of myosin heads in a layer of thickness fib per unit of cross-section area. The series elastic element is here assumed to have a linear constitutive equation of elasticity modulus Es. Note that a nonlinear hyperelastic behavior could be considered, at the price of having to deal with the dimensionless extension es as an additional internal variable. Nevertheless, in physiological conditions the extension es remains small and a linear behavior is adequate [START_REF] Caruel | Stochastic modeling of chemical-mechanical coupling in striated muscles[END_REF]. As regards viscosity, we here incorporate a simple component of viscous modulus ν 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 ∂Ψ v ∂ ė = η ė.

We can now summarize the 3D equations as

                                   Ω 0 ρ 0 ÿ • w dΩ + Ω 0 Σ : dye • w dΩ = Γ 0 N t N • w dΓ, ∀w ∈ V ad with Σ = ∂Ψ ∂e + η ė + T fib (1 + 2τ • e • τ ) 1 2 τ ⊗ τ T fib = ν ėc + Ťc = Eses with Ťc(x, t) = ρ surf da s + s - a(x, s, t) ∂w 1 ∂s (s, t)ds ∂a ∂t (x, s, t) = k + (s) 1 -a(x, s, t) -k -(s)a(x, s, t) -fib ėc(x, t) ∂a ∂s (x, s, t) (33a) (33b) (33c)
where (33c) is based on the Huxley'57 model with sliding velocity vc = fib ėc. This velocity is independent of the microscopic variable s, which justifies our above study.

Note, however, the dependency of a(x, s, t) on 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 ẏ as an admissible virtual displacement field in (33a). We thus get

dK dt + Ω 0 Σ : ė dΩ = P ext ,
where K = 1 2 Ω 0 ρ 0 | ẏ| 2 dΩ stands for the kinetic energy and

P ext = Γ 0 N t N • ẏ dΓ is
the power of external forces. Then, we decompose

Ω 0 Σ : ė dΩ = d dt Ω 0 Ψ dΩ + Ω 0 η|ė| 2 dΩ + Ω 0 T fib τ • ė • τ (1 + 2τ • e • τ ) 1 2 dΩ,
where the first term is associated with the stored hyperelastic energy, the second term is a macroscopic viscous dissipation, and the last term -denoted P fib -is the power of internal forces in the sarcomere. Then, using the rheological rules we find

P fib = Ω 0 T fib τ • ė • τ (1 + 2τ • e • τ ) 1 2 dΩ = Ω 0 T fib ėfib dΩ = Ω 0
Ts ės + Tc ėc dΩ =

Ω 0 Eses ės + ν ė2 c + Ťc ėc dΩ = d dt Ω 0 Es 2 e 2 s dΩ (1)
+

Ω 0 ν ė2 c dΩ (2)
+

Ω 0 Ťc ėc dΩ (3)
, 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

d dt K + Ω 0 Ψ + Es 2 e 2 s + ρvF dΩ = P ext - Ω 0 η|ė| 2 + ν ė2 c dΩ - Ω 0 ρvT Ṡprod dΩ + Ω 0 ρv Ė dΩ, ( 34
)
where ρv = ρ surf / fib is the density of myosin head per unit volume in the reference configuration, and where we recall that F is the internal free energy of the bridges introduced in (13), Ṡprod is the entropy production term defined in (22) corresponding to energy dissipation associated with chemical transitions, and Ė = µ T J2 as defined in (12b) corresponds to the input flux in energy provided by ATP hydrolysis.

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.

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 δs and with the convention s 0 = s -and s = s + . We then choose an upwind implicit scheme, that for the sake of simplicity we only present for a positive sliding velocity vc, with a natural extension to negative sliding velocities vc 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 0 = 0 and with the natural condition a 0 i ∈ [0, 1]. The discretization scheme then reads

   a n+1 i -a n i δt = k +,i d n+1 i -k -,i a n+1 i -vc a n+1 i -a n+1 i-1 δs , ∀i ∈ [[1, ]] d n+1 i = 1 -a n+1 i , ∀i ∈ [[1, ]] (35)
with the definition

k +,i = k + (s -+ iδs), i ∈ [[1, ]], k -,i = k -(s -+ iδs), i ∈ [[1, -1]], k -, = 2k -, -1 . (36a) (36b) (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 n 0 = a n . Note that, with this choice, we do not strictly have a n 0 = a n = 0. This property is only satisfied approximately, or asymptotically when the spatial discretization length goes to zero.

Defining α = vcδt/δs and k i = k +,i + k -,i , the scheme can be written in a matrix form on the state vector

a n = a n 1 . . . a n T Da n+1 = a n + δtk + , where D =       1 + δtk 1 + α -α -α 1 + δtk 2 + α . . . . . . -α 1 + δtk + α       , k + =       k +,1 k +,2 . . . k +,       .

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 [START_REF] Richtmyer | Difference methods for initial-value problems[END_REF]).

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 ∀n, ∀i ∈ [[1, ]], a n i ∈ [0, 1].
To obtain this property, we need D to preserve the positivityi.e. for a ∈ R , Da ≥ 0 ⇒ a ≥ 0 (where we use the convention that a vector is positive if all its coefficients a i are positive). Let us then take a ∈ R such that Da ≥ 0. We have ∀i ∈

[[1, ]] 1 + δt k i + α a i -αa i-1 ≥ 0, (37)
with the boundary condition a 0 = a . Multiplying (37) by α l-1 for i = 1 and by

α -i i-1 j=1 1 + k j + α for i ∈ [[2, ]],
and summing, we obtain

i=2 α -i i j=1 1 + k j + α a i -α -i+1 i-1 j=1 1 + k j + α a i-1 + α l-1 1 + k 1 + α a 1 -α l a l ≥ 0 ⇔ j=1 1 + k j + α a + i=3 α -i+1 i-1 j=1 1 + k j + α a i-1 - i=2 α -i+1 i-1 j=1 1 + k j + α a i-1 + α l-1 1 + k 1 + α a 1 -α l a l ≥ 0.
Noting that the middle term is a telescoping series, we obtain as expected

j=1 1 + k j + α -α ≥0 a ≥ 0.
Then, recursively from (37), we get ∀i ∈ [[1, ]], a i ≥ 0, which shows that the matrix operator D preserves the positivity. Knowing that the initial condition and the transition rates are positive, we obtain ∀n ≥ 0, ∀i ∈ [ [1, ]], a n i ≥ 0. Writing the numerical scheme for the variable 1 -a n i from (35), we similarly obtain that ∀n ≥ 0

, ∀i ∈ [[1, ]], a n i ≤ 1.
Therefore, the proposed numerical scheme preserves the adequate positivity and boundedness.

Consistency -Let us here denote a a sufficiently regular solution of [START_REF] Gonzales | Exact energy and momentum conserving algorithm for general models in nonlinear elasticity[END_REF]. Note that this solution satisfies a Dirichlet boundary condition on one side of the simulation interval (i.e. a(s -, t) = 0) and that we showed in (10) that lim s→s + a(s, t) = 0. We denote by a the vector of the values of a at the spatial discretization points at time n δt, a n = a(δs, n δt) . . . a( δs, n δt)

T . We define, as usual, the convergence error by

e n = a n -a n .
We have De n+1 = e n -δt η n , where the consistency error η n is given by

η n = 1 δt Da n+1 -a n + δtk + .
Evaluating the continuous equation ( 7) at s = s -+ iδs and t = (n + 1)δt, we obtain

∂ā ∂t s -+ iδs, (n + 1)δt = k +,i -k i ān+1 i -vc ∂ā ∂s
s -+ iδs, (n + 1)δt , and this also holds for i = with the finite numerical value chosen for k , due to the fact that ā(s + ) = 0 and k -(s)ā(s) tends to zero in s + . We directly infer, ∀i ∈ [[1, ]],

η n i = ā s -+ iδs, (n + 1)δt -ā s -+ iδs, nδt δt - ∂ā ∂t s -+ iδs, (n + 1)δt +vc ā s -+ iδs, (n + 1)δt -ā s -+ (i -1)δs, (n + 1)δt δs - ∂ā ∂s s -+ iδs, (n + 1)δt = O(δs + δ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 [START_REF] Godunov | A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics[END_REF].

a n+1 i -a n i δt a n+1 i + vc a n+1 i -a n+1 i-1 δs a n+1 i = -k i (a n+1 i ) 2 .
Using the identity

-ab = 1 2 a -b 2 -1 2 a 2 -1 2 b 2 , we obtain 1 2δt (a n+1 i ) 2 -(a n i ) 2 + 1 2δt a n+1 i -a n i 2 = - vc 2δs (a n+1 i ) 2 -(a n+1 i-1 ) 2 - vc 2δs a n+1 i -a n+1 i-1 2 -k i (a n+1 i ) 2 .
Summing over i, and using the periodic boundary conditions (a n+1 ) 2 = (a n+1 0

) 2 , we find D -1 a n 2 ≤ a n 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 = D -1 n e 0 -δt n-1 k=0 D -1 n-k-1 η k ,
so that there exist C, T > 0 such that

e n 2 ≤ Cδt n-1 k=0
(δs + δt) = CT δs + δt .

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 [START_REF] Hill | Free Energy Transduction And Biochemical Cycle Kinetics[END_REF]. 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

w 1,i = w 1 (s -+ iδs), ∀i ∈ [[1, l -1]], w 1,l = 2w 1,l-1 ,
the energy is discretized as

U n = δs da i=1 w 1,i a n i + w 0 d n i .
Then, defining the fluxes as

J n 1,i = k 1,i d n i -k -1,i a n i , J n 2,i = k 2,i a n i -k -2,i d n i ,
we have ).

U n+1 -U n δt = δs da i=1 w 1,i a n+1 i -a n i δt + w 0 d n+1 i -d n i δt = δs da i=1 w 1,i J n+1 1,i -J n+1 2,i -vc a n+1 i -a n+1 i-1 δs + w 0 J n+1 2,i -J n+1 1,i + vc a n+1 i -a n+1 i-1 δs = δs da i=1 J n+1 1,i w 1,i -w 0 + J n+1 2,i w 0 -w 1,i -µ T + µ T J n+1 2,i - vcδs da -1 i=1 w 1,i -w 1,
(

) 38
With the definition of the discrete force

τ n+1 c = δs da i=1 w 1,i+1 -w 1,i δs a n+1 i ,
where we define w 1, +1 = w 1,1 , and using the periodicity of the solution, (38) becomes

U n+1 -U n δt = vcτ n+1 c + δs da i=1 J n+1 1,i w 1,i -w 0 + J n+1 2,i w 0 -w 1,i -µ T + µ T J n+1 2,i
.

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

U n+1 -U n δt = W n+1 -W n δt + Q n+1 -Q n δt + E n+1 -E n δt , (39)
with

W n+1 -W n δt = vcτ n+1 c , Q n+1 -Q n δt = δs da i=1 J n+1 1,i w 1,i -w 0 + J n+1 2,i w 0 -w 1,i -µ T , E n+1 -E n δt = µ T δs da i=1 J n+1 2,i .

Second principle

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

S n = -k B δs da i=1 a n i ln a n i + d n i ln d n i ,
and the free energy F n = U n -T S n , which can be rewritten as

F n = δs da i=1 µ n 1,i a n i + µ n 0,i d n i ,
by introducing the discrete chemical potentials

µ n 1,i = w 1,i + k B T ln a n i , µ n 0,i = w 0 + k B T ln d n i .
We then rewrite the previous calculation in a manner that closely follows the calculation in the continuous case. We have

F n+1 -F n δt = 1 δt δs da i=1 µ n+1 1,i a n+1 i -µ n 1,i a n i + µ n+1 0,i d n+1 i -µ n 0,i d n i = δs da i=1 µ n+1 1,i a n+1 i -a n i δt + a n i µ n+1 1,i -µ n 1,i δt + µ n+1 0,i d n+1 i -d n i δt + d n i µ n+1 0,i -µ n 0,i δt .
Hence,

F n+1 -F n δt = δs da i=1 µ n+1 1,i J n+1 1,i -J n+1 2,i -vc a n+1 i -a n+1 i-1 δs + a n i w 1,i + k B T ln(a n+1 i ) -w 1,i + k B T ln(a n i ) δt + µ n+1 0,i J n+1 2,i -J n+1 1,i + vc a n+1 i -a n+1 i-1 δs + d n i w 0 + k B T ln(d n+1 i ) -w 0 + k B T ln(d n i ) δt .
Developing the expression of the chemical potentials, performing Abel transformations and using the periodic boundary conditions, we obtain

F n+1 -F n δt = vcτ n+1 c + δs da i=1 J n+1 1,i µ n+1 1,i -µ n+1 0,i + J n+1 2,i µ n+1 0,i -µ n+1 1,i -µ T + µ T J n+1 2,i + vc da k B T i=1 a n+1 i ln(a n+1 i+1 ) -ln(a n+1 i ) + vc da k B T i=1 d n+1 i ln(d n+1 i+1 ) -ln(d n+1 i ) + k B T δs da i=1 a n i ln(a n+1 i ) -ln(a n i ) δt + k B T δs da i=1 d n i ln(d n+1 i ) -ln(d n i ) δt .
Since x → ln x is a concave function, we have ln a j -ln a i ≤ ln (a i ) a j -a i so that

F n+1 -F n δt ≤ vcτ n+1 c + δs da i=1 J n+1 1,i µ n+1 1,i -µ n+1 0,i + J n+1 2,i µ n+1 0,i -µ n+1 1,i -µ T + µ T J n+1 2,i + vc da k B T i=1 a n+1 i+1 -a n+1 i + vc da k B T i=1 d n+1 i+1 -d n+1 i + k B T δs da δt i=1 a n+1 i -a n i + k B T δs da δt i=1 d n+1 i -d n i .
Using the fact that the scheme imposes ∀n, ∀i ∈ [[1, ]], a n i + d n i = 1, the sums vanish two by two. We finally find

F n+1 -F n δt ≤ vcτ n+1 c δs da i=1 J n+1 1,i µ n+1 1,i -µ n+1 0,i + J n+1 2,i µ n+1 0,i -µ n+1 1,i -µ T + µ T J n+1 2,i
.

As a point-wise evaluation of the continuous expression ( 21), we have ∀i ∈ [[1, ]]

J n+1 1,i µ n+1 1,i -µ n+1 0,i ≤ 0, J n+1 2,i µ n+1 0,i -µ n+1 1,i -µ T ≤ 0.
Hence, in our case where vc > 0, we finally obtain

F n+1 -F n δt -vcτ n+1 c -µ T δs da i=1 J n+1 2,i ≤ 0. (40)
To relate the decrease of the free energy to the creation of entropy, we first note

F n+1 -F n δt -vcτ n+1 c -µ T δs da i=1 J n+1 2,i = F n+1 -F n δt - U n+1 -U n δt - Q n+1 -Q n δt = -T S n+1 -S n δt + Q n+1 -Q n δt . (41)
Comparing (41) with the formal expression of the second principle [START_REF] Mirrahimi | A homogenization approach for the motion of motor proteins[END_REF], we define the discrete entropy creation by

S n+1 prod -S n prod δt = S n+1 -S n δt - 1 T Q n+1 -Q n δt ≥ 0. (42)
Note that the entropy creation is formally given by

S n+1 prod -S n prod δt = - 1 T δs da i=1 J n+1 1,i µ n+1 1,i -µ n+1 0,i + J n+1 2,i µ n+1 0,i -µ n+1 1,i -µ T + vck B T da i=1 ln(d n+1 i ) -ln(a n+1 i ) a n+1 i -a n+1 i-1 + k B T δs da i=1 a n i ln(a n+1 i ) -ln(a n i ) δt + k B T δs da i=1 d n i ln(d n+1 i ) -ln(d n i ) δt .
and finally combining (40),( 41) and ( 42) we have the time-discrete counterpart of the free energy balance (23)

F n+1 -F n δt = vcτ n+1 c + µ T δs da i=1 J n+1 2,i -T S n+1 prod -S n prod δt ≤ vcτ n+1 c + µ T δs da i=1 J n+1 2,i . (43)

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 vc, we discretize the system (6) with the following implicit upwind numerical scheme

               x n+1 q,i+j -x n q,i+j δt = p|Epq∈G J n+1 pq,i+j -vc x n+1 q,i+j -x n+1 q,i+j -1 δs ,∀Vq ∈ Va, ∀i ∈ [[1, ]], j ∈ Z x n+1 q,i -x n q,i δt = p|Epq∈G j∈Z J n+1 pq,i+j -vc x n+1 q,i -x n+1 q,i-1 δs ,∀Vq ∈ V d , i ∈ [[1, ]] (44)
with the notation x n q,i+j = xq(s -+(i+j )δs, nδt) for Vq ∈ Va, x n q,i = xq(s -+iδs, nδt)

for Vq ∈ V d and J n pq,i+j = Jpq(s -+ (i + j )δs, nδt). Note that xq for Vq ∈ V d has periodic boundary conditions and thus x q,i+j = x q,i , ∀j ∈ Z, and we keep j ∈ Z, albeit in practice we will bound the attachment zone, introducing a boundary consistency error.

Mass conservation -Defining the total quantity of matter

m n i = Vq∈Va j∈Z x n q,i+j + Vq∈Vd x n q,i ,
-which has the periodicity m n 0 = m n -and using the scheme (44), we find the classical implicit transport equation

m n+1 i -m n i δt = -vc m n+1 i -m n+1 i-1
δs , hence, we retrieve, as in the Huxley'57 model, the conservation of matter

Vq∈Va j∈Z

x n q,i+j + Vq∈Vd x n q,i = 1 ∀n, ∀i ∈ [[1, ]]. (45)
First principle -Concerning the energy balance, the internal energy is now defined as

U n = δs da i=1 Vq∈Va j∈Z w q,i+j x n q,i+j + Vq∈Vd wqx n q,i .
and we find this time, with the notation Wpq,i+j = Wpq(s

-+ (i + j )δs), that U n+1 -U n δt = δs da i=1 Epq∈G j∈Z J n+1 pq,i+j Wpq,i+j + µ T i=1 j∈Z J 25,i+j + J 34,i+j - vcδs da i=1 Vq∈Va j∈Z w q,i+j x n+1 q,i+j -x n+1 q,i+j -1 δs - vcδs da i=1 Vq∈Vd wq x n+1 q,i -x n+1 q,i-1
δs .

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

τ n+1 c = δs da i=1 Vq∈Va j∈Z
x q,i+j w n+1 q,i+j +1 -w n+1 q,i+j δs ,

we obtain the discrete first principle

U n+1 -U n δt = W n+1 -W n δt + Q n+1 -Q n δt + E n+1 -E n δt , with W n+1 -W n δt = vcτ n+1 c , Q n+1 -Q n δt = δs da i=1 Epq∈G j∈Z J n+1 pq,i+j Wpq,i+j , E n+1 -E n δt = µ T δs da i=1 j∈Z
J 25,i+j + J 34,i+j .

Second principle -We now define the discrete entropy as

S n = -k B δs da i=1 Vq∈Va j∈Z
x n q,i+j ln x n q,i+j +

Vq∈Vd

x n q,i ln x n q,i , and the free energy F n = U n -T S n , that we rewrite as

F n = δs da i=1 Vq∈Va j∈Z
µ n q,i+j x n q,i+j + Vq∈Vd µ n q,i x n q,i , using discrete chemical potentials

µ n q,i+j = w q,i+j + k B T ln x n q,i+j , Vq ∈ Va, µ n p,i = wq + k B T ln x n q,i , Vq ∈ V .
The discrete time derivative of the free energy is

F n+1 -F n δt = δs da i=1 Vq∈Va j∈Z µ n+1 q,i+j
x n+1 q,i+j -x n q,i+j δt

+ x n q,i+j µ n+1 q,i+j -µ n q,i+j δt + i=1 Vq∈Vd µ n+1 q,i x n+1 q,i -x n q,i δt + x n q,i µ n+1 q,i -µ n q,i δt .
Using the notation M n pq,i+j = µ n q,i+j -µ n p,i+j -µ T δ pq={25,34} and the results of the first principle, we get

k B T ln x n+1 q,i+j x n+1 q,i+j -x n+1 q,i+j -1 δs -vc δs da i=1 Vq∈Vd k B T ln x n+1 q,i x n+1 q,i -x n+1 q,i-1 δs + k B T δs da i=1 Vq∈Va j∈Z x n q,i+j ln(x n+1 q,i+j ) -ln(x n q,i+j ) δt + k B T δs da i=1 Vq∈Vd x n q,i ln(x n+1 q,i ) -ln(x n q,i ) δt + vcτ n+1 c
.

Performing an Abel transformation on the transport terms, we obtain

F n+1 -F n δt = δs da i=1 Epq∈G j∈Z J n+1 pq,i+j M n+1 pq,i+j + µ T i=1 j∈Z J 25,i+j + J 34,i+j + vc da k B T i=1 Vq∈Va j∈Z x n+1
q,i+j ln(x n+1 q,i+j +1 ) -ln(x n+1 q,i+j )

+ vc da k B T i=1 Vq∈Vd x n+1
q,i ln(x n+1 q,i+1 ) -ln(x n+1 q,i )

+ k B T δs da i=1 Vq∈Va j∈Z
x n q,i+j ln(x n+1 q,i+j ) -ln(x n q,i+j )

δt

+ k B T δs da i=1 Vq∈Vd
x n q,i ln(x n+1 q,i ) -ln(x n q,i )

δt + vcτ n+1 c .
Using again that ln x j -ln x i ≤ ln (x i ) x j -x i , the detailed balance and the conservation of matter (45), we finally obtain

F n+1 -F n δt -vcτ n+1 c -µ T i=1 j∈Z J 25,i+j + J 34,i+j ≤ 0. (46)
As in Section (3.4) this property is equivalent to

S n+1 prod -S n prod δt = S n+1 -S n δt - 1 T Q n+1 -Q n δt ≥ 0. (47)
Moreover, the entropy production is in fact given by

S n+1 prod -S n prod δt = - 1 T δs da i=1 Epq∈G j∈Z J n+1 pq,i+j M n+1 pq,i+j - vck B T da i=1 Vq∈Va j∈Z ln(x n+1 q,i+j ) x n+1 q,i+j -x n+1 q,i+j -1 - vck B T da i=1 Vq∈Vd ln(x n+1 q,i ) x n+1 q,i -x n+1 q,i-1 + k B T δs da i=1 Vq∈Va j∈Z x n q,i+j ln(x n+1 q,i+j ) -ln(x n q,i+j ) δt + k B T δs da i=1 Vq∈Vd
x n q,i ln(x n+1 q,i ) -ln(x n q,i ) 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 [START_REF] Chapelle | Energy-preserving muscle tissue model: formulation and compatible discretizations[END_REF] -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

v n+ 1 2 = v n+1 +v n 2
for any variable v and the notation n + 1 2 will allow to indicate when we depart from this classical rule. First, as recommended in [START_REF] Gonzales | Exact energy and momentum conserving algorithm for general models in nonlinear elasticity[END_REF][START_REF] Tallec | Energy conservation in fluid structure interactions[END_REF][START_REF] Hauret | Energy controlling time integration methods for nonlinear elastodynamics and low velocity impact[END_REF][START_REF] Chapelle | Energy-preserving muscle tissue model: formulation and compatible discretizations[END_REF], we will consider the following non-standard -albeit classical -mid-point quantities

e n+ 1 2 = e(y n+ 1 2 ), ėn+ 1 2 = e n+1 -e n δt ,
dy e n+ 1 2 • w = 1 2 ∇ w + (∇ w) T + (∇ y n+ 1 2 ) T • ∇ w + (∇ w) T • ∇ y n+ 1 2 ,
and a passive hyperelastic stress law discretization that includes an energy correction term, namely ∂Ψ ∂e

n+ 1 2 = ∂Ψ ∂e (e n+ 1 2 ) + Ψ (e n+1 ) -Ψ (e n ) δt - ∂Ψ ∂e (e n+ 1 2 ) : ėn+ 1 2 ėn+ 1 2 ėn+ 1 2 : ėn+ 1 2 . (49)
Then, we propose the following discretization of (33)

                                                   y n+1 -y n δt = v n+ 1 2 = v n+1 + v n 2 Ω 0 ρ 0 v n+1 -v n δt • w dΩ + Ω 0 Σ n+ 1 2 : dye n+ 1 2 • w dΩ = P n+ 1 2 ext (w), ∀w ∈ V ad with Σ n+ 1 2 = ∂Ψ ∂e n+ 1 2 + η e n+1 -e n δt + T n+ 1 2 fib (1 + 2τ • e n • τ ) 1 2 τ ⊗ τ T n+ 1 2 fib = ν e n+1 c -e n c δt + Ť n+1 c = Ese n+ 1 2 s with Ť n+1 c = ρ surf δs da i=1 w 1,i+1 -w 1,i δs a n+1 i a n+1 i -a n i δt = k +,i (1 -a n+1 i ) -k -,i a n+1 i -fib e n+1 c -e n c δt a n+1 i -a n+1 i-1 δs (50a) (50b) (50c) (50d)
Note here that (50c) and (50d) are defined at each quadrature point 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 [START_REF] Chapelle | Energy-preserving muscle tissue model: formulation and compatible discretizations[END_REF] by considering the mid-point velocity v n+ 1 2 as an admissible displacement field and recalling that

dy e n+ 1 2 • v n+ 1 2 = e n+1 -e n δ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

K n+1 -K n δt + Ω 0 Ψ n+1 -Ψ n δt dΩ + 1 δt Ω 0 T n+ 1 2 fib τ • (e n+1 -e n ) • τ (1 + 2τ • e n • τ ) 1 2 dΩ = P n+ 1 2 ext - Ω 0 η |e n+1 -e n | 2 δt 2 dΩ. The function x → (1 + 2x) 1 2 is concave of derivative (1 + 2x) -1 2 , hence τ • e n+1 • τ -τ • e n • τ (1 + 2τ • e n • τ ) 1 2 ≥ (1 + 2τ • e n+1 • τ ) 1 2 -(1 + 2τ • e n • τ ) 1 2 = e n+1
fib -e n fib .

Therefore we have, recalling that we have defined es such that e fib = es + ec, We now incorporate the Huxley'57 discrete-time free energy balance (43) -and we could proceed identically with the other models using (48) -to get

P fib = 1 δt Ω 0 T n+ 1 2 fib τ • (e n+1 -e n ) • τ (1 + 2τ • e n+1 • τ ) 1 2 dΩ ≥ Ω 0 T n+ 1
K n+1 -K n δt + Ω 0 Ψ n+1 -Ψ n δt + Es |e n+1 s | 2 -|e n s | 2 2δt + ρv F n+1 -F n δt dΩ ≤ P n+ 1 2 ext - Ω 0 η |e n+1 -e n | 2 δt 2 + ν (e n+1 c -e n c ) 2 δt 2 + T ρv S n+1 prod -S n prod δt dΩ + Ω 0 ρvµ T δs da i=1 J n+1 2,i dΩ, (51)
which is the discrete-time counterpart of the Clausius-Duhem relation (34), with here an inequality only due to numerical dissipation and consistent in δ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 Ť n+1 c is redefined by weighing the intermediate states over all sub-iterations.

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.

Huxley'57 model

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

w 1 (s) = κw(s -s * ) 2 + αw s + -s + αw s -s -, k 1 (s) = kmax exp -λ 1 (s -s) 8 + α k1 (s + -s) 2 + α k1 (s -s -) 2 , k 2 (s) = k mid -k mid -k min exp -λ 2 (s -s) 8 + α k2 (s + -s) 2 + α k2 (s -s -) 2 .
(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. 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.

s - 0 s + 0 500 1,000 w 1 w 0 w 0 -µ T Energy (zJ)

Model parameters µ

T 100 zJ s + 20 nm s - -20 nm s * 0 s 9 nm s 5 nm κ w 1.1 pN nm -1 k max 41.3 × 10 -3 ms -1 k min 10 × 10 -3 ms -1 k mid 30 k max λ 1 6.21 × 10 -5 nm -8 λ 2 3 λ 1 α w 1 2 κ w (s -s) 2 α k1 k mid /λ 1 α k2 k mid
Table 1: Model parameters used in the simulations with the Huxley'57 model.

We consider two simulation cases. First, we simulate the tension rise in isometric conditions (vc = 0). Then, we compute the muscle response in contraction at constant shortening velocity (vc < 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

U n+1 -U n δt - W n+1 -W n δt + Q n+1 -Q n δt + E n+1 -E n δt , and
F n+1 -F n δt - W n+1 -W n δt - E n+1 -E n δt .
The results for both simulation cases is presented in Figure 8. We notice that the first expression is ten orders of magnitude smaller that the individual terms that compose it (see Figures 9 and10), 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 (vc = 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 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 (vc < 0) starting from the isometric steadystate 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 [START_REF] Hill | The heat of shortening and the dynamic constants of muscle[END_REF] (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.

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 [START_REF] Piazzesi | A cross-bridge model that is able to explain mechanical and energetic properties of shortening muscle[END_REF]. We use modified transition rates to ensure that detachment rates diverge at The results are presented in Figure 11. They match the results presented in the original paper [START_REF] Piazzesi | A cross-bridge model that is able to explain mechanical and energetic properties of shortening muscle[END_REF], which shows the consistency of our approach with the original model, hence completed with thermodynamic balances.

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

Declarations

Availability of data and material. Nothing applicable.

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

Funding. Nothing applicable.

Fig. 1 :

 1 Fig. 1: Lymn-Taylor cycle representation. Each stage of the cycle can be seen as a change of chemical state.

Fig. 2 :

 2 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.

Fig. 3 :

 3 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.

Fig. 5 :

 5 Fig.5: Fiber rheology combining a 1D active element (upper branch) and a 3D viscohyperelastic element (lower branch).

L 2 -

 2 Stability -We now analyze the L 2 -stability by justifying the 2 -stability of the operator D -1 . Multiplying (35) without the source term by a n+1 i , we have

F

6

 6 Discretization of the macroscopic model coupling

Fig. 6 :Fig. 7 :

 67 Fig. 6: Energy levels parametrization for the Huxley'57 model.

Fig. 8 :

 8 Fig. 8: Validation of the discrete thermodynamical balances in two test cases for the Huxley'57 model.

Fig. 9 :

 9 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 τ 0 .

Fig. 10 :

 10 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 .

Fig. 11 :

 11 Fig. 11: Tension evolution and thermodynamic per myosin head in a length step experiment simulation with Piazzesi-Lombardi'95 model. The isometric tension is denoted by τ 0 .

  

	∂xq ∂t	(s+jda, t) = p|Epq∈G Jpq(s+jda, t) -vc	∂xq ∂s	(s+jda, t), ∀Vq ∈ Va, j ∈ Z,
	∂xq ∂t	(s, t) = p|Epq∈Ga j∈Z	Jpq(s+jda, t) + p|Epq∈G d Jpq(s, t) -vc	∂xq ∂s

  )2.4 Extension to multi-state, multi-site modelsLet 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 ∀s ∈ [s -, s + ] d dt

						xq(s + jda, t) +	xq(s, t) = 0.	(24)
						Vq∈Va j∈Z	Vq∈Vd
	First principle -We define the energy as
					s +
		U(t) =	1 da	s -	Vq∈Va j∈Z	xq(s + jda, t)wq(s + jda) +
	The time derivative reads, after integrating by parts the transport term
						s +
		d dt	U(t) =	1 da
	s We then obtain
	d dt	U(t) = µ

Vq∈Vd xq(s, t)wq ds. -Epq∈G j∈Z Jpq(s + jda, t)Wpq(s + jda) ds + vcτc(t), with Wpq = wq -wp and the active force defined as τc(t) = 1 da s + s -Vq∈Va j∈Z xq(s + jda, t) ∂wq ∂s (s + jda) ds. T J 25 (t) + J 34 (t) + vcτc(t) + 1 da s + s -Epq∈G j∈Z Jpq(s + jda, t) Wpq(s + jda) ds,

Acknowledgements. Nothing applicable.

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.