http://orcid.org/0000-0002-4207-1400
Figures
Abstract
We propose four novel mathematical models, describing the microscopic mechanisms of force generation in the cardiac muscle tissue, which are suitable for multiscale numerical simulations of cardiac electromechanics. Such models are based on a biophysically accurate representation of the regulatory and contractile proteins in the sarcomeres. Our models, unlike most of the sarcomere dynamics models that are available in the literature and that feature a comparable richness of detail, do not require the time-consuming Monte Carlo method for their numerical approximation. Conversely, the models that we propose only require the solution of a system of PDEs and/or ODEs (the most reduced of the four only involving 20 ODEs), thus entailing a significant computational efficiency. By focusing on the two models that feature the best trade-off between detail of description and identifiability of parameters, we propose a pipeline to calibrate such parameters starting from experimental measurements available in literature. Thanks to this pipeline, we calibrate these models for room-temperature rat and for body-temperature human cells. We show, by means of numerical simulations, that the proposed models correctly predict the main features of force generation, including the steady-state force-calcium and force-length relationships, the length-dependent prolongation of twitches and increase of peak force, the force-velocity relationship. Moreover, they correctly reproduce the Frank-Starling effect, when employed in multiscale 3D numerical simulation of cardiac electromechanics.
Author summary
Computer-based numerical simulations of the heart are increasingly assuming a recognized role in the context of computational medicine and cardiology. They are based on mathematical models describing the different physical phenomena occurring during an heartbeat. Among these models, a pivotal role is played by those describing how cardiomyocytes—the cardiac muscle cells—produce active force, driven by changes in calcium concentration. However, due to the intrinsic complexity of these subcellular mechanisms, the computational cost associated with the solution of cardiac active force models is often prohibitive. For this reason, phenomenological models are typically used in place of biophysically detailed ones in organ-scale simulations. In this paper, we propose some new biophysically detailed mathematical models of cardiac force generation. Our models are rigorously derived on the basis of physically motivated assumptions that allow to drastically reduce the computational cost associated to their resolution, making them suitable for organ-scale numerical simulations, without renouncing to their biophysical detail.
Citation: Regazzoni F, Dedè L, Quarteroni A (2020) Biophysically detailed mathematical models of multiscale cardiac active mechanics. PLoS Comput Biol 16(10): e1008294. https://doi.org/10.1371/journal.pcbi.1008294
Editor: Daniel A. Beard, University of Michigan, UNITED STATES
Received: April 20, 2020; Accepted: August 27, 2020; Published: October 7, 2020
Copyright: © 2020 Regazzoni et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: The codes implementing the proposed models are available in an open-source repository at https://github.com/FrancescoRegazzoni/cardiac-activation. The datasets accompanying the article are available at https://doi.org/10.5281/zenodo.3992553.
Funding: This project has received funding from the European Research Council (ERC, https://erc.europa.eu/) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 740132, iHEART - An Integrated Heart Model for the simulation of the cardiac function, P.I. Prof. A. Quarteroni). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
This is a PLOS Computational Biology Software paper.
Introduction
Cardiovascular mathematical and numerical models are increasingly used, with a twofold role [1–6]. On the one hand, realistic and detailed in silico models of the heart can deepen the understanding of its function, help the interpretation of experimental observations and explain the delicate links between the organ-level emergent phenomena and the underlying biophysical mechanisms. On the other hand, patient-specific numerical simulations, which are increasingly becoming available, can provide clinicians with valuable quantitative information to improve patient care and to support decision-making procedures.
The construction of an integrated mathematical and numerical model of cardiac electromechanics (EM) is however a remarkably arduous task. This is mainly due to the multiphysics (due to the interplay of biochemistry, electricity, solid mechanics, fluid dynamics) and multiscale nature of the heart: characteristic spatial scales range from nanometers to centimeters and the temporal ones from microseconds to seconds. This makes it difficult to devise computationally efficient and accurate algorithms for a plurality of mathematical models featuring a broad degree of details [5, 7–11].
The contrasting needs between model accuracy and computational efficiency of numerical simulations is mainly due to the multiscale nature of the heart, for which the mechanical work enabling the macroscopic motion of the organ is fueled by the energy consumed at the microscale by subcellular mechanisms. The generation of active force takes place inside sarcomeres and involves a complex chain of chemical and mechanical reactions. This mechanism can be split into two steps, that we sketch in Fig 1. First, a ionic signal (specifically, an increase of calcium ions concentration) triggers the so-called regulatory units, protein complexes consisting of troponin and tropomyosin, that act as on-off switches for the muscle contraction. Then, when the regulatory units are activated, the actin and myosin proteins are free to interact and form the so-called crossbridges, molecular motors that generate an active force by consuming the chemical energy stored in ATP [12, 13].
Microscopic force generation includes many regulatory mechanisms, forming the subcellular basis of organ-level phenomena, such as the Frank-Starling effect [13]. Hence, if a microscale mathematical model of force generation is used in a multiscale setting to build an integrated organ-level EM model, then it should be able to reproduce the above-mentioned mechanisms.
In the past decades, several efforts have been dedicated to the construction of mathematical models describing the complex dynamics of the processes taking place in sarcomeres [14–29]. However, because of the intrinsic complexity of the phenomenon of force generation, huge computational costs are associated with the numerical approximation of such models, thus limiting their application within multiscale EM simulations. Despite several attempts to capture the fundamental mechanisms underlying the force generation phenomenon into a tractable number of equations [16, 17, 30–34], the existing organ-level cardiac mathematical models rely on two alternative strategies to describe microscopic force generation.
- Phenomenological models (see e.g. [18, 23, 26, 35, 36]) are built by fitting the measured data with simple laws, chosen by the modeler. However, the parameters characterizing phenomenological models often lack a clear physical interpretation; moreover, the noisy nature and deficiency of data coming from the subcellular contractile units and the intrinsic difficulties in measuring sarcomeres under the conditions occurring during an heartbeat hamper the predictive power of such models [4].
- Biophysically detailed models are based on an accurate description of the proteins involved in the force generation process and are derived from physics first principles. However, their numerical solution, because of their complexity, is typically obtained by means of a Monte Carlo (MC) approximation (see e.g. [19, 24, 25]). The MC method is in fact inefficient, featuring a huge computational cost, both in terms of time and memory storage. Indeed, to accurately approximate the solution of a single heartbeat for a single myofilament, tens of hours of computational time may be required; see e.g. [37].
The purpose of this paper is to develop a biophysically detailed model for active force generation, that explicitly describes the fundamental ingredients of the force generation apparatus, yet with a tractable computational cost, so that it is suitable for multiscale EM simulations.
Paper outline
This paper is structured as follows. First, we recall the main features of the force generation phenomenon in cardiomyocytes and the main difficulties encountered in the construction of mathematical models describing the associated mechanisms. In the section Models, we present the models proposed in this paper and we describe the strategy employed for their calibration. Then, in the section Results and Discussion, we show some numerical results obtained with the proposed models, including filament-level 0D simulations and multiscale 3D cardiac EM simulations. Finally, we provide some concluding remarks. In Table 1 we provide a list of the abbreviation used throughout this paper.
Microscale models of cardiac contraction
We recall the microscopical mechanisms by which active force is generated in the cardiac tissue and we highlight the difficulties, rooted in the their intrinsic complexity, in describing such phenomena with a tractable number of equations. We also review the main contributions available in literature.
Sarcomeres are cylindrically-shaped, 2 μm length units made of a regular arrangement 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).
In the next sections, first we deal with the activation of the thin filament, involving the troponin-tropomyosin regulatory units (RUs). Then, we address the crossbridge (XB) cycling and finally we consider the full-sarcomere dynamics (see Fig 1).
Modeling the thin filament activation.
The activation of the thin filament is mainly regulated by two variables, namely the intracellular calcium ions concentration ([Ca2+]i) and the sarcomere length (SL). The experimental signature of the regulation mechanisms is given by the steady-steady relationships between calcium, sarcomere length and generated force.
The force-calcium relationship (see Fig 2) features a sigmoidal shape, well described by the Hill equation [40–42]:
(1)
where
is the maximum tension (tension at saturating calcium levels), EC50 is the half maximal effective concentration (i.e. the calcium concentration producing half maximal force) and nH is the Hill coefficient. The experimentally measured force-calcium curves in the cardiac tissue feature a steep slope in correspondence of half activation (Hill coefficient greater than one) [41–43], thus revealing the presence of cooperative effects. Despite several explanations have been proposed [16, 44–50], the most likely hypothesis lies in the end-to-end interactions of Tm units [51–53].
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 and the 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 the increase of extension of the single-overlap zone (i.e. the region of the sarcomere where the MF filament a single AF), a well-assessed explanation for the second effect (known as length-dependent activation, LDA) has not yet been found [44, 50, 53–60].
The earliest attempts to model the calcium-driven regulation of the muscular contractile system date back to the 1990s [16, 30, 31, 61–63]. Those models rely on the formalism of continuous-time Markov Chains (CTMC), also known as Markov Jump processes (see e.g. [64]), to model the transitions between the different configurations assumed by the proteins involved in the force regulation process. In those models, the necessity of representing the end-to-end interactions of Tm units dictates a spatially-explicit representation of the RUs. Indeed, mean-field models, where only a single representative RU is considered, fail to correctly predict the cooperative activation and the resulting steep force-calcium curves [17, 44], unless phenomenological laws are introduced, as discussed below. However, the number of degrees of freedom of the CTMC increases exponentially with the number of RUs represented. This hinders the possibility of numerically approximating the solution of the Forward Kolmogorov Equation, also known as Master Equation in natural sciences, which describes the time evolution of the probabilities associated with the states of a stochastic process [65]. As a matter of fact, the Forward Kolmogorov Equation associated with this CTMC would have a number of variables that is exponential in the number of RUs, resulting in as many as 1020 or more variables [37]. For this reason, spatially-explicit models require a MC approximation for their numerical resolution, thus resulting in very large computational costs [33, 37, 44].
To avoid an explicit representation of end-to-end interactions, phenomenological models, where the transition rates are set as a nonlinear functions of the calcium concentration, have been proposed [18, 23, 26, 35, 36]. These models are however based on phenomenological laws. Alternatively, to overcome the large computational cost induced by the MC method without renouncing to represent end-to-end interactions (unlike in phenomenological models), several attempts to capture cooperative phenomena by means of numerically tractable ODE systems have been done in literature. In [17] an analytical solution is derived for the steady-state. In [32], a periodicity assumption is used to reduce the number of unknowns. In [33] each RU is considered independently of each others, while end-to-end interactions are accounted for by fitting the parameters of an integro-differential system with memory from a collection of simulations. In [34], the states of the CTMC are grouped by the number of unblocked RUs and a MC sampling technique is used to estimate the average free energy of each group and, thus, the transition rates within groups. For further details on these modeling attempts, the interested reader can refer to [37, 66].
Modeling the crossbridge dynamics.
Active force is generated by XBs by the cyclical attachment and detachment of myosin heads (MHs) to actin binding sites (BSs). When MHs are in their attached configuration, 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 as Lymn-Taylor cycle [67], features a wide range of time scales (nearly from 1 to 100 ms) [12, 28, 68]; hence, also the response of the force generation apparatus to external stimuli is characterized by different time scales. Indeed, when a fast step in force (respectively, in length) is applied to an isometrically contracted muscle fiber, three different phases can be observed [20, 28, 68–70]. First, an instantaneous elastic response occurs along the so-called T1-L1 curve (see Fig 2), whose slope corresponds to the stiffness of the attached myosin proteins. Then, a fast transient (2-3 ms) occurs, and the sarcomere reaches a length (respectively, a tension) belonging to the T2-L2 curve (see Fig 2). Such second phase corresponds to a rearrangement of MHs within their pre- and post-power-stroke configuration. Finally, with a time scale on nearly 100 ms, the muscle fiber reaches a steady-state regime, characterized (in case the step is applied by controlling the force) by a constant shortening (or lengthening) velocity. The relationship between the steady-state velocity and the muscle tension constitutes the so-called force-velocity curve (see Fig 2), firstly measured by Hill in [71], and it is characterized by a finite value of velocity (denoted by vmax) for which the generated tension is zero [12, 13]. The experimental measurements of the T1-L1, the T2-L2 and the force-velocity curves are invariant after normalization of the tension by its isometric value, denoted by . This reveals that the underlying mechanisms are related to the XB dynamics, while they are independent of the thin filament regulation, whose effect is simply that of tuning the number of recruitable XBs [12, 13, 70].
The attachment-detachment process of MHs has been described accordingly with the formalism of the Huxley model [14] (that we denote by H57 model), where the population of MHs is described by the distribution density of the distortion of attached XB. The time evolution of such distribution is driven by a PDE, where a convection term 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 destruction of XBs [22, 72–74]. In order to capture the separation between the fastest time scales (i.e. between the first two phases following a fast step either in force or in length), an explicit representation of the power-stroke must be included in the model, by introducing a multistable discrete [15, 75] or continuous [20, 28, 29, 69, 76] degree of freedom, representing the angular position of the rotating MH.
Modeling the full sarcomere dynamics.
In the past two decades, several models describing the generation of active force in the 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 lies in the spatial dependence of the cooperativity phenomenon, crucial to reproduce the calcium dependence of muscle activation. As a matter of fact, an explicit representation of spatial-dependent cooperative mechanisms dramatically increases the computational complexity of activation models, even more so when such models are coupled with models describing XB cycling. When the interactions between BSs and MHs are considered, indeed, one must face the further difficulty of tracking which BS faces which MH when the filaments mutually slide. The attempt of capturing such spatially dependent phenomena in a compact system of ODEs is the common thread of most of the literature on sarcomere modeling (see e.g. [17, 22, 32–34, 36, 37, 77–79]). We remark that computational efficiency is a major issue when sarcomere models 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 points of the computational domain (typically at each nodal point of the computational mesh). Nonetheless, most of biophysically detailed full-sarcomere models rely on the time-consuming MC method for their numerical approximation [19, 24, 25, 33, 80].
Models
In this section, we propose four different microscale models of active force generation in the cardiac tissue. These models are derived from a biophysically detailed CTMC, accurately describing the dynamics of the regulatory and contractile proteins. We present a strategy to derive a set of equations, with a dramatically smaller number of variables than the Forward Kolmogorov Equation, describing the evolution of the biophysically detailed CTMC. Our strategy is based on physically motivated assumptions, aimed at neglecting second-order interactions among the proteins, focusing only on the main interactions, so that the variables describing the stochastic processes associated with the states of the proteins can be partially decoupled, similarly to what we have done in [37]. This results in a drastic reduction in the size of models. Moreover, we change the classical MH-centered description of XBs (see e.g. [15, 29, 72, 75]), in favor of a BS-centered one. This prevents the necessity of tracking the mutual sliding between the filaments, still without the need of introducing any simplifying assumption.
As in most of RUs models (see Introduction), we describe Tn and Tm by discrete states. Moreover, based on the experimental evidence that cooperativity is linked to RUs end-to-end interactions [51–53], we include nearest-neighbor interactions among RUs with the formalism of the model of [17].
Concerning the modeling of XBs, we are here interested in developing a model of cardiomyocytes contraction in the heart, which is characterized by slower time-scales than the fast time-scale of the power-stroke. This suggests that the level of detail that best suits the application to cardiac EM does not require to explicitly represent the power-stroke [81]. In [29], indeed, the authors showed that, if the time-scales of interest are slower than the time-scale of the power-stroke, the detailed models including the power-stroke reduce to H57-like models, where only the attachment-detachment process of XBs is explicitly represented. Therefore, we model the XB dynamics as a two-states process, within the H57 formalism, where the attachment-detachment rates depend on the myosin arm distortion.
Model setup
We consider a pair of interacting myofilaments and we denote by NA the number of RUs located on an AF and by NM the number of MHs located on half MF. To identify a RU we employ the index , while to identify MHs we employ the index
. The CTMC model describing the dynamics of the RUs and of the MHs is sketched in Fig 3.
Each RU is described by a 4-state model (top left), whose dynamics is affected by the state nearest-neighboring RUs. An example of nearest-neighbor interactions is shown in the bottom-left box, where the notation for the transition rates kT,i is illustrated by the orange arrows. MHs are described as 2-state elements, whose transition rates are affected by the XB elongation, the sliding velocity between myofilaments and the permissivity state of the RU.
Each RU is composed by a Tn unit and by a Tm unit, respectively associated with the variables and
. Tn can be either unbound (
) or bound (
) to calcium (we write
and
, respectively). On the other hand, Tm can be either in non-permissive (
) or in permissive (
) configuration (we write
and
, respectively). In our model, the calcium binding kinetics is affected by the state of Tm. Hence, when the i-th Tm unit is non-permissive (i.e.
), we denote the binding and unbinding calcium rates by
and
, respectively; conversely, when
, we denote the binding and unbinding calcium rates by
and
. Similarly, the kinetics of Tm is affected by the calcium-binding state of the corresponding Tn unit. Moreover, because of the Tm end-to-end interactions, the Tm transition rates also depend on the state of the nearest-neighboring Tm units. Hence, the transition rates from the non-permissive to the permissive states and vice versa—given the state of Ti−1, Ti+1 and Ci—are respectively denoted by
and
. To better clarify the notation, an example is shown in Fig 3 (bottom-left box). In the example,
denotes the transition rate for i-th Tm unit from the non-permissive to the permissive state, when the nearest-neighboring units are both permissive and the associated Tn unit is bound to calcium. Concerning the Tm units located at the ends of the filaments, for which a neighbor is missing, the latter is assumed to be in state
. We exclude any feedback from XBs on the dynamics of the RUs, as recent experimental evidence suggests that this kind of feedback is not present [49, 50].
Each myosin arm is modeled as a linear spring with stiffness kXB. The attachment and detachment rates of XBs depend on the distance between the MH resting position and the BS, denoted by x, and on the relative velocity between the myofilaments, denoted by . For simplicity, we define as v(t) = 2vhs(t)/SL0 the normalized shortening velocity, where SL0 denotes a reference sarcomere length. Moreover, to model the calcium-driven regulation, the XB attachment and detachment rates depend on the state of the corresponding Tm unit. Hence, when
, we denote the XB binding and unbinding rates by
and
, respectively; conversely, when
, we denote the XB binding and unbinding rates by
and
. In particular, we set
since new XBs can form only if the corresponding Tm unit is permissive. Clearly, a XB can form only if neither the BS nor the MH is already attached to another site. Moreover, the attachment rates
and
are zero sufficiently far from x = 0.
In order to describe the state of XBs, we introduce the variables ,
and
, respectively denoting the state of actin BSs, the state of MHs and the displacement of attached XBs. Specifically, when the i-th actin BS is attached to the j-th MH we write
. Similarly, when the j-th MH is attached to the i-th actin BS we write
. Moreover, we write
whenever the i-th actin BS is attached to a MH with displacement x. Clearly, these variables are redundant, as we have:
where we have denoted by dij(t) the distance between the i-th BS and the j-th MH.
To summarize, the CTMC is described by the following stochastic processes, for ,
and t ≥ 0:
(2)
We remark that we denote the detached state by rather than
, because the latter notation is employed to denote the case when the i-th actin BS is attached with displacement x = 0.
The total force exerted by the pair of interacting half MF and AF is given by the sum of the force generated by each attached XB. Therefore, the expected value of the force is given by:
where FXB(x) denotes the force exerted by an attached MH with distortion x and where we set by convention FXB(∅) = 0. Here and in what follows, we denote by
the expected value of a random variable.
The size of the CTMC (2), that is the number of its states, is overwhelming. As a matter of fact, each RU can be in four possible states (,
,
and
), and the corresponding BS can be either unbound or bound to one of the NM MHs. In conclusion, the number of states of the CTMC is
. Thus, the numerical solution of the associated Forward Kolmogorov Equation is unaffordable, as it would feature as many variables as the number of states of the CTCM [37]. To overcome this inconvenient, we introduce some physical motivated assumptions that allow to partially decouple the dynamics of the stochastic processes, thus yielding a much reduced set of equations.
Models equations
In this section, we present equations describing the evolution of the stochastic processes of Eq (2). A detailed derivation of these equations is provided in the Supporting Information (S1 Appendix).
Thin filament regulation.
Due to the lack of feedback from XBs to RUs, it is possible to write an equation describing the evolution of the stochastic processes and
independently of the stochastic processes associated with the XBs (while the converse clearly does not hold). Similarly to [37], we focus on the joint probabilities of triplets of consecutive RUs. Hence, we consider the following variables, where i = 2, …, NA − 1,
and
:
(3)
where
denotes the probability of an event. For instance,
denotes the probability that the triplet centered in the i-th unit has the Tm units in states
,
and
and the Tn units in states
,
and
respectively, as shown in Fig 4.
The arrows illustrate the meaning of the notation.
For each i = 2, …, NA − 1, we have 64 variables written in the form (3), corresponding to as many states of the triplet. The dynamics of each variable is determined by 6 possible forward and backward transitions, as depicted in Fig 5. However, the transition rates associated with the Tm units at the edge of the triplet cannot be computed from the variables of Eq (3), as they depend on the state of a Tm unit outside the triplet. For instance, the rate of the transition from to
depends on the state of
, which does not belong to the triplet. Nonetheless, under a suitable hypothesis, the transition rates between Tm being in permissive and non-permissive state can be defined as:
(4)
where the symbol ° recalls that the corresponding unit has an arbitrary state. In Eq (4) and in what follows, we use the notation
,
,
and
to denote opposite states. For instance, if
, then
.
For example, this figure shows the transitions of the configuration associated with the state variable , with the corresponding transition rates. The transition rates computed thanks to Ass. (H1) are highlighted with a colored background.
Eq (4) is rigorously derived in the Supporting Information (S1 Appendix). The derivation is rather technical and is based upon the following assumption:
(H1)
where, given three events A, B, C we write A ⫫ B|C if A and B are conditionally independent given C (i.e.
). Assumption (H1) states that distant RUs are conditionally independent given the states of the intermediate RUs. From a modeling viewpoint, this means that the interaction of far units is mediated by the intermediate ones, which is coherent with the short-range nature of end-to-end interactions.
In conclusion, we obtain the following system of ODEs, for t ≥ 0 and for any i = 2, …, NA − 1, and
:
(5)
endowed with suitable initial conditions.
The permissivity of a given RU is defined as its probability of being in permissive state (i.e. ) and can be obtained from the variables
as:
Crossbridge dynamics.
Similarly to the H57 model, we introduce distribution density functions tracking the elongation of attached XBs. However, since in our CTMC the XB transition rates depend on the state of the associates Tm unit, we split attached XBs into two families: those associated with a non-permissive Tm and those associated with a permissive one. Hence, we define the following variables, for , corresponding to the probability density that the i-th BS is attached to a MH with displacement x and that the associated RU is in a given permissivity state:
(6)
where the symbol
denotes a probability density function.
We notice that we make here the choice of tracking the XBs from the point of view of the BSs, rather than of the MHs, as it is traditionally done in literature [15, 29, 72, 75]. This change of perspective has the significant advantage that it does not require to track which RU faces which MH at each time. Indeed, each BS and each RU, being located on the same filament, rigidly move with respect of each others and, thus, each BS is always associated with the same RU.
In the Supporting Information (S1 Appendix) we derive the following system of PDEs, describing the time evolution of and
:
(7)
endowed with suitable initial conditions, where we define:
(8)
The sink and source terms in Eq (7) account for the fluxes among the two groups (XBs with non-permissive Tm and with permissive Tm). The terms Pi and (1 − Pi) represent the maximum possible proportion of attached BSs in each group. The term DM, defined as the distance between two consecutive MHs, appears because in this setting and
are, from a dimensional point of view, the inverse of length units (they are probability densities), whereas the variables of the H57 model are dimensionless. We remark that, differently than the H57 model, Eq (7) is referred to BSs rather than to MHs.
The derivation of Eq (7), presented in the Supporting Information (S1 Appendix), is based on two assumptions. First, we assume that the state of a BS is conditionally independent of the state of surrounding RUs, given the permissivity state of the associated RU. This is coherent with the physics of the model, as the only feature of the RUs that directly affects the XBs binding rates is the permissivity state of Tm. In mathematical terms, this assumption reads:
(H2)
Second, we assume that the shortening velocity v is never so large to convect attached BSs within the range of attachment of a different MH. In mathematical terms, we assume that:
(H3)
The latter assumption allows to decouple the dynamics of the different units. We notice that all the models belonging to the family of the H57 model are based on assumptions analogous to Ass. (H3), without which the H57 equation cannot be derived.
By combining Eq (5) with Eq (7), 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 is described by a PDE system). In this model, the expected value of the force exerted by the whole half filament can be determined as follows:
(9)
Distribution-moments equations.
When the XB attachment-detachment transition rates assume special forms, the PDE system of Eq (7) can be reduced to a more compact system of ODEs, by following a general strategy in statistical physics, already used for H57-like models [22, 77, 78]. Specifically, under suitable hypotheses on the transition rates, the distributions of the elongation of attached XBs (i.e. and
) can be fully characterized by their first two moments that we denote by
,
and by
,
, respectively. More precisely, we define, for
, for
and for p = 0, 1:
(10)
We notice that the zero order moment (respectively,
) can be interpreted as the probability that the i-th BS is attached and associated to a non-permissive (respectively, permissive) RU. Moreover, the ratio
(respectively,
) corresponds to the expected value of the distortion (normalized with respect to SL0/2) of the XB attached to the i-th RU, provided that the corresponding RU is in non-permissive (respectively, permissive) state. The physical meaning of the variables
and
is illustrated in Fig 6.
The variable (respectively,
) corresponds to the fraction of permissive (respectively, non-permissive) BSs to which a MH is bound, while the ratio
(respectively,
) corresponds to the average XB distortion within the permissive (respectively, non-permissive) attached BSs.
Let us assume that the total attachment-detachment rate is independent of the XB distortion (i.e. there exist functions and
, for
, such that
and
for any
). Under this assumption, as shown in the Supporting Information (S1 Appendix), we get the following distribution-moments equations:
(11)
By assuming a linear spring hypothesis for the tension generated by attached XBs (i.e. FXB(x) = kXB x), the expected value of the force of half filament is given by:
(12)
where we have defined:
(13)
Since the active force generated by half MF is proportional to μ1(t), there exists some constant aXB (with the dimension of a pressure) such that the macroscopic active tension can be written as Ta(t) = aXB μ1(t). In the following, we denote by SE-ODE model the one obtained by combining Eq (5) with Eq (11) (where ODE denotes the fact that the XB dynamics is described by an ODE system).
Mean-field approximation.
The models proposed so far are based on an explicit spatial description of the physical arrangements of proteins along the myofilaments. The spatial description allows to model the cooperativity mechanism (linked to the nearest-neighbor interactions within RUs) and the SL related effects on the force generation machinery (linked to the filament overlapping). However, the first phenomenon, despite being spatially dependent, is based on interactions of local type; the effect of the second phenomenon, in turn, largely depends on the size of the single-overlap zone, that is a scalar quantity non dependent on the spatial variable. Based on the above considerations, we propose a mean-field approximation of the spatially dependent CTMC presented above, where the nearest-neighbor interaction are captured as a local effect, and the effect of SL is modeled in function of the size of the single-overlap zone.
This mean-field model is based on the assumption that the single-overlap zone can be considered as infinitely long. Such approximation is reasonable as far as the effect of the edges can be neglected (the validity of such approximation will be discussed in the Conclusions). A direct consequence of this assumption is the invariance by translation of the joint distribution of RUs. In other terms, the variables defined in Eq (3) coincide for each i. In addition, we further reduce the number of variables by tracking only the state of the Tn unit of the central RU of the triplet (this further reduction is made possible by the fact that we never have to track the behavior at the boundaries of the filaments, as we will see in what follows). We define thus the following variables, for
and
:
(14)
A visual representation of these variables is provided in Fig 7. We notice that the variables παβδ,η(t) are well-defined thanks to the translational invariance of the distribution of RUs. Moreover, the transition rates and
for the units in the single-overlap region do not depend on i. Hence, we will denote them simply as
and
.
The arrows illustrate the meaning of the notation.
The dynamics of the variables παβδ,η(t) involves 4 possible forward and backward transitions (see Fig 8). Similarly to the spatially-explicit model, the transitions rates associated with the edge Tm units cannot be determined without additional assumptions. Therefore, we introduce the following assumption:
(H4)
The transition rates computed thanks to Ass. (H4) are highlighted with a colored background.
Assumption (H4), similarly to (H1), states the conditional independence of far units, given the states of the intermediate ones. In other terms, by Assumption (H4) we neglect the long-range interactions between Tm units, which is a secondary effect compared to short-range (i.e. end-to-end) interactions. We remark that Ass. (H4) and (H1) are two different mathematical translations of the same physical concept. The difference is due to the different state representation done in the two models: while the state of the SE model comprises three Tn units, the one of the MF model comprises only one Tn unit.
In this way we obtain (as proved in the Supporting Information (S1 Appendix)), the following ODE model, valid for t ≥ 0 and for any and
:
(15)
where:
The permissivity of a RU in the single-overlap zone, defined as (such that the i-th RU belongs to the single-overlap zone), can be obtained as:
By similar arguments, it follows that also the joint distribution of the stochastic processes associated with XB formation enjoys the translational invariance property and, consequently, the following variables are well defined, as the right-hand sides are independent of the index i (for i belonging to the single-overlap zone):
(16)
By proceeding as before, we get the following model:
(17)
where we have defined:
(18)
The expected value of the force exerted by the whole half filament can be obtained as follows:
(19)
where the single-overlap ratio χso denotes the fraction of the AF filament in the single-overlap zone:
(20)
LA being the length of the AF, LM the length of the MF and LH the length of the bare zone (see Fig 9). We notice that we are here assuming that the relative sliding between the filaments is such that χso slowly varies, so that we can neglect the effects linked to the state transitions taking place at the border of the single-overlap zone. The combination of Eqs (15) and (17) gives a model for the full-sarcomere dynamics, which we denote as the MF-PDE model (where MF stands for mean-field).
The thick filament (MF) is represented in red and two thin filaments (AF) are represented in blue (top). The origin of the frame of reference is the left side of the reference AF. The functions χSF and χM are also represented (bottom).
Moreover, under the assumption that the total attachment-detachment rate does not depend on the XB elongation (i.e. there exist two functions and
such that
and
for any
), we can derive the following distribution-moment equation:
(21)
endowed with suitable initial conditions and where we define, for
:
(22)
The force exerted by half thick filament is then given by:
(23)
where
(24)
for p = 0, 1 and, hence, the tissue level active tension is Ta(t) = aXB μ1(t). Finally, by combining Eqs (15) and (21) we obtain a model that we denote as MF-ODE model.
List of proposed models.
Table 2 provides a recap of the different models proposed in this paper. For each model, we report the assumptions underlying its derivation. The two models derived within the distribution-moments formalism (SE-ODE and MF-ODE) also require that the sum of 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 modeling choice.
For future reference, we assign a name to each model (SE stands for spatially-explicit,MF stands for mean-field). In the second column we report the number of ODEs and PDEs included in each model as a function of NA and NM and we specify the resulting values in the case NA = 32, NM = 18. In the “Assumptions” column,m.f. stands for mean-field assumption.
Table 2 contains four different models describing the same biological phenomenon. A natural question is when each model is to be preferred with respect to the others and which are its advantages and disadvantages. The modeler should make two choices:
- ODE versus PDE models. The two PDE models leave a great freedom to the modeler in the choice of the XB transition rate functions
,
,
and
. Hence, if the modeler has knowledge of the specific form of these functions (or if he wants to test different choices), PDE models are to be preferred. Otherwise, the two ODE models allow to avoid having to define the precise form of the functions
,
,
and
, only requiring to set a few scalar parameters, that can be easily calibrated from macroscopic measurements, as we show later.
If the interest of the modeler is oriented towards the microscopical dynamics of XBs, PDE models allow to simulate the precise distribution of XB strains, whereas ODE models only describe its first two moments. The greater detail of PDE models, however, comes at the price of larger computational costs. Therefore, in the context of multiscale simulations such as cardiac electromechanics—where the model should be solved simultaneously in a large number of points—the ODE models are computationally more attractive than PDE models. - SE (spatially-explicit) versus MF (mean-field) models. The two MF models are derived from the corresponding SE models by considering a single triplet of RUs rather than a whole AF. Hence, the former models neglect the profile that the states of RUs assume alongside the AF. As we show later, this profile—specifically the behavior near the end-points of the single-overlap zone—may play a role in the SL-induced change in calcium sensitivity (the so-called LDA), but the sources of this phenomenon have not been not fully understood yet [53, 57, 58, 60]. Clearly, the MF models are much lighter than the SE ones. Hence, the former should be preferred when computational cost is a major issue, such as in large-scale simulations of cardiac electromechanics. On the other hand, SE models are to be preferred if the modeler is interested in studying the activation profile alongside the AFs (e.g. to investigate its role in the onset of the LDA phenomenon).
Definition of transition rates
We have shown that, under some physically motivated assumptions, the CTMC presented above can be described by different systems of ODEs and/or PDEs. The models listed in Table 2 are thus valid independently of the specific choice of the transition rates (with the only exception of the models SE-ODE and MF-ODE that require that the sum of the detachment and attachment rate is independent of the XB distorsion). In this section, we present and motivate the specific choice of transition rates that we will adopt in the rest of this paper.
In the spatially-explicit models introduced above, the transition rates possibly depend on the location of the RUs in the myofilaments. This allows to account for the myofilaments overlap (single overlap, duble overlap, no overlap). To facilitate the identification of the overlap region corresponding to a given RU from its index , we introduce the following functions, which define a smooth transition between the regions (see Fig 9):
where we have defined:
as the coordinates (with respect to the end of the AF closer to the center of the sarcomere) of the left and right ends of the MF (yLM, yRM), of the beginning of the single-overlap region (ySF) and of the i-th RU (yi). Hence, we have χM(SL, i) ≃ 1 if the i-th RU faces the considered half MF and χSF(SL, i) ≃ 1 if the i-th RU is in the single filament region (no overlap with other AFs occurs).
RUs transition rates.
The RU dynamics is determined by the eight rates associated with the forward and backward transitions ,
,
and
. The transition rates are affected by [Ca2+]i (that enhances in a multiplicative way the transition
), the filament overlap and the state of the nearest-neighboring Tm units (for the latter interaction we adopt the cooperative interactions proposed in [17]). We start by considering the single-overlap zone, where we adopt the transition rates of the model of [17]. We show below that the transition rates of [17] are, however, rather general, as they are based on just a couple of assumptions. We keep the notation consistent with [17] to allow for comparisons.
We call and, without loss of generality, we set
, where the constant μ allows to differentiate the two rates. Experiments carried out with protein isoforms from different species highlight that there is no apparent variation in the transition
in different combinations of Tn subunits and Tm [18]. We assume thus that the transition rates for
and for
coincide, and we set
. Conversely, we allow the reverse transition rates to depend on the state of the associated Tm. Concerning the transitions involving Tm, we assume that the calcium binding state of Tn affects the transition rate of
for the associated Tm, but not the reverse rate. Therefore, we set
, where n(α, δ) ∈ {0, 1, 2} denotes the number of permissive states among α and δ, as proposed in [17]. The physical interpretation of the constant γ corresponds to
, where kB is the Boltzmann constant, T the absolute temperature and ΔE denotes the energetic gain of the configuration of neighboring units in the same state (i.e.
or
) with respect to that with different states (i.e.
or
). See [66] for more details in this regard.
Then, without loss of generality we denote , where the constant Q allows to differentiate the forward and backward transition rates. The only transition rate left is given, to satisfy the detailed-balance consistency with the other rates, by
.
We remark that, due to the definition of n(α, δ) as the number of neighboring units in permissive state, the units located at the ends of the filament cannot have n = 2. Coherently with this, in Eq (4), the state of the missing neighboring RUs is set to .
In conclusion, the transition rates are determined by the five parameters Q, μ, kd, koff, kbasic (plus the parameter γ that regulates the amount of cooperativity), resulting from the eight free parameters constrained by the two assumptions ( not affected by Tm,
not affected by Tn) and by the detailed-balance consistency.
Concerning the dependence on the filament overlap, we assume that the only transition affected by filament overlap is , that is prevented in the central zone of the sarcomere, where the two AFs meet [82]. Specifically, we set, for
and for
:
(25)
The resulting 4-states CTMC associated with each RU is represented in Fig 10.
The terms depending on the intracellular calcium concentration [Ca2+]i are highlighted in red; terms depending on the state of neighbouring RUs (i.e. depending on n) are highlighted in green; terms depending on the position of the RU and the current sarcomere elongation are highlighted in blue.
XBs transition rates.
On the basis of the results of [81], we work under the hypothesis that the total attachment-detachment rate is independent of the XB elongation. In this case, the models SE-ODE and MF-ODE can be used in place of the more computationally expansive counterparts SE-PDE and MF-PDE, which involve the solution of a PDE system. As a matter of fact, in [81] we have shown that the introduction of such hypothesis significantly reduces the number of parameters to be fitted by experiments, still preserving the capability of the models of reproducing a wide range of experimental characterizations. Moreover, we also make the reasonable assumption that the sliding velocity only affects the detachment rate.
As already mentioned, the transition rates may be affected by the mutual arrangement of the filaments. Specifically, we assume that binding is possible only in the single-overlap region and when Tm is in the state . In other terms, we set:
(26)
We assume that the unbinding rate can take different values inside and outside the single-overlap region. Hence, we set, for and
:
(27)
Moreover, it is well motivated that out of the single-overlap zone, the detachment rates are not affected by the state of Tm (i.e. ) and that the detachment rate when Tm is in state
is not affected by the filament overlap (i.e.
). In summary, we have:
for some constant r0 and function q, such that q(0) = 0. Moreover, since we focus on the 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.
We notice that, from (10) and (26), it follows, for p = 0, 1:
(28)
Therefore, the parameters to calibrate are the scalars ,
, r0, α and aXB, to link the microscopic force with the macroscopic active tension.
Parameters calibration
We developed a pipeline to calibrate the parameters of the models proposed in this paper from measurements typically available from experiments. Our strategy is based on the observation that different experimental setups involve different time scales and different aspects of the force generation phenomenon. This allows to decouple the role of parameters associated with different phenomena and with different time scales, thus calibrating them in sequential manner. The calibration pipeline is made of three building blocks: the calibration of the XB rates, the calibration of the RU rates ruling the steady-state solution and, finally, the calibration of the RU rates ruling the kinetics of force development and relaxation. We report below the main steps of the calibration pipeline. A detailed description of the different steps is available in the Supporting Information (S2 Appendix).
Calibration of the XBs rates.
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 will be clarified later. In [81] we have shown that the parameters of the distribution-moments equations describing the XB dynamics can be calibrated to fit the steady state force, the shape of the force-velocity relationship (see Fig 2b) and the slope of the tension-elongation curve following a fast step (see Fig 2c). We remark that the experimental setups associated with these measurements are are such that the thin filament activation machinery can be considered in steady-state. This observation is crucial since it allows to decouple the calibration of the parameters involved in the thin filament regulation from the calibration of the parameters involved in XBs cycling.
In conclusion, once the parameters of the thin filament activation model has been calibrated, we have at our disposal an automatic procedure to calibrate the remaining parameters. For this reason, we first setup such calibration procedure for the parameters associated with XB cycling (i.e. (11) or (21)) and, successively, we calibrate the parameters associated with RU activation (i.e. (5) or (15)), so that we can directly see the effect of changes of such parameters on the resulting force (the remaining parameters are automatically adjusted).
Calibration of the RUs rates (steady-state).
The steady-state solution of the thin filament activation models (i.e. (5) and (15)) only depends on the ratio between the pairs of opposite transition rates (e.g. the ratio ). Therefore, the six parameters can be split into two groups: the first group (Q, μ, kd and γ) determines the steady-state solution, while the second group (koff, kbasic) only affects the kinetics of the model (that is to say how fast the transients are). This allows to calibrate first the parameters of the first group, and, only successively, those of the second group.
The fingerprint of the steady-state solution of the RU model is the force-calcium relationship (see Fig 2a). Hence, we tune the parameters Q, μ, kd and γ to fit this curve. Specifically, kd mainly acts on the calcium sensitivity (i.e. EC50), γ on the apparent cooperativity (i.e. nH), while Q and μ affect cooperativity, calcium sensitivity, the asymmetry of the force-calcium relationship below and above EC50 [17] and the SL-driven regulation on calcium sensitivity (in the SE-ODE model).
As we show in the section Results, the force-calcium curves obtained with the SE-ODE model exhibit the SL-induced change in calcium sensitivity observed in experiments (LDA, see Introduction). We remark that we are here able to reproduce the LDA without any phenomenological SL-dependent tuning of the parameter, as done, e.g. in [18, 24, 25, 33]. The LDA emerges from the SE-ODE model in a spontaneous way, as a consequence of the spatial-dependent interaction between the RUs (see [66] for a discussion on this topic). Conversely, with the MF-ODE model it is not possible to reproduce the LDA by simply acting on the model parameters. Indeed, the only effect of SL in the model is 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 mechanism reproducing LDA in the SE-ODE model is indeed intrinsically linked to the explicit spatial representation of the myofilaments [66]. Therefore, in the mean-field model MF-ODE, without an explicit spatial representation, we phenomenologically tune the calcium sensitivity kd in function of SL, by setting
(29)
where
.
Calibration of the RUs rates (kinetics).
To complete the calibration of the SE-ODE and MF-ODE models, we only need to set the parameters kbasic and koff, ruling the rapidity at which the transitions and
take place, respectively. Despite the fact that, at this stage, we need to calibrate just two parameters, this reveals some difficulties, mainly related to the following two aspects. First, the interplay between the two parameters is tight and their roles cannot be easily decoupled [18, 83, 84]. This results in a poor identifiability of the parameters: different combinations of parameters give similar results in terms of force transients. This issue has been reported also by [85], while calibrating the models of [23] and [18]. Additionally, the force transients predicted by the model are very sensitive to the calcium transient at input (this is a typical feature of activation models, see e.g. [85]). Therefore, since the experimentally measured calcium transients are much affected by noise (see e.g. [43, 85, 86], a calibration based on the best fit of the model response to experimentally measured calcium transients should be performed with care.
Based on the former remarks, we calibrate the parameters kbasic and koff by the following procedure. We consider force transients experimentally measured during isometric twitches and synthetic calcium transients fitted from experimentally measured ones. Then, we perform a MC sampling of the parameters kbasic and koff within prescribed ranges, and we select those values for which the force transients predicted by the model best fit the experimental ones.
Calibration from rat and human experimental data.
Due to the lack of a sufficiently large set of measurements from human cells at body temperature [26] to adequately fit all the model parameters, to calibrate the SE-ODE and MF-ODE models for body-temperature human cardiomyocytes we proceed as follows. First, we calibrate the model parameters from rat experiments at room temperature (for which available data are more abundant) and then we adjust the parameters that are reasonably affected by the two varying factors (i.e. inter-species variability and temperature) to fit the available data from human cells as body temperature. We compensate in this way for the data deficiency. We notice that we work under the hypothesis that inter-species variability does not affect the fundamental processes of tissue activation and force generation, but, since different species express different isoforms of the same protein, it can be encompassed by changing the parameters of the same mathematical model (see [85] for a detailed discussion on this topic).
Specifically, different species mainly differ in their calcium-sensitivity (i.e. kd) and in the kinetics (different species feature highly different heart rates), while temperature mainly affects the kinetics [12, 87, 88]. By exploiting the decoupling of the parameters of the RUs model ruling the steady-state relationships from those ruling the kinetics, we first focus on the steady-state, and we adjust kd to reflect the higher calcium sensitivity of human cells compared to rat [18, 26, 85]. Then, we re-calibrate the parameters koff and kbasic based on the kinetics of human force transients experimentally measured from human cells at body temperature.
We provide in Table 3 the full list of parameters for both species (room-temperature rat and body-temperature human) and for both models (SE-ODE and MF-ODE). Finally, we report in Table 4 the geometrical constants describing the size of the myofilament components that we use in the following.
Results and discussion
We show the results of numerical simulations obtained with the SE-ODE and MF-ODE models, performed under different experimental settings. The numerical schemes employed to approximate the solutions of the models here proposed are described in the Supporting Information (S3 Appendix). The codes implementing the models proposed in this paper and used to produce the results presented in this section are freely available in the following online repository: https://github.com/FrancescoRegazzoni/cardiac-activation.
With this implementation, the computational cost of the SE-ODE model is of nearly 12 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 of physical time on the same computer platform. Moreover, all the data used to generate the figures of this paper are available in the static repository [89]: https://doi.org/10.5281/zenodo.3992553.
Steady-state results
First, we consider steady-state solutions. To numerically obtain the steady-state curves, we fix a level of [Ca2+]i and SL and we let the model reach the equilibrium solution.
Force-calcium relationship.
We report in Fig 11 the force-calcium curves obtained with the SE-ODE and MF-ODE models calibrated from room-temperature rat data, compared with the experimental data used for the calibration. We are able to well fit the main features of the curves, including the characteristic S-shape, the plateau forces at both SL, the significant cooperativity typical of the cardiac tissue and the SL-induced change in calcium sensitivity. In Fig 12 we report the steady-state curves obtained with the sets of parameters calibrated for human body temperature cells. Also in this case, the curves reproduce the main experimentally observed features reported in the Introduction.
Fig 13 shows the dependence of the Hill coefficient nH and of the half-activating calcium concentration EC50 on the sarcomere length SL. We notice that, while the MF-ODE model produces an Hill coefficient that is independent of SL (the reason is that the role of SL on activation is just that of shifting and rescaling the curves, thus leaving nH unaffected), the SE-ODE model predicts a small increase of nH with SL. Both the results are equally acceptable since there is no common agreement on whether the Hill coefficient depend on SL or not [40–42, 90, 91]. Both the models correctly predict for both species an increase of EC50 as SL decreases. The relationship is approximately linear in the typical working range of SL (as experimentally observed, e.g., in [41]), while, for small values of SL, the SE-ODE model produces a faster decrease of sensitivity.
Force-length relationship.
Fig 14 shows the ascending limb of the steady-state force-length relationship. For both the SE-ODE and MF-ODE models we observe a change of slope for saturating calcium concentration around 1.65 μM, coherently with the experimental observations [42, 90, 92, 93]. Moreover, both models predict the observed change of convexity of the force-length curves at different calcium levels [40, 42, 90].
Isometric twitches
The predicted isometric twitches obtained with the room-temperature rat models are compared with the experimental data used for their calibration in Fig 15. The calcium transient here employed is obtained by fitting the experimental transient of [86] with the synthetic curve reported in the Supporting Information (S2 Appendix).
Similarly, we show in Fig 16 the tension transients obtained by simulating isometric twitches giving as input to the human models the calcium transient of the ToR-ORd model [39]. We notice that both models predict the tension-dependent prolongation of the twitch time [43, 86, 94], as it can be seen from the normalized traces reported in the bottom lines of the figures. We remark that, despite recent measurements on rabbit cells show that part of the increase of twitch force is linked to a larger calcium release under stretch [59], in this paper—due to the lack of experimental data showing a similar effect in human cells—we employ for simplicity the same calcium transient for all lengths.
In order to quantitatively assess the effects of SL on twitches, we report in Fig 17 the dependence on SL of the peak force () and of some synthetic indicators of the kinetics. Specifically, we consider the time-to-peak TTP (defined as the time separating the beginning of the stimulus and the tension peak) and the relaxation times RT50 and RT90 (defined as the time needed to accomplish 50% and 90% of relaxation, respectively). We notice that both models feature a kinetics with characteristic times matching those obtained experimentally. Moreover, both the peak tension and the characteristic times feature the expected increasing behavior with respect to SL. Remarkably, in the SE-ODE model, this feature spontaneously emerges from the model. Conversely, in the MF-ODE model, the prolongation of TTP and of the relaxation times is to be ascribed to the change of calcium sensitivity with respect to SL. Furthermore, it is remarkable that, although in both the cases the kinetic constants koff and kbasic have been calibrated to fit the twitch kinetics for a given SL, the simulations yields a good experimental match also for other values of SL (this happens despite the SL-dependent effect on calcium sensitivity is calibrated under a different experimental setting, i.e. steady-state conditions).
When possible, model results are compared with experimental data from [95].
Force-velocity relationship
Fig 18 shows the force-velocity relationship predicted with the rat models, compared with the experimental data used for the calibration. The human model yields similar results. For both the SE-ODE and the MF-ODE model, the numerical results correctly reproduce the experimentally observed convex profile, with the force reaching zero in correspondence of a finite value of velocity, the so-called maximum shortening velocity (see Introduction). Moreover, the value of is not significantly affected by the level of activation (that is to say, by the values of [Ca2+]i and SL), as well as the curvature of the curve. This is also coherent with the experimental observations [70]. The good agreement with the experimental data serves as a validation for the automatic calibration procedure presented above.
Fast force transients
We consider the response to fast steps predicted by the SE-ODE and the MF-ODE models. With this aim, we set a fixed value for the calcium concentration and sarcomere length (we set [Ca2+]i = 1.2 μM and SL = 2.2 μm, but the results are not significantly 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 small time interval Δt = 200 μs, and we plot the tension at the end of the step as a function of the step length ΔL.
We show in Fig 19 the results obtained with the rat models, with a comparison with experimental data (similar results are obtained with the human models). The good match between the simulation results and the experimental measurements provide a further validation of the calibration procedure.
Solid line: model result; circles: T2-L2 experimental data from [70].
Multiscale cardiac electromechanics
Finally, we consider a multiscale cardiac EM model of the left ventricle (LV), that we describe in the following. Further details on the EM model and its numerical discretization can be found in [66]. For the sake of brevity, we only show the results obtained by considering the MF-ODE model for human body-temperature cardiomyocytes.
We employ a computational domain Ω0 derived from the Zygote heart model, representing the 50th percentile of a 21 years old healthy caucasian male [96], in which we generate the fibers and sheets distribution by the rule-based algorithm proposed in [97]. We model the electrophysiological activity of the heart tissue by means of the monodomain equation [98, 99], coupled with the TTP06 model for the ionic activity [38]. We introduce a deformation map and we define the displacement vector as d(X, t) = φ(X, t) − X. The mechanical behavior of the myocardium is described by means of the balance of momentum equation for the displacement d [100, 101], where we model the passive properties of the cardiac tissue through the quasi-incompressible exponential constitutive law of [102]. To account for the presence of the pericardium, we employ the generalized Robin boundary conditions of [7, 103] on the epicarial boundary. On the LV base we adopt the energy-consistent boundary condition that we proposed in [104], accounting for the effect of the neglected part of the domain on the artificial boundary of the LV base.
Within a multiscale setting, we describe the force generation mechanisms at the microscale by means of the MF-ODE model, proposed in this paper. The intracellular calcium concentration ([Ca2+]i) is provided in each point of the computational domain by the TTP06 model. The local sarcomere length, in turn, is assumed to be proportional to the tissue stretch in the direction of muscel fibers f0, that is we set , where SL0 denotes the sarcomere length at rest,
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
.
The output of the MF-ODE model, namely the active tension magnitude field Ta, provides the link from the microscopic to the macroscopic level, where the effect of the microscopically generated active force is given by the following active stress tensor [66, 105]:
Finally, the 3D LV model is coupled to a 0D model of blood circulation, consisting of a two-elements Windkessel model [106].
For the spatial discretization of the equations involved in the EM model, we employ piece-wise linear Finite Elements of a tetrahedral computational mesh with 354 ⋅ 1010 cells. For the discretization of time derivatives, we use first-order finite difference schemes [107] with the implicit-explicit (IMEX) scheme described in [66]. The coupling of the different core models is obtained by means of the segregated strategy presented in [8].
The results of the numerical simulation are provided in Fig 20, where we show the deformation of the LV at different time steps of an heartbeat, and in Fig 21, showing the time evolution of the variables involved in the force generation phenomenon and of the LV pressure and volume predicted by the multiscale EM model. In Fig 22 we show the pressure-volume loops obtained for different preloads, i.e. by varying the end-diastolic pressure pED. Our model correctly predicts the increased stroke volume following a raise in the preload. This phenomenon, known as Frank-Starling effect, represents a self-regulatory mechanism of fundamental importance, as it allows the cardiac output to be synchronized with the venous return [13]. The microscopic driver of 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 force generation.
Top row: full geometry. Middle row: half domain (top view). Bottom row: half domain (frontal view).
Time evolution of [Ca2+]i, SL, v and Ta (minimum, maximum and average over the computational domain) and of left ventricle pressure and volume.
Conclusions
In this paper we have derived four different models for the generation of active force in the cardiac muscle tissue. Such models are based on a biophysically detailed description of the microscopic mechanisms of regulation and activation of the contractile proteins. Still, their numerical realization features a contained computational cost. Indeed, their numerical resolution does not require the computationally expensive MC method.
The main difficulties to be addressed in the derivation of full-sarcomere models concern (i) the spatial correlation of the states of the RUs due to the nearest-neighbor interactions of Tm, which hinders a straightforward decoupling of the adjacent RUs, and (ii) the mutual filament sliding, that changes which RU regulates which XB from time to time. Similarly to what we did in [37], we address the first problem by introducing a conditional independence assumption for far RUs, given the state of intermediate RUs (this is coherent with the local nature of nearest-neighboring interactions). This assumption dramatically reduces the number of equations needed to describe the evolution of the stochastic processes: while in the original model of [17] the number of variables of the CTMC increases exponentially with the number of RUs, 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 latter can be considered independently of that of the former. Moreover, we depart from the traditional MHs-centered representation of XBs, in favour of a BSs-centered point of view. Thanks to this change of perspective, we derived a set of equations describing the XB dynamics without the need to track the mutual position of the RUs and the MHs.
Under the hypothesis that the total attachment-detachment rate is independent of the myosin arm stretch (as done in [22, 78]), the PDE system describing the XBs can be replaced by a system of ODEs. We remark that this is not a simplificatory assumption, like the conditional independence assumptions mentioned before, but rather a feature of the specific modeling choice for the transition rates describing the attachment-detachment process.
We have also presented a class of models (MF-PDE and MF-ODE), such that the 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-field models differ from the mean-field models of [16, 30, 31, 61, 79]. The latter, indeed, consider a single RU, instead of a triplet. In this manner, the short-range spatial correlation, 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 remarkably good agreement between model predictions as experimental measurements, in particular in the reproduction of the highly cooperative steady-state force-calcium curves. We have then calibrated the SE-ODE and the MF-ODE models for room-temperature rat intact cardiomyocytes and, later, body-temperature human intact cardiomyocytes.
The SE-ODE model predicts the so-called LDA (the increment of calcium sensitivity when the sarcomere length increases), phenomenon whose microscopic source is still debated [44, 50, 53, 56–58, 60]. Interestingly, in our SE-ODE model, the LDA spontaneously emerges without any phenomenological tuning the parameters in dependence of SL, thanks to the explicit spatial representation of the spatially-dependent nearest-neighborhood interactions among RUs. In particular, we believe that this is linked to the RUs located at the end-points of the single-overlap zone, characterized by a low probability of being in the state (because they have at most one neighbor in state
). Since nearest-neighbor interactions propagate the low probability of
towards the center of the filament, small values of SL are penalized with respect to large values of SL, the role of such border effect being more pronounced (further details on this topic can be found in [66]). Therefore, if the hypothesis that the RUs located at the end-point of the single-overlap zone behave as if the outer neighboring units are in state
is accepted, then the SE-ODE model provides a possible explanation for LDA. Conversely, if this hypothesis in not accepted, then the effect of the edges can be neglected and the mean-field approximation underlying the MF-ODE model is well motivated. In conclusion, according to the hypothesis made on the behavior of 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 based on a microscopically detailed representation of the regulatory and contractile proteins, where phenomenological modeling choices are only introduced for phenomena whose underlying mechanisms are not fully understood [53, 57, 58, 60]. The models proposed in this paper could be enriched to investigate alternative hypotheses for the LDA (such as the force-dependent recruitment of MHs from an “off” state, in which the interaction with the BSs is prevented [59]).
The results of the numerical simulations showed that our models can reproduce the main features of the experimental characterizations of muscle contraction associated with the time scales of interest of cardiac EM (that is to say, the time scales larger than that of the power-stroke), including the steady-state relationship between calcium and force and between sarcomere length and force, the kinetics of activation of relaxation, the length-dependent twitches prolongation and peak force increase, and the force-velocity relationship. Additionally, we plan to study the response of the models to changes in the heart rate, which has been reported to have an impact on calcium sensitivity, developed force and twitch kinetics [108–110]. However, this investigation requires a model for calcium dynamics capable of accounting for the effects of the heart rate on the calcium transients [108].
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-ODE model, able of reproducing realistic pressure-volume loops. Moreover, our multiscale EM model is capable of reproducing—at the macroscale—the Frank-Starling self-regulation mechanism, in virtue of the length-dependent effects characterizing—at the microscale—the force-generation model. This macroscopic effect, emerging from a microscopic phenomenon, can be regarded as a proof of concept for potential uses of the models proposed in this paper in the context of 3D numerical simulations. The latter, indeed, can help to investigate the links between microscopic mechanisms and organ-level phenomena and to elucidate the relationships intercurrent between the microscopic variables and the macroscopic ones. Furthermore, the employment of a biophysically detailed model in organ-level simulations might reveal the links between cardiomyopathies and their cellular or molecular basis.
Supporting information
S1 Appendix. Models derivation.
This appendix contains the mathematical details about the formal derivation of the models proposed in this paper.
https://doi.org/10.1371/journal.pcbi.1008294.s001
(PDF)
S2 Appendix. Parameters calibration.
This appendix provides further details about the calibration pipeline of the proposed models.
https://doi.org/10.1371/journal.pcbi.1008294.s002
(PDF)
S3 Appendix. Numerical schemes.
This appendix provides details about the numerical schemes employed to approximate the solution of the models proposed in this paper.
https://doi.org/10.1371/journal.pcbi.1008294.s003
(PDF)
References
- 1. Fishman GI, Chugh SS, DiMarco JP, Albert CM, Anderson ME, Bonow RO, et al. Sudden cardiac death prediction and prevention: report from a National Heart, Lung, and Blood Institute and Heart Rhythm Society Workshop. Circulation. 2010;122(22):2335–2348. pmid:21147730
- 2. Smith N, Nickerson D, Crampin E, Hunter P. Multiscale computational modelling of the heart. Acta Numerica. 2004;13:371–431.
- 3. Crampin EJ, Halstead M, Hunter P, Nielsen P, Noble D, Smith N, et al. Computational physiology and the physiome project. Experimental Physiology. 2004;89(1):1–26. pmid:15109205
- 4. Fink M, Niederer SA, Cherry EM, Fenton FH, Koivumäki JT, Seemann G, et al. Cardiac cell modelling: observations from the heart of the cardiac physiome project. Progress in Biophysics and Molecular Biology. 2011;104(1):2–21. pmid:20303361
- 5. Chabiniok R, Wang VY, Hadjicharalambous M, Asner L, Lee J, Sermesant M, et al. Multiphysics and multiscale modelling, data–model fusion and integration of organ physiology in the clinic: ventricular cardiac mechanics. Interface Focus. 2016;6(2):20150083. pmid:27051509
- 6. Quarteroni A, Lassila T, Rossi S, Ruiz-Baier R. Integrated Heart–Coupling multiscale and multiphysics models for the simulation of the cardiac function. Computer Methods in Applied Mechanics and Engineering. 2017;314:345–407.
- 7. Gerbi A, Dedè L, Quarteroni A. A monolithic algorithm for the simulation of cardiac electromechanics in the human left ventricle. Mathematics in Engineering. 2018;1(1):1–37.
- 8.
Dedè L, Gerbi A, Quarteroni A. Segregated algorithms for the numerical simulation of cardiac electromechanics in the left human ventricle. In: Ambrosi D, Ciarletta P, editors. The Mathematics of Mechanobiology. Springer; 2020. p. 81–116.
- 9.
Quarteroni A, Dede’ L, Manzoni A, Vergara C. Mathematical modelling of the human cardiovascular system: data, numerical approximation, clinical applications. Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press; 2019.
- 10. Salvador M, Dedè L, Quarteroni A. An intergrid transfer operator using radial basis functions with application to cardiac electromechanics. Computational Mechanics. 2020;66:491–511.
- 11. Azzolin L, Dedè L, Gerbi A, Quarteroni A. Effect of fibre orientation and bulk modulus on the electromechanical modelling of the human ventricles. Mathematics in Engineering. 2020;2(4):614–638.
- 12.
Bers D. Excitation-contraction coupling and cardiac contractile force. vol. 237. Springer Science & Business Media; 2001.
- 13.
Katz AM. Physiology of the heart. Lippincott Williams & Wilkins; 2010.
- 14. Huxley AF. Muscle structure and theories of contraction. Progress in Biophysics and Biophysical Chemistry. 1957;7:255–318. pmid:13485191
- 15. Huxley AF, Simmons RM. Proposed mechanism of force generation in striated muscle. Nature. 1971;233(5321):533–538. pmid:4939977
- 16. Rice JJ, Winslow RL, Hunter WC. Comparison of putative cooperative mechanisms in cardiac muscle: length dependence and dynamic responses. American Journal of Physiology-Heart and Circulatory Physiology. 1999;276(5):H1734–H1754.
- 17. Rice JJ, Stolovitzky G, Tu Y, de Tombe PP. Ising model of cardiac thin filament activation with nearest-neighbor cooperative interactions. Biophysical Journal. 2003;84(2):897–909. pmid:12547772
- 18. Niederer SA, Hunter PJ, Smith NP. A quantitative analysis of cardiac myocyte relaxation: a simulation study. Biophysical Journal. 2006;90(5):1697–1722. pmid:16339881
- 19. Hussan J, de Tombe PP, Rice JJ. A spatially detailed myofilament model as a basis for large-scale biological simulations. IBM Journal of Research and Development. 2006;50(6):583–600.
- 20. Marcucci L, Truskinovsky L. Muscle contraction: A mechanical perspective. The European Physical Journal E. 2010;32(4):411–418.
- 21.
Caruel M. Mechanics of Fast Force Recovery in striated muscles. PhD Thesis (Ecole Polytechnique); 2011.
- 22. Chapelle D, Le Tallec P, Moireau P, Sorine M. Energy-preserving muscle tissue model: formulation and compatible discretizations. International Journal for Multiscale Computational Engineering. 2012;10(2).
- 23. Land S, Niederer SA, Aronsen JM, Espe EK, Zhang L, Louch WE, et al. An analysis of deformation-dependent electromechanical coupling in the mouse heart. The Journal of Physiology. 2012;590(18):4553–4569. pmid:22615436
- 24. Washio T, Okada J, Takahashi A, Yoneda K, Kadooka Y, Sugiura S, et al. Multiscale heart simulation with cooperative stochastic cross-bridge dynamics and cellular structures. Multiscale Modeling & Simulation. 2013;11(4):965–999.
- 25. Washio T, Yoneda K, Okada J, Kariya T, Sugiura S, Hisada T. Ventricular fiber optimization utilizing the branching structure. International Journal for Numerical Methods in Biomedical Engineering. 2015;. pmid:26453026
- 26. Land S, Park-Holohan S, Smith NP, dos Remedios CG, Kentish JC, Niederer SA. A model of cardiac contraction based on novel measurements of tension development in human cardiomyocytes. Journal of Molecular and Cellular Cardiology. 2017;106:68–83. pmid:28392437
- 27. Riccobelli D, Ambrosi D. Activation of a muscle as a mapping of stress–strain curves. Extreme Mechanics Letters. 2019;28:37–42.
- 28. Caruel M, Truskinovsky L. Physics of muscle contraction. Reports on Progress in Physics. 2018;81(3):036602. pmid:28649969
- 29. Caruel M, Moireau P, Chapelle D. Stochastic modeling of chemical–mechanical coupling in striated muscles. Biomechanics and Modeling in Mechanobiology. 2019;18(3):563–587. pmid:30607642
- 30. Sachse FB, Glänzel KG, Seemann G. Modeling of protein interactions involved in cardiac tension development. International Journal of Bifurcation and Chaos. 2003;13(12):3561–3578.
- 31. Razumova MV, Bukatina AE, Campbell KB. Stiffness-distortion sarcomere model for muscle simulation. Journal of Applied Physiology. 1999;87(5):1861–1876. pmid:10562631
- 32. Campbell SG, Lionetti FV, Campbell KS, McCulloch AD. Coupling of adjacent tropomyosins enhances cross-bridge-mediated cooperative activation in a Markov model of the cardiac thin filament. Biophysical Journal. 2010;98(10):2254–2264. pmid:20483334
- 33. Washio T, Okada J, Sugiura S, Hisada T. Approximation for cooperative interactions of a spatially-detailed cardiac sarcomere model. Cellular and Molecular Bioengineering. 2012;5(1):113–126. pmid:22448201
- 34. Land S, Niederer SA. A spatially detailed model of isometric contraction based on competitive binding of troponin I explains cooperative interactions between tropomyosin and crossbridges. PLoS Computational Biology. 2015;11(8):e1004376. pmid:26262582
- 35. Hunter PJ, McCulloch AD, Ter Keurs HEDJ. Modelling the mechanical properties of cardiac muscle. Progress in Biophysics and Molecular Biology. 1998;69(2):289–331. pmid:9785944
- 36. Rice JJ, Wang F, Bers DM, de Tombe PP. Approximate model of cooperative activation and crossbridge cycling in cardiac muscle using ordinary differential equations. Biophysical Journal. 2008;95(5):2368–2390. pmid:18234826
- 37. Regazzoni F, Dedè L, Quarteroni A. Active contraction of cardiac cells: a reduced model for sarcomere dynamics with cooperative interactions. Biomechanics and Modeling in Mechanobiology. 2018;17:1663–1686. pmid:30003434
- 38. Ten Tusscher KH, Panfilov AV. Alternans and spiral breakup in a human ventricular tissue model. American Journal of Physiology-Heart and Circulatory Physiology. 2006;291(3):H1088–H1100. pmid:16565318
- 39. Tomek J, Bueno-Orovio A, Passini E, Zhou X, Minchole A, Britton O, et al. Development, calibration, and validation of a novel human ventricular myocyte model in health, disease, and drug block. Elife. 2019;8:e48890. pmid:31868580
- 40. Kentish JC, ter Keurs HEDJ, Ricciardi L, Bucx JJJ, Noble MIM. Comparison between the sarcomere length-force relations of intact and skinned trabeculae from rat right ventricle. Influence of calcium concentrations on these relations. Circulation Research. 1986;58(6):755–768. pmid:3719928
- 41. Dobesh DP, Konhilas JP, de Tombe PP. Cooperative activation in cardiac muscle: impact of sarcomere length. American Journal of Physiology-Heart and Circulatory Physiology. 2002;51(3):H1055.
- 42. Ter Keurs HEDJ, Shinozaki T, Zhang YM, Zhang ML, Wakayama Y, Sugai Y, et al. Sarcomere mechanics in uniform and non-uniform cardiac muscle: a link between pump function and arrhythmias. Progress in biophysics and molecular biology. 2008;97(2-3):312–331. pmid:18394686
- 43. Backx PH, Gao W, Azan-Backx MD, Marbán E. The relationship between contractile force and intracellular [Ca2+] in intact rat cardiac trabeculae. The Journal of General Physiology. 1995;105(1):1–19. pmid:7730787
- 44. Rice JJ, de Tombe PP. Approaches to modeling crossbridges and calcium-dependent activation in cardiac muscle. Progress in Biophysics and Molecular Biology. 2004;85(2):179–195. pmid:15142743
- 45. Dupuis LJ, Lumens J, Arts T, Delhaas T. Mechano-chemical Interactions in Cardiac Sarcomere Contraction: A Computational Modeling Study. PLoS Computational Biology. 2016;12(10):e1005126. pmid:27716775
- 46. Gordon A, Regnier M, Homsher E. Skeletal and cardiac muscle contractile activation: tropomyosin “rocks and rolls”. Physiology. 2001;16(2):49–55.
- 47. Fitzsimons DP, Patel JR, Moss RL. Cross-bridge interaction kinetics in rat myocardium are accelerated by strong binding of myosin to the thin filament. The Journal of Physiology. 2001;530(2):263–272. pmid:11208974
- 48. Craig R, Lehman W. Crossbridge and tropomyosin positions observed in native, interacting thick and thin filaments. Journal of Molecular Biology. 2001;311(5):1027–1036. pmid:11531337
- 49. Sun YB, Lou F, Irving M. Calcium-and myosin-dependent changes in troponin structure during activation of heart muscle. The Journal of Physiology. 2009;587(1):155–163. pmid:19015190
- 50. Farman GP, Allen EJ, Schoenfelt KQ, Backx PH, De Tombe PP. The role of thin filament cooperativity in cardiac length-dependent calcium activation. Biophysical Journal. 2010;99(9):2978–2986. pmid:21044595
- 51. Heeley D, Smillie L, Lohmeier-Vogel E. Effects of deletion of tropomyosin overlap on regulated actomyosin subfragment 1 ATPase. Biochemical Journal. 1989;258(3):831–836. pmid:2525026
- 52. Pan B, Gordon A, Luo Z. Removal of tropomyosin overlap modifies cooperative binding of myosin S-1 to reconstituted thin filaments of rabbit striated muscle. Journal of Biological Chemistry. 1989;264(15):8495–8498. pmid:2722785
- 53. Sequeira V, van der Velden J. The Frank–Starling Law: a jigsaw of titin proportions. Biophysical Reviews. 2017;9(3):259–267. pmid:28639137
- 54. Ter Keurs H, Rijnsburger WH, Van Heuningen R, Nagelsmit MJ. Tension development and sarcomere length in rat cardiac trabeculae. Evidence of length-dependent activation. Circulation Research. 1980;46(5):703–714. pmid:7363419
- 55. Fitzsimons DP, Moss RL. Strong binding of myosin modulates length-dependent Ca2+ activation of rat ventricular myocytes. Circulation Research. 1998;83(6):602–607. pmid:9742055
- 56. Pearson JT, Shirai M, Tsuchimochi H, Schwenke DO, Ishida T, Kangawa K, et al. Effects of sustained length-dependent activation on in situ cross-bridge dynamics in rat hearts. Biophysical Journal. 2007;93(12):4319–4329. pmid:17766361
- 57. Ait-Mou Y, Hsu K, Farman GP, Kumar M, Greaser ML, Irving TC, et al. Titin strain contributes to the Frank–Starling law of the heart by structural rearrangements of both thin-and thick-filament proteins. Proceedings of the National Academy of Sciences. 2016;113(8):2306–2311.
- 58. de Tombe PP, ter Keurs HE. Cardiac muscle mechanics: sarcomere length matters. Journal of Molecular and Cellular Cardiology. 2016;91:148–150. pmid:26678623
- 59. Campbell KS, Janssen PM, Campbell SG. Force-dependent recruitment from the myosin off state contributes to length-dependent activation. Biophysical journal. 2018;115(3):543–553. pmid:30054031
- 60. Niederer SA, Campbell KS, Campbell SG. A short history of the development of mathematical models of cardiac mechanics. Journal of molecular and cellular cardiology. 2019;127:11–19. pmid:30503754
- 61. Landesberg A, Sideman S. Coupling calcium binding to troponin C and cross-bridge cycling in skinned cardiac cells. American Journal of Physiology-Heart and Circulatory Physiology. 1994;266(3):H1260–H1271.
- 62. Zou G, Phillips GN Jr. A cellular automaton model for the regulatory behavior of muscle thin filaments. Biophysical Journal. 1994;67(1):11–28. pmid:7918978
- 63. Dobrunz LE, Backx PH, Yue DT. Steady-state [Ca2+] i-force relationship in intact twitching cardiac muscle: direct evidence for modulation by isoproterenol and EMD 53998. Biophysical Journal. 1995;69(1):189–201. pmid:7669896
- 64.
Norris JR. Markov Chains. 2. Cambridge University Press; 1998.
- 65.
Bailey NTJ. The elements of stochastic processes with applications to the natural sciences. vol. 25. John Wiley & Sons; 1990.
- 66.
Regazzoni F. Mathematical modeling and Machine Learning for the numerical simulation of cardiac electromechanics. PhD Thesis (Politecnico di Milano); 2020.
- 67. Lymn R, Taylor EW. Mechanism of adenosine triphosphate hydrolysis by actomyosin. Biochemistry. 1971;10(25):4617–4624. pmid:4258719
- 68.
Keener JP, Sneyd J. Mathematical physiology. vol. 1. Springer; 2009.
- 69. Marcucci L, Truskinovsky L. Mechanics of the power stroke in myosin II. Physical Review E. 2010;81(5):051915.
- 70. Caremani M, Pinzauti F, Reconditi M, Piazzesi G, Stienen GJ, Lombardi V, et al. Size and speed of the working stroke of cardiac myosin in situ. Proceedings of the National Academy of Sciences. 2016;113(13):3675–3680.
- 71. Hill AV. The heat of shortening and the dynamic constants of muscle. Proceedings of the Royal Society of London B: Biological Sciences. 1938;126(843):136–195.
- 72. Huxely AF. Muscle structure and theories of contraction. Prog Biophys Biophys Chem. 1957; p. 255–318.
- 73. Kimmig F, Chapelle D, Moireau P. Thermodynamic properties of muscle contraction models and associated discrete-time principles. Advanced Modeling and Simulation in Engineering Sciences. 2019;6(1):6.
- 74.
Kimmig F, Caruel M, Moireau P, Chapelle D. Activation-contraction coupling in a multiscale heart model. In: Proceedings of CMBE 2019 (volume 1); 2019. p. 96–99.
- 75. Smith D, Geeves MA, Sleep J, Mijailovich SM. Towards a unified theory of muscle contraction. I: foundations. Annals of Biomedical Engineering. 2008;36(10):1624–1640. pmid:18642081
- 76.
Kimmig F. Multi-scale modeling of muscle contraction—From stochastic dynamics of molecular motors to continuum mechanics. PhD Thesis (Université Paris-Saclay); 2019.
- 77. Zahalak GI. A distribution-moment approximation for kinetic theories of muscular contraction. Mathematical Biosciences. 1981;55(1-2):89–114.
- 78.
Bestel J, Clément F, Sorine M. A biomechanical model of muscle contraction. In: International conference on medical image computing and computer-assisted intervention. Springer; 2001. p. 1159–1161.
- 79.
Sachse FB. Computational cardiology: modeling of anatomy, electrophysiology, and mechanics (lecture notes in Computer Science). Secaucus, NJ, USA: Springer-Verlag New York, Inc.; 2004.
- 80. Sugiura S, Washio T, Hatano A, Okada J, Watanabe H, Hisada T. Multi-scale simulations of cardiac electrophysiology and mechanics using the University of Tokyo heart simulator. Progress in Biophysics and Molecular Biology. 2012;110(2):380–389. pmid:22828714
- 81. Regazzoni F, Dedè L, Quarteroni A. Active force generation in cardiac muscle cells: mathematical modeling and numerical simulation of the actin-myosin interaction. Vietnam Journal of Mathematics. 2020;.
- 82. Gordon A, Huxley AF, Julian F. The variation in isometric tension with sarcomere length in vertebrate muscle fibres. The Journal of Physiology. 1966;184(1):170–192. pmid:5921536
- 83. Poggesi C, Tesi C, Stehle R. Sarcomeric determinants of striated muscle relaxation kinetics. Pflügers Archiv. 2005;449(6):505–517. pmid:15750836
- 84. Piroddi N, Belus A, Scellini B, Tesi C, Giunti G, Cerbai E, et al. Tension generation and relaxation in single myofibrils from human atrial and ventricular myocardium. Pflügers Archiv-European Journal of Physiology. 2007;454(1):63–73. pmid:17123098
- 85. Tøndel K, Land S, Niederer SA, Smith NP. Quantifying inter-species differences in contractile function through biophysical modelling. The Journal of Physiology. 2015;593(5):1083–1111. pmid:25480801
- 86. Janssen PM, de Tombe PP. Uncontrolled sarcomere shortening increases intracellular Ca2+ transient in rat cardiac trabeculae. American Journal of Physiology—Heart and Circulatory Physiology. 1997;272(4):H1892–H1897.
- 87. Harrison SM, Bers DM. Influence of temperature on the calcium sensitivity of the myofilaments of skinned ventricular muscle from the rabbit. The Journal of General Physiology. 1989;93(3):411–428. pmid:2703821
- 88. Janssen PML, Stull LB, Marbán E. Myofilament properties comprise the rate-limiting step for cardiac relaxation at body temperature in the rat. American Journal of Physiology-Heart and Circulatory Physiology. 2002;282(2):H499–H507. pmid:11788397
- 89.
Regazzoni F, Dedè L, Quarteroni A. PLOS Comput. Biol. “Biophysically detailed mathematical models of multiscale cardiac active mechanics”: datasets; 2020. Available from: https://doi.org/10.5281/zenodo.3992553.
- 90. Ter Keurs HEDJ, Hollander EH, ter Keurs MHC. The effect of sarcomere length on the force–cytosolic [Ca2+] relationship in intact rat cardiac trabeculae. Skeletal muscle mechanics: from mechanics to function Wiley, New York. 2000; p. 53–70.
- 91. Jelicks LA, Gupta RK. Intracellular Free Magnesium and High Energy Phosphates in the Perfused Normotensive and Spontaneously Hypertensive Rat Heart A 31P NMR Study. American Journal of Hypertension. 1991;4(2_Pt_1):131–136. pmid:2021444
- 92.
Kentish J, Ter Keurs H, Allen D. The contribution of myofibrillar properties to the sarcomere length-force relationship of cardiac muscle. In: Starling’s law of the heart revisited. Springer; 1988. p. 1–17.
- 93.
Keurs H, Noble MI. Starling’s law of the heart revisited. vol. 89. Springer Science & Business Media; 2012.
- 94. Janssen PM, Hunter WC. Force, not sarcomere length, correlates with prolongation of isosarcometric contraction. American Journal of Physiology—Heart and Circulatory Physiology. 1995;269(2):H676–H685.
- 95. Chung JH, Martin BL, Canan BD, Elnakish MT, Milani-Nejad N, Saad NS, et al. Etiology-dependent impairment of relaxation kinetics in right ventricular end-stage failing human myocardium. Journal of molecular and cellular cardiology. 2018;121:81–93. pmid:29981798
- 96.
Zygote 3D models; 2019. Available from: https://www.zygote.com/.
- 97. Bayer JD, Blake RC, Plank G, Trayanova NA. A novel rule-based algorithm for assigning myocardial fiber orientation to computational heart models. Annals of Biomedical Engineering. 2012;40(10):2243–2254. pmid:22648575
- 98.
Colli Franzone P, Pavarino LF, Savaré G. Computational electrocardiology: mathematical and numerical modeling. In: Complex systems in biomedicine. Springer; 2006. p. 187–241.
- 99.
Colli Franzone P, Pavarino LF, Scacchi S. Mathematical cardiac electrophysiology. vol. 13. Springer; 2014.
- 100.
Ogden RW. Non-linear elastic deformations. Courier Corporation; 1997.
- 101.
Antman SS. Nonlinear problems of elasticity. Springer; 1995.
- 102. Usyk TP, LeGrice IJ, McCulloch AD. Computational model of three-dimensional cardiac electromechanics. Computing and Visualization in Science. 2002;4(4):249–257.
- 103. Pfaller MR, Hörmann JM, Weigl M, Nagler A, Chabiniok R, Bertoglio C, et al. The importance of the pericardium for cardiac biomechanics: From physiology to computational modeling. Biomechanics and Modeling in Mechanobiology. 2019;18(2):503–529. pmid:30535650
- 104. Regazzoni F, Dedè L, Quarteroni A. Machine learning of multiscale active force generation models for the efficient simulation of cardiac electromechanics. Computer Methods in Applied Mechanics and Engineering. 2020;370:113268.
- 105. Giantesio G, Musesti A, Riccobelli D. A comparison between active strain and active stress in transversely isotropic hyperelastic materials. Journal of Elasticity. 2019; p. 1–20.
- 106. Westerhof N, Lankhaar JW, Westerhof BE. The arterial Windkessel. Medical & Biological Engineering & Computing. 2009;47(2):131–141.
- 107.
Quarteroni A, Valli A. Numerical approximation of partial differential equations. vol. 23. Springer Science & Business Media; 2008.
- 108. Coppini R, Ferrantini C, Yao L, Fan P, Del Lungo M, Stillitano F, et al. Late sodium current inhibition reverses electromechanical dysfunction in human hypertrophic cardiomyopathy. Circulation. 2013;127(5):575–584. pmid:23271797
- 109. Chung JH, Biesiadecki BJ, Ziolo MT, Davis JP, Janssen PM. Myofilament calcium sensitivity: role in regulation of in vivo cardiac contraction and relaxation. Frontiers in physiology. 2016;7:562. pmid:28018228
- 110. Chung JH, Milani-Nejad N, Davis JP, Weisleder N, Whitson BA, Mohler PJ, et al. Impact of heart rate on cross-bridge cycling kinetics in failing and nonfailing human myocardium. American Journal of Physiology-Heart and Circulatory Physiology. 2019;317(3):H640–H647. pmid:31347914