# Varying thin filament activation in the framework of the Huxley'57 model

François Kimmig, Matthieu Caruel, Dominique Chapelle

► To cite this version:

François Kimmig, Matthieu Caruel, Dominique Chapelle. Varying thin filament activation in the framework of the Huxley'57 model. International Journal for Numerical Methods in Biomedical Engineering, 2022, 10.1002/cnm.3655. hal-03788978v2

HAL Id: hal-03788978

https://hal.science/hal-03788978v2

Submitted on 24 Oct 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Varying thin filament activation in the framework of the Huxley'57 model

François Kimmig¹,²*, Matthieu Caruel³, and Dominique Chapelle²,¹

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

²Inria, Palaiseau, France

³Univ Paris Est Creteil, Univ Gustave Eiffel, CNRS, UMR 8208, MSME, Créteil, France

*corresponding author: francois.kimmig@inria.fr

## Abstract

Muscle contraction is triggered by the activation of the actin sites of the thin filament by calcium ions. It results that the thin filament activation level varies over time. Moreover, this activation process is also used as a regulation mechanism of the developed force. Our objective is to build a model of varying actin site activation level within the classical Huxley'57 two-state framework. This new model is obtained as an enhancement of a previously proposed formulation of the varying thick filament activation within the same framework [1]. We assume that the state of an actin site depends on whether it is activated and whether it forms a cross-bridge with the associated myosin head, which results in four possible states. The transitions between the actin site states are controlled by the global actin sites activation level and the dynamics of these transitions is coupled with the attachment-detachment process. A preliminary calibration of the model with experimental twitch contraction data obtained at varying sarcomere lengths is performed.

**Keywords** – Cardiac modeling; Mathematical modeling; Thick and thin filament activation; Huxley'57 model

## Introduction

Cardiac muscle contraction originates from the interaction between two types of aligned filaments located at the microscale of the tissue: myosin filaments (also called thick filaments) and actin filaments (also called thin filaments). These two types of filaments, structured in an array, are the main components of the sarcomeres, which are the elementary contractile units of the muscle fibers. The myosin proteins forming the thick filaments have heads that cyclically bind to the actin sites of the adjacent thin filaments, forming the so-called *cross-bridges*. Once formed, cross-bridges exert an active traction force between the filaments. At rest, the attachment of the myosin heads is blocked by the troponin-tropomyosin complex. The electrical activation wave travelling through the muscle tissue triggers the release of calcium ions inside the cell. They bind to the troponin-tropomyosin complex, which results in a conformation change ultimately unblocking the actin sites. We call this transition the thin filament activation. The myosin heads can now attach and the force rises. Shortly after the release, calcium ions are taken off the cytosol, which contributes to the detachment of the calcium ions from the thin filament. Subsequently, the myosin heads detach from the actin sites, which finally leads to the tissue relaxation. This sequence of activation-deactivation of the actin sites constitutes a so-called *twitch contraction* in *ex vivo* experimental preparations mirroring the *in vivo* heartbeat. The force magnitude and the dynamics of the twitch contraction depend on the amount of calcium released in the cell (see Section 1). We call *level of calcium supply* the ratio between the amount of calcium actually released in the intracellular medium and the maximal amount that could be possibly released – that thus corresponds to the maximal force that can be generated.

In parallel to this activation-deactivation by calcium, the degree of activation of thin filament is also controlled by the sarcomere extension [2, 3, 4, 5]. This effect is called the length dependent activation
<figure><img src="image_2.png" /><figcaption>Figure 1: Presentation of our novel modeling framework. It aims to couple the actin-myosin interaction (c) and the thin and thick filament activations (d) in a mathematically consistent manner. (a) Structure of a sarcomere (~µm). Using the symmetries in the sarcomere structure, its description can be made considering a half-myosin filament. (b) Actin sites and myosin heads interact with each other to form cross-bridges (cycling rate ~10 ms). (e) Thick filament regulation. As a result of sarcomere extension, a myosin head transitions between a state where it is not available for attachment (indexed by 0) to a state in which it is available for attachment (indexed by 1). (f) Thin filament regulation. As a result of sarcomere extension and calcium supply, an actin site transitions between a permissive state (denoted by *P*) in which it is available for myosin attachment and a non-permissive state (denoted by *N*) in which is not available for myosin attachment.</figcaption></figure>

(LDA). Recent experimental results show that the length-dependent activation is an intrinsic mechanism [2, 5], meaning in particular that it operates without feedbacks from the formation of cross-bridges. The thick filament is also subjected to an activation process varying the number of myosin heads that are available for attachment controlled by the sarcomere extension [6, 4]. For both filaments, an increase in sarcomere extension is associated with an increase of the activation level, and therefore eventually with an increase of the developed force.

In addition to this intrinsic regulation by the sarcomere length, the thin filament is also subjected to extrinsic regulation mechanisms driven by the neuroendocrine system *in vivo*, that affect the way the thin filament responds to variations of the sarcomere length and that modulate the supply of calcium ions. As a starting point, we choose to focus in this paper on modeling the muscle behavior observed with *ex vivo* samples. These samples are separated from the neuroendocrine system, hence the extrinsic regulation mechanisms will not be considered in this work. Our model will thus aim at the mechanisms by which activation varies with the sarcomere length and with the level of calcium supply.

The mechanism by which activation varies with the sarcomere length and with the level of calcium supply has been widely investigated using mechanical tests in various activation conditions, possibly coupled with structural analysis by X-ray crystallography [1, 3, 7, 9, 11], but remains only partially understood [3, 4]. In particular, although it is now known that it is intrinsic, the exact origin of length-dependent activation of the thin filament remains subject to debate. The complete role of titin in this process is not fully elucidated either. In parallel, a controversy is still open about the implication in the thick filament regulation of the myosin head transition between the on- and off-states – in particular in the light of the dynamical aspects of this transition.

The investigation of the active force regulation can also be considered from a modeling point of view by formulating a simulation tool where the force is seen as a proxy to investigate the thick and thin filament activation processes. Most models that have been developed to that purpose are phenomenological, and derive from the seminal work of [12, 13, 14, 15] (see [16] for a review).

In the original papers [12, 13, 14, 15, 17], the actin sites are grouped into regulatory units (made of the seven consecutive actin sites that are covered by a single troponin-tropomyosin complex) and it is assumed that there is only one attachment site for the myosin heads per regulatory unit (RU). Having thus paired myosin heads and regulatory units, these formulations propose continuous-time mean-field Markov models mixing the description of a representative myosin head (attached or not) and the description of a representative regulatory unit (linked with calcium or not, activated or not, associated with a myosin
head or not). In these models, the transition rates between the model states can be modulated by the state itself to represent the cooperativity effects – i.e. the fact that the activation of some regulatory units may have a positive feedback effect on the other regulatory units, favoring their activation.

A subsequent major contribution is the model of J.J. Rice and co-authors [18], which explicitly introduces spatial extension of the cooperativity processes into the modeling framework, while focusing on the description of the thin filament alone. To account for the cooperativity effect observed experimentally, they propose to link the value of the transition rates between the states of individual sites to the configuration of the first adjacent regulatory units, with an increased probability of being activated when these units are already activated. With this new modeling ingredient, the thin filament cannot be modeled by a representative regulatory unit anymore. Instead, all regulatory units must be considered. This model was extended to incorporate a simplified myosin attachment-detachment process that, in particular, does not take into account the strain of the cross-bridges [19]. Variations have been proposed to incorporate specific mechanisms and test their ability to capture the observed cooperativity effects, such as the coupling of adjacent tropomyosins [20]. A difficulty that arises with these models is the large number of possible states. Indeed, a thin filament is made of 26 regulatory units [18], which leads to a number of possible states of $4^{26} = 4.5 \times 10^{15}$. A first method to overcome this difficulty (proposed in the original paper) consists in assuming that the thin filament is periodic, although it does not correspond to the actual structure. This assumption allows an analytical resolution of the system's steady-state. Several other approaches have then been proposed to tackle the issue. They include simulating the whole system using high performance computing [21] or a reduction of the model again with a single representative regulatory unit but with adjusted transitions rates so that cooperativity effects can be reproduced [22].

This modeling framework has further been refined by assuming that the level of overlap between thick and thin filaments impacts the transition rates in the thin filament model [23]. This hypothesis can lead to realistic predictions, which supports the idea that the length-dependent activation originates from mechanisms associated with the level of overlap between the thin and the thick filaments. Recently, a novel and powerful technique for the efficient simulation of the latter model by introducing wise simplifying assumptions has been proposed along with an improvement of the model formalization [24].

In addition to the main family of models, alternative approaches have been proposed to model the thin filament activation based on a mechano-chemical description. They propose to consider the proteins as a single long flexible chain whose deformations are associated with an elastic energy [25, 26, 27, 28] or to modulate the properties of the actin site with the mechanical force carried by the thin filament [29]. These approaches may seem more appealing than the phenomenological approach because they are related to the actual structure of the thin filament. However, since the origin of the length-dependent activation remains partially unknown, the establishment of non-phenomenological models and their evaluation with respect to phenomenological ones is difficult.

A drawback of the above approaches is the impossibility to use as actin-myosin interaction model the complete partial differential equation based model proposed by A.H. Huxley, which includes the effect of the relative actin and myosin filaments sliding on the states densities evolution [30]. The inclusion of the actin-myosin dynamics with a contribution of filaments sliding is important to model the sarcomere for two reasons. First, filaments sliding is a key phenomenon of the cardiac muscle behavior, in particular – in the context of heart modeling – in the blood ejection phase when the left ventricle shrinks. Secondly, there does not exist any experimental protocol capable of measuring the time transient variations of the thin filament activation – i.e. the conformational changes of the thin filament that allow myosin heads to bind to the actin sites – happening in a twitch contraction. Therefore, models have to be calibrated from force transient measurements, which, in part, depends on the relative filament sliding. There is thus the additional dynamics between the regulatory units activation and the myosin heads binding that has to be considered to avoid calibrating the models on biased data.

The aim of this paper is precisely to set the framework for combining detailed actin-myosin interaction models with filament activation models.

Many models of the actin-myosin interaction including the filament sliding contribution have been developed based on the work of A.F. Huxley [30] proposing various levels of refinement [31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47]. Our goal in this paper is to link models targeting the level of the thick and thin filament activation with the family of the Huxley'57 model in a rigorous way – which has never been done before in a general framework – in order to build a complete contraction model that takes into account both the actin-myosin interaction dynamics and the activation-regulation mechanisms
of the contraction. The point is not to bring new modeling elements for actin-myosin interaction or the
activation processes themselves but instead to build a framework able to couple them in a strong manner,
the thick and thin filament activations being strain-dependent. In particular, the intracellular calcium
concentration [Ca²⁺]ᵢ is neither a parameter nor a variable of our model, and we use as input the
filament activations. The developed modeling framework is illustrated in Figure 1. It is developed as
an enhancement of a previously proposed formulation of the varying thick filament activation within the
same framework [1]. We validate the potential of our approach by comparing the simulation results with
experimental twitch contractions signals using a state-of-the-art Huxley'57 model and surrogate models
of the thick and thin filament activation calibrated to phenomenologically reproduce some features of
the existing dedicated models. This approach could then pave the way for finely investigating potential
regulation underlying mechanisms and ultimately allowing to obtain a complete model of the sarcomere
that could be used in organ simulations.

Several works have already followed the path of embedding the thin filament activation in the Hux-
ley'57 model family, but they use a phenomenological approach, typically representing the degree of
filament activation as a multiplicative modulation of the attachment rate [48, 49, 1] or as an *ad hoc*
modification directly put in factor of the active force [50]. We propose here to overcome such limitations
with a formulation where thick and thin filament activations and the Huxley'57 model are explicitly
coupled.

The closest contribution in the literature to our work is the recent paper of F. Regazzoni and co-
authors [51]. Both approaches combine the thick and thin filaments regulation with the actin-myosin in-
teraction and use similar ground assumptions. They nevertheless differ in several major points. F. Regaz-
zoni and co-authors assume a particular mechanism for the thick filament regulation and focuses its
description on the actin site while our approach is compatible with many thick filament regulation
mechanisms, and it focuses on the myosin head. Moreover, in the modeling process, they introduce
assumptions along the derivation of the model while we consider a priori assumptions.

Another recent paper tackles a similar issue of coupling the activation of the thin filament and
the actin-myosin interaction [52]. This work introduces new modeling bricks in the form of additional
chemical states and provides a comprehensive description of the whole activation-contraction coupling.
In our work, we focus on a subpart of this coupling from the thin filament activation to the contraction.
To do so, we complement the seminal Huxley'57 partial differential equation to integrate the thick and
thin filament activation, while properly handling mass conservation. Moreover, our approach has the
advantage of being embedded into a multiscale framework paving the way for organ level simulations
[53].

This paper is organized as follows. Section 1 presents the experimental data that support our modeling ingredient choices and with which the model calibration is performed. In Section 2, we first present a review of the Huxley'57 model and its previously proposed extension that accounts for the thick filament activation [1]. Then, we introduce our description of the thin filament, and we derive our further enhancement of these models from the conservation of matter. In Section 3, we exhibit a possible calibration of our model that we validate with experimental data. Sections 4 and 5 then present the discussion of our results and our conclusions, respectively.

# 1 Physiological review

We recall here the essential properties of the thin filament activation that we will need for the development
of our model.

## 1.1 Activation of the thin filament

The baseline evolution of the thin filament activation (the conformation change of the thin filament
enabling myosin heads attachment) is a transient process. It cannot be directly measured experimentally
and must therefore be inferred from the force transient measurements through models. Note that the
sarcomere extension being an important parameter affecting the muscle contraction (see Section 1.2),
it must be controlled in the experimental setup. The *ex vivo* experimental results presented in this
paper are obtained in sarcomere length controlled conditions. While the muscle contracts, the sarcomere
<figure><img src="image_3.png" /><figcaption>Figure 2: (a) Comparison of the steady-state force (black symbols) and the twitch peak force (blue symbols) in high level of calcium supply. (Black triangle) intact tetanized rat cardiac cell at 26°C [59], (square) skinned rat cardiac cells at 22-24°C [60], (black diamond) skinned rat cardiac cells at 15°C [61], (blue circles) intact rat cardiac cells at 25°C [56], (blue squares) intact rat cardiac cells at 25°C [55]. In steady-state conditions, the variation of the force in maximal calcium activation conditions originates solely from the variation of the thick filament activation – i.e. the variations of the myosin heads availability. The normalized force thus allows to calibrate the function $n_0(e_c)$ (red line), which represents the ratio of available myosin heads (see also Section 2.1). (b) Variation of the thin filament activation with respect to the sarcomere extension for various concentrations of intracellular calcium obtained from intact tetanized rat cell force measurement and using the function $n_0(e_c)$ presented in (a). In the experimental conditions several levels of extracellular calcium concentration are considered. The resulting measured free intracellular concentration is reported near the curves. Approximate data fits are displayed with solid lines to enhance readability. Data from [59].</figcaption></figure>

length is measured and a feedback loop adjusts the length of the contracting muscle sample so that the sarcomere length remains constant [54, 55].

## 1.2 Regulation mechanisms

In this paper, we focus on the intrinsic regulation mechanisms, which can be observed on *ex vivo* samples. These regulation mechanisms operate through the sarcomere length and the supply of calcium.

The sarcomere length regulation affects both the thick and thin filaments activation but leaves the properties of the cross-bridges unchanged and in particular does not affect the force generated per cross-bridge [56, 57, 58]. For both the thick and the thin filaments, an increase in sarcomere extension results in an increase of the activation level, hence of the maximal force developed by cardiac muscles. The calcium supply affects the activation of the thin filament only and is positively linked to the developed active force.

### 1.2.1 Modulation of the maximum force

The effects of sarcomere length and amount of calcium supply on the maximum forces produced under steady state conditions and during twitch contraction are illustrated in Figure 2. We there define the sarcomere extension $e_c = SL/(2\ell_{hs}) - 1$, where SL denotes the sarcomere length, and $\ell_{hs}$ is the reference half-sarcomere length taken equal to 0.925 µm for rat cardiac muscle [54].

**Steady-state conditions** We first consider the steady-state regime, in which the thin filament activation remains at a constant level to discard the dynamical effects. These non-physiological conditions are obtained experimentally with skinned cells and tetanized intact cells. In these conditions, the force is proportional to the levels of thick and thin filaments activation.

To extract information on the thin filament activation from the force experimental data, the contribution of the varying thick filament activation must be eliminated. This contribution is identified by considering the maximum force obtained in steady state contraction at saturating calcium concentration,
see Figure 2(a, black symbols). In such conditions, the thin filament is maximally activated at all sarcomere lengths, and the variation of the maximum force is solely due to the thick filament regulation. We then identify the thick filament dependency of the activation level at maximum calcium concentration by the function $n_0(e_c)$ illustrated in Figure 2(a, red line).

In Figure 2(b, black lines and symbols), we show the steady state tension obtained at various non-saturating calcium concentration and at different sarcomere length. The curves are normalized by the function $n_0(e_c)$ to illustrate highlight the effect of the thin filament activation level.

These data show that the thin filament activation level always increases with the level of calcium supply, but in a way that depends on the sarcomere length. The origin of this effect remains partially unknown [4], but there are now evidences that it is a mechanism intrinsic to the thin filament [5].

**Non-steady-state conditions** Physiological cardiac muscle contractions are series of twitches, in which force generation and relaxation phases alternate. During a twitch, the peak force remains below the maximum steady-state force at all tested sarcomere lengths and calcium concentrations, see blue lines and symbols in Figure 2(a). This effect is due to the sarcomere length and calcium concentration dependence of the dynamics of the activation itself: the relaxation process starts before the steady-state activation level is reached.

### 1.2.2 Modulation of the force transients

The sarcomere length and the amount of calcium supply also have an impact on the time evolution of the thin filament activation as assessed by the observed variations of the force transient evolution, see Figure 3 (dashed lines). The traces are obtained from twitch contraction at imposed sarcomere length and two different calcium ions concentration. The time taken to reach the peak force is not strongly affected by the sarcomere length, but the relaxation duration is increased with increasing sarcomere lengths. These data also show that the level of calcium activation only slightly affects the time dependency of the force (after rescaling to its maximum value) in particular with a delay of the peak time at higher levels of calcium supply, but this feature is not observed in all experimental data, see [62]. This result suggests that the activation dynamics is marginally dependent on the level of calcium supply.

<figure><img src="image_4.png" /><figcaption>Figure 3: Comparison of simulated twitch contractions at various sarcomere extension with experimental data (from top to bottom, $e_c$ = {0.19, 0.162, 0.135, 0.108, 0.081, 0.054, 0.027} or SL = {2.2 µm, 2.15 µm, 2.1 µm, 2.05 µm, 2 µm, 1.95 µm, 1.9 µm}) (a) High contractility conditions, (b) low contractility conditions. The experimental data are obtained with a strict sarcomere length control throughout the experiment [55].</figcaption></figure>

As mentioned in Section 1.2, the dynamics of the actin-myosin interactions is left unchanged by sarcomere length variations and the level of calcium supply [56, 57, 58].

By applying rapid (5 ms) preconditioning length steps immediately before a twitch contraction and observing no changes in the time-course of the force, it has been showed that the regulation by the sarcomere length operates on a timescale much faster than both the thin filament activation (timescale
of 100 ms) and the myosin motors attachment-detachment dynamics (the shortest timescale involved being of 5 ms, see [47]) [63].

We conclude from this experimental analysis that both the thick filament and the thin filament activation variations due to changes in the sarcomere length can be considered to be instantaneous as a first approximation.

## 2 Model presentation

<figure><img src="image_5.png" /><figcaption>Figure 4: Representation of the main model ingredients. (a) Parametrization of the Huxley'57 model. The myosin and actin filaments are assumed to be rigid and slide past each other at a relative velocity $\dot{x}_c$. The myosin head can interact only with its closest actin site. The variable $s$ denotes the distance to this nearest actin site. The distance $s$ is thus defined in an interval $[s^-, s^+]$ of width $d_a$. The myosin head is presented in the detached state (black) and in the attached state (gray). (b) Definition of the transition rates between the attached and detached states. (c) The ratios of available myosin heads in the thick filament and of activated actin sites in the thin filament are given by the quantities $n_0$ and $n_a$, respectively. An increase in $n_0$ is associated with an increase in the number of available myosin heads. An increase in $n_a$ is associated with an increase in the number of activated actin sites. (d) For each myosin head, the interval in which it can reach actin sites is represented in blue. Due to the difference in length of periodicity between myosin heads and actin sites, some actin sites are thus not paired with any myosin head (actin site represented in gray).</figcaption></figure>
## 2.1 Actin-myosin interaction and thick filament activation

We start our model presentation by recalling the classical framework used to describe the actin-myosin interaction. This framework is based on the original Huxley’s 57 model [30], which considers that all myosin heads and all actin sites are activated. This framework also considers that the thick and thin filaments are rigid in the overlap region, which implies that the cross-bridges operate in parallel.

Considering the group of myosin heads located at a distance $s$ from their nearest actin site, the ratio of these heads that are attached at time $t$ is denoted by $P_1(s, t)$. The actin sites accessible for the myosin heads are supposed to be regularly located along the thin filament with a spatial period $d_a$, and therefore the value of the parameter $s$ can vary in a possibly non-symmetric interval $[s^-, s^+]$ with $s^+ - s^- = d_a$, see Figure 4(a). The dynamics of the system is given by

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

where $\dot{x}_c$ is the relative sliding velocity between the actin and myosin filaments (taken positive when the sarcomere length increases), $f$ and $g$ are the attachment and detachment rates, respectively, see Figure 4(b). The number of myosin heads per unit surface in a cross-section of half-sarcomere thickness ($\ell_{hs}$) is denoted by $\rho_{surf}$, and each myosin develops a force $\frac{dw_1}{ds}$, where $w_1$ is the free energy of the attached state. The active stress developed in a muscle cross-section of thickness $\ell_{hs}$ is thus given by [64, 47]

$$ T_c = \frac{\rho_{surf}}{d_a} \int_{s^{-}}^{s^{+}} \frac{dw_1}{ds}(s) P_1(s, t) ds. $$

The availability of the myosin heads is then introduced by considering two pools of heads [1]: those that are available for attachment and those that are not. Each head is characterized by an additional discrete parameter $\gamma$ describing the belonging to one of the two pools. The parameter $\gamma$ takes the value 1 in the pool of the available heads and the value 0 in the other pool. The fraction of heads that are available is denoted by $n_0$, as in Figure 2. In this work, $n_0$ will only depend on the sarcomere extension $e_c$ but could take a more general form depending on the thick filament activation model considered.

To preserve mathematical consistency, the function $n_0$ must belong to $[0,1]$ and we will consider that the heads that are not available can still attach—albeit with a much lower rate. In each pool, the probability of being attached is then denoted by $P_1(s, t, \gamma)$, and the attachment and detachment rates in each pool are denoted by $f_\gamma$ and $g_\gamma$, respectively. Defining

$$ |x|_{+} = \begin{cases} x, & \text{if } x \ge 0, \\ 0, & \text{otherwise,} \end{cases} \quad \text{and} \quad |x|_{-} = \begin{cases} -x, & \text{if } x \le 0, \\ 0, & \text{otherwise,} \end{cases} $$

the dynamics of the system is now given by (see [1] for the details)

$$ \left\{ \begin{aligned} \partial_t P_1(s, t, 1) + \dot{x}_c \partial_s P_1(s, t, 1) + \frac{|\dot{n}_0|_+}{n_0} [P_1(s, t, 1) - P_1(s, t, 0)] &= f_1(s)(1 - P_1(s, t, 1)) - g_1(s)P_1(s, t, 1), \\ \partial_t P_1(s, t, 0) + \dot{x}_c \partial_s P_1(s, t, 0) + \frac{|\dot{n}_0|_-}{1 - n_0} [P_1(s, t, 0) - P_1(s, t, 1)] &= f_0(s)(1 - P_1(s, t, 0)) - g_0(s)P_1(s, t, 0). \end{aligned} \right. $$

One can note that additional terms, compared with (1), appear in the equation to account for the transfer of myosin heads from one pool to the other as the thick filament activation function $n_0$ varies. The active tension is finally given as the weighted average between the two pools, hence we have

$$ T_c = \frac{\rho_{surf}}{d_a} \int_{s^{-}}^{s^{+}} \frac{dw_1}{ds}(s) \left[ n_0 P_1(s, t, 1) + (1 - n_0) P_1(s, t, 0) \right] ds. \quad (2) $$
## 2.2 Thin filament activation

We now extend this model to rigorously incorporate the variation of the thin filament activation.

### 2.2.1 Geometry of the system

A first step is to describe the geometry of our system. The longitudinal periodicity of the myosin heads along the thick filament is 43 nm while the periodicity along the helix direction is 14.3 nm [65]. Similarly, the thin filament also has two types of periodicity: the longitudinal periodicity 38.5 nm corresponding to the length of the regulatory units, and the 5.5 nm periodicity along the helix direction [65].

The thick and thin filament are organized in a pseudo-crystalline structure where a myosin head could possibly interact with several actin filaments and, conversely, an actin filament could be linked to multiple myosin filaments. In these conditions, it is not straightforward to define what are the relevant periodicity lengths that should be considered when describing a group of myosin heads and actin sites interacting with each other.

In this work, we consider that the periodicity of the myosin heads is 43 nm and that of the actin sites is $d_a = 38.5$ nm. Since the periodicity length is larger for the myosin heads, some actin sites are not paired with any myosin head, see Figure 4(d). These actin sites do not interact with myosin heads, and their states can thus be easily described: they are necessarily not part of a cross-bridge and their probability of being activated is equal to the global thin filament activation level. For this reason, we choose to center our description on the myosin heads, and it follows that the number of considered actin sites is the same as the number of myosin heads, although there are more actin sites in the real physical system.

If the relevant periodicity length were larger for actin sites than for myosin heads, the point of view of the description could be moved to the actin sites as in [51] but the modeling principles would not change.

The signed distance from a myosin head rest position to its nearest actin site is still denoted by $s$. This relation pairs a myosin head and an actin site. The signed distance between the actin site and its paired myosin head is denoted by $\bar{s}$ and satisfies the relation $\bar{s} = -s$ (see Figure 4(d)). Since we only considered actin sites that are paired with a myosin head, we can express every quantity as a function of the distance $s$.

### 2.2.2 Actin sites states and fluxes between them

We assume that each actin site has four possible states (see Figure 5). It is either

* non-activated and non-occupied by a myosin head,

* activated and non-occupied,

* non-activated but occupied by a myosin head,

* or activated and occupied.

For the subset of actin sites located at distance $s$ of their associated myosin head, we denote the respective ratios of each state at time $t$: $\hat{n}_a(s, t, \gamma)$, $\bar{n}_a(s, t, \gamma)$, $\tilde{n}_a(s, t, \gamma)$ and $\tilde{n}_a(s, t, \gamma)$. Naturally, these ratios must satisfy the condition

$$ \hat{n}_a(s, t, \gamma) + \bar{n}_a(s, t, \gamma) + \tilde{n}_a(s, t, \gamma) + \tilde{n}_a(s, t, \gamma) = 1, \forall t, \forall s \in [-s^-, s^+], \forall \gamma \in \{0, 1\}. $$

Actin sites can transition from one state to another. We denote the flux of actin sites from state A to state B normalized by the total number of myosin heads by $J_{A \to B}$ (we recall that in our model the description is made from the point of view of the myosin head and each myosin head is associated with a single actin site). We represent the four states of our model and the fluxes between them in Figure 5.

The ratio of activated actin sites is defined by $n_a = \bar{n}_a(s, t, \gamma) + \tilde{n}_a(s, t, \gamma) \forall \gamma \in \{0, 1\}$, which we consider as a prescribed parameter in our system. The ratio $n_a$ must belong to $|0,1|$ for mathematical consistency.

We assume that this parameter is homogeneous over the whole filaments and may thus depend explicitly on time and on other macroscale variables.
<figure><img src="image_6.png" /><figcaption>Figure 5: Actin model representation. The ratio associated with each state is indicated in blue.</figcaption></figure>

We further assume that when $n_a(t)$ varies, the changes occur indistinctly between occupied and non-occupied actin sites. With this assumption, we have

$$J_{na,xb \to a,xb}(s, t, \gamma) = \frac{|\dot{n}_a|_+}{\tilde{n}_a(s, t, \gamma) + \hat{n}_a(s, t, \gamma)} \tilde{n}_a(s, t, \gamma),$$

$$J_{a,xb \to na,xb}(s, t, \gamma) = \frac{|\dot{n}_a|_-}{\tilde{n}_a(s, t, \gamma) + \bar{n}_a(s, t, \gamma)} \tilde{n}_a(s, t, \gamma),$$

$$J_{na,nxb \to a,nxb}(s, t, \gamma) = \frac{|\dot{n}_a|_+}{\tilde{n}_a(s, t, \gamma) + \hat{n}_a(s, t, \gamma)} \hat{n}_a(s, t, \gamma),$$

$$J_{a,nxb \to na,nxb}(s, t, \gamma) = \frac{|\dot{n}_a|_-}{\tilde{n}_a(s, t, \gamma) + \bar{n}_a(s, t, \gamma)} \bar{n}_a(s, t, \gamma).$$

The remaining transition rates are related to the attachment-detachment process. Similarly to the fact that it is possible for a myosin head to attach when it is in the pool $\gamma = 0$, we still do not assume that it is impossible to bind on a non-activated actin site. However, the rate of this transition is lower than when the actin site is activated. We define the attachment rates with an activated actin site by $f_{\gamma,a}$ and by $f_{\gamma,na}$ when the actin site is not activated.

We now consider the flux of actin sites changing from an activated non-occupied state to an activated occupied state, which is equal to the flux of myosin heads attaching to activated actin sites. Myosin heads that can undergo this transition are a subgroup of the detached heads (the detached heads represent a ratio of $1 - P_1$), that are paired with an activated actin site. The ratio of myosin heads that, among the detached heads, are paired with an activated actin site is denoted by $r_{a|nxb}$. We thus have

$$J_{a,nxb \to a,xb}(s, t, \gamma) = f_{\gamma,a}(s)(1 - P_1(s, t, \gamma))r_{a|nxb} = f_{\gamma,a}(s)(1 - P_1(s, t, \gamma))\frac{\bar{n}_a(s, t, \gamma)}{\bar{n}_a(s, t, \gamma) + \hat{n}_a(s, t, \gamma)}.$$

Similarly, we obtain

$$J_{na,nxb \to na,xb}(s, t, \gamma) = f_{\gamma,na}(s)(1 - P_1(s, t, \gamma))r_{na|nxb} = f_{\gamma,na}(s)(1 - P_1(s, t, \gamma))\frac{\hat{n}_a(s, t, \gamma)}{\bar{n}_a(s, t, \gamma) + \hat{n}_a(s, t, \gamma)}.$$

Using the same reasoning for the actin sites paired with attached myosin heads (the attached myosin heads represent a ratio of $P_1$), we have

$$J_{a,xb \to a,nxb}(s, t, \gamma) = g_{\gamma,a}(s)P_1(s, t, \gamma)r_{a|xb} = g_{\gamma,a}(s)P_1(s, t, \gamma)\frac{\tilde{n}_a(s, t, \gamma)}{\tilde{n}_a(s, t, \gamma) + \tilde{n}_a(s, t, \gamma)},$$

$$J_{na,xb \to na,nxb}(s, t, \gamma) = g_{\gamma,na}(s)P_1(s, t, \gamma)r_{na|xb} = g_{\gamma,na}(s)P_1(s, t, \gamma)\frac{\tilde{n}_a(s, t, \gamma)}{\tilde{n}_a(s, t, \gamma) + \tilde{n}_a(s, t, \gamma)}.$$

### 2.2.3 Conservation laws and active tension

Having considered the fluxes between the different states, the dynamics of each population is given by the following system (the explicit variable dependency is omitted here in the right-hand side for the sake
of readability):

$$ \left\{ \begin{aligned} \frac{d}{dt} P_1(s, t, \gamma) &= f_{\gamma,a} \frac{\bar{n}_a}{\bar{n}_a + \hat{n}_a} (1 - P_1) + f_{\gamma,na} \frac{\hat{n}_a}{\bar{n}_a + \hat{n}_a} (1 - P_1) \\ &\quad - g_{\gamma,a} P_1 \frac{\tilde{n}_a}{\tilde{n}_a + \tilde{n}_a} - g_{\gamma,na} P_1 \frac{\tilde{n}_a}{\tilde{n}_a + \tilde{n}_a}, \end{aligned} \right. \quad (3a) $$

$$ \frac{d}{dt} \bar{n}_a(s, t, \gamma) = -f_{\gamma,a}(1-P_1) \frac{\bar{n}_a}{\bar{n}_a + \hat{n}_a} + g_{\gamma,a} P_1 \frac{\tilde{n}_a}{\tilde{n}_a + \tilde{n}_a} + \frac{|\dot{n}_a|_+}{\tilde{n}_a + \hat{n}_a} \hat{n}_a - \frac{|\dot{n}_a|_-}{\tilde{n}_a + \bar{n}_a} \bar{n}_a, \quad (3b) $$

$$ \frac{d}{dt} \tilde{n}_a(s, t, \gamma) = f_{\gamma,a}(1-P_1) \frac{\bar{n}_a}{\bar{n}_a + \hat{n}_a} - g_{\gamma,a} P_1 \frac{\tilde{n}_a}{\tilde{n}_a + \tilde{n}_a} + \frac{|\dot{n}_a|_+}{\tilde{n}_a + \hat{n}_a} \tilde{n}_a - \frac{|\dot{n}_a|_-}{\tilde{n}_a + \bar{n}_a} \tilde{n}_a, \quad (3c) $$

$$ \frac{d}{dt} \tilde{n}_a(s, t, \gamma) = f_{\gamma,na}(1-P_1) \frac{\hat{n}_a}{\bar{n}_a + \hat{n}_a} - g_{\gamma,na} P_1 \frac{\tilde{n}_a}{\tilde{n}_a + \tilde{n}_a} + \frac{|\dot{n}_a|_-}{\tilde{n}_a + \bar{n}_a} \tilde{n}_a - \frac{|\dot{n}_a|_+}{\tilde{n}_a + \hat{n}_a} \tilde{n}_a, \quad (3d) $$

$$ \frac{d}{dt} \hat{n}_a(s, t, \gamma) = -f_{\gamma,na}(1-P_1) \frac{\hat{n}_a}{\bar{n}_a + \hat{n}_a} + g_{\gamma,na} P_1 \frac{\tilde{n}_a}{\tilde{n}_a + \tilde{n}_a} + \frac{|\dot{n}_a|_-}{\tilde{n}_a + \bar{n}_a} \bar{n}_a - \frac{|\dot{n}_a|_+}{\tilde{n}_a + \hat{n}_a} \hat{n}_a. \quad (3e) $$

From this system, we obtain the desired relations by direct summation

$$ \hat{n}_a(s, t, \gamma) + \bar{n}_a(s, t, \gamma) + \tilde{n}_a(s, t, \gamma) + \tilde{n}_a(s, t, \gamma) = 1, \quad (4a) $$

$$ \tilde{n}_a(s, t, \gamma) + \bar{n}_a(s, t, \gamma) = n_a, \quad (4b) $$

$$ \tilde{n}_a(s, t, \gamma) + \hat{n}_a(s, t, \gamma) = P_1(s, t, \gamma). \quad (4c) $$

The first one is the conservation of actin sites; the second one corresponds to the definition that the ratio of activated actin sites (occupied or not) is given by $n_a$; and the third one reflects the fact that the ratio of attached head and the ratio of occupied actin sites is the same.

Using Equations (4), the system (3) can be reduced to

$$ \frac{d}{dt} P_1(s, t, \gamma) = f_{\gamma,a} \bar{n}_a + f_{\gamma,na} (1 - P_1 - \bar{n}_a) - g_{\gamma,a} (n_a - \bar{n}_a) - g_{\gamma,na} (\bar{n}_a - n_a + P_1), \quad (5a) $$

$$ \frac{d}{dt} \bar{n}_a(s, t, \gamma) = \frac{|\dot{n}_a|_+}{1 - n_a} (1 - P_1 - \bar{n}_a) + g_{\gamma,a} (n_a - \bar{n}_a) - \frac{|\dot{n}_a|_-}{n_a} \bar{n}_a - f_{\gamma,a} \bar{n}_a, \quad (5b) $$

$$ \hat{n}_a(s, t, \gamma) = 1 - \bar{n}_a - P_1, \quad (5c) $$

$$ \tilde{n}_a(s, t, \gamma) = n_a - \bar{n}_a, \quad (5d) $$

$$ \tilde{n}_a(s, t, \gamma) = P_1 - n_a - \bar{n}_a. \quad (5e) $$

The derivative used here is a total time derivative, meaning that it describes the time variations following a group of myosin heads. To fully establish the system dynamics, these total time derivatives need to be made explicit. The system is described from the point of view of the myosin heads, and therefore the pool exchange terms are written as in [1] and, since we have $\frac{ds}{dt} = -\frac{d\bar{s}}{dt} = \dot{x}_c$, the total time derivatives in (5) are given by

$$ \left\{ \begin{aligned} \frac{d}{dt} P_1(s, t, 1) &= \partial_t P_1(s, t, 1) + \dot{x}_c \partial_s P_1(s, t, 1) + \frac{|\dot{n}_0|_+}{n_0} [P_1(s, t, 1) - P_1(s, t, 0)], \\ \frac{d}{dt} P_1(s, t, 0) &= \partial_t P_1(s, t, 0) + \dot{x}_c \partial_s P_1(s, t, 0) + \frac{|\dot{n}_0|_-}{1 - n_0} [P_1(s, t, 0) - P_1(s, t, 1)], \\ \frac{d}{dt} \bar{n}_a(s, t, 1) &= \partial_t \bar{n}_a(s, t, 1) - \dot{x}_c \partial_s \bar{n}_a(s, t, 1) + \frac{|\dot{n}_0|_+}{n_0} [\bar{n}_a(s, t, 1) - \bar{n}_a(s, t, 0)], \\ \frac{d}{dt} \bar{n}_a(s, t, 0) &= \partial_t \bar{n}_a(s, t, 0) - \dot{x}_c \partial_s \bar{n}_a(s, t, 0) + \frac{|\dot{n}_0|_-}{1 - n_0} [\bar{n}_a(s, t, 0) - \bar{n}_a(s, t, 1)]. \end{aligned} \right. $$

The active force is still given by (2)

$$ T_c = \frac{\rho_{\text{surf}}}{d_a} \int_{s^{-}}^{s^{+}} \frac{d\bar{w}_1}{ds}(s) \left[ n_0 P_1(s, t, 1) + (1 - n_0) P_1(s, t, 0) \right] ds. $$
### 2.2.4 Example of a steady state isometric contraction

We consider as an example the steady-state regime of our model in isometric condition, i.e. we assume that $\dot{n}_0 = 0$, $\dot{n}_a = 0$ and $\dot{x}_c = 0$. The system (5) yields

$$ \left| \begin{array}{l} \bar{n}_a(s, \gamma) = \frac{g_{\gamma,a}(s)}{f_{\gamma,a}(s) + g_{\gamma,a}(s)} n_a, \\ P_1(s, \gamma) = n_a \frac{f_{\gamma,a}(s)}{f_{\gamma,a}(s) + g_{\gamma,a}(s)} + (1 - n_a) \frac{f_{\gamma,na}(s)}{f_{\gamma,na}(s) + g_{\gamma,na}(s)} \end{array} \right. $$

The global probability of being attached is then given by

$$ \begin{aligned} P_1(s) &= n_0 P_1(s, 1) + (1 - n_0) P_1(s, 0), \\ &= n_0 \left[ n_a \frac{f_{1,a}(s)}{f_{1,a}(s) + g_{1,a}(s)} + (1 - n_a) \frac{f_{1,na}(s)}{f_{1,na}(s) + g_{1,na}(s)} \right] \\ &\quad + (1 - n_0) \left[ n_a \frac{f_{0,a}(s)}{f_{0,a}(s) + g_{0,a}(s)} + (1 - n_a) \frac{f_{0,na}(s)}{f_{0,na}(s) + g_{0,na}(s)} \right]. \end{aligned} $$

The probability of being attached is the weighted average of the probability of being attached in each configuration (actin site activated or not, myosin head attached or not). With the natural approximations

$$ \frac{f_{0,a}}{g_{0,a}} \ll 1, \quad \frac{f_{0,na}}{g_{0,na}} \ll 1, \quad \frac{f_{1,na}}{g_{1,na}} \ll 1, $$

we obtain

$$ P_1(s) = n_0 n_a \frac{f_{1,a}(s)}{f_{1,a}(s) + g_{1,a}(s)}, $$

and the ratio of attached heads is finally given by

$$ n_{\text{att}} = \frac{1}{d_a} \int_{s^{-}}^{s^{+}} P_1(s) \, ds = n_0 n_a \int_{s^{-}}^{s^{+}} \frac{f_{1,a}(s)}{f_{1,a}(s) g_{1,a}(s)} \, ds. $$

The number of attached heads is therefore proportional to the level of thick filament activation and the level of thin filament activation, which is consistent with our modeling assumptions.

## 3 Model calibration

To validate our proposed model extension that takes into account the variations of the thin filament activation level, we need to specify the form of the thin filament activation function $n_a$. Our goal here is to give a preliminary calibration of the newly introduced modeling framework. To do so, we build surrogate models of the thick and thin filament activations, i.e. we prescribe an analytical form for the functions $n_0$ and $n_a$ such that they capture the basic features expected for the thick and thin filament activation models, respectively. In a further stage, these surrogate models should be replaced by specific predictive models.

To account for the available physiological data reviewed in Section 1, $n_a$ should be a function of the time, the sarcomere extension and the level of calcium supply. In experiments, the extracellular or intracellular calcium concentrations may be controlled. Since the interaction between the thin filament and the calcium ions are not modeled here, we cannot use these quantities as variables.

Instead, we introduce the parameter $C$ to represent the level of calcium supply, which we have defined as the ratio between the amount of calcium actually released in the cell and the maximal amount that could be possibly released. The parameter $C$ thus takes its values between 0 and 1: when $C = 0$, there is no calcium supply; when $C = 1$, the calcium supply is maximal. This parameter can be related in a bijective manner to the contractility $\Sigma_0$, which is the relevant clinical quantity of interest. This relation is given as follows: the contractility is defined as the maximal tension that can be developed for a given level of calcium supply – this tension is naturally obtained at the optimal sarcomere extension and in steady-state conditions.
Based on these considerations, we assume that the thin filament activation can be multiplicatively decomposed as

$$n_a(e_c, t, C) = n_{a,\infty}(e_c, C)n_{a,t}(e_c, t)$$

where $n_{a,\infty}$ accounts for the steady-state variation of the thin filament activation with the sarcomere extension and the level of calcium supply – scaled in such a way that it tends to 1 when $C$ becomes large – and where $n_{a,t}$ accounts for the transient evolution, which is affected by the sarcomere extension. The time transient function $n_{a,t}(e_c, t)$ is itself decomposed as $n_{a,t}(e_c, t) = n_{a,peak}(e_c)n_{a,t,norm}(e_c, t)$, where $n_{a,t,norm}$ represents the shape of the transient activation normalized by the peak value and $n_{a,peak}$ modulates peak activation level and accounts for the fact that the thin filament may not have time to reach its complete activation. The function $n_{a,peak}(e_c)$ can then be calibrated based on experimental data, see Figure 7(d). To avoid the overfitting of the twitch contraction data, we choose to have a time evolution that does not depend on the level of calcium supply. The details of the surrogate model functions definition are provided in Appendix A.

We choose the transition rates and the attached free energy level as calibrated in [47]. They are presented in Figure 6. Moreover, we choose the modeling assumptions

$$ \begin{cases} f_{\gamma,na} = \frac{1}{500} f_{\gamma,a} \quad \forall \gamma \in \{0,1\}, & f_{1,a} = \frac{1}{100} f_{0,a}, & f_{1,na} = \frac{1}{100} f_{0,na}, \\ g_{1,a} = g_{1,na} = g_{0,a} = g_{0,na}. \end{cases} $$

The functions $n_{a,\infty}$, $n_{a,t,norm}$ and $n_{a,peak}$ defined with these parameters along with the relation between the level of calcium supply $C$ and the contractility $\Sigma_0$ are presented in Figure 7.

<figure><img src="image_7.png" /><figcaption>Figure 6: (a) Attached state free energy level $w_1$. (b) Transition rates. The derivation of these quantities functions is presented in [47].</figcaption></figure>

We validate our model by comparing its outputs with the experimental data presented in Figure 2 (steady-state) and Figure 3 (transient). As a direct consequence of the calibration of the function $n_{a,\infty}$, the model in steady-state condition is able to match the observed relations between the developed force and the sarcomere extension at different level of calcium supply, see Figure 8.

We then simulate twitch contractions in “high” and “low” contractility conditions (we use a level of calcium supply $C$ of 0.95 and 0.46, respectively) in isometric conditions, i.e. with $\dot{x}_c = 0$, see Figure 3.

At both high and low level of calcium supply, the peak force obtained at the largest sarcomere extension ($e_c = 0.19$) is consistent with the value of the contractility. At such a sarcomere extension, we have $n_0 \approx n_{a,peak} \approx 1$. The contractility is thus directly controlled by the function $n_{a,\infty}$ corresponding to this extension, see Figure 7(a, light blue).

In high contractility conditions, we have $n_{a,\infty} \approx 1$ for $e_c > 0$, see Figure 7(a). Since all the tested sarcomere extensions are positive, the level of the peak force attained during the twitch contraction is essentially prescribed by the functions $n_0$ and $n_{a,peak}$. The former is calibrated beforehand, based on steady state data obtained with tetanized samples, see Figure 2. The latter is assumed to be independent of the contractility, see Figure 7(d), and adjusted precisely to match the peak force in the high contractility conditions, hence the good match with the observed peak force shown in Figure 3. Notice that the time to reach the peak force is the same for all sarcomere extensions as prescribed by the definition of
<figure><img src="image_8.png" /><figcaption>Figure 7: (a) Steady-state thin filament activation function $n_{a,\infty}$ for various sarcomere extensions. (b) Relation between the contractility $\Sigma_0$ and the level of calcium supply $C$. (c) Normalized thin filament activation function $n_{a,t,\text{norm}}$ for various sarcomere extensions. (d) Function $n_{a,\text{peak}}$ compared with experimental data. (Closed symbols) Peak twitch forces in high calcium supply conditions for cardiac rat cells normalized by the steady-state force predicted by the model – with the prior calibration of $n_0$ and $n_{a,\infty}$ – in the corresponding calcium activation conditions ($C = 0.95$); (circle) data from [56], (square) data from [55].</figcaption></figure>
<figure><img src="image_9.png" /></figure>

Figure 8: Comparison between the model steady-state force predictions at various levels of calcium supply as a function of the sarcomere extension and experimental data from [59] (black symbols). Each symbol corresponds to a single level of extracellular calcium concentration. The colored lines represent the active force normalized by its value at both maximum sarcomere extension and calcium supply, for various levels of calcium supply C. Note that the values of C chosen in the simulations are adjusted to match the data.

the function $n_{a,t,norm}$, see Figure 2(a). The relaxation dynamics is also well reproduced thanks to the dedicated calibration of $n_{a,t,norm}$ on the sarcomere stretch $e_c$, see Figure 7(c).

In the low contractility condition $C = 0.46$, the function $n_{a,∞}$ now has a significant dependency on $e_c$. Since $n_{a,peak}$ is already calibrated to reproduce the data at high contractility (Figure 3(a)), the simulated decay of the peak force with the sarcomere extension is larger than the observed one, especially at short sarcomere extension (see Figure 3(b)), due to the rapid decay of $n_{a,peak}$ with $e_c$ for $e_c < 0.1$, see Figure 8.

Note that no confidence intervals are provided in reference [55] from which the experimental data shown in Figure 3 are taken. However, some statistical estimates can be inferred for the late systolic phase duration, defined as the duration between the time at which the transient force reaches its peak and that at which it has decreased half-way from the peak value. More specifically, the standard deviation of the measured late systolic phase duration appears to be 38 ms and 16 ms in high and low contractility conditions, respectively. In comparison, the standard deviation of the differences in late systolic phase duration between the model predictions and the experiments are 9 ms and 6 ms in high and low contractility conditions, respectively, which shows that the model predictions are statistically relevant.

# 4 Discussion

## 4.1 Limits of the calibration

As observed above, the decrease of the peak force with the sarcomere extension is stronger in the simulations than in the data in low contractility condition, see Figure 8. The results could be improved by adjusting the value of $n_{a,∞}$ in the corresponding sarcomere extension interval ($e_c \approx 0.02 - 0.2$). However, this change would be rather artificial as the data used for this calibration are obtained with a range of sarcomere extensions between -0.15 and 0.05 whereas the twitch contractions are performed at sarcomere extensions between 0.02 and 0.19, i.e. in a region where the function $n_{a,∞}$ extrapolates the data. Of course this problem could easily be solved by making either $n_{a,peak}$ or $n_{a,t,norm}$ depend on C.

A reason why we did not pursue this fine-tuning of the calibration is that the data in steady-state conditions and the experimental twitch contractions have not been obtained in the same study, and there may thus be some intrinsic incompatibilities between them due to different experimental conditions.
It is also worth noting that the calibration presented in this paper only focuses on experimental data obtained from rat muscle samples at sub-physiological temperatures. To extend the scope of this work, confrontation with data obtained with human cardiomyocytes obtained at body temperature should be conducted [66], but these data are not available yet for the experimental setup required to perform the model validation.

In this theoretical work, we restrict our calibration to rat data to benefit from the vast number of available data sets and a lower dispersion of the data.

Moreover, our calibration considers that the actin-myosin interaction is not altered by variations of sarcomere length, as supported by experimental facts [56, 57, 58] (although some alternative points of view can be found in the literature, which account for the reduction of force due to the thin filament double overlap with a reduction of the effective myosin attachment rate [67]). In our modeling framework, the observed variation of the force due to sarcomere length variation is therefore solely attributed to the modulation of the thick and thin filaments activation. In particular, for the thick filament, myosin heads located in the thin filament double overlap double region are considered unactivated since they would on average not contribute to the force. This effect, possibly superimposed with a transition between the on- and off-states, is embedded in the calibration of the activation function $n_0$.

Finally, our model does not consider the intracellular calcium concentration [$Ca^{2+}]_i$ as a model parameter or variable. As a consequence, when using experimental data to calibrate the model, all variability with the sarcomere length is projected onto the filaments activation levels. It has been shown experimentally that the peak [$Ca^{2+}]_i$ is not affected by sarcomere length and so is thus our calibrated thin filament activation level at the peak of the twitch contraction [68, 69]. However, in the relaxation phase, variations of the order of 10% may appear in the case of a 7 times reduction of the force [69]. This difference would result in small variations of the force (through the thin filament activation level), which are not captured in our calibration. Nevertheless, the main variation of the twitch force magnitude is thus related to the thin and thick filament regulation, which validates the approach pursued in our paper.

## 4.2 Compared influences of cross-bridge and calcium activation dynamics on the twitch dynamics

<figure><img src="image_10.png" /><figcaption>Figure 9: Comparison of the time evolution of the thin filament activation function $n_a$ (blue lines) and the active tension $T_c$ (black lines). The simulations are performed in isometric conditions ($\dot{x}_c = 0$) with $e_c = 0.108$ and $C = 0.33$. We choose two sets of transitions rates. (a) Transition rates based on the force-velocity calibration from [47]. (b) Transition rates divided by a factor six with respect to (a). These rates are similar to those proposed by [70] that are calibrated on the tension development dynamics.</figcaption></figure>

We now compare the effects of the dynamics of the thin filament activation and the dynamics of actin-myosin interaction on the tension development. There is no consensus on the calibration of the transition rates of the Huxley’s77 model family [47]. These rates can be calibrated to match the isotonic sliding velocity data (force-velocity curve) or to match the rate of tension redevelopment after a quick sarcomere shortening at a constant level of activation [70].
The twitch contractions obtained with both calibration approaches are shown in Figure 9. If the calibration favors the match with the force-velocity relation – which corresponds to the choice made in Figure 3 – the transition rates are six times higher than if the data from tension redevelopment are used for calibration. In that case, the tension trace closely follows the activation trace (see Figure 9(a)) showing that the contraction dynamics is mostly controlled by activation. Here, the detailed modeling of the actin-myosin interaction and in particular of the filament sliding only brings a small contribution to the active force dynamics.

Keeping the same calibration of the thin filament activation $n_a$ (the point is not to match experimental data), we show the tension dynamics resulting from a six-fold reduction in the attachment-detachment rates in Figure 9(b). Not surprisingly, the delay between the thin filament activation and the active tension is longer when the attachment-detachment rate is slower. In this case, the interplay between the two dynamics, and thus the underlying models, has to be considered.

# 5 Conclusions

In this paper, we proposed a framework to rigorously incorporate the varying actin sites activation level into the classical Huxley'57 model family. This framework may be coupled with a separately defined thin filament activation model or, as done in this work, used as a standalone model with the time-dependent thin filament activation being provided as a model input. The model is calibrated on experimental twitch contractions obtained in controlled sarcomere length conditions and a preliminary calibration is performed. However, additional validations are required to fully assess the capabilities of the model, in particular in conditions where the sarcomere extension varies during the contraction as it is the case in a heartbeat.

Our novel framework opens new possibilities to develop thin filament activation models by allowing to couple them with an actin-myosin interaction model, therefore enabling a rigorous calibration. Moreover, in the context of heart simulation, our model may also prove to be an important element as it bridges the electrophysiological part and the mechanical part of the description. Once thoroughly calibrated and coupled within the heart simulation framework, our model may allow to investigate the relation between the force and calcium concentration measured *in vitro* and that required *in vivo* in physiological conditions. If a discrepancy appears (possibly due to the absence of neuroendocrine regulation in the model apart from its effects on the level of calcium supply), it will encourage the development of new experimental protocols to understand this discrepancy and allow the development of associated models.

Moreover, our modeling framework has the full capability to be embedded in a micro-macro coupled formulation of the muscle description. To perform this coupling, the sarcomere extension $e_c$ must be related to the deformation of the complete tissue at the macro-level [49, 53]. This paves the way for the use of this novel modeling framework as a basis for modeling the active muscle contraction for cardiomyocyte or complete heart simulations.

**Ethical statement** The authors declare no conflict of interest. They consciously assure that the material presented in the manuscript is the authors' own original work and reflects the authors' own research and analysis in a truthful and complete manner.

# A Appendix

In this section, we provide the detailed definition of the functions used to build the surrogate thin filament activation model. The model parameters used for the definition of these functions are presented in Table 1.

Getting our inspiration from experiments on skinned cells showing that the force varies with a sigmoid shape with respect to the level of calcium supply [60, 61], we choose a family of sigmoid functions to represent the variation of the force with respect to the contractility, parametrized by the sarcomere
extension. We define

$$n_{a,\infty}(e_c, C) = \frac{(C)^{n_{\text{Ha}}}}{(C)^{n_{\text{Ha}}} + (C_{50}(e_c))^{n_{\text{Ha}}}}, \quad \text{with}$$

$$C_{50}(e_c) = \begin{cases} \frac{C_{50}^{\text{M}} - C_{50}^{\text{L}}}{e_c - e_{c,50}^{\text{L}}}(C_{50}^{1} - e_{c,50}^{\text{L}}) + C_{50}^{\text{L}}, & \text{if } e_c < e_{c,50}^{1}, \\ C_{50}^{1}\varphi_1(e_c) + \bar{C}_{50}^{1}\varphi_2(e_c) + C_{50}^{2}\varphi_3(e_c) + \bar{C}_{50}^{2}\varphi_4(e_c) + C_{50}^{3}\varphi_5(e_c), & \text{if } e_{c,50}^{1} \le e_c \le e_{c,50}^{2}, \\ \frac{C_{50}^{\text{H}} - C_{50}^{\text{M}}}{e_{c,50}^{\text{H}} - e_{c,50}^{\text{M}}}(C_{50}^{2} - e_{c,50}^{\text{M}}) + C_{50}^{\text{M}}, & \text{if } e_{c,50}^{2} \le e_c, \end{cases}$$

where $C_{50}$ represents the effect of the sarcomere extension on the thin filament activation. It is chosen as a regularized piecewise linear function. The regularization interpolation functions $\varphi_i$ and interpolation coefficients are given by

$$\left\{ \begin{aligned} \varphi_1(e_c) &= \frac{(e_c - e_{c,50}^2)^2}{(e_c - e_{c,50}^2)^2} \left( 2 \frac{e_c - e_{c,50}^1}{e_{c,50}^1 - e_{c,50}^2} + 1 \right), \\ \varphi_2(e_c) &= \frac{(e_c - e_{c,50}^2)^2}{(e_{c,50}^1 - e_{c,50}^2)^2} (e_c - e_{c,50}^1), \\ \varphi_3(e_c) &= \frac{(e_c - e_{c,50}^1)^2}{(e_{c,50}^1 - e_{c,50}^2)^2} \left( 3 - 2 \frac{e_c - e_{c,50}^1}{e_{c,50}^1 - e_{c,50}^2} \right), \\ \varphi_4(e_c) &= \frac{(e_c - e_{c,50}^1)^2}{(e_{c,50}^1 - e_{c,50}^2)^2} (e_c - e_{c,50}^2), \\ \varphi_5(e_c) &= \varphi_1(e_c) \cdot \varphi_3(e_c), \end{aligned} \right.$$

and

$$\left\{ \begin{aligned} C_{50}^{1} &= \frac{C_{50}^{\text{M}} - C_{50}^{\text{L}}}{e_{c,50}^{\text{M}} - e_{c,50}^{\text{L}}}(C_{50}^{1} - e_{c,50}^{\text{L}}) + C_{50}^{\text{L}}, \\ C_{50}^{2} &= \frac{C_{50}^{\text{H}} - C_{50}^{\text{M}}}{e_{c,50}^{\text{H}} - e_{c,50}^{\text{M}}}(C_{50}^{2} - e_{c,50}^{\text{M}}) + C_{50}^{\text{M}}, \\ \bar{C}_{50}^{1} &= \frac{C_{50}^{\text{M}} - C_{50}^{\text{L}}}{e_{c,50}^{\text{M}} - e_{c,50}^{\text{L}}}, \\ \bar{C}_{50}^{2} &= \frac{C_{50}^{\text{H}} - C_{50}^{\text{M}}}{e_{c,50}^{\text{H}} - e_{c,50}^{\text{M}}}, \end{aligned} \right.$$
<table><caption>Table 1: Model parameters for the thin filament activation function.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>n<sub>Ha</sub></td><td>5.5</td><td>τ<sub>rise</sub></td><td>0.14 s</td></tr><tr><td>C<sup>H</sup><sub>50</sub></td><td>0.35</td><td>τ<sup>H</sup><sub>relax</sub></td><td>0.38 s</td></tr><tr><td>C<sup>M</sup><sub>50</sub></td><td>0.40</td><td>τ<sup>L</sup><sub>relax</sub></td><td>0.15 s</td></tr><tr><td>C<sup>L</sup><sub>50</sub></td><td>0.60</td><td>e<sup>H</sup><sub>c,relax</sub></td><td>0.18</td></tr><tr><td>e<sup>H</sup><sub>c,50</sub></td><td>0.252</td><td>e<sup>L</sup><sub>c,relax</sub></td><td>6.5 × 10<sup>-3</sup></td></tr><tr><td>e<sup>M</sup><sub>c,50</sub></td><td>0.05</td><td>t<sub>0</sub></td><td>0.117 s</td></tr><tr><td>e<sup>L</sup><sub>c,50</sub></td><td>-0.177</td><td>e<sup>L</sup><sub>c,t</sub></td><td>0.028</td></tr><tr><td>e<sup>1</sup><sub>c,50</sub></td><td>0</td><td>e<sup>M</sup><sub>c,t</sub></td><td>0.1895</td></tr><tr><td>e<sup>2</sup><sub>c,50</sub></td><td>0.1</td><td>n<sup>M</sup><sub>a,peak</sub></td><td>1</td></tr><tr><td>C<sup>3</sup><sub>50</sub></td><td>0</td><td>n<sup>L</sup><sub>a,peak</sub></td><td>0.50</td></tr><tr><td></td><td></td><td>n<sub>a,0</sub></td><td>1 × 10<sup>-3</sup></td></tr></tbody></table>

respectively. The transient part of the thin filament activation function is defined by

$$n_{a,t}(e_c, t) = n_{a,\text{peak}}(e_c)n_{a,t,\text{norm}}(e_c, t), \text{ with}$$

$$n_{a,t,\text{norm}}(e_c, t) = \begin{cases} \frac{1 - n_{a,0}}{2} \left( \sin \left( \frac{\pi}{\tau_{\text{rise}}} (t - t_0) - \frac{\pi}{2} \right) + 1 \right) + n_{a,0}, & \text{if } t_0 \le t \le t_0 + \tau_{\text{rise}}, \\ \frac{1 - n_{a,0}}{2} \left( \sin \left( \frac{\pi}{\tau_{\text{relax}}(e_c)} (t - t_0 - \tau_{\text{rise}}) + \frac{\pi}{2} \right) + 1 \right) + n_{a,0}, & \text{if } t_0 + \tau_{\text{rise}} \le t \le t_0 + \tau_{\text{rise}} + \tau_{\text{relax}}(e_c), \\ 0, & \text{otherwise,} \end{cases}$$

$$\tau_{\text{relax}}(e_c) = \frac{\tau_{\text{relax}}^{\text{H}} - \tau_{\text{relax}}^{\text{L}}}{e_{c,\text{relax}}^{\text{H}} - e_{c,\text{relax}}^{\text{L}}} (e_c - e_{c,\text{relax}}^{\text{L}}) + \tau_{\text{relax}}^{\text{L}},$$

$$n_{a,\text{peak}}(e_c) = \begin{cases} \frac{n_{a,\text{peak}}^{\text{M}} - n_{a,\text{peak}}^{\text{L}}}{e_{c,t}^{\text{M}} - e_{c,t}^{\text{L}}} (e_c - e_{c,t}^{\text{M}}) + n_{a,\text{peak}}^{\text{L}}, & \text{if } e_c < e_{c,t}^{\text{M}}, \\ n_{a,\text{peak}}^{\text{M}}, & \text{if } e_c \ge e_{c,t}^{\text{M}}, \end{cases}$$

where $\tau_{\text{rise}}$ is the constant activation rising time, $\tau_{\text{relax}}$ accounts for the varying relaxation duration, $n_{a,\text{peak}}$ represents the fact that the thin filament may not have time to reach the complete activation state for non-maximal sarcomere extensions and $n_{a,0}$, satisfying $0 < n_{a,0} \ll 1$, ensures that $n_a$ is never equal to zero, which is forbidden in the system (5).

## References

[1] F. Kimmig, P. Moireau, and D. Chapelle, “Hierarchical modeling of length-dependent force generation in muscles and associated thermodynamically-consistent numerical schemes.,” *Computational Mechanics*, vol. 68, no. 4, pp. 885–920, 2021.

[2] G. P. Farman, E. J. Allen, K. Q. Schoenfelt, P. H. Backx, and P. P. de Tombe, “The role of thin filament cooperativity in cardiac length-dependent calcium activation.,” *Biophysical Journal*, vol. 99, no. 9, pp. 2978–2986, 2010.

[3] P. de Tombe, R. Mateja, K. Tachampa, Y. Mou, G. Farman, and T. Irving, “Myofilament length dependent activation.,” *Journal of Molecular and Cellular Cardiology*, vol. 48, no. 5, pp. 851–858, 2010.
[4] P. de Tombe and H. ter Keurs, “Cardiac muscle mechanics: sarcomere length matters.,” *Journal of Molecular and Cellular Cardiology*, vol. 91, no. C, pp. 148–150, 2016.

[5] Y. Ait-Mou, K. Hsu, G. Farman, M. Kumar, M. Greaser, T. Irving, and P. de Tombe, “Titin strain contributes to the Frank–Starling law of the heart by structural rearrangements of both thin- and thick-filament proteins.,” *Proceedings of the National Academy of Sciences*, vol. 113, no. 8, pp. 2306–2311, 2016.

[6] D. G. Allen and J. C. Kentish, “The cellular basis of the length-tension relation in cardiac muscle,” *Journal of Molecular and Cellular Cardiology*, vol. 17, no. 9, pp. 821–840, 1985.

[7] J. P. Konhilas, T. C. Irving, and P. P. de Tombe, “Myofilament calcium sensitivity in skinned rat cardiac trabeculae.,” *Circulation Research*, vol. 90, no. 1, pp. 59–65, 2002.

[8] M. Reconditi, E. Brunello, L. Fusi, M. Linari, M. F. Martinez, V. Lombardi, M. Irving, and G. Piazzesi, “Sarcomere-length dependence of myosin filament structure in skeletal muscle fibres of the frog.,” *The Journal of Physiology*, vol. 592, no. 5, pp. 1119–1137, 2014.

[9] M. Reconditi, M. Caremani, F. Pinzauti, J. D. Powers, T. Narayanan, G. J. M. Stienen, M. Linari, V. Lombardi, and G. Piazzesi, “Myosin filament activation in the heart is tuned to the mechanical task.,” *Proceedings of the National Academy of Sciences*, vol. 114, no. 12, pp. 3240–3245, 2017.

[10] G. Piazzesi, M. Caremani, M. Linari, M. Reconditi, and V. Lombardi, “Thick filament mechano-sensing in skeletal and cardiac muscles: a common mechanism able to adapt the energetic cost of the contraction to the task.,” *Frontiers in Physiology*, vol. 9, 2018.

[11] M. Caremani, F. Pinzauti, J. D. Powers, S. Governali, T. Narayanan, G. J. M. Stienen, M. Reconditi, M. Linari, V. Lombardi, and G. Piazzesi, “Inotropic interventions do not change the resting state of myosin motors during cardiac diastole.,” *The Journal of General Physiology*, vol. 151, no. 1, pp. 53–65, 2019.

[12] A. Landesberg and S. Sideman, “Coupling calcium binding to troponin C and cross-bridge cycling in skinned cardiac cells.,” *The American Journal of Physiology*, vol. 266, no. 3 Pt 2, pp. H1260–71, 1994.

[13] A. Landesberg and S. Sideman, “Mechanical regulation of cardiac muscle by coupling calcium kinetics with cross-bridge cycling: a dynamic model.,” *The American Journal of Physiology*, vol. 267, no. 2 Pt 2, pp. H779–95, 1994.

[14] J. Rice, R. Winslow, and W. Hunter, “Comparison of putative cooperative mechanisms in cardiac muscle: length dependence and dynamic responses.,” *AJP: Heart and Circulatory Physiology*, vol. 276, no. 5, pp. H1734–H1754, 1999.

[15] M. V. Razumova, A. E. Bukatina, and K. B. Campbell, “Different myofilament nearest-neighbor interactions have distinctive effects on contractile behavior.,” *Biophysical Journal*, vol. 78, no. 6, pp. 3120–3137, 2000.

[16] N. Trayanova and J. Rice, “Cardiac electromechanical models: from cell to organ.,” *Frontiers in Physiology*, vol. 2, p. 43, 2011.

[17] Y. Yaniv, W. C. Stanley, G. M. Saidel, M. E. Cabrera, and A. Landesberg, “The role of Ca2+ in coupling cardiac metabolism with regulation of contraction: in silico modeling.,” *Annals of the New York Academy of Sciences*, vol. 1123, no. 1, pp. 69–78, 2008.

[18] J. Rice, G. Stolovitzky, Y. Tu, and P. de Tombe, “Ising model of cardiac thin filament activation with nearest-neighbor cooperative interactions.,” *Biophysical Journal*, vol. 84, no. 2, pp. 897–909, 2003.

[19] J. Rice and P. de Tombe, “Approaches to modeling crossbridges and calcium-dependent activation in cardiac muscle.,” *Progress in Biophysics and Molecular Biology*, vol. 85, no. 2-3, pp. 179–195, 2004.
[20] S. G. Campbell, F. V. Lionetti, K. S. Campbell, and A. D. McCulloch, "Coupling of adjacent tropomyosins enhances cross-bridge-mediated cooperative activation in a Markov model of the cardiac thin filament.," *Biophysical Journal*, vol. 98, no. 10, pp. 2254–2264, 2010.

[21] J. Hussan, P. de Tombe, and J. Rice, "A spatially detailed myofilament model as a basis for large-scale biological simulations.," *IBM Journal of Research and Development*, vol. 50, no. 6, pp. 583–600, 2006.

[22] J. Rice, F. Wang, D. Bers, and P. de Tombe, "Approximate model of cooperative activation and crossbridge cycling in cardiac muscle using ordinary differential equations.," *Biophysical Journal*, vol. 95, no. 5, pp. 2368–2390, 2008.

[23] T. Washio, J.-I. Okada, S. Sugiura, and T. Hisada, "Approximation for cooperative interactions of a spatially-detailed cardiac sarcomere model.," *Cellular and Molecular Bioengineering*, vol. 5, no. 1, pp. 113–126, 2011.

[24] F. Regazzoni, L. Dedè, and A. Quarteroni, "Active contraction of cardiac cells: a reduced model for sarcomere dynamics with cooperative interactions.," *Biomechanics and Modeling in Mechanobiology*, vol. 17, no. 6, pp. 1663–1686, 2018.

[25] M. Geeves, H. Griffiths, S. Mijailovich, and D. Smith, "Cooperative [Ca2+]-dependent regulation of the rate of myosin binding to actin: solution data and the tropomyosin chain model.," *Biophysical Journal*, vol. 100, no. 11, pp. 2679–2687, 2011.

[26] S. Mijailovich, O. Kayser-Herold, X. Li, H. Griffiths, and M. Geeves, "Cooperative regulation of myosin-S1 binding to actin filaments by a continuous flexible Tm–Tn chain.," *European Biophysics Journal*, vol. 41, no. 12, pp. 1015–1032, 2012.

[27] N. Metalnikova and A. Tsaturyan, "A mechanistic model of Ca regulation of thin filaments in cardiac muscle.," *Biophysical Journal*, vol. 105, no. 4, pp. 941–950, 2013.

[28] S. Land and S. Niederer, "A spatially detailed myofilament model as a basis for large-scale biological simulations.," *PLOS Computational Biology*, vol. 11, no. 8, pp. e1004376–28, 2015.

[29] L. J. Dupuis, J. Lumens, T. Arts, and T. Delhaas, "Mechano-chemical interactions in cardiac sarcomere contraction: a computational modeling study.," *PLOS Computational Biology*, vol. 12, no. 10, pp. e1005126–20, 2016.

[30] A. F. Huxley, "Muscle structures and theories of contraction.," *Progress in Biophysics and Biophysical Chemistry*, vol. 7, pp. 255–318, 1957.

[31] A. F. Huxley and R. M. Simmons, "Proposed mechanism of force generation in striated muscle.," *Nature*, vol. 233, no. 5321, pp. 533–538, 1971.

[32] T. L. Hill, "Theoretical formalism for the sliding filament model of contraction of striated muscle Part I.," *Progress in Biophysics and Molecular Biology*, vol. 28, pp. 267–340, 1974.

[33] T. L. Hill, "Theoretical formalism for the sliding filament model of contraction of striated muscle part II.," *Progress in Biophysics and Molecular Biology*, vol. 29, pp. 105–159, 1976.

[34] E. Eisenberg and T. L. Hill, "A cross-bridge model of muscle contraction.," *Progress in Biophysics and Molecular Biology*, vol. 33, no. 1, pp. 55–82, 1978.

[35] E. Eisenberg, T. L. Hill, and Y. Chen, "Cross-bridge model of muscle contraction. Quantitative analysis.," *Biophysical Journal*, vol. 29, no. 2, pp. 195–227, 1980.

[36] G. Piazzesi and V. Lombardi, "A cross-bridge model that is able to explain mechanical and energetic properties of shortening muscle.," *Biophysical Journal*, vol. 68, no. 5, pp. 1966–1979, 1995.

[37] G. I. Zahalak, "The two-state cross-bridge model of muscle is an asymptotic limit of multi-state models.," *Journal of Theoretical Biology*, vol. 204, no. 1, pp. 67–82, 2000.
[38] D. A. Smith, M. A. Geeves, J. Sleep, and S. M. Mijailovich, “Towards a unified theory of muscle contraction. I: foundations.,” *Annals of Biomedical Engineering*, vol. 36, no. 10, pp. 1624–1640, 2008.

[39] L. Marcucci and L. Truskinovsky, “Mechanics of the power stroke in myosin II.,” *Physical Review E*, vol. 81, no. 5, pp. 051915–8, 2010.

[40] M. Caruel, J.-M. Allain, and L. Truskinovsky, “Muscle as a metamaterial operating near a critical point.,” *Physical Review Letters*, vol. 110, no. 24, p. 248103, 2013.

[41] R. Sheshka and L. Truskinovsky, “Power-stroke-driven actomyosin contractility.,” *Physical Review E*, vol. 89, no. 1, p. 012708, 2014.

[42] M. Caremani, L. Melli, M. Dolfi, V. Lombardi, and M. Linari, “Force and number of myosin motors during muscle shortening and the coupling with the release of the ATP hydrolysis products.,” *The Journal of Physiology*, vol. 593, no. 15, pp. 3313–3332, 2015.

[43] A. Månsson, D. Rassier, and G. Tsiavaliaris, “Poorly understood aspects of striated muscle contraction.,” *BioMed Research International*, vol. 2015, 2015.

[44] S. M. Mijailovich, O. Kayser-Herold, B. Stojanovic, D. Nedic, T. C. Irving, and M. A. Geeves, “Three-dimensional stochastic model of actin–myosin binding in the sarcomere lattice.,” *The Journal of General Physiology*, vol. 148, no. 6, pp. 459–488, 2016.

[45] M. Caruel, P. Moireau, and D. Chapelle, “Stochastic modeling of chemical-mechanical coupling in striated muscles.,” *Biomechanics and Modeling in Mechanobiology*, vol. 18, no. 3, pp. 563–587, 2019.

[46] A. Månsson, M. Persson, N. Shalabi, and D. E. Rassier, “Nonlinear actomyosin elasticity in muscle?,” *Biophysical Journal*, vol. 116, no. 2, pp. 330–346, 2019.

[47] F. Kimmig and M. Caruel, “Hierarchical modeling of force generation in cardiac muscle.,” *Biomechanics and Modeling in Mechanobiology*, vol. 19, no. 6, pp. 2567–2601, 2020.

[48] G. I. Zahalak and I. Motabarzadeh, “A re-examination of calcium activation in the Huxley cross-bridge model.,” *Journal of Biomechanical Engineering*, vol. 119, no. 1, pp. 20–29, 1997.

[49] D. Chapelle, P. Moireau, M. Sorine, and P. Le Tallec, “Energy-preserving muscle tissue model: formulation and compatible discretizations.,” *International Journal for Multiscale Computational Engineering*, vol. 10, no. 2, pp. 189–211, 2012.

[50] G. I. Zahalak, “A distribution-moment approximation for kinetic theories of muscular contraction.,” *Mathematical Biosciences*, vol. 55, no. 1-2, pp. 89–114, 1981.

[51] F. Regazzoni, L. Dedè, and A. Quarteroni, “Biophysically detailed mathematical models of multiscale cardiac active mechanics.,” *PLOS Computational Biology*, vol. 16, no. 10, p. e1008294, 2020.

[52] S. M. Mijailovich, M. Prodanovic, C. Poggesi, M. A. Geeves, and M. Regnier, “Multiscale modeling of twitch contractions in cardiac trabeculae.,” *The Journal of General Physiology*, vol. 153, no. 3, p. e202012604, 2021.

[53] F. Kimmig, D. Chapelle, and P. Moireau, “Thermodynamic properties of muscle contraction models and associated discrete-time principles.,” *Advanced Modeling and Simulation in Engineering Sciences*, vol. 6, no. 1, p. 6, 2019.

[54] H. E. D. J. Ter Keurs, W. H. Rijnsburger, R. Van Heuningen, and M. J. Nagelsmit, “Tension Development and Sarcomere Length in Rat Cardiac Trabeculae: Evidence of Length-Dependent Activation.,” in *Cardiac Dynamics*, pp. 25–36, Dordrecht: Springer Netherlands, 1980.

[55] P. M. Janssen and W. C. Hunter, “Force, not sarcomere length, correlates with prolongation of isosarcometric contraction.,” *AJP: Heart and Circulatory Physiology*, vol. 269, no. 2, pp. H676–H685, 1995.
[56] M. Caremani, F. Pinzauti, M. Reconditi, G. Piazzesi, G. J. M. Stienen, V. Lombardi, and M. Linari, "Size and speed of the working stroke of cardiac myosin in situ.," *Proceedings of the National Academy of Sciences*, vol. 113, no. 13, pp. 3675-3680, 2016.

[57] D. Amiad Pavlov and A. Landesberg, "The cross-bridge dynamics is determined by two length-independent kinetics: implications on muscle economy and Frank-Starling law.," *Journal of Molecular and Cellular Cardiology*, vol. 90, pp. 94-101, 2016.

[58] F. Pinzauti, I. Pertici, M. Reconditi, T. Narayanan, G. J. M. Stienen, G. Piazzesi, V. Lombardi, M. Linari, and M. Caremani, "The force and stiffness of myosin motors in the isometric twitch of a cardiac trabecula and the effect of the extracellular calcium concentration.," *The Journal of Physiology*, vol. 596, no. 13, pp. 2581-2596, 2018.

[59] H. E. D. J. Ter Keurs, T. Shinozaki, Y. M. Zhang, M. L. Zhang, Y. Wakayama, Y. Sugai, Y. Kagaya, M. Miura, P. A. Boyden, B. D. M. Stuyvers, and A. Landesberg, "Sarcomere mechanics in uniform and non-uniform cardiac muscle: A link between pump function and arrhythmias.," *Progress in Biophysics and Molecular Biology*, vol. 97, no. 2-3, pp. 312-331, 2008.

[60] J. C. Kentish, H. Ter Keurs, and L. Ricciardi, "Comparison between the sarcomere length-force relations of intact and skinned trabeculae from rat right ventricle. Influence of calcium concentrations on these relations.," *Circulation Research*, 1986.

[61] D. P. Dobesh, J. P. Konhilas, and P. P. de Tombe, "Cooperative activation in cardiac muscle: impact of sarcomere length.," *AJP: Heart and Circulatory Physiology*, vol. 282, no. 3, pp. H1055-H1062, 2002.

[62] R. Van Heuningen, W. H. Rijnsburger, and H. E. ter Keurs, "Sarcomere length control in striated muscle.," *The American Journal of Physiology*, vol. 242, no. 3, pp. H411-20, 1982.

[63] R. D. Mateja and P. P. de Tombe, "Myofilament length-dependent activation develops within 5&nbsp;ms in Guinea-pig myocardium-," *Biophysical Journal*, vol. 103, no. 1, pp. L13-L15, 2012.

[64] T. L. Hill, *Free Energy Transduction in Biology*. Academic press, 1977.

[65] R. Craig and R. Padrón, "Molecular structure of the sarcomere.," *Myology*, vol. 3, pp. 129-144, 2004.

[66] S. Land, S.-J. Park-Holohan, N. Smith, C. dos Remedios, J. Kentish, and S. Niederer, "A model of cardiac contraction based on novel measurements of tension development in human cardiomyocytes.," *Journal of Molecular and Cellular Cardiology*, vol. 106, pp. 68-83, 2017.

[67] S. Mijailovich, B. Stojanovic, D. Nedic, M. Svicevic, M. Geeves, T. Irving, and H. Granzier, "Nebulin and titin modulate cross-bridge cycling and length-dependent calcium sensitivity," *The Journal of General Physiology*, vol. 151, no. 5, pp. 680-704, 2019.

[68] D. G. Allen and S. Kurihara, "The effects of muscle length on intracellular calcium transients in mammalian cardiac muscle.," *The Journal of Physiology*, vol. 327, June 1982.

[69] P. H. Backx and H. E. ter Keurs, "Fluorescent properties of rat cardiac trabeculae microinjected with fura-2 salt," *AJP: Heart and Circulatory Physiology*, vol. 264, pp. H1098-H1110, Apr. 1993.

[70] P. de Tombe and G. Stienen, "Impact of temperature on cross-bridge cycling kinetics in rat myocardium.," *The Journal of Physiology*, vol. 584, no. 2, pp. 591-600, 2007.