2 Macroscale: contraction at the organ level
In this chapter, we present the framework for modeling the heart contraction at the macroscale. We will here focus on the continuum mechanics formulation to make clear where the active component of the model enters the description. The modeling of the internal variables that represent the underlying active processes will be presented in the next chapter.
The contributions summarized in this chapter are from the work of Caruel et al. “Dimensional reductions of a cardiac model for effective validation and calibration” (2014) and Kimmig and Caruel “Hierarchical modeling of force generation in cardiac muscle” (2020). The general approach used in these two references is quite representative of the literature.
2.1 Background: mechanical models of the heart
2.1.1 A composite material
From the macroscale point of view, the cardiac muscle tissue can be represented as a composite material consisting in a hyperelastic polymer matrix with embeded active stress fibers, see Figure 2.1. The fibers’ orientation vary continuously from +60º at the endocardium to -60º at the epicardium1. Locally parallel fibers are themselves grouped in fascicles of 5 to 10 individuals surrounded by sheetlets of collagen, whose orientation also vary from the endocardium to the pericardium.2
Macroscale models are formulated using the laws of continuum mechanics which assume the existence of an elementary representative microscopic volume. For the muscle tissue, the usual elementary representative volume is constituted of a hyper-viscoelastic material representing the extracellular matrix, in which a 1D contractile fiber—characterized by its local orientation \(\underline{\tau}_{}\)—is embeded, see Figure 2.2 (a).3 This representation assumes the possibility to define a preferred fiber direction at each material point.
3 In this presentation we neglect the sheetlets organization and refer to Tueni, Allain, and Genet “On the structural origin of the anisotropy in the myocardium: Multiscale modeling and analysis” (2023) for a presentation of a model that takes it into account.
2.1.2 Active strain decomposition
We start the model presentation by mentioning that our approach is based on the so-called “active stress” assumption. The alternative “active strain” approach is based on the multiplicative decomposition of the deformation gradient tensor \[ \underline{\underline{F}}= \underline{\underline{F}}_{{\,p}}\cdot\underline{\underline{F}}_{\,a}, \] where \(\underline{\underline{F}}_{\,a}\) and \(\underline{\underline{F}}_{{\,p}}\) represent the active and passive deformations, respectively. The theory is developped in the following references: Nardinocchi and Teresi, “On the Active Response of Soft Living Tissues” (2007); Nobile, Quarteroni, and Ruiz-Baier, “An active strain electromechanical model for cardiac tissue: Active strain in cardiac electromechanics” (2012); Göktepe, Menzel, and Kuhl, “The generalized Hill model: A kinematic approach towards active muscle contraction” (2014). As an example, Colorado-Cervantes et al. “Patient-specific modeling of left ventricle mechanics” (2022) proposed the following expression for the active strain \[ \underline{\underline{F}}_{\,a}= \lambda_{\|}\, \underline{\tau}\otimes\underline{\tau}+ \frac{1}{\sqrt{\smash[b]{\lambda_{\|}}}}(\underline{\underline{\mathbb{I}}} - \underline{\tau}\otimes\underline{\tau}), \] where \(\lambda_{\|}\) measures the contraction along the fibers. Models of this type then provide a phenomenological relation between the parameter \(\lambda_{\|}\) and an activation field linked for instance to the transmembrane electric potential4.
We chose to follow a different approach based on a formulation that naturally describes the internal microscopic kinematics of a fiber and that can be closely linked to the actual molecular force generation process. This approach can thus benefit from a large set of experimental data targeting these microscale processes and from which different model ingredients can be specifically calibrated.
2.1.3 Model of the contracting muscle tissue
The proposed macroscale muscle tissue model originates from the work of Chapelle et al. “An energy-preserving muscle tissue model: Formulation and compatible discretizations” (2012) It is based on the Hill-Maxwell rheology represented in Figure 2.2 (b). This model was used in Caruel et al., “Dimensional reductions of a cardiac model for effective validation and calibration” (2014), and has been refined in Kimmig, Chapelle, and Moireau, “Thermodynamic properties of muscle contraction models and associated discrete-time principles” (2019). We here summarize the latter formulation.
The tissue model is decomposed of two parallel branches: a 3D branch representing the extracellular matrix and a 1D branch representing the fiber. The Green-Lagrange deformation tensors of the extracellular matrix and the contractile fiber are denoted by \(\underline{\underline{e}}_{{\,p}}\) and \(\underline{\underline{e}}_{{\,a}}\), respectively. Since the passive and active branch are in parallel, the global Green-Lagrange tensor reads \(\underline{\underline{e}}= \underline{\underline{e}}_{{\,p}}= \underline{\underline{e}}_{{\,a}}\). The fiber is usually viewed as a one dimensional homogenous medium where the global deformation can be mapped to the deformation of single half-sarcomere. Denoting by \(\ell_{\text{fib}}\) the characteristic length of a half sarcomere and by \(\delta\ell_{\text{fib}}\) its variation, we define the total extension of the fiber by \[ e_{\text{fib}}= \frac{\delta\ell_{\text{fib}}}{\ell_{\text{fib}}} = (1 + 2 \underline{\tau}\cdot\underline{\underline{e}}\cdot\underline{\tau})^{\tfrac{1}{2}}-1. \]
A mechanical geometrical representation of a single half-sarcomere that have been validated by numerous in situ mechanical experiments consists in a series connexion of two elements:5 (i) a passive element representing the elastic deformation of the myofilaments and other non-contractile sarcomeric structure, like the Z-disks, the M-line and titin,6 and (ii) an active element representing the array of myosin molecular motors interacting with the surrounding actin filaments, see Figure 2.2 (b).
5 see Pertici, Caremani, and Reconditi, “A mechanical model of the half-sarcomere which includes the contribution of titin” (2019); Caruel and Truskinovsky, “Physics of muscle contraction” (2018) and references therein. The kinematic decomposition is based on the work of Ford, Huxley, and Simmons “The relation between stiffness and filament overlap in stimulated frog muscle fibres.” (1981).
6 for more details about the sarcomeric proteins and their role in the contraction mechanism, see Chapter 5.
Following this assumption, we can decompose the extension \(\delta\ell_{\text{fib}}\) into \[ \delta\ell_{\text{fib}}= \delta\ell_{c}+ \delta\ell_{s}, \] where \(\delta\ell_{c}\) and \(\delta\ell_{s}\) denote the extensions of the active and passive series elements, respectively. We further define the associated deformations \(e_{c}= \delta\ell_{c}/\ell_{\text{fib}}\) and \(e_{s}= \delta\ell_{s}/\ell_{\text{fib}}\), which verify \(e_{\text{fib}}= e_{c}+ e_{s}\).
We denote by \(T_{\text{fib}}\) the active force in the direction of the fiber per unit area of transverse cross-section considered in the reference configuration. The stress \(T_{\text{fib}}\) combines the force produced by the contractile machinery \(T_{c}\) and a viscous drag that is considered linear, so that \[ T_{\text{fib}}= T_{c}+ \nu\dot{e}_{c}= \kappa_{s}e_{s}. \tag{2.1}\] The first Piola-Kirchhoff active stress tensor \(T_{\text{fib}}\) can finally be mapped to the second Piola-Kirchhoff stress tensor \[ \underline{\underline{\Sigma}}_{{\,a}}= \frac{T_{\text{fib}}}{(1+2\underline{\tau}\cdot\underline{\underline{e}}\cdot\underline{\tau})^{\frac{1}{2}}}\underline{\tau}\otimes\underline{\tau}. \]
The extracellular matrix is endowed with a 3D visco-hyperelastic constitutive law derived from a hyperelastic potential \(\mathcal{W}_{p}\) and a viscous pseudo-potential \(\mathcal{W}_{v}\). The resulting passive second Piola-Kirchhoff stress tensor is \[ \underline{\underline{\Sigma}}_{{\,p}}= \frac{\partial\mathcal{W}_{p}}{\partial\underline{\underline{e}}} + \frac{\partial\mathcal{W}_{v}}{\partial\dot{\underline{\underline{e}}}}. \] Finally, the total stress tensor aggregates the stresses in the active and passive branches \[ \underline{\underline{\Sigma}}= \underline{\underline{\Sigma}}_{{\,p}}+ \underline{\underline{\Sigma}}_{{\,a}}= \frac{\partial\mathcal{W}_{p}}{\partial\underline{\underline{e}}} + \frac{\partial\mathcal{W}_{v}}{\partial\dot{\underline{\underline{e}}}} + \frac{T_{c}+ \nu\dot{e}_{c}}{(1+2\,\underline{\tau}\cdot\underline{\underline{e}}\cdot\underline{\tau})^{\frac{1}{2}}}\underline{\tau}\otimes\underline{\tau}. \tag{2.2}\]
The principle of vitual work can then be written for any admissible displacement field \(\underline{w}_{}\) on the reference configuration \(\Omega_{0}\) with boundary \(S_{0}\): \[ \int_{\Omega_{0}}\rho_{0}\,\ddot{\underline{y}}\cdot\underline{w}\,\mathrm{d}\Omega + \int_{\Omega_{0}}\underline{\underline{\Sigma}}:\mathrm{d}_{\underline{y}}\,\underline{\underline{e}}\cdot\underline{w}\,\mathrm{d}\Omega_{0} = \int_{S_{0}} \underline{t}\cdot\underline{w}\, \mathrm{d}S, \tag{2.3}\] where, \(\rho_{0}\) is the mass density of the tissue in the reference configuration, \(y_{}\) is the displacement field, \(\underline{t}_{}\) is a force field actin on the boundary \(S_{0}\) of the domain \(\Omega_{0}\) and \[ \mathrm{d}_{\underline{y}}\,\underline{\underline{e}}\cdot\underline{w}= \frac{1}{2}\left[ \underline{\underline{\nabla}}\,\underline{w}+ (\underline{\underline{\nabla}}\,\underline{w})^{\top} + (\underline{\underline{\nabla}}\,\underline{y})^{\top}\cdot\underline{\underline{\nabla}}\,\underline{w} + (\underline{\underline{\nabla}}\,\underline{w})^{\top}\cdot\underline{\underline{\nabla}}\,\underline{y} \right]. \] In cardiac simulation the external load results typically from the intraventricular blood pressure and contact forces at the pericardium. To reproduce a cardiac cycle, Equation 2.3 have to be supplemented with an external circulation model and opening/closing valve laws that define the relashionship between the internal pressure and the blood outflow.7
The coupling between the macroscopic behavior and the process of force generation by the molecular motors is contained in the active stress \(T_{c}\) appearing in Equation 2.2, which will depend on local internal variable characterizing the state of the force production machinery. The usual hypothesis is to consider that the macroscopic stress \(T_{c}\) can be obtained by a simple rescaling of the force produced by a single molecular motor \(\tau_{c}\) in a mean-field-like approach: \[ T_{c}= \rho_{\text{m}}\tau_{c}, \] where \(\rho_{\text{m}}\) is the number of myosin molecular motors in a layer of thickness \(\ell_{\text{fib}}\) per unit cross-sectional area.
We delay the presentation of the molecular motors models behind the active force \(\tau_{c}\) to Chapter 3, but note here that the current framework directly connect the macroscale continnum mechanics balance laws with the nanoscale force generation dynamics. In this sence, the presented model, like most of the existing muscle tissue models is not multiscale, indeed. This aspect is further discussed in Chapter 4 and Chapter 5.
2.2 Contribution: model reduction for validation and calibration
To illustrate the model summarized in Section 2.1.3, we present the results obtained by Caruel et al. “Dimensional reductions of a cardiac model for effective validation and calibration” (2014). In this work we introduced two reduced models for the purpose of fast calibration and simulation.
8 For details about the formulation inhomogenous 1D model and some illustrations, we refer to Caruel et al. “Dimensional reductions of a cardiac model for effective validation and calibration” (2014)
The first model is a unidimensional representation of a muscle fiber aiming at reproducing uniaxial experiments. The formulation assumes a displacement field in the direction of the fiber. With this assumption, the principle of virtual work (Equation 2.3) reduced to a one dimensional partial differential equation and eventually an ordinary differential equation (0D model) if a uniform deformation is assumed. This uniform formulation is usually used for simulating the experiment performed on fibers. We will present some results of this approach in Chapter 3.8.
We here present another reduced model, that was formulated for fast simulation of the left ventricule cardiac cycle, using a spherical representation of that ventricle. Originally formulated by Caruel et al., this reduction was rewritten by Manganotti et al. “Coupling reduced-order blood flow and cardiac models through energy-consistent strategies: Modeling and discretization” (2021) to take into account the kinematics introduced by Kimmig, Chapelle, and Moireau “Thermodynamic properties of muscle contraction models and associated discrete-time principles” (2019), see Section 2.1.3.
The cardiac ventricle is represented by a thick sphere with radius \(R\) and thickness \(d\), see Figure 2.3. In the reference configuration, \(R= R_{0}\) and \(d=d_{0}\). We denote by \(y\) the uniform radial displacement and \(V\) the volume of the cavity such that
\[ R(y) = R_{0} + y,\, d(y) = d_{0}\left(1+\frac{y}{R_{0}}\right)^{-2} \text{ and } V(y) = \frac{4}{3}\pi\left(R(y) - \frac{d(y)}{2}\right)^{3}, \]
The passive component of the stress (extracellular matrix) can be derived from the hyperelastic potential \[ \mathcal{W}_{p}(\underline{\underline{e}}) = C_{0}\exp\big[C_{1}(J_{1}(\underline{\underline{e}})-3)^{2}\big] + C_{2}\exp\big[C_{3}(J_{4}(\underline{\underline{e}})-1)^{2}\big] \tag{2.4}\] that depends on the reduced invariant of the Green-Lagrange strain tensor \(J_{1}\) and \(J_4\);9, on the viscous potential \(\mathcal{W}_{v}(\underline{\underline{e}}) = \eta\,\dot{\underline{\underline{e}}}:\dot{\underline{\underline{e}}}\) and on four material parameters: \(C_{1},\,C_{2},\,C_{3}\) and \(C_{4}\).
9 Denoting by \(\underline{\underline{C}}\) the rigth Cauchy-Green deformation tensor we have \(J_{1} = I_{1}I_{3}^{-\frac{1}{3}}\) and \(J_{4} = I_{4}I_{3}^{-\frac{1}{3}}\) where \(I_{1} = \mathrm{tr}(\underline{\underline{C}})\) and \(I_{3} = \mathrm{det}(\underline{\underline{C}})\), and \(J_{4} = \underline{\tau}\cdot \underline{\underline{C}}\cdot\underline{\tau}\)
Combined with the sphere kinematics, Equation 2.3 and Equation 2.1 become \[ \left\{ \begin{aligned} &\rho_{0}|\Omega_{0}|\ddot{y} + \frac{|\Omega_{0}|}{R_{0}}\kappa_{s}\left(\frac{y}{R_{0}-e_{c}}\right) + \frac{\partial\mathcal{W}_{p}}{\partial y}(y) + \Lambda_{v}(y,\dot{y}) = P_{v}\frac{\partial V}{\partial y},\\ & \kappa_{s}\left(\frac{y}{R_{0}} - e_{c}\right) - \nu\dot{e}_{c}= \tau_{c}, \end{aligned} \right. \] where \(|\Omega_{0|}\) is the volume of the cavity in the reference configuration, \(P_{v}\) is the ventricular pressure and \[ \begin{aligned} &\mathcal{W}_{p}(y) = |\Omega_{0}|\left[ C_{0}\,e^{ C_{1}\left[ 2\left(1+\frac{y}{R_{0}}\right)^{2} + \left(1+\frac{y}{R_{0}}\right)^{4} - 3 \right]^{2} } + C_{2}\,e^{ C_{3}\left[ \left(1+\frac{y}{R_{0}}\right)^{2}-1 \right]^{2} } \right]\\ & \Lambda(y,\dot{y}) = 2\eta |\Omega_{0}| \left(1+\frac{y}{R_{0}}\right)^{2} \left[ 1 + 2\left(1+\frac{y}{R_{0}}\right)^{-12} \right]\dot{y}. \end{aligned} \]
To compute the ventricular pressure \(P_{v}\), the 0D model is coupled to a Windkessel representation of the blood circulation, which closes the system of equations. In this framework, the circulation system is analog to an electric circuit containing diodes, representing the valves, resistances representing the head loss in arteries and capacitances representing the elasticity of the arteries.10 The application of the classical Kirchhoff laws leads the closure of the system of equation driving the dynamics of the cavity contractile cycle.
10 see Manganotti et al., “Coupling reduced-order blood flow and cardiac models through energy-consistent strategies: Modeling and discretization” (2021), for an example of a two-stage Winkessel model
The parameters of the model are:
- the material paramters \(C_{i}\), characterizing the hyperelastic potential (2.4),
- the activation of the force generation process, here represented by the active forece \(\tau_{c}\) (see Chapter 3 for more details),
- the atrium pressure, which is a given function of time if one considers only the left ventricle,
- the characteristics of the circulation model.
A caridac cycle simulation performed using the 0D model is illustrated in Figure 2.4. The model was calibrated using experimental data obtained from mechanical test on heart papillary muscles mechanical tests.11 In this model, the active stress \(\tau_{c}\) is computed using a reduced model of the actin-myosin interaction that will be discussed in Chapter 3.
2.3 Remaining challenges and future work
Simulation platforms that produce physiologically relevant macroscale simulations of the heart contraction are now available.12 Such simulations in patient specific geometries are costly, but strategies have been developed to accelerate that process, using surrogate models based on asymptotic methods or AI assisted methods.13 Now these models can be used also for predicting the effect of drugs and mutations in the context of developed cardiomyopathies and large consortium are being set to develop clinical applications.
12 Chabiniok et al., “Multiphysics and multiscale modelling, datamodel fusion and integration of organ physiology in the clinic” (2016); Regazzoni, Dedè, and Quarteroni, “Biophysically detailed mathematical models of multiscale cardiac active mechanics” (2020); Filipovic et al., “SILICOFCM platform, multiscale modeling of left ventricle from echocardiographic images and drug influence for cardiomyopathy disease” (2022); Sugiura et al., “UT-Heart: A Finite Element Model Designed for the Multiscale and Multiphysics Integration of our Knowledge on the Human Heart” (2022).
13 Kimmig, Moireau, and Chapelle, “Hierarchical modeling of length-dependent force generation in cardiac muscles and associated thermodynamically-consistent numerical schemes” (2021); Regazzoni et al., “A machine learning method for real-time numerical simulations of cardiac electromechanics” (2022); Milićević et al., “Huxley muscle model surrogates for high-speed multi-scale simulations of cardiac contraction” (2022).
The most recent adavances in cardiac research concerns the activation and regulation pathways. We will report in Chapter 5 that the cytoskeletal proteins linking the muscle contractile unit together in the tissue play a major role in these fundamental processes. These mechanisms thus operate at a scale lying between the tissue scale and the molecular motors scale. The problem is that the large majority of available heart models are formulated using similar hypotheses as the ones used in this chapter. Two hypotheses have to be reconsidered if one wants to incorporate a realistic description of the intrinsic contraction and regulation pathways into a heart simulation framework:
- The contractile force in the macroscale mechanical formulation is obtained by a direct rescaling of the molecular motor force,
- The muscle fibers are locally homogenous.
At the time this manuscript is being writted, there seems to be no mechanical model available that provides convincing representation of the intermediate scales, bridging the nanoscale and the macroscale. This issue is explained in more details in Chapter 5.