Electro-Mechanical Whole-Heart Digital Twins: A Fully Coupled Multi-Physics Approach

Abstract

Mathematical models of the human heart are evolving to become a cornerstone of precision medicine and support clinical decision making by providing a powerful tool to understand the mechanisms underlying pathophysiological conditions. In this study, we present a detailed mathematical description of a fully coupled multi-scale model of the human heart, including electrophysiology, mechanics, and a closed-loop model of circulation. State-of-the-art models based on human physiology are used to describe membrane kinetics, excitation-contraction coupling and active tension generation in the atria and the ventricles. Furthermore, we highlight ways to adapt this framework to patient specific measurements to build digital twins. The validity of the model is demonstrated through simulations on a personalized whole heart geometry based on magnetic resonance imaging data of a healthy volunteer. Additionally, the fully coupled model was employed to evaluate the effects of a typical atrial ablation scar on the cardiovascular system. With this work, we provide an adaptable multi-scale model that allows a comprehensive personalization from ion channels to the organ level enabling digital twin modeling.

1. Introduction

Cardiovascular diseases are the biggest contributors to the mortality and morbidity in the European Union, affecting millions of people every year [1]. While diagnostic tools and therapeutic options continuously improve and more data become available to researchers, the treatment of diseases, such as ischemic heart disease or atrial fibrillation remains difficult due to the highly complex nature of the human heart and cardiovascular system. Evidently, the complex mechanisms underlying these pathophysiological conditions are notoriously difficult to evaluate in the clinical environment due to ethical and technical limitations. Computational models of the human heart have the ability to avoid these limitations and provide a valuable tool for clinical research and practice [2]. Already today, these models can improve diagnosis [3], risk stratification [4], therapy planning [5,6], and intraprocedural support [7]. Furthermore, the vision of providing therapies customized to each individual patient heavily relies on mechanistic models to build a digital twin based on our knowledge of physiology and the fundamental laws of physics on the one hand and measured characteristics of the individual patient on the other hand [8]. An essential component of such a mechanistic model is the representation of the cardiac electrophysiology (EP) to reproduce the depolarization and repolarization sequence with a reaction-diffusion model. On the microscopic or cellular level, the reaction part was first described by Hodgkin and Huxley [9], while the diffusion on the macroscopic level can be modeled as described in [10]. Another essential component of the mechanistic model is cardiac mechanics (M) to describe the deformation, and the interplay between these systems described by myofilament models to drive the contraction and relaxation of the myocardium. The latter depends on the loading conditions imposed by the circulatory system, which is most accurately described by a fluid-structure interaction (FSI) problem, and the tissue that is surrounding the heart. Due to the complexity of this multi-physics problem, only a few cardiac simulation frameworks have been proposed that solve the complete EP-M-FSI system [11,12,13,14]. The most common way is to replace the FSI part by lumped parameter models of the circulatory system. More specifically, most simulations of the coupled EP-M problem focus on (bi-)ventricular models of a single heart beat and therefore only require isolated ventricular pre- and afterload models such as a 3- or 4-element Windkessel. To observe changes in cardiac function over multiple heartbeats, closed-loop models [15,16,17,18] are necessary and need to be coupled to the mechanical system [19,20,21,22,23,24,25]. EP-M models of the whole heart have been present for a few years now and have made impressive progress (e.g., [26,27,28,29,30]). However, only a few included the active contraction of the atria [31,32,33]. Land et al. [34] developed a model of active contraction for the atria based on human tissue preparations. This model can now be incorporated into coupled tissue level simulations, as was previously done for the ventricles [35,36].

Another aspect of such a complex model is the parameterization of all its components. With personalized heart models in mind, this is ideally done using patient-specific measurements. However, the majority of parameters can only be measured by invasive procedures, estimated indirectly, or not at all due technical limitations. Hence, efficient workflows are necessary to gain as much information as possible and incorporate them into digital twin models. Anatomically accurate heart models based on imaging data are created using (semi-)automatic workflows [37,38,39]. To make these heart models comparable to each other, universal coordinate systems have been proposed for the atria [40] and the ventricles [41,42]. Furthermore, we can build on already existing pipelines for the personalization of passive mechanical behavior of the heart [43,44], ventricular afterload models [45], and cardiac electrophysiology based on electrocardiograms [46,47,48,49] or electroanatomical mapping [50,51].

In this study, we present a detailed mathematical description of a fully coupled multi-scale model of the human heart building on previous work in cardiac mechanics [25,31,43,52,53,54,55] and cardiac electrophysiology [49,56,57,58,59,60]. We parameterize and use state-of-the-art myofilament models based on human physiology to describe membrane kinetics, excitation-contraction coupling and active tension generation in the atria and the ventricles. To solve the coupled EP-M problem, we apply a staggered scheme in time where the monodomain equation and the deformation are solved sequentially. Additionally, the proposed electro-mechanical whole-heart model is coupled to a novel 0D closed-loop model of the cardiovascular system. Here, we use a quasi Newton method to update the pressure values in all four chambers and reach convergence in fewer iterations compared to standard Newton methods. Furthermore, we briefly touch on the possibilities of such models to be tailored to patient-specific measurements and demonstrate their capabilities by simulating the complete heart cycle using an anatomical model of the heart of a healthy volunteer. Subsequently, the simulation data are evaluated and compared to magnetic resonance imaging (MRI) data with a focus on left ventricular function and atrio-ventricular plane displacement (AVPD). Finally, a typical radio-frequency ablation (RFA) scar is introduced into the left atrium and simulation results are compared to the healthy reference case. Hence, we demonstrate not only the possibility to adapt the model to pathological scenarios but the capability of the whole-heart model to respond to local changes in loading conditions.

2. Methods

2.1. Four-Chamber Heart Model

The geometric model used in this study is based on magnetic resonance imaging (MRI) data of a 32-year-old healthy volunteer provided by Heidelberg University Hospital. The data were acquired using a 1.5 T MR tomography system (Philips Medical Systems) and consist of a static whole heart image at diastasis, as well as time resolved MRI in 2-, 3-, and 4-chamber long axis view and 12 time resolved short axis slices with a 10 mm spacing. Endocardial and epicardial boundaries of the myocardial wall were segmented and labeled manually from the static whole heart image to create the atrial (${\mathsf{\Omega }}_{\mathrm{A}}$) and ventricular (${\mathsf{\Omega }}_{\mathrm{V}}$) domains. As in Fritz et al. [31], the heart is surrounded by a pericardial layer (${\mathsf{\Omega }}_{\mathrm{Pericardium}}$) representing the surrounding tissue in which the heart is embedded. Additionally, the tetrahedral mesh was extended by volumetric representations of the inferior and superior vena cava, the pulmonary veins, pulmonary artery, ascending aorta (further summarized in ${\mathsf{\Omega }}_{\mathrm{MajorVessels}}$), the mitral and tricuspidal valve (${\mathsf{\Omega }}_{\mathrm{Valves}}$), and adipose tissue (${\mathsf{\Omega }}_{\mathrm{AdiposeTissue}}$) around the base of the heart. Previously, the openings of the pulmonary veins and the venae cavae were fixed directly, which constrained atrial motion. This limitation was ameliorated by instead fixing the terminal ends of the major arteries and veins to allow for a free contraction of the atria [61]. Adipose tissue was added to the free space between the epicardial surface and the pericardium to obtain a continuous contact surface [32]. Two different discretizations are used: (1) the mechanical reference domain ${\mathsf{\Omega }}_{\mathrm{M}}={\mathsf{\Omega }}_{\mathrm{V}}\cup {\mathsf{\Omega }}_{\mathrm{A}}\cup {\mathsf{\Omega }}_{\mathrm{Valves}}\cup {\mathsf{\Omega }}_{\mathrm{AdiposeTissue}}\cup {\mathsf{\Omega }}_{\mathrm{MajorVessels}}\cup {\mathsf{\Omega }}_{\mathrm{Pericardium}}$, as shown in Figure 1b; (2) the electrophysiological reference domain ${\mathsf{\Omega }}_{\mathrm{EP}}={\mathsf{\Omega }}_{\mathrm{V}}\cup {\mathsf{\Omega }}_{\mathrm{A}}$ as a subset of ${\mathsf{\Omega }}_{\mathrm{M}}$.

2.2. Cardiac Elasto-Mechanics

We consider the bounded reference domain ${\mathsf{\Omega }}_{\mathrm{M}}\subset {\mathbb{R}}^3$ for the human heart, which is deformed at time $t\in [0,T]$ into the current configuration ${\mathsf{\Omega }}_{\mathrm{M}}^{\mathit{\phi }}=\mathit{\phi }(t,{\mathsf{\Omega }}_{\mathrm{M}})$ by the deformation vector

$$\mathit{\phi }:[0,T]\times {\mathsf{\Omega }}_{\mathrm{M}}⟶{\mathbb{R}}^3,\mathit{\phi }(t,{\displaystyle \mathbf{x}})={\displaystyle \mathbf{x}}+{\displaystyle \mathbf{u}}(t,{\displaystyle \mathbf{x}}),$$

where $\mathbf{u}:[0,T]\times {\mathsf{\Omega }}_{\mathrm{M}}⟶{\mathbb{R}}^3$ is the displacement at a given position $\mathbf{x}\in {\mathsf{\Omega }}_{\mathrm{M}}$ over time. The deformation on the continuum body ${\mathsf{\Omega }}_{\mathrm{M}}$ is characterized by the axioms of conservation of mass, linear momentum and angular momentum. Following ([62], Chapter 3), the displacement $\mathbf{u}$ is determined by the solution of the Navier–Cauchy equations for $t\in (0,T)$

$$\begin{array}{cc}\hfill {\varrho }_0{{\displaystyle \mathbf{a}}}^{\mathit{\phi }}-{div}^{\mathit{\phi }}{{\displaystyle \mathbf{T}}}^{\mathit{\phi }}& ={{\displaystyle \mathbf{f}}}^{\mathit{\phi }}\mathrm{in}{\mathsf{\Omega }}_{\mathrm{M}}^{\mathit{\phi }},\hfill \end{array}$$

(1a)

$$\begin{array}{cc}\hfill {{\displaystyle \mathbf{T}}}^{\mathit{\phi }}{{\displaystyle \mathbf{n}}}^{\mathit{\phi }}& ={{\displaystyle \mathbf{g}}}^{\mathit{\phi }}\mathrm{on}{\mathsf{\Gamma }}_{\mathrm{N}}^{\mathit{\phi }},\hfill \end{array}$$

(1b)

with reference density ${\varrho }_0$, acceleration ${\mathbf{a}}^{\mathit{\phi }}$ and the Cauchy stress ${\mathbf{T}}^{\mathit{\phi }}$. The vector-valued function ${\mathbf{g}}^{\mathit{\phi }}$ represents applied boundary forces on all Neumann boundaries ${\mathsf{\Gamma }}_{\mathrm{N}}^{\mathit{\phi }}$. The PDE (1) describes the deformation on the unknown domain ${\mathsf{\Omega }}_{\mathrm{M}}^{\mathit{\phi }}$. We therefore pull back the displacement onto the reference domain using the techniques provided in [62,63]: Denoting the Fréchet derivative in space by $\mathrm{D}$, we introduce the deformation gradient

$${\displaystyle \mathbf{F}}:=\mathrm{D}\mathit{\phi }:[0,T]\times {\mathsf{\Omega }}_{\mathrm{M}}⟶{\mathbb{R}}^{3\times 3},{\displaystyle \mathbf{F}}(t,{\displaystyle \mathbf{x}})={\displaystyle \mathbf{I}}+\mathrm{D}{\displaystyle \mathbf{u}}(t,{\displaystyle \mathbf{x}}),$$

and define the symmetric second Piola–Kirchhoff stress tensor

$${\displaystyle \mathbf{S}}({\displaystyle \mathbf{x}},{\displaystyle \mathbf{F}}):=det\left({\displaystyle \mathbf{F}}\right){{\displaystyle \mathbf{F}}}^{-1}{{\displaystyle \mathbf{T}}}^{\mathit{\phi }}{{\displaystyle \mathbf{F}}}^{-\top },$$

to obtain the following equations in the reference configuration for $t\in (0,T)$

$$\begin{array}{cc}\hfill {\varrho }_0{\partial }_t^2{\displaystyle \mathbf{u}}-div\left({\displaystyle \mathbf{F}}{\displaystyle \mathbf{S}}\right)& ={\displaystyle \mathbf{f}}\mathrm{in}{\mathsf{\Omega }}_{\mathrm{M}},\hfill \\ \hfill {\displaystyle \mathbf{F}}{\displaystyle \mathbf{S}}{\displaystyle \mathbf{n}}& ={\displaystyle \mathbf{g}}\mathrm{on}{\mathsf{\Gamma }}_N.\hfill \end{array}$$

The passive material properties of the cardiac tissue are described as follows: we assume that the material is hyperelastic, i.e., a stored energy function $W=W\left(\mathbf{F}\right)$ exists, such that the stress response is given by

$${\displaystyle \mathbf{S}}\left({\displaystyle \mathbf{F}}\right)={{\displaystyle \mathbf{F}}}^{-1}\frac{\mathrm{D}W\left(\mathbf{F}\right)}{\mathrm{D}\mathbf{F}}.$$

As a consequence of objectivity, the material is frame-indifferent, i.e., there exist representations $\tilde{W}=\tilde{W}\left(\mathbf{C}\right)$ and $\hat{W}=\hat{W}\left(\mathbf{E}\right)$, such that

$$\begin{array}{cccccc}\hfill W\left({\displaystyle \mathbf{F}}\right)& =\tilde{W}\left({\displaystyle \mathbf{C}}\right)=\hat{W}\left({\displaystyle \mathbf{E}}\right),\hfill & \hfill {{\displaystyle \mathbf{F}}}^{-1}\frac{\mathrm{D}W\left(\mathbf{F}\right)}{\mathrm{D}\mathbf{F}}& =2\frac{\mathrm{D}\tilde{W}\left(\mathbf{C}\right)}{\mathrm{D}\mathbf{C}}=\frac{\mathrm{D}\hat{W}\left(\mathbf{E}\right)}{\mathrm{D}\mathbf{E}},\hfill & \hfill {\displaystyle \mathbf{C}}=& {{\displaystyle \mathbf{F}}}^{\top }{\displaystyle \mathbf{F}},{\displaystyle \mathbf{E}}=\frac{1}{2}({{\displaystyle \mathbf{F}}}^{\top }{\displaystyle \mathbf{F}}-{\displaystyle \mathbf{I}}).\hfill \end{array}$$

The model is complemented by additional boundary conditions (Figure 1a), namely a Dirichlet boundary on ${\mathsf{\Gamma }}_{\mathrm{D}}$ on the distal end of the major vessels ${\mathsf{\Omega }}_{\mathrm{MajorVessels}}$ and the outer surface of the pericardium ${\mathsf{\Omega }}_{\mathrm{Pericardium}}$. On the inner surface of the pericardium ${\mathsf{\Gamma }}_{\mathrm{P}}$, a contact boundary condition for the deformation is used (see Section 2.2.1). The Neumann boundary ${\mathsf{\Gamma }}_{\mathrm{N}}$ is given by the inner surfaces of the cardiac chambers ${\mathsf{\Gamma }}_{\mathrm{LV}}\cup {\mathsf{\Gamma }}_{\mathrm{RV}}\cup {\mathsf{\Gamma }}_{\mathrm{LA}}\cup {\mathsf{\Gamma }}_{\mathrm{RA}}$, where a pressure

$${\displaystyle \mathbf{p}}=\left(p_{\mathrm{LV}},p_{\mathrm{RV}},p_{\mathrm{LA}},p_{\mathrm{RA}}\right)\in {\mathbb{R}}^4$$

(2)

from the cardiovascular system (see Section 2.2.2) is applied.

With these assumptions, we obtain the equations of cardiac elasticity in the reference configuration for $t\in (0,T)$

$$\begin{array}{cccc}\hfill {\varrho }_0{\partial }_t^2{\displaystyle \mathbf{u}}-div\left({\displaystyle \mathbf{F}}{\displaystyle \mathbf{S}}({\displaystyle \mathbf{x}},{\displaystyle \mathbf{F}})\right)& ={\displaystyle \mathbf{0}}\hfill & & \mathrm{in}{\mathsf{\Omega }}_{\mathrm{M}},\hfill \end{array}$$

(3a)

$$\begin{array}{cccc}\hfill {\displaystyle \mathbf{u}}& ={\displaystyle \mathbf{0}}\hfill & & \mathrm{on}{\mathsf{\Gamma }}_{\mathrm{D}},\hfill \end{array}$$

(3b)

$$\begin{array}{cccc}\hfill \left[\mathit{\phi }\right]& ={\displaystyle \mathbf{0}}\hfill & & \mathrm{on}{\mathsf{\Gamma }}_{\mathrm{P}},\hfill \end{array}$$

(3c)

$$\begin{array}{cccc}\hfill {\displaystyle \mathbf{F}}{\displaystyle \mathbf{S}}({\displaystyle \mathbf{x}},{\displaystyle \mathbf{F}}){\displaystyle \mathbf{n}}& =p_C{\displaystyle \mathbf{n}}\hfill & & \mathrm{on}{\mathsf{\Gamma }}_C,C\in \left\{\mathrm{LV},\mathrm{RV},\mathrm{LA},\mathrm{RA}\right\}.\hfill \end{array}$$

(3d)

The cardiac muscle tissue is nearly incompressible, i.e., $J:=det\left(\mathbf{F}\right)\approx 1$. Therefore, we include, in the elastic energy, a penalty function $W_{\kappa }$ with bulk modulus $\kappa \gg 1$, which satisfies the conditions

$$W_{\kappa }\left(1\right)=0,\underset{J\to \infty }{lim}{W}_{\kappa }\left(J\right)=\infty ,\underset{J\to 0}{lim}{W}_{\kappa }\left(J\right)=\infty .$$

For our applications, we use two different hyperelastic materials in subsets of ${\mathsf{\Omega }}_{\mathrm{M}}$. For the active part of the heart composed of cardiomyocytes, we use the anisotropic model of Guccione et al. [64]:

$$W_{\mathrm{G}}={\widehat{W}}_{\mathrm{G}}\left({\displaystyle \mathbf{E}}\right)=\frac{{\mu }_{\mathrm{G}}}{2}exp\left[Q\left(\overline{{\displaystyle \mathbf{E}}}\right)-1\right]+\frac{\kappa }{2}{(J-1)}^2,\overline{{\displaystyle \mathbf{E}}}=\frac{1}{2}({\displaystyle \mathbf{Q}}{{\displaystyle \mathbf{F}}}^{\top }{\displaystyle \mathbf{F}}{{\displaystyle \mathbf{Q}}}^{\top }-{\displaystyle \mathbf{I}}),$$

(4)

where ${\mu }_{\mathrm{G}}>0$ is the shear modulus, $W_{\mathrm{G},\kappa }\left(J\right)={(J-1)}^2$, and $\mathbf{Q}=(\mathbf{f}\mid \mathbf{s}\mid \mathbf{n})$ is the transformation resulting from the orthogonal fiber, sheet and sheet-normal directions

$$\begin{array}{ccccccc}\hfill {\displaystyle \mathbf{f}}& :{\mathsf{\Omega }}_{\mathrm{M}}⟶S^2,\hfill & \hfill {\displaystyle \mathbf{s}}& :{\mathsf{\Omega }}_{\mathrm{M}}⟶S^2,\hfill & \hfill {\displaystyle \mathbf{n}}& :{\mathsf{\Omega }}_{\mathrm{M}}⟶S^2,\hfill & \hfill \end{array}$$

(5)

and $Q:{\mathbb{R}}^{3\times 3}\to \mathbb{R}$ is a scalar function

$$Q\left(\overline{\mathbf{E}}\right)=b_{\mathbf{f}}{E}_{11}^2+b_{\mathbf{s}}(E_{22}^2+E_{33}^2+E_{23}^2+E_{32}^2)+b_{\mathbf{f}\mathbf{s}}(E_{12}^2+E_{21}^2+E_{13}^2+E_{31}^2),\overline{\mathbf{E}}={\left(E_{ij}\right)}_{i,j=1,\dots ,3},$$

with scaling parameters $b_{\mathbf{f}}$, $b_{\mathbf{s}}$, $b_{\mathbf{f}\mathbf{s}}$. Purely passive tissue is characterized by a Neo–Hookean material

$$W_{\mathrm{NH}}={\tilde{W}}_{\mathrm{NH}}\left({\displaystyle \mathbf{C}}\right)={\mu }_{\mathrm{NH}}(tr\left({\displaystyle \mathbf{C}}\right)-3)-{\mu }_{\mathrm{NH}}ln(det{\displaystyle \mathbf{C}})+\frac{\kappa }{2}{ln}^2(det{\displaystyle \mathbf{C}}),$$

with shear modulus ${\mu }_{\mathrm{NH}}>0$ and $W_{\mathrm{NH},\kappa }=\frac{1}{2}{ln}^2\left(J\right)$.

In addition to the passive material behavior, we model the active material response to electrical activation by an active stress approach. Given the deformed fiber directions ${\mathbf{l}}^{\mathit{\phi }}=\mathbf{F}\mathbf{l}$, $\mathbf{l}=\mathbf{f},\mathbf{s},\mathbf{n}$ and the active tension $T_{\mathrm{tot}}$ generated by a cellular tension evolution model depending on $(\mathbf{c},\mathbf{q})$ and determined by (10c) and (12), we calculate the active stress response by

$${{\displaystyle \mathbf{T}}}_{\mathrm{a}}=T_{\mathrm{tot}}({\displaystyle \mathbf{c}},{\displaystyle \mathbf{q}})\frac{{\mathbf{f}}^{\mathit{\phi }}\otimes {\mathbf{f}}^{\mathit{\phi }}}{\parallel {\mathbf{f}}^{\mathit{\phi }}\parallel }.$$

(6)

We use the additive decomposition $\mathbf{S}={\mathbf{F}}^{-1}{\mathrm{D}}_{\mathbf{F}}W\left(\mathbf{F}\right)+J{\mathbf{F}}^{-1}{\mathbf{T}}_{\mathrm{a}}{\mathbf{F}}^{-\top }$ to compute the stress.

2.2.1. Contact Boundary Conditions

Cardiac deformation is physiologically restricted by the pericardium and the surrounding tissue. In particular, it has been shown that including the interaction between ventricles, atria, and the pericardium during cardiac simulations allows one to better reproduce atrioventricular plane displacement (AVPD) [31,32] and thus cardiac motion. Realistic implementations of the interaction between the heart and surrounding tissue were first introduced by Fritz et al. [31]. They assumed frictionless and permanent contact between the epicardium and the pericardium using a penalty formulation. Other models describe these interactions using normal Robin boundary conditions on the epicardium with constant spring stiffness [26,32] or spatially varying spring stiffness [61] informed by image-derived motion. In this study, the sliding contact problem in Fritz et al. [31] between the epicardium and the inner pericardial surface (3c) is used by defining a contact energy

$$W_{\mathrm{P}}=\epsilon {\int }_{{\mathsf{\Gamma }}_{\mathrm{P}}}{g}_{\mathrm{P}}^2\mathrm{d}{\displaystyle \mathbf{a}},$$

with the penalty parameter $\epsilon >0$ and the penalty functional

$$g_{\mathrm{P}}(t,{\displaystyle \mathbf{x}})=({{\displaystyle \mathbf{x}}}^{\perp }-\mathit{\phi }(t,{\displaystyle \mathbf{x}}))·{\displaystyle \mathbf{n}},$$

where $\mathbf{n}$ is the surface normal of the epicardium and ${\mathbf{x}}^{\perp }$ the projection of the deformed point onto the pericardium. This contact energy is then added to the hyperelastic energy functional.

2.2.2. Closed-Loop Circulatory Model

The heart’s main purpose is to effectively pump blood throughout the body. Therefore, it is important to define realistic hemodynamic boundary conditions for a computer model of cardiac mechanics. Most accurately, this is done by solving an FSI problem. To date, FSI models were typically restricted to a single chamber due to their computational complexity [11,12,13,14]. Since here, we are not interested in a detailed distribution of pressure and flow throughout all cardiac chambers, an alternative solution is to use a lumped parameter model of the circulatory system. For simulations spanning multiple heart beats, a closed-loop model is required to preserve the total blood volume in the circulatory system. Several closed-loop models of the circulatory system have been proposed [15,18,23,24,65] and all of them share some of their structure. However, in all cases, the circulatory system was only coupled to the finite element model of one or both of the ventricles, while the atria were represented by a time-varying elastance. To the best of our knowledge, there is no published work on coupling all four heart chambers with the circulatory model.

In the deformed configuration, the ventricular volume is computed from the reference domains for each chamber $\mathrm{C}$ and the corresponding deformed chamber volumes are then denoted by $|{\mathsf{\Omega }}_{\mathrm{LV}}^{\mathit{\phi }}|,|{\mathsf{\Omega }}_{\mathrm{RV}}^{\mathit{\phi }}|,|{\mathsf{\Omega }}_{\mathrm{LA}}^{\mathit{\phi }}|,|{\mathsf{\Omega }}_{\mathrm{RA}}^{\mathit{\phi }}|$ with

$$|{\mathsf{\Omega }}_{\mathrm{C}}^{\mathit{\phi }}|={\int }_{{\mathsf{\Omega }}_{\mathrm{C}}^{\mathit{\phi }}}1\mathrm{d}{\displaystyle \mathbf{x}}={\int }_{{\mathsf{\Omega }}_{\mathrm{C}}}J\mathrm{d}{\displaystyle \mathbf{x}},\mathrm{C}\in \left\{\mathrm{LV},\mathrm{RV},\mathrm{LA},\mathrm{RA}\right\}.$$

We approximate the chamber volumes on the discrete geometry by summing up the volumes of all tetrahedrons constructed by a surface triangle $K\subset {\mathsf{\Gamma }}_{\mathrm{C}}$ and an inner point ${\mathbf{x}}_{\mathrm{C}}$. We assume that each chamber ${\mathsf{\Omega }}_{\mathrm{C}}^{\mathit{\phi }}$ is star-shaped around the point ${\mathbf{x}}_{\mathrm{C}}$. By denoting the vectors of the vertices of K by ${\mathbf{a}}_K$, ${\mathbf{b}}_K$, ${\mathbf{c}}_K$, we obtain

$$|{\mathsf{\Omega }}_{\mathrm{C}}^{\mathit{\phi }}| \approx \sum _{K\subset \partial {\mathsf{\Omega }}_{\mathrm{C}}}\frac{1}{6}(\mathit{\phi }(t,{{\displaystyle \mathbf{a}}}_K)-\mathit{\phi }(t,{{\displaystyle \mathbf{x}}}_{\mathrm{C}}))·\left((\mathit{\phi }(t,{{\displaystyle \mathbf{b}}}_K)-\mathit{\phi }(t,{{\displaystyle \mathbf{a}}}_K))\times (\mathit{\phi }(t,{{\displaystyle \mathbf{c}}}_K)-\mathit{\phi }(t,{{\displaystyle \mathbf{a}}}_K))\right).$$

To simulate the response of the circulatory system given the discrete volumes, we reinterpret the system as a series of transmissions [66] as shown in Figure 2. Our goal is to find a pressure vector (2) in such a way, that the chamber volumes of the circulatory model coincide with the discrete volumes of the mechanics simulation [20]. The circulatory chamber and vessel volumes $\mathbf{v}:[0,T]\to {\mathbb{R}}^8$ are denoted by

$${\displaystyle \mathbf{v}}=(v_{\mathrm{LV}},v_{\mathrm{RV}},v_{\mathrm{LA}},v_{\mathrm{RA}},v_{\mathrm{SysVen}},v_{\mathrm{SysArt}},v_{\mathrm{PulVen}},v_{\mathrm{PulArt}}).$$

A series of ODEs is solved to find a solution $\mathbf{v}$, to obtain

$$v_{\mathrm{LV}}=|{\mathsf{\Omega }}_{\mathrm{LV}}^{\mathit{\phi }}|,v_{\mathrm{RV}}=|{\mathsf{\Omega }}_{\mathrm{RV}}^{\mathit{\phi }}|,v_{\mathrm{LA}}=|{\mathsf{\Omega }}_{\mathrm{LA}}^{\mathit{\phi }}|,v_{\mathrm{RA}}=\left|{\mathsf{\Omega }}_{\mathrm{RA}}^{\mathit{\phi }}\right|.$$

Let $\mathbf{z}:[0,T]\to {\mathbb{R}}^{12}$ be the intermediate variables

$${\displaystyle \mathbf{z}}=\left(p_{\mathrm{SysVen}},p_{\mathrm{SysArt}},p_{\mathrm{PulVen}},p_{\mathrm{PulArt}},Q_{\mathrm{SysArt}},Q_{\mathrm{SysPer}},Q_{\mathrm{SysVen}},Q_{\mathrm{Rav}}{Q}_{\mathrm{PulArt}},Q_{\mathrm{PulPer}},Q_{\mathrm{PulVen}},Q_{\mathrm{Lav}}\right),$$

describing the systemic venous, systemic arterial, pulmonary venous and pulmonary arterial pressures and flows. Then the closed-loop equations read

$$\begin{array}{cc}\hfill {\partial }_t{\displaystyle \mathbf{v}}& ={{\displaystyle \mathbf{G}}}_{\mathbf{v}}({\displaystyle \mathbf{p}},{\displaystyle \mathbf{z}})\mathrm{in}(0,T),\hfill \end{array}$$

(7a)

depending on

$$\begin{array}{cc}\hfill {\displaystyle \mathbf{z}}& ={{\displaystyle \mathbf{G}}}_{\mathbf{z}}({\displaystyle \mathbf{p}},{\displaystyle \mathbf{v}}).\hfill \end{array}$$

(7b)

The relations $\mathbf{z}={\mathbf{G}}_{\mathbf{z}}(\mathbf{p},\mathbf{v})$ in (7b) are given by the flows

$$\begin{array}{cccc}\hfill Q_{\mathrm{SysArt}}& =max\left\{\frac{p_{\mathrm{LV}}-\frac{v_{\mathrm{SysArt}}}{C_{\mathrm{SysArt}}}}{R_{\mathrm{SysArtValve}}+R_{\mathrm{SysArt}}},0\right\},\hfill & \hfill Q_{\mathrm{SysPer}}& =\frac{\frac{v_{\mathrm{SysArt}}}{C_{\mathrm{SysArt}}}-p_{\mathrm{SysVen}}}{R_{\mathrm{SysPer}}},\hfill \\ \hfill Q_{\mathrm{SysVen}}& =\frac{p_{\mathrm{SysVen}}-p_{\mathrm{RA}}}{R_{\mathrm{SysVen}}},\hfill & \hfill Q_{\mathrm{Rav}}& =max\left\{\frac{p_{\mathrm{RA}}-p_{\mathrm{RV}}}{R_{\mathrm{RavValve}}},0\right\},\hfill \\ \hfill Q_{\mathrm{PulArt}}& =max\left\{\frac{p_{\mathrm{RV}}-\frac{v_{\mathrm{PulArt}}}{C_{\mathrm{PulArt}}}}{R_{\mathrm{PulArtValve}}+R_{\mathrm{PulArt}}},0\right\},\hfill & \hfill Q_{\mathrm{PulPer}}& =\frac{\frac{v_{\mathrm{PulArt}}}{C_{\mathrm{PulArt}}}-p_{\mathrm{PulVen}}}{R_{\mathrm{PulPer}}},\hfill \\ \hfill Q_{\mathrm{PulVen}}& =\frac{p_{\mathrm{PulVen}}-p_{\mathrm{LA}}}{R_{\mathrm{PulVen}}},\hfill & \hfill Q_{\mathrm{Lav}}& =max\left\{\frac{p_{\mathrm{LA}}-p_{\mathrm{LV}}}{R_{\mathrm{LavValve}}},0\right\},\hfill \end{array}$$

and the circulatory pressures

$$\begin{array}{cccc}\hfill p_{\mathrm{SysArt}}& =p_{\mathrm{LV}}-R_{\mathrm{SysArtValve}}{Q}_{\mathrm{SysArt}},\hfill & \hfill p_{\mathrm{SysVen}}& =\frac{v_{\mathrm{SysVen}}}{C_{\mathrm{SysVen}}},\hfill \\ \hfill p_{\mathrm{PulArt}}& =p_{\mathrm{RV}}-R_{\mathrm{PulArtValve}}{Q}_{\mathrm{PulArt}},\hfill & \hfill p_{\mathrm{PulVen}}& =\frac{v_{\mathrm{PulVen}}}{C_{\mathrm{PulVen}}}.\hfill \end{array}$$

The evolution of $\mathbf{v}$ in (7a) is determined by

$$\begin{array}{cccc}\hfill {\partial }_t{v}_{\mathrm{RV}}& =Q_{\mathrm{Rav}}-Q_{\mathrm{PulArt}},\hfill & \hfill {\partial }_t{v}_{\mathrm{SysArt}}& =Q_{\mathrm{SysArt}}-Q_{\mathrm{SysPer}},\hfill \\ \hfill {\partial }_t{v}_{\mathrm{LV}}& =Q_{\mathrm{Lav}}-Q_{\mathrm{SysArt}},\hfill & \hfill {\partial }_t{v}_{\mathrm{SysVen}}& =Q_{\mathrm{SysPer}}-Q_{\mathrm{SysVen}},\hfill \\ \hfill {\partial }_t{v}_{\mathrm{RA}}& =Q_{\mathrm{SysVen}}-Q_{\mathrm{Rav}},\hfill & \hfill {\partial }_t{v}_{\mathrm{PulArt}}& =Q_{\mathrm{PulArt}}-Q_{\mathrm{PulPer}},\hfill \\ \hfill {\partial }_t{v}_{\mathrm{LA}}& =Q_{\mathrm{PulVen}}-Q_{\mathrm{Lav}},\hfill & \hfill {\partial }_t{v}_{\mathrm{PulVen}}& =Q_{\mathrm{PulPer}}-Q_{\mathrm{PulVen}}.\hfill \end{array}$$

The pressure values for the boundary condition (3d) are computed iteratively at each time step. Within one time step $t_n$ of the mechanical problem, we first find an approximated pressure ${\mathbf{p}}^{n,0}$ through the closed-loop model. During the first time steps, we choose ${\mathbf{p}}^{n,0}$ as a fixed increment of $1\mathrm{Pa}$ from the previous time step. After the fifth time step, we extrapolate ${\mathbf{p}}^{n,0}$ using a 4th order Adams-Bashforth scheme:

$${{\displaystyle \mathbf{p}}}^{n,0}={{\displaystyle \mathbf{p}}}^{n-1}+\Delta t_n\left(\frac{55}{24}\frac{\Delta {\mathbf{p}}^{n-1}}{\Delta t_{n-1}}-\frac{59}{24}\frac{\Delta {\mathbf{p}}^{n-2}}{\Delta t_{n-2}}+\frac{37}{24}\frac{\Delta {\mathbf{p}}^{n-3}}{\Delta t_{n-3}}-\frac{9}{24}\frac{\Delta {\mathbf{p}}^{n-4}}{\Delta t_{n-4}}\right),$$

(8)

with $\frac{\Delta {\mathbf{p}}^k}{\Delta t_k}=\frac{{\mathbf{p}}^k-{\mathbf{p}}^{k-1}}{t_k-t_{k-1}}$, $k\in \mathbb{N}$. We then start the iteration to calculate the correct pressure: during each iteration cycle i, we update the mechanical and circulatory models for the given pressure boundary ${\mathbf{p}}^{n,i}$. The contraction is calculated as in Section 2.5. The ODEs of the circulatory model are updated with a Runge–Kutta 4 scheme. We compare the residual $\mathbf{r}:{\mathbb{R}}^4\to {\mathbb{R}}^4$ for each chamber

$${\displaystyle \mathbf{r}}\left({{\displaystyle \mathbf{p}}}^{n,i}\right)={(r_{\mathrm{LV}},r_{\mathrm{RV}},r_{\mathrm{LA}},r_{\mathrm{RA}})}^{\top }$$

with $r_{\mathrm{C}}=\left|{\mathsf{\Omega }}_{\mathrm{C}}\right|-v_{\mathrm{C}}$ for $\mathrm{C}\in \left\{\mathrm{LV},\mathrm{RV},\mathrm{LA},\mathrm{RA}\right\}$. If the absolute values of the residuals $|r_{\mathrm{C}}|$ are below a threshold ${\epsilon }_{\mathbf{p}}$, the pressure ${\mathbf{p}}^{n,i}$ is accepted and we move to the next time step. Otherwise, we update $\mathbf{p}$ by a quasi-Newton method [67]:

$${{\displaystyle \mathbf{p}}}^{n,i}={{\displaystyle \mathbf{p}}}^{n,i-1}-{{\displaystyle \mathbf{C}}}_i^{-1}{\displaystyle \mathbf{r}}\left({{\displaystyle \mathbf{p}}}^{n,i}\right),$$

(9)

where ${\mathbf{C}}_i$ is the compliance matrix determined by

$${{\displaystyle \mathbf{C}}}_i^{-1}={{\displaystyle \mathbf{C}}}_{i-1}^{-1}+\left({\Delta {\displaystyle \mathbf{p}}}^{n,i}-{{\displaystyle \mathbf{C}}}_{i-1}^{-1}{{\displaystyle \Delta \mathbf{r}}}^{n,i}\right)\frac{{{\Delta \mathbf{p}}^{n,i}}^{\top }{\mathbf{C}}_{i-1}^{-1}}{{{\Delta \mathbf{p}}^{n,i}}^{\top }{\mathbf{C}}_{i-1}^{-1}{\Delta \mathbf{r}}^{n,i}},{{\displaystyle \mathbf{C}}}_0^{-1}={\displaystyle \mathbf{I}}$$

with ${\Delta \mathbf{p}}^{n,i}={\mathbf{p}}^{n,i}-{\mathbf{p}}^{n,i-1}$, ${\Delta \mathbf{r}}^{n,i}=\mathbf{r}\left({\mathbf{p}}^{n,i}\right)-\mathbf{r}\left({\mathbf{p}}^{n,i-1}\right)$. Updating the compliance matrix by applying a quasi-Newton method is an improvement compared to the modified Newton method proposed by Kerckhoffs et al. [20], which makes it possible to reach convergence in fewer iterations.

2.3. Cardiac Electrical Activity

On the reference domain ${\mathsf{\Omega }}_{\mathrm{EP}}$, the cardiac electric activity can be described by the monodomain model as explained by Keener and Sneyd [10]. This reaction-diffusion equation includes the transmembrane voltage $V_{\mathrm{m}}:[0,T]\times {\mathsf{\Omega }}_{\mathrm{EP}}\to \mathbb{R}$ coupled to ODEs describing the vector of gating variables $\mathbf{w}:[0,T]\times {\mathsf{\Omega }}_{\mathrm{EP}}\to {[0,1]}^{n_{\mathbf{w}}}$ and the vector of intracellular ion concentrations $\mathbf{c}:[0,T]\times {\mathsf{\Omega }}_{\mathrm{EP}}\to {\mathbb{R}}_{+}^{n_{\mathbf{c}}}$, so that the system in $(0,T)\times {\mathsf{\Omega }}_{\mathrm{EP}}$ is given by

$$\begin{array}{cc}\hfill \beta C_{\mathrm{m}}{\partial }_t{V}_{\mathrm{m}}& =\nabla ·\sigma \nabla V_{\mathrm{m}}-\beta I_{\mathrm{ion}}(V_{\mathrm{m}},{\displaystyle \mathbf{w}},{\displaystyle \mathbf{c}})+\beta I_{\mathrm{ext}},\hfill \end{array}$$

(10a)

$$\begin{array}{cc}\hfill {\partial }_t{\displaystyle \mathbf{w}}& ={{\displaystyle \mathbf{G}}}_{\mathbf{w}}(V_{\mathrm{m}},{\displaystyle \mathbf{w}},{\displaystyle \mathbf{c}}),\hfill \end{array}$$

(10b)

$$\begin{array}{cc}\hfill {\partial }_t{\displaystyle \mathbf{c}}& ={{\displaystyle \mathbf{G}}}_{\mathbf{c}}(V_{\mathrm{m}},{\displaystyle \mathbf{w}},{\displaystyle \mathbf{c}}),\hfill \end{array}$$

(10c)

depending on the conductivity tensor $\sigma :{\mathsf{\Omega }}_{\mathrm{EP}}\to {\mathbb{R}}^{3\times 3}$, the external stimulus $I_{\mathrm{ext}}:[0,T]\times {\mathsf{\Omega }}_{\mathrm{EP}}\to \mathbb{R}$, the total ionic current $I_{\mathrm{ion}}:\mathbb{R}\times {[0,1]}^{n_{\mathbf{w}}}\times {\mathbb{R}}_{+}^{n_{\mathbf{c}}}\to \mathbb{R}$, the functions ${\mathbf{G}}_{\mathbf{w}}:\mathbb{R}\times {[0,1]}^{n_{\mathbf{w}}}\times {\mathbb{R}}_{+}^{n_{\mathbf{c}}}\to {\mathbb{R}}^{n_{\mathbf{w}}}$ and ${\mathbf{G}}_{\mathbf{c}}:\mathbb{R}\times {[0,1]}^{n_{\mathbf{w}}}\times {\mathbb{R}}_{+}^{n_{\mathbf{c}}}\to {\mathbb{R}}^{n_{\mathbf{c}}}$, the cellular surface-to-volume ratio $\beta \in {\mathbb{R}}_{+}$ and the membrane capacitance $C_{\mathrm{m}}\in {\mathbb{R}}_{+}$. The system (10) is complemented by initial values at $t=0$ for all $\mathbf{x}\in {\mathsf{\Omega }}_{\mathrm{EP}}$

$$\begin{array}{c}\hfill V_{\mathrm{m}}(0,{\displaystyle \mathbf{x}})=V_{\mathrm{m}}^0\left({\displaystyle \mathbf{x}}\right),{\displaystyle \mathbf{w}}(0,{\displaystyle \mathbf{x}})={{\displaystyle \mathbf{w}}}^0\left({\displaystyle \mathbf{x}}\right),{\displaystyle \mathbf{c}}(0,{\displaystyle \mathbf{x}})={{\displaystyle \mathbf{c}}}^0\left({\displaystyle \mathbf{x}}\right),\end{array}$$

(10d)

and homogeneous Neumann boundary conditions

$$\begin{array}{c}\hfill \sigma \nabla V_{\mathrm{m}}·{\displaystyle \mathbf{n}}=0\mathrm{on}(0,T)\times \partial {\mathsf{\Omega }}_{\mathrm{EP}}.\end{array}$$

(10e)

The dimensions $n_{\mathbf{w}}$ and $n_{\mathbf{c}}$, the ionic current $I_{\mathrm{ion}}(V_{\mathrm{m}},\mathbf{w},\mathbf{c})$, the vector functions ${\mathbf{G}}_{\mathbf{w}}(V_{\mathrm{m}},\mathbf{w},\mathbf{c})$ and ${\mathbf{G}}_{\mathbf{c}}(V_{\mathrm{m}},\mathbf{w},\mathbf{c})$ are specified by the choice of the cell model. According to the Hodgkin–Huxley model [9], the total ionic current can be written as

$$\begin{array}{c}\hfill I_{\mathrm{ion}}(V_{\mathrm{m}},{\displaystyle \mathbf{w}},{\displaystyle \mathbf{c}})=\sum _{i=1}^k{g}_i\left(\prod _{j=1}^{n_{\mathbf{w}}}{w}_j^{p_{i,j}}\right)\left(V_{\mathrm{m}}-E_i\right)+H(V_{\mathrm{m}},{\displaystyle \mathbf{w}},{\displaystyle \mathbf{c}}),\end{array}$$

with the maximal conductance $g_i$ and the reversal potential $E_i$ for each ion channel i. More complex models add terms to the classical Hodgkin-Huxley model for additional ion channels. Those are of the form $I=g\left({\displaystyle \prod _{j=1}^{n_{\mathbf{w}}}}{w}_j^{p_j}\right)(V_{\mathrm{m}}-E)$ and increase the number k or are summarized in the function $H(V_{\mathrm{m}},\mathbf{w},\mathbf{c})$. The components of ${\mathbf{G}}_{\mathbf{w}}(V_{\mathrm{m}},\mathbf{w},\mathbf{c})$ are typically given by

$$\begin{array}{c}\hfill G_{w_j}(V_{\mathrm{m}},w_j)={\alpha }_j\left(V_{\mathrm{m}}\right)\left(1-w_j\right)-{\beta }_j\left(V_{\mathrm{m}}\right)w_j\mathrm{for}j=1,\dots ,n_{\mathbf{w}},\end{array}$$

where for $j=1,\dots ,n_{\mathbf{w}}$ the positive functions ${\alpha }_j\left(V_{\mathrm{m}}\right)$ and ${\beta }_j\left(V\right)$ are fitted to the experimental data measured for the specific gate j. Due to different physiological properties on the cellular level for the atria and the ventricles, different models should be chosen on the domains ${\mathsf{\Omega }}_{\mathrm{A}}$ and ${\mathsf{\Omega }}_{\mathrm{V}}$. The atrial and ventricular cell model used for simulations presented in this work are specified in Section 2.4.1.

The anisotropic conductivity tensor

$$\sigma \left({\displaystyle \mathbf{x}}\right)={\sigma }_{\mathrm{f}}{\displaystyle \mathbf{f}}\otimes {\displaystyle \mathbf{f}}+{\sigma }_{\mathrm{s}}{\displaystyle \mathbf{s}}\otimes {\displaystyle \mathbf{s}}+{\sigma }_{\mathrm{n}}{\displaystyle \mathbf{n}}\otimes {\displaystyle \mathbf{n}}\in {\mathbb{R}}_{\mathrm{sym}}^{3\times 3},$$

depends on the fiber, sheet, and sheet-normal direction (5) with conductivity parameters ${\sigma }_{\mathrm{f}}\ge {\sigma }_{\mathrm{s}}\ge {\sigma }_{\mathrm{n}}\ge 0$. The depolarization waves are initiated by the stimuli $I_{\mathrm{ext}}(t,\mathbf{x})$ in the stimulation area ${\mathsf{\Omega }}_{\mathrm{stim}}\subset {\mathsf{\Omega }}_{\mathrm{EP}}$, and we set

$$\begin{array}{c}\hfill I_{\mathrm{ext}}(t,{\displaystyle \mathbf{x}})=\left\{\begin{array}{cc}{A}_j\hfill & t\in (t_j,t_j+{\tau }_j),{\displaystyle \mathbf{x}}\in {\mathsf{\Omega }}_{\mathrm{stim}}\hfill \\ 0\hfill & \mathrm{else}\hfill \end{array}\right.\end{array}$$

(11)

depending on the starting time $t_j\ge 0$, duration ${\tau }_j>0$, and amplitude $A_j>0$.

2.4. Electro-Mechanical Coupling Mechanisms

2.4.1. Cellular Electro-Mechanical Model

The mathematical description of electrophysiological activity on the cellular level is represented by the term $I_{\mathrm{ion}}$ in Equation (10a) coupled to the ODEs (10b) and (10c). In general, Equation (10b,c) can represent any model that describes gating mechanisms and the transport of ions. Here, the models by Courtemanche et al. [68] and O’Hara et al. [69] are used for the atria and the ventricles, respectively. For the model by O’Hara et al., the modifications to the h- and j-gate as proposed by [70,71] were implemented.

Active tension development for $\mathbf{q}:[0,T]\times {\mathsf{\Omega }}_{\mathrm{EP}}\to {\mathbb{R}}^{n_{\mathbf{q}}}$ is represented by a recent model for human cardiac contraction proposed by Land et al. [72]. It consists of a system of ODEs summarized as

$${\partial }_t{\displaystyle \mathbf{q}}={{\displaystyle \mathbf{G}}}_{\mathbf{q}}({\displaystyle \mathbf{q}},{\displaystyle \mathbf{c}},{\gamma }_{\mathbf{f}},\dot{{\gamma }_{\mathbf{f}}})\mathrm{in}(0,T)\times {\mathsf{\Omega }}_{\mathrm{EP}},$$

(12)

where ${\gamma }_{\mathbf{f}}=\sqrt{{\left(\mathbf{F}\mathbf{f}\right)}^{\top }\left(\mathbf{F}\mathbf{f}\right)}$ is the stretch in fiber direction. Bidirectional coupling between the models of electrophysiology and active stress is established by replacing the algebraic formulation of the troponin buffer in the Courtemanche et al. and O’Hara et al. models by the evolution of calcium bound to troponin from the Land et al. model. Therefore, the ion concentration $\mathbf{c}$ depends on the stretch ${\gamma }_{\mathbf{f}}$ and Equation (10c) extends to

$${\partial }_t{\displaystyle \mathbf{c}}={{\displaystyle \mathbf{G}}}_{\mathbf{c}}(V_{\mathrm{m}},{\displaystyle \mathbf{w}},{\displaystyle \mathbf{c}},{\gamma }_{\mathbf{f}})\mathrm{in}(0,T)\times {\mathsf{\Omega }}_{\mathrm{EP}}.$$

(13)

The components of (13) for the intracellular concentration of calcium ions ${\left[{\mathrm{Ca}}^{2+}\right]}_{\mathrm{i}}$ change to

$$\begin{array}{cc}\hfill {\partial }_t{\left[{\mathrm{Ca}}^{2+}\right]}_{\mathrm{i}}& =\frac{1}{1+\frac{{\left[\mathrm{CMDN}\right]}_{\max }{K}_{\mathrm{CMDN}}}{{({\left[{\mathrm{Ca}}^{2+}\right]}_{\mathrm{i}}+K_{\mathrm{CMDN}})}^2}}\hfill \\ \hfill & ·\left[\frac{2I_{\mathrm{NaCa}}-I_{\mathrm{p},\mathrm{Ca}}-I_{\mathrm{Ca},\mathrm{L}}-I_{\mathrm{b},\mathrm{Ca}}}{2F{V}_{\mathrm{i}}}+\frac{V_{\mathrm{up}}(I_{\mathrm{up},\mathrm{leak}}-I_{\mathrm{up}})+I_{\mathrm{rel}}{V}_{\mathrm{rel}}}{V_{\mathrm{i}}}-{\partial }_t{\left[{\mathrm{Ca}}^{2+}\right]}_{\mathrm{TRPN}}\right],\hfill \end{array}$$

(14)

for the model by Courtemanche et al. and

$$\begin{array}{cc}\hfill {\partial }_t{\left[{\mathrm{Ca}}^{2+}\right]}_{\mathrm{i}}& =\frac{1}{1+\frac{{\left[\mathrm{CMDN}\right]}_{\max }{K}_{\mathrm{CMDN}}}{{({\left[{\mathrm{Ca}}^{2+}\right]}_{\mathrm{i}}+K_{\mathrm{CMDN}})}^2}}\hfill \\ \hfill & \left[-(I_{\mathrm{pCa}}+I_{\mathrm{Cab}}-2I_{\mathrm{NaCa}})\frac{A_{\mathrm{cap}}}{2F{v}_{\mathrm{myo}}}-J\frac{v_{\mathrm{nsr}}}{v_{\mathrm{myo}}}+J_{\mathrm{Ca}}\frac{v_{\mathrm{ss}}}{v_{\mathrm{myo}}}-{\partial }_t{\left[{\mathrm{Ca}}^{2+}\right]}_{\mathrm{TRPN}}\right],\hfill \end{array}$$

(15)

for the model by O’Hara et al. with

$$\begin{array}{cc}\hfill {\partial }_t{\left[{\mathrm{Ca}}^{2+}\right]}_{\mathrm{CaTRPN}}& ={\left[\mathrm{TRPN}\right]}_{\max }·{\partial }_t\mathrm{CaTRPN},\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill {\partial }_t\mathrm{CaTRPN}& =k_{\mathrm{CaTRPN}}\left[{\left(\frac{{\left[{\mathrm{Ca}}^{2+}\right]}_{\mathrm{i}}}{{\left[{\mathrm{Ca}}^{2+}\right]}_{\mathrm{T}50}\left({\gamma }_{\mathbf{f}}\right)}\right)}^{n_{\mathrm{CaTRPN}}}(1-\mathrm{CaTRPN})-\mathrm{CaTRPN}\right].\hfill \end{array}$$

(17)

CaTRPN is a component of $\mathbf{q}$ and represents the fraction of troponin C units with calcium bound to its regulatory binding site in the Land et al. model, ${\left[\mathrm{TRPN}\right]}_{\max }$ is the maximum concentration of troponin in the myoplasm, and ${\left[{\mathrm{Ca}}^{2+}\right]}_{\mathrm{T}50}\left({\gamma }_{\mathbf{f}}\right)$ is the length-dependent sensitivity to intra-cellular calcium.

The above modifications are introduced in [35,36] for the O’Hara et al. model. Parameters or equations not defined specifically are adopted from the original models and given in detail in the Supplementary Material.

Originally, the active stress model by Land et al. was developed and driven by experimentally measured human calcium transients (CaT) for the atria [34] and ventricles [72], respectively. Due to differences in the CaTs of the two models used in this study compared to the ones used by Land et al., a re-parameterization of the tension model is necessary to achieve physiologically correct active tension development. Based on the work of [34,36,72], we manually adjusted the parameters of the model to match key features of active tension to isometric twitch data of intact human tissue preparations from literature.

2.4.2. Mechano-Electric Feedback

As the conductivity tensor $\sigma$ is oriented along the myocyte orientation $\mathbf{f}$, $\mathbf{s}$, $\mathbf{n}$, Equation (10a) has to be evaluated on the deformed geometry ${\mathsf{\Omega }}^{\mathit{\phi }}$ for $t\in (0,T)$

$$\beta C_{\mathrm{m}}{\partial }_t^{\mathit{\phi }}{V}_{\mathrm{m}}={\nabla }^{\mathit{\phi }}·\sigma {\nabla }^{\mathit{\phi }}{V}_{\mathrm{m}}-\beta I_{\mathrm{ion}}(V_{\mathrm{m}},{\displaystyle \mathbf{w}},{\displaystyle \mathbf{c}})+\beta I_{\mathrm{ext}}\mathrm{in}{\mathsf{\Omega }}_{\mathrm{EP}}^{\mathit{\phi }},$$

(18)

where the derivatives are dependent on the deformed variable ${\mathbf{x}}^{\mathit{\phi }}=\mathit{\phi }(t,\mathbf{x})$. Applying the Piola transform on this equation, the corresponding Lagrange formulation reads

$$\beta C_{\mathrm{m}}{\partial }_t{V}_{\mathrm{m}}=\nabla ·(\left({{\displaystyle \mathbf{F}}}^{-1}\sigma {{\displaystyle \mathbf{F}}}^{-\top }\right)\nabla V_{\mathrm{m}})-\beta I_{\mathrm{ion}}(V_{\mathrm{m}},{\displaystyle \mathbf{w}},{\displaystyle \mathbf{c}})+\beta I_{\mathrm{ext}}\mathrm{in}(0,T)\times {\mathsf{\Omega }}_{\mathrm{EP}}.$$

(19)

Furthermore, the existence of stretch-activated currents $I_{\mathrm{SAC}}$ was first confirmed by Guharay and Sachs [73]. Different models were developed to describe the current $I_{\mathrm{SAC}}(V_{\mathrm{m}},{\gamma }_{\mathbf{f}})$ carried by these channels [74,75,76,77,78,79]. In this case, $I_{\mathrm{SAC}}$ is added to $I_{\mathrm{ion}}$, so that the ion current $I_{\mathrm{ion}}=I_{\mathrm{ion}}(V_{\mathrm{m}},\mathbf{w},\mathbf{c},{\gamma }_{\mathbf{f}})$ also depends on the stretch. For the simulations presented in Section 4, the deformation was not considered for the calculation of electrophysiology and no stretch-activated channels were incorporated.

2.5. Electro-Mechanical Coupling Algorithm

The fully coupled whole heart model used for the numerical simulations in Section 4 is described by the following PDE-ODE system

$$\begin{array}{cccc}\hfill \beta C_{\mathrm{m}}{\partial }_t{V}_{\mathrm{m}}& =\nabla ·\sigma \nabla V_{\mathrm{m}}-\beta I_{\mathrm{ion}}(V_{\mathrm{m}},{\displaystyle \mathbf{w}},{\displaystyle \mathbf{c}},{\gamma }_{\mathbf{f}})+\beta I_{\mathrm{ext}}\hfill & \hfill \mathrm{in}& (0,T)\times {\mathsf{\Omega }}_{\mathrm{EP}},\hfill \end{array}$$

(20a)

$$\begin{array}{cccc}\hfill {\partial }_t{\displaystyle \mathbf{w}}& ={{\displaystyle \mathbf{G}}}_{\mathbf{w}}(V_{\mathrm{m}},{\displaystyle \mathbf{w}},{\displaystyle \mathbf{c}})\hfill & \hfill \mathrm{in}& (0,T)\times {\mathsf{\Omega }}_{\mathrm{EP}},\hfill \end{array}$$

(20b)

$$\begin{array}{cccc}\hfill {\partial }_t{\displaystyle \mathbf{c}}& ={{\displaystyle \mathbf{G}}}_{\mathbf{c}}(V_{\mathrm{m}},{\displaystyle \mathbf{w}},{\displaystyle \mathbf{c}},{\gamma }_{\mathbf{f}})\hfill & \hfill \mathrm{in}& (0,T)\times {\mathsf{\Omega }}_{\mathrm{EP}},\hfill \end{array}$$

(20c)

$$\begin{array}{cccc}\hfill {\partial }_t{\displaystyle \mathbf{q}}& ={{\displaystyle \mathbf{G}}}_{\mathbf{q}}({\displaystyle \mathbf{q}},{\displaystyle \mathbf{c}},{\gamma }_{\mathbf{f}},{\dot{\gamma }}_{\mathbf{f}})\hfill & \hfill \mathrm{in}& (0,T)\times {\mathsf{\Omega }}_{\mathrm{EP}},\hfill \end{array}$$

(20d)

$$\begin{array}{cccc}\hfill {\displaystyle \mathbf{S}}({\displaystyle \mathbf{x}},{\displaystyle \mathbf{F}})& ={\mathrm{D}}_{\mathbf{F}}W\left({\displaystyle \mathbf{F}}\right)+{{\displaystyle \mathbf{S}}}_{\mathrm{a}}\hfill & \hfill \mathrm{in}& (0,T)\times {\mathsf{\Omega }}_{\mathrm{M}},\hfill \end{array}$$

(20e)

$$\begin{array}{cccc}\hfill {\varrho }_0{\partial }_t^2{\displaystyle \mathbf{u}}& =div\left({\displaystyle \mathbf{F}}{\displaystyle \mathbf{S}}\right({\displaystyle \mathbf{x}},{\displaystyle \mathbf{F}}\left)\right)\hfill & \hfill \mathrm{in}& (0,T)\times {\mathsf{\Omega }}_{\mathrm{M}},\hfill \end{array}$$

(20f)

$$\begin{array}{cccc}\hfill {\partial }_t{\displaystyle \mathbf{v}}& ={{\displaystyle \mathbf{G}}}_{\mathbf{v}}({\displaystyle \mathbf{p}},{\displaystyle \mathbf{z}})\hfill & \hfill \mathrm{in}& (0,T),\hfill \end{array}$$

(20g)

with boundary conditions

$$\begin{array}{cccc}\hfill \sigma \nabla V_{\mathrm{m}}·{\displaystyle \mathbf{n}}& =0\hfill & \hfill \mathrm{on}& (0,T)\times \partial {\mathsf{\Omega }}_{\mathrm{EP}},\hfill \end{array}$$

(20h)

$$\begin{array}{cccc}\hfill {\displaystyle \mathbf{u}}& ={\displaystyle \mathbf{0}}\hfill & \hfill \mathrm{on}& (0,T)\times {\mathsf{\Gamma }}_{\mathrm{D}},\hfill \end{array}$$

(20i)

$$\begin{array}{cccc}\hfill \left[\mathit{\phi }\right]& ={\displaystyle \mathbf{0}}\hfill & \hfill \mathrm{on}& (0,T)\times {\mathsf{\Gamma }}_{\mathrm{P}},\hfill \end{array}$$

(20j)

$$\begin{array}{cccc}\hfill {\displaystyle \mathbf{F}}{\displaystyle \mathbf{S}}({\displaystyle \mathbf{x}},{\displaystyle \mathbf{F}}){\displaystyle \mathbf{n}}& =p_{\mathrm{C}}{\displaystyle \mathbf{n}}\hfill & \hfill \mathrm{on}& (0,T)\times {\mathsf{\Gamma }}_{\mathrm{C}},\mathrm{C}\in \left\{\mathrm{LV},\mathrm{RV},\mathrm{LA},\mathrm{RA}\right\}.\hfill \end{array}$$

(20k)

This is complemented by the evaluation of the pressure vector $\mathbf{p}$ and the intermediate variables $\mathbf{z}$ as described in Section 2.2.2. As the cardiac electrophysiology is computed on ${\mathsf{\Omega }}_{\mathrm{EP}}$ and the elasto-mechanical equations are solved on ${\mathsf{\Omega }}_{\mathrm{M}}$, a mapping between the corresponding discrete meshes is required to couple the problem. This mapping from the vertices of the mesh of ${\mathsf{\Omega }}_{\mathrm{EP}}$ to the quadrature points of the finite element disctetization of ${\mathsf{\Omega }}_{\mathrm{M}}$ is denoted as ${\mathcal{M}}_{\mathrm{EP},\mathrm{M}}$ and the mapping vice versa as ${\mathcal{M}}_{\mathrm{M},\mathrm{EP}}$. Both mappings are realized via linear interpolation. In space, the PDEs in ${\mathsf{\Omega }}_{\mathrm{M}}$ and ${\mathsf{\Omega }}_{\mathrm{EP}}$ are discretized with the finite element method, using linear conforming tetrahedral elements for the electrophysiology and quadratic elements for the elasto-mechanic part.

In time, we employ a staggered approach. First, the equations for the electrophysiology defined on ${\mathsf{\Omega }}_{\mathrm{EP}}$ are solved, updating $V_{\mathrm{m}}$, $\mathbf{w}$, $\mathbf{c}$, $\mathbf{q}$. This is realized via a first order operator splitting method as proposed by Sundnes et al. [80] where (10) is split in the reaction and diffusion part. The space-independent ODE system for the reaction is solved in time by explicit integration methods. The linear PDE modeling the diffusion is solved by a Crank-Nicolson scheme in time. The details are given in the supplement. The time step $\Delta t_{\mathrm{EP}}=10$ µs is used, so that the monodomain equation is computed $\Delta t_{\mathrm{M}}/\Delta t_{\mathrm{EP}}$ times, before the mechanics is updated. Then, $T_{\mathrm{tot}}(\mathbf{c},\mathbf{q})$ is mapped via ${\mathcal{M}}_{\mathrm{M},\mathrm{EP}}$ to ${\mathsf{\Omega }}_{\mathrm{M}}$ to evaluate the active stress response (6). Now, the equations for elasto-mechanics are solved on ${\mathsf{\Omega }}_{\mathrm{M}}$. Starting with the update of the pressure vector $\mathbf{p}$ as explained in Section 2.2.2, a new displacement $\mathbf{u}$ is computed with a Newmark $\beta$-scheme, the details of which are again given in the Supplementary Material. Then, the chamber volumes of the discrete geometry are computed, cf. Section 2.2.2. Additionally, the volumes of the circulatory system $\mathbf{v}$ are updated with a fourth order Runge–Kutta scheme and time step $\Delta t_{\mathbf{p}}=0.1$ ms. Then, the new vector $\mathbf{z}$ is computed. If the difference between the computed volumes is larger than a tolerance value, the procedure starts with an update of $\mathbf{p}$ again. Else, the elasto-mechanics are finished for the time step and the stretch ${\gamma }_{\mathbf{f}}$ is mapped with ${\mathcal{M}}_{\mathrm{M},\mathrm{EP}}$ to ${\mathsf{\Omega }}_{\mathrm{EP}}$. The new stretch is used for the next update of the electrophysiology. The overall scheme is illustrated in Figure 3.

3. Patient-Specific Simulation and Evaluation

3.1. Personalizing Electro-Mechanical Whole Heart Models: Building Digital Twins

Personalized in silico cardiology is evolving to become an important component of therapy planing and starts to inform the decision making process throughout the clinical workflow by maximizing the value of available clinical data [8]. Building a digital twin can comprise different aspects including statistical and mechanistic models. In the following, a brief overview of the important aspects of cardiac modeling is given and we highlight strategies on how to incorporate personalized information in our modeling framework.

3.1.1. Cardiac Anatomy

The first step to a personalized in silico model is to create an anatomically accurate representation of the heart in a discretized finite element model. Medical imaging data plays an important role in this step, since it provides the individuals shape of the heart and structural information such as the location of a scar, fiber architecture, and fibrosis in a non-invasive procedure [81,82]. In the case of whole heart models, a semi-automatic workflow has been proposed by Strocchi et al. [39]. However, a manual adjustment of the discretized model is still necessary and multiple tools are available to make this process more efficient (e.g., meshtool [83], vmtk [84]). Since anatomy is highly individualized, it is necessary to ensure the comparability between datasets. To solve this problem, a standardized 17 sector map of the ventricles can be used [85]. Similarly, universal coordinate systems have been proposed for the ventricles [41,42] and atria [40].

3.1.2. Fiber Orientation

Recovering patient-specific fiber and sheet orientation is possible through diffusion tensor MRI [86]. However, this kind of imaging data are only available from ex vivo measurements and is not acquired in vivo and thus, Laplace–Dirichlet rule-based (LDRB) methods are used to define fiber and sheet orientations in the ventricles [87,88,89] and atria [58,90]. Nevertheless, information about the transmural distribution of fiber and sheet angles from observations [91,92] can be used to parametrize LDRB methods to specify fiber and sheet angles at the endocardium and epicardium (${\alpha }_{\mathrm{endo}}$,${\alpha }_{\mathrm{epi}}$, ${\beta }_{\mathrm{endo}}$,${\beta }_{\mathrm{epi}}$).

3.1.3. Passive Stress

MRI data are typically acquired during the early diastolic state of the heart cycle when movement and chamber pressure are minimal. However, residual stress in the tissue cannot be measured with standard imaging techniques and the pressure inside the heart’s cavities can only be measured by invasive procedures. Therefore, a pressure and stress-free reference configuration of the heart has to be estimated to accurately model the biomechanical diastolic function. Methods applied in the context of cardiac mechanics try to find a stress-free reference configuration by using fixed-point iterations, as originally proposed by Sellier et al. [93]. Parameters for the constitutive model are mostly chosen manually [94] by utilizing the empirical end-diastolic pressure-volume relation (EDPVR) described by Klotz et al. [95] as a fitting target. In Kovacheva et al. [43], the constitutive parameters were determined by solving an optimization problem using a gradient-free method: the distance between the simulated EDPVR and the one proposed by Klotz et al. was minimized together with an additional condition imposed on the unloaded volume. The latter volume relation was again proposed by Klotz et al. However, values from different sources contradict each other [96]. Recently, Marx et al. [44] presented an automated method to match passive material parameters to the shape of a patient-specific EDPVR and to find an unloaded reference configuration of the heart simultaneously. Their method requires a mesh generated from clinical imaging data and at least one data point of the EDPVR.

3.1.4. Active Stress

Personalizing the systolic part of the heart cycle is done by adjusting the characteristics of the active tension model. This is done until the deformation of the heart is similar to the one observed in MRI data or until the ejected blood volume of the recorded heart beat is reached. Personalizing other characteristics of active tension development, such as the rate of contraction is more difficult in biophysically motivated models, due to the high number of parameters, which are typically not well constrained by the available experimental data [72]. Therefore, most personalized electro-mechanical models utilize a simplified model to drive contraction [33,61]. These models only comprise a small set of parameters that are linked to distinct characteristics (e.g., peak tension, rate of contraction and relaxation) of tension development, making them easier to fit to experimental data. As a consequence, these simpler models miss the important connection to calcium and are driven by activation times only.

3.1.5. Electrophysiology

Myocardial tissue in the atria and the ventricles shows regionally varying anisotropic properties with regard to conduction velocity (CV) [97]. Considering that CV is not a parameter of the momodomain equation, we cannot adjust it directly. However, assuming a planar wavefront propagation, Costa et al. [98] found that CV is proportionally related to the conductivity $\sigma$ and the surface-to-volume ratio $\beta$ by $\mathrm{CV}\propto \sqrt{\sigma /\beta }$. Since CV is much more likely to be measured in clinical scenarios than conductivity (e.g., [51,99]), one way to personalize the conduction model is to match CV. Costa et al. [98] realized this and proposed an automated parameterization strategy to find optimal conductivity values for a given CV. This strategy should be applied to each patient-specific mesh due to the dependency of CV on the spatial discretization [100] and the cellular models used in the study.

Mostly, atrial and ventricular activation is initiated by an externally applied stimulus $I_{\mathrm{ext}}$ mimicking the effect of the specialized cardiac conduction system at the sinus node and the Purkinje muscle junction. Timing and position of the stimulus are crucial for a proper depolarization pattern and orchestrated mechanical activation of the heart. In the atria, a stimulation at the sinus node is enough to initiate a heart beat. However, it is not trivial to properly activate the ventricles. This is due to the conduction of the depolarization wave through the atrio-ventricular node (AVN), the His–Purkinje system (HPS), and consequently the ventricular myocardium. Especially the HPS is highly individualized and rarely modeled explicitly [101] but rather implicitly by capturing its effect at the Purkinje muscle junctions [49,57,102]. Respective activation sequences can be personalized by using in-vivo measurements like electrocardiograms (ECG) or body surface potential maps (BSPM) to estimate sites of earliest activation [49,102,103]. Alternatively, the HPS can be modeled using N fascicular sites for the stimulus combined with a thin, fast conducting endocardial layer [104]. The conduction delay introduced by the AVN is typically represented by a predefined time delay between the onset of atrial and ventricular depolarization or by a model representing the AVN function [105].

3.2. Experimental Setup

Two simulations were conducted: one representing a normal heart beat as a reference and a second including a scar in the left atrium as a result of a typical radio-frequency ablation therapy for atrial fibrillation.

3.2.1. Parameterization

To capture the highly anisotropic characteristics of myocardial tissue, fiber and sheet architecture was incorporated into the model. For the ventricles, a method based on Bayer et al. [87] was applied. The original algorithm was adapted (https://github.com/KIT-IBT/LDRB_Fibers (accessed on 23 May 2021), release v0.1) [106] to eliminate a discontinuity of the fibers in the free walls and to yield a fiber rotation that is proportional to the transmural Laplace solution. As shown in Figure 4, fiber angles change linearly from ${\alpha }_{\mathrm{endo}}={60}^{\circ }$ at the endocardium to ${\alpha }_{\mathrm{epi}}=-{60}^{\circ }$ at the epicardium, as suggested by Streeter et al. [91]. Furthermore, the sheet angles change from ${\beta }_{\mathrm{endo}}=-{65}^{\circ }$ to ${\beta }_{\mathrm{epi}}={25}^{\circ }$ at the endo- and epicardium, respectively [87]. The atrial fiber and sheet orientation in this work are based on 22 seed points that define anatomical landmarks used to calculate numerous paths which split the atria in regions with distinct fiber orientation [58,59].

Using the method proposed by Costa et al. [98], $\sigma$ was optimized to achieve CVs that resulted in a realistic activation sequence. To simplify atrial activation, only three different conductivities were assigned to differentiate atrial bulk tissue, fast conducting bundles, and slow conducting regions [107]. The algorithm by Wachter et al. [58] was used to label the fast conducting regions including the crista terminalis, the pectinate muscles, Bachmann’s bundle, middle and upper posterior inter-atrial connection, the coronary sinus, and the scar in the left atrium (Figure 5).

They were assumed to be transversely isotropic with CVs of $2.25$ m/s in longitudinal direction and $1.45$ m/s in transverse direction. Slow regions were set to be isotropic with a CV of $0.65$ m/s, which is a commonly used simplification of the structural reality [107]. Atrial bulk tissue used the same value for a CV of $1.45$ m/s in longitudinal direction and $0.65$ m/s in transverse direction.

Ventricular myocardium was assumed to be orthotropic with CVs of $0.87$ m/s, $0.6$ m/s, $0.4$ m/s in longitudinal, transversal, and normal directions respectively. Additionally, a thin, fast-conducting endocardial layer was labeled using the consistent ventricular coordinates (Cobiveco) [42]. The CV in this layer was three times higher than in the rest of the ventricle to represent the effect of the HPS. All parameters used in the whole heart simulations are listed in

Table 1. Overview of electrophysiological parameters for the whole heart model.

Parameter	Value	Unit	Description
$({\sigma }_{\mathrm{f}},{\sigma }_{\mathrm{s}},{\sigma }_{\mathrm{n}})$	$(0.47,0.27,0.15)$	S/m	conductivities in ventricular bulk tissue
$({\sigma }_{\mathrm{f}},{\sigma }_{\mathrm{s}},{\sigma }_{\mathrm{n}})$	$(3.25,2.75,2.25)$	S/m	conductivities in ventricular fast conducting layer
$({\sigma }_{\mathrm{f}},{\sigma }_{\mathrm{s}},{\sigma }_{\mathrm{n}})$	$(0.99,0.26,0.26)$	S/m	conductivities in atrial bulk tissue
$({\sigma }_{\mathrm{f}},{\sigma }_{\mathrm{s}},{\sigma }_{\mathrm{n}})$	$(2.35,0.99,0.99)$	S/m	conductivities in atrial fast conducting regions
$({\sigma }_{\mathrm{f}},{\sigma }_{\mathrm{s}},{\sigma }_{\mathrm{n}})$	$(0.26,0.26,0.26)$	S/m	conductivities in atrial slow conducting regions
$({\sigma }_{\mathrm{f}},{\sigma }_{\mathrm{s}},{\sigma }_{\mathrm{n}})$	$({10}^{-12},{10}^{-12},{10}^{-12})$	S/m	conductivities in scar tissue
$\beta$	140,000	1/m	membrane surface-to-volume ratio
$C_{\mathrm{m}}$	$0.01$	F/m${}^2$	membrane capacitance
AV-delay	$0.160$	s	atrio-ventricular conduction delay
BCL	1	s	basic cycle length (=$1/\mathrm{heartrate}$)

.

The depolarization wave was initiated by a stimulus $I_{\mathrm{ext}}$ with a cycle length of $\mathrm{BCL}=1$ s at the position of the most common sinus node exit site at the junction of the right atrial appendage and the superior vena cava [108]. Conduction between the left and right atrium was only possible via Bachmann’s bundle, a middle and upper posterior inter-atrial connection, and the coronary sinus. Otherwise, the tissue was isolated in the middle of the septal wall through duplicated vertices. 160 ms after the atrial stimulus, the ventricles were excited at five distinct positions. These five positions represent common sites of earliest activation in the ventricles and were chosen based on observations by Durrer et al. [109] and simulation results obtained in [104]. The activation pattern includes three sites of earliest activation in the LV (mid-posterior inferior ${\mathbf{x}}_{\mathrm{LV},\mathrm{mpi}}$, mid-posterior superior ${\mathbf{x}}_{\mathrm{LV},\mathrm{mps}}$, basal anterior paraseptal ${\mathbf{x}}_{\mathrm{LV},\mathrm{bap}}$) and two in the RV (septal ${\mathbf{x}}_{\mathrm{RV},\mathrm{s}}$ and free wall ${\mathbf{x}}_{\mathrm{RV},\mathrm{fw}}$). Their position is given in terms of universal coordinates in

Table 2. Sites of earliest activation in terms of the coordinate system Cobiveco [42]: a apicobasal; m transmural; r rotational; v transventricular.

	Root Point	Extent
$\mathbf{x}$	${\mathbf{x}}_{\mathbf{root}}=\{\mathit{a},\mathit{m},\mathit{r},\mathit{v}\}$	${\mathbf{\delta }}_{\mathbf{m}}$	${\mathbf{\delta }}_{\mathbf{rad}}$
${\mathbf{x}}_{\mathrm{LV},\mathrm{mps}}$	$\{0.65,1,0.99,0\}$	$0.05$	3 mm
${\mathbf{x}}_{\mathrm{LV},\mathrm{mpi}}$	$\{0.55,1,0.175,0\}$	$0.05$	3 mm
${\mathbf{x}}_{\mathrm{LV},\mathrm{bap}}$	$\{0.9,1,0.57,0\}$	$0.05$	3 mm
${\mathbf{x}}_{\mathrm{RV},\mathrm{s}}$	$\{0.4,1,0.85,1\}$	$0.05$	3 mm
${\mathbf{x}}_{\mathrm{RV},\mathrm{fw}}$	$\{0.45,1,0.15,1\}$	$0.1$	3 mm

and the resulting activation times are shown in Figure 6.

In the mechanical model, myocardial tissue in the atria and ventricles was modeled as a transversely isotropic material as defined in Equation (4). These parameters were determined using the method proposed by Kovacheva et al. [43] to match the empirical EDPVR of Klotz et al. [95]. Purely passive tissue was modeled as an isotropic Neo-Hookean solid using different material parameters. All mechanical parameters are listed in

Table 3. Overview of passive mechanical parameters in $W_{\mathrm{G}},W_{\mathrm{NH}}$ for the whole heart model.

		Parameters
Domain	Model	$\mathit{\mu }$ (Pa)	${\mathit{b}}_{\mathit{f}}$	${\mathit{b}}_{\mathit{s}}$	${\mathit{b}}_{\mathit{f}\mathit{s}}$	$\mathit{\kappa }$ (Pa)	${\mathit{\rho }}_0$(kg/m${}^2$)
${\mathsf{\Omega }}_{\mathrm{V}}$	Guccione	$325.56$	$11.01$	$4.4$	$7.71$	${10}^6$	1082
${\mathsf{\Omega }}_{\mathrm{A}}$	Guccione	$325.56$	$11.01$	$4.4$	$7.71$	${10}^6$	1082
${\mathsf{\Omega }}_{\mathrm{Valves}}$	Neo-Hooke	${10}^5$	-	-	-	${10}^3$	1082
${\mathsf{\Omega }}_{\mathrm{AdiposeTissue}}$	Neo-Hooke	3725	-	-	-	${10}^3$	1082
${\mathsf{\Omega }}_{\mathrm{MajorVessels}}$	Neo-Hooke	${10}^4$	-	-	-	${10}^3$	1082
${\mathsf{\Omega }}_{\mathrm{Pericardium}}$	Neo-Hooke	${10}^4$	-	-	-	${10}^3$	1082

.

For our particular subject, no measured data were available that can be related to circulatory system parameters. Therefore, all values were chosen based on models found in literature that share some of their structure with the model proposed in Section 2.2.2. The list of parameter values and corresponding references can be found in

Table 4. Overview of circulatory system parameters in Section 2.2.2 for the whole heart model.

Parameter	Value	Unit	Description	Ref.
$R_{\mathrm{SysArtValve}}$	0.006	$\mathrm{mmHg}·\mathrm{s}·{\mathrm{mL}}^{-1}$	aortic valve resistance	[110,111]
$R_{\mathrm{SysArt}}$	0.03	$\mathrm{mmHg}·\mathrm{s}·{\mathrm{mL}}^{-1}$	systemic arterial resistance	[112,113]
$C_{\mathrm{SysArt}}$	3.0	$\mathrm{mL}·{\mathrm{mmHg}}^{-1}$	systemic arterial compliance	[112,113,114]
$V_{\mathrm{SysArtUnstr}}$	800.0	mL	unstressed systemic arterial volume	[110]
$R_{\mathrm{SysPer}}$	0.6	$\mathrm{mmHg}·\mathrm{s}·{\mathrm{mL}}^{-1}$	systemic peripheral resistance	[112,113]
$R_{\mathrm{SysVen}}$	0.03	$\mathrm{mmHg}·\mathrm{s}·{\mathrm{mL}}^{-1}$	systemic venous resistance	[115,116]
$C_{\mathrm{SysVen}}$	150.0	$\mathrm{mL}·{\mathrm{mmHg}}^{-1}$	systemic venous compliance	[110,111,116]
$V_{\mathrm{SysVenUnstr}}$	2850.0	mL	unstressed systemic venous resistance	[110,115]
$R_{\mathrm{RavValve}}$	0.003	$\mathrm{mmHg}·\mathrm{s}·{\mathrm{mL}}^{-1}$	tricuspid valve resistance	[20,111]
$R_{\mathrm{PulArtValve}}$	0.003	$\mathrm{mmHg}·\mathrm{s}·{\mathrm{mL}}^{-1}$	pulmonary valve resistance	[110]
$R_{\mathrm{PulArt}}$	0.02	$\mathrm{mmHg}·\mathrm{s}·{\mathrm{mL}}^{-1}$	pulmonary arterial resistance	[117,118]
$C_{\mathrm{PulArt}}$	10.0	$\mathrm{mL}·{\mathrm{mmHg}}^{-1}$	pulmonary arterial compliance	[117,119]
$V_{\mathrm{PulArtUnstr}}$	150.0	mL	unstressed pulmonary arterial volume	[110]
$R_{\mathrm{PulPer}}$	0.07	$\mathrm{mmHg}·\mathrm{s}·{\mathrm{mL}}^{-1}$	pulmonary peripheral resistance	[120,121]
$R_{\mathrm{PulVen}}$	0.03	$\mathrm{mmHg}·\mathrm{s}·{\mathrm{mL}}^{-1}$	pulmonary venous resistance	[20]
$C_{\mathrm{PulVen}}$	15.0	$\mathrm{mL}·{\mathrm{mmHg}}^{-1}$	pulmonary venous compliance	[119]
$V_{\mathrm{PulVenUnstr}}$	200.0	mL	unstressed pulmonary venous volume	[110]
$R_{\mathrm{LavValve}}$	0.003	$\mathrm{mmHg}·\mathrm{s}·{\mathrm{mL}}^{-1}$	mitral valve resistance	[20]

.

The radio-frequency ablation scar only affected electrophysiological properties in this study. Specifically, the scar was assumed to be non-conducting and consists of a single lesion around all pulmonary veins isolating the whole roof of the left atrium plus a linear lesion towards the atrioventricular plane on the posterior side (Figure 5).

3.2.2. Initialization

As a first step, the cellular models were paced to a limit cycle at a cycle length of 1 s for a total of 1000 cycles. In this step, stretch ${\gamma }_{\mathbf{f}}$ and stretch rate ${\dot{\gamma }}_{\mathbf{f}}$ were set to 1 and 0 respectively. The resulting values of the state variables $\{\mathbf{w},\mathbf{c},\mathbf{q}\}$ were assigned to each vertex of ${\mathsf{\Omega }}_{\mathrm{EP}}$ as initial values. These values can be found in the Supplementary Material.

A backward displacement method [93,122] was used to find the unloaded configuration of ${\mathsf{\Omega }}_{\mathrm{M}}$ for all four cavities simultaneously. As no measured pressure data was available, population-based values for the diastatic pressure were assumed ($p_{\mathrm{LV}}=p_{\mathrm{LA}}=8.25$ mmHg, $p_{\mathrm{RV}}=p_{\mathrm{RA}}=3.5$ mmHg). The algorithm was stopped manually after 3 iterations with a residual norm of 0.0052 m and unloaded volumes $v_{\mathrm{RV}}=49.8\%,v_{\mathrm{LV}}=52.2\%,$ $v_{\mathrm{RA}}=55.0\%,$ $v_{\mathrm{LA}}=55.2\%$ with respect to the reference configuration. Afterwards, the unloaded configuration was inflated with the same pressure to pre-stress the tissue. Both, the unloaded and pre-stressed geometries are shown in the Supplementary Material together with the simulated EDPVR. To reduce the amount of heart beats that need to be run to reach a stable limit cycle of the circulatory system, purely mechanical simulations were done to find more suitable initial conditions for the circulatory model. The resulting initial values from which the coupled reference simulations was started are shown in

Table 5. Initial conditions for the circulatory system model.

Parameter	Value	Unit	Description
$V_{\mathrm{tot}}$	5700.0	mL	total volume
$V_{\mathrm{SysArt}}$	981.1396	mL	systemic arterial volume
$V_{\mathrm{PulArt}}$	303.7683	mL	pulmonary arterial volume
$V_{\mathrm{PulVen}}$	349.6759	mL	pulmonary venous volume
$p_{\mathrm{LV}}$	8.0246	mmHg	left ventricular pressure
$p_{\mathrm{LA}}$	8.2061	mmHg	left atrial pressure
$p_{\mathrm{RV}}$	5.8073	mmHg	right ventricular pressure
$p_{\mathrm{RA}}$	5.8071	mmHg	right atrial pressure

.

3.2.3. Evaluation

To evaluate the simulation results, the deformation of the heart was analyzed in terms of atrio-ventricular valve displacement (AVPD), left ventricular blood volume, and intra-cavitary pressure-volume relationships. For the left ventricular volume and the AVPD, the in silico results were compared to time-resolved MR imaging data that were analyzed using the freely available software Segment (version 3.2 R8531, [123]). AVPD was measured by projecting the mean displacement of the surrounding tissue of the mitral- and tricuspidvalve onto the apico-basal heart axis. Pressure-volume relationships are provided as output by the circulatory system model. Since the circulatory model is a closed loop, all simulations were run for multiple cycles until the difference in stroke volumes of the left and right ventricle was less than 1 mL.

4. Results

4.1. Cellular Electro-Mechanical Model

The calcium transients (CaT) of the cell models differ from the experimental traces used by Land et al. to calibrate their model regarding key biomarkers such as time to peak of the calcium transient (TPCaT) and relaxation times to 50% and 90% decay from the peak calcium concentration (RT50, RT90). Therefore, a re-parameterization of the tension model by Land et al. was necessary after introducing the bidirectional coupling into the models of Courtemanche et al. and O’Hara et al.

For the atrial model, the parameters from Land et al. [34] were used. In tissue simulations this choice of parameters caused unphysiological behavior due to the introduction of length-dependent effects. An unphysiological rise in tension in the atria was observed during the contraction of the ventricles. The passive increase in volume and the atrio-ventricular plane displacement (AVPD) resulted in the atria being stretched by the ventricles, thus prolonging the time of contraction of the atria. Land et al. [34] encountered the same problem and suggested faster atrial contraction and crossbridge cycling rates as a solution. The CaT of the Courtemanche et al. electrophysiological model already promoted a faster contraction with faster TPCaT compared to physiological values. Therefore, calcium cycling rates were left unchanged. Instead, the change in calcium sensitivity with respect to ${\gamma }_{\mathbf{f}}$ was reduced to ${\beta }_1=-0.5$ and the half-activation point was increased to ${\left[{\mathrm{Ca}}^{2+}\right]}_{\mathrm{T}50}({\gamma }_{\mathbf{f}}=1)=1.05$ µM. The re-parameterized atrial electro-mechanical model (Courtemanche–Land opt.) showed a time to peak tension TPT $=73.5$ ms and RT50 $=78.1$ ms after a total of 1000 cycles at a basic cycle length (BCL) of 1 s (Figure 7), which is in close agreement with atrial human tissue preparation data (

Table 6. Literature values for time to peak of calcium transient (TPCaT) and active tension development (TPT) as well as relaxation time to 50%, 90% and 95% respectively (RT50, RT90, RT95) from human tissue preparations. The list of ventricular values were originally compiled in [36].

	Calcium Transient			Active Tension
Tissue	TPCaT (ms)	RT50 (ms)	RT90 (ms)	TPT (ms)	RT50 (ms)	RT95 (ms)	Ref.
Ventricle	$47.8\pm 10.0$	$151.1\pm 89.2$	$315.6\pm 161.2$	-	-	-	[124]
Ventricle	-	-	-	$165\pm 7$	$116\pm 6$	$334\pm 43$	[125]
Ventricle	-	-	-	$157\pm 10$	$117\pm 8$	$477\pm 31$	[126]
Ventricle	-	-	-	$235.0\pm 13.4$	$153\pm 71$	$309\pm 13.7$	[127]
Ventricle	-	-	-	$151.0\pm 6.1$	$98.0\pm 7.7$	$173.0\pm 10.7$	[127]
Atria	$52.5\pm 3.1$	$177.5\pm 9.0$	-	$109.6\pm 3.6$	$110.2\pm 84.0$	-	[128]
Atria	-	-	-	$85.0\pm 5.5$	$66.1\pm 5.9$	-	[129]
Atria	-	-	-	$88.3\pm 2.5$	$73.3\pm 1.7$	-	[129]

).

As a baseline for the ventricular model, the parameterization of [72] is used. First, the modifications suggested by Margara et al. [36] to the Hill coefficient of cooperative activation and the tropomyosin rate constant were adopted. Additionally, the same value for the half-activation point ${\left[{\mathrm{Ca}}^{2+}\right]}_{\mathrm{T}50}({\gamma }_{\mathbf{f}}=1)=1.05$ µM was used as for the atrial model. This reduced diastolic resting tension, which otherwise exceeded 2 kPa in tissue simulations in pre-stressed (${\gamma }_{\mathbf{f}}>1$) conditions. The CaT and active tension development of the final ventricular model is shown in Figure 8. Compared to the original model by Land et al. [72], the optimized model (OHaraRudy+Land opt.) showed more physiological behavior compared to experimental data after stimulating the model for 1000 cycles at a BCL of 1 s (TPT $=154.8$ ms, RT50 $=121.7$ ms, RT95 $=275.0$ ms). To achieve deformations comparable with MRI data, the parameter $T_{\mathrm{ref}}$ in the OHaraRudy+Land opt. model was set to 480 kPa in multi-scale simulations. All adjusted parameters are given in

Table 7. Optimized parameters for the atrial and ventricular tension development model. The original values are taken from [34,72]. Values not listed in this table are unchanged from the original publications.

	Atria		Ventricle
Parameter	Original Value	Optimized	Original Value	Optimized
$k_{\mathrm{u}}$	1/ms	1/ms	1/ms	0.04/ms
$n_{\mathrm{Tm}}$	5	5	5	2.4
${\left[{\mathrm{Ca}}^{2+}\right]}_{\mathrm{T}50}^{\mathrm{ref}}$	0.86 µM	1.05 µM	0.805 µM	1.05 µM
${\beta }_1$	−2.4	−0.5	−2.4	−2.4

.

4.2. Electro-Mechanical Whole-Heart Simulation

The final discretizations were created using Gmsh [130]. ${\mathsf{\Omega }}_{\mathrm{M}}$ consists of 26,681 quadratic tetrahedral elements with an average edge length of 7.2 mm, yielding 48,780 nodes with 135,897 degrees of freedom. Since electrophysiological simulations require a finer spatial discretization [100], a subset of ${\mathsf{\Omega }}_{\mathrm{M}}$ was rediscretized to yield ${\mathsf{\Omega }}_{\mathrm{EP}}$ with 7,434,101 linear tetrahedral elements and an average edge length of 0.6 mm. All simulations were run on a 2019 Apple iMac™ using 8 MPI processes. A single heart beat took between 20 to 24 h to compute.

Starting from the pre-stressed reference geometry as described in Section 3.2, the reference heartbeat was simulated using the geometry shown in Figure 1. After eight beats, the simulation reached a stable limit cycle with a stroke volume (SV) difference between the two ventricles of less than 1 mL. As shown in Figure 9, the four major phases of the cardiac cycle were reproduced faithfully. After initiating the heart beat with an electrical stimulus at the location of the sinus node (Figure 9a), the atria contract and ventricular end-diastole later reached 180 ms (Figure 9b). Followed by a short time of isovolumetric contraction, the aortic and pulmonary valves open and the ejection of blood from the ventricles begins. Meanwhile, the atria relax and passively fill with blood (Figure 9c). During isovolumetric relaxation (Figure 9d), all valves close and the pressure reduces. As soon as the pressure is low enough, the mitral and tricuspid valves open and the ventricles fill with blood.

The results of the circulatory system are shown in Figure 10 for the last simulated heartbeat. Atrial PV-loops exhibit the typical figure-of-eight shape, with the A-wave and the V-wave both present. Additionally, a short spike in pressure with little change in volume (C-wave) can be observed in both atria immediately after the onset of ventricular contraction. With ${\mathrm{SV}}_{\mathrm{LV}}=126.16$ mL and ${\mathrm{SV}}_{\mathrm{RV}}=125.24$ mL, ventricular stroke volumes are close to the values derived from the imaging data (${\mathrm{SV}}_{\mathrm{LV}}^{\mathrm{MRI}}=132$ mL, ${\mathrm{SV}}_{\mathrm{RV}}^{\mathrm{MRI}}=129$ mL). However, systolic blood pressure is elevated in both the left ($p_{\mathrm{LV}}=154.6$ mmHg) and the right ventricle ($p_{\mathrm{RV}}=48.7$ mmHg). Flow through the mitral and tricuspid valve is observed during early diastolic filling (E-wave) and during atrial contraction (A-wave). With an E/A ratio of 1.11 in the left ventricle and 0.85 in the right ventricle, both are considered to function normally.

The LV and RV volume at the beginning of the last heart beat demonstrates a discrepancy with respect to the segmented volumes from the MRI data and the initial volumes in the pre-stressed geometry. The volume reduced from 200.5 mL to 179.4 mL in the LV and increased from 188.5 mL to 224.8 mL in the RV during equilibration to a stable limit cycle. Therefore, LV volume was normalized for both simulated and MRI results for comparison (Figure 11). In both control and ablation cases, atrial systole contributed 11–13% to the end-diastolic ventricular blood volume. Major differences can be observed during ventricular systole and diastole. With a peak ejection rate (PER) of ${\mathrm{PER}}_{\mathrm{LV}}=2094.8$ mL/s, it took only 160 ms until end-systole was reached. The ${\mathrm{PER}}_{\mathrm{LV}}^{\mathrm{MRI}}=882$ mL/s in the MRI data was markedly lower resulting in a time to end-systole of 398 ms. The opposite behavior was observed during ventricular relaxation. Peak filling rate (PFR) in the MRI data was more than two times higher than in the simulation (${\mathrm{PFR}}_{\mathrm{LV}}^{\mathrm{MRI}}=1027$ mL/s and ${\mathrm{PFR}}_{\mathrm{LV}}=408.3$ mL/s). Hence, ventricular relaxation in the model was significantly slower than in the MRI data, taking 660 ms and 289 ms, respectively.

AVPD was measured relative to the position of the valves during diastasis and is shown in Figure 12. Negative values are associated with the movement of the valves towards the apex. During the contraction of the atria, both the mitral and the tricuspid valve were pulled up by 5 mm and 7 mm, respectively. Peak AVDP was reached at end-systole with a displacement of $-10.7$ mm for the mitral valve and $-11.9$ mm for the tricuspid valve. However, before moving towards the apex, the valves moved 2–3 mm in the opposite direction, which is associated with the bulging of the atrioventricular valves into the atria. This behavior could not be observed in the MRI data, since the temporal resolution was too low to capture it. Compared to the MRI data, the tricuspid valve is displaced 5 mm less during systole. However, valve displacement during atrial contraction is of the same magnitude. In contrast, the mitral valve is displaced 1.5 mm more during atrial contraction compared to the MRI data and differs only by $-2.1$ mm during systole.

Starting from the last cycle of the control heart beat simulation, an ablation scar was introduced in the left atrium. It took three heart cycles for the volumes in the circulatory system to adapt to the new conditions and to reach a limit cycle with a stroke volume difference of less than 1 mL. As shown in Figure 13, the scar isolates the whole roof of the left atrium from electrical signals. Since all the pathways between the left and right atrium are still intact, there is no measurable difference with regard to the electrical activation in the rest of the atrium. Due to the isolated area in the left atrium, a significant proportion of tissue is not developing any active tension hence reducing the overall contractile function of the left atrium. Less blood is pumped out of the LA into the LV, reducing the LV end-diastolic volume by 2.5 mL to 200.4 mL. This can be observed by a decreased flow through the mitral valve during atrial contraction and consequently by an increase of the E/A ratio to 1.46. Right atrial and ventricular pumping efficiency are unaffected by the ablation scar. Additionally, a decreased AVPD during atrial contraction is observed.

5. Discussion

5.1. Bidirectional Coupling between the Mechanical and Electrophysiological Systems

Bidirectional coupling or strong coupling is required to simulate physiological behavior of the heart including mechano-electric feedback (MEF). Electrophysiological properties of cells are temporarily altered due to mechanical deformation. Currently, three major mechanisms are known to contribute to MEF: (1) the deformation of cardiac tissue during contraction and relaxation changes the electrical gradient $\sigma \nabla V_{\mathrm{m}}$; (2) the transmembrane voltage $V_{\mathrm{m}}$ is directly modified through stretch activated currents (SACs) $I_{\mathrm{SAC}}(V_{\mathrm{m}},{\gamma }_{\mathbf{f}})$; (3) the stretch in the cells changes the binding affinity of calcium to troponin C, thus altering the intracellular levels of calcium. For the numerical experiments described in Section 3.2, we did not consider the effect of deformation on the electrical gradient. Its effect should be marginal during normal sinus rhythm, since large deformations only occur at times when the electrical gradient is small, i.e., during the action potential plateau phase or repolarization. However, this effect can alter the electrical activity of the heart, especially at shorter cycle lengths or reentry, and should be considered in such scenarios [131]. While SACs have been used extensively in simulation studies, the underlying mechanisms that relate stretch to changes in ion channel conductance are not well understood and the majority of published models of SACs are not well constrained by experimental data [74,75,76,77,78,79]. Since we use human electro-mechanical cell models and the parameterization of SACs is out of the scope of this study, we decided to exclude them from our model until more data from humans are available. Nevertheless, both of the currently missing effects can be easily added to the presented framework.

As illustrated in Figure 7 and Figure 8, the cellular electro-mechanical combination of models proposed in this study yields physiologically accurate active tension in line with data from human tissue preparations. However, when used in a multi-scale environment, it still has room for improvement as reflected in the time course of LV volume (Figure 11). The global behavior of the LV does not entirely reflect the contractile dynamics observed in MRI data. During relaxation, the reduction in stress and increase in ventricular volume is too slow. As a result, no clear diastasis phase can be observed in our model. Furthermore, the rate of contraction during systole is higher and the model reaches end-systole more than twice as fast as the MRI data suggests. This behavior can be observed in other simulation studies using the Land et al. tension model [72,132] or other tension models based on sliding filament theory [31,133], suggesting that not all mechanisms of cardiac contraction are fully represented in the most recent models. A simulation study by Campbell suggests that an acceleration of the relaxation can be achieved by introducing interfilamentory movement resulting from compliance in the sarcomere [134]. Possibly, a similar effect could result from an increase in myocardial stiffness and should be investigated further.

A limitation of our electro-mechanical model is that the dependence on stretch rate ${\gamma }_{\mathbf{f}}$ in Equation (12) had to be neglected by setting ${\dot{\gamma }}_{\mathbf{f}}=0$ for all numerical experiments. Adding stretch rate dependence introduced strong unphysical oscillations in the multi-scale model resulting in an unstable numerical scheme. Regazzoni et al. [135] identified the source of these instabilities to be an inconsistent treatment of macroscopic and microscopic strains in staggered solution strategies of deformation in conjunction with an active stress model with a stretch rate dependence. This condition applies to the staggered approach used in this work (Section 2.5) and should be solved in the future either by implementing the proposed stabilization term in [135] or solving the whole system monolithically.

5.2. Circulatory System

As demonstrated in the numerical experiments, the closed-loop lumped-parameter model proposed in this manuscript is strongly coupled to the deformation problem and can reproduce major features of the human circulatory system. In the atria, the typical figure-of-eight shape is observed in the pressure-volume (PV) loops consisting of an A- and a V-loop corresponding to the atrial and ventricular contraction respectively. Additionally, a short increase in pressure in early systole known as the C-wave can be observed. Physiologically, this is the result of the atrioventricular valve’s cusps bulging into the atria during ventricular isovolumetric contraction. In the model, this observation was made possible by including physical representations of the valves in the mechanical mesh and was visible in the AVPD as well (Figure 12). In the ventricular PV-loops, the four major phases are present: (1) isovolumetric contraction; (2) ejection of blood from the ventricles; (3) isovolumetric relaxation; (4) passive filing with blood. Except for LV volume, no data related to the circulatory system were available from the healthy volunteer. Therefore, the parameterization of the lumped parameter model is motivated entirely by values obtained from literature and might not be fully representative to our specific case. The systolic pressures of the LV and the RV are too high for a healthy heart and are likely the result of the steep increase in stress due to the tension model. An optimized tension model could potentially address this issue and properly reflect systolic and diastolic behavior of the ventricles.

Using a closed-loop system makes it possible to study the influence of local changes in the cardiac tissue on the whole system. This was demonstrated by introducing an RFA scar into the left atrium (LA). The steady state of the circulatory system adapted to the decreased contractile function of the LA in a matter of 3 heart cycles and indicates that the general pumping efficiency of the ventricles is not heavily affected by the scar. Scar tissue typically shows an increased stiffness due to the higher collagen content compared to myocardial tissue. This could contribute to higher local stresses on the tissue in the atria and increasingly affect ventricular function. We did not consider this in our simulation setup and it should be included in further investigations of such scenarios. Studies like this can only be done with a whole heart model including the atria, the ventricles, and a circulatory system with a description of all necessary compartments. However, such a model is computationally very expensive and renders a patient-specific optimization of haemodynamic parameters challenging, especially since most of the parameters are unavailable from human measurements due to their invasive nature.

5.3. Numerical Considerations

The simulation framework proposed in this study was previously verified in N-version benchmarks that investigated the spatio-temporal convergence of the electrophysiological system [100] and the mechanical system [136], respectively.

For the solution of the monodomain equation, a small spatial resolution is required to avoid spatial undersampling artifacts leading to a reduction in conduction velocity or even conduction block. With an average spatial discretization of 0.6 mm of the electrophysiological domain ${\mathsf{\Omega }}_{\mathrm{EP}}$, we use a coarser mesh than all values tested in [100]. Recently, Woodworth et al. published in [137] a convergence study for the monodomain problem coupled to the cell model of ten Tusscher et al. [138]. The spatial discretization needed to be in the convergence region is out of scope of the technical requirements for the simulations of this work. However, to match prescribed velocities [132] for the electrical activation pattern the conductivity values are increased for the whole heart simulations. As shown in Figure 6 and Figure 13, the resulting activation sequence matches normal patterns observed in humans [107,109].

The benchmark of the mechanical deformation problem only considered a simplified approach, which cannot be directly translated to a whole heart simulation scenario. However, it includes important aspects of cardiac function, such as the active contraction of the heart muscle, pressure boundary conditions, and a complex distribution of fiber orientations. We extended this benchmark problem to include different spatial resolutions and finite elements of a different order. The results are presented in the Supplementary Material and show negligible deviations in the solution for a similar spatial discretization as the left ventricle in our whole heart model as long as quadratic elements are used. Quadratic elements additionally reduce the effect of volumetric locking that is introduced by the penalty formulation used to enforce quasi-incompressability. For linear tetrahedral elements, volumetric locking is a major concern and would require a substantially smaller spatial discretization to achieve the large deformations required in cardiac mechanics.

To the best of our knowledge, there is currently no extensive existence result nor convergence study available for the fully coupled electro-mechanical problem. Some results were first presented by Andreianov et al. [139] and an extensive proof is given in [140], where an active strain decomposition was assumed in contrast to the active stress formulation used in this work. Moreover, the equations of finite elasticity were linearized and assumed to be isotropic. The authors mention that the proof is extensible to anisotropic and polyconvex material formulations, for which there are some existence results ([63], Chapter 7). Unfortunately, the material formulation and boundary conditions used in this work do not satisfy the necessary assumptions.

First, ideas to establish a benchmark for the coupled system were proposed by Quarteroni et al. [11] and Santiago et al. [12]. Their results implicate that the change in conduction velocity due the spatial discretization errors in the solution of the monodomain equation are the most prominent and cause a temporal shift in mechanical activation. The magnitude of the deformation itself was not affected for the tested scenarios. This supports the approach applied in this study to use different spatial resolutions for the mechanical and the electrophysiological problems. However, using two domains requires projection techniques to transfer values between these domains. To solve this, linear interpolation using shape functions was used. Extrapolation of values is avoided by ensuring that the boundaries of the domains are identical.

6. Conclusions

In this work, we propose a framework for the fully coupled cardiac electro-mechanical problem with a detailed description of appropriate boundary conditions such as a lumped parameter model of the human circulatory system and a contact handling that replicates the effects of the tissue surrounding the heart.

We provide parameterizations of a fully coupled excitation contraction model for cells of the atrial and ventricular myocardium. Both the intracellular calcium transient and the tension development match data of human tissue preparations from literature. Furthermore, the validity of the model is demonstrated through a simulation on a personalized heart geometry created from magnetic resonance imaging data of a healthy volunteer. Coupling the 0D lumped parameter model of human circulation to all four chambers of the 3D electromechanical model enables a faithful reproduction of the major phases of the cardiac cycle as well as the characteristic figure of eight shape in the atrial PV-loops and flow patterns observed in clinical practice. Including the influence of the surrounding tissue into the model yields an atrioventricular plane displacement close to that observed in the MRI data. Finally, we demonstrate the potential of the model to include pathological changes and predict effects of therapeutic interventions by introducing an ablation lesion in the left atrium which led to changes in the activation and contraction of the left atrium. As a result of the changed loading conditions, the end-diastolic volume, and thus, the cardiac output of the left ventricle was reduced as well.

The adaptability of the framework allows comprehensive personalization of the model and paves the way towards digital twin modeling.