A monolithic algorithm for the simulation of cardiac electromechanics in the human left ventricle

1.   Introduction

The heart plays the crucial role of pumping blood into the circulatory system by providing all sort of vital substances to the cells. Moreover, cardiovascular related diseases represent the leading causes of death in the whole world ^[60]. While advancements in medical practice are continuously improving patients conditions and diseases outcomes, recent progresses in mathematical modeling and computational mechanics allow to perform realistic cardiovascular numerical simulations ^{[38,65,76,77,84,97,105,107]}, thus potentially providing medical doctors and clinicians with valuable diagnostic and predictive tools. Moreover, clinical data and the application of image segmentation techniques to Magnetic Resonance Imaging (MRI) and Computed Tomography (CT) scans feed, as inputs, the mathematical models, thus allowing numerical simulations in subject- and patient-specific frameworks ^[23,55]; in addition, uncertainty quantification techniques allow to estimate the models parameters and data, and to cope with their variability ^[30,59].

The mathematical modeling of the heart involves several challenges intrinsically related to the complexity of its function ^[17,20,33]. A satisfactory model must be able to describe a wide range of different processes, such as the evolution of the transmembrane potential in the myocardium, the deformation caused by the muscles contraction, and the dynamics of the blood in the heart chambers and through the valves ^{[31,71,76,77,96,98,99]}. All these processes feature different temporal and spatial scales, yielding the so-called tyranny of scales ^[8]. A satisfactory knowledge of the models for isolated processes - the electrophysiology, the active and the passive mechanics - which we refer to as "single core models", is nowadays quite established; however, further theoretical and numerical studies are still necessary to better understand their mutual interactions ^{[18,22,30,82]}. In this work, we use state-of-the-art models in passive myocardial tissue modeling – the Holzapfel-Ogden model ^[46] – together with the active strain formulation ^[2,3] in combination with a recently proposed model for the transmurally heterogeneous thickening of the myocardium ^[7]; the latter is for the first time used in this paper for the integrated electromechanics modeling of the left ventricle (LV). The link between electrophysiology and mechanics is finally established through a model describing the shortening of the myocardial fibers ^[88], which is in turn triggered by a change in the ionic concentrations in the cardiac cells.

From the numerical viewpoint, we discretize in space the continuous single core models by means of the Finite Element Method (FEM) with both linear and quadratic elements, while the time discretization is carried on using Backward Differentiation Formulas (BDFs) of orders 1 and 2 ^[78]. BDFs are used both within implicit and semi-implicit schemes, the latter consisting in equal-order extrapolation of the unknowns in the nonlinear terms^[16,37]; however, the semi-implicit scheme poses strong limitations on the choice of the timestep size with respect to the implicit scheme, even if the latter yields a nonlinear problem, which we approximate by means of the Newton method. Integrating the discrete single core models leads to the formulation of a monolithic algebraic coupled problem. Numerical coupling of cardiac electromechanics is however typically performed by means of segregated or staggered approaches, for which the electrophysiology and mechanics solvers are not applied simultaneously ^{[4,5,18,38,54,84,107]}. It is however known that partitioned schemes (segregated or staggered) do not generally guarantee unconditional stability ^[15], an issue not concerning instead the monolithic scheme.

Monolithic schemes for cardiac electromechanics of the heart were firstly proposed and analyzed in ^[24,25], where only the contraction phase was simulated and the blood pressure at the endocardium (the LV inner wall) was neglected. More recently, a monolithic scheme including a 0D model for the vascular system, heart valve, and atrial chamber model was proposed in ^[44]. However, other than being less explored, monolithic schemes have been used so far only for cardiac electromechanics within the framework of the active stress formulation and only for idealized or simplified LV geometries. Instead, in this paper we propose and successfully use a monolithic algorithm for the numerical simulation of cardiac electromechanics in the framework of the active strain formulation; in addition, we simulate full cardiac cycles (one heartbeat) for both idealized and patient-specific geometries. Essential components of our framework are: modeling FSI at the LV endocardium wall by means of $0D$ models for the pressure variable along the heartbeat (^[30,84,107]) and the prestress technique, which we apply to the patient-specific LVs to estimate the internal stresses of the myocardium at the initial time of the simulation ^[47,100].

We solve the large, sparse block linear system stemming from the monolithic strategy by means of the GMRES method and, in order to speed-up the convergence of the linear solver, we employ a novel (right) preconditioner ^[89]. Our preconditioning strategy extends the one proposed in ^[27] for FSI problems and is based on the factorization of a block Gauss-Seidel preconditioner which exposes the multiphysics nature of the integrated problem, thus allowing to adapt its action on each single factor of the preconditioner. This strategy can be easily extended to a wide range of integrated multiphysics problems: the factors are indeed independently preconditioned by using algebraic multigrid ^[11] or domain decomposition ^[14,104] preconditioners tailored by exploiting physics-specific information (such as the number of equations) for each block; this local information would be instead lost by using a global monolithic preconditioner.

We use our monolithic electromechanics solver to perform numerical simulations in the High Performance Computing framework of the whole cardiac cycle, both for idealized and patient-specific LV geometries, and we analyze the numerical results in terms of clinically relevant indicators; specifically, we produce the so-called pressure-volume loops in order to assess the LV function and its properties.

The paper is organized as follows: in Sec. 2 we introduce the mathematical models for the electrophysiology, the mechanics and the activation of the myocardium; we then integrate them thus obtaining a continuous integrated model. In Sec. 3 we carry on the space and time numerical approximations of the single core models; in Sec. 4 we introduce the preconditioner used to solve the monolithic linear system arising after the time and space discretizations. In Sec. 5 we report and discuss the numerical results obtained with the proposed methods, and we finally draw our conclusions in Sec. 6.

2.   Mathematical models

The myocardium is a complex tissue composed of cardiomyocytes, organized in fibers and laminar collagen sheets ^[94], characterising the LV orthotropic internal structure. Therefore, the material properties of the LV are strongly dependent on the direction of fibers and sheets. First, the electrical conductivity of the cardiomyocytes is much larger in the longitudinal direction than transversally ^[52,75]; hence, at the macroscopic level, the transmembrane electric potential travels faster along the fibers. Secondly, it has been experimentally observed in ^[46] that the response of internal stresses to an external load significantly varies when measured among different directions since the myocardium is stiffer along the fibers. Finally, the contraction of the LV is made possible thanks to an active force acting on the cardiomyocytes lined up along the fibers direction ^[86]. In order to mathematically define fibers and sheets, we identify a local frame of reference by defining the mutually orthonormal vector fields ${\bf f}_0$ (fibers), ${\bf s}_0$ (sheets) and ${\bf n}_0$ (aka normals), where ${\bf n}_0$ is such that ${\bf n}_0$, ${\bf f}_0$, and ${\bf n}_0$ are mutually orthogonal, as depicted in Figure 1.

2.1.   Electrophysiology

Electrophysiology models take into account the electrochemical reactions occuring in the myocardium ^[21] triggered by an electric impulse originated at the sino-atrial node and then conveyed through the Purkinje fibers ^[20]. Such signal causes a quick depolarization of the LV cardiomyocytes, meaning that the transmembrane potential (i.e. the difference in electric potential between the interior and the exterior of a cell) changes sign in few microseconds. The potential drives a change in the concentration of different ionic species flowing through the so-called ionic gates located on the cellular membrane; the concentration of these ionic species, in turn, influences the potential thus slowly repolarizing the cells. The continuous interaction between the ions concentration and the potential causes a cascade effect for which a fast traveling wave known as action potential propagates in the whole myocardium ^[57].

In this work, we consider the monodomain equation for the description of the evolution of the cellular transmembrane potential $V$, a nonlinear diffusion-reaction equation obtained by homogenization of the bidomain equations ^{[20,48,75,92]}. The latter is indeed a richer, but more complicated model than the monodomain equation which is required for pathological conditions. In physiological conditions, the monodomain model is adequate and reads:

Here, $\Omega_0$ is the reference computational domain, a smooth and bounded subset of $\mathbb R^3$ represented e.g. by the configuration of the LV at the end of the diastolic phase, and $T>0$ is the final time. The parameters $\chi$ and $C_m\in \mathbb{R}^+$ are the ratio of membrane surface with respect to the volume and the membrane capacitance, respectively. In order to take into account the anisotropic electrical conductance ^[90], we define the diffusion tensor - a second order tensor in $\mathbb R^3$ - as ${\bf D}_m = \sigma_t \, {\bf I} + (\sigma_l - \sigma_t) \, {\bf f}_0 \otimes {\bf f}_0,$ where $\sigma_t, \sigma_l \in \mathbb{R}^+$ are the electric conductivities in the directions transversal and longitudinal with respect to the fibers, respectively. The current geometry displacement ${\bf d} = {\bf X}-{\bf x}$, with ${\bf d} \, :\, \Omega_0 \longrightarrow \mathbb R^3$, also influences the diffusive term since the second order strain tensor $\displaystyle{{\bf F} = {\bf I}+\frac{\partial {\bf d}}{\partial {\bf X}}}$ and $J = \det({\bf F})$, where ${\bf X}$ and ${\bf x}$ are the reference and deformed coordinates, respectively; for $J>0$, the second order tensor $\left(J {\bf F}^{-1} {\bf D}_m {\bf F}^{-T}\right)$ has uniform elliptic properties. The function $I^{app}(t)$ represents an externally applied current, which stands for the electric stimulus injected at the endocardium by the terminal fibers of the Purkinje network; for our purposes, we consider it as a source term which triggers the electrophysiological activity. The nonlinear term $I^{ion}(V, {\bf w})$ in Eq. 2.1 depends on the ionic variable ${\bf w}$ and is peculiar of the ionic model at hand. In the Hodgkin-Huxley formalism ^[45], a ionic model takes the form:

where the unknowns $w_i \in [0, 1]$ represent the ions concentration and/or the fraction of open ionic channels on the cellular membrane. Among the manychoices proposed in literature - see e.g. ^{[1,57,58,64,102]} – we use the Bueno-Orovio minimal model ^[12] for its simplicity (for which $N_I = 3$); nonetheless, this model is capable of capturing the main features of the electrophysiology in healthy myocardial tissues. The system of differential equations modeling the electrophysiology hence reads:

since $\chi = C_m = 1$ as prescribed in ^[12], where the monodomain equation is presented in dimensionless form; moreover, $\displaystyle I^{ion}(V, {\bf w} ) = \hspace{-1.0em} \sum_{\substack{q \in \{fi, so, si\}} } \hspace{-1.0em} I^q(V, {\bf w})$. Albeit a system of ODEs, the ionic model is defined in $\Omega_0 \times (0, T)$ since the dependence on $V$ indirectly introduces a dependence on the space independent variable. The terms of the ionic model are:

and $\displaystyle{I^{si}(V, w_2, w_3) = - \frac{H_{V_2}(V)}{\tau_{si}}w_2 w_3}$; moreover, $\tau^-_i(V) = H_{V^-_i}(V)(\tau"_i-\tau'_i) + \tau'_i$ for $i = 1, 2$, and $\tau_\eta(V) = H_{\eta}(V)(\tau"_\eta-\tau'_\eta) + \tau'_\eta$, with $H_\eta \in \left\{H_{V_2}, H_{V_o}, H_{V_{so}}^{k_{so}}\right\}$ for $\eta \in\{3, o, so\}$, respectively. Here $H_a(z)$ is the Heaviside function centered in $a\in\mathbb{R}$ and its smooth approximation is $\displaystyle{H_a^{\varepsilon}(z) = \frac{1 + \text{tanh}(\epsilon(z-a))}{2}}$, with $\epsilon \in \mathbb{R}^+$.

2.2.   Tissue passive mechanics

An adequate mechanical model for the description of the myocardium's displacement must account for the tissue's complex behavior. Firstly, the internal stresses induced by a prescribed deformation are highly anisotropic ^[41] and in our formulation depend on the directions ${\bf f}_0$, ${\bf s}_0$, and ${\bf n}_0$. We model the compressibility of the tissue through a nearly-incompressible formulation by weakly penalizing large volumetric variations ^[95] since with this approach small volumetric changes are allowed. Finally, the contraction of the muscle due to the electrophysiological activity must be taken into account: with this aim, we introduce the auxiliary dimensionless variable $\gamma_f$, which represents the relative stretching (or elongation) of the fibers. The model that we use to describe the evolution of $\gamma_f$ will be detailed in Sec. 2.3 as it constitutes the link between the electrophysiology and the mechanics.

We recall the momentum conservation equation in the reference configuration $\Omega_0$, endowed with boundary and initial conditions, in the unknown displacement variable ${\bf d}$ ^[66]:

Here $\Gamma_0^\eta$, with $\eta = \{epi, base\}$, represents the subsets of the boundary corresponding to the epicardium and the base of the myocardium as depicted in Figure 2. We denote by $\displaystyle{{\bf q}\left({\bf d}, \frac{\partial {\bf d}}{\partial t} \right) = ({\bf N} \otimes {\bf N}) \left( K_{\perp}^\eta {\bf d} + C_{\perp}^\eta \frac{\partial {\bf d}}{\partial t} \right) + ({\bf I} - {\bf N} \otimes {\bf N}) \left( K_{\parallel}^\eta {\bf d} + C_{\parallel}^\eta \frac{\partial {\bf d}}{\partial t} \right)}$, a vector in $\mathbb R^3$, the zero-th order term of the Robin boundary condition and by $K_{\perp}^\eta, K_{\parallel}^\eta, C_{\perp}^\eta, C_{\parallel}^\eta \in\mathbb{R}^+$ its parameters: the symbols $\perp$ and $\parallel$ identify either a parameter relative to the normal or the tangential direction at the corresponding boundary subset $\gamma_0^\eta$, respectively. $p^{endo}(t)$ is the external load applied by the fluid at the endocardium wall which, at this stage of the model description only, we assume to be prescribed. ${\bf N}$ is the outward directed unit vector normal to the boundary; ${\bf d}_0$ and $\dot{{\bf d}_0}$ are the initial data. The information related to the mechanical behaviour is embedded in the nonlinear Piola-Kirchhoff strain tensor ${\bf P} = {\bf P}({\bf d}, \gamma_f)$ - a second order tensor in $\mathbb R^3$ - which also must incorporate the active properties of the muscle. In order to characterize ${\bf P}$, we first define the right Cauchy-Green tensor as ${\bf C} = {\bf F}^T {\bf F}$, where ${\bf F} = {\bf I} + \nabla_0 {\bf d}$ is the strain tensor, and we introduce the strain energy function $\mathcal{W}({\bf C}):\mathbb{R}^{3 \times 3} \longrightarrow \mathbb{R}$ which relates the strain energy of the material to the strain tensor. The tissue mechanical properties are taken into account through the strain energy function: under the hyperelasticity assumption, the latter is differentiated with respect to the deformation tensor in order to obtain the Piola-Kirchhoff strain tensor, i.e.:

A very popular model for the myocardial tissue is the Holzapfel-Ogden ^[46] one. This is obtained by considering different contributions and by taking into account the anisotropic nature of the muscle:

where $\mathcal{I}_1 = \text{tr} \, {\bf C}$, $\mathcal{I}_{4f} = {\bf C :} \, {\bf f}_0 \otimes {\bf f}_0 = {\bf f} \cdot {\bf f}$, $\mathcal{I}_{4s} = {\bf C :} \, {\bf s}_0 \otimes {\bf s}_0 = {\bf s} \cdot {\bf s}$, and $\mathcal{I}_{8fs} = {\bf C :} \, {\bf f}_0 \otimes {\bf s}_0 = {\bf f} \cdot {\bf s}$ are the invariants of the tensor ${\bf C}$; the parameters $a_k, b_k$ are fitted from experimental data ^[46]. The function $\left<y\right> = y\, H_0(y)$ indicates the positive part of $y$ and its role consists in switching off the contributions of the fibers and sheets to $\widetilde{\mathcal W}$ when the material is under compression along these directions.

The mechanical model should also account for the volumetric change to which the myocardium undergoes during the cardiac cycle. It has been observed in ^[19,113] that albeit this change is moderate, still it significantly ranges from 2% to 15%. For these reasons we use a nearly-incompressible formulation, which allows for small volume variations ^[28]. We multiplicatively decompose the deformation gradient ${\bf F}$ into the isochoric and the volumetric parts as:

where $J = \det({\bf F}_v) = \det({\bf F}),$ being $\det(\overline{{\bf F}}) = 1,$ and we weakly enforce the incompressibility constraint by adding to the strain energy function $\widetilde{\mathcal{W}}$ a convex term $\mathcal{W}_{vol}(J)$ such that $J = 1$ is its global minimum; in this manner, large variations in volume are penalized. We choose:

such that the larger is the bulk modulus $B\in\mathbb{R}^+$, the "stronger" is the enforcement of the incompressibility constraint. Following ^[93], we evaluate the isotropic term $\mathcal{W}_1$ in $\mathcal{\overline{I}}_1 = \text{tr} (\overline{{\bf F}}^T \overline{{\bf F}}) = J^{-\frac{2}{3}} \mathcal{I}_1$ (instead of $\mathcal{I}_1$). The energy function hence reads:

We now proceed by modeling the active behavior of the myocardium driven by the stretching of the fibers by means of the active strain approach ^[2,3,63,86]. A virtual intermediate state $\widehat{\Omega}$, representing the active part of the deformation between the reference domain $\Omega_0$ and the deformed one $\Omega$ (see Figure 3), is introduced. The domain $\widehat{\Omega}$ is reached from $\Omega_0$ by applying a prescribed active transformation (which we will specify later) represented by the tensor ${\bf F}_A$. On the other hand, the material's elastic response to the prescribed active transformation is embedded in the tensor ${\bf F}_E$ and finally transforms $\widehat{\Omega}$ into $\Omega$. Mathematically, this approach requires a decomposition in the form ${\bf F} = {\bf F}_E {\bf F}_A = {\bf F}_v \overline{{\bf F}}_E {\bf F}_A = J^{\frac{1}{3}} \overline{{\bf F}}_E {\bf F}_A$, where we also take into account the factorization (2.5) for the term ${\bf F}_E$. Finally, we define ${\bf P}$ in $\Omega_0$ with respect to the total displacement ${\bf d}$ of the tissue by applying a pull-back to the stress computed in the intermediate state $\widehat{\Omega}$, i.e.:

In Sec. 2.3, we will provide the explicit form of the tensor ${\bf F}_A$ under the requirement of symmetry and identity of its determinant. Under these assumptions, the final form of the tensor ${\bf P}$ reads:

where ${\bf f}_E = {\bf F}_E {\bf f}_0$, ${\bf f}_A = {\bf F}_A^{-1} {\bf f}_0$, ${\bf s}_E = {\bf F}_E {\bf s}_0$, and ${\bf s}_A = {\bf F}_A^{-1} {\bf s}_0$.

2.3.   Mechanical activation

The activation model of the myocardium represents the link between the electrophysiology and the tissue (passive) mechanics. Instead of cellular models describing the complex dynamic taking place inside the sarcomeres ^[81], we exploit a phenomenological model for the local shortening of the fibers $\gamma_f$ at the macroscopic level. This model has been proposed in ^[88] and further developed in ^[84], and assumes that the evolution of $\gamma_f$ is due to the concentration of calcium ions $[Ca^{2+}]$ and by a feedback from the mechanics (through the variable ${\bf d}$):

Here, $g(c) = \widehat{\mu}_A c^2$, while the active force $\Phi(c, \gamma_f, {\bf d})$ reads $\Phi(c, \gamma_f, {\bf d}) = \alpha H_{c_0}(c) \, (c-c_0)^2 R_{FL}(I_{4f}) + \sum \limits_{j = 1}^5 (-1)^j(j+1)(j+2)I_{4f}\gamma_f^j$. With respect to the formulation ^[84] we added the diffusive term $\varepsilon \Delta \gamma_f$ in (2.8) to yield a model in the form of a PDE. While this is not strictly motivated by physical considerations, it can be interpreted as the upscaling of the microscopic activation at the macroscopic continuum level of the tissue. Moreover, this choice yields a more regular solution $\gamma_f$ in terms of the space variable ${\bf X}$, from which the numerical approximation will also benefit.

The contraction of the tissue is triggered by the calcium concentration exceeding a threshold value $c_0$. The parameters $\alpha$ and $\widehat{\mu}_A$ represent quantities to be properly tuned for the case under consideration, while $R_{FL}$ is the sarcomere force-length relationship ^[39] of cardiac cells which we represent as truncated Fourier series $R_{FL}(z) = \left( \frac{d_0}{2}+ \sum \limits_{n = 1}^3 \left[d_n \text{sin}(n l_0 z^{1/2}) + e_n \text{cos}(n l_0 z^{1/2}) \right] \right) \chi_{\left[SL_{min}, SL_{max}\right]}(z^{1/2})$. The calcium concentration is not explicitly represented in the set of variables ${\bf w}$ in the Bueno-Orovio model; however, the variable $w_3$ acts as a generic ion concentration, for which it can be interpreted as $c = [Ca^{2+}]$; this choice has been already made in literature e.g. in ^[84,85].

Being $\gamma_f$ the solution of (2.8), we choose the following orthotropic form for the tensor ${\bf F}_A$ ^{[3,6,73,84,85]}:

which is symmetric and we set $\gamma_n = \gamma_n(\gamma_f), \, \gamma_s = \gamma_s(\gamma_f, \gamma_n)$. In analogy with $\gamma_f$, these functions represent the local shortening (or elongation) of the tissue in the directions ${\bf s}_0$ and ${\bf n}_0$, respectively. Following ^[7], we define $\gamma_n$ such that, according to experimental observations ^[67], the thickening of the ventricle's walls is transversely non-homogeneous:

Indeed, the latter depends on an independent variable $\lambda$ which represents the transmural coordinate, taking the value $\lambda_{endo}$ at the endocardium and the value $\lambda_{epi}$ at the epicardium. As detailed in ^[7] - where the results of numerical simulations are compared and validated against experimental data for strains in the canine LV provided in ^[67] - this model is able to capture the transmural heterogeneity of the thickening of the LV. The function $k'(\lambda)$ is defined as:

Finally, in order to have $\det \left( {\bf F}_A \right) = 1$, we set $\displaystyle{\gamma_s = \frac{1}{(1+\gamma_f)(1+\gamma_n)} - 1}$.

2.4.   The coupled model: electromechanics

Writing together Eqs.(2.3), (2.4), and (2.8), the coupled electromechanics problem finally reads:

We remark that at this stage the pressure $p^{endo}(t)$ is still prescribed: however, later in Sec. 4.1, we will consider $p^{endo}(t)$ as an additional unknown of the coupled problem, being the solution of 0D fluid problems for each of the phases of the cardiac cycle.

2.5.   Prestress

A common problem which has to be considered in biological modeling involving a fluid-structure interface is that the reference geometry $\Omega_0$, acquired from medical images at the telediastole (the phase immediately before the systole), does not necessarily correspond to a stress-free configuration. This because the blood exerts a pressure on the endocardium walls $p^{endo}(t)$ which, at each time instant of the heartbeat, is larger than zero, taking values ranging approximately from a minimum of 5 mmHg to a maximum of 120 mmHg in healthy individuals. In turn, this implies that solving problem (2.4) with a physiological endocardial pressure $p^{endo}>0$ would give rise to non-physiological displacements as the internal stresses are not in equilibrium with the intraventricular blood's pressure of LV. This issue is particularly evident in the case of patient-specific settings since, unlike the case of idealized geometries, the actual computational geometry is obtained from medical images and hence the reference domain $\Omega_0$ cannot be arbitrarily chosen. Typically, the LV geometry is acquired at the so-called telediastole and the electromechanics (or mechanics) model initiated at this stage. The corresponding pressure is indicated with $\overline{p}^{endo}$ and the stressed LV configuration is determined in these conditions. To our knowledge, two strategies have been proposed in literature to address this issue.

Pressure preload (^{[30,84,101,103]}): the reference geometry $\Omega_0$ is loaded with the prescribed pressure $\overline{p}^{endo}$. This is done by solving the steady variant of the mechanics problem (2.4) while gradually increasing the pressure until the desired value $\overline{p}^{endo}$. The displacement field so obtained is then used as an initial datum ${\bf d}_0$ for the unsteady problem.

Pressure prestress (^[47,100]): we compute internal stresses distribution such that the reference geometry is in equilibrium with the blood pressure $\overline{p}^{endo}$. An additive decomposition of the stress tensor $\widetilde{{\bf P}} = {\bf P}({\bf d}) + {\bf P}_0$ is operated, where the prestress tensor ${\bf P}_0$ is determined to ensure a null displacement ${\bf d}_0$ in correspondance of the assigned pressure $\overline{p}^{endo}$.

We adopt the pressure prestress approach since in the preload one the procedure returns a configuration which might be significantly different with respect to the initial geometry as shown in ^[47]. Indeed, we already know the loaded configuration from medical images. In order to compute ${\bf P}_0$ according to this approach, we proceed by adapting the method proposed in ^[47] to our model by first defining the following mechanical problem:

Eq. (2.11) is the steady, passive version of problem (2.4) with the additive decomposition of the strain tensor. Moreover, we set:

where $S \in \mathbb{N}$ is the number of steps of the continuation method that we exploit to gradually increase the pressure value. Then, the fixed point Algorithm 1 is applied to compute ${\bf P}_0$. First, we compute the approximation $\widetilde{{\bf P}}_0 = \,$PRESTRESS$(100, 10^{-2}, \overline{p}^{endo}, {\bf 0})$ and finally we set ${\bf P}_0 = \,$PRESTRESS$(1, 10^{-5}, \overline{p}^{endo}, \widetilde{{\bf P}}_0)$. The additional step is performed to require a smaller tolerance when the pressure has already reached a closer value of the tensor $\widetilde{{\bf P}}_0$ to the target. We observe that, while $\frac{\|{\bf P}\left({\bf d}_k^m, {\bf I}\right)\|_{\infty}}{\|{\bf P}^m_{0, k}\|_{\infty}}\longrightarrow 0$ for $m\longrightarrow+\infty$ in Algorithm 1, the displacement ${\bf d}^{m+1}_k$ does not converge to ${\bf 0}$ but to a vector which we denote by $\widehat{{\bf d}}$. However, we observe that the quantity $\|\widehat{{\bf d}}\|_{\infty}$ is negligible with respect to the endocardial walls thickness, and that this initial displacement ensures that the prestress is in equilibrium with the pressure at the epicardium. Therefore, we decide to set ${\bf d}_0 = \widehat{{\bf d}}$ in (2.10).

3.   Numerical discretization

In this section we describe the numerical approximation of the electromechanics problem (2.10). The FEM ^[78] is used for the space discretization of the PDEs, while BDF ^[78] are used for the time discretization.

3.1.   Space discretization

We consider a mesh composed of tetrahedrons $\mathcal{T}_h$, with $h$ representing the maximum size of the elements $K \in \mathcal{T}_h$, such that $\cup_{K \in \mathcal{T}_h} K = \Omega_0$; the mesh elements are pairwise disjoint and their union $\Omega_0 \subset \mathbb{R}^3$ is the region of the space identified by the myocardium in the telediastolic phase of the heartbeat. In this work, we consider two LV geometries, that is an idealized prolate ellipsoid (as is often done in literature ^[30,40,84]) and a patient-specific geometry extracted through a segmentation procedure from 3D MRI images^*, as we will describe in Sec. 5.1. We will approximate each of the single core models in the computational domain with a space discretization based on the FEM. We will denote by $N^{\text{dof}}_V, N^{\text{dof}}_{{\bf w}}, N^{\text{dof}}_{{\bf d}},$ and $N^{\text{dof}}_{\gamma_f}$ the number of Degrees of Freedom (DoF) (i.e. the size of the discretized single core model) for the potential, ionic variables, displacement, and fiber shortening, respectively. The underlying number of nodes determined by the mesh $\mathcal{T}_h$ and the degree $r$ is indicated as $N_h$. Hence, we have that, for scalar PDEs as the electrophysiology, $N^{\text{dof}}_V = N_h$.

^*The MRI images are provided by Prof. J. Schwitter (Chief physician at the Centre Hospitalier Universitaire Vaudois CHUV, Lausanne) and Dr. P. Masci in the framework of the collaboration CMCS@EPFL-CHUV hospital.

3.1.1.   Monodomain equation

We use the FEM for the spatial approximation of the monodomain equation (2.1). We first introduce the finite dimensional space $\mathcal{X}^r_h = \left\{ v \in C^0(\overline{\Omega}_0) \; : \; v|_{K} \in \mathbb{P}^r(K) \; \forall K \in \mathcal{T}_h \right\}$, where $\mathbb{P}^r(K)$ is the set of polynomials of degree smaller than or equal to $r$ in the element $K$. Moreover, we denote by $\left\{ \psi_i \right\}_{i = 1}^{N^{\text{dof}}_V}$ a basis of $\mathcal{X}^r_h$ with $N^{\text{dof}}_V = \text{dim}\left(\mathcal{X}^r_h\right)$. The projection of the solution $V(t)$ on $\mathcal{X}^r_h$ can hence be written as $V_h(t) = \sum_{j = 1}^{N^{\text{dof}}_V} V_j(t) \psi_j$ and the weak semidiscrete formulation of the problem reads: given ${\bf w}_h(t)$ and ${\bf d}_h(t)$, find, for all $t \in (0, T)$, $V_h(t) \in \mathcal{X}^r_h$ such that

with $V_h(0) = \sum_{j = 1}^{N^{\text{dof}}_V} \left(V_0, \psi_j\right)_{L^2(\Omega_0)} \psi_j$. By setting ${\bf V}_h(t) = \left\{ V_j(t) \right\}_{j = 1}^{N^{\text{dof}}_V}$ and ${\bf V}_{0, h} = \left\{\left(V_0, \psi_j\right)_{\mathcal{H}} \right\}_{j = 1}^{N^{\text{dof}}_V}$, we rewrite Eq. (3.1) as a system of nonlinear ODEs:

where $\mathbb{M}_{ij} = \int_{\Omega_0} \psi_j \psi_i \, d\Omega_0$, $\mathbb{K}_{ij}({\bf d}_h) = \int_{\Omega_0} (J {\bf F}_h^{-1} {\bf D}_m {\bf F}_h^{-T} \nabla \psi_j) \cdot \nabla \psi_i \, d\Omega_0$, and

We will discuss in Sec. 3.1.3 two different strategies for the approximation of the nonlinear term ${\bf I}_{ion}({\bf V}_h, {\bf w}_h)$. Moreover, we will use a lumped mass matrix $\mathbb{M}^L$ in place of $\mathbb{M}$ in Eq. (3.2) in order to avoid numerical instabilities, a known numerical issue that may occur in the case of problems with "traveling waves" due to the dominance of zero order terms over the second order ones ^[13].

3.1.2.   Ionic model

The ionic model (2.2) is a system of ODEs that indirectly depends on the space variable through the transmembrane potential $V$. By denoting with $\left\{{\bf x}_j\right\}_{j = 1}^{N^{\text{dof}}_{{\bf w}}}$ the set of the degrees of freedom, we write the equations of the model evaluated at each of these points and we denote the value of the $l$-th ionic variable in ${\bf x}_j$ by $w^l_j(t)$. Similarly, we write $V_j(t) = V_h({\bf x}_j, t)$, and finally rearrange the unknowns in the vector ${\bf w}_h(t)$ in the following fashion:

The problem thus obtained reads: given $V_h(t)$, find, for all $t \in (0, T)$, ${\bf w}_h(t)$ such that

with $w^l_j(0) = w^l_0$, for $l = 1, \dots, N_I, \; j = 1, \dots, N^{\text{dof}}_{{\bf w}}$. Following Eq. (3.4), system (3.5) can be conveniently rewritten in algebraic form as

where the block matrix $\mathbb{U}$ and the block vector ${\bf Q}$ are defined as $\left(\mathbb{U}(V_h)\right)_{mm} = \alpha_l(V_j) - \beta_l(V_j)$ and $\left({\bf Q}(V_h)\right)_m = \alpha_l(V_j)w_l^{\infty}(V_j)$, respectively, where $m = l\, N^{\text{dof}}_{{\bf w}} +j$, for $l = 1, \dots, N_I, \; j = 1, \dots, N^{\text{dof}}_{{\bf w}}$.

3.1.3.   Ionic currents

As anticipated in Sec. 3.1.1, different choices can be made for the approximation of the term ${\bf I}_{ion}({\bf V}_h, {\bf w}_h)$ in Eq. (3.1). By denoting with $\left\{{\bf x}^K_s\right\}_{s = 1}^{N_q}$ and $\left\{\omega^K_s\right\}_{s = 1}^{N_q}$ the quadrature nodes and the corresponding weights in a generic mesh element $K \in \mathcal{T}_h$, we consider two approaches which are investigated, for instance, in ^{[53,68,69,70]}.

State Variable Interpolation (SVI): the variables $V_h$ and ${\bf w}_h$ in Eq. (3.3) are evaluated at the quadrature nodes

This approach corresponds to the standard FEM approximation.

Ionic Currents Interpolation (ICI): for each element, $I^q$ is evaluated at the degrees of freedom of $V$ and its value at the quadrature nodes is obtained by interpolation through the same Lagrangian polynomial basis of degree $r$ used for the FE space.

The two approaches are equivalent if $I^q$ is linear in $V_h$ and ${\bf w}_h$. Both the approaches tend to overestimate the conduction velocities if the computational mesh is not "sufficiently" fine ^[62], but yield convergent results under h-refinement. The SVI approach corresponds to a standard FE approximation, while the ICI one yields a faster assembly of the ionic currents term ^[53,68,84] but introduces an error in the evaluation of ${\bf I}_{ion}({\bf V}_h, {\bf w}_h)$. In our case, as we will describe in Sec. 5, the overall time required for the assembly of the terms of the system arising from the ful discretization of Eq. (2.10) is mostly determined by the construction of the mechanics terms, while the assembly of the ionic currents term does not represent a computational bottleneck; we hence use the SVI method since the advantage of using ICI here is negligible.

3.1.4.   Active and passive mechanics

As previously done for the monodomain equation, we use the FEM to approximate the momentum equation (2.4). We denote by $[\mathcal{X}^r_h]^3$ the finite dimensional subspace of vector valued functions and by $\{\mathit{\boldsymbol{\psi}}_i\}_{i = 1}^{N^{dof}_{{\bf d}}}$ its basis. The Galerkin formulation of Eq. (2.4) reads: given $\gamma_{f, h}(t)$, find, for all $t \in (0, T)$, ${\bf d}_h(t) \in [\mathcal{X}^r_h]^3$ such that

with ${\bf d}_h(0) = \sum_{\eta = 1}^{N^{\text{dof}}_{{\bf d}}} \left({\bf d}_0, \mathit{\boldsymbol{\psi}}_j\right)_{[L^2(\Omega_0)]^3} \mathit{\boldsymbol{\psi}}_j$ and $\dot{{\bf d}}_h(0) = \sum_{\eta = 1}^{N^{\text{dof}}_{\dot{{\bf d}}}} \left(\dot{ {\bf d}}_0, \mathit{\boldsymbol{\psi}}_j\right)_{[L^2(\Omega_0)]^3} \mathit{\boldsymbol{\psi}}_j$.

We rewrite it in algebraic form as:

where $\mathbb{F}_{ij} = \sum_{\substack{\eta \in \{epi, base\}} }\int_{\Gamma_0^\eta} \left( C_{\perp}^\eta ({\bf N} \otimes {\bf N}) + C_{\parallel}^\eta ({\bf I} - {\bf N} \otimes {\bf N}) \right) \mathit{\boldsymbol{\psi}}_j \cdot \mathit{\boldsymbol{\psi}}_i \, d\Gamma_0$ and $\mathbb{G}_{ij} = \sum_{\substack{\eta \in \{epi, base\}} } \int_{\Gamma_0^\eta} \left( K_{\perp}^\eta ({\bf N} \otimes {\bf N}) + K_{\parallel}^\eta ({\bf I} - {\bf N} \otimes {\bf N}) \right) \mathit{\boldsymbol{\psi}}_j \cdot \mathit{\boldsymbol{\psi}}_i \, d\Gamma_0$; moreover, we have ${\bf S}_i({\bf d}_h(t), \gamma_{f, h}(t)) = \int_{\Omega_0} {\bf P}({\bf d}_h, \gamma_{f, h}) : \nabla_0 \mathit{\boldsymbol{\psi}}_i \, d\Omega_0$.

Finally, we once again use FEM to discretize in space the equation for the unknown $\gamma_f$: given $c_h(t)$ and ${\bf d}_h(t)$, find, for all $t \in (0, T)$, $\gamma_{f, h}(t) \in \mathcal{X}_h$ such that

with $\gamma_{f, h}(0) = 0$. We highlight that, with respect to Eq. (2.8), we divided each term by $g(c_h)$. Thus, we obtain the following system of ODEs:

3.2.   Time discretization

After having applied the space discretization to the single core models, we obtain the semi-discretized formulation of problem (2.10) in the form of a nonlinear system of ODEs. We denote by ${\bf z} = ({\bf z}_{{\bf w}}, \, {\bf z}_V, \, {\bf z}_{\gamma_f}, \, {\bf z}_{{\bf d}})^T$ the block vector containing the unknowns of each semi-discrete single core problem, where the ordering of the unknowns will be motivated in Sec.3.2.2, and we write:

where $\displaystyle{\mathcal{M} = \text{diag}\left(\frac{\text{d}}{\text{dt}}, \frac{\text{d}}{\text{dt}}, \frac{\text{d}}{\text{dt}}, \frac{\text{d}^2}{\text{dt}^2}\right)}$ is a differential operator which applies a first order time derivative to the ionic variables, the transmembrane potential and the fibers shortening (motivating the division by $g(c_h)$ of the corresponding equation), while a second order time derivative to the displacement. In order to obtain a fully discretized formulation, we use the BDF for the time approximation of Eq. (3.6). Hence, we write:

where $\Delta t = \frac{T}{N_T}$ is the timestep, $N_T$ being the number of timesteps, while the parameters $\vartheta^{(Ⅰ)}_k, \vartheta^{(Ⅱ)}_k, k = 0, \dots, \sigma$ depend on the order $\sigma$ of the BDF scheme. We will in particular use BDF of order $\sigma = 1$ and 2, and will consider two schemes for the time discretization of the monolithic problem: a fully implicit and a semi-implicit one.

In the implicit case, we obtain the following nonlinear system:

with ${\bf z}^k$ assigned for $k = 0, \dots, \sigma$ and we set for simplicity $\dot{{\bf d}}_{0, h} = {\bf 0}$. Problem (3.7) is solved by exploiting the Newton method ^[78] at each timestep.

In the semi-implicit case, on the other hand, we extrapolate the variables in the nonlinear term ${\bf A}({\bf z}^{n+1})$ by means of the Newton-Gregory backward polynomials ^[16] - as is done, e.g., for the Navier-Stokes equations in ^[37] - thus yielding a linear system at each timestep. The extrapolated variables are evaluated as a linear combination of the variables at previous timesteps, with an approximation of the same order $\sigma$ of the BDF scheme:

We then avoid the nonlinear terms by partially evaluating ${\bf A}$ in the extrapolated variable ${\bf z}^*$, i.e. we approximate the nonlinear term as

By recalling Eq.(3.7), we hence obtain a system in the form:

with ${\bf z}^k$ assigned for $k = 0, \dots, \sigma$ and $\widetilde{{\bf b}}^{n+1} = {\bf h}^{n+1}-{\bf g}^{n+1}$.

3.2.1.   The implicit scheme

By applying the Newton method ^[78] to (3.7), we iteratively solve for $n = \sigma, \dots, N_T-1$ the linear problem

for $k = 0, \dots$, until some convergence criterion is met. $\mathbb{J}^{n+1}_k$ is the Jacobian matrix of ${\bf A}$ evaluated in the ${\bf z}^{n+1}_k$ and is endowed with the following block structure:

while the residual is defined as ${\bf r}_k^{n+1} = {\bf b}^{n+1} - {\bf A}({\bf z}_k^{n+1})$.

3.2.2.   The semi-implicit scheme

In this case, the block structure of the matrix associated to linear system (3.8) to be solved is sparser with respect to the implicit one, as $\mathbb{A}$ is a block lower triangular matrix endowed with only one extra diagonal block for each $n = \sigma, \dots, N_T-1$:

The matrix (3.10) is block lower triangular thanks to the ordering of the variables in the electrophysiology model (3.6): indeed, by swapping the ionic and the electric blocks we would have obtained a block upper triangular matrix. The semi-implicit scheme has the clear advantage of avoiding, at each timestep, the solution of a nonlinear problem. Moreover, both the assembly of the system matrix $\mathbb{A}({\bf z}^*)$ and the solution of the linear system (3.8) require smaller computational time than in the implicit case.

On the other hand, semi-implicit schemes may impose strong limitations on the timestep size $\Delta t$ to ensure stability. This is indeed what we observe from our numerical tests: the maximum timestep size that we used in order to ensure stability for the semi-implicit scheme is at least one order of magnitude smaller than the one that used for the implicit one, for which accuracy is the main factor driving the choice of $\Delta t$. This, in turn, makes the computational cost of the semi-implicit scheme much larger than the implicit one, especially when considering the simulation of a full heartbeat. By analyzing the behavior of the linear solver (which we will detail in Sec. 4) and the residual of the system (3.8), we observe that numerical instabilities are driven by the mechanics block: we conclude that the strong nonlinearity of the Piola-Kirchhoff tensor (2.7) and the presence of the exponential terms in particular yield evaluations of the stress tensor ${\bf P}({\bf d}_h^*, \gamma^*_{f, h})$ in the extrapolated variables that are not "sufficiently" close to ${\bf P}({\bf d}_h^{n+1}, \gamma^{n+1}_{f, h})$, unless the timestep is significantly "small". Indeed, instability does not appear in numerical simulations of the isolated electrophysiology problem, even using a larger timestep with respect to the one used for the monolithic problem, thusconfirming that the bottleneck in the choice of $\Delta t$ is due to the mechanics core model.

We will further investigate this issue in the future; however, for the next numerical simulations of this paper we will consider only the implicit scheme.

4.   Linear solver: the preconditioning strategy

We rewrite, for the sake of simplicity, the linear system associated to a single Newton iteration for the implicit scheme (3.9) in the following form:

We use the GMRES method ^[89] for the solution of the linear problem (4.1) as $\mathbb{J}$ is non-symmetric and the distribution of the eigenvalues of its spectrum is not known a priori.

The choice of the preconditioner is critical in order to ensure the convergence of the linear solver; while this is true in general, it is even more relevant in the case of monolithic multiphysics problems ^[51]. Indeed, using a black-box algebraic preconditioner for problem (4.1) the information related to each differential problem associated to a single core model would be neglected; we instead consider a strategy exploiting such information at the block level, that is we create a preconditioner that exploits the "physics" of the coupled problem at the level of the block structure. Following ^[26,27,35] for FSI problems, we propose a block Gauss-Seidel preconditioner $\mathcal{P}$ obtained by dropping the upper triangular blocks of matrix $\mathbb{J}$, i.e. the off-diagonal blocks, thus obtaining:

$\mathcal{P}$ can then be factorized into four model-specific nonsingular matrices, namely $\mathcal{P}_{ion}$, $\mathcal{P}_{pot}$, $\mathcal{P}_{act}$, and $\mathcal{P}_{mec}$ corresponding to the ionic, the potential, the activation, and the mechanics single core models, respectively:

This decomposition can also be interpreted as an inexact factorization of matrix $\mathbb{J}$. The preconditioned version of problem (4.1) then reads:

Since each diagonal block $J_{ii}$ appears in a distinct factor $\mathcal{P}_i$, for $i = 1, \dots, 4$, then physics-specific ad-hoc preconditioners can be efficiently used to approximate its inverse. Indeed, we define the symbolic operation $\delta {\bf z} = \mathcal{P}^{-1}{\bf y}$ in (4.2) as the application of the following sequential steps

the solution of each linear system in (4.3) is carried out by using again the GMRES method and by exploiting Algebraic Multigrid (AMG) preconditioners ^[11]. This new preconditioner $\mathcal{P}$ can be regarded as a generalization of the FaCSI preconditioner for FSI problems proposed in ^[27] (see also ^[26,35]).

4.1.   Cardiac cycle

For our simulations we will consider a full heartbeat, by taking the conventional duration of $T = 0.8$ s. With this aim, we take into account for the interaction of the endocardium with the blood by modeling the pressure $p^{endo}$ as in Eq. (2.10). Before introducing the models used to describe such pressure, we first recall that the heartbeat is conventionally split into four phases (see Figure 4):

(1) an isovolumic contraction phase in which the endocardial pressure $p^{endo}(t)$ increases from the End Diastolic Pressure (EDP) to the value measured in the aorta (about 95 mmHg). This increment is driven by the early stages of the LV contraction;

(2) an ejection phase characterized by a decrement in the ventricular volume $V^{endo}(t)$ due to the contraction of the LV forcing the blood to flow through the aortic valve;

(3) an isovolumic relaxation phase during which $p^{endo}(t)$ decreases as a consequence of the LV early relaxation;

(4) a filling phase starting with the opening of the mitral valve causing an increment of $V^{endo}(t)$ until $p^{endo}(t)$ reaches the EDP value.

When necessary, we compute the volume $V^{endo}(t)$ at time $n$ by exploiting the formula $V^{endo, n} = \int_{\Gamma_0^{endo}} J({\bf d}_h^n)\mathit{\boldsymbol{\xi}}\cdot{\bf F}^{-T}({\bf d}_h^n){\bf N} \, d \Gamma_0$, which is rigorously derived in ^[84] and where $\mathit{\boldsymbol{\xi}}$ is a vector directed as the centerline of the LV.

We then model the endocardial pressure $p^{endo}(t)$ with different 0D models, following ^[30,84,107], for each of the aforementioned phases (in the following, we drop the "$endo$" superscript for simplicity):

1) Isovolumic phase I: At each timestep, we iteratively solve problem (2.10) by updating at each subiteration the pressure in the following manner:

with $p^{n+1}_0 = p^n$, $V^{n+1}_0 = V^n$, until the condition $\frac{|V^{n+1}_k-V^n|}{V^n} < \varepsilon$ is satisfied. The parameter $C_p < 0$ has to be "sufficiently" large in order for the fixed point algorithm to converge. This phase ends when the pressure $p^{n+1}_{k+1}$ reaches the value $\overline{p} = 95$ mmHg.

2) Ejection phase: A two-element Windkessel model ^[110] of the form:

is solved in the pressure variable with a BDF scheme of order $\sigma = 1$ while the term $\frac{dV}{dt}$ is approximated at time $n+1$ as $\frac{V^{n}-V^{n-1}}{\Delta t}$, for simplicity. $T'$ and $T"$ are the initial and final time of this phase, respectively, while the parameters $C, R>0$ represent the capacitance and resistance of the equivalent electric circuit. This phase ends when the (initially negative) term $\frac{dV}{dt}$ changes sign.

3) Isovolumic phase II: Modeled as the first isovolumic phase. This phase ends when the pressure reaches the minimal value of $\overline{\overline{p}} = 5$ mmHg.

4) Filling phase: The pressure is linearly increased at each timestep so that it reaches the EDP value at the final time $T$. We are aware that this is not fully coherent with the real behavior of the diastolic phase, but a model for the muscle and the function of the left atrium are not taken into account in this work. Future work should better address this points to represent the full heartbeat in a detailed fashion.

Even if (4.5) is strongly coupled with (2.10) during the ejection phase, for simplicity we solve the two problems in a segregated fashion.

5.   Numerical simulations

We now describe the geometries used for the numerical simulation of the integrated electromechanics problem. As anticipated in Sec. 3, we use both an idealized prolate geometry and a patient-specific geometry segmented from MRI images as detailed in Sec. 5.1. Moreover, in Sec. 5.2, we outline the procedure that is exploited to define the fibers and the sheets fields on the different meshes. We conclude this section by showing and analyzing the results of the simulations.

5.1.   Patient-specific geometry segmentation

Image segmentation is the process of locating regions of interest (ROI) in the form of a subset of pixels ^[42]. In biomedical applications, this accounts to assign different flags to regions containing different types of tissues and/or fluids. This result, depending on the properties of the biological material which is to be segmented, can often be achieved through semi-automatic or automatic procedures (see e.g. ^[32,49,106] for arteries and blood vessels and ^[10,109] for the Purkinje network), exploiting a large set of different techniques. However, since the development of algorithms for the segmentation of the chambers of the heart ^[72,114,115] is beyond the scope of this work, we used a manual procedure based on the brightness of the pixels of the aforementioned 3D MRI images. We first apply a threshold filter to select the smallest set of pixels containing the whole myocardium; then, we manually remove artifacts and features that we decide to neglect such as the papillary muscles and the upper part of the LV. Finally, a Gaussian smoothing procedure ^[111] was performed to improve the quality of the mesh. See Figure 5.

The image was taken at the end of the diastolic phase, that is when the LV reaches its maximum enlargement. The internal volume of the patient-specific myocardium at this stage of the heartbeat measures approximately 95 ml.

5.2.   Fibers and sheets distribution

Unlike the geometry, the fibers and sheets field distribution in the myocardium may not be extracted from medical images such as MRIs unless special procedures are applied ^[80]. For this reason, several mathematically rule-based definition of the fields have been used in literature ^{[38,56,61,86]}, which try to approximate their orientation. At the epicardium and at the endocardium, the fiber direction is tangential to the boundary, while the sheet direction belongs to the plane identified by the normal and the ventricle centerline. Then, in the most general case, angles $\alpha_{endo}$, $\alpha_{epi}$, $\beta_{endo}$, and $\beta_{epi}$, representing the inclination of the fibers and the sheets with respect to the base plane, are assigned. Finally, the direction of fibers and sheets inside the myocardium is determined by a transmurally linear mapping. A first study of the influence of the angles on both the conductivity and the deformation was carried out in ^[30]. Since this kind of analysis goes beyond the scope of this work, we will consider for both the idealized and the patient-specific geometries the algorithm proposed in ^[112] and further developed in ^[84], thereforse setting $\alpha_{endo} = -60^{\circ}$, $\alpha_{epi} = +60^{\circ}$, $\beta_{endo} = \beta_{epi} = 0^{\circ}$. In Figure 6, we show the fields obtained by applying the algorithm to a set of subsequently refined ideal meshes, while in Figure 7 the same is done for the patient-specific mesh.

5.3.   Numerical results

We setup a total of five numerical simulations, each for a single heartbeat, using the four meshes shown in Figures 6-7. The monodomain, the activation, and the mechanics equations are approximated with $\mathbb{P}^r$ elements; the ionic model equations, on the other hand, are solved in the degrees of freedom determined by the discretization of the monodomain equations (e.g. the vertices of the tetrahedrons if $r = 1$). We set $r = 1$ for each mesh in order to limit the size of the monolithic system and, additionally, $r = 2$ for the coarse mesh only. This yields a linear system of size $M = 8\times N_h$. A second order BDF scheme ($\sigma = 2$) is considered for the time discretization in all cases to ensure A-stability while maximizing the convergence order ^[78]. The information on the meshes and on the numerical simulations is summarized in Tables 1 and 2, respectively.

We use LifeV^†, an open-source finite element library for the solution of problems described by PDEs in a High Performance Computing framework, to implement the numerical methods. All the computations were carried out using Piz Daint, a Cray XC50/XC40 supercomputer installed at the Swiss National Supercomputing Center (CSCS)^‡.

^†https://cmcsforge.epfl.ch

^‡http://www.cscs.ch

As detailed in Sec. 4 we exploit the Newton method for the solution of the nonlinear system (3.7) by setting a maximum number of 5 iterations per timestep. Moreover we use the relative residual as stopping criterion, with a tolerance equal to $10^{-7}$. At each Newton iteration (3.9), as detailed in Sec. 4, we use the GMRES method for the inversion of the blocks $J_{ii}$, $i = 1, \dots, 4$ together with AMG preconditioners by taking into account the number of underlying PDEs of each block; we exploit the ML package ^[36] of the Trilinos library ^[43] by performing three sweeps of the Gauss-Seidel algorithm for pre- and post-smoothing, while the solution on the coarsest level is obtained through an LU factorization. We consider the relative residual as stopping criterion for the GMRES method, with a tolerance equal to $10^{-8}$.

For simplicity, we initiate the electric signal propagation by applying at the same time 2-milliseconds long stimuli at three distinct points on the endocardium; we remark that, even if this pacing strategy is not fully realistic since the electrical signal delivered by the Purkinje network activates asynchronously hundreds of cells ^[83,108], it provides physically meaningful results. The other parameters used for the simulations are reported in Table 3.

We compare in Figure 8 the transmembrane potential obtained using $r = 1$ elements on the three idealized meshes at different times. As expected, the conduction velocities are overestimated when the mesh size is not sufficiently fine ^[9,29]. The same conclusion can be drawn by observing the plots in Figure 9, where the transmembrane potential sampled at the apex is represented against time: the delay of the action potential (the quick depolarization of the tissue followed by a slow repolarization ^[21]) in the finer meshes is due to the longer time needed for the wave to reach the apex. We display in Figure 10 the transmembrane potential for the patient-specific case at the same instants.

We observe from Figure 9 a difference of several milliseconds in the complete activation of the myocardium, which is significant in the electrophysiology characteristic time scales. However, because of the multiscale nature of the integrated problem, the mechanics is not affected by such delay as can be observed in Figure 11. Here, the only remarkable difference between the three results is the underestimation of the displacement magnitude in the simulation with the coarsest mesh (a). In this case, the number of degrees of freedom is not sufficient to either fullfill the nearly-incompressible constraint nor to represent the large deformations observable in the other two cases. The deformation of the patient-specific mesh (d) at the same times is shown in Figure 12.

In both the idealized and the patient-specific cases, a significant thickening of the myocardium walls takes place, which is in accordance with experimental observations ^[79]. The model, with the set of parameters in Table 3, produces however a significant rotation of the LV: a recent work ^[74] suggests that this behavior is related to the choice of the incompressibility constraint, the bulk modulus $B$ magnitude, and the boundary conditions. As a first step towards the investigation of this aspect, we report in Figure 13 the muscle obtained with different values of the bulk modulus $B$, while all the other parameters are kept unchanged with respect to Table 3 Larger values of $B$ correspond to a larger torsion, more accentuated and uniform along the ventricle's walls; moreover, the thickening of the myocardium is larger when the bulk modulus is smaller as expected, since the incompressibility condition becomes weaker. However, the imposition of a stronger volumetric constraint through $B$ negatively affects the conditioning of the system matrix $\mathbb{J}$ in (4.1), and forces the linear solver to perform more iterations in order to reach convergence. We will further focus on this issue in future works.

In order to better evaluate the ability to describe the LV activity of the integrated numerical model for a full cardiac cycle, pressure-volume (pV) loops are also represented. A comparison with in-vivo measurements would be pointless in the case of the ideal geometry, however we stress the fact that the pV-loops reported in Figure 14 are in both qualitative and quantitative agreement with those observed experimentally as e.g. in ^[87]; moreover, the maximal endocardial pressure, which is reached during the ejection phase, is approximately 110 mmHg and hence lays well within the physiological range of average individuals ^[91]. The difference between the four cases in Figure 14 accounts to a maximum variation of about 2 % in the minimal ventricle volume during the second isochoric phase.

The pV-loop for the patient-specific case is reported in Figure 15. Experimental measurements of the pressure and of the volume for the patient under study were not available, therefore we did not carry on a quantitative comparison of our results against patient-specific data. However, we once again observe that the values of the pV-loop are in line with real data ^[87].

We then report in Table 4, the timestep used for each simulation, the average number of Newton iterations, the average number of fixed point iterations for each of the two isochoric phases, and the average time spent for different operations during a single timestep. The overall time required for the complete simulations is also reported. Table 4 clearly shows that finer meshes require larger computational resources for a full heartbeat simulation ($T^{tot}$), and this is mostly due to the longer time required for the solution of a single linear system. On the other hand, we notice that the proposed preconditioner ensures that the convergence of the GMRES method is obtained with a relatively small number of iterations $\overline{N}^G$ even with large linear system sizes.

Finally, we test the strong scalability of the proposed monolithic solver for cardiac electromechanics. With this aim, we further refine the ideal mesh of Figure 6, thus obtaining a finer one featuring 4'629'817 vertices and 26'624'000 tetrahedra, which is more adequate for the strong scalability test. We then set $\Delta t = 2 \times 10^{-4}$ s and solve the electromechanics problem with the implicit scheme up to the final time $T = 10^{-2}$ s. As preconditioner for the GMRES linear solver, we use the block Gauss-Seidel one $\mathcal P$ outlined in Sec. 4. However, for this test we use an Additive Schwarz preconditioner for preconditioning of the mechanics core block as, in our experience, it shows better scalability properties that the Algebraic Multigrid one ^[34]. For our test, we use 800, 1'600, and 3'200 CPUs. Since with this mesh the total number of unknowns of the monolithic discrete problem (3.9) amounts to 37'038'536, the number of unknowns attributed to each CPU is about 46'298, 23'149, and 11'574, respectively^§. We remark however that our preconditioning strategy applies sequentially to each core model according to Eq.(4.3); hence, the true number of unknowns per CPU is smaller and amounts, e.g., to about 5'787, 2'894, and 1'447 for the monodomain model, respectively. We report in Fig. 16 the mean number of Newton $\overline{N}^N$ and $\overline{N}^G$ GMRES iterations, the average time $\overline{T}^P$ spent in assembling and applying the preconditioner, and the total wall time $T^W$ of the simulation. We observe that, while $\overline{N}^N$ and $\overline{N}^G$ are kept almost constant with the number of CPUs, the total wall time $T^W$ first decreases but then increases when the largest number of unknowns per CPU is achieved. This is due to the fact that communications among CPUs becomes dominating, the application of the preconditioner being the main responsible as highlighted by $\overline{T}^P$.

^§Performing the strong scalability test with a number of CPUs smaller than 800 was not feasible with this mesh; indeed, memory limitations occurred due the excessive number of unknowns per single CPU.

6.   Conclusion

In this work we proposed a novel monolithic solver for the simulation of the integrated electromechanics problem for the LV of the human heart. We coupled the monodomain equation, the ionic minimal model, the activation model for the fibers contraction, and the myocardial mechanics in the active strain framework with a transmurally variable activation. The interaction of the myocardium with the blood inside the ventricle was taken into account by prestressing and by solving additional 0D problems for the fluid.

We considered both implicit and semi-implicit BDF schemes for our monolithic solver together with a preconditioning strategy based on block Gauss-Seidel preconditioner which exploits the multiphysics structure of the coupled problem. The strong limitations on the timestep size in the semi-implicit case are such that the implicit scheme is by far the most efficient for numerical simulations lasting for one or multiple heartbeats. We solve the coupled problem in the High Performance Computing framework and we performed a scalability analysis for a benchmark test involving more than 37'000'000 unknowns.

The results obtained by simulating a complete heartbeat with both ideal and patient-specific geometries are qualitatively in agreement with physiological values measured in healthy patients. Moreover, the comparison of the simulations on the idealized geometry with different mesh sizes shows that the proposed numerical model is able to capture the physical solution with satisfactory accuracy even in the case of relatively coarse meshes. From our simulations we observe that choosing a mesh with maximum size $h$ of about 2 mm is sufficient to correctly reproduce the contraction of the LV, while negligible differences to the corresponding pV-loops are observable when using finer meshes.

Acknowledgments

This research was partially supported by the Swiss Platform for Advanced Scientific Computing (PASC, project "Integrative HPC Framework for Coupled Cardiac Simulations"). We also gratefully acknowledge the Swiss National Supercomputing Center (CSCS) for providing the CPU resources for the numerical simulations under project ID s635. Finally, we acknowledge the support from the ERC Advanced Grant iHEART, “An Integrated Heart Model for the simulation of the cardiac function”, 2017-2022, P.I.A. Quarteroni (ERC-2016-ADG, project ID: 740132).

We acknowledge Prof. J. Schwitter and Dr. P. Masci (CHUV) for providing the MRI images used in this work and for the constructive discussions on the topic from a physician's perspective. We also thank Prof. P. Tozzi (CHUV) for the invaluable insights in the functioning of the human heart. Finally, we acknowledge both Dr. D. Forti and MER S. Deparis (EPFL) for the useful discussions and their precious help with the LifeV library.

Conflict of interest

The authors declare no conflict of interest.