Biophysically detailed mathematical models of multiscale cardiac active mechanics
BBiophysically detailed mathematical models ofmultiscale cardiac active mechanics
F. Regazzoni ∗ , L. Ded`e , A. Quarteroni MOX - Dipartimento di Matematica, Politecnico di Milano,P.zza Leonardo da Vinci 32, 20133 Milano, Italy Mathematics Institute, ´Ecole Polytechnique F´ed´erale de Lausanne,Av. Piccard, CH-1015 Lausanne, Switzerland (Professor Emeritus)
Abstract
We propose four novel mathematical models, describing the microscopic mech-anisms of force generation in the cardiac muscle tissue, which are suitable formultiscale numerical simulations of cardiac electromechanics. Such models arebased on a biophysically accurate representation of the regulatory and contractileproteins in the sarcomeres. Our models, unlike most of the sarcomere dynamicsmodels that are available in the literature and that feature a comparable rich-ness of detail, do not require the time-consuming Monte Carlo method for theirnumerical approximation. Conversely, the models that we propose only requirethe solution of a system of PDEs and/or ODEs (the most reduced of the fouronly involving 20 ODEs), thus entailing a significant computational efficiency.By focusing on the two models that feature the best trade-off between detail ofdescription and identifiability of parameters, we propose a pipeline to calibratesuch parameters starting from experimental measurements available in litera-ture. Thanks to this pipeline, we calibrate these models for room-temperaturerat and for body-temperature human cells. We show, by means of numericalsimulations, that the proposed models correctly predict the main features offorce generation, including the steady-state force-calcium and force-length re-lationships, the tension-dependent prolongation of twitches, the force-velocityrelationship. Moreover, they correctly reproduce the Frank-Starling effect, whenemployed in multiscale 3D numerical simulation of cardiac electromechanics.
Keywords
In silico models, Cardiac modeling, Active stress, Sarcomeres, Multiscalemodeling
Cardiovascular mathematical and numerical models are increasingly used, with atwofold role [15, 20, 26, 27, 72, 90]. On the one hand, realistic and detailed in silicomodels of the heart can deepen the understanding of its function, help the interpreta-tion of experimental observations and explain the delicate links between the organ-levelemergent phenomena and the underlying biophysical mechanisms. On the other hand,patient-specific numerical simulations, which are increasingly becoming available, can ∗ Corresponding author. Email address: [email protected] l e c t r o p h y s i o l o g y M e c h a n i c s Tropomyosin permissivenon-permissive
Troponin-C
Ca-boundCa-unbound
Crossbridges (actin+myosin)
Lymn – Taylor cycle
Thin filament activation Crossbridge cycling
Strain,Rate of strain
Sec. 2.1 Sec. 2.2Sec. 2.3
Regulatory Unit (RU)
Figure 1: Representation of the different stages of the force generation mechanism andof the sections where they are discussed in this paper.provide clinicians with valuable quantitative information to improve patient care andto support decision-making procedures.The construction of an integrated mathematical and numerical model of cardiacelectromechanics (EM) is however a remarkably arduous task. This is mainly dueto the multiphysics (due to the interplay of biochemistry, electricity, solid mechanics,fluid dynamics) and multiscale nature of the heart: characteristic spatial scales rangefrom nanometers to centimeters and the temporal ones from microseconds to seconds.This makes it difficult to devise computationally efficient and accurate algorithms fora plurality of mathematical models featuring a broad degree of details [3, 15, 21, 30,74, 87].The contrasting needs between model accuracy and computational efficiency ofnumerical simulations is mainly due to the multiscale nature of the heart, for whichthe mechanical work enabling the macroscopic motion of the organ is fueled by theenergy consumed at the microscale by subcellular mechanisms. The generation ofactive force takes place inside sarcomeres and involves a complex chain of chemicaland mechanical reactions. This mechanism can be split into two steps, that we sketchin Fig. 1. First, a ionic signal (specifically, an increase of calcium ions concentration)triggers the so-called regulatory units , protein complexes consisting of troponin andtropomyosin, that act as on-off switches for the muscle contraction. Then, when theregulatory units are activated, the actin and myosin proteins are free to interact andform the so-called crossbridges , molecular motors that generate an active force byconsuming the chemical energy stored in ATP [7, 47].Microscopic force generation includes many regulatory mechanisms, forming thesubcellular basis of organ-level phenomena, such as the Frank-Starling effect [47].Hence, if a microscale mathematical model of force generation is used in a multi-scale setting to build an integrated organ-level EM model, then it should be able toreproduce the above-mentioned mechanisms.In the past decades, several efforts have been dedicated to the construction ofmathematical models describing the complex dynamics of the processes taking placein sarcomeres [12–14, 16, 39, 41, 42, 56, 57, 61, 62, 80, 82, 84, 101, 102]. How-ever, because of the intrinsic complexity of the phenomenon of force generation, huge2omputational costs are associated with the numerical approximation of such mod-els, thus limiting their application within multiscale EM simulations. Despite severalattempts to capture the fundamental mechanisms underlying the force generation phe-nomenon into a tractable number of equations [10, 55, 75, 82, 84, 86, 100], the existingorgan-level cardiac mathematical models rely on two alternative strategies to describemicroscopic force generation. • Phenomenological models (see e.g. [38, 56, 57, 62, 83]) are built by fittingthe measured data with simple laws, chosen by the modeler. However, theparameters characterizing phenomenological models often lack a clear physicalinterpretation; moreover, the noisy nature and deficiency of data coming from thesubcellular contractile units and the intrinsic difficulties in measuring sarcomeresunder the conditions occurring during an heartbeat hamper the predictive powerof such models [26]. • Biophysically detailed models are based on an accurate description of the pro-teins involved in the force generation process and are derived from physics firstprinciples. However, their numerical solution, because of their complexity, istypically obtained by means of a Monte Carlo (MC) approximation (see e.g.[39, 101, 102]). The MC method is in fact inefficient, featuring a huge compu-tational cost, both in terms of time and memory storage. Indeed, to accuratelyapproximate the solution of a single heartbeat for a single myofilament, tens ofhours of computational time may be required; see e.g. [76].The purpose of this paper is to develop a biophysically detailed model for active forcegeneration, that explicitly describes the fundamental ingredients of the force genera-tion apparatus, yet featuring a tractable computational cost, so that it is suitable formultiscale EM simulations.
This paper is structured as follows. In Sec. 2 we recall the main features of the forcegeneration phenomenon in cardiomyocytes and the main difficulties encountered inthe construction of mathematical models describing the associated mechanisms. InSec. 3 we present the models proposed in this paper and in Sec. 4 we describe thestrategy employed for their calibration. Then, in Sec. 5, we show some numericalresults obtained with the proposed models, including filament-level 0D simulationsand multiscale 3D cardiac EM simulations. We provide some final remarks in Sec. 6.In Tab. 1 we provide a list of the abbreviations used throughout this paper.
We recall the microscopical mechanisms by which active force is generated in thecardiac tissue and we highlight the difficulties, rooted in the their intrinsic complexity,in describing such phenomena with a tractable number of equations. We also reviewthe main contributions available in the literature.Sarcomeres are cylindrically-shaped, 2 µ m length units made of a regular arrange-ment of thick and thin filaments. The former, also known as myosin filaments (MF),are mainly made of the protein myosin, while the latter are made of actin, troponin(Tn) and tropomyosin (Tm) and are also called actin filaments (AF).3 natomical terms Mathematical models LV Left ventricle H57 Model of [42]Tm Tropomyosin TTP06 Model of [93]Tn TroponinMH Myosin head
Others abbreviations
XB Crossbridge ODE Ordinary differential equationRU Regulatory unit PDE Partial differential equationBS Binding site CTMC Continuous-time Markov ChainMF Myosin (thick) filament MC Monte CarloAF Actin (thin) filamentEM ElectromechanicsTable 1: List of abbreviations of this paper.In Sec. 2.1 we deal with the activation of the thin filament, involving the troponin-tropomyosin regulatory units (RUs). Then, in Sec. 2.2, we address the crossbridge(XB) cycling and finally, in Sec. 2.3, we consider the full-sarcomere dynamics (seeFig. 1).
The activation of the thin filament is regulated by two factors, namely the intracellularcalcium ions concentration ([Ca ] i ) and the sarcomere length ( SL ). The experimentalsignature of the regulation mechanisms is given by the steady-steady relationshipsbetween calcium, sarcomere length and generated force.The force-calcium relationship (see Fig. 2) features a sigmoidal shape, well de-scribed by the Hill equation [22, 50, 96]: T isoa = T maxa (cid:16) EC [Ca ] i (cid:17) n H , (1)where T maxa is the maximum tension (tension at saturating calcium levels), EC isthe half maximal effective concentration (i.e. the calcium concentration producinghalf maximal force) and n H is the Hill coefficient. The experimentally measuredforce-calcium curves in the cardiac tissue feature a steep slope in correspondence ofhalf activation (Hill coefficient greater than one) [4, 22, 96], thus revealing the presenceof cooperative effects. Despite several explanations have been proposed [19, 24, 25,29, 33, 81, 82, 92], the most likely hypothesis lies in the end-to-end interactions of Tmunits [35, 66, 88].An increase of SL has a two-fold effect on the steady-state tension (see Fig. 2):the plateau force (i.e. the tension associated with saturating calcium) increases andthe calcium-sensitivity is enhanced (i.e. the sigmoidal curves are left-ward shifted).Whereas the explanation of the first effect is commonly-agreed to be linked to theincrease of extension of the single-overlap zone (i.e. the region of the sarcomere wherethe MF filament a single AF), a well-assessed explanation for the the second effect(known as length-dependent activation, LDA) has not yet been found [1, 25, 28, 63,67, 81, 88, 94, 97].The earliest attempts to model the calcium-driven regulation of the muscular con-tractile system date back to the 1990s [23, 58, 75, 82, 86, 105]. Those models rely4 ncreasing SL (a) (b) (c) Figure 2: Representation of the steady-state force-calcium relationship (a), the force-velocity relationship (b) and tension-elongation curves after a fast transient (c).on the formalism of continuous-time Markov Chains (CTMC), also known as MarkovJump processes (see e.g. [64]), to model the transitions between the different config-urations assumed by the proteins involved in the force regulation process. In thosemodels, the necessity of representing the end-to-end interactions of Tm units dictatesa spatially-explicit representation of the RUs. Indeed, mean-field models, where onlya single representative RU is considered, fail to correctly predict the cooperative acti-vation and the resulting steep force-calcium curves [81, 84]. However, the number ofdegrees of freedom of the CTMC increases exponentially with the number of RUs rep-resented. This hinders the possibility of numerically approximating the solution of theForward Kolmogorov Equation, also known as Master Equation in natural sciences,which describes the time evolution of the probabilities associated with the states of astochastic process [5]. As a matter of fact, the Forward Kolmogorov Equation asso-ciated with this CTMC would have a number of variables that is exponential in thenumber of RUs, resulting in as many as 10 or more variables [76]. For this reason,spatially-explicit models require a MC approximation for their numerical resolution,thus resulting in very large computational costs [76, 81, 100]. Active force is generated by XBs by the cyclical attachment and detachment of myosinheads (MHs) to actin binding sites (BSs). When MHs are in their attached configura-tion, they rotate towards the center of the sarcomere, performing the so-called power-stroke, thus pulling the AF along the same direction. Such cyclical path, known asLymn-Taylor cycle [59], features a wide range of time scales (nearly from 1 to 100 ms)[7, 14, 48]; hence, also the response of the force generation apparatus to externalstimuli is characterized by different time scales. Indeed, when a fast step in force(respectively, in length) is applied to an isometrically contracted muscle fiber, threedifferent phases can be observed [11, 14, 48, 60, 61]. First, an instantaneous elasticresponse occurs along the so-called T - L curve (see Fig. 2), whose slope correspondsto the stiffness of the attached myosin proteins. Then, a fast transient (2 − T - L curve (see Fig. 2). Such second phase corresponds to a rearrangement of MHswithin their pre- and post-power-stroke configuration. Finally, with a time scale onnearly 100 ms, the muscle fiber reaches a steady-state regime, characterized (in casethe step is applied by controlling the force) by a constant shortening (or lengthening)velocity. The relationship between the steady-state velocity and the muscle tensionconstitutes the so-called force-velocity curve (see Fig. 2), firstly measured by Hill in[36], and it is characterized by a finite value of velocity (denoted by v max ) for whichthe generated tension is zero [7, 47]. The experimental measurements of the T - L ,5he T - L and the force-velocity curves are invariant after normalization of the tensionby its isometric value, denoted by T isoa . This reveals that the underlying mechanismsare related to the XB dynamics, while they are independent of the thin filament reg-ulation, whose effect is simply that of tuning the number of recruitable XBs [7, 11,47].The attachment-detachment process of MHs has been described accordingly withthe formalism of the Huxley model [42] (that we denote by H57 model), where thepopulation of MHs is described by the distribution density of the distortion of attachedXB. The time evolution of such distribution is driven by a PDE, where a convectionterm accounts for the mutual sliding between filaments, and a source and a sink term(whose rates depend on the XB distortion) account for the creation and destructionof XBs [16, 40, 53, 54]. In order to capture the separation between the fastest timescales (i.e. between the first two phases following a fast step either in force or inlength), an explicit representation of the power-stroke must be included in the model,by introducing a multistable discrete [41, 89] or continuous [13, 14, 52, 60, 61] degreeof freedom, representing the angular position of the rotating MH. In the past two decades, several models describing the generation of active force inthe cardiac tissue, including both the calcium-driven regulation and the XB cycling,have been proposed. The main challenge faced in the development of such models liesin the spatial dependence of the cooperativity phenomenon, crucial to reproduce thecalcium dependence of muscle activation (see Sec. 2.1). As a matter of fact, an explicitrepresentation of spatial-dependent cooperative mechanisms dramatically increasesthe computational complexity of activation models, even more so when such modelsare coupled with models described XB cycling. When the interactions between BSsand MHs is considered, indeed, one must face the further difficulty of tracking whichBS faces which MH when the filaments mutually slide. The attempt of capturing suchspatially dependent phenomena in a compact system of ODEs is the common threadof most of the literature on sarcomere modeling (see e.g. [8, 10, 16, 55, 76, 83–85,100, 104]). We remark that computational efficiency is a major issue when sarcomeremodels are employed in multiscale simulation, such as cardiac EM. In this case, indeed,the microscale force generation model needs to be simultaneously solved in many pointsof the computational domain (typically at each nodal point of the computationalmesh). Nonetheless, most of biophysically detailed full-sarcomere models rely on thetime-consuming MC method for their numerical approximation [39, 91, 100–102].
In this section, we propose four different microscale models of active force generationin the cardiac tissue. These models are derived from a biophysically detailed CTMC,accurately describing the dynamics of the regulatory and contractile proteins. Wepresent a strategy to derive a set of equations, with a dramatically smaller numberof variables than the Forward Kolmogorov Equation, describing the evolution of thebiophysically detailed CTMC. Our strategy is based on physically motivated assump-tions, aimed at neglecting second-order interactions among the proteins, focusing onlyon the main interactions, so that the variables describing the stochastic processes as-sociated with the states of the proteins can be partially decoupled, similarly to what6e have done in [76]. This results in a drastic reduction in the size of models. More-over, we change the classical MH-centered description of XBs (see e.g. [13, 40, 41,89]), in favor of a BS-centered one. This prevents the necessity of tracking the mutualsliding between the filaments, still without the need of introducing any simplifyingassumption.As in most of RUs models (see Sec. 2.1), we describe Tn and Tm by discrete states.Moreover, based on the experimental evidence that cooperativity is linked to RUs end-to-end interactions [35, 66, 88], we include nearest-neighbor interactions among RUswith the formalism of the model of [84].Concerning the modeling of XBs, we are here interested in developing a model ofcardiomyocytes contraction in the heart, which is characterized by slower time-scalesthan the fast time-scale of the power-stroke. This suggests that the level of detailthat best suits the application to cardiac EM does not require to explicitly representthe power-stroke [77]. In [13], indeed, the authors showed that, if the time-scalesof interest are slower than the time-scale of the power-stroke, the detailed modelsincluding the power-stroke reduce to H57-like models, where only the attachment-detachment process of XBs is explicitly represented. Therefore, we model the XBdynamics as a two-states process, within the H57 formalism, where the attachment-detachment rates depend on the myosin arm distortion.
We focus on half MF, where N M MHs are located, and one AF, regulated by N A RUs(see Fig. 3). We place a reference system at the end of the AF closer to the center ofthe sarcomere and we consider the following smoothed geometrical factors, for i ∈ I A : χ M ( SL, i ) = 12 tanh (cid:18) y i − y LM ε (cid:19) + 12 tanh (cid:18) − y i − y RM ε (cid:19) ,χ SF ( SL, i ) = 12 (cid:18) (cid:18) y i − y SF ε (cid:19)(cid:19) , where y LM = (2 L A − SL + L H ) / , y SF = 2 L A − SL,y RM = (2 L A − SL + L M ) / , y i = L A N A ( i − . . The geometrical factors are defined so that we have χ M ( SL, i ) (cid:39) i -th RU facesthe considered half MF and χ SF ( SL, i ) (cid:39) i -th RU is in the single filamentregion (no overlap with other AFs occurs). Let d ij ( t ) be the distance between the i -thactin BS and the j -th MH at time t . For some constant d we have: d ij ( t ) = D A i − D M j + SL ( t )2 − d . (2)We denote by v hs ( t ) = − ddt SL ( t ) / v ( t ) = 2 v hs ( t ) /SL the normalized shortening velocity, where SL denotes a referencesarcomere length.In the model that we propose, each RU is characterized by the state of Tn (boundto calcium or not) and Tm (permissive or not). Moreover, each RU corresponds to a7 𝑆𝐹 𝑦 𝐿𝑀 𝝌 𝑴 𝝌 𝑺𝑭 𝑦 𝑦 … 𝑦 𝑅𝑀 𝐿𝑀 𝐿𝐴 𝑆𝐿(𝑡)𝐿𝐻 actin filament half myosin filament RU MH Figure 3: Sketch of the sarcomere model. The thick filament (MF) is represented inred and two thin filaments (AF) are represented in blue (top). The origin of the frameof reference is the left side of the reference AF. The functions χ SF and χ M are alsorepresented (bottom).BS to which the MHs can bind. Obviously, each MH can bind to a single BS at a timeand, conversely, each BS can have at most a single MH attached. We consider thereforethe following stochastic processes, for i ∈ I A := { , . . . , N A } , j ∈ I M := { , . . . , N M } and t ≥ C ti = (cid:40) B if the i -th Tn is bound to calcium, U else; T ti = (cid:40) P if the i -th Tm is permissive, N else; A ti = (cid:40) j if the i -th actin BS is attached to the j -th MH,0 if the i -th actin BS is not attached to any MH; M tj = (cid:40) i if the j -th MH is attached to the i -th actin BS,0 if the j -th MH is not attached to any actin BS; ,Z ti = (cid:40) x if the i -th actin BS is attached to a MH with displacement x, ∅ if the i -th actin BS is not attached to any MH. (3)We remark that we denote the detached state by Z ti = ∅ rather than Z ti = 0, becausethe latter notation is employed to denote the case when the i -th actin BS is attachedwith displacement x = 0. We notice that the last three processes are redundant, aswe have, for any i , j and t :( A ti = j ) ⇐⇒ ( M tj = i ) ⇐⇒ ( Z ti = d ij ( t )) . In the following, we will use the notation N = P , P = N , U = B and B = U to denoteopposite states.We assume, coherently with the hypotheses of the model in [84], that the dynamicsof Tn is affected by the state of the corresponding Tm, and that the dynamics of Tm8s affected by the Tn belonging to the same RUs and by the state of Tm in the nearestneighboring RUs, by a cooperativity mechanism. We exclude any feedback from XBson the dynamics of the RUs, as recent experimental evidence supports that this kindof feedback is not present [25, 92]. Conversely, the attachment-detachment rates of theXBs may depend on the distance between the MH and the BS, the relative velocitybetween the filaments and the permissivity state of the RU. Moreover, transition ratescan be affected by the mutual arrangements between the filaments: we allow thus fora dependence of the transition rates on the index corresponding to each unit. Weconsider therefore the following transition rates for the stochastic processes in (3),where δ ∈ {B , U} , α, β, η ∈ {P , N } : k δδ | βC,i = lim ∆ t → t P (cid:2) C t +∆ ti = δ | ( C i , T i ) t = ( δ, β ) (cid:3) ,k ββ | α · η,δT,i = lim ∆ t → t P (cid:2) T t +∆ ti = β | ( T i − , T i , T i +1 , C i ) t = ( α, β, η, δ ) (cid:3) ,f iα ( x, v ( t )) = lim ∆ t → t P (cid:2) Z t +∆ ti = x | Z ti = ∅ , T ti = α, ∃ j ∈ I M : d tij = x + v hs ∆ t, M tj = 0 (cid:3) ,g iα ( x, v ( t )) = lim ∆ t → t P (cid:2) Z t +∆ ti = ∅ | Z ti = x, T ti = α (cid:3) , (4)where v ( t ) plays the role of independent variable. In the definition of f iα , the eventsconditioning the probability ensure that, at time t , the i -th BS is not attached andthat there exists a non-attached MH at distance x + v hs ∆ t (so that at time t + ∆ t thedistance is reduced to x ).As mentioned earlier, the transition rates may be affected by the mutual arrange-ment of the filaments. Specifically, we assume that binding is possible only in thesingle-overlap region and that it may be characterized by different rates inside andoutside the single-overlap region. Hence, we have, for α ∈ {N , P} and i ∈ I A : f iα ( x, v ) = f α ( x, v ) χ M ( SL, i ) χ SF ( SL, i ) ,g iα ( x, v ) = g α ( x, v ) χ M ( SL, i ) χ SF ( SL, i ) + ˜ g α ( x, v ) (1 − χ M ( SL, i ) χ SF ( SL, i )) . (5)For the sake of simplicity, we shorten the notations as: f ij P = f i P ( d ij ( t ) , v ( t )) , f ij N = f i N ( d ij ( t ) , v ( t )) ,g ij P = g i P ( d ij ( t ) , v ( t )) , g ij N = g i N ( d ij ( t ) , v ( t )) . (6)Finally, the expected value of the force exerted by the considered pair of interactinghalf MF and AF is given by: F hf ( t ) = (cid:88) i ∈I A E (cid:2) F XB ( Z ti ) (cid:3) , where F XB ( x ) denotes the force exerted by an attached MH with distortion x andwhere we set by convention F XB ( ∅ ) = 0. In this section, we present several equations describing the evolution of the stochasticprocesses of Eq. (3). Such equations are rigorously derived from a number of physicallymotivated assumptions. For the sake of clarity, however, we first present the modelequations and we later list the underlying assumption in Sec. 3.3. A detailed derivationof the equations, with the mathematical details, is provided in App. A.9 .2.1 Thin filament regulation
Due to the lack of feedback from XBs to RUs, it is possible to write an equationdescribing the evolution of the stochastic processes C ti and T ti independently of thestochastic processes associated with the XBs. Similarly to what we have done in [76],we focus on the joint probabilities of triplets of consecutive RUs. Hence, we considerthe following collection of probabilities, where i = 2 , . . . , N A − ϑ, η, λ ∈ {U , B} and α, β, δ ∈ {N , P} : π αβδ,ϑηλi ( t ) := P (cid:2) ( T i − , T i , T i +1 ) t = ( α, β, δ ) , ( C i − , C i , C i +1 ) t = ( ϑ, η, λ ) (cid:3) . (7)We show in App. A.1 that, under suitable assumptions, the time evolution of suchvariables is described by the following system of ODEs, for t ≥ i =2 , . . . , N A − ϑ, η, λ ∈ {U , B} and α, β, δ ∈ {N , P} : ddt π αβδ,ϑηλi = (cid:101) k αα |◦ · β, ◦ ϑηT,i − π αβδ,ϑηλi − (cid:101) k αα |◦ · β, ◦ ϑηT,i − π αβδ,ϑηλi + k ββ | α · δ,ηT,i π αβδ,ϑηλi − k ββ | α · δ,ηT,i π αβδ,ϑηλi + (cid:101) k δδ | β · ◦ ,ηλ ◦ T,i +1 π αβδ,ϑηλi − (cid:101) k δδ | β · ◦ ,ηλ ◦ T,i +1 π αβδ,ϑηλi + k ϑϑ | αC,i − π αβδ,ϑηλi − k ϑϑ | αC,i − π αβδ,ϑηλi + k ηη | βC,i π αβδ,ϑηλi − k ηη | βC,i π αβδ,ϑηλi + k λλ | δC,i +1 π αβδ,ϑηλi − k λλ | δC,i +1 π αβδ,ϑηλi , (8)endowed with suitable initial conditions and where we have defined the followingtransition rates (in which the symbol ◦ denotes that the corresponding unit can be inan arbitrary state): (cid:101) k αα |◦ · β, ◦ ϑηT,i := (cid:80) ξ,ζ k αα | ξ · β,ϑT,i π ξαβ,ζϑηi (cid:80) ξ,ζ π ξαβ,ζϑηi for i = 2 , . . . , N A − ,k αα |N · β,ϑT,i for i = 1; (cid:101) k δδ | β · ◦ ,ηλ ◦ T,i := (cid:80) ξ,ζ k δδ | β · ξ,λT,i π βδξ,ηλζi (cid:80) ξ,ζ π βδξ,ηλζi for i = 2 , . . . , N A − ,k δδ | β · N ,λT,i for i = N A . (9)The permissivity of a given RU is defined as its probability of being in permissivestate (i.e. P i ( t ) = P [ T ti = P ]) and can be obtained from the variables π αβδ,ϑηλi ( t ) as: P i ( t ) = (cid:88) β,δ ∈{N , P} (cid:88) ϑ,η,λ ∈{U , B} π P βδ,ϑηλ ( t ) for i = 1 , (cid:88) α,δ ∈{N , P} (cid:88) ϑ,η,λ ∈{U , B} π α P δ,ϑηλi ( t ) for i = 2 , . . . , N A − , (cid:88) α,β ∈{N , P} (cid:88) ϑ,η,λ ∈{U , B} π αβ P ,ϑηλN A − ( t ) for i = N A . .2.2 Crossbridge dynamics In order to account for the imperfections in the sarcomere lattice, we consider the valueof d (see Eq. (2)) as a random variable rather than a constant. Hence, we assume that,given a BS in front of the MF, the probability that the closest MH is located at distance x is uniform for x ∈ [0 , D M ). We denote by ρ M := f (cid:2) ∃ j ∈ Z : d tij = x (cid:3) = D M − theMH linear density, where the symbol f denotes the probability density function, thatis: P (cid:2) ∃ j ∈ Z : d tij = x ∈ ( a, b ) (cid:3) = (cid:90) ba f (cid:2) ∃ j ∈ Z : d tij = x (cid:3) dx. We define the following variables, for i ∈ I A , corresponding to the probability densitythat the i -th BS is attached to a MH with displacement x and that the associated RUis in a given permissivity state: n i, P ( x, t ) = f (cid:2) Z ti = x, T ti = P (cid:3) ,n i, N ( x, t ) = f (cid:2) Z ti = x, T ti = N (cid:3) . (10)We notice that we make here the choice of tracking the XBs from the point of view ofthe BSs, rather than of the MHs, as it is traditionally done in literature [13, 40, 41,89]. This change of perspective has the significant advantage that it does not requireto track which RU faces which MH at each time. Indeed, each BS and each RU, beinglocated on the same filament, rigidly move with respect of each others and, thus, eachBS is always associated with the same RU.In App. A.3 we derive the following system of PDEs, describing the time evolutionof n i, P and n i, P : ∂n i, P ∂t − v hs ∂n i, P ∂x = ( ρ M P i − n i, P ) f i P − g i P n i, P − (cid:101) k PN T,i n i, P + (cid:101) k N P
T,i n i, N x ∈ R , t ≥ , i ∈ I A ,∂n i, N ∂t − v hs ∂n i, N ∂x = ( ρ M (1 − P i ) − n i, N ) f i N − g i N n i, N − (cid:101) k N P
T,i n i, N + (cid:101) k PN T,i n i, P x ∈ R , t ≥ , i ∈ I A , (11)endowed with suitable initial conditions, where we define: (cid:101) k N P
T,i := (cid:80) α,δ,ϑ,η,λ k N P| α · δ,ηT,i π α N δ,ϑηλi − P i , (cid:101) k PN T,i := (cid:80) α,δ,ϑ,η,λ k PN | α · δ,ηT,i π α P δ,ϑηλi P i . (12)We notice that Eq. (11) can be interpreted as a generalized H57 model, where theXB population is split into two groups, according to the permissivity state of thecorresponding RU. The sink and source terms in Eq. (11) account for the fluxes amongthe two populations. The terms P i and (1 − P i ) represent the maximum possibleproportion of attached BSs in each group, while the factor ρ M appears because inthis setting n i, P and n i, N are, from a dimensional point of view, the inverse of lengthunits (they are probability densities), whereas the variables of the H57 model aredimensionless. We remark that, differently than the H57 model, Eq. (11) is referredto BSs rather than to MHs. 11y combining Eq. (8) with Eq. (11), describing the dynamics of RUs and XBs,respectively, we obtain a model that we denote as the SE-PDE model (where SE stands for spatially-explicit, while PDE denotes the fact that the XB dynamics isdescribed by a PDE system). The expected value of the force exerted by the wholehalf filament can be determined as follows: F hf ( t ) = (cid:88) i ∈I A (cid:90) + ∞−∞ F XB ( x ) ( n i, P ( x, t ) + n i, N ( x, t )) dx. (13) When the XB attachment-detachment transition rates assume special forms, the PDEsystem of Eq. (11) can be reduced to a more compact system of ODEs, by followinga general strategy in statistical physics, already used for H57-like models [8, 16, 104].Specifically, let us assume that the total attachment-detachment rate is independentof the XB distortion (i.e. there exist functions r i P ( v ) and r i N ( v ), for i ∈ I A , such that r i P ( v ) = f i P ( x, v ) + g i P ( x, v ) and r i N ( v ) = f i N ( x, v ) + g i N ( x, v ) for any x ∈ R ). Then,we define, for α ∈ {N , P} , for ψ ∈ { f i P , f i N , g i P , g i N } and for p ∈ N : µ pi,α ( t ) := (cid:90) + ∞−∞ (cid:18) xSL / (cid:19) p n i,α ( x, t ) dx,µ pψ ( v ) := (cid:90) + ∞−∞ (cid:18) xSL / (cid:19) p ψ ( x, v ) dxD M . (14)We notice that µ i, N ( t ) (respectively, µ i, P ( t )) can be interpreted as the probability thatthe i -th BS is attached and associated to a non-permissive (respectively, permissive)RU. Moreover, µ i, N ( t ) /µ i, N ( t ) (respectively, µ i, P ( t ) /µ i, P ( t )) corresponds to the ex-pected value of the distortion (normalized with respect to SL /
2) of the XB attachedto the i -th RU, provided that the corresponding RU is in non-permissive (respectively,permissive) state.Under these assumptions, as shown in App. A.4, we get the following distribution-moments equations: ddt µ i, P = − (cid:16) r i P + (cid:101) k PN T,i (cid:17) µ i, P + (cid:101) k N P
T,i µ i, N + P i µ f i P t ≥ , i ∈ I A ,ddt µ i, N = − (cid:16) r i N + (cid:101) k N P
T,i (cid:17) µ i, N + (cid:101) k PN T,i µ i, P + (1 − P i ) µ f i N t ≥ , i ∈ I A ,ddt µ i, P + v µ i, P = − (cid:16) r i P + (cid:101) k PN T,i (cid:17) µ i, P + (cid:101) k N P
T,i µ i, N + P i µ f i P t ≥ , i ∈ I A ,ddt µ i, N + v µ i, N = − (cid:16) r i N + (cid:101) k N P
T,i (cid:17) µ i, N + (cid:101) k PN T,i µ i, P + (1 − P i ) µ f i N t ≥ , i ∈ I A . (15)We notice that, thanks to (5), we have: µ pf iα = µ pf α χ M ( SL, i ) χ SF ( SL, i ) ,µ pg iα = µ pg α χ M ( SL, i ) χ SF ( SL, i ) + µ p ˜ g α (1 − χ M ( SL, i ) χ SF ( SL, i )) . (16)By assuming a linear spring hypothesis for the tension generated by attached XBs(i.e. F XB ( x ) = k XB x ), the expected value of the force of half filament is given by: F hf ( t ) = k XB SL N A µ ( t ) , (17)12here we have defined: µ p ( t ) = 1 N A N A (cid:88) i =1 (cid:16) µ pi, P + µ pi, N (cid:17) . (18)Since the active force generated by half MF is proportional to µ ( t ), there exists someconstant a XB (with the dimension of a pressure) such that the macroscopic activetension can be written as T a ( t ) = a XB µ ( t ). In the following, we denote by SE-ODEmodel the one obtained by combining Eq. (8) with Eq. (15).
In the previous sections, we proposed a model of tension generation in the muscle tissuebased on an explicit spatial description of the physical arrangements of proteins alongthe myofilaments. The spatial description allows to model the cooperativity mech-anism (linked to the nearest-neighbor interactions within RUs) and the SL relatedeffects on the force generation machinery (linked to the filament overlapping). How-ever, the first phenomenon, despite being spatially dependent, is based on interactionsof local type; the effect of the second phenomenon, in turn, largely depends on thesize of the single-overlap zone, that is a scalar quantity non dependent on the spatialvariable. Based on the above considerations, we propose a mean-field approximationof the spatially dependent CTMC presented in Sec. 3.1, where the nearest-neighboringinteraction are captured as a local effect, and the effect of SL is modeled in functionof the size of the single-overlap zone.This mean-field model is based on the assumption that the single-overlap zone canbe considered as infinitely long. Such approximation is reasonable as far as the effectof the edges can be neglected (the validity of such approximation will be discussedin Sec. 6). A direct consequence of this assumption is the invariance by translationof the joint distribution of RUs. This means that for any set of indices I ⊂ Z and I ⊂ Z and for any collection of states α i ∈ {N , P} (for i ∈ I ) and β i ∈ {U , B} (for i ∈ I ), the joint distribution of the states of the corresponding RUs is not affected ifthe RUs are translated by a count of k ∈ Z units: P (cid:34)(cid:32) (cid:92) i ∈I T ti = α i (cid:33) ∩ (cid:32) (cid:92) i ∈I C ti = β i (cid:33)(cid:35) = P (cid:34)(cid:32) (cid:92) i ∈I T ti + k = α i (cid:33) ∩ (cid:32) (cid:92) i ∈I C ti + k = β i (cid:33)(cid:35) . It follows that the variables π αβδ,ϑηλi ( t ) defined in Eq. (7) coincide for each i . Inaddition, we further reduce the number of variables by tracking only the state of theTn of the central RU of the triplet (this further reduction is made possible by the factthat we never have to track the behavior at the boundaries of the filaments, as we willsee in what follows). We define thus the following variables, for α, β, δ ∈ {N , P} and η ∈ {U , B} : π αβδ,η ( t ) := P (cid:2) ( T i − , T i , T i +1 ) t = ( α, β, δ ) , C ti = η (cid:3) . (19)We notice that the variables π αβδ,η ( t ) are well-defined thanks to the translational in-variance of the distribution of RUs. Moreover, the transition rates k δδ | βC,i and k ββ | α · η,δT,i for the units in the single-overlap region do not depend on i . Hence, we will denotethem simply as k δδ | βC and k ββ | α · η,δT .The time evolution of the variables π αβδ,η ( t ) is described, under suitable assump-tions (as proved in App. A.2), by the following ODE model, valid for t ≥ α, β, δ ∈ {N , P} and η ∈ {U , B} : ddt π αβδ,η = (cid:101) k αα |◦ · β, ◦ T π αβδ,η − (cid:101) k αα |◦ · β, ◦ T π αβδ,η + k ββ | α · δ,ηT π αβδ,η − k ββ | α · δ,ηT π αβδ,η + (cid:101) k δδ | β · ◦ , ◦ T π αβδ,η − (cid:101) k δδ | β · ◦ , ◦ T π αβδ,η + k ηη | βC π αβδ,η − k ηη | βC π αβδ,η , (20)where: (cid:101) k αα |◦ · β, ◦ T := (cid:80) ξ,ζ k αα | ξ · β,ζT π ξαβ,ζ (cid:80) ξ,ζ π ξαβ,ζ , (cid:101) k δδ | β · ◦ , ◦ T := (cid:80) ξ,ζ k δδ | β · ξ,ζT π βδξ,ζ (cid:80) ξ,ζ π βδξ,ζ . The permissivity of a RU in the single-overlap zone, defined as P ( t ) = P [ T ti = P ](such that the i -th RU belongs to the single-overlap zone), can be obtained as: P ( t ) = (cid:88) α,δ ∈{N , P} (cid:88) η ∈{U , B} π α P δ,η ( t ) . By similar arguments, it follows that also the joint distribution of the stochastic pro-cesses associated with XB formation enjoy the translational invariance property and,consequently, the following variables are well defined, as the right-hand sides are in-dependent of the index i (for i belonging to the single-overlap zone): n P ( x, t ) = f (cid:2) Z ti = x, T ti = P (cid:3) ,n N ( x, t ) = f (cid:2) Z ti = x, T ti = N (cid:3) . (21)By proceeding as in Sec. 3.2.2, we get the following model: ∂n P ∂t − v hs ∂n P ∂x = ( ρ M P − n P ) f P − g P n P − (cid:101) k PN T n P + (cid:101) k N P T n N x ∈ R , t ≥ ,∂n N ∂t − v hs ∂n N ∂x = ( ρ M (1 − P ) − n N ) f N − g N n N − (cid:101) k N P T n N + (cid:101) k PN T n P x ∈ R , t ≥ , (22)where we have defined: (cid:101) k N P T ( t ) := (cid:80) α,δ,η k N P| α · δ,ηT π α N δ,η ( t )1 − P ( t ) , (cid:101) k PN T ( t ) := (cid:80) α,δ,η k PN | α · δ,ηT π α P δ,η ( t ) P ( t ) . (23)The expected value of the force exerted by the whole half filament can be obtained asfollows: F hf ( t ) = N A χ so ( SL ( t )) (cid:90) + ∞−∞ F XB ( x ) ( n P ( x, t ) + n N ( x, t )) dx, (24)14here the single-overlap ratio χ so denotes the fraction of the AF filament in the single-overlap zone: χ so ( SL ) := SL ≤ L A , SL − L A ) L M − L H if L A < SL ≤ L M ,SL + L M − L A L M − L H if L M < SL ≤ L A − L H , L A − L H < SL ≤ L A + L H ,L M + 2 L A − SLL M − L H if 2 L A + L H < SL ≤ L A + L M , SL > L A + L M . (25)We notice that we are here assuming that the relative sliding between the filamentsis such that χ so slowly varies, so that we can neglect the effects linked to the statetransitions taking place at the border of the single-overlap zone. The combination ofEqs. (20) and (22) gives a model for the full-sarcomere dynamics, that we denote asthe MF-PDE model (where MF stands for mean-field).Moreover, under the assumption that the total attachment-detachment rate doesnot depend on the XB elongation (i.e. there exist two functions r P ( v ) and r N ( v ) suchthat r P ( v ) = f P ( x, v ) + g P ( x, v ) and r N ( v ) = f N ( x, v ) + g N ( x, v ) for any x ∈ R ), wecan derive the following distribution-moment equation: ddt µ P = − (cid:16) r P ( v ) + (cid:101) k PN T (cid:17) µ P + (cid:101) k N P T µ N + P µ f P t ≥ ,ddt µ N = − (cid:16) r N ( v ) + (cid:101) k N P T (cid:17) µ N + (cid:101) k PN T µ P + (1 − P ) µ f N t ≥ ,ddt µ P + v µ P = − (cid:16) r P ( v ) + (cid:101) k PN T (cid:17) µ P + (cid:101) k N P T µ N + P µ f P t ≥ ,ddt µ N + v µ N = − (cid:16) r N ( v ) + (cid:101) k N P T (cid:17) µ N + (cid:101) k PN T µ P + (1 − P ) µ f N t ≥ , (26)endowed with suitable initial conditions and where we define, for α ∈ {N , P} : µ pα ( t ) := (cid:90) + ∞−∞ (cid:18) xSL / (cid:19) p n α ( x, t ) dx. (27)The force exerted by half thick filament is then given by: F hf ( t ) = k XB SL N A µ ( t ) , (28)where µ p ( t ) := χ so ( SL ( t )) [ µ p P ( t ) + µ p N ( t )] , (29)for p = 0 , T a ( t ) = a XB µ ( t ). Finally, bycombining Eqs. (20) and (26) we obtain a model that we denote as MF-ODE model . The models presented in Sec. 3.2 are derived based on a few assumptions, that arediscussed in this section. Let us introduce the following assumptions:15I1) ( T i +1 , C i ) t ⊥⊥ ( T i − , C i − ) t | ( T i − , T i ) t for i = 3 , . . . , N A − T i − , C i ) t ⊥⊥ ( T i +2 , C i +1 ) t | ( T i +1 , T i ) t for i = 2 , . . . , N A − T i +1 , C i +1 ) t ⊥⊥ T ti − | ( T i − , T i , C i − , C i ) t for i = 3 , . . . , N A − T i − , C i − ) t ⊥⊥ T ti +2 | ( T i +1 , T i , C i +1 , C i ) t for i = 2 , . . . , N A − A ti ⊥⊥ ( T i − , T i +1 , C i ) t | T ti for i = 2 , . . . , N A − A, B, C we write A ⊥⊥ B | C if A and B are conditionallyindependent given C (i.e. P [ A ∩ B | C ] = P [ A | C ] P [ B | C ]).Assumptions (I1) and (I2) state that distant RUs are conditionally independentgiven the states of the intermediate RUs. From a modeling viewpoint, this means thatthe interaction of far units is mediated by the intermediate ones, which is coherentwith the short-range nature of end-to-end interactions. Assumption (I1) is used for thederivation of the spatially-explicit models, while assumption (I2) is used for mean-fieldmodels.Assumption (I3), on the other hand, entails that the state of a BS is conditionallyindependent of the state of surrounding RUs, given the permissivity state of the asso-ciated RU. This is also coherent with the assumptions underlying the model, as theonly feature of the RUs that directly affects the XBs binding rates is the permissivitystate of Tm.Moreover, we consider the following dual assumptions:(A1) f ij P (cid:54) = 0 = ⇒ A h (cid:54) = j ∀ h (cid:54) = i ;(A2) f ij P (cid:54) = 0 = ⇒ M k (cid:54) = i ∀ k (cid:54) = j .Assumption (A1) states that, whenever a XB can form, the MH cannot be involvedin a XB with a farther BS. Suppose that the support of f is contained in the interval[ x , x + h ]. Then, this is equivalent to say that, if d ij ∈ [ x , x + h ], the XBs betweenthe couples ( i − , j ) and ( i + 1 , j ), which feature displacements d ij − D A and d ij + D A respectively, cannot exist. This condition is automatically fulfilled if XBs are presentonly for displacements in the interval ( − D A + x + h, D A + x ), which has width2 D A − h . The interval consists in the support of f , with width h , surrounded by twobands of width D A − h .On the other hand, assumption (A2) states that, whenever a XB can form, theBS cannot be involved in a XB with a farther MH. By similar considerations, it turnsout that this hypothesis is satisfied if XBs are present only in the range ( − D M + x + h, D M + x ). Since D M > D A , assumption (A1) is stronger than (A2).Assumptions (A1)-(A2) allow to decouple the dynamics of the different units. Theirvalidity is justified when the shortening velocity is relatively small, whereas, for largevelocities, the XB displacements may be convected outside the region ( − D A + x + h, D A + x ). Figure 4 provides a visual representation of assumptions (A1)-(A2). Wenotice that all the models belonging to the family of the H57 model are based on theseassumptions, without which the H57 equation cannot be derived.In App. A we report the detailed derivation of the models proposed in this paper,based on the above presented assumptions. Table 2 provides a recap of the differentmodels, together with the assumptions underlying each model. The two models de-rived by the distribution-moments (SE-ODE and MF-ODE) also require that the sumof the attachment and detachment rates is independent of x (we write f + g ⊥⊥ x ).We notice that this is not a simplificatory assumption, but rather a specific modelingchoice that allows to write the model with the distribution-moments formalism.16 🗸 ✗ AFMF ✗ 🗸 ✗ AFMF
Assumption (A1) Assumption (A2)
Figure 4: Representation of assumptions (A1)-(A2). According to assumption (A1)(respectively, assumption (A2)) when a BS-MH pair is within the XB formation range,then the adjacent BSs (respectively, MHs) cannot be bound to the considered MH(respectively, BS).
Model name (Equations) Number of ODEsNumber of PDEs Assumptions Modelingchoices
SE-PDE (8)-(11) ( N A − N A = 1280= 64 (I2),(I3),(A1) SE-ODE (8)-(15) ( N A − + 4 N A - = 1408 (I2),(I3),(A1) f + g ⊥⊥ x MF-PDE (20)-(22) 2 MF-ODE (20)-(26) 2 + 4- = 20 (I1),(I3),(A1), m.f. f + g ⊥⊥ x Table 2: List of the models proposed in Sec. 3. For future reference, we assign a nameto each model ( SE stands for spatially-explicit, MF stands for mean-field). In thesecond column we report the number of ODEs and PDEs included in each model asa function of N A and N M and we specify the resulting values in the case N A = 32, N M = 18. In the “Assumptions” column, m.f. stands for mean-field assumption .17 .4 Transition rates definition In Sec. 3.2 we have shown that, under some physically motivated assumptions, theCTMC presented in Sec. 3.1 can be described by different systems of ODEs and/orPDEs. The models listed in Tab. 2 are thus valid independently of the specific choiceof the transition rates (with the only exception of the models SE-ODE and MF-ODEthat require that the sum of the detachment and attachment rate is independent ofthe XB distorsion). In this section, we present and motivate the specific choice oftransition rates that we will adopt in the rest of this paper.
The RU dynamics is determined by the eight rates associated with the forward andbackward transitions
UN − (cid:42)(cid:41) − BN , UP − (cid:42)(cid:41) − BP , UN − (cid:42)(cid:41) − UP and
BN − (cid:42)(cid:41) − BP . Thetransition rates are affected by [Ca ] i (that enhances in a multiplicative way thetransition U → B ), the filament overlap and the state of the nearest-neighboring Tmunits (for the latter interaction we adopt the cooperative interactions proposed in [84]).We start by considering the single-overlap zone, where we adopt the transition rates ofthe model of [84]. We show below that the transition rates of [84] are, however, rathergeneral, as they are based on just a couple of assumptions. We keep the notationconsistent with [84] to allow for comparisons.We call k BU|N C := k off and, without loss of generality, we set k BU|P C := k off /µ ,where the constant µ allows to differentiate the two rates. Experiments carried outwith protein isoforms from different species highlight that there is no apparent vari-ation in the transition U → B in different combinations of Tn subunits and Tm [62].We assume thus that the transition rates for
UN → BN and for
UP → BP coin-cide, and we set k UB|N C = k UB|P C := k off /k d [Ca ] i . Conversely, we allow the reversetransition rates to depend on the state of the associated Tm. Concerning the tran-sitions involving Tm, we assume that the calcium binding state of Tn affects thetransition rate of N → P for the associated Tm, but not the reverse rate. Therefore,we set k PN | α · δ, U T = k PN | α · δ, B T = k basic γ − n ( α,δ ) , where n ( α, δ ) ∈ { , , } denotesthe number of permissive states among α and δ , as proposed in [84]. Then, withoutloss of generality we denote k N P| α · δ, B T = Q k basic γ n ( α,δ ) , where the constant Q al-lows to differentiate the forward and backward transition rates ∗ . The only transitionrate left is given, to satisfy the detailed-balance consistency with the other rates, by k N P| α · δ, U T = Q/µ k basic γ n ( α,δ ) .We remark that, due to the definition of n ( α, δ ) as the number of neighboring unitsin permissive state, the units located at the ends of the filament cannot have n = 2.Coherently with this, in Eq. (9), the state of the missing neighboring RUs is set to N .In conclusion, the transition rates are determined by the five parameters Q , µ , k d , k off , k basic (plus the parameter γ that regulates the amount of cooperativity), resultingfrom the eight free parameters constrained by the two assumptions ( U → B not affectedby Tm,
P → N not affected by Tn) and by the detailed-balance consistency.Concerning the dependence on the filament overlap, we assume that the only tran-sition affected by filament overlap is
N → P , that is prevented in the central zone ofthe sarcomere, where the two AFs meet [32]. Specifically, we set, for η ∈ {U , B} and ∗ The physical interpretation of the constant γ corresponds to γ = exp(2 ∆ Ek B T ), where k B isthe Boltzmann constant, T the absolute temperature and ∆ E denotes the energetic gain of theconfiguration of neighboring units in the same state (i.e. N - N or P - P ) with respect to that withdifferent states (i.e. N - P or P - N ). See [79] for more details in this regard. 𝒩 ℬ𝒩ℬ𝒫𝒰𝒫 𝑘 𝑜𝑓𝑓 𝛾 − 𝑛 𝑘 𝑏 𝑎 𝑠 𝑖 𝑐 𝛾 − 𝑛 𝑘 𝑏 𝑎 𝑠 𝑖 𝑐 𝛾 𝑛 𝑄 𝑘 𝑏 𝑎 𝑠 𝑖 𝑐 𝜇 𝜒 𝑆 𝐹 ( 𝑆 𝐿 , 𝑖 ) 𝑘 𝑜𝑓𝑓 /𝜇 𝑘 𝑜𝑓𝑓 /𝑘 𝑑 [ 𝐶𝑎 ] 𝑖 𝑘 𝑜𝑓𝑓 /𝑘 𝑑 [ 𝐶𝑎 ] 𝑖 𝛾 𝑛 𝑄 𝑘 𝑏 𝑎 𝑠 𝑖 𝑐 𝜒 𝑆 𝐹 ( 𝑆 𝐿 , 𝑖 ) Figure 5: The proposed four states Markov model describing each RU. The termsdepending on the intracellular calcium concentration [Ca ] i are highlighted in red;terms depending on the state of neighbouring RUs (i.e. depending on n ) are high-lighted in blue; terms depending on the position of the RU and the current sarcomereelongation are highlighted in green.for α, δ ∈ {N , P} : k N P| α · δ,ηT,i = χ SF ( SL, i ) k N P| α · δ,ηT . (30)The resulting 4-states CTMC associated with each RU is represented in Fig. 5. On the basis of the results of [77], we work under the hypothesis that the totalattachment-detachment rate is independent of the XB elongation. In this case, themodels SE-ODE and MF-ODE can be used in place of the more computationallyexpansive counterparts SE-PDE and MF-PDE, which involve the solution of a PDEsystem. As a matter of fact, in [77] we have shown that the introduction of such hy-pothesis significantly reduces the number of parameters to be fitted by experiments,still preserving the capability of the models of reproducing a wide range of experi-mental characterizations. Moreover, we also make the reasonable assumption that thesliding velocity only affects the detachment rate.We assume that XB attachment is only possible when the corresponding Tm isin the state P (thus f N ≡ g N ≡ ˜ g P )and that the detachment rate when Tm is in state N is not affected by the filamentoverlap (i.e. ˜ g N ≡ g N ). In summary, if we assume that the total transition rate inthe single-overlap zone does not depend on the Tm state, we have: f P ( x, v ) + g P ( x, v ) = g N ( x, v ) = ˜ g N ( x, v ) = ˜ g P ( x, v ) = r + q ( v ) , for some constant r and function q , such that q (0) = 0. Moreover, since we focus onthe small velocity regimes, the function q is well characterized by its behavior around v = 0. Hence, we set q ( v ) = α | v | , in order to ease the calibration process. Therefore,the parameters to calibrate are µ f P , µ f P , r , α and a XB , to link the microscopic forcewith the macroscopic active tension. 19arameter Value Units Parameter Value Units SL . µ m ε . µ m L A . µ m N M
18 - L M . µ m N A
32 - L H . µ mTable 3: List of geometrical constants. To calibrate the models, we will restrict ourselves to experimental measurements fromthe literature coming from intact cardiac cells since the skinning procedure alters ina significant (and only partially understood) manner the activation and force genera-tion apparatus [4, 50, 96]. Moreover, thanks to the technique of flura-2 fluorescence,it is nowadays possible to measure the intracellular calcium concentration withoutdepriving the cell of its membrane, and it is also possible to inhibit the sarcoplas-mic reticulum calcium uptake by cyclopiazonic acid, so that the calcium level can becontrolled without the need of skinning the cells [95, 96].However, only few data are available from human cells, compared to animal cells,and most of experiments are performed at room temperature. Since the availabledata for human cells at body temperature are not sufficient to adequately constrainthe parameters of our models, we proceed as follows. First, we calibrate the modelparameters from rat experiments at room temperature and than we adjust the pa-rameters that are reasonably affected by the two varying factors (i.e. inter-speciesvariability and temperature) to fit the available data from human cells as body tem-perature. We compensate in this way for the data deficiency. We notice that we workunder the hypothesis that inter-species variability does not affect the fundamentalprocesses of tissue activation and force generation, but, since different species expressdifferent isoforms of the same protein, it can be encompassed by changing the pa-rameters of the same mathematical model (see [98] for a detailed discussion on thistopic).We report in Table 3 the geometrical constants describing the size of the myofila-ment components that we use in the following.
Even if the thin-filament activation precedes the XB cycling from a logical viewpoint,we start by illustrating the calibration procedure for the latter part. The reason willbe clarified later.In [77] we have shown that the parameters of the distribution-moments equationsdescribing the XB dynamics can be calibrated starting from five quantities, that aredescribed in the following. Under isometric conditions, we consider the isometrictension T isoa := a XB µ and the fraction of attached XBs µ := µ , where µ and µ are the steady-state solution for v = 0. The force-velocity relationship is characterizedby the maximum shortening velocity v max , and by v , its intersection with the axis T a = 0 of the tangent of the curve in isometric conditions, defined as: v := − (cid:18) ∂T a ( v ) /T isoa ∂v (cid:12)(cid:12)(cid:12)(cid:12) v =0 (cid:19) − , T a ( v ) denotes the steady-state tension for velocity v . Finally, as discussed in[77], the response to fast inputs is characterized, in the small-velocity regime, by thenormalized slope of the tension-elongation curve after a fast step. Such quantity isdefined as follows: by applying a step in length ∆ L to an isometrically contractedmuscle in a small time interval ∆ t , we define by T a (∆ t ) the tension recorded after thestep is applied, and we define the normalized stiffness as:˜ k := − ∂T a (∆ t ) /T isoa ∂ ∆ L (cid:12)(cid:12)(cid:12)(cid:12) ∆ L =0 . As discussed in [77], when the small-velocity regimes is considered, the quantity ˜ k corresponds to the slope of the T - L relationship.In conclusion, as shown in [77], by acting on the five parameters µ f P , µ f P , r , α and a XB we can fit experimentally measured values concerning the isometric solution( T isoa and µ ), the force-velocity relationship ( v max and v ) and the fast transientsresponse (˜ k ). All the above mentioned experimental setups are such that the thinfilament activation machinery can be considered in steady-state. Indeed, [Ca ] i isconstant in all the cases and, concerning SL : it is also constant under isometricconditions; constant shortening experiments are typically performed in the plateauregion of the force-length relationship, and thus the effect of changes in SL is irrelevant;fast transient experiments are carried out at a time-scale such that the activationvariables can be considered constant, since they are characterized by a much slowerdynamics [7, 47]. Therefore, in these cases, the values of P i , (cid:101) k PN T,i and (cid:101) k N P
T,i can beconsidered as fixed in Eq. (15) (and similarly in Eq. (26)). This observation is crucialsince it allows to decouple the calibration of the parameters involved in the thinfilament activation from the calibration of the parameters involved in XBs cycling.In fact, the five parameters of the XBs cycling model can be calibrated from the fiveabove mentioned experimental values.We notice that, while for the models considered in [77] the relationship betweenthe five parameters and the five experimentally measured values can be analyticallyinverted, in this case we find the values of the parameters with a numerical strategy.Specifically, to find the steady-state solution we solve Eq. (15) by setting to zero thetime derivative terms; we consider the exact solution after the fast transient (the linearODE system (15) can be solved analytically); we approximate the derivative appearingin the definition of v and ˜ k by finite differences. Finally, we solve, for the fiveparameters, the nonlinear system of equations linking the five measured values withthe parameters themselves. With this aim we employ the Newton-Raphson method,by approximating the Jacobian matrix by means of finite differences.Therefore, once the thin filament activation model (8) (or (20)) has been calibrated,we have at our disposal an automatic procedure to calibrate the remaining modelparameters. For this reason, we first setup such calibration procedure for the modelparameters associated with the XBs cycling (i.e. Eq. (15) or (26)) and, successively,we calibrate the model parameters associated with the RUs activation (i.e. Eq. (8)or (20)), so that we can directly see the effect of changes of such parameters on theresulting force since the remaining parameters are automatically adjusted.The experimental data used to calibrate the model are reported in Table 4, togetherwith a reference to the source in literature. As we mentioned at the beginning of Sec. 4,we employ data coming from room-temperature intact cardiac rat cell, apart from µ (acquired from skeletal frog muscle), for which we did not find measurements fromcardiac muscles. However, as shown in [77], this variable only affects the value of themicroscopic variables (i.e. µ pi,α ), but not that of the predicted active tension T a .21arameter Value Units Reference T isoa
120 kPa [95] µ .
22 - [9] v max − [11] v − [11]˜ k
66 - [11]Table 4: List of the experimental data used for model calibration.We notice that the constants v max , v and ˜ k are normalized with respect to T isoa and are thus valid for a wide range of activation levels (see Sec. 2.2). Conversely, thevalue of T isoa is associated with a SL in the plateau region and to saturating calciumconcentration. Therefore, when we calibrate the parameters we set [Ca ] i = 10 µ Mand SL = 2 . µ m. The steady-state solution of the activation models (8) and (20) only depends on theratio between the pairs of opposite transition rates (e.g. k N P| α · δ, B T /k PN | α · δ, B T = Q γ n ( α,δ ) − ). Therefore, the six parameters can be split into two groups: the firstgroup ( Q , µ , k d and γ ) determines the steady-state solution, while the second group( k off , k basic ) only affects the kinetics of the model (that is to say how fast the tran-sients are). This allows to calibrate first the parameters of the first group, and, onlysuccessively, those of the second group. The steady-state characterization of the muscle tissue activation is represented by theforce-calcium and force-length relationships (see Sec. 2.1), whose main features arethe behavior for SL in the plateau region (characterized by the tension for saturatingcalcium T isoa , the calcium sensitivity EC and the cooperativity coefficient n H ) andthe effect of SL (on the saturating tension T isoa and on the calcium sensitivity EC ).The tension for saturating calcium concentrations T isoa in the plateau region of SL is automatically fitted, thanks to Sec. 4.1. The effect of k d is that of shiftingthe force-calcium curves with respect to the log [Ca ] i axes since it only appears incombination with [Ca ] i in the model equations. Therefore, the value of k d can beeasily calibrated to match the experimental data as it only affects EC . The effectof γ , on the other hand, is that of tuning the amount of cooperativity (indeed, itacts on n H ). The role of the remaining parameters ( Q and µ ) is more involved andcannot be easily decoupled, as they affect the cooperativity, calcium sensitivity, andthe asymmetry of the force-calcium relationship below and above EC [84]. Moreover,in the SE-ODE case, they also act on the SL -driven regulation on calcium sensitivity.In the following we set µ = 10, coherently with the fact that the experimentallymeasured dissociation rate of Tn from calcium varies of one order of magnitude indifferent combinations [62]. For the SE-ODE model, we set γ , Q and k d , to fit thesteady-state force-calcium measurements of [95] (referred to the two different valuesof SL of 1 .
85 and 2 . µ m) from intact rat cardiac cells at room temperature. Theresulting steady-state curves are reported in Fig. 6. We are able to well fit the mainfeatures of the curves, including the characteristic S-shape, the plateau forces at both22 -1 [Ca ] i [ M] T a [ k P a ] SE-ODE model -1 [Ca ] i [ M] T a [ k P a ] MF-ODE model
SL = 1.85 m (model)SL = 2.15 m (model) SL = 1.85 m (exp)SL = 2.15 m (exp)
Figure 6: Steady-state force-calcium curves obtained with the SE-ODE (left) andMF-ODE (right) models with the optimal parameters values reported in Tab. 7 for SL = 1 . µ m and SL = 2 . µ m, compared to experimental data from intact ratcardiac cells at room temperature, from [95]. SL , the significant cooperativity typical of the cardiac tissue and the SL -inducedchange in calcium sensitivity (LDA, see Sec. 2.1). We remark that we are here ableto reproduce the LDA without any phenomenological SL -dependent tuning of theparameter, as done, e.g. in [62, 100–102]. The LDA emerges from the SE-ODE modelin a spontaneous way, as a consequence of the spatial-dependent interaction betweenthe RUs (see [79] for a discussion on this topic).Conversely, with the MF-ODE model it is not possible to reproduce the LDA bysimply acting on the model parameters. Indeed, the only effect of SL in the modelis to multiplicatively tune the generated force by the factor χ so ( SL ( t )). Therefore,no SL induced effect on the calcium sensitivity can be achieved. The mechanismreproducing LDA in the SE-ODE model is indeed intrinsically linked to the explicitspatial representation of the myofilaments [79]. Therefore, in the mean-field modelMF-ODE, without an explicit spatial representation, we phenomenologically tune thecalcium sensitivity k d in function of SL , by setting k d ( t ) = k d + α k d ( SL ( t ) − SL k d ) , (31)where SL k d = 2 . µ m. In Fig. 6, we show that, with this modification, we can fit theexperimental data also with the MF-ODE model. To complete the calibration of the SE-ODE and MF-ODE models, we only need toset the parameters k basic and k off , ruling the rapidity at which the transitions N − (cid:42)(cid:41) − P and
U − (cid:42)(cid:41) − B take place, respectively. Despite the fact that, at this stage, we needto calibrate just two parameters, this reveals some difficulties, mainly related to thefollowing two aspects. First, the interplay between the two parameters is tight andtheir roles cannot be easily decoupled [62, 69, 70]. This results is a poor identifiabilityof the parameters: different combinations of parameters give similar results in termsof force transients. This issue has been reported also by [98], while calibrating themodels of [57] and [62]. Additionally, the force transients predicted by the model are23ery sensitive to the calcium transient at input (this is a typical feature of activationmodels, see e.g. [98]). Therefore, since the experimentally measured calcium transientsare much affected by noise (see e.g. [4, 44, 98], a calibration based on the best fit ofthe model response to experimentally measure calcium transients should be performedwith care.Based on the former remarks, we calibrate the parameters k basic and k off by thefollowing procedure. We consider the isometric twitches of intact rat cardiac musclefibers reported in [43] for different values of SL , and with [Ca ] o = 1 . ] i ( t ) = c t < t c ,c + c max − c β (cid:20) e − t − tc τc − e − t − tc τc (cid:21) t ≥ t c , (32)where β = (cid:18) τ c τ c (cid:19) − (cid:16) τc τc − (cid:17) − − (cid:18) τ c τ c (cid:19) − (cid:16) − τc τc (cid:17) − , with the constants values c max = 1 . µ M, t c = 0 .
05 s, τ c = 0 .
02 s, τ c = 0 .
19 s.Then, we sample the candidate parameters space ( k basic , k off ) ∈ [0 ,
80 s − ] × [0 ,
300 s − ]by a MC strategy, for each parameters value we compute the tension transients corre-sponding to SL = 1 .
90, 2 .
05 and 2 . µ m and we compare them with the experimen-tal measurements from [43]. We consider the following discrepancy metrics, where T i, moda ( t ) denotes the tension predicted by the model for the i -th value of SL and T i, expa ( t ) denotes the experimentally measured one: E L := (cid:118)(cid:117)(cid:117)(cid:116) (cid:88) i =1 (cid:90) T (cid:12)(cid:12)(cid:12) T i, moda ( t ) − T i, expa ( t ) (cid:12)(cid:12)(cid:12) dt,E peak := (cid:118)(cid:117)(cid:117)(cid:116) (cid:88) i =1 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) sup t ∈ [0 ,T ] T i, moda ( t ) − sup t ∈ [0 ,T ] T i, expa ( t ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . The first metric coincides with the L error over time, while the second one evaluatesthe error of the predicted force peak. We notice indeed that the parameters k basic and k off also determine the force peak attained during a sarcomere twitch: the most rapidthe tissue activation is, the more the force-calcium curve stays close to the steady-state curve and thus it reaches higher force values before the relaxation begins. Ascriterion to select the best parameters values, we consider as overall metric a weightedcombination between the two, given by E tot = ( E L + w E ) / , where we set w peak = 5.The obtained values of the discrepancy metric E tot in the parameters space for boththe SE-ODE and the MF-ODE models are reported in Fig. 7. We notice that the levelcurves do not clearly identify an optimal pair ( k basic , k off ), rather these exhibit a wideregion in the parameters space producing very similar results. Given the uncertaintyin the measurements of both force and, mostly, calcium, it makes no sense to selectthe best parameters by blindly selecting the pair that realizes the smaller discrepancyfrom experimental results. Therefore, we supplement the results of Fig. 7 with directmeasurements of calcium binding rates to Tn, showing that k on = k off /k d lies between24 F-ODESE-ODE
Figure 7: Discrepancy metric E tot in the parameters space for intact, room-temperature rat cells, obtained with the SE-ODE model (left) and with the MF-ODEmodel (right).50 and 200 µ M − s − [62]. On this basis, we restrict the region of candidate valuesand we select the parameters reported in Tab. 7.The predicted isometric twitches obtained with the selected values of the param-eters are compared with the experimental data in Fig. 8. We notice here that theMF-ODE model accomplishes a better fit of experimental data than the SE-ODEmodel. This is an effect of the phenomenological tuning of k d of Eq. (31), that allowsfor a significant increase of calcium sensitivity and, consequently, of twitch duration,for higher values of SL . Nonetheless, also the SE-ODE model predicts, even if toa lower extent, the experimentally observed prolongation of twitches at higher SL ,without any phenomenological tuning of the calcium sensitivity.We notice that there is room for improvement in the calibration of the kineticparameters k basic and k off , which could be better constrained in presence of moreabundant and more reliable experimental data and when a deeper understanding onthe determinants of the kinetics of activation and relaxation will be available. Never-theless, the calibration of k basic and k off for the rat model does not affect the qualityof the human model, since those two parameters are the only ones to be completelyre-calibrated for the human model.The full list of the parameters, including those associated with the XB cycling,calibrated by the procedure presented in Sec. 4.1, is provided in Tab. 7. Due to the lack of a sufficiently large set of measurements from human cells at bodytemperature [56] to fit all the model parameters, to calibrate the SE-ODE and MF-ODE models for body-temperature human cardiomyocytes we start from the cor-responding sets of parameters for room-temperature rat cells and we adapt thoseparameters that are reasonably affected by the change in species and temperature.Specifically, different species mainly differ in their calcium-sensitivity (i.e. k d ) and inthe kinetics (different species feature highly different heart rates), while temperaturemainly affects the kinetics [7, 34, 45].By exploiting the decoupling of the parameters of the RUs model ruling the steady-25 t [s] T a [ k P a ] MF-ODE model [Ca ] i [ M] T a [ k P a ] t [s] T a [ k P a ] SE-ODE model [Ca ] i [ M] T a [ k P a ] SL = 1.90 m (model)SL = 1.90 m (exp) SL = 2.05 m (model)SL = 2.05 m (exp) SL = 2.20 m (model)SL = 2.20 m (exp)
Figure 8: Force transients (top line) and phase-loops (bottom line) in isometrictwitches, for different values of SL , predicted by the SE-ODE model (left column)and MF-ODE model (right column), in comparison with the experimental measure-ments from intact rat cardiac cells taken from [43].Species Temperature Preparation SL ( µ m) EC ( µ M) ReferenceRat Room Skinned 2.15 3.51 [22]Rat Room Intact 2.15 0.68 [95]Human Body Skinned 2.00 1.72 [56]Human Body Skinned 2.20 1.56 [56]Table 5: List of the experimental values used to calibrate the calcium sensitivity forthe human models at body temperature.state relationships from those ruling the kinetics (see Sec. 4.2), we first focus on thesteady-state, and we adjust k d to reflect the higher calcium sensitivity of human cellscompared to rat [56, 62, 98]. For this purpose, we employ the data reported in [56],which, however, refer to skinned cells. In order to estimate the effect of skinning on k d , we compare the calcium sensitivity measured for room-temperature rat cardiaccells in skinned [22] and intact preparations [95] at SL = 2 . µ m and we assume thatthe same relationship holds for skinned versus intact, body-temperature human cells.Finally, we rescale the values of k d to reflect the estimated change in calcium sensitivitybetween intact, body-temperature human cells and intact, room-temperature rat cells,obtaining the values reported in Tab. 7. The experimental data used in such procedureare listed in Tab. 5.Since the RUs kinetics may depend on both the species and the temperature, were-calibrate the parameters k off and k basic on the base of the output metrics reportedin Tab. 6 (taken from [98]), referred to body-temperature human cells. The outputmetric comprises the peak tension T peaka , the time-to-peak T T P (defined as the time26arameter Lower bound Upper bound Units T peaka
40 45 kPa
T T P
147 172 s RT
109 125 s RT
291 377 sTable 6: Output metrics used to calibrate the kinetics of the body-temperature humanmodels. Data are taken from [98]. The range of values of T peaka has been reduced withthe aim of more strongly forcing the models to match the desired output. MF-ODESE-ODE
Figure 9: Overall distance, in the parameters space, of the output metrics from thetarget values ranges indicated in Tab. 6 for intact, body-temperature human cells,obtained with the SE-ODE model (left) and with the MF-ODE model (right).separating the beginning of the stimulus and the tension peak) and the relaxationtimes RT and RT (defined as the time needed to accomplish 50% and 95% ofrelaxation, respectively). Since, at the best of our knowledge, detailed calcium tran-sients measurements for intact human cells at body temperature are not available, weemploy the synthetic calcium transient predicted by the ten Tusscher-Panfilov 2006ionic model [93], that we denote as the TTP06 model. In Fig. 9 we show the land-scape, in the parameters space ( k basic , k off ) of the overall distance of the output metricsfrom the ranges reported in Tab. 6. Similarly to the rat models, large parameters re-gions feature similar values. Therefore, we proceed as for the calibration of rat cells,by reducing the parameters space based on direct experimental measurements (seeSec. 4.2). In this case, we manage to select a unique couple of values suitable for boththe SE-ODE and the MF-ODE model (see Tab. 7).Finally, for the calibration of the parameters ruling the XBs cycling, we use thesame values of Tab. 4. Therefore, since the calibration depends on the parameterspreviously set for the RUs activation, the resulting values of the parameters are slightlydifferent. We provide in Tab. 7 the full list of parameters for both species (room-temperature rat and body-temperature human) and for both models (SE-ODE andMF-ODE). 27E-ODE MF-ODEParameter Units Rat,room temp. Human,body temp. Rat,room temp. Human,body temp. RU steady-state µ - 10 10 10 10 γ - 20 20 12 12 Q - 3 3 2 2 k d µ M 1 .
347 0 .
615 0 .
775 0 . α k d µ M µ m − − . − . RU kinetics k off s −
110 121 135 121 k basic s − . . . . XB cycling µ f P s − .
421 58 .
192 33 .
161 33 . µ f P s − .
390 1 .
384 0 .
789 0 . r s − .
958 134 .
075 134 .
182 134 . α - 25 .
113 25 .
129 25 .
149 25 . Micro-macro upscaling a XB MPa 22 .
941 22 .
931 22 .
919 22 . We show the results of numerical simulations obtained with the SE-ODE and MF-ODE models, performed under different experimental settings. The numerical schemesemployed to approximate the solutions of the models here proposed are described inApp. B. The codes implementing the models proposed in this paper and used toproduce the results presented in this section are freely available in the following onlinerepository: https://github.com/FrancescoRegazzoni/cardiac-activation
With this implementation, the computational cost of the SE-ODE model is of nearly12 s of simulation for 1 s of physical time on a single-core standard laptop. The MF-ODE model, instead, requires nearly 1 s of computational time to simulate 1 s ofphysical time on the same computer platform.
First, we consider steady-state solutions. To numerically obtain the steady-statecurves, we fix a level of [Ca ] i and SL and we let the model reach the equilibriumsolution. We show in Figs. 10 and 11 the force-calcium curves predicted by the SE-ODE andMF-ODE models for rat and human cells, respectively. In all the cases, the curvesreproduce the main experimentally observed features reported in Sec. 2.1. Figure 1228 -1 [Ca ] i [ M] T a [ k P a ] SE-ODE model -1 [Ca ] i [ M] T a [ k P a ] MF-ODE model
SL [ m]
Figure 10: Steady-state force-calcium relationship at different SL obtained with theSE-ODE (left) and the MF-ODE (right) models for intact, room-temperature ratcardiomyocytes. -1 [Ca ] i [ M] T a [ k P a ] SE-ODE model -1 [Ca ] i [ M] T a [ k P a ] MF-ODE model
SL [ m]
Figure 11: Steady-state force-calcium relationship at different SL obtained with theSE-ODE (left) and the MF-ODE (right) models for intact, body-temperature humancardiomyocytes.shows the dependence of the Hill coefficient n H and of the half-activating calciumconcentration EC on the sarcomere length SL . We notice that, while the MF-ODEmodel produces an Hill coefficient that is independent of SL (the reason is that therole of SL on activation is just that of shifting and rescaling the curves, thus leaving n H unaffected), the SE-ODE model predicts a small increase of n H with SL . Boththe results are equally acceptable since there is no common agreement on whether theHill coefficient depend on SL or not [22, 46, 50, 95, 96]. Both the models correctlypredict for both species an increase of EC as SL decreases. The relationship isapproximately linear in the typical working range of SL (as experimentally observed,e.g., in [22]), while, for small values of SL , the SE-ODE model produces a fasterdecrease of sensitivity. Figures 13 and 14 show the ascending limb of the steady-state force-length relation-ships. For both the SE-ODE and MF-ODE models we observe a change of slope forsaturating calcium concentration around 1 . µ m, coherently with the experimentalobservations [49, 51, 95, 96]. Moreover, both models predict the change of convexityof the force-length curves at different calcium levels [50, 95, 96].29 .6 1.8 2 2.2 SL [ m] n H [ - ] SL [ m] E C [ M ] SE-ODE model (room-temperature rat)MF-ODE model (room-temperature rat) SE-ODE model (body-temperature human)MF-ODE model (body-temperature human)
Figure 12: Dependence of the Hill coefficient n H and of the half-activating calciumconcentration EC on the sarcomere length SL , for the SE-ODE (blue lines) and MF-ODE (red lines) model, calibrated for intact, room-temperature rat cardiomyocytes(dashed lines) and for intact, body-temperature human cardiomyocytes (solid lines). SL [ m] T a [ k P a ] SE-ODE model [Ca ] i [ M] SL [ m] T a [ k P a ] MF-ODE model [Ca ] i [ M] Figure 13: Steady-state force-length relationship at different [Ca ] i obtained withthe SE-ODE (left) and the MF-ODE (right) models for intact, room-temperature ratcardiomyocytes. SL [ m] T a [ k P a ] SE-ODE model [Ca ] i [ M] SL [ m] T a [ k P a ] MF-ODE model [Ca ] i [ M] Figure 14: Steady-state force-length relationship at different [Ca ] i obtained with theSE-ODE (left) and the MF-ODE (right) models for intact, body-temperature humancardiomyocytes. 30 t [s] T a [ k P a ] SE-ODE model t [s] T a / T a m a x [ - ] SL [ m] t [s] T a [ k P a ] MF-ODE model t [s] T a / T a m a x [ - ] SL [ m]
Figure 15: Tension transients during isometric twitches at different SL obtained withthe SE-ODE (left) and the MF-ODE (right) models for intact, body-temperaturehuman cardiomyocytes. We show in Fig. 15 the tension transients obtained by simulating isometric twitchesgiving as input to the human models the calcium transient of the TTP06 model (wehave already reported similar numerical results for the rat models in Fig. 8). We noticethat both models predict the tension-dependent prolongation of the relaxation time[4, 43, 44], as it can be seen from the normalized traces reported in the bottom linesof the figures. Moreover, we report in Fig. 16 the dependence of the time indicators
T T P , RT and RT on SL . Coherently with the experimental measurements, the T T P is nearly independent of SL , while RT shows an increasing trend with SL . SL [ m] [ s ] SE-ODE model
TTP RT RT SL [ m] [ s ] MF-ODE model
TTP RT RT Figure 16: Metrics of activation and relaxation kinetics as function of SL duringisometric twitches obtained with the SE-ODE (left) and the MF-ODE (right) modelsfor intact, body-temperature human cardiomyocytes.31 T a /T aiso [-] -20246810 v [ h s / s ] SE-ODE model
Ca = 1.0 M, SL = 1.9 m (exp)Ca = 1.0 M, SL = 2.2 m (exp)Ca = 2.5 M, SL = 2.2 m (exp) Ca = 1.0 M, SL = 1.9 m (model)Ca = 1.0 M, SL = 2.2 m (model)Ca = 2.5 M, SL = 2.2 m (model) T a /T aiso [-] -20246810 v [ h s / s ] MF-ODE model
Ca = 1.0 M, SL = 1.9 m (exp)Ca = 1.0 M, SL = 2.2 m (exp)Ca = 2.5 M, SL = 2.2 m (exp) Ca = 1.0 M, SL = 1.9 m (model)Ca = 1.0 M, SL = 2.2 m (model)Ca = 2.5 M, SL = 2.2 m (model)
Figure 17: Normalized force-velocity relationships for different combinations of [Ca ] i and SL obtained with the SE-ODE (left) and the MF-ODE (right) models for intact,room-temperature rat cardiomyocytes in comparison with experimental measurementsfrom [11]. Figures 17 and 18 show the force-velocity relationship predicted with the models cal-ibrated from rat and human data, respectively. In both the cases, and for both theSE-ODE and the MF-ODE model, the numerical results correctly reproduce the ex-perimentally observed convex profile, with the force reaching zero in correspondenceof a finite value of velocity, the so-called maximum shortening velocity (see Sec. 2.2).Moreover, the value of v maxhs is not significantly affected by the level of activation (thatis to say, by the values of [Ca ] i and SL ), as well as the curvature of the curve. Thisis also coherent with the experimental observations [11].In Fig. 17 we also compare the force-velocity curves obtained by the models withthe experimental data used to calibrate the rat models (from [11]), thus validatingthat the automatic calibration procedure presented in Sec. 4.1 has been successful. We consider the response to fast steps predicted by the SE-ODE and the MF-ODEmodels. With this aim, we set a fixed value for the calcium concentration and sar-comere length (we set [Ca ] i = 1 . µ M and SL = 2 . µ m, but the results are notsignificantly affected by this choice) and we let the system reach the steady-state.Then, we apply a length step, by applying a constant shortening velocity in a smalltime interval ∆ t = 200 µ s, and we plot the tension at the end of the step as a functionof the step length ∆ L .We show in Figs. 19 and 20 the results obtained for the rat and the human models,respectively, with a comparison for the first case with experimental data. The goodmatch between the simulation results and the experimental measurements provide afurther validation of the calibration procedure.32 T a /T aiso [-] -20246810 v [ h s / s ] SE-ODE model
Ca = 0.4 M, SL = 2.0 mCa = 0.5 M, SL = 2.0 mCa = 0.6 M, SL = 2.0 m Ca = 0.4 M, SL = 2.2 mCa = 0.5 M, SL = 2.2 mCa = 0.6 M, SL = 2.2 m T a /T aiso [-] -20246810 v [ h s / s ] MF-ODE model
Ca = 0.4 M, SL = 2.0 mCa = 0.5 M, SL = 2.0 mCa = 0.6 M, SL = 2.0 m Ca = 0.4 M, SL = 2.2 mCa = 0.5 M, SL = 2.2 mCa = 0.6 M, SL = 2.2 m
Figure 18: Normalized force-velocity relationships for different combinations of [Ca ] i and SL obtained with the SE-ODE (left) and the MF-ODE (right) models for intact,body-temperature human cardiomyocytes. -10 -8 -6 -4 -2 0 L [nm/hs] T / T i s o [ - ] SE-ODE model experimental model -10 -8 -6 -4 -2 0
L [nm/hs] T / T i s o [ - ] MF-ODE model experimental model
Figure 19: Normalized force after the application of a fast length step ∆ L for in-tact, room-temperature rat cardiomyocytes. Solid line: model result; circles: T - L experimental data from [11]. -10 -8 -6 -4 -2 0 L [nm/hs] T / T i s o [ - ] SE-ODE model model -10 -8 -6 -4 -2 0
L [nm/hs] T / T i s o [ - ] MF-ODE model model
Figure 20: Normalized force after the application of a fast length step ∆ L for intact,body-temperature human cardiomyocytes.33 .5 Multiscale cardiac electromechanics Finally, we consider a multiscale cardiac EM model of the left ventricle (LV), thatwe describe in the following. Further details on the EM model and its numericaldiscretization can be found in [79]. For the sake of brevity, we only show the resultsobtained by considering the MF-ODE model for human body-temperature cardiomy-ocytes.We employ a computational domain Ω derived from the Zygote heart model, rep-resenting the 50 th percentile of a 21 years old healthy caucasian male [106], in whichwe generate the fibers and sheets distribution by the rule-based algorithm proposedin [6]. We model the electrophysiological activity of the heart tissue by means of themonodomain equation [17, 18], coupled with the TTP06 model for the ionic activity[93]. We introduce a deformation map ϕ : Ω × [0 , T ] → R and we define the displace-ment vector as d ( X , t ) = ϕ ( X , t ) − X . The mechanical behavior of the myocardiumis described by means of the balance of momentum equation for the displacement d [2, 65], where we model the passive properties of the cardiac tissue through the quasi-incompressible exponential constitutive law of [99]. To account for the presence of thepericardium, we employ the generalized Robin boundary conditions of [30, 68] on theepicarial boundary. On the LV base we adopt the boundary conditions derived in [78],accounting for the effect of the neglected part of the domain on the artificial boundaryof the LV base.Within a multiscale setting, we describe the force generation mechanisms at themicroscale by means of the MF-ODE model, proposed in this paper. The intracel-lular calcium concentration ([Ca ] i ) is provided in each point of the computationaldomain by the TTP06 model. The local sarcomere length, in turn, is assumed to beproportional to the tissue stretch in the direction of muscel fibers f , that is we set SL = SL (cid:112) I ,f , where SL denotes the sarcomere length at rest, I ,f = Ff · Ff and F = I + ∇ d denotes the deformation gradient. The last input of the MF-ODE model,the normalized shortening velocity, is obtained by definition as v = − ∂∂t SL/SL .The output of the MF-ODE model, namely the active tension magnitude field T a ,provides the link from the microscopic to the macroscopic level, where the effect ofthe microscopically generated active force is given by the following active stress tensor[31, 79]: P act = T a Ff ⊗ f (cid:107) Ff (cid:107) . Finally, the 3D LV model is coupled to a 0D model of blood circulation, consisting ofa two-elements Windkessel model [103].For the spatial discretization of the equations involved in the EM model, we em-ploy piece-wise linear Finite Elements of a tetrahedral computational mesh. For thediscretization of time derivatives, we use first-order finite difference schemes [73] withthe implicit-explicit (IMEX) scheme described in [79]. The coupling of the differentcore models is obtained by means of the segregated strategy presented in [21].The results of the numerical simulation are provided in Fig. 21, where we show thedeformation of the LV at different time steps of an heartbeat, and in Fig. 22, showingthe time evolution of the variables involved in the force generation phenomenon andof the LV pressure and volume predicted by the multiscale EM model. In Fig. 23 weshow the pressure-volume loops obtained for different preloads, i.e. by varying theend-diastolic pressure p ED . Our model correctly predicts the increased stroke volumefollowing a raise in the preload. This phenomenon, known as Frank-Starling effect,represents a self-regulatory mechanism of fundamental importance, as it allows the34 = 0 s t = 0 .
10 s t = 0 .
15 s t = 0 .
40 s t = 0 .
60 s
Figure 21: LV multiscale EM: deformed geometry and magnitude of displacement d atdifferent times. Top row: full geometry. Middle row: half domain (top view). Bottomrow: half domain (frontal view).cardiac output to be synchronized with the venous return [47]. The microscopic driverof the Frank-Starling effect is the length-dependent mechanisms of force generation.We remark that, in our EM model, the Frank-Starling is not artificially imposed,but it spontaneously emerges from the properties of the microscopic model of forcegeneration. In this paper we have derived four different models for the generation of active forcein the cardiac muscle tissue. Such models are based on a biophysically detailed de-scription of the microscopic mechanisms of regulation and activation of the contractileproteins. Still, their numerical realization features a contained computational cost. In-deed, their numerical resolution does not require the computationally expensive MCmethod.The main difficulties to be addressed in the derivation of full-sarcomere modelsconcern (i) the spatial correlation of the states of the RUs due to the nearest-neighborinteractions of Tm, which hinders a straightforward decoupling of the adjacent RUs,and (ii) the mutual filament sliding, that changes which RU regulates which XB from35 . . . . . [ C a + ] i [ µ M ] minavgmax 0 0 . . . . . S L [ µ m ] minavgmax 0 0 . . . . − − v [ s − ] minavgmax0 0 . . . . T a [ k P a ] minavgmax 0 0 . . . . P r e ss u r e [ mm H g ] . . . . V o l u m e [ m L ] Figure 22: LV multiscale EM. Time evolution of [Ca ] i , SL , v and T a (minimum,maximum and average over the computational domain) and of left ventricle pressureand volume.
40 60 80 100 120 140020406080100 Volume [mL] P r e ss u r e [ mm H g ] p ED = 7 . p ED = 10 mmHg p ED = 12 . p ED = 15 mmHg Figure 23: Pressure-volume loops obtained with different preloads (reported in leg-end). 36ime to time. Similarly to what we did in [76], we address the first problem byintroducing a conditional independence assumption for far RUs, given the state ofintermediate RUs (this is coherent with the local nature of nearest-neighboring in-teractions). This assumption dramatically reduces the number of equations neededto describe the evolution of the stochastic processes: while in the original model of[84] the number of variables of the CTMC increases exponentially with the number ofRUs, in our model the number of variables grows linearly with the number of RUs.Since no feedback from the XBs to the RUs is foreseen, the dynamics of the lattercan be considered independently of that of the former. Moreover, we depart from thetraditional MHs-centered representation of XBs, in favour of a BSs-centered point ofview. Thanks to this change of perspective, we derived a set of equations describingthe XB dynamics without the need to track the mutual position of the RUs and theMHs.Under the hypothesis that the total attachment-detachment rate is independentof the myosin arm stretch (as done in [8, 16]), the PDE system describing the XBscan be replaced by a system of ODEs. We remark that this is not a simplificatoryassumption, like the conditional independence assumptions mentioned before, butrather a feature of the specific modeling choice for the transition rates describing theattachment-detachment process.We have also presented a class of models (MF-PDE and MF-ODE), such thatthe myofilaments overlap is not explicitly described, but is rather replaced by a mean-field description of a single representative RUs triplet. We remark that such mean-fieldmodels differ from the mean-field models of [58, 75, 82, 85, 86]. The latter, indeed,consider a single RU, instead of a triplet. In this manner, the short-range spatialcorrelation, responsible of cooperativity, cannot be captured. Conversely, the mean-field triplet framework here proposed, thanks to the local nature of cooperativity,allows to capture the effect of nearest-neighbor interactions, as testified by the re-markably good agreement between model predictions as experimental measurements,in particular in the reproduction of the highly cooperative steady-state force-calciumcurves (see Sec. 5.1). We have then calibrated the SE-ODE and the MF-ODE modelsfor room-temperature rat intact cardiomyocytes and, later, body-temperature humanintact cardiomyocytes.The SE-ODE model predicts the so-called LDA (the increment of calcium sensi-tivity when the sarcomere stretch lies within the physiological range), thanks to theexplicit spatial representation of the spatially-dependent nearest-neighborhood inter-actions among RUs. We remark that there is not a common agreement in the scientificcommunity about the microscopical phenomena producing the above mentioned effect[1, 25, 28, 63, 67, 81, 88, 94, 97]. Therefore, we believe that our results are remarkablein this sense, because the phenomena of LDA spontaneously emerges without any phe-nomenological tuning the parameters in dependence of SL . This is due to the effectthat the RUs located at the end-points of the single-overlap zone, characterized by alow probability of being in the P state, have on the whole filament (further details onthis topic can be found in [79]). Therefore, if the hypothesis that the RUs located atthe end-point of the single-overlap zone behave as if the outer neighboring units arein state N is accepted, then the SE-ODE model provides a possible explanation forLDA. Conversely, if this hypothesis in not accepted, then the effect of the edges canbe neglected, so that the mean-field approximation underlying the MF-ODE modelis well motivated. In conclusion, according to the hypothesis made on the behaviorof the RUs near the end-points of the single-overlap zone, the SE-ODE and the MF-ODE models represent two alternative descriptions of the sarcomere dynamics based37n a microscopically detailed representation of the regulatory and contractile proteins,where phenomenological modeling choices are only introduced for phenomena whoseunderlying mechanisms are not fully understood [1, 63, 88, 97].The results of the numerical simulations showed that our models can reproduce themain features of the experimental characterizations of muscle contraction associatedwith the time scales of interest of cardiac EM (that is to say, the time scales larger thanthat of the power-stroke), including the steady-state relationship between calcium andforce and between sarcomere length and force, the kinetics of activation of relaxation,the force-dependent twitches prolongation, and the force-velocity relationship. Finally,we have presented the results of multiscale EM numerical simulations in a human LV,obtained by modeling the microscopic generation of active force through the MF-ODEmodel, able of reproducing realistic pressure-volume loops. Moreover, our multiscaleEM model is capable of reproducing – at the macroscale – the Frank-Starling self-regulation mechanism, in virtue of the length-dependent effects characterizing – atthe microscale – the force-generation model. Acknowledgements
This project has received funding from the European Research Council (ERC) un-der the European Union’s Horizon 2020 research and innovation programme (grantagreement No 740132, iHEART - An Integrated Heart Model for the simulation of thecardiac function, P.I. Prof. A. Quarteroni).
A Models derivation
We report the detailed derivation of the models proposed in Sec. 3. The derivation isbased on the assumptions discussed in Sec. 3.3. These assumptions allow to neglectsecond-order interactions among the stochastic processes, so that the variables canbe partially decoupled, thus leading to drastic reductions in the size of model. Suchstrategy is illustrated in the following proposition.
Proposition 1.
Let (Ω , A , P ) be a probability space. Let A, B, C ⊂ Ω and let D be afinite partition of Ω , such that:(H1) A ⊥⊥ B | C, D ∀ D ∈ D ;(H2) B ⊥⊥ D | C ∀ D ∈ D .Then, we have: P [ A | B, C ] = (cid:80) D ∈D P [ A | C, D ] P [ C, D ] P [ C ] = P [ A | C ] . roof. We have: P [ A | B, C ] = P [ A, B, C ] P [ B, C ] = (cid:80) D ∈D P [ A, B, C, D ] P [ B, C ] = (cid:80) D ∈D P [ A | B, C, D ] P [ B, C, D ] P [ B, C ]= (cid:80) D ∈D P [ A | B, C, D ] P [ B | C, D ] P [ C, D ] P [ B, C ] . From (H1), it follows: P [ A | B, C, D ] = P [ A | C, D ] . Moreover, in virtue of (H2), we have: P [ B | C, D ] = P [ B | C ] = P [ B, C ] / P [ C ] . By substituting into the above equation, the thesis follows.In the following, we will use several times the result of Prop. 1, where (H1) is amodeling choice on the dynamics of the system and (H2) is a simplifying assumption.Specifically, A is the target event, whose probability is the aim of the computation. Inmany situations, we know the joint probability of C and B , whereas the probability of A can be obtained by the joint probability of C and a different event D . Proposition 1allows to pass from B to D , by assuming that the knowledge of B does not provideany further information when C and D are known. A.1 Derivation of Eq. (8)
Le us consider the time increment ∆ t and let us compute the probability π αβδ,ϑηλi ( t + ∆ t ).In virtue of the Bayes formula [64], we have: π αβδ,ϑηλi ( t + ∆ t ) ∆ t → ∼ P (cid:2) T t +∆ ti − = α | ( T i − , T i , T i +1 ) t = ( α, β, δ ) , ( C i − , C i , C i +1 ) t = ( ϑ, η, λ ) (cid:3) π αβδ,ϑηλi ( t )+ P (cid:2) T t +∆ ti = β | ( T i − , T i , T i +1 ) t = ( α, β, δ ) , ( C i − , C i , C i +1 ) t = ( ϑ, η, λ ) (cid:3) π αβδ,ϑηλi ( t )+ P (cid:2) T t +∆ ti +1 = δ | ( T i − , T i , T i +1 ) t = ( α, β, δ ) , ( C i − , C i , C i +1 ) t = ( ϑ, η, λ ) (cid:3) π αβδ,ϑηλi ( t )+ P (cid:2) C t +∆ ti − = ϑ | ( T i − , T i , T i +1 ) t = ( α, β, δ ) , ( C i − , C i , C i +1 ) t = ( ϑ, η, λ ) (cid:3) π αβδ,ϑηλi ( t )+ P (cid:2) C t +∆ ti = η | ( T i − , T i , T i +1 ) t = ( α, β, δ ) , ( C i − , C i , C i +1 ) t = ( ϑ, η, λ ) (cid:3) π αβδ,ϑηλi ( t )+ P (cid:2) C t +∆ ti +1 = λ | ( T i − , T i , T i +1 ) t = ( α, β, δ ) , ( C i − , C i , C i +1 ) t = ( ϑ, η, λ ) (cid:3) π αβδ,ϑηλi ( t )+ P (cid:2) ( T i − , T i , T i +1 ) t +∆ t = ( α, β, δ ) , ( C i − , C i , C i +1 ) t +∆ t = ( ϑ, η, λ ) | ( T i − , T i , T i +1 ) t = ( α, β, δ ) , ( C i − , C i , C i +1 ) t = ( ϑ, η, λ ) (cid:3) π αβδ,ϑηλi ( t ) , where, by definition, we have: P (cid:2) T t +∆ ti = β | ( T i − , T i , T i +1 ) t = ( α, β, δ ) , ( C i − , C i , C i +1 ) t = ( ϑ, η, λ ) (cid:3) ∆ t → ∼ k ββ | α · δ,ηT,i ∆ t, and P (cid:2) C t +∆ ti − = ϑ | ( T i − , T i , T i +1 ) t = ( α, β, δ ) , ( C i − , C i , C i +1 ) t = ( ϑ, η, λ ) (cid:3) ∆ t → ∼ k ϑϑ | αC,i − ∆ t, P (cid:2) C t +∆ ti = η | ( T i − , T i , T i +1 ) t = ( α, β, δ ) , ( C i − , C i , C i +1 ) t = ( ϑ, η, λ ) (cid:3) ∆ t → ∼ k ηη | βC,i ∆ t, P (cid:2) C t +∆ ti +1 = λ | ( T i − , T i , T i +1 ) t = ( α, β, δ ) , ( C i − , C i , C i +1 ) t = ( ϑ, η, λ ) (cid:3) ∆ t → ∼ k λλ | δC,i +1 ∆ t.
39y adopting assumption (I2) and applying Prop. 1 for A = ( T t +∆ ti − = α ), B = ( T ti +1 = δ, C ti +1 = λ ), C = (( T i − , T i ) t = ( α, β ) , ( C i − , C i ) t = ( ϑ, η )) and D = { ( T ti − = ξ, C ti − = ζ ) } ξ,ζ , we have: P (cid:2) T t +∆ ti − = α | ( T i − , T i , T i +1 ) t = ( α, β, δ ) , ( C i − , C i , C i +1 ) t = ( ϑ, η, λ ) (cid:3) = (cid:88) ξ,ζ π ξαβ,ζϑηi − ( t ) − (cid:88) ξ,ζ P (cid:2) T t +∆ ti − = α | ( T i − , T i − , T i ) t = ( ξ, α, β ) , ( C i − , C i − , C i ) t = ( ζ, ϑ, η ) (cid:3) π ξαβ,ζϑηi − ( t )= (cid:80) ξ,ζ k αα | ξ · β,ϑT,i π ξαβ,ζϑηi − ( t ) (cid:80) ξ,ζ π ξαβ,ζϑηi − ( t ) ∆ t + o (∆ t ) , and similarly for the term related to T t +∆ ti +1 . In conclusion, by taking the limit ∆ t → A.2 Derivation of Eq. (20)
Similarly to Sec. A.1, we consider a finite time increment ∆ t and we write: π αβδ,η ( t + ∆ t ) ∆ t → ∼ P (cid:2) T t +∆ ti − = α | ( T i − , T i , T i +1 ) t = ( α, β, δ ) , C ti = η (cid:3) π αβδ,η ( t )+ P (cid:2) T t +∆ ti = β | ( T i − , T i , T i +1 ) t = ( α, β, δ ) , C ti = η (cid:3) π αβδ,η ( t )+ P (cid:2) T t +∆ ti +1 = δ | ( T i − , T i , T i +1 ) t = ( α, β, δ ) , C ti = η (cid:3) π αβδ,η ( t )+ P (cid:2) C t +∆ ti = η | ( T i − , T i , T i +1 ) t = ( α, β, δ ) , C ti = η (cid:3) π αβδ,η ( t )+ P (cid:2) ( T i − , T i , T i +1 ) t +∆ t = ( α, β, δ ) , C t +∆ ti = η | ( T i − , T i , T i +1 ) t = ( α, β, δ ) , C ti = η (cid:3) π αβδ,η ( t ) , where, by definition of the transition rates, it holds: P (cid:2) T t +∆ ti = β | ( T i − , T i , T i +1 ) t = ( α, β, δ ) , C ti = η (cid:3) ∆ t → ∼ k ββ | α · δ,ηT ∆ t, and P (cid:2) C t +∆ ti = η | ( T i − , T i , T i +1 ) t = ( α, β, δ ) , C ti = η (cid:3) ∆ t → ∼ k ηη | βC ∆ t. By adopting assumption (I1), Prop. 1 for A = ( T t +∆ ti − = η ), B = ( T ti +1 = δ, C ti = η ), C = (( T i − , T i ) t = ( α, β )) and D = { ( T ti − = ξ, C ti − = ζ ) } ξ,ζ leads to: P (cid:2) T t +∆ ti − = α | ( T i − , T i , T i +1 ) t = ( α, β, δ ) , C ti = η (cid:3) = (cid:80) ξ,ζ P (cid:2) T t +∆ ti − = α | ( T i − , T i − , T i ) t = ( ξ, α, β ) , C ti − = ζ ) (cid:3) π ξαβ,ζi − ( t ) (cid:80) ξ,ζ π ξαβ,ζi − ( t ) ∆ t → ∼ (cid:101) k αα |◦ · β, ◦ T ∆ t. In conclusion, by letting ∆ t →
0, we get Eq. (20).40 .3 Derivation of Eqs. (11) and (22)
Considering for instance the variable n i, P ( x, t ) (similar calculations can be carried outfor n i, N ( x, t )), we have: f (cid:2) ( Z i , T i ) t +∆ t = ( x − v hs ( t )∆ t, P ) (cid:3) ∆ t → ∼ f (cid:2) ( Z i , T i ) t +∆ t = ( x − v hs ( t )∆ t, P ) | ( Z i , T i ) t = ( ∅ , P ) (cid:3) P (cid:2) ( Z i , T i ) t = ( ∅ , P ) (cid:3) + f (cid:2) ( Z i , T i ) t +∆ t = ( x − v hs ( t )∆ t, P ) | ( Z i , T i ) t = ( x, N ) (cid:3) f (cid:2) ( Z i , T i ) t = ( x, N ) (cid:3) + f (cid:2) ( Z i , T i ) t +∆ t = ( x − v hs ( t )∆ t, P ) | ( Z i , T i ) t = ( x, P ) (cid:3) f (cid:2) ( Z i , T i ) t = ( x, P ) (cid:3) . Thanks to Prop. 1, by taking A = (( Z i , T i ) t +∆ t = ( x − v hs ( t )∆ t, P )), B = ( Z ti = x ), C = ( T ti = N ) and D = { ( T i − , T i +1 , C i ) t = ( α, η, δ ) } α,η,δ , assumption (I3) leads to f (cid:2) ( Z i , T i ) t +∆ t = ( x − v hs ( t )∆ t, P ) | ( Z i , T i ) t = ( x, N ) (cid:3) ∆ t → ∼ (cid:101) k N P
T,i ∆ t. where we have defined: (cid:101) k N P
T,i := (cid:80) α,η,δ k N P| α · η,δT,i P [( T i − , T i , T i +1 , C i ) t = ( α, N , η, δ )] P [ T ti = N ] , (cid:101) k PN T,i := (cid:80) α,η,δ k PN | α · η,δT,i P [( T i − , T i , T i +1 , C i ) t = ( α, P , η, δ )] P [ T ti = P ] . We notice that the transition rates (cid:101) k N P
T,i and (cid:101) k PN T,i can be obtained from the variables π αβδ,ϑηλi as in Eq. (12). Moreover, we have: P (cid:2) ( Z i , T i ) t +∆ t = ( x − v hs ( t )∆ t, P ) | ( Z i , T i ) t = ( x, P ) (cid:3) ∆ t → ∼ − P (cid:2) ( Z i , T i ) t +∆ t = ( ∅ , P ) | ( Z i , T i ) t = ( x, P ) (cid:3) − P (cid:2) ( Z i , T i ) t +∆ t = ( x − v hs ( t )∆ t, N ) | ( Z i , T i ) t = ( x, P ) (cid:3) ∆ t → ∼ − ∆ t (cid:16) g i P ( x, v ( t )) − (cid:101) k PN T,i (cid:17) , where we have applied once again assumption (I3). Concerning the XB formationterm, we have:( F ) := f (cid:2) ( Z i , T i ) t +∆ t = ( x − v hs ( t )∆ t, P ) | ( Z i , T i ) t = ( ∅ , P ) (cid:3) P (cid:2) ( Z i , T i ) t = ( ∅ , P ) (cid:3) = f (cid:2) ( Z i , T i ) t +∆ t = ( x − v hs ( t )∆ t, P ) , ( Z i , T i ) t = ( ∅ , P ) (cid:3) = f (cid:2) ( Z i , T i ) t +∆ t = ( x − v hs ( t )∆ t, P ) , ( Z i , T i ) t = ( ∅ , P ) , ∃ j ∈ I M : d tij = x, M tj = 0 (cid:3) + f (cid:2) ( Z i , T i ) t +∆ t = ( x − v hs ( t )∆ t, P ) , ( Z i , T i ) t = ( ∅ , P ) , ∃ j ∈ I M : d tij = x, M tj (cid:54) = 0 (cid:3) + f (cid:2) ( Z i , T i ) t +∆ t = ( x − v hs ( t )∆ t, P ) , ( Z i , T i ) t = ( ∅ , P ) , ∃ j ∈ Z \ I M : d tij = x (cid:3) . The last two terms are at least of second order in ∆ t for ∆ t →
0, while the first termgives: f (cid:2) ( Z i , T i ) t +∆ t = ( x − v hs ( t )∆ t, P ) , ( Z i , T i ) t = ( ∅ , P ) , ∃ j ∈ I M : d tij = x, M tj = 0 (cid:3) = P (cid:2) ( Z i , T i ) t +∆ t = ( x − v hs ( t )∆ t, P ) | ( Z i , T i ) t = ( ∅ , P ) , ∃ j ∈ I M : d tij = x, M tj = 0 (cid:3) f (cid:2) ( Z i , T i ) t = ( ∅ , P ) , ∃ j ∈ Z : d tij = x, M tj = 0 (cid:3) ∆ t → ∼ f i P ( x, v hs ( t )) f (cid:2) ( Z i , T i ) t = ( ∅ , P ) , ∃ j ∈ Z : d tij = x, M tj = 0 (cid:3) ∆ t ;41he remaining two terms are null. Thus:( F ) ∼ f i P ( x, v ( t )) ∆ t f (cid:2) ( Z i , T i ) t = ( ∅ , P ) , ∃ j ∈ I M : d tij = x, M tj = 0 (cid:3) . By assumption (A1), for any i and x such that f i P ( x, v ( t )) (cid:54) = 0, the event ( M tj = 0)for j s.t. d tij = x implies the event ( Z ti = ∅ ), thus: f (cid:2) ( Z i , T i ) t = ( ∅ , P ) , ∃ j ∈ I M : d tij = x, M tj = 0 (cid:3) = f (cid:2) ( Z i , T i ) t = ( ∅ , P ) , ∃ j ∈ I M : d tij = x (cid:3) = ( f (cid:2) T ti = P , ∃ j ∈ I M : d tij = x (cid:3) − (cid:88) k f (cid:2) ( Z i , T i ) t = ( x + kD M , P ) , ∃ j ∈ I M : d tij = x (cid:3) ) , since a BS can be only attached with displacements that are multiple of D M . Moreover,we recall that the RU dynamics is independent of the interaction with XBs and thusof d (see Sec. 3.2.1) and that for i and x such that f i P ( x, v ( t )) (cid:54) = 0 the events ( ∃ j ∈I M : d tij = x ) and ( ∃ j ∈ Z : d tij = x ) coincide. Therefore, we have (on the support of f i P ): f (cid:2) T ti = P , ∃ j ∈ I M d tij = x (cid:3) = P (cid:2) T ti = P (cid:3) f (cid:2) ∃ j ∈ Z : d tij = x (cid:3) . In addition, since ( Z i = x + kD M ) implies ( ∃ j ∈ Z : d tij = x ), on the support of f i P itholds true: f (cid:2) ( Z i , T i ) t = ( x + kD M , P ) , ∃ j ∈ I M : d tij = x (cid:3) = f (cid:2) ( Z i , T i ) t = ( x + kD M , P ) (cid:3) . Since assumption (A1) implies (A2), the unique nonzero term of the sum is k = 0 andthus:( F ) ∼ = f i P ( x, v ( t ))∆ t ( P (cid:2) T ti = P (cid:3) f (cid:2) ∃ j ∈ Z : d tij = x (cid:3) − f (cid:2) ( Z i , T i ) t = ( x, P ) (cid:3) ) . Finally, we divide everything by ∆ t we let ∆ t → n i, P ( x − v hs ( t )∆ t, t + ∆ t ) − n i, P ( x, t )∆ t = n i, P ( x − v hs ( t )∆ t, t + ∆ t ) − n i, P ( x − v hs ( t )∆ t, t )∆ t + n i, P ( x − v hs ( t )∆ t, t ) − n i, P ( x, t )∆ t v hs ( t ) v hs ( t ) → ∂n i, P ∂t ( x, t ) − v hs ( t ) ∂n i, P ∂x ( x, t ) . We get in such a way Eq. (11). Moreover, the expected value of the force exerted bythe whole half filament is given by: F hf ( t ) = (cid:88) i (cid:90) + ∞−∞ F XB ( x ) f (cid:2) Z ti = x (cid:3) dx = (cid:88) i (cid:90) + ∞−∞ F XB ( x ) (cid:0) f (cid:2) Z ti = x, T ti = P (cid:3) + f (cid:2) Z ti = x, T ti = N (cid:3)(cid:1) dx = (cid:88) i (cid:90) + ∞−∞ F XB ( x ) ( n i, P ( x, t ) + n i, N ( x, t )) dx. On the other hand, Eq. (22) can be derived similarly to Eq. (11), by dropping thedependence on the RU index i . 42 .4 Derivation of Eqs. (15) and (26) By following [104], we multiply Eq. (11) by ( xSL / ) p , for p = 0 ,
1, and we integratewith respect to x over the real line. Thanks to the fact that, for x → ±∞ , thedistributions n i,α are definitively equal to zero, for α ∈ {N , P} and for i ∈ I A , theconvective terms give raise to the following terms. For p = 0, we have: (cid:90) + ∞−∞ v hs ∂n i,α ∂x dx = [ n i,α ] + ∞−∞ = 0 . On the other hand, for p = 1, we have: (cid:90) + ∞−∞ xSL / v hs ∂n i,α ∂x dx = − (cid:90) + ∞−∞ v hs SL / n i, P dx + v hs (cid:20) xSL / n i,α (cid:21) + ∞−∞ = − v µ i,α ( t ) . Hence, simple calculations lead to Eqs. (15) and (26).
B Numerical schemes
We provide details about the numerical schemes employed to approximate the solutionof the models proposed in this paper.
B.1 RU dynamics
Since the equations describing the evolution of the RUs are independent of the vari-ables describing the XB states, their solution can be approximated independently ofthat of the models describing the XB dynamics. Specifically, we adopt a ForwardEuler scheme with a time step size of 2 . · − s. Further details on the properties ofthe numerical solution obtained by applying this scheme can be found in [79]. B.2 XB dynamics
Concerning the equations describing the XB dynamics, we notice that Eqs. (15) and(26) can be written in the following form: (cid:40) ˙ x ( t ) = A ( t ) x ( t ) + r ( t ) t ∈ (0 , T ] , x (0) = x , (33)where x ( t ) is the vector of the variables describing the states of XBs, while A ( t ) and r ( t ) are respectively a time-dependent matrix and vector, determined by the input v ( t ) and by the RUs states π αβδ,ϑηλi ( t ) (or, for mean-field models, π αβδ,η ( t )).In order to approximate the solution of Eq. (33), we consider a subdivision 0 = t < t < · · · < t M = T of the time interval [0 , T ] with time step size ∆ t and wedenote by x ( k ) ≈ x ( t k ) the approximated solution at time t k . Due to the linearity ofEq. (33), we consider the following exponential scheme [37]: x (0) = x , x ( k ) ∞ = − A − ( t k ) r ( t k ) for k ≥ , x ( k ) = x ( k ) ∞ + e ∆ t A ( t k ) ( x ( k − − x ( k ) ∞ ) for k ≥ . (34)43ue to the implicit nature of the scheme of Eq. (34), that entails better stabilityproperties than the explicit scheme used for the RUs equations [71], we solve it witha larger time step size than the one used for the RUs model (∆ t = 1 · − s), by afirst-order time splitting scheme. References [1] Y. Ait-Mou, K. Hsu, G. P. Farman, M. Kumar, M. L. Greaser, T. C. Irving,and P. P. de Tombe. “Titin strain contributes to the Frank–Starling law of theheart by structural rearrangements of both thin-and thick-filament proteins”.In:
Proceedings of the National Academy of Sciences
Nonlinear problems of elasticity . Springer, 1995.[3] L. Azzolin, L. Ded`e, A. Gerbi, and A. Quarteroni. “Effect of fibre orientationand bulk value on the electromechanical modelling of the human ventricles”.In:
Mathematics in Engineering 2020 (to appear) (2020).[4] P. Backx, W. Gao, M. Azan-Backx, and E. Marb´an. “The relationship betweencontractile force and intracellular [Ca2+] in intact rat cardiac trabeculae.” In:
The Journal of General Physiology
The elements of stochastic processes with applications to the naturalsciences . Vol. 25. John Wiley & Sons, 1990.[6] J. D. Bayer, R. C. Blake, G. Plank, and N. A. Trayanova. “A novel rule-basedalgorithm for assigning myocardial fiber orientation to computational heartmodels”. In:
Annals of Biomedical Engineering
Excitation-contraction coupling and cardiac contractile force . Vol. 237.Springer Science & Business Media, 2001.[8] J. Bestel, F. Cl´ement, and M. Sorine. “A biomechanical model of muscle con-traction”. In:
International conference on medical image computing and computer-assisted intervention . Springer. 2001, pp. 1159–1161.[9] E. Brunello, M. Caremani, L. Melli, M. Linari, M. Fernandez-Martinez, T.Narayanan, M. Irving, G. Piazzesi, V. Lombardi, and M. Reconditi. “The con-tributions of filaments and cross-bridges to sarcomere compliance in skeletalmuscle”. In:
The Journal of Physiology
Biophysical Journal
Proceedings of the National Academy of Sciences
Biomechanics and Modeling inMechanobiology
Reportson Progress in Physics
InterfaceFocus
InternationalJournal for Multiscale Computational Engineering
Complex systems inbiomedicine . Springer, 2006, pp. 187–241.[18] P. Colli Franzone, L. F. Pavarino, and S. Scacchi.
Mathematical cardiac elec-trophysiology . Vol. 13. Springer, 2014.[19] R. Craig and W. Lehman. “Crossbridge and tropomyosin positions observed innative, interacting thick and thin filaments”. In:
Journal of Molecular Biology
Exper-imental Physiology
TheMathematics of Mechanobiology . Ed. by D. Ambrosi and P. Ciarletta. Springer,2020 (to appear).[22] D. Dobesh, J. Konhilas, and P. de Tombe. “Cooperative activation in cardiacmuscle: impact of sarcomere length”. In:
American Journal of Physiology-Heartand Circulatory Physiology
Biophysical Journal
PLoS Computational Biology
Biophysical Journal
Progress in Biophysics and Molecular Biology
Circulation
Circulation Re-search
The Journal of Physiology
Mathematicsin Engineering
Journal of Elasticity (2019), pp. 1–20.[32] A. Gordon, A. F. Huxley, and F. Julian. “The variation in isometric tensionwith sarcomere length in vertebrate muscle fibres”. In:
The Journal of Physi-ology
Physiology
The Journal of General Physiology
Biochemical Jour-nal
Pro-ceedings of the Royal Society of London B: Biological Sciences
Acta Numer-ica
19 (2010), pp. 209–286.[38] P. Hunter, A. McCulloch, and H. Ter Keurs. “Modelling the mechanical prop-erties of cardiac muscle”. In:
Progress in Biophysics and Molecular Biology
IBM Journal of Researchand Development
Prog BiophysBiophys Chem (1957), pp. 255–318.[41] A. Huxley and R. Simmons. “Proposed mechanism of force generation in stri-ated muscle”. In:
Nature
Progress inBiophysics and Biophysical Chemistry
American Journal of Physiology -Heart and Circulatory Physiology
American Journalof Physiology - Heart and Circulatory Physiology
American Journal of Physiology-Heart and Circulatory Physiology
American Journal of Hypertension
Physiology of the heart . Lippincott Williams & Wilkins, 2010.[48] J. Keener and J. Sneyd.
Mathematical physiology . Vol. 1. Springer, 2009.[49] J. Kentish, H. Ter Keurs, and D. Allen. “The contribution of myofibrillar prop-erties to the sarcomere length-force relationship of cardiac muscle”. In:
Starlingslaw of the heart revisited . Springer, 1988, pp. 1–17.[50] J. Kentish, H. ter Keurs, L. Ricciardi, J. Bucx, and M. Noble. “Comparisonbetween the sarcomere length-force relations of intact and skinned trabeculaefrom rat right ventricle. Influence of calcium concentrations on these relations.”In:
Circulation Research
Starlings law of the heart revisited . Vol. 89. SpringerScience & Business Media, 2012.[52] F. Kimmig. “Multi-scale modeling of muscle contraction - From stochastic dy-namics of molecular motors to continuum mechanics”. PhD thesis. Universit´eParis-Saclay, 2019.[53] F. Kimmig, M. Caruel, P. Moireau, and D. Chapelle. “Activation-contractioncoupling in a multiscale heart model”. In:
Proceedings of CMBE 2019 (volume1) . 2019, pp. 96–99.[54] F. Kimmig, D. Chapelle, and P. Moireau. “Thermodynamic properties of mus-cle contraction models and associated discrete-time principles”. In:
AdvancedModeling and Simulation in Engineering Sciences
PLoS Computational Biology
Journal of Molecular andCellular Cardiology
106 (2017), pp. 68–83.[57] S. Land, S. A. Niederer, J. M. Aronsen, E. K. Espe, L. Zhang, W. E. Louch,I. Sjaastad, O. M. Sejersted, and N. P. Smith. “An analysis of deformation-dependent electromechanical coupling in the mouse heart”. In:
The Journal ofPhysiology
American Journal ofPhysiology-Heart and Circulatory Physiology
Biochemistry
Physical Review E
The European Physical Journal E
Biophysical Journal
Journal ofmolecular and cellular cardiology
127 (2019), pp. 11–19.[64] J. Norris.
Markov chains . 2. Cambridge University Press, 1998.[65] R. W. Ogden.
Non-linear elastic deformations . Courier Corporation, 1997.[66] B. Pan, A. Gordon, and Z. Luo. “Removal of tropomyosin overlap modifiescooperative binding of myosin S-1 to reconstituted thin filaments of rabbitstriated muscle.” In:
Journal of Biological Chemistry
Biophysical Journal
Biomechanics and Modelingin Mechanobiology
P߬ugers Archiv-European Journal ofPhysiology
P߬ugers Archiv
Numerical mathematics . Vol. 37. SpringerScience & Business Media, 2010.[72] A. Quarteroni, T. Lassila, S. Rossi, and R. Ruiz-Baier. “Integrated Heart–Coupling multiscale and multiphysics models for the simulation of the cardiacfunction”. In:
Computer Methods in Applied Mechanics and Engineering
Numerical approximation of partial differentialequations . Vol. 23. Springer Science & Business Media, 2008.[74] A. Quarteroni, L. Dede’, A. Manzoni, and C. Vergara.
Mathematical modellingof the human cardiovascular system: data, numerical approximation, clinicalapplications . Cambridge Monographs on Applied and Computational Mathe-matics. Cambridge University Press, 2019.4875] M. Razumova, A. Bukatina, and K. Campbell. “Stiffness-distortion sarcomeremodel for muscle simulation”. In:
Journal of Applied Physiology
Biomechanics and Modeling in Mechanobiology (2018), pp. 1–24.[77] F. Regazzoni, L. Ded`e, and A. Quarteroni. “Active force generation in cardiacmuscle cells: mathematical modeling and numerical simulation of the actin-myosin interaction”. In:
MOX Report 2019/45, Politecnico di Milano (2019).[78] F. Regazzoni, L. Ded`e, and A. Quarteroni. “Machine learning of multiscaleactive force generation models for the efficient simulation of cardiac electrome-chanics”. In:
MOX Report 33/2019, Politecnico di Milano (2019).[79] F. Regazzoni. “Mathematical modeling and Machine Learning for the numericalsimulation of cardiac electromechanics”. PhD thesis. Politecnico di Milano,2020.[80] D. Riccobelli and D. Ambrosi. “Activation of a muscle as a mapping of stressstraincurves”. en. In:
Extreme Mechanics Letters
28 (Apr. 2019), pp. 37–42.[81] J. Rice and P. de Tombe. “Approaches to modeling crossbridges and calcium-dependent activation in cardiac muscle”. In:
Progress in Biophysics and Molec-ular Biology
American Journal of Physiology-Heart and Circulatory Physiology
Biophysical Journal
Bio-physical Journal
Computational cardiology: modeling of anatomy, electrophysiology,and mechanics (lecture notes in Computer Science) . Secaucus, NJ, USA: Springer-Verlag New York, Inc., 2004.[86] F. Sachse, K. Gl¨anzel, and G. Seemann. “Modeling of protein interactions in-volved in cardiac tension development”. In:
International Journal of Bifurcationand Chaos
MOXReport 2019/36, Politecnico di Milano (2019).[88] V. Sequeira and J. van der Velden. “The Frank–Starling Law: a jigsaw of titinproportions”. In:
Biophysical Reviews
Annals of Biomedical Engi-neering
Acta Numerica
13 (2004), pp. 371–431.4991] S. Sugiura, T. Washio, A. Hatano, J. Okada, H. Watanabe, and T. Hisada.“Multi-scale simulations of cardiac electrophysiology and mechanics using theUniversity of Tokyo heart simulator”. In:
Progress in Biophysics and MolecularBiology
The Journal ofPhysiology
American Journal of Physiology-Heartand Circulatory Physiology
Circulation Research
Skeletal muscle mechanics: from mechanics to function. Wiley, New York (2000), pp. 53–70.[96] H. Ter Keurs, T. Shinozaki, Y. Zhang, M. Zhang, Y. Wakayama, Y. Sugai,Y. Kagaya, M. Miura, P. Boyden, B. Stuyvers, et al. “Sarcomere mechanics inuniform and non-uniform cardiac muscle: a link between pump function andarrhythmias”. In:
Progress in biophysics and molecular biology
Journal of Molecular and Cellular Cardiology
91 (2016),pp. 148–150.[98] K. Tøndel, S. Land, S. A. Niederer, and N. P. Smith. “Quantifying inter-speciesdifferences in contractile function through biophysical modelling”. In:
The Jour-nal of Physiology
Computing and Visualizationin Science
Cellularand Molecular Bioengineering
Multiscale Modeling & Simulation
InternationalJournal for Numerical Methods in Biomedical Engineering (2015).[103] N. Westerhof, J.-W. Lankhaar, and B. E. Westerhof. “The arterial Windkessel”.In:
Medical & Biological Engineering & Computing
Mathematical Biosciences
Biophysical Journal
Zygote 3D models . 2019. url :