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.
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 epicardium (Streeter Jr. and Bassett 1966). 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 (LeGrice et al. 1995; Nielles-Vallespin et al. 2017).
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).1 This representation assumes the possibility to define a preferred fiber direction at each material point.
1 In this presentation we neglect the sheetlets organization and refer to Tueni, Allain, and Genet (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}}_{\raise 0.3em {\,p}}\cdot\underline{\underline{F}}_{\raise 0.3em {\,a}}, \]
where \(\underline{\underline{F}}_{\raise 0.3em {\,a}}\) and \(\underline{\underline{F}}_{\raise 0.3em {\,p}}\) represent the active and passive deformations, respectively. The theory is developped in the following references: Nardinocchi and Teresi (2007); Nobile, Quarteroni, and Ruiz-Baier (2012); Göktepe, Menzel, and Kuhl (2014). As an example, Colorado-Cervantes et al. (2022) proposed the following expression for the active strain \[ \underline{\underline{F}}_{\raise 0.3em {\,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 potential (Göktepe, Menzel, and Kuhl 2014).
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. (2012). It is based on the Hill-Maxwell rheology represented in Figure 2.2 (b). This model was used in (Caruel et al. 2014), and has been refined by Kimmig, Chapelle, and Moireau (2019). We here summarize this latter formulation.
The models decomposed into 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}}_{\raise 0.3em {\,p}}\) and \(\underline{\underline{e}}_{\raise 0.3em {\,a}}\), respectively. Since the passive and active branch are in parallel, the global Green-Lagrange tensor reads \(\underline{\underline{e}}= \underline{\underline{e}}_{\raise 0.3em {\,p}}= \underline{\underline{e}}_{\raise 0.3em {\,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:2 (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,3 and (ii) an active element representing the array of myosin molecular motors interacting with the surrounding actin filaments, see Figure 2.2 (b).
2 see (Pertici, Caremani, and Reconditi 2019), (Caruel and Truskinovsky 2018) and references therein. The kinematic decomposition is based on the work of Ford, Huxley, and Simmons (1981).
3 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}}_{\raise 0.3em {\,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}}_{\raise 0.3em {\,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}}_{\raise 0.3em {\,p}}+ \underline{\underline{\Sigma}}_{\raise 0.3em {\,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, see Caruel et al. (2014) for more details.
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 are 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 in (Caruel et al. 2014). In this work we introduced two reduced models for the purpose of fast calibration and simulation.
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. For details about the formulation inhomogenous 1D model and some illustrations, we refer to (Caruel et al. 2014).
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 in by Caruel et al. (2014), this reduction was rewritten by Manganotti et al. (2021) to take into account the kinematics introduced by Kimmig, Chapelle, and Moireau (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,\quad 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\) (Caruel et al. 2014); and the viscous potential \(\mathcal{W}_{v}(\underline{\underline{e}}) = \eta\,\dot{\underline{\underline{e}}}:\dot{\underline{\underline{e}}}\).
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, see (Manganotti et al. 2021) for an example of a two-stage Winkessel model. The application of the classical Kirchhoff laws leads the closure of the system of equation driving the dynamics of the cavity contractile cycle.
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, see Caruel et al. (2014). 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 (Chabiniok et al. 2016; Regazzoni, Dedè, and Quarteroni 2020; Filipovic et al. 2022; Sugiura et al. 2022). 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 Milićević et al. (2022). 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.
The most recent adavances in cardiac research concerns the activation and regulation pathways. We will report in Chapter 5 that the structural 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.