Modeling actin-myosin interaction: beyond the Huxley–Hill framework

. Contractile force in muscle tissue is produced by myosin molecular motors that bind and pull on specific sites located on surrounding actin filaments. The classical framework to model this active system was set by the landmark works of A.F. Huxley and T.L. Hill. This framework is built on the central assumption that the relevant quantity for the model parametrization is the myosin head reference position. In this paper, we present an alternative formulation that allows to take into account the current position of the myosin head as the main model parameter. The actin-myosin system is described as a Markov process combining Langevin drift-diffusion and Poisson jumps dynamics. We show that the corresponding system of Stochastic Differential Equation is well-posed and derive its Partial Differential Equation analog in order to obtain the thermodynamic balance laws. We finally show that by applying standard elimination procedures, a modified version of the original Huxley-Hill framework can be obtained as a reduced version of our model. Theoretical results are supported by numerical simulations where the model outputs are compared to benchmark experimental data.


Muscle contraction
Muscle contraction is produced at the nanoscale by the concerted action of myosin molecular motors, converting metabolic energy into mechanical work.These motors are enzymatic proteins called myosins that catalyze the hydrolysis of ATP 1 in the presence of another protein, called f-actin.While interacting with f-actin, individual myosin motors can temporarily form a force generating unit called a cross-bridge.In the muscle tissue, at the microscale, myosins form thick filaments running in the longitudinal direction of the muscle fibers, see Figure 1.1 (a).F-actin also takes the form of filaments made out of polymerized g-actin (globular actin) which runs parallel to the myosin filaments.
The metabolic energy extracted from the hydrolysis of ATP by the motors is transformed into a mechanical working stroke that pulls on the actin filament.Both filaments are oriented to build a contractile unit as the myosin-actin interaction produces their relative sliding in opposite directions, see Figure 1.1 (a).At the mesoscale, antagonists contractile units are stacked in the transverse direction and form a sarcomere, whose length is of the order of 2 µm, see Figure 1.1 (b).Sarcomere themselves are arranged in series to form bundles of fibrils inside the muscle cells, see Figure 1.1 (c).Thanks to this highly organized structure, the metabolic activity of the actin-myosin molecular system results in the macroscopic contraction.At the core of the contraction mechanism is thus the interaction between a myosin protein and its neighboring actin monomers.This interaction takes the form of a cyclic succession of structural and chemical transitions during which a myosin binds actin to form a force producing a cross-bridge, and then disengage while hydrolyzing one ATP molecule, see [19] for more details on the molecular mechanism.The constant supply of fresh ATP molecules sustain a net flux around this cycle, which keeps the system out of equilibrium.

Modeling the actin-myosin interaction
The first physiologically relevant theory of actin-myosin interaction was formulated in 1957 by A.F. Huxley [20].The relative sliding of the actin and myosin filament was assumed to be powered by protrusions of the myosin filament, the heads of the myosin proteins.These protrusions are capable of attaching to specific binding sites positioned periodically on the actin filament, within fixed windows of length d, see Figure 1.2 (a).The force is produced by a large population of such heads, and it is assumed that the overall distribution of the binding site is uniform on the fixed window.
A fundamental hypothesis is that a given head can interact only with the nearest binding site, whose distance relative to the anchor point of the myosin head is denoted by s(t).A representative myosin head is then modeled as a linear elastic spring existing in two states: detached or attached to actin.The attachment and detachment events are modeled as a "chemical" jump process with rates that depend on the distance to the nearest binding site.
While the myosin is bound, the generated tension is proportional to the distance between the head of the protein, which is located at the binding site, and its anchor in the myosin filament.The dynamics of the system takes the form of a transport partial differential equation (PDE) for the population of attached and detached myosins, with the transport term accounting for the relative sliding velocity of the actin and myosin filament ( ẋc , see Figure 1.2 (a)) and sources and sinks terms accounting for the transition between the two states.
The parameters of the model to be calibrated are the energies of the two states and the transition rates.The aim of the model being to generate positive mechanical work, the basic calibration principle is to ensure that a given motor spends on average most of its time in a tensed configuration when it is attached, as represented in Figure 1.2 (a).separating the myosin anchor from its nearest binding site.The distance between two consecutive sites is denoted by d.After attachment, the formed cross-bridge is represented as a spring of length s(t).(b) The h-model considers the position of the detached myosin head X t relative to the same anchor point.The nearest binding site is now defined relative to the position of the head, and parametrized by the variable h t .This modeling framework has been formalized and generalized by T.L. Hill and co-workers [11,12,16,17,18].The generalization consisted in formulating the dynamics equations for an arbitrary number of chemical-mechanical states and arbitrary transitions between them, and in formalizing the fundamental mathematical requirements that ensure the compatibility of the model with conservation laws, in particular the thermodynamic principles.
Since then, plethora of such chemical-mechanical models have been proposed to account for a constantly increasing body of experimental data.They enable, for instance, interactions with multiple binding sites, see [3,8,21,28,29,31,35,39,40,41,42].Another refinement consists in using generalized non-linear energy potentials for the elastic spring representation of the myosin head [24,30,34].
In this paper, we start by presenting the Huxley-Hill framework and some of its fundamental properties, see Section 2.
Second, we propose in Section 3 an alternative approach by reconsidering the initial assumptions of the Huxley-Hill framework.In particular, we propose a parametrization of the model by the current position of the myosin head, which is closer to the actual molecular configuration.Our new model, called h-model, is based on the following assumptions see Figure 1.2 (b).
(1) A myosin head is parametrized by the position X t of its head with respect to its anchor point on the myosin filament at time t.
(2) The position of the nearest actin binding site is not defined with respect to the anchor point anymore, but instead with respect to the position of the head X t .This relative position is denoted by h t .
(3) A given myosin head can bind to the nearest actin site only, but the fact that the head can extend beyond the inter-site distance d allows interactions with a site that is not necessarily the closest to the anchor point.
Situations where several myosin heads compete to bind to the same actin site can theoretically occur.Unlike a recent model [37], our framework does not account for these situations, which will not occur in practice when a physiological calibration is used.
We formulate the h-model as a continuous-time stochastic Markov process involving a continuous overdamped Langevin dynamics for the position of the head, and a discrete Poisson dynamics for the state (attached or detached).
Throughout the paper, we use two approaches to describe the system: a stochastic approach based on Stochastic Differential Equations (SDE in short) and a deterministic approach based on Partial Differential Equations (PDE in short).The latter is used to draw a parallel between our proposed framework and the Huxley-Hill model.In particular, we show that the Huxley-Hill model can in fact be viewed as an analog of a reduced version of the h-model (the h-reduced model, see Section 4).
For all models (Huxley-Hill, h-model and h-reduced model), we study the compatibility with the thermodynamic principles.
Finally, numerical illustrations of the h-model and its reduced version are presented in Section 5, where we compare the steady state properties of the models with benchmark experimental data that characterize the mechanical performance of cardiac cells.

Notations and conventions
• In this article, elements of the h-model (fraction of myosin heads and model parameters) are denoted by non-overlined quantities.Overlined quantities correspond to reduced quantities associated with the h-reduced model.
• Quantities indexed by t are stochastic variables.
• Energy functions are denoted by w.
• Poisson random measures are denoted by N (dt, du).
• The one-dimensional torus of length d is denoted by T d .The related canonical projection is π : R → T d .We will use a coordinate map • Jump rates for the Huxley-Hill model are denoted by k.
• Jump rates for the h-model are denoted by K.
• Jump rates for the h-reduced model are denoted by k.
• The notation a ∨ b is used for the maximum of two real numbers a and b.
• The space of continuous functions R n → R is denoted by C(R n , R), and the space of differentiable functions with continuous partial derivatives by C 1 (R n , R).
• The space of bounded differentiable functions with bounded continuous partial derivatives is denoted by The Huxley-Hill model considers the dynamics of the attachment state α t ∈ {0, 1}, which takes the form of a jump process with two jumps associated with the rates f and g.(b) As the myosin head hydrolyzes one ATP molecule per cycle, the energies of the detached state at the beginning and at the end of the cycle differ by the amount of metabolic energy µ T , harvested through this reaction.Hence, the cycle is better represented by two "chemical" reactions, each having a forward rate (solid lines) and a reverse rate (dashed lines).(c) Each state α is associated with an energy function w α , and the metabolic input takes the form of a downward energy shift µ T .The bold line in (b) illustrates a possible evolution of a myosin head in this energy landscape as the position of the nearest actin site slides at a negative velocity (contraction).Each segment of this evolution, labelled with Q, W and E, can be associated with an element of the first principle of thermodynamics, namely heat exchange, work, and energy influx, respectively.

The framework
The paradigmatic Huxley-Hill model [20] sets a 1D framework to describe the dynamics of a population of independent myosin motors located on a representative rigid myosin filament, interacting with a single rigid actin filament, see Figure 1.1 (a).A representative myosin motor is shown in Figure 1.2 (a).
The binding sites are regularly positioned on the actin filament, the distance between consecutive sites being denoted by d.At all times, the myosin head has the possibility to bind only to the nearest actin site, whose position with respect to the anchor point is parametrized by the variable s.
Hence, the binding to the nearest site occurs necessarily within a restricted fixed window of size d, here assumed to be (−d/2, +d/2) for simplicity.The current actin site is thus at most d/2 far from the base point located at 0. This site is assumed to move with a known velocity ẋc (t) at time t ( ẋc < 0 for contraction).Therefore, during contraction, a site enters the window (s = d/2), at the time the previous one leaves (s = −d/2).

Model formulation
A discrete variable α t ∈ {0, 1} is introduced to characterize the attachment state of the myosin head.This variable takes the value α t = 1 if the head is attached, and the value α t = 0 if the head is detached, see Figure 2.1(a).The dynamics of α t is a jump process with two jumps, whose rates are continuous non-negative functions of s ∈ (−d/2, +d/2).We define • the attachment transition 0 → 1 with rate f (s); • the detachment transition 1 → 0 with rate g(s).
The Huxley-Hill model considers the population of myosin heads having their closest attachment site located at s ∈ (−d/2, +d/2).The function (t, s) → P 1 (t, s), represents the fraction of this population that is attached at time t.Similarly, the function (t, s) → P 0 (t, s) := 1 − P 1 (s, t) represents the fraction of this population that is detached at time t.The fraction of attached myosin heads verifies the following PDE (2.1a) for any initial condition s → P ini 1 (s) ∈ [0, 1] that satisfies P ini 1 (±d/2) = 0.The latter boundary condition is the direct consequence of the hypothesis that a given myosin head can interact only with its nearest binding site, which must be less than d/2 away from the head anchor point.Attachment to sites located at s = ±d/2 is not allowed.To satisfy this condition, the following constraint is imposed on the attachment rates and detachment rate g must ensure that the heads detach when nearing the boundaries.Explicit conditions on the rate functions will be exposed in the well-posedness analysis below.

Well-posedness of the Huxley-Hill model
To the best of our knowledge, the conditions required for the well-posedness of (2.1) had never been made explicit in the muscle contraction modeling literature before.We establish the following well-posedness result.
Proposition 2.1 (Well-posedness).Given T > 0, we consider the following hypotheses. • The main interest of this result is to make explicit the mathematical conditions on rates f and g that must be imposed for (2.1) to be valid.The following proof is based on the method of characteristics: looking at some specific curves t → s(t) with ṡ(t) = ẋc (t), it is possible to obtain the explicit expression (2.5) below.An analogous proof would work if ẋc (t) < 0 for every t, instead of ẋc (t) > 0.
Proof.For any s in (−d/2, +d/2) and t > 0, we define the backward characteristic curve rooted at (t, s) by Since ẋc (v) > 0, u → S(u, t, s) is an increasing function.Furthermore, we define the hitting time As explained above, the Huxley-Hill model postulates that the attached myosin head behaves as a spring, transmitting force between the two filaments.To link the dynamic equation (2.1a) to the force produced by the myosin head, it is necessary to model the energetic properties of a myosin head.Let s → w 1 (s) be a C 1 function on (−d/2, +d/2), that models the energy of the myosin head in the α = 1 (attached) configuration.For instance, if the attached head is modeled as an elastic spring, the energy w 1 (s) is a quadratic function of s, see Figure 2.1 (c).Similarly, w 0 stands for the energy of the myosin head in the α = 0 (detached) configuration.In the classical formulation [20], w 0 is a constant, but the following results remain valid if w 0 is a function of s.
Myosin motors are ATPases, capable of converting the energy provided by ATP hydrolysis into mechanical work.The concentration of ATP molecules is considered buffered in the cell, so that the system can be maintained out of equilibrium.The chemical potential µ T of ATP is introduced as a constant parameter, representing this energy supply.
The recruitment of an ATP molecule occurs when the myosin head detaches from the actin site.This supply modifies the energy of the detached state, so that the energy of the head is lower at the end of the cycle than at the beginning, see Figure 2.1 (b,c).Hence, the detachment transition rate g cannot be considered, using the language of solution chemistry, as the reverse of the attachment transition f .Following the formalism proposed by T.L. Hill, the cycle is thus better represented as a sequence (α = 0) → (α = 1) → (α = 0), where the two transitions can be viewed as a two different "reactions", each being characterized by a forward and a reverse jump (see Figure 2.1 (b)): We define accordingly • a jump 0 → 1 with rate k rev 1→0 (s) called reverse detachment.
The above functions can be paired to recover the total transition rates, f and g: Since the rates are non-negative functions of s, the condition (2.2) translates into We will now show how this decomposition can be used to ensure the compatibility of the model with the first and second principles of thermodynamics.Definition 2.2 (First principle).Formally, the thermodynamic properties of the Huxley-Hill model are computed as averages over the whole population of heads.We start by defining the quantities associated with the first principles of thermodynamics: internal energy, mechanical power, energy influx and heat dissipation.
• The internal energy of the system is defined as the average value of 1 α=0 w 0 + 1 α=1 w 1 (s) over the spatial window, i.e.
where the coefficient 1/d translates the assumed uniform distribution of the binding sites in the interval (−d/2, +d/2).
• The instantaneous power of external efforts is defined as where the force generated by the attached motors is defined as the averaged gradient of w 1 over the window: ∂ s w 1 (s)P 1 (t, s) ds.
• Each actin-myosin interaction cycle consumes one molecule of ATP [27], which is recruited upon detachment of the myosin head from actin.Each ATP molecule providing the energy µ T , the energy supply can be defined as • The heat dissipation is defined as This definition is tailored to guarantee the compatibility of the model with the first principle of which appears as a direct consequence of (2.1a).The first principle is illustrated on Figure 2.1(c), where the thick line represents a possible evolution of a myosin head as its site slides ( ẋc < 0).The two vertical jumps correspond, from right to left, to the attachment and the detachment, respectively.The corresponding (net) heat exchanges are the first and the second terms of (2.7), respectively.The work is associated with the part of the evolution along the energy landscape w 1 , and the energy input is represented by the vertical shift of the energy level w 0 .

Definition 2.3 (Second principle)
. We now define the thermodynamic quantities associated with the second principle of thermodynamics.
• For α ∈ {0, 1}, the chemical potential of state α at time t is defined as ), where T is the absolute temperature and k B is the Boltzmann constant.
• The free energy of the system is defined as the average value of chemical potentials, i.e.
• The system is said to satisfy the second principle of thermodynamics if d dt If this holds, the entropy production is defined as the non-negative quantity There is no a priori guarantee that the model will satisfy (2.9) for any choice of the rate functions k 0→1 , k rev 0→1 , k 1→0 , and k rev 1→0 .A reasonable assumption is based on the classical thermochemistry rule that the quotient of direct and reverse rates of a chemical reaction scales as the exponential of the energy difference between the states.In this case, the rates are said to verify the so-called detailed balance condition.This law is for instance directly fulfilled if individual rates follow the Arrhenius law, which stipulates that the inverses of rates scale as the exponential of the energy barrier between the states.
We recover that the conditions of detailed balance between the forward and reverse rates are sufficient to ensure that the model satisfies the second principle (2.9) [18,Chapter 5].

Lemma 2.4. Assuming the detailed balance relations
) the system satisfies the second principle of thermodynamics according to Definition 2.9.
In [20], the model was directly formulated with the total rates f and g (see Equation (2.6)).Imposing a detailed balance relation directly on f and g would also make the model compatible with the second principle.However, for the system to be maintained out of equilibrium and to produce mechanical work, the energy landscapes before attachment and after detachment must differ, see Figure 2.1.Therefore, the detailed balance between f and g should be broken, see also [23].Splitting these two jumps into four jumps, that verify pairwise detailed balance conditions, is a way to fulfill this condition, while ensuring the compatibility of the model with the second principle.Furthermore, these relations constrain the choice of w 0 and w 1 , according to the well-posedness conditions of Proposition 2.1.
Proof.Using that P 0 + P 1 = 1, the time-derivative of Adapting the computations in the proof of Lemma B.1 to (2.1a) yields  [4,11,12,35].Our framework could be extended to these settings, if suitable detailed balance conditions are provided.

Towards a stochastic formalism
In this section, we reformulate the Huxley-Hill model (2.1) using a stochastic description of the jump process.One advantage of this formulation compared to the PDE formulation is to enable Monte-Carlo simulations, which can be numerically more efficient than PDE discretization.
Another advantage is that the stochastic formulation gives access to the behavior of individual myosin heads (instead of a whole population), which can be interpreted in terms of a nanoscale molecular mechanism.This formulation is also more appropriate to model in vitro experimental setups where only a few motors are involved [34].
Eventually, the stochastic formulation sets the framework for presenting the h-model, see Section 3. As a preliminary, we consider the following toy model, which is the analogous of the Huxley-Hill PDE system extended to the full real line: for some Lipschitz continuous functions f , g : R → R + .
Lemma 2.5.Given a C 1 initial condition P ini 1 , Equation (2.11) has a unique solution P in Proof.The coefficients f and g being Lipschitz, this is a classical consequence of the method of characteristics, which moreover provides the explicit formula using the same notations as (2.5).The well-posedness issue is easier to handle here than for (2.1), because there is no explosion on rates anymore, nor boundary condition.□ We define the Markov process (α t ) t≥0 that models the state of a given myosin head.This process shall depend on the location s(t) of the closest actin site, which can now be anywhere in the infinite window (−∞, +∞).
The position of this site is represented by the C 1 b curve t → s(t) with velocity ẋc (t), which is assumed to be a given parameter.To build the desired process, we consider a filtered probability space (Ω, F, (F t ) t≥0 , P).On this space, we define two independent adapted Poisson random measures on R + × R + , N 0→1 (dt, du) and N 1→0 (dt, du), whose intensity measure dt ⊗ du is the Lebesgue measure on R + × R + .We then consider the right-continuous with left limits (càdlàg) Markov process (α t ) t≥0 defined by (2.12) For a practical approach, this formulation is further explained and illustrated in the presentation of the algorithm used to simulate it, see Appendix A. A complete presentation can be found in [1,22,44].
Proof.The validity of (2.12) is a classical result because of the affine bounds that we assumed on the coefficients f and g: the construction of the process can be found in e.g.[ We could now proceed to the restriction of the toy model (2.11) to a finite window (−d/2, +d/2), alike the original Huxley-Hill PDE (2.1a).Since the functions f , g should satisfy the growth assumptions required by the well-posedness result 2.1, the rates f = f and g = g appearing in (2.11) would not have linear growth anymore.In this context, it is not clear that the SDE (2.12) is still well-posed: the solution could explode in finite time.
Instead of performing a SDE well-posedness analysis in this non-bounded setting, which would be interesting in itself, we propose an alternative stochastic formulation of the system, for which a well-posedness result can be obtained.

The h-model
The h-model is built as an alternative to the Huxley-Hill model that avoids its well-posedness difficulties in the SDE formulation.Apart from this mathematical consideration, the motivation behind the h-model is to set a framework where the position of the nearest actin site is defined with respect to the position X t of the tip of the myosin head instead of the myosin anchor, see Figure 1.2.The dynamics of the new variable X t is driven by the gradient of the detached state potential, and by thermal fluctuations.
Furthermore, the h-model does not assume that the maximal extension of the head is half the inter-site distance, this fact being experimentally uncertain.Instead, the tip X t can explore the whole real line, and in this respect the h-model radically differs from the Huxley-Hill model.
The h-model parametrization is set in Section 3.1.In Sections 3.2 and 3.3 we present its SDE and PDE formulations respectively, before deriving the associated thermodynamic balances in Section 3.4.Finally, in Section 3.5, we show that the Huxley-Hill model can be viewed as a heuristic limit of the h-model.

Model parametrisation
As in Section 2.5, the h-model uses a SDE to describe the attachment variable α t of the myosin head at time t, depending on the location of the closest actin site.
A first stochastic equation describes the dynamics of the attachment variable α t as in (2.12).A second equation, that is coupled to the first one, is introduced to describe the dynamics of the position X t ∈ R of the myosin head, as it was done in the model developed in [6].In addition to jump terms, this latter dynamics is driven by a Brownian motion which models thermal fluctuations; it is thus a diffusion with jumps.As in the Huxley-Hill framework, the head can only bind to a single actin site at each time, but this site is now the closest one to the head tip X t , and not the closest one to the myosin anchor anymore.
The actin sites are separated by a distance d, so that either there is a unique actin site within the window [X t − d/2, X t + d/2] and this site is in the interior (X t − d/2, X t + d/2), or there are exactly two actin sites in [X t − d/2, X t + d/2], and they are located at X t ± d/2.To exclude this singular case, we will require, as previously, that attachment rates vanish at s = X t ± d/2.
To take this constraint into account, we consider the distance between the myosin head tip and its closest action site as a continuous variable h t in a one-dimensional torus of width d centered at X t , see Figure 1

SDE formulation
As in Section 2, let (Ω, F, (F t ) t≥0 , P) be a filtered probability space.We require Ω to be large enough to provide four independent Poisson measures, N 0→1 (dt, du) and N rev 1→0 (dt, du) on R + × R + with intensity measure dt ⊗ du, together with N 1→0 (dt, du, dh) and N rev 0→1 (dt, du, dh) on R + × R + × T with intensity measure dt ⊗ du ⊗ dh.The measure dh on the torus is inherited from the Lebesgue measure on [−d/2, +d/2].Let (B t ) t≥0 be a Brownian motion independent of the Poisson measures.We then search for a {0, 1} × R × T d -valued Markov process (α t , X t , h t ) t≥0 , as the càdlàg strong solution of the SDE system When α t = 0, X t is driven by the Brownian noise and the gradient of w 0 , w 0 : R → R being a given C 2 b function which models the energy of an unattached myosin tip.Similarly, we introduce a C 2 b function w 1 : R → R, which models the energy of an attached myosin tip.As previously, T is the ambient temperature and k B denotes the Boltzmann constant, η being a constant damping coefficient.
The h-model considers two possible 0 → 1 jumps for α t (using the Poisson measures N 0→1 and N rev 1→0 ), together with two 1 → 0 jumps (using the Poisson measures N 1→0 and N rev 0→1 ).During an attachment event 0 → 1, (X t − , h t − ) jumps to the new value (X t , π(0)) such that X t = X t − + h t − .From Equation (3.1), X t − + h t − is the location of the closest actin site at t − .These attachment jumps (x, h) → (x + h ⋆ , π(0)) for (X t , h t ) occur at rates K 0→1 (x, h) and K rev 1→0 (x, h).To guarantee that the attachment occurs only on the nearest site, we ensure that X t cannot jump over a distance larger than d/2 by imposing 3) Similarly, during a detachment event 1 → 0, the couple (X t − , h t − ) with h ⋆ t − = 0 jumps to a new value (X t , h t ), with X t + h ⋆ t = X t − .The closest actin site at time t is the one located at X t − , to which the tip was attached before the jump.A jump 1 → 0 must select a new value for h t , hence the need for a h-component in the jump measures N 1→0 and N rev 0→1 .These detachment jumps (x, π(0)) → (x − h ⋆ , h) for (X t , h t ) occur at rates K 1→0 (x, h) and K rev 0→1 (x, h).As in Section 2.4, we want to link the rates of direct and reverse jumps using some detailed balance conditions like (2.10a)-(2.10b).We thus require that We emphasize that in (3.4) the variable h stands for the value before the jump, while it stands for the value after the jump in (3.5).From the above relations, the detachment rates K 1→0 and K rev 1→0 can be recovered from the knowledge of the energy functions w 0 and w 1 , together with the attachment rates K 0→1 and K rev 1→0 .We require K 0→1 and K rev 1→0 to be bounded and Lipschitz continuous.Since w 0 and w 1 are C 2 b , K 1→0 and K rev 1→0 are bounded and Lipschitz continuous too.As in Section 2.4, the reverse jumps with rates K rev 0→1 and K rev 1→0 will allow the model to be thermodynamically consistent, see Lemma 3.4 below.Lemma 3.1 (Well-posedness).Under the above assumptions on rates, the SDE system (3.2) has a strong càdlàg solution (α t , X t , h t ) t≥0 which enjoys pathwise uniqueness.
The proof of this well-posedness result can be done using classical fixed-point methods for jump-diffusion SDEs.The standard method is given in [22, Chapter 4.9] in the case of jumpdiffusion processes in R d .In our case the rates K 0→1 , K rev 0→1 , K 1→0 , K rev 1→0 , and drift terms ∂ x w 0 , ẋc , are bounded and Lipschitz continuous (in the present situation, the linear growth of jump rates is actually sufficient), hence the fixed-point method can be applied for (α t , X t , h t ) to each term of our large SDE system.

PDE related system
The above section has built a stochastic version of the h-model, building the process (α t , X t , h t ) t≥0 as the solution of the SDE system (3.2).We now write the PDE counterpart of (3.2), describing the law of (α t , X t , h t ) at time t as the solution of a PDE system.From classical results on jump-diffusion processes (using Ito's formula and martingale compensation see e.g.[22,Chapter 4.9] or [1, Theorem 4.4.7]), the law of (α t , X t , h t ) is expected to solve the PDE system below in the (weak) measure sense.However, for ease of reading, the following PDE system is written in strong from.
The first three lines of (3.6a) correspond to the diffusion and drift terms which drive dX t and dh t in (3.2).The cross differential term −2η −1 k B T ∂ x ∂ h p(t, 0, x, h) comes from the fact that X t and h t are driven by the same Brownian noise.Equation (3.6b) can be seen as a transport equation, for which (t, x, h) → p(t, 0, x, h) is a source term.If (3.6) admitted a unique solution in the sense of distributions, then one could establish a link with (3.2) analogous to Lemma 2.6.However, such a uniqueness result is not guaranteed here.If (t, x, h) → p(t, 0, x, h) is given and known to be C 1 in the t and x variables, then the method of characteristics yields using the notations of (2.5).Inserting (3.7) into (3.6a)gives a large PDE on (t, x, h) → p(t, 0, x, h).
For the sake of completeness, we mention a well-posedness result for this latter equation.

Lemma 3.2 (PDE well-posedness).
The PDE on (t, x, h) → p(t, 0, x, h) obtained by inserting (3.7) into (3.6a)has a unique solution in Proof.The change of variable s = π(x) + h in T d allows us to write this large PDE as where Q is a non local operator which depends on p in a bounded linear way, and such that Q(t, x, s, p) remains a linear in p Lipschitz perturbation of the heat equation.Therefore, the fixedpoint method in the proof of [33,Theorem 3.13] can be applied here.An alternative variational approach is given in [32, Chapter 2].
□ The remaining of this paper aims to derive qualitative properties of the h-model, in a less formal way.In particular, we will assume that the solution of the PDE system (3.6) is a bounded smooth function.

Thermodynamic balances
Thermodynamic properties can now be studied as in Section 2.4.The analogous of Definition 2.2 is the following definition.

Definition 3.3 (First principle)
. We here define the quantities associated with the first principles of thermodynamics: internal energy, mechanical power, energy influx and heat dissipation.
• The internal energy of the system reads • As previously, the instantaneous power of external efforts reads • The energy supplied by ATP consumption reads which depends on the net detachment flux • The heat dissipation is defined as Comparing to (2.7), there are two new terms involving w 0 : these terms account for the overdamped Langevin dynamics on x.This definition is tailored to guarantee the compatibility of the model with the first principle of thermodynamics which appears as a direct consequence of (3.6).For the second principle, Definition 2.3 applies formally to the h-model.
As in Lemma 2.4, we now show that the detailed balance conditions (3.4)-(3.5)are sufficient to ensure that the model satisfies the second principle (2.9).

Lemma 3.4. Assuming the detailed balance relations (3.4) and (3.5), the h-model satisfies the second principle of thermodynamics according to Definition 2.3.
Proof.Chemical potentials being defined as µ(t, 0, x, h) = w 0 (x) + k B T ln p(t, 0, x, h) and µ(t, 1, x) = w 1 (x) + k B T ln p(t, 1, x), the free energy reads As previously, the time derivative of the free energy can be put in the form and straightforward computations show that this requires As in the proof of Lemma B.1, Equations (3.4) and (3.5) then imply that Ṡprod (t) is nonnegative.□

The Huxley-Hill model as a heuristic limit of the h-model
From the h-model, we expect to recover the Huxley-Hill system (2.1a) in some suitable limit.
For that purpose, we re-introduce the framework in which the Huxley-Hill is defined; namely, assume formally that X t almost surely stays within the finite window (−d/2, +d/2).This should be possible by allowing rates to explode, with no guarantee that, in this non-Lipschitz situation, (3.2)-(3.6)would still be meaningful.
To recover a function defined on [−d/2, +d/2], we identify h to the coordinate h ⋆ ∈ [−d/2, +d/2), requiring periodic boundary conditions in h for p.To alleviate notations, h ⋆ will be kept written as h in the sequel.The Huxley-Hill variable s of the Huxley-Hill model can then be introduced as s = x + h.We perform the change of variable x = s − h, and define p(t, 0, s, h) = p(t, 0, s − h, h) and p(t, 1, s) = p(t, 1, s), for s in [−d/2, +d/2].From (3.6), we obtain (3.9) We now define the total rates

Model reduction
To further investigate the relation between the h-model and the seminal Huxley-Hill model, we put forward a formal model reduction that provides an extension of (2.1a) in which the assumption of a restricted window for binding is relaxed.Therefore, in this derived h-reduced model, the myosin head position evolves on the full real line and the head can potentially bind any actin site.However, this model retains the hypothesis that only one possible binding site is available for attachment at each time, which implies that the myosin head position discontinuity allowed during an attachment event is at most d/2.Like the Huxley-Hill model, the obtained PDE only involves the variables (t, s), s being the location of the closest actin site to the myosin head.
The derivation of the h-reduced model is performed in Section 4.1, and its form is compared to the original Huxley-Hill formulation (2.1a) in Section 4.2.Section 4.3 then proposes a formal writing of the h-reduced model using SDEs.The thermodynamic properties of the h-reduced model are finally presented in Section 4.4.In the following, we do not perform the model derivation in a completely rigorous mathematical manner.Instead, we aim at a heuristic result, whose validity can then be tested numerically (see Section 5).

Derivation of the h-reduced model
Starting from the PDE h-model (3.6), the model derivation is performed in three steps: (1) Change of variable.We first parametrize h ∈ T d by its coordinate h ⋆ , supplemented with appropriate boundary conditions.We then replace the location x of the myosin head by the position s = x + h of its closest actin.
(2) Integration over h.We integrate Equation (3.6) over h: this yields an exact but nonclosed form equation on the ratio of myosin head in state (0, s).
(3) Adiabatic elimination of h.To obtain a closed-form formulation, we assume that the dynamics of h t is much faster than the dynamics of s t .The validity of this latter assumption is investigated in Section 5.
The functions s → p(t, 0, s), p(t, 1, s) thus quantify the ratio of myosin heads in state (0, s) and (1, s), respectively.Integrating (3.9) in h yields the exact relation which is non-closed form as it still depends on p(t, 0, s, h).Thanks to the periodicity relation (4.3), the first three lines become:

) Adiabatic elimination of h.
To obtain a closed-form of (4.4) in the variable s, we now eliminate the variable h by assuming that p(t, 0, s, h) = p(t, 0, s)p th 0 (h|s) ( where p th 0 stands for the thermalized density of h given s.This density is assumed to take the canonical Boltzmann-Gibbs form with the partition function With this assumption, we consider that the dynamics of h t is much faster than the dynamics of s t , which is driven by the model parameter ẋc , and the attachment-detachment dynamics.The variable h t is thus assumed to be a fast variable.As in Section 3.5, we define the total transition rates Finally, the thermal equilibrium assumption (4.6) cancels the diffusion terms in (4.4), giving the following system (4.9) This system meets the objectives of the model reduction and defines the h-reduced model.It generalizes the Huxley-Hill system (2.1a) to the full real line.

Comparison between the h-reduced
appearing on the first line of (4.9).This term is reminiscent of the fact a myosin head does not discriminate sites that are equally far from it.
When the actin sites translate with positive speed ẋc (t) > 0, the population of myosin heads located at distance d/2 to these sites changes.Sites located at s are no longer the closest ones for myosin heads that were located at x = s − d/2, hence the loss of the related fraction of heads and the term − ẋc (t)p(t, 0, s, d/2).Simultaneously, myosin heads located at x = s+d/2 now have their nearest actin site located at s, which corresponds to the addition of the related fraction of heads, hence the term ẋc (t)p(t, 0, s, −d/2).After application of the periodicity relation (4.3) and the adiabatic elimination of the variable h, we obtain the correction term (4.10).Similarly, the same result is obtained in the case where ẋc (t) < 0. Remark 4.1 (Loss of the periodicity relation).Under the closure assumption (4.6), the periodicity relation (4.3) is not verified anymore.Indeed, the reduced system (4.9) is only an approximation of the exact marginal relation (4.4).The periodicity relation could be preserved through the approximation process, if the periodicity relation satisfied the closure assumption, i.e. if we had the relation which is not satisfied.The loss of this periodicity relation will impinge on the thermodynamic balances in Section 4.4.

Relation between p(t, α, s) and P α (t, s).
The densities p(t, α, s) of (4.9) differ from the global ratios P 0 (t, s) and P 1 (t, s) of detached and attached heads knowing s appearing in (2.1a).Indeed, the variable s t of the h-reduced model is reminiscent of the myosin head position X t at time t and is thereby a random variable, while s(t) is a deterministic parameter in (2.12).The ratios P 0 and P 1 can however be computed from p(t, 0, s) and p(t, 1, s) using where p(t, s) := p(t, 0, s) + p(t, 1, s).Summing the equations in (4.9), we obtain Notice that the function (t, s) → p(t, s) cannot be computed from the latter equation, because it is not in closed-form.Therefore, it is not possible to directly write an evolution PDE for P 0 (t, s) and P 1 (t, s).However, P 0 and P 1 can be computed from (4.12), after solving (4.9).

Stochastic differential equation representation of the reduced model
The correction term (4.10) corresponds to a mass transfer from s to s − d.As for (2.11), it is possible to build a SDE representation for (4.9).To simplify the presentation, we impose that ẋc (t) > 0 for every t.We then consider a filtered probability space (Ω, F, (F t ) t≥0 , P), and three independent adapted Poisson random measures on R + × R + , N 0→1 (dt, du), N 1→0 (dt, du) and N corr (dt, du), with intensity measures dt ⊗ du.We then define a Markov process (α t , s t ) t≥0 as the càdlàg strong solution of the SDE system The additional jump on s t stands for the aforementioned mass transfer from s to s − d.Similarly to Lemma 2.6, one can prove a well-posedness result for (4.13) and identify the law of (α t , s t ) as the solution of (4.9).

Reduced thermodynamic balances
In this section, we look at the thermodynamic properties of the reduced model (4.9).Our purpose is to give an insight on the way the thermodynamic balances presented in Sections 2.4 and 3.4 for the Huxley-Hill model and the h-model, respectively, are modified when considering the h-reduced model.Therefore, we do not use the structure of Sections 2.4 and 3.4, and we rather focus on the changes at the level of each balance.

Detailed balance conditions
The energy functions involved in the detailed balance conditions (3.4)- (3.5) for the h-model are w 0 and w 1 , w 0 being a function of x = s − h ⋆ .In the h-reduced model, h ⋆ has been thermalized and s is considered as the leading coordinate.It is thus natural to replace w 0 by the free energy w 0 (s) := −k B T ln Z 0 (s), the attached-state energy being w 1 (s) := w 1 (s).The compatibility of w 0 and w 1 with the jump rates of the h-reduced model must now be checked.Introducing the change of variable x = s−h ⋆ into the detailed balance relations (3.4)-(3.5),we obtain To recover the rates of the h-reduced model, we integrate over h.Using (4.6)-(4.7),let us define the reduced jump rates as This yields and similarly s).These relations give the desired detailed balance conditions with respect to the reduced energy w 0 , which must be satisfied by the rates of the h-reduced model.

Energy balance
We now compute the time-derivative of the internal energy U(t) := R w 0 (s)p(t, 0, s) + w 1 (s)p(t, 1, s) ds, to recover the first principle of thermodynamics.The terms are the same as in (2.8), up to an extra contribution due to the correction term, which reads: The first principle eventually reads is the power required for the displacement of s t , the closest actin site to X t , in the energy landscapes w 1 and w 0 .In the detached state, s t is no longuer subject to a classical springlike effort with potential energy w 0 , because of the correction term in (4.9).This explains the quadratic correction to the energy w 0 (s) appearing in (4.14): less power is needed, because the correction term prevents the head from escaping too far from the origin.This latter behavior is reminiscent of the dynamics of the head X t , which was pulled back to the origin like a spring when escaping too far.If the periodicity relation (4.11) were true for the closure relation (4.6), the second line in the computation of d dt U(t) would be zero.This latter energy contribution is thus an artificial extra term, due to the approximation error.

Free energy balance
As in Section 2.4, chemical potentials can be defined as µ(t, 0, s) := w 0 (s) + k B T ln p(t, 0, s) and µ(t, 1, s) := w 1 (s) + k B T ln p(t, 1, s).The free energy rewrites  (4.6) introduces an additional contribution, which would be zero in the exact setting.This contribution is not always negative, hence the second principle of thermodynamics is not satisfied anymore.To satisfy the second principle, one could add a reverse jump which compensates the jump induced by the correction term.However, the physical interpretation of such a reverse jump would be unclear.

Numerical illustrations
We now present a calibration of our newly introduced hierarchy of models and show its capability to reproduce key indicators of the cardiac muscle contraction physiology.Our target for the calibration is two categories of fundamental physiological modes of contraction: the isometric contraction (contraction mode in which the relative sliding of the actin and myosin filaments is prevented in the experimental setup) and the steady-state contraction at constant shortening velocity.These two modes of contraction have historically been used to characterize the behavior of muscles [15] and are key features of the physiology at the heart level in the isovolumetric contraction phase and ejection phase, respectively.More precisely, we aim to reproduce the ratio of attached myosin heads and the force per attached myosin head during isometric contraction and the force-velocity curve in steady-state conditions (where the force τ c is normalized by the isometric force τ c,0 , see Figure 5.2).

Model calibration
To calibrate the models, the following parameters need to be provided: the actin inter-site distance d, the attached and detached state energy potential w 1 (x) and w 0 (x) and the attachment rates K 0→1 (x, h) and K rev 1→0 (x, h).Note that the detachment rates K rev 0→1 (x, h) and K 1→0 (x, h) result from the attachment rates and the energy potentials through the detailed balance equations (3.4) and (3.5).
We describe here the main principles that have guided the calibration process.A detailed presentation of the chosen model parameters is given in Table D.1 and illustrated in Figure 5.1 for the h-model and in Figure 5.3 for the h-reduced model.
First, a structural choice must be made.The actin filament is composed of a double helix of actin monomers that is covered by tropomyosin proteins whose length corresponds to the helix periodicity and which drive the muscle activation by uncovering the actin sites [26].This gives rise to two paradigms that are both used in the literature for the choice of the inter-site distance: the actin filament helix periodicity (∼40 nm) [6,12,24,28,39,41] (this paradigm is selected in this paper) or the size of the actin monomer (∼5.5 nm) [4,35].
Then, the chosen isometric indicators (the ratio of attached heads ňa and the force per attached head τc ) bring constraints on the transitions rates and the attached energy potential.These indicators for the h-model are given by ňa where p ∞ (1, x) = lim t→∞ p(t, 1, x) in isometric conditions ( ẋc = 0).The ratio of attached heads depends solely on p ∞ .It is mostly driven by the attachment-detachment process, which is characterized by the associated rates.The force per attached head depends on the attached myosin head nonlinear stiffness and on the stretch of the attached myosin heads, and therefore on transition rates, which determine the favored extension levels of attached myosin heads.Finally, the shape of the force-velocity curve will add constraints on the transition rates and in particular the detachment rate.Indeed, the variation of the force with the sliding velocity in a steady-state contraction results from the competition between the variation of the force per attached head and the ratio of attached heads.The former results from variations of the attached myosin extensions due to filament sliding and the attachment-detachment process, which allows counter-productive myosin heads to detach and re-attach in a position where they now contribute positively to the force.The latter solely depends on the attachment and detachment rates.
The attachment rate is set so that myosin heads can attach when x lies in an interval of positive values and when they are close to their nearest actin site, i.e. when h belongs to a limited interval in absolute value, see Figure 5.1 (b).
Following recent works [6,24,34], we choose a regularized piecewise quadratic potential for the attached potential w 1 , see Figure 5.1 (a).This shape is reminiscent of the power-stroke  [5].The force per attached head is derived from the data given in [36], considering that the average macroscopic force scale linearly with the ratio of attached heads.

Isometric indicators
Experiments h-model h-reduced model Ratio of attached heads 0.15 0.151 0.145 Force per attached head 6.14 pN 6.39 pN 6.19 pN conformational change, which operates at very short timescale [21] and which is associated with a double-well energy landscape [6,24,34].Our lump representation is thus valid as long as this conformational change can be considered fast compared to other processes, in particular the attachment-detachment processes.This is reasonably the case in classical force-velocity experiments.However, without this additional dynamics, the model cannot be used to reproduce the fast transient dynamics, see [24].
Fewer experimental data are available to guide the calibration of the detached potential w 0 .We thus choose the same form as the attached potential (piecewise quadratic, reminiscent of the underlying double-well energy landscape of the nano-level myosin conformations) and set the energy wells such that the detached myosin heads are attracted towards the region where they can attach.Even if uncertain, this representation seems compatible with recent structural and molecular dynamics studies [2].The detachment rate is then set to impose that myosins detach when they operate in compression and to control the ratio of attached heads at the physiological sliding velocities.In the calibration process, it is the attachment rate K rev 1→0 that is prescribed, the detachment rate K 1→0 being obtained through the detailed balance relation (3.5), see Figure 5.1 (c).

h-model
We first present the results obtained with the h-model.The model equations are simulated with the scheme presented in Appendix A. The isometric indicators (obtained with ẋc = 0) are summarized in Table 5.1.The force-velocity curve is presented in Figure 5.2.As it has already been noticed in the literature [14,23,24], the force slightly increases for shortening velocities that are small in absolute value compared to the isometric force level [24].This results from the fact that, in this range of velocities, the ratio of attached head increases and compensates the decrease of the force per attached head [24].This range of shortening velocities is not physiologically relevant in a cardiac muscle contraction, and instabilities are observed experimentally [9,10,13,43].The shape of the force-velocity curve in this region has thus not been a focus in the calibration process.
The physiological indicators simulated with the h-model show a good agreement with the experiment data.Since most of the model parameters are functions, various calibrations would give similar results.The chosen rates are regularized piecewise constant functions to keep the number of scalar parameter minimal.With this choice of function we found that the leeway for the parameter values is relatively restrained showing that the model is well constrained.

h-reduced model
To assess the impact of the modeling assumptions used to derive the h-reduced model (see Section 4.1), we compare the model predictions against that of the h-model.The h-reduced 0 0.2 0.4 0.6 0. where p ∞ (1, s) = lim t→∞ p t (1, s) in isometric conditions ( ẋc = 0).Table 5.1 presents a comparison of the physiological indicators.Both models provide predictions of the isometric indicators that are close to the experimental data.
The h-reduced model is based on the assumption that the variable h of the h-model can be considered at thermal equilibrium.Its dynamics is thus assumed to be much faster than that of the actin filament, which corresponds to the dynamics of the variable s.To challenge this assumption, we simulate the force velocity curves of both models considering two values for the internal viscosity η.At higher internal viscosity values, the dynamics of h will be slowed down and we expect the fundamental assumption to fault.At a given velocity, the force developed by the h-model should be reduced because of a reduction of the myosin head cycling rate.Indeed, due to the higher viscosity, myosin heads need a longer time to reach region where they can re-attach while the h-reduced model does not depend on the viscosity value.
The results are presented in Figure 5.4.For the physiological value of the internal viscosity (that we denote by η ref ), the force velocity curve for the h-model and the h-reduced model display little differences.This validates the h-reduced model fundamental assumption when the timescale of interest is the one involved in the force-velocity curve, which is mainly driven by the attachment-detachment dynamics.However, we see that with an internal viscosity 50 times highers, the two force velocity curves differ showing that the fundamental assumption would not be applicable.Further investigation on the validity of the h-reduced model would be necessary if the model targetted finer timescales such as the timescale of the power stroke.fundamental assumption that the variable h can be abiabatically eliminated is not valid anymore, and the two models differ, the h-model producing less force than the h-reduced model.

Conclusion
In this work, we introduced a novel modeling framework of the actin-myosin interaction, which represents an alternative to the seminal Huxley-Hill framework.A jump-diffusion stochastic model is proposed.It describes the myosin configuration by its attachment state, the extension of the myosin neck and the distance to the nearest actin site relatively to the current myosin tip position.Unlike Huxley-Hill formulations, no restriction of the myosin extension in the detached state is imposed, while still assuming that the head can attach to only one site at a time.
Interpolation functions

Figure 1 . 1 .
Figure 1.1.Physiological context.(a) The muscle contraction is generated by the interaction of myosin molecular motors with actin filaments in so-called contractile units.(b) Antagonist contractile units are stacked in parallel to form sarcomeres. (c) Sarcomeres are assembled in series to form the fibrils of the muscle cells.

Figure 1 . 2 .
Figure 1.2.(a) The Huxley-Hill model is parametrized by the distance s(t)separating the myosin anchor from its nearest binding site.The distance between two consecutive sites is denoted by d.After attachment, the formed cross-bridge is represented as a spring of length s(t).(b) The h-model considers the position of the detached myosin head X t relative to the same anchor point.The nearest binding site is now defined relative to the position of the head, and parametrized by the variable h t .
.2.Let T d denote the one-dimensional torus built by identifying the extremities of the interval [−d/2, +d/2].This torus is associated with the projection map π : [−d/2, +d/2] → T d , that we extend in a map π : R → T d .Together with π, we consider the coordinate map

Figure 5 . 1 .
Figure 5.1.h-model parameter functions.(a) Attached energy potential w 1 and detached energy potential w 0 .(b) Attachment rate.(c) Detachment rate.Expressions and parameters values can be found in Table D.1.

Figure 5 . 2 .
Figure 5.2.Evaluation of the calibrated h-model prediction against experimental data.Experimental data are obtained with rat cardiac intact muscles cells around 25 • C, see [25, Figure4(a)] for more details.Note that for low shortening velocities ( ẋc < 0.5 µm s −1 ), experimental data display instabilities[9,10,13,43], which are classically the subject of dedicated experimental investigations.Reproducing closely the experimental data in this regime is out of the scope of this paper.

Force production and Hill's thermodynamic formalism
Proposition 2.1 and Lemma 2.4 set the Huxley-Hill framework.As mentioned in the introduction, the original two-state model is not sufficient to reproduce all available experimental data, and other models considering several intermediate attached and detached states have been proposed in the literature [1,,formula for jump processes which is proven in e.g.[38, Proposition (4.6)] or[1,  Lemma 4.4.5].The bounds on coefficients ensure that the local martingales in the Ito formula are true martingales having zero expectancy.□ 22, Chapter 2.3].Equation (2.11) having a unique solution in the distribution sense, Equation (2.13) follows from the +d/2).Note that an alternative coordinate map T d → (−d/2, +d/2] could equivalently be considered.The position s t ∈ R of the closest actin site to The location s t of the closest actin site is now a random variable which can be discontinuous, depending on the displacements of X t on the full real line.On the contrary, h t is a continuous variable confined in a compact domain, making the SDE analysis simpler.The assumption of a finite window is thus replaced by a bound on the distance X t can jump at once, which cannot exceed d/2. model (4.9) and the Huxley-Hill model(2.1a) 1, s)p(t, 1, s)ds.

Table 5 . 1 .
Experimental and simulation results values of isometric physiological indicators: ratio of attached heads and force per attached head.The experimental values are obtained from rat cardiac isometric experiments.The ratio of attached head is taken from