Towards a Computational Framework for Modeling the Impact of Aortic Coarctations Upon Left Ventricular Load

Computational fluid dynamics (CFD) models of blood flow in the left ventricle (LV) and aorta are important tools for analyzing the mechanistic links between myocardial deformation and flow patterns. Typically, the use of image-based kinematic CFD models prevails in applications such as predicting the acute response to interventions which alter LV afterload conditions. However, such models are limited in their ability to analyze any impacts upon LV load or key biomarkers known to be implicated in driving remodeling processes as LV function is not accounted for in a mechanistic sense. This study addresses these limitations by reporting on progress made toward a novel electro-mechano-fluidic (EMF) model that represents the entire physics of LV electromechanics (EM) based on first principles. A biophysically detailed finite element (FE) model of LV EM was coupled with a FE-based CFD solver for moving domains using an arbitrary Eulerian-Lagrangian (ALE) formulation. Two clinical cases of patients suffering from aortic coarctations (CoA) were built and parameterized based on clinical data under pre-treatment conditions. For one patient case simulations under post-treatment conditions after geometric repair of CoA by a virtual stenting procedure were compared against pre-treatment results. Numerical stability of the approach was demonstrated by analyzing mesh quality and solver performance under the significantly large deformations of the LV blood pool. Further, computational tractability and compatibility with clinical time scales were investigated by performing strong scaling benchmarks up to 1536 compute cores. The overall cost of the entire workflow for building, fitting and executing EMF simulations was comparable to those reported for image-based kinematic models, suggesting that EMF models show potential of evolving into a viable clinical research tool.

1. Introduction

CFD models of blood flow in the LV and aorta are important tools for analyzing the mechanistic links between myocardial deformation and flow patterns. Typically, such models are either driven by prescribed flow profiles measured in the LV outflow tract or the aortic root (Goubergrits et al., 2013; Ralovich et al., 2015; Andersson et al., 2017), or by image-based kinematic models (Doenst et al., 2009; Schenkel et al., 2009; Mihalef et al., 2011; Seo et al., 2013; Chnafa et al., 2014; Su et al., 2016) built from segmentation of 4D medical imaging datasets. While such models have proven to be valuable for analyzing the hemodynamic status quo of a patient or for predicting changes in hemodynamics in the aorta secondary to intervention such as aortic valve repair (Kelm et al., 2017) or stenting of a coarctation (Goubergrits et al., 2015), they are inherently limited in their ability to assess cardiac function as the biophysics driving myocardial activation and deformation is not taken into consideration in the model formulation. EMF models that capture the entire physics of a heartbeat based on first principles show promise to overcome this limitation (Crozier et al., 2016a) by rendering feasible the assessment of all essential myocardial parameters, which are known to be key factors driving ventricular remodeling and disease progression. Thus EMF models may offer, in principal, the potential of predicting longer term outcomes beyond changes in the acute response to therapies.

However, due to a number of factors such as the inherent complexity of multiphysics models, the large-scale motion and complex deformation of the myocardial walls as well as the significant computational burden, these models pose substantial methodological challenges. For LV EMF models and similar applications, methods to overcome the problem of large-scale deformations can be roughly classified into two categories: ALE formulations using a moving fluid mesh (Tang et al., 2008, 2010; Nordsletten et al., 2011; Vázquez et al., 2015; de Vecchi et al., 2016) and immersed boundary (IB) methods (Vigmond et al., 2008a; Seo and Mittal, 2013; Choi et al., 2015). While ALE formulations often rely on severe simplifications or automatic remeshing strategies (Long et al., 2013), IB methods are more versatile as the moving wall of the ventricle is not explicitly tracked. However, IBs and all related non-boundary-fitting methods have a reduced accuracy for the solution near the fluid-solid structure interface due to interpolation errors, pose severe challenges on the implementation, and additional degrees of freedom have to be introduced on interface cut elements, which all contributes to significantly higher computational costs (van Loon et al., 2007).

In this study, we report on the progress made toward a novel EMF model of the human LV that is entirely based on first principles and that copes with significantly large defomations, i.e., ejection fractions (EFs) beyond 60%, without requiring remeshing or IB principles. Validated in silico models taken from a recent clinical modeling study where a cohort of in silico EM LV and aorta models of patients suffering from aortic valve disease (AVD) and/or CoA (Augustin et al., 2016a) were built, served as kinematic driver to a computational model of hemodynamics in the LV cavity and aorta. A hybrid two stage modeling approach was adopted with regard to hemodynamics. First, the afterload imposed by the circulatory system onto the LV was represented by a lumped model of afterload and coupled to an EM model of LV and aorta to compute LV kinematics. Subsequently a full-blown CFD model with moving domain boundaries based on an ALE formulation was unidirectionally or weakly coupled to the EM model using the kinematics of its endocardial surface as input. We show validation results for two selected clinical CoA cases under pre-treatment conditions and compare pre-treatment and post-treatment simulation results for one patient case in which the CoA was geometrically repaired by a virtual stenting procedure. Further, we demonstrate numerical feasibility of the implemented approach by analyzing changes in mesh quality and its impact upon solver performance under the significantly large deformations of the LV blood pool mesh and also provide strong scaling benchmarking results for a range of 96–1,536 compute cores. The overall cost of the entire workflow for building, fitting and execution of EMF simulations is ≈ 48 h which is comparable to plain image-based kinematic driver models (Mittal et al., 2016).

2. Methods

The methodology to develop a coupled model of cardiac and cardiovascular hemodynamics based on an ALE formulation is structured as follows.

 i) We begin in section 2.1 by describing MRI data acquisition and anatomical FE model generation of the LV and aorta for two patients suffering from CoA.

ii) Then, a brief summary of all model components is given comprising an electrophysiology (EP) model to drive electrical activation and repolarization (section 2.2.1); an EM model describing passive biomechanics as well as the generation of active stresses (section 2.2.2); afterload models to provide appropriate boundary conditions on the LV endocardium during the ejection phase (section 2.2.3); and a CFD model with moving domain boundaries representing blood flow in the LV and aorta during ejection. The EM and CFD model are weakly coupled in a forward fluid structure interaction (FSI) framework, where the EM model is used as a kinematic driver to move the fluid domain (section 2.3).

iii) The solution procedure and software implementation details are outlined in section 2.4.

iv) Finally, procedures implemented for the patient-specific parameterization of the major model components is described in section 2.5.

2.1. Clinical Data Acquisition and Model Generation

Hemodynamic data of two patients with clinical indication for catheterization due to CoA—all preceding a cardiac magnetic resonance study—were acquired before and after CoA treatment by stent implant, see Table 1. CoA treatment indicators included an echocardiographic measured, peak systolic pressure gradient across the stenotic region of > 20 mmHg and/or arterial hypertension. The study was approved by the institutional research ethics committee following the ethical guidelines of the 1975 Declaration of Helsinki. Written informed consent was obtained from the participants' guardians. Acquired data are summarized in Table 1.

TABLE 1

Table 1. CoA patient characteristics from MRI and invasive catheter pressure recordings including end-diastolic volume (EDV), end-systolic volume (ESV), stroke volume (SV), ejection fraction (EF), heart rate (HR), cardiac output (CO), diastolic and systolic pressures recorded in the aorta or estimated from cuff measurements (P_ao/cuff,dia and P_{ao/cuff, sys}), mean arterial pressure (MAP) computed from pressure recorded invasively in the aorta or estimated from P_cuff,dia and P_cuff,sys, and aortic valve open pressure P_open determined from invasive pressure recordings.

2.1.1. MRI Acquisition and Post Processing

MR imaging was done with a whole body 1.5 Tesla MR scanner Achieva R 3.2.2.0 using a five-element cardiac phased-array coil (Philips Medical System, Best, Netherlands). Three MRI sequences were used further in our study: (i) flow-sensitive four-dimensional (4D) velocity-encoded magnetic resonance imaging (4D VEC-MRI), (ii) three-dimensional (3D) anatomical imaging of the whole heart (3DWH) during diastasis, and (iii) 4D gapless short axis Cine MRI.

4D VEC-MRI of the thorax was performed using an anisotropic 4D segmented k-space phase contrast gradient echo sequence. Retrospective electrocardiographic gating without navigator gating of respiratory motion in order to minimize acquisition time was used. Sequence parameters were: acquired voxel 2.5 × 2.5 × 2.5 mm; reconstructed voxel 1.7 × 1.7 × 2.5 mm; repetition time 3.5 ms; echo time 2.2 ms; flip angle 5°; 25 reconstructed cardiac phases; number of signal averages 1; High velocity encoding (3–6 m/s) in all three directions was used in order to avoid phase wraps in the presence of coarctation and associated secondary flow. Flow measurements were completed with automatic correction of concomitant phase errors. Postprocessing for analysis of flow rates across the aortic valve was carried out with GTFlow 1.6.8 software¹ (Gyrotools, Zurich, Switzerland).

The 3DWH exemplary sequence parameters were: acquired voxel 0.66 × 0.66 × 3.2 mm; reconstructed voxel 0.66 × 0.66 × 1.6 mm; repetition time 4.0 ms; echo time 2.0 ms; flip angle 90°; and number of signal averages 3.

Short axes Cine imaging data were acquired with sequence parameters: 16 slices, with an acquisition resolution of 0.86 × 0.86 × 6.0 mm, repetition time 4.24 ms, echo time 2.12 ms, flip angle 60° and 25 automatically reconstructed cardiac phases which were used to determine LV volume traces. The non-compact myocardium as well as papillary muscles were counted toward blood pool volume.

MRI based pressure mapping allowing to assess non-invasively the relative pressures in a vessel by solving Pressure Poisson equation (PPE) was done with MevisFlow². Briefly, the PPE can be derived from the Navier–Stokes equations by taking the divergence of the momentum Equation (26), see Gresho and Sani (1987) and Krittian et al. (2012) for more details. For more details we refer to Krittian et al. (2012). The processing and analysis pipeline of the pressure mapping consists of the following four steps.

 i) Semi-automatic segmentation (labeling) of the aortic domain from 3DWH data generating 3D mask of the aorta.

ii) Background phase correction and phase-unwrapping of the 4D VEC-MRI data and generation of a sequence of volumetric velocity vector fields.

iii) Coarse semi-automatic segmentation of the aorta based on magnitude and phase contrast of the 4D VEC-MRI data and registration with 3DWH based mask of the aorta.

iv) Solving (PPE) at each time step having 4D VEC-MRI data as input. Furthermore, a 5 % mask size reduction is applied in order to avoid numerical inconsistencies close to the vessel wall as suggested earlier (Meier et al., 2010).

Relative pressure maps are represented with zero pressure located at the center of the CoA (narrowest location). 3D mask based on 3DWH data was used due to its better spatial resolution compared to 4D VEC-MRI data. Correction of velocity data (step ii) was done in order to minimize noise and aliasing artifacts originating from multiple sources.

2.1.2. Invasive Catheter Recordings

During catheterization, pressure was recorded over the cardiac cycle in the ascending aorta and the LV before treatment and repeated in the ascending aorta after an interventional treatment procedure was performed. Pressures were recorded simultaneously at three predefined locations (LV, ascending aorta, and descending aorta) and the femoral artery during catheterization. Patients were sedated by intravenous administration of a bolus of midazolam (0.1–0.2 mg/kg, max. 5 mg), followed by a bolus of propofol (1–2 mg/kg, as needed) and continuous infusion of propofol (1–2 mg/kg, as needed). Pressure measurements were taken with senior cardiologists present. Pigtail catheters (Cordis, Warren, NJ, USA) of 5-6F were connected to pressure transducers (Becton-Dickinson, Franklin Lakes, NJ, USA). Routinely, patients received balloon angioplasty with or without additional placement of a stent in order to treat a given stenosis by removing the narrowing of the vessel and thus the pressure gradient. To reduce duration of catheterization, pressures were measured post-treatment only in the ascending aorta. The Schwarzer hemodynamic analysis system (Schwarzer, Heilsbronn, Germany) was used to amplify, acquire, and analyze pressure signals.

2.1.3. Anatomical FE Model Generation

Multi-label segmentation of the LV myocardium, LV blood pool, left atrium (LA) and aortic cavities was done at the DHZB using 3DWH data and the ZIB Amira software³ (Stalling et al., 2005). The segmentations were smoothed and upsampled to a 0.1 mm isotropic resolution using a variational smoothing method (Crozier et al., 2016a). The resulting high resolution multi-label segmentation was meshed using CGAL⁴ (The CGAL Project, 2017), giving a global mesh ${\Omega }_{\text{s},\text{total}}^0$ consisting of tetrahedral elements. Here, (•)⁰ denotes the mechanical reference configuration at end-diastolic pressure. The mesh was subdivided into various subdomains corresponding to predefined labels which are summarized in Equation (3). We write

$\begin{array}{ll}{\Omega }_{\text{s},\text{total}}^0={\displaystyle \underset{i\in I}{\cup }{\Omega }_{\text{s},i}^0},& (1)\end{array}$

with the index set

$\begin{array}{ll}I:=\{\text{lv},\text{ao},\text{cushion},\text{av},\text{mv},\text{lvbp},\text{aobp}\},& (2)\end{array}$

see Figures 1E–G for illustration. The elements in the index set are abbreviations for the following labels

$\begin{array}{l}\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}lv }\leftrightarrow \text{Myocardium},\\ \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}ao }\leftrightarrow \text{Aortic wall},\\ \text{cushion }\leftrightarrow \text{Elastic cushion},\\ \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}av }\leftrightarrow \text{Aortic valve},\\ \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}mv }\leftrightarrow \text{Mitral valve},\\ \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}lvbp }\leftrightarrow \text{Left ventricular bloodpool},\\ \text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}aobp }\leftrightarrow \text{Aortic bloodpool}.& (3)\end{array}$

With this, we define the following submeshes

$\begin{array}{ll}{\Omega }_{\text{s}}^0≔{\Omega }_{\text{s},\text{total}}^0\backslash \left({\Omega }_{\text{s},\text{lbvp}}^0\cup {\Omega }_{\text{s},\text{aobp}}^0\right),& (4)\end{array}$

$\begin{array}{ll}{\Omega }_{\text{s},\text{bp}}^0={\tilde{\Omega }}_f^0≔{\Omega }_{\text{s},\text{av}}^0\cup {\Omega }_{\text{s},\text{lvbp}}^0\cup {\Omega }_{\text{s},\text{aobp}}^0,& (5)\end{array}$

where ${\Omega }_{\text{s}}^0$ is the solid domain and ${\Omega }_{\text{s},\text{bp}}^0$ is the unsmoothed blood pool domain used for extracting a smoothed CFD mesh, see Figures 1E,F. For later use, we define the following surfaces

$\begin{array}{ll}{\Gamma }_{\text{s},\text{N}}^0≔\partial \left(\left({\Omega }_{s,\text{lv}}^0\cup {\Omega }_{\text{s},\text{av}}^0\cup {\Omega }_{\text{s},\text{mv}}^0\right)\cap {\Omega }_{\text{s},\text{lvbp}}^0\right),& (6)\end{array}$

$\begin{array}{ll}{\Gamma }_{\text{s},\text{H}}^0≔\partial {\Omega }_{\text{s}}^0\backslash \left({\Gamma }_{\text{s},\text{N}}^0\cup {\Gamma }_{\text{s},\text{D}}^0\right),& (7)\end{array}$

$\begin{array}{ll}{\Gamma }_{\text{s},\text{bp}}^0≔\partial {\Omega }_{\text{s},\text{bp}}^0\backslash {\Gamma }_{\text{s},\text{D}}^0,& (8)\end{array}$

where ${\Gamma }_{\text{s},\text{D}}^0$ denote the cutoff faces as indicated by blue lines in Figure 1; ${\Gamma }_{\text{s},\text{N}}^0$ are surfaces subject to pressure; and ${\Gamma }_{\text{s},\text{H}}^0$ are surfaces with homogeneous Neumann boundary conditions. In order to avoid numerical difficulties with non-smooth, jagged boundaries, the surface of the mechanical blood pool domain ${\Gamma }_{\text{s},\text{bp}}^0$ was extracted and smoothed using the VMTK toolbox⁵ (Antiga et al., 2008). The smoothed surface, ${\Gamma }_{\text{f},\text{wall}}^0$, was used to define the boundary of the fluid domain reference configuration, ${\Omega }_{\text{f}}^0$, for volumetric FE meshing using ANSYS ICEM CFD⁶. Refined boundary layers were included in this process to better resolve sharp gradients in the vicinity of ${\Gamma }_{\text{f},\text{wall}}^0$ occurring during simulation of hemodynamics. The various processing stages for building EM and CFD models are illustrated in Figures 1, 4, respectively.

FIGURE 1

Figure 1. Mechanics model generation: Starting from a patient specific MRI scan (A) a segmentation was performed (B) which was then upsampled and smoothed (C). Myocardial fibers were generated in the tissue according to Bayer et al. (2012) (D). A labeled FE geometry ${\Omega }_{\text{s},\text{total}}^0$ including the blood pool was generated (G). The geometry has been sliced to reveal the blood pool and valves and has been color coded according to the labels defined in Equation (3). Boundaries ${\Gamma }_{\text{s},\text{D}}^0$ used for prescribing homogeneous Dirichlet boundary conditions are sketched as blue curves. From this mesh the EM submesh ${\Omega }_{\text{s}}^0$ (E) and the unsmoothed blood pool (F) were extracted. Boundary ${\Gamma }_{\text{s},\text{N}}^0$ was used to prescribe pressure boundary conditions inside the LV and ${\Gamma }_{\text{s},\text{bp}}^0$ is the surface of the blood pool.

2.2. Electromechanical Model

2.2.1. Electrophysiology of the LV

A recently developed reaction-eikonal (R-E) model (Neic et al., 2017) was employed to generate electrical activation sequences which serve as a trigger for active stress generation in cardiac tissue. The hybrid R-E model combines a standard reaction-diffusion (R-D) model based on the monodomain equation with an eikonal model. Briefly, the eikonal equation is given as

$\begin{array}{ll}\{\begin{array}{c}\sqrt{{\nabla }_{\text{X}}{t}_{\text{a}}^{\top }\text{V}{\nabla }_{\text{X}}{t}_{\text{a}}}=1\text{\hspace{1em}in }{\Omega }_{\text{s},\text{lv}}^0,\\ t_{\text{a}}=t_0\text{ on }{\Gamma }_{\text{s},\ast }^0,\end{array}& (9)\end{array}$

where (∇_X) is the gradient with respect to the end-diastolic reference configuration ${\Omega }_{\text{s},\text{lv}}^0$; t_a is a positive function describing the wavefront arrival time at location $\mathbf{\text{X}}\in {\Omega }_{\text{s},\text{lv}}^0$; and t₀ are initial activations at locations ${\Gamma }_{\text{s},*}^0\subseteq {\Gamma }_{\text{s},\text{N}}^0$. The symmetric positive definite 3 × 3 tensor V(X) holds the squared velocities (v_f(X), v_s(X), v_n(X)) associated to the tissue's eigenaxes, referred to as fiber, f₀, sheet, s₀, and sheet normal, n₀, orientations. The arrival time function t_a(X) was subsequently used in a modified monodomain R-D model given as

$\begin{array}{ll}\beta C_{\text{m}}\frac{\partial V_{\text{m}}}{\partial t}={\nabla }_{\text{X}}·{\sigma }_{\text{i}}{\nabla }_{\text{X}}{V}_{\text{m}}+I_{\text{foot}}-\beta I_{\text{ion}},& (10)\end{array}$

where an arrival time dependent foot current, I_foot(t_a), was added which is designed to mimic subthreshold electrotonic currents to produce a physiological foot of the action potential. The key advantage of the R-E model is its ability to compute activation sequences at much coarser spatial resolutions that are not afflicted by the spatial undersampling artifacts leading to conduction slowing or even numerical conduction block as it is observed in standard R-D models. Ventricular EP was represented by the tenTusscher–Noble–Noble–Panfilov model of the human ventricular myocyte (ten Tusscher et al., 2004). As indicated in Equations (9, 10), activation sequences and electrical source distribution in the LV were computed in its end-diastolic configuration ${\Omega }_{\text{s},\text{lv}}^0$, that is, any effects of deformation upon electrotonic currents remained unaccounted for.

2.2.2. Active and Passive Mechanics in the LV and Aorta

The deformation of the heart is governed by imposed external loads such as pressure in the cavities or from surrounding tissue and active stresses intrinsically generated during contraction. Tissue properties of the LV myocardium and the aorta are characterized as a hyperelastic, nearly incompressible, anisotropic material with a non-linear stress-strain relationship. Mechanical deformation was described by Cauchy's equation of motion under stationary equilibrium assumptions leading to a quasi-static boundary value problem

$\begin{array}{ll}-{\nabla }_{\text{X}}·\text{FS}({\text{d}}_{\text{s}},t)=0\text{\hspace{1em}in }{\Omega }_{\text{s}}^0,& (11)\end{array}$

for t ∈ [0, T], where d_s is the unknown displacement; F is the deformation gradient; S is the second Piola–Kirchhoff stress tensor; and (∇_X ·) denotes the divergence operator in the Lagrange reference configuration. Homogeneous Dirichlet boundary conditions

$\begin{array}{ll}{\text{d}}_{\text{s}}=0\text{\hspace{1em}on\hspace{1em}}{\Gamma }_{\text{s},\text{D}}^0,& (12)\end{array}$

homogeneous Neumann boundary conditions

$\begin{array}{ll}\begin{array}{lll}\mathit{F}\mathit{S}({\overrightarrow{d}}_{\text{s}},t)=\overrightarrow{0}\hfill & \hfill & \text{on }{\Gamma }_{\text{s},\text{H}}^0,\hfill \end{array}& (13)\end{array}$

and inhomogeneous Neumann boundary conditions

$\begin{array}{ll}\text{FS}({\text{d}}_{\text{s}},t){\text{n}}_{\text{s,0}}=p(t)J{\text{F}}^{-\top }({\text{d}}_{\text{s}},t){\text{n}}_{\text{s},\text{0}}\text{\hspace{1em}on\hspace{1em}}{\Gamma }_{\text{s},\text{N}}^0& (14)\end{array}$

were imposed, where n_s,0 is the outward unit normal vector; p(t) is the pressure; and J = det F. For sake of clarity, boundary conditions are illustrated in Figure 1C.

The total stress S was additively decomposed according to

$\begin{array}{ll}\text{S}={\text{ S}}_{\text{pas}}+{\text{ S}}_{\text{act}},& (15)\end{array}$

where S_pas and S_act refer to the passive and active stresses, respectively. Passive stresses were modeled based on the constitutive equation

$\begin{array}{ll}{\text{S}}_{\text{pas}}=2\frac{\partial \Psi (\text{C})}{\partial \text{C}}& (16)\end{array}$

given a hyper-elastic strain-energy function Ψ and the right Cauchy–Green strain tensor C = F^⊤ F. Two different strain-energy functions were used for characterizing passive mechanical behavior in the LV and the aorta. In the LV, where the underlying mesh ${\Omega }_{\text{s},\text{lv}}^0$ and fiber orientations (f₀, s₀, n₀) are the same as for the EP model, section 2.2.1, the transversely isotropic constitutive relation

$\begin{array}{ll}{\Psi }_{\text{Guc}}(\text{C})=\frac{\kappa }{2}{\left(\log J\right)}^2+\frac{C_{\text{Guc}}}{2}\left[\exp ({Q})-1\right].& (17)\end{array}$

by Guccione et al. (1995) was employed. Here, the term in the exponent is

$\begin{array}{l}{Q}=b_{\text{f}}{\left({\text{f}}_0·\overline{\text{E}}{\text{f}}_0\right)}^2+b_{\text{t}}\left[{({\text{s}}_0·{\overline{\text{E}}\text{s}}_0)}^2+{({\text{n}}_0·{\overline{\text{E}}\text{n}}_0)}^2+2{({\text{s}}_0·{\overline{\text{E}}\text{n}}_0)}^2\right]\\ +2b_{\text{fs}}\left[{({\text{f}}_0·{\overline{\text{E}}\text{s}}_0)}^2+{({\text{f}}_0·{\overline{\text{E}}\text{n}}_0)}^2\right]& (18)\end{array}$

and $\overline{\mathbf{\text{E}}}=\frac{1}{2}(\overline{\mathit{\text{C}}}-\mathbf{\text{I}})$ is the modified isochoric Green–Lagrange strain tensor, where $\overline{\text{C}}:=J^{-2/3\text{C}}$. Default values of b_f = 18.48, b_t = 3.58, and b_fs = 1.627 were used. The parameter C_Guc was varied for the different cases, see Table 2. In the aorta ${\Omega }_{\text{s},\text{ao}}^0$, unlike in previous studies (Augustin et al., 2014), we refrained from assigning fiber structures, since our efforts were primarily focused on modeling the biomechanics of the LV and, to a lesser degree, the aorta. Thus, in absence of information on structural anisotropy, an isotropic model due to Demiray (1972) was used

$\begin{array}{ll}{\Psi }_{\text{Dem}}(\text{C}):=\frac{\kappa }{2}{\left(\log J\right)}^2+\frac{a}{2b}\left\{\exp \left[b(\text{tr}(\overline{\text{C}})-3)\right]-1\right\}.& (19)\end{array}$

The parameter $\tilde{C}=\frac{a}{2b}$ was chosen such that $\tilde{C}$ = 3,000 kPa in the aorta, $\tilde{C}$ = 30,000 kPa for valves, and $\tilde{C}=300kPa$ for the elastic cushion. The bulk modulus κ, which serves as a penalty parameter to enforce nearly incompressible material behavior, was chosen as κ = 650 kPa in both Equations (17, 19). For the elastic cushion a value of κ = 100 kPa was used.

TABLE 2

Table 2. Fitted parameters for EM Model.

A simplified phenomenological contractile model was used to represent active stress generation (Niederer S. A. et al., 2011). Owing to its small number of parameters and its direct relation to clinically measurable quantities such as peak pressure, p_lv, and the maximum rate of rise of pressure, dp_lv/dt_max, this model is fairly easy to fit and thus very suitable for being used in clinical EM modeling studies. Briefly, the active stress transient is given by

$\begin{array}{l}{S}_{\text{a}}(t,\lambda )=S_{\text{peak}}\varphi (\lambda ){\tanh }^2\left(\frac{t_{\text{s}}}{{\tau }_{\text{c}}}\right){\tanh }^2\left(\frac{t_{\text{dur}}-t_{\text{s}}}{{\tau }_{\text{r}}}\right),\\ \text{\hspace{1em}\hspace{1em}for }0<t_{\text{s}}<t_{\text{dur}},& (20)\end{array}$

with

$\begin{array}{ll}\varphi =\tanh (\text{ld}(\lambda -{\lambda }_0)),{\tau }_{\text{c}}={\tau }_{{\text{c}}_0}+{\text{ld}}_{\text{up}}(1-\varphi ),t_{\text{s}}=t-t_{\text{a}}-t_{\text{emd}}& (21)\end{array}$

and t_s is the onset of contraction; ϕ(λ) is a non-linear length-dependent function in which λ is the fiber stretch and λ₀ is the lower limit of fiber stretch below which no further active tension is generated; t_a is the local activation time from Equation (9); t_emd is the EM delay between the onsets of electrical depolarization and active stress generation; S_peak is the peak isometric tension; t_dur is the duration of active stress transient; τ_c is time constant of contraction; τ_c0 is the baseline time constant of contraction; ld_up is the length-dependence of τ_c; τ_r is the time constant of relaxation; and ld is the degree of length dependence. Thus, active stresses in this simplified model are only length-dependent, but dependence on fiber velocity, $\stackrel{.}{\lambda }$, is ignored. Unlinke in previous studies (Niederer S. A. et al., 2011) we set the nonlinear length-dependent function ϕ(λ) = 1 for the whole simulation. The active stress tensor in the reference configuration ${\Omega }_{\text{s},\text{lv}}^0$ induced in fiber direction f₀ is defined as

$\begin{array}{ll}{\text{S}}_{\text{a}}=S_{\text{a}}{\left({\text{f}}_0·{\text{Cf}}_0\right)}^{-1}{\text{f}}_0\otimes {\text{f}}_0,& (22)\end{array}$

with S_a defined in Equation (20). This active stress involves a scaling by ${\lambda }^2={\mathbf{\text{f}}}_0·\mathbf{\text{C}}{\mathbf{\text{f}}}_0$, see Pathmanathan and Whiteley (2009) for details.

2.2.3. Mechanical and Hemodynamic Afterload Models

Hydrostatic pressures in the LV, p_lv, and the proximal aorta, p_ao, were modeled using a 3-element Windkessel model (Westerhof et al., 1971), and the system of PDEs (11) was linked to this lumped model of the arterial system, see Figure 2. The models were coupled by a diode (aortic valve) which opens at the end of the isovolumetric contraction (IVC) phase when the pressure in the LV cavity, p_lv, exceeds the pressure in the proximal aorta, p_ao, and closes at the end of ejection when p_lv drops below p_ao and the flow q_lv starts to reverse. In its open state the aortic valve was modeled as a linear resistor, R_av, in series with the characteristic impedance of the aorta, Z_c. During ejection, the pressure in the LV was then computed by the Windkessel equation

$\begin{array}{ll}\frac{\text{d}{p}_{\text{lv}}}{\text{d}t}=\frac{1}{C}\left(1+\frac{Z_{\text{c}}+R_{\text{av}}}{R}\right)q_{\text{lv}}+(Z_{\text{c}}+R_{\text{av}})\frac{\text{d}{q}_{\text{lv}}}{\text{d}t}-\frac{1}{RC}{p}_{\text{lv}},& (23)\end{array}$

which predicts the rate of change of pressure in the LV as a function of flow q_lv out of the LV into the aorta. The resistor R represents peripheral arterial resistance placed in parallel with a capacitor C, representing vascular compliance.

FIGURE 2

Figure 2. Lumped circuit representation of the coupled EM PDE model of the LV with the cardiovascular system. The time-varying compliance of the LV is represented as a PDE model which was coupled through the aortic valve (R_av) to a 3-element Windkessel model representing aortic impedance, Z_c, and peripheral arterial compliance, C, and resistance, R, during ejection, and through the mitral valve (R_mv) to a constant pressure p_la in the left atrium during filling. Negative flows −q_la and −q_lv mean the respective cavity is ejecting, while positive flow means cavity is being filled.

A similar form of Equation (23) was also used to estimate the pressure in the aorta, p_ao. In this case, there is no additional resistance due to an outlet valve and hence R_av is omitted. Balancing of the PDE (11) and the ODE (23) was achieved by recasting Equation (11) as a saddle point problem, see Gurev et al. (2015) and Hirschvogel et al. (2017).

For CFD simulations, hydrostatic pressures at artificial aortic fluid outlets, were modeled using a similar 3-element Windkessel model as in Equation (23) that was rewritten in the form of the following differential algebraic equations for outlet i

$\begin{array}{ll}{C}_i\frac{\text{d}{p}_{\text{d},i}}{\text{d}t}+\frac{p_{\text{d},i}}{R_i}=q_i,& (24)\end{array}$

$\begin{array}{ll}{p}_{\text{wk},i}=Z_i{q}_i+p_{\text{d},i},& (25)\end{array}$

see Fouchet-Incaux (2014) and Bertoglio et al. (2017) for more details. During ejection the Windkessel pressure p_wk at an outlet was then applied as an outflow boundary condition for the fluid flow model, see section 2.5.5. In Equations (24, 25), C_i represents compliance, Z_i impedence, and R_i resistance of the peripheral arteries for the respective aortic outlet and q_i denotes the flux through this outlet. Fitting of the parameters involved will be discussed in section 2.5.5.

2.3. Fluid Flow Model

Human blood in larger vessels such as the LV or the aorta complies with the assumptions of an incompressible, isothermal, Newtonian and single-phase liquid (Nichols et al., 2011). Let ${\Omega }_{\text{f}}⊊{ℝ}^3$ denote the fluid domain, then the evolution of flow is governed by the incompressible Navier–Stokes equations

$\begin{array}{ll}\begin{array}{ll}{\rho }_{\text{f}}(\frac{\partial }{\partial t}{u}_{\text{f}}+u_{\text{f}}\cdot {\nabla }_x{u}_{\text{f}})-{\nabla }_x\cdot {\sigma }_{\text{f}}(u_{\text{f}},p_{\text{f}})=0\hfill & \text{in }{\Omega }_{\text{f}}\hfill \end{array},& (26)\end{array}$

$\begin{array}{ll}\begin{array}{ll}{\nabla }_x\cdot u_{\text{f}}=0\hfill & \text{in\hspace{0.17em}}\hfill \end{array}{\Omega }_{\text{f}},& (27)\end{array}$

$\begin{array}{rc}{u}_{\text{f}}=0\text{\hspace{1em}\hspace{1em}on\hspace{0.17em}\hspace{0.17em}}{\Gamma }_{\text{noslip}},& (28)\end{array}$

$\begin{array}{rc}{u}_{\text{f}}=g_{\text{f}}\text{ on }{\Gamma }_{\text{inflow}}\text{,}& (29)\end{array}$

$\begin{array}{ll}\begin{array}{ll}{\sigma }_{\text{f}}{n}_{\text{f}}-{\rho }_{\text{f}}\beta \left(u_{\text{f}}\cdot n_{\text{f}}\right)\_u_{\text{f}}=h_{\text{f}}\hfill & \text{on\hspace{0.17em}}{\Gamma }_{\text{outflow}}\hfill \end{array},& (30)\end{array}$

$\begin{array}{ll}{u_{\text{f}}|}_{t=0}=u_0,& (31)\end{array}$

where u_f denotes fluid velocity in m/s; p_f is fluid pressure in Pa; ρ_f is the density of blood, given as 1.060 kg/m³; σ_f is the fluid stress tensor in, Pa, defined as $-p_{\text{f}}\text{I}+{\mu }_{\text{f}}({\nabla }_{\mathbf{\text{x}}}{\mathbf{\text{u}}}_{\text{f}}+{\nabla }_{\mathbf{\text{x}}}{\mathbf{\text{u}}}_{\text{f}}^{\top })$, with dynamic viscosity of blood μ_f given as 0.004 Pa s; g_f, in m/s is a velocity inlet; p_wk, in Pa, is the Windkessel pressure solution to Equations (24, 25); u₀, in m/s, refers to the initial condition; n_f is the outward unit normal of the fluid domain; and (∇_x) is the gradient and (∇_x·) is the divergence operator in the fluid domain Ω_f. The sets Γ_noslip, Γ_inflow, and Γ_outflow denote the complementary subsets of Γ_f: = ∂Ω_f and we assume that |Γ_outflow| > 0. Note that Equation (29) is given only for the sake of completeness but was not used in this study, as the inflow of blood into the aorta is driven by the motion of the LV thus avoiding the need for prescribing an inflow profile as it is necessary in models which consider the aorta in isolation. For p_wk ≡ 0, boundary condition Equation (30) is referred to as directional do-nothing boundary condition, see Esmaily Moghadam et al. (2011) and Braack et al. (2014), and the term

$\begin{array}{ll}\left(u_{\text{f}}\cdot n_{\text{f}}\right)\_:=\frac{1}{2}\left(u_{\text{f}}\cdot n_{\text{f}}-|u_{\text{f}}\cdot n_{\text{f}}\text{|}\right)& (32)\end{array}$

is added for backflow stabilization. A value of $\beta >\frac{1}{2}$ was assumed to guarantee stability of the system. However, in practical applications values of $\beta \le \frac{1}{2}$ were also used without causing numerical issues, see Esmaily Moghadam et al. (2011). In presence of multiple outlets outflow boundary conditions as given in Equation (30) were prescribed at each of the outlets.

2.3.1. Extension to Moving Geometries

For time-dependent fluid domains, i.e., ${\Omega }_{\mathrm{\text{f}}}={\Omega }_{\mathrm{\text{f}}}^t$, Equations (26–31) need to be modified to account for the domain movement. This requires the linking of the equations governing fluid dynamics—posed in an Eulerian coordinate frame—with the structural mechanics equations—posed in a Lagrangian reference frame. This is achieved by using the ALE formulation which combines both Lagrangian and Eulerian formulation in a generalized description, see Bazilevs et al. (2013, section 1.3) and Hirt et al. (1974). Similar to structural mechanics, a reference fluid configuration ${\Omega }_{\mathrm{\text{f}}}^0⊊{ℝ}^3$ is used which we identify with the mesh been generated at end-diastolic state, see section 2.1.3. The coordinate system of the Eulerian frame is denoted by x and the reference coordinate system is denoted by X. Their relation is given by the ALE mapping x = X + d_f(t, X). Here, d_f(t, X) refers to an arbitrary, not necessarily physical, displacement of points to track the deformation of the fluid domain. Using this ALE mapping the time-dependent moving fluid domain is represented as

$\begin{array}{ll}{\Omega }_{\text{f}}^t≔\left\{x:x=X+d_{\text{f}}(t,X),\forall X\in {\Omega }_{\text{f}}^0\right\}.& (33)\end{array}$

Further, we define the fluid domain velocity w_f as

$\begin{array}{ll}{w}_{\text{f}}≔\frac{\partial }{\partial t}{\text{d}}_{\text{f}}{|}_{\text{X}},& (34)\end{array}$

where $\frac{\partial }{\partial t}(·){}_{|\mathbf{\text{X}}}$ is the derivative with respect to t with X being fixed, and the moving interface between fluid and solid domain as

$\begin{array}{ll}{\Gamma }_{\text{f},\text{mov}}^t≔\partial {\Omega }_{\text{f}}^t\backslash {\displaystyle \underset{i=1}{\overset{n_{\text{outlets}}}{\cup }}{\Gamma }_{\text{f},\text{outflow},i}^t},& (35)\end{array}$

where ${\Gamma }_{\text{f},\text{outflow},i}^t$ are the individual aortic outlets. The fluid displacement at this point remains unknown and will be specified in section 2.3.3. Combining these concepts, an ALE description of the Navier–Stokes equations can be derived, see e.g., Bazilevs et al. (2013) and Förster et al. (2006),

$\begin{array}{ll}{\rho }_{\text{f}}\left(\frac{\partial }{\partial t}{u}_{\text{f}}{|}_X+\left(u_{\text{f}}-w_{\text{f}}\right)\cdot {\nabla }_x{u}_{\text{f}}\right)-{\nabla }_x\cdot {\sigma }_{\text{f}}(u_{\text{f}},p_{\text{f}})=0\text{ \hspace{1em}\hspace{1em}\hspace{1em}on }{\Omega }_{\text{f}}^t,& (36)\end{array}$

$\begin{array}{ll}{\nabla }_x·u_{\text{f}}=0\text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}on }{\Omega }_{\text{f}}^t,& (37)\end{array}$

$\begin{array}{ll}{u}_{\text{f}}=g_{\text{mov}}\text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}on }{\Gamma }_{\text{f},\text{mov}}^t,& (38)\end{array}$

$\begin{array}{ll}{\sigma }_{\text{f}}(u_{\text{f}},p_{\text{f}})n_{\text{f}}-{\rho }_{\text{f}}\beta ((u_{\text{f}}-w_{\text{f}})·n_{\text{f}})\_u_{\text{f}}=-p_{\text{wk},\text{i}}{n}_{\text{f}}\text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}on each }{\Gamma }_{\text{f},\text{outflow},i}^t,& (39)\end{array}$

$\begin{array}{ll}{u}_{\text{f}}{|}_{t=0}=u_0\text{\hspace{1em}\hspace{1em}\hspace{1em}in }{\Omega }_{\text{f}}^0.& (40)\end{array}$

Along ${\Gamma }_{\text{f},\text{mov}}^t$ we imposed equality between fluid velocity and the velocity of the moving surfaces. Boundary condition (Equation 39) is the ALE equivalent of the outflow stabilization in Equation (30), see Bazilevs et al. (2013, section 8.4.2.3). Details on how domain movement and velocity were chosen in our application will be discussed later in sections 2.3.3 and 2.5.5.

2.3.2. Variational Formulation of the Navier–Stokes Equations

Following Bazilevs et al. (2007), Bazilevs et al. (2013), and Pauli and Behr (2017), the discrete variational formulation of the ALE Equations (36)–(40) can be stated in the following abstract form: find ${\mathbf{\text{u}}}_{\text{f}}^{\text{h}}\in {\left[{{S}}_{\text{h},\mathbf{\text{g}}}^1({{T}}_{\text{N}})\right]}^3,p_{\text{f}}^{\text{h}}\in {{S}}_{\text{h}}^1({{T}}_{\text{N}})$ such that for all ${\mathbf{\text{v}}}^{\text{h}}\in {\left[{{S}}_{h,\mathbf{\text{0}}}^1({{T}}_{\text{N}})\right]}^3$ and for all $q^{\text{h}}\in {{S}}_{\text{h}}^1({{T}}_{\text{N}})$

$\begin{array}{ll}{A}_{\text{NS}}(v^{\text{h}},q^{\text{h}};u_{\text{f}}^{\text{h}},p_{\text{f}}^{\text{h}})+S_{\text{VMS}}(v^{\text{h}},q^{\text{h}};u_{\text{f}}^{\text{h}},p_{\text{f}}^{\text{h}})=F_{\text{NS}}(v^{\text{h}}),& (41)\end{array}$

with the classical bilinear form of the Navier–Stokes equations

$\begin{array}{ll}\begin{array}{ll}{A}_{\text{NS}}(v^{\text{h}},q^{\text{h}};u_{\text{f}}^{\text{h}},p_{\text{f}}^{\text{h}})≔\hfill & {\rho }_{\text{f}}{\displaystyle \underset{{\Omega }_{\text{f}}^t}{\int }{v}^{\text{h}}}·\left(\frac{\partial }{\partial t}{u}_{\text{f}}^{\text{h}}+\left(u_{\text{f}}^{\text{h}}-w_{\text{f}}^{\text{h}}\right)\cdot {\nabla }_x{u}_{\text{f}}^{\text{h}}\right)\text{d}x\hfill \\ \hfill & +{\displaystyle \underset{{\Omega }_{\text{f}}^t}{\int }\epsilon }(v^{\text{h}}):{\sigma }_{\text{f}}(u_{\text{f}}^{\text{h}},p_{\text{f}}^{\text{h}})\text{d}x\hfill \\ \hfill & +{\displaystyle \underset{{\Omega }_{\text{f}}^t}{\int }{q}^{\text{h}}}{\nabla }_x·u_{\text{f}}^{\text{h}}\text{ d}x-{\rho }_{\text{f}}\beta {\displaystyle \sum _{i=1}^{n_{\text{outlets}}}}\hfill \\ \hfill & {\displaystyle \underset{{\Gamma }_{\text{f},\text{outflow},i}^t}{\int }\text{}}\text{}((u_{\text{f}}^{\text{h}}-w_{\text{f}}^{\text{h}})\cdot n_{\text{f}})\_v^{\text{h}}\cdot u_{\text{f}}^{\text{h}}{\text{d}}_{s_x},\hfill \end{array}& (42)\end{array}$

the bilinear form S_VMS, which is explained later in Equation (45), and the right-hand side contribution

$\begin{array}{ll}{F}_{\text{NS}}(v^{\text{h}})≔-{\displaystyle \sum _{i=1}^{n_{\text{outlets}}}{p}_{\text{wk},\text{i}}}{\displaystyle \underset{{\Gamma }_{\text{f},\text{outflow},i}^t}{\int }\text{}}\text{}{v}^{\text{h}}·n_{\text{f}}{\text{d}}_{s_x}.& (43)\end{array}$

In Equation (42), ε is the strain-rate tensor and ${\mathbf{\text{w}}}_{\text{f}}^{\text{h}}$ is the discrete counterpart of the fluid domain velocity w_f, i.e.,

$\begin{array}{ll}{w}_{\text{f}}^{\text{h}}(t^{n+1},X)=\frac{d_{\text{f}}(t^{n+1},X)-d_{\text{f}}(t^n,X)}{\Delta t}.& (44)\end{array}$

The FE function space ${{S}}_{\text{h},*}^1({{T}}_{\text{N}})$ is the conformal trial space of piecewise linear, globally continuous basis functions over a decomposition ${{T}}_{\text{N}}$ of ${\Omega }_{\text{f}}^t$ into N simplicial elements constrained by v^h = * on essential boundaries. The FE function space ${{S}}_{\text{h}}^1({{T}}_{\text{N}})$ denotes the same space without constraints. For further details we refer to Brenner and Scott (2007) and Steinbach (2007).

From a mathematical point of view, the Navier–Stokes equation can be seen as a multidimensional convection–diffusion equation with pressure acting as a Lagrangian multiplier of the incompressibility constraint. In the common case where velocity and pressure are retained as unknowns, as above, the Ladyzhenskaya–Babuška–Brezzi (LBB) condition has to be satisfied by the velocity and pressure spaces (Donea and Huerta, 2003). A violation of the LBB condition may lead to pressure oscillations. Stabilization techniques allowing the circumvention of the LBB condition exist and have been extensively studied (see for example Hughes et al., 1986; Franca and Hughes, 1988; Douglas and Wang, 1989; Bochev et al., 2006). However, with increasing Reynolds number the Navier–Stokes equations become convection dominated. This requires increasingly finer mesh resolutions to accurately resolve finer flow details which, eventually, renders numerical solution in this form computationally intractable. As a remedy, one can resort to using turbulence models. In particular, in this study the residual based variational multiscale turbulence model (RBVMS), see Hughes (1995), Bazilevs et al. (2007), Bazilevs et al. (2013), and Pauli and Behr (2017) was employed which acts as a stabilization and a turbulence model. The underlying main idea is to split the unknown solution into resolvable (coarse) and unresolvable (fine) scales by the FE approximation, where the finer scale details are taken into account based on element residuals. For details on the derivation we refer to elsewhere (Bazilevs et al., 2007). The term S_VMS in Equation (41) denotes the bilinear form of the RBVMS formulation and reads as

$\begin{array}{ll}\begin{array}{l}{S}_{\text{VMS}}(v^{\text{h}},q^{\text{h}};u_{\text{f}}^{\text{h}},p_{\text{f}}^{\text{h}})≔\hfill \\ \frac{1}{{\rho }_{\text{f}}}{\displaystyle \sum _{l=1}^{n_{\text{el}}}{\displaystyle \underset{{\tau }_{\ell }}{\int }{\tau }_{\text{MOM}}}}\left({\rho }_{\text{f}}\left(u_{\text{f}}^{\text{h}}-w_{\text{f}}^{\text{h}}\right)\cdot {\nabla }_x{v}^{\text{h}}+q^{\text{h}}\right)\cdot r_{\text{MOM}}(u_{\text{f}}^{\text{h}},p_{\text{f}}^{\text{h}})\text{d}x\hfill \\ +{\displaystyle \sum _{l=1}^{n_{\text{el}}}{\displaystyle \underset{{\tau }_{\ell }}{\int }{\tau }_{\text{CONT}}}}{\nabla }_x·v^{\text{h}}{\nabla }_x·u_{\text{f}}^{\text{h}}\text{d}x\hfill \\ -{\displaystyle \sum _{l=1}^{n_{\text{el}}}{\displaystyle \underset{{\tau }_{\ell }}{\int }{\tau }_{\text{MOM}}}}{v}^{\text{h}}\cdot \left({\nabla }_x{u}_{\text{f}}^{\text{h}}{r}_{\text{MOM}}(u_{\text{f}}^{\text{h}},p_{\text{f}}^{\text{h}})\right)\text{d}x\hfill \\ -\frac{1}{{\rho }_{\text{f}}}{\displaystyle \sum _{l=1}^{n_{\text{el}}}{\displaystyle \underset{{\tau }_{\ell }}{\int }{\tau }_{\text{MOM}}^2}}\epsilon (v^{\text{h}}):(r_{\text{MOM}}(u_{\text{f}}^{\text{h}},p_{\text{f}}^{\text{h}})\otimes r_{\text{MOM}}(u_{\text{f}}^{\text{h}},p_{\text{f}}^{\text{h}}))\text{dx},\hfill \end{array}& (45)\end{array}$

where the vector r_MOM is defined as

$\begin{array}{ll}{r}_{\text{MOM}}(u_{\text{f}}^{\text{h}},p_{\text{f}}^{\text{h}})≔{\rho }_{\text{f}}\left(\frac{\partial }{\partial t}{u}_{\text{f}}^{\text{h}}+\left(u_{\text{f}}^{\text{h}}-w_{\text{f}}^{\text{h}}\right)\cdot {\nabla }_x{u}_{\text{f}}^{\text{h}}\right)-{\nabla }_x·{\sigma }_{\text{f}}(u_{\text{f}}^{\text{h}},p_{\text{f}}^{\text{h}}).& (46)\end{array}$

The definition of the parameters τ_MOM, τ_CONT according to Pauli and Behr (2017) is given by

$\begin{array}{ll}{\tau }_{\text{MOM}}≔\min \left\{{\left(\frac{4}{\Delta t^2}+(u_{\text{f}}^{\text{h}}-w_{\text{f}}^{\text{h}})\cdot \mathbf{\text{G}}(u_{\text{f}}^{\text{h}}-w_{\text{f}}^{\text{h}})\right)}^{-\frac{1}{2}},\frac{{\rho }_{\text{f}}{C}_{\text{M}}}{{\mu }_{\text{f}}\sqrt{\mathbf{\text{G}}:\mathbf{\text{G}}}}\right\},& (47)\end{array}$

with Δt being the time step size and $\mathbf{\text{G}}:={\frac{\partial \xi }{\partial \mathbf{\text{x}}}}^{\top }\mathbf{\text{K}}\frac{\partial \xi }{\partial \mathbf{\text{x}}}$, where $\frac{\partial \xi }{\partial \mathbf{\text{x}}}$ denotes the Jacobian of the mapping from a physical FE to the reference FE, the tensor K is defined as

$\begin{array}{ll}\mathbf{\text{K}}≔\frac{1}{2\sqrt[3]{2}}\left(\begin{array}{ccc}3& -1& -1\\ -1& 3& -1\\ -1& -1& 3\end{array}\right)& (48)\end{array}$

and the constant C_M = 0.0285. Further, the stabilization parameter τ_CONT is defined as

$\begin{array}{ll}{\tau }_{\text{CONT}}≔\frac{1}{{\tau }_{\text{MOM}}{g}_{\text{f}}·g_{\text{f}}},& (49)\end{array}$

$\begin{array}{ll}{g}_{\text{f},i}≔{\displaystyle \sum _{j=1}^3{\left(\frac{\partial \xi }{\partial x}\right)}_{ji}}.& (50)\end{array}$

2.3.3. Em-based Kinematic Driver Model

Displacements computed with the EM model were used to prescribe the kinematics of the blood pool mesh which in turn was used for simulating hemodynamics in the CFD model. This was achieved by imposing ${\mathbf{\text{g}}}_{\text{mov}}=\frac{\partial }{\partial t}{\mathbf{\text{d}}}_{\text{s}}$ in Equation (38). Since the surface of the reference CFD blood pool mesh, $\partial {\Omega }_{\text{f}}^0$, is not conformal with the surface of the reference EM blood pool mesh, ${\Omega }_{\text{s},\text{bp}}^0$, and the overlap of the two surfaces is imperfect due to smoothing of $\partial {\Omega }_{\text{f}}^0$ and remeshing of ${\Omega }_{\text{f}}^0$, a direct transfer of displacements between the two surfaces is not readily feasible. As a remedy, we proceeded as follows. After solving the EM problem the subset of displacements ${\widetilde{\mathbf{\text{d}}}}_{\text{s}}$ that form the endocardial interface with the blood pool, ${\Gamma }_{\text{s},\text{bp}}^0$, were extracted from the solution d_s defined at ${\Omega }_{\text{s}}^0$. Since the mesh interface between ${\Omega }_{\text{s}}^0$ and ${\Omega }_{\text{s},\text{bp}}^0$ is conformal the extracted displacements can be applied as inhomogeneous time-varying Dirichlet boundary conditions to the blood pool mesh ${\Omega }_{\text{s},\text{bp}}^0$ to solve a linear elastic problem given as

$\begin{array}{ll}-{\nabla }_X·\sigma (d_{\text{s}}(t))=0\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}in }{\Omega }_{s,\text{bp}}^0,& (51)\end{array}$

$\begin{array}{ll}{d}_{\text{s}}(t)={\tilde{d}}_{\text{s}}(t)\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}on }\partial {\Omega }_{s,\text{bp}}^0,& (52)\end{array}$

where stress and strain tensor are

$\begin{array}{ll}\sigma (d_{\text{s}})≔\frac{E}{1+\nu }\left(\frac{\nu }{1-2\nu }{\nabla }_X·d_{\text{s}}\mathbf{\text{I}}+\epsilon (d_{\text{s}})\right),& (53)\end{array}$

$\begin{array}{ll}\epsilon (d_{\text{s}})≔\frac{1}{2}\left({\nabla }_X{d}_{\text{s}}+{({\nabla }_X{d}_{\text{s}})}^{\top }\right),& (54)\end{array}$

the constant E is Young's modulus in kPa and the constant ν is Poisson's ratio which is dimensionless in the range of [−1, 0.5). Combining the solutions d_s computed for ${\Omega }_{\text{s}}^0$ and ${\Omega }_{\text{s},\text{bp}}^0$ yields displacements d_s for ${\Omega }_{\text{s},\text{total}}^0$. Since $\partial {\Omega }_{\text{f}}^0$ is fully embedded in this domain, ${\Omega }_{\text{s},\text{total}}^0$ ${\Omega }_{\text{s},\text{total}}^0$ can be used as a hanging background mesh for interpolating displacements onto the blood pool mesh, ${\Omega }_{\text{f}}^0$, used for CFD simulations. However, for reasons of mesh quality, interpolation is solely applied on the boundary ${\Omega }_{\text{f}}^0$ itself, and to find the interior displacement field the exact same linear elastic problem 51–54, is solved for d_f instead of d_s.

In both patient cases studied, ejection fractions were large leading to a substantial deformation of the blood pool mesh ${\Omega }_{\text{f}}^t$. To maintain mesh quality under such large deformations the parameters E and ν governing stiffness and incompressibility of the material were altered accordingly. Initially, a fixed E₀ and ν₀ was chosen while the subsequent modification of E and ν was guided by a combination of the two following strategies.

 i) Quality based stiffening: For each element τ_ℓ in the fluid mesh a tetrahedral quality indicator κ(τ_ℓ) based on the movement from the previous time step was calculated, see Freitag and Knupp (2002) and Kanchi and Masud (2007), and rescaled such that for elements of good quality κ is close to 1, while for elements with poor quality κ tends toward infinity. Eventually, the parameter E was multiplied by κ within each element.

ii) ν-Volume based stiffening: For larger deformation elements in the fluid mesh may collapse or even invert, yielding a zero or negative volume. When solving Equations (51)–(54), the current element volumes were tracked and a volume ratio relative to an undeformed reference element was computed as $\frac{\left|{\tau }_{\ell }\right|}{\left|{\widehat{\tau }}_{\ell }\right|}$. For ratios below a predefined critical value the parameter ν was set close to 0.5 to make this element nearly incompressible.

2.4. Numerical Solution

Spatio-temporal discretization of all PDEs and the solution of the arising systems of equations relied upon the Cardiac Arrhythmia Research Package (CARP), see Vigmond et al. (2003). Numerical details on FE discretization (Rocha et al., 2011) and solution of EP (Vigmond et al., 2008b; Neic et al., 2012, 2017) and EM (Augustin et al., 2016b) have been discussed in detail elsewhere. FE discretization and solution of the Navier–Stokes equations were implemented recently using the same numerical framework which was extended to account for non-linear saddle-point problems arising from the discretized CFD equations.

Two time discretization schemes were implemented and compared for the applications in mind, and a computationally cheap semi-implicit scheme, modified from Forti (2016, section 1.4.2), showed similar results to the more expensive fully-implicit generalized-α method (Jansen et al., 2000). Hence, all results in section 3 were obtained using the semi-implicit scheme; to advance from time step tⁿ to tⁿ⁺¹, only a linear block system needs to be solved, where each block depends on data from the previous time step only. Solvers for the block system were taken from the PETSc library (Balay et al., 1997, 2016a,b). We used a right preconditoned flexible GMRES method with PETSc fieldsplit preconditioning (Silvester et al., 2001; Elman et al., 2008) which in turn uses BoomerAMG (Van Emden and Yang, 2002) to approximate sub-block inverses. While the time step size for mechanics and CFD was the same, Δt_mech = Δt_CFD = 0.5 ms, it was significantly smaller for EP, where Δt_EP = 25 µs.

The implementation of the CFD solvers has been subjected to various validation procedures against standard CFD benchmarks (Schäfer et al., 1996). All simulations were executed at the national HPC computing facility ARCHER in the United Kingdom using 384 and 768 cores for EM and CFD simulations, respectively.

2.5. Model Parameterization

2.5.1. Electrophysiology

Electrical activation sequences were indirectly parameterized using the QRS complex of a given patient's ECG as guidance. Unlike in previous studies (Augustin et al., 2016a), we refrained from a detailed parameterization which aimed at reproducing the QRS complex of the ECG for a given patient by finding appropriate locations and timings for the main fascicles of the cardiac conduction system in the LV. Rather, default locations and timings were used which yielded a total activation time within the physiological range.

2.5.2. Passive Biomechanics

The LV myocardium was characterized as a hyperelastic, nearly incompressible, transversely isotropic material with a nonlinear stress–strain relationship (Guccione et al., 1995). Orthotropic material axes were aligned with the local fiber, sheet and sheet normal directions. To remove rigid body motion, homogeneous displacement boundary conditions were applied by fixing the terminal rims of the clipped brachiocephalic, left common carotid and left subclavian arteries as well as the clipped rim of the aorta descendens, see Figure 1. The model was stabilized by resting the LV apex on an elastic cushion of which the bottom face was rigidly anchored also by applying homogeneous displacement boundary conditions.

The constitutive model was fitted to recorded clinical data as previously reported with minor modifications (Augustin et al., 2016a). The passive biomechanical model governed by the strain-energy function given in Equation (17) was fitted to approximate the end-diastolic pressure-volume relation (EDPVR). Due to limitations in the recorded data we refrained from directly fitting the model to the recorded pressure and volume data. Rather, only one data pair—EDV and end-diastolic pressure (EDP)—was used to fit the stress-free residual volume to the empiric Klotz relation (Klotz et al., 2007) by adjusting the isotropic scaling parameter C_Guc in Equation (17). As the model anatomy was built from a segmented 3DWH MRI scan—acquired during diastasis—the FE model was inflated to increase the volume of the cavity by the difference between the volume at mid diastasis and the EDV. Using the end-diastolic geometry, default material parameters and the recorded EDP, an initial guess of the stress-free reference configuration was computed by unloading the model using a backward displacement method (Sellier, 2011; Bols et al., 2013; Krishnamurthy et al., 2013). The unloading procedure was repeated with varying trial material parameters, C_Guc, until the difference between the unstressed LV volume of the model and the prediction of the Klotz relation was less than 5 %.

2.5.3. Active Stresses

Parameters of the active stress model were fitted during IVC and ejection phase. During IVC the LV volume was held constant (Gurev et al., 2015) and the parameters of the active stress given in Equation (20) rate of contraction, τ_c, and peak active stress, S_peak, were manually adjusted to fit the maximum rate of rise of pressure, (dP/dt)_max, and peak pressure, p_lv.

2.5.4. Afterload

When the LV pressure p_lv exceeded the aortic pressure, p_ao, ejection was initiated by connecting the LV model with the lumped 3-element Windkessel model (Westerhof et al., 1971). Volume traces recorded from a given patient during ejection were used as input to compute aortic pressure traces by solving Equation (23). Both types of data were not recorded simultaneously as volume traces were computed from Cine MRI scans and pressure traces were recorded later invasively by catheterization. Volume and pressure traces were synchronized in time by aligning the onset of ejection of the volume trace V_lv(t) with the instant of opening of the aortic valve in the pressure trace p_ao(t). In those cases where heart rates were markedly different between the two measurements, volume traces were scaled in time to adjust LV ejection time (LVET) to the duration of ejection in the pressure traces, that is, the time elapsed between opening and closing of the aortic valve as these two instants in time were clearly identifiable in all traces p_ao(t), see Figure 3. Moreover, volume traces were offset to ensure that the model volume based on the segmentation of the 3DWH scan acquired during diastasis matched up with the Cine-MRI based volume trace at mid diastasis. The parameter space of the Windkessel model comprising characteristic impedance of the aorta, Z_c, as well as resistance, R, and compliance, C, of the arterial system was sampled using a recently developed stochastic sampling approach (Crozier et al., 2016b).

FIGURE 3

Figure 3. (A) Invasive clinical recordings from cases 28-Pre and 44-Pre. Top: Recorded aortic pressure P_ao (black curve) and recorded LV pressure P_LV (blue curve). Marked with dashed lines are Systolic pressure P_sys, mean arterial pressure MAP, and diastolic pressure P_dia; Center: Volume change in the LV, V_LV, in red ranging from end-diastolic volume EDV to end-systolic volume ESV. Bottom: LV flow Q_LV in orange with marked peak flow Q_peak. (B) Comparison of EM simulations and clinical data. Upper part shows a comparison of the LV model in end-diastolic (colored opaquley blue) and end-systolic configuration (colored by displacement). Lower part shows comparison of clinical (colored blue) and simulated PV loops (colored red). The dashed orange curve shows the ideal Klotz curve, while the green curve shows the simulated Klotz curve, with volume of stress-free unloaded configuration marked as V₀.

Numerous box constraints were used to constrain the search space of parameter sweeps. In particular, we used reported measurements in humans to define the mean values and restricted the search space for each parameter to fall within ± 20% around the mean. Due to high frequency errors introduced by the pressure transducer we refrained from computing norms ||p_{ao, meas} − p_{ao, fit}|| to quantify the deviations of fitted from measured pressure and opted for manual selection using three criteria, aortic peak pressure, p_ao, closing pressure of aortic valve and exponential decay of p_ao during diastole. For the sake of fitting Z_c we assumed p_ao ≈ p_lv since transvalvular pressure gradients in all patients were very minor.

2.5.5. CFD Boundary Conditions

The validated EM models yield the time-dependent displacement fields, d_s, which were transferred onto the fluid domain to drive simulations of blood flow in LV and aorta as described in section 2.3.3 yielding d_f(t, x) defined on the whole CFD mesh. Figure 4G shows a summary of the boundary conditions. On the boundary ${\Gamma }_{\text{f},\text{mov}}^t$ a Dirichlet boundary condition enforcing the mesh velocity ${\mathbf{\text{w}}}_{\text{f}}^{\text{h}}$ is applied. On each aortic outlet Γ_{f, outflow, i}(t) a 3-Element Windkessel model as described in section 2.2.3 is attached. Further, the stabilization parameter β in Equation (39) was set to 0.2. Estimation of the input parameters for the hemodynamical Windkessel equations relied on an extension of the simple hydraulic analog of Ohm's law. Given the patient specific MAP, CO, and a percentage α_i of total CO running through the outlet the resistance R_i was estimated as

$\begin{array}{ll}{R}_i\approx \frac{\text{MAP}}{{\alpha }_i\text{CO}}.& (55)\end{array}$

The percentages α_i were obtained either by measurement or by applying Murray's law (Murray, 1926). The impedances Z_i were chosen as 5 % of R_i, and the compliances C_i were chosen such that R_iC_i ≈ 1, 000ms. To keep the semi-implicit character of the CFD system the Windkessel equations were solved with a semi-implicit backward Euler method using the flow $q_i^n$ through the aortic outlet, from the previous time step as input.

FIGURE 4

Figure 4. Processing workflow used for generating blood pool FE models: (A,B) Elements labeled as blood pool or valve were extracted from the mesh used for EM modeling. (C) Surfaces of extracted meshes were smoothed to avoid numerical instabilities due to reentrant corners resulting from a jagged surface. A closeup view of the smoothing effect is displayed in the upper right. The smoothed surface is then used as input for the fluid mesh generation. (D,E) Comparison of the smoothed and unsmoothed blood pool mesh immersed in the original EM mesh. (F) Closeup view of the generated boundary layer mesh. (G) Boundary conditions used for CFD. Moving wall boundary ${\Gamma }_{\text{f},\text{mov}}^t$ colored in orange, outlet boundaries ${\Gamma }_{\text{f},\text{outflow},i}^t$ colored in blue with attached illustration of the 3-element Windkessel models.

3. Results

3.1. Building Electromechanical Kinematic Driver Models

Using a previously developed automated workflow (Crozier et al., 2016a), anatomical FE models of LV and aorta were built for patient cases 28-Pre and 44-Pre based on segmented imaging data acquired under pre-treatment conditions. Figure 1 illustrates the key processing steps and the resulting FE model for case 28-Pre. For the case 28-Pre the CoA was repaired by a virtual dilatation procedure applied to the segmented image data with the aim to restore normal cross sectional areas. Subsequently, a new FE mesh was generated referred to as 28-Post, which was essentially identical to 28-Pre, with the only difference being the anatomical adjustment of the CoA in the aortic arch to the target post-treatment anatomy after stenting, see Figure 5.

FIGURE 5

Figure 5. CoA anatomy of case 28 before and after virtual stenting procedure. CoA location is indicated with a red circle.

Passive biomechanical properties, afterload and active stress models of cases 28-Pre and 44-Pre were parameterized using clinically recorded pressure and volume data under pre-treatment conditions, see Figure 3A. The fitted final parameters used are summarized in Table 2. The goodness of fit of both integrated EM models was verified by standard PV loop analysis as shown in Figure 3B. Results of a quantitative comparison with clinically derived metrics including EF, EDV and ESV, CO, and peak systolic pressure are summarized in Table 3.

TABLE 3

Table 3. Comparison of clinical indicators and indicators computed from simulation for the EM models.

3.2. Blood Pool FE Modeling for CFD

Conformal FE blood pool meshes were extracted from EM FE meshes, surfaces were smoothed and used for volumetric remeshing with increased spatial resolution including boundary layers. The corresponding workflow is illustrated in Figure 4.

Kinematics of the EM model were transferred to the CFD blood pool mesh and the result is illustrated in terms of displacements d_s, d_f in Figure 6II. Due to the large EF of about 65 % for both 28-Pre and 44-Pre, the blood pool underwent a significant deformation. However, using a combination of element quality and ν-Volume based stiffening with an initial Young's Modulus E₀ = 100 kPa and Poisson's ratio ν₀ = 0.3, sufficient element quality was preserved throughout the entire ejection phase and numerical instabilities could be avoided. Figure 6I shows the 80th-percentile of bad element quality against the number of linear iterations required for convergence for the 28-Pre case. The quality of elements was calculated with the same quality inidcator (Freitag and Knupp, 2002; Kanchi and Masud, 2007) as described in section 2.3.3 but was rescaled to the interval [0, 1], with the best element quality being 0 and the worst element quality being 1. The modest increase in iteration numbers of the iterative preconditioned GMRES solver provides indirect evidence of sufficiently preserved mesh quality (see Figure 6). Spatially, most lower quality elements were located in the CFD boundary layer.

FIGURE 6

Figure 6. (I) shows quality analysis for case 28-Pre. Spatial locations of elements of poor quality > 0.8 (in red) are shown at the top for different snapshots of deformation (green lines in graph). The graph below shows linear iterations per time step (in blue) and percentage of elements with poor quality > 0.8 (in red). (II) shows the processing stages of kinematic transfer for the 28-Pre case at maximum displacement. (A) Displacement d_s on EM mesh ${\Omega }_{\text{s}}^0$. (B) Displacement d_s extended to conformal EM blood pool mesh ${\Omega }_{\text{s},\text{total}}^0$ which serves as hanging background mesh for the kinematic transfer onto the CFD blood pool mesh ${\Omega }_{\text{f}}^0$. (C) Displacement d_s on ${\Omega }_{\text{s}}^0$ superimposed with fluid mesh displacement d_f on ${\Omega }_{\text{f}}^0$.

3.3. Numerical CFD Benchmarks

The implementation of the Navier–Stokes solver was verified by solving a set of standardized benchmark problems, see Schäfer et al. (1996). Computational performance was evaluated by performing strong scaling experiments by repeating the post-treatment hemodynamics simulation of case 28-Post with varying numbers of cores ranging from 96 to 1.536. Details on computational complexity and costs are summarized in Table 4. For temporal discretization a time step of Δt = 0.5 ms was used to simulate the ejection phase lasting for 208 ms. The overall discrete system comprised 5,177,056 degrees of freedom, which was solved over 416 time steps. Strong scaling results are summarized in Figure 7. Efficient strong scaling behavior was observed up to 768 cores with parallel efficiency slowly degrading from 100 % at 96 cores down to 55 % at 768 cores. Scalability stalled when doubling the core count to 1,536 which reduced the degrees of freedom per parallel partition down to 3,386. Parallel efficiency dropped to 27 % which is attributed due to the unfavorable ratio between local compute work and communication.

TABLE 4

Table 4. Discretization details for the studied cases.

FIGURE 7

Figure 7. Results of strong scaling benchmark based on case 28-Post with 5.2 million overall degrees of freedom. T_AvgSolv is the total solving time divided by the total amount of linear iterations per simulation run.

3.4. Simulating Cardiac and Cardiovascular Hemodynamics

Hemodynamics in the LV and aorta was simulated using the EM simulations as a kinematic driver. Flow rates through various aortic cross sections and outflow orifices were calculated as the integral over measured fluxes through the cross-sectional plane for both 4D VEC MRI and simulated flow data. At locations of interest which were ε_DSC, ε_BCA, ε_LCA, and ε_LSCA denoting cross sections in the aorta descendens and the orifices of brachocephalic, left carotid and left subclavian artery, respectively, relative flows were computed from 4D VEC MRI data as fractions α_i expressed in percent of the total peak flow through the aorta ascendens as determined over the plane ε_ASC. For those planes of interest where measurements were not feasible due to noise, flow percentages were estimated based on Murray's law. Flow curves during ejection at selected cross sections are shown in Figures 8A,E. MAP and computed mean flow through each outlet orifice were used to determine the parameters of the coupled Windkessel models of afterload in Equations (24, 25), see Table 5. In the 28-Pre case this resulted in flow splits of α_i ≈ 23, 51.3, 12.83, and 12.83% whereas in the 44-Pre case the flow split ratios were α_i ≈ 5.68, 57.45, and 34.01% for ε_DSC, ε_BCA, ε_LCA, and ε_LSCA, respectively.

TABLE 5

Table 5. Windkessel parameters for each outlet of cases 28-Pre and 44-Pre.

FIGURE 8

Figure 8. CFD results. (A,E) show the given clinical measurements for flow through different planes. The planes are depicted in (B,F). (B,F) also depict the pressure along the centerlines at peak flow conditions at t = 167ms and t = 142ms respectively. (C) shows velocity streamlines at peak flow. (D) shows the relative pressure map from the Pressure–Poisson mapping used for validating the pressure drop in our simulations. (G,H) show velocity streamlines at peak flow and t = 200ms for case 44-Pre.

For the CFD analysis a time step of Δt = 0.5 ms was used. The ejection phases of the EM simulations were chosen as time horizons for the CFD simulation which lasted from t = 90 ms to t = 302 ms in the 28-Pre case and from t = 70 ms to t = 329 ms in the 44-Pre case, yielding 424 and 518 time steps, respectively. The Windkessel parameters for each outlet, calculated as described in section 2.5.5, are summarized in Table 5. Pressure p_f along the centerline s_c and fluxes through the planes ε_DSC, ε_LSC, ε_BCA, and ε_ASC were computed at the instant of peak flow in the aorta ascendens and compared against measured data, which were pressures derived from Pressure–Poisson mapping (see Figure 8D) and 4D VEC MRI fluxes. For case 28-Pre pressure drops were calculated from the pressure values on the intersection of the centerline and ε_DSC, ε_ASC respectively. Further, we calculated the average pressure over the aforementioned planes as well. Both ways yielded a simulated pressure drop across the CoA of ≈ 29.2 mmHg which agreed well with the clinically estimated pressure drop of ≈ 30 mmHg. Furthermore, we calculated the flux through the various planes and compared them against the clinically estimated fluxes. A quantitative comparison of fluxes is given in Table 6. Figures 8C,G,H show velocity profiles at peak flow condtions. Figures 8B,F show the pressure along the centerlines, the velocity field ${\overrightarrow{v}}_{\text{f}}$ through the plane ε_ASC, and the position of all planes used for evaluating fluxes. Supplementary Materials 1, 2 contain videos of the time evolution of the velocity distribution for cases 28-Pre and 44-Pre.

TABLE 6

Table 6. Comparison of clincal estimated flow rates and simulated flow rates through the various planes for cases 28-Pre and 44-Pre.

3.5. Post-treatment Simulations

Simulations of case 28-Pre were repeated on geometry of case 28-Post using almost the same set of parameters, see Table 2. Only S_peak was slightly adjusted, which resulted in a better peak pressure value in the LV. The geometry of case 28-Post was almost identical to case 28-Pre with the only exception being the virtual repair of CoA anatomy. In this scenario only pre- and post-treatment simulations were compared to evaluate their relative differences in terms of pressure and flow velocities. Figure 9 shows results. Pressure drops were calculated as in section 3.4 for both scenarios. For 28-Pre we calculated a pressure drop of ≈ 29.2 mmHg while for 28-Post a pressure drop of ≈ 14.15 mmHg was calculated.

FIGURE 9

Figure 9. Comparison of cases 28-Pre and 28-Post. Shown on the left are the pressures along the centerline at peak flow. Depicted in the middle are the slices used for calculating the pressure drops. Shown on the right are velocity streamlines at peak flow.

4. Discussion

In this study, we report on the progress made toward a novel EMF model of the human LV that is entirely based on first principles and as such, in principle, is able to represent all cause-effect relationships with full biophysical detail. Unlike in the majority of cardiac CFD studies where the use of image-based kinematic driver models prevails, EM LV and aorta models of CoA patients were employed to serve as a kinematic driver to a computational model of hemodynamics in the LV cavity and aorta. A hybrid two stage modeling approach was adopted with regard to hemodynamics where EM and CFD model are executed sequentially. First, in the EM simulations the afterload imposed by the circulatory system upon the LV was represented by a lumped model to compute LV kinematics. These EM models were carefully fitted to available clinical data to replicate important clinical metrics characterizing hemodynamic and biomechanical work performed by the LV (Gsell et al., under review). In a subsequent step, a full-blown ALE-based CFD model with moving domain boundaries was unidirectionally or weakly coupled to the EM model. The motion of the fluid domain was driven by the kinematics of the EM model. Kinematics was transferred from EM mesh onto the CFD blood pool mesh by generating a combined kinematic model comprising LV, valve, aortic structure and a conformal blood pool mesh which served as a hanging background mesh for interpolation. The higher resolution blood pool CFD mesh with refined boundary layers was fully immersed in the EM background mesh. Kinematics was transferred by interpolation only onto the surface of the CFD blood pool mesh and extended into the volume of the blood pool by solving a linear solid mechanics problem.

We show validation results for two selected clinical CoA cases under pre-treatment conditions and compare between pre-treatment and post-treatment for one patient case in which the CoA was anatomically modified by a virtual stenting procedure. Further, we demonstrate numerical tractability of the implemented approach by providing strong scaling benchmark results. The overall cost of the entire work flow for building, fitting and execution of EMF simulations is comparable to plain image-based kinematic driver models (Mittal et al., 2016), suggesting that the proposed methodology may be, in principle, compatible with clinical time scales.

4.1. Biomechanical Modeling vs. Image-Based Kinematics

Modalities such as CMR and Cardiac CT on the other hand, provide excellent spatial resolution. CMR has an in-plane resolution of 1.5 × 1.5 mm, but more limited through-plane resolution (typically about 8 mm) while CT is capable of isotropic spatial resolution on the sub millimeter scale (≈ 0.5 mm) and clear delineation of trabeculae and lumen boundaries. CMR has the advantage of higher temporal resolution (30–50 ms) while temporal resolution in CT depends on the scanning system (50–200 ms). This is orders-of-magnitude lower than the temporal resolution required for the flow simulation (≈ 1, 000 phases per cardiac cycle) and appropriate interpolation methods need to be employed to create CFD-ready models. This stage of model generation has been very difficult to automate, and remains the biggest bottleneck for patient-specific cardiac flow modeling. Compared to pure image-based kinematic approaches our model is able to compute, e.g., the spatio-temporal distribution of wall stresses, power density, the length of diastolic intervals available for myocardial perfusion, O₂ consumption, and metabolic supply/demand ratios. The variations of all these parameters in response to a changed afterload and many other biomarkers of physiological interest can be derived, which is not feasible with image-based models.

4.2. Kinematic Transfer to CFD Blood Pool Model

Both patients modeled in this study featured healthy EFs of > 60 %, that is, EF was ≈ 65 % in both cases. At a such high EFs the wall motion of the LV is significant, leading to substantial reductions in the LV blood pool volume. IB methods (Vigmond et al., 2008a; Seo and Mittal, 2013; Choi et al., 2015) are known to be more convenient to cope with the large deformation of the CFD blood pool (Quarteroni et al., 2017). IB methods and other non-boundary-fitting methods rely on a fixed fluid mesh and the moving wall of the ventricle is not explicitly tracked. The coupling between the CFD mesh and the structure is performed via Dirac Delta functions (IB) or Lagrange multipliers (fictitious domain methods) and is usually realized by introducing additional degrees of freedom on interface cut elements. While mesh generation is only necessary prior to computation fixed mesh methods typically require adaptive mesh refinement or modifications (Wang and Liu, 2004) to obtain reasonable accuracy for the solution near the fluid-solid interface.

In contrast, ALE algorithms capture the fluid-solid interface more accurately, are in general stable and easy to implement, no extra degrees of freedoms are introduced, and computational costs are low in comparison (Tallec and Mouro, 2001; van Loon et al., 2007). However, it is often assumed that unstructured FE approaches, as implemented in this paper, critically depend on automatic remeshing strategies (Long et al., 2013) to keep mesh quality within acceptable bounds (Mittal et al., 2016). Our study demonstrates that this may not necessarily be the case. While the mesh quality decreased with deformation over the course of ejection, the linear elastic deformation of the CFD blood pool mesh combined with the quality-based stiffening approach prevented the degeneration of any elements. The number of elements in which element quality degraded noticeably was very small. As illustrated in Figure 6, virtually all elements of reduced quality were located in the higher resolution boundary layer of the CFD blood pool mesh. According to the element quality metric used, an element quality of 1 refers to a fully degenerated element of zero volume. Despite the significant compression of the blood pool mesh, not a single element was deformed to this degree. Even when applying a stricter threshold where element quality is deemed poor if the quality indicator is > 0.8, which is not critical from a numerical point of view, the number of elements in this range remained small with < 0.8 % (Figure 6). The worst element quality observed in the entire mesh was 0.9994. Using a threshold of > 0.95 where element quality may be sufficiently poor to impact more notably on solver performance, only 24 out of 2,506,987 elements were found. Nonetheless, an increase in number of linear iterations required for convergence was observed which is likely to be linked to the gradual degradation of element quality. The number of iterations per solver step increased from around ≈ 17 iterations during early ejection up to ≈ 80 iterations during late ejection. While the more than fourfold increase in linear iterations negatively impacted overall solver performance and rendered simulations computationally more expensive, the complexity of automatic remeshing was avoided. We consider this a pivotal importance as automatic remeshing in combination with a MPI parallel FE solver is definitely feasible, but highly non-trivial to implement robustly and efficiently.

4.3. Computational Feasibility

Computational feasibility of human scale cardiac simulations by using strongly scalable numerical implementations has been demonstrated previously for electrophysiology (Niederer S. et al., 2011) and mechanics (Augustin et al., 2016b). More recently, we reported on a novel reaction-eikonal model which reduces the cost of EM simulations significantly by alleviating constraints imposed by reaction-diffusion models upon mesh resolution (Neic et al., 2017). In this study, this recent reaction-eikonal approach was used for simulating EM using the same FE grid with an average resolution of ≈ 1 mm for both EP and mechanics. Such lower resolutions suffice for solving for mechanics with sufficient accuracy (Land et al., 2015). The overall reduction in terms of nodes and degrees of freedom reduces the compute cost substantially, rendering simulations in desktop environments feasible. Using 96 cores, EM simulations of a full cardiac cycle only lasted ≈ 180 min which facilitated sufficiently short simulation cycles for efficient model fitting. The entire workflow for building and parameterizing one patient-specific EM model is feasible within a day.

Owing to the higher resolution of the blood pool mesh and the presences of a refined boundary layer the number of nodes and degrees of freedom were higher than for EM simulations, around 350,000/1,500,000 nodes/degrees of freedom for case 28-Pre and 400,000/1,700,000 nodes/degrees of freedom for case 44-Pre, respectively. To assess strong scaling properties of our CFD solver implementation, the resolution was further increased to 1,300,000/5,000,000 nodes/degrees of freedom for case 28-Post to cover a wider range of core counts. Strong scaling efficiency leveled off when doubling from 768 to 1,536 cores. Local compute load with 1,536 was 900/2,600 nodes/degrees of freedom per core. The patient simulations were performed using 384 cores, resulting in a load per core of about 900/2,700 nodes/dofs, respectively. At these resolutions CFD simulations were executed in ≈ 40 min, suggesting that compatibility with clinical time frames will be achievable.

4.4. Limitations

In the presented modeling approach numerous simplifying assumptions were made which may affect the biophysical fidelity of the model. In particular, while the aorta was taken into account as a solid structure in the EM simulations, its biomechanical description was simplified by assuming isotropic behavior, that is, the fibrous organization of aortic walls remained unaccounted for (Augustin et al., 2014). Further, as our main focus was on the EM of the LV and, to a much lesser degree, on the aorta, the aortic lumen remained unpressurized and, in absence of distensibility measurements of the aortic wall, parameters of the passive biomechanics model used for the aortic wall were not fitted. Thus the model of the aorta does not respond to the rise in pressure during ejection with an adequate distension ΔV of its lumen. In the CFD simulations ΔV ≈ 0 translates into a stiff aorta of low compliance which may cause a bias toward overestimation of the computed pressure fields. Further, the influence of the aortic valve upon blood flow was not taken into account. Rather, it was assumed that with the start of ejection the aortic valve is in its full open configuration, which allows blood flow over the entire orifice area and in which the valve does not influence the blood flow out of the LV in a significant way. Since only CoA patients were modeled which showed no indications of AVD this simplifying assumption may be well justified.

A potential main strength of the presented modeling approach—the ability to predict the biomechanical response of the LV to changed flow patterns in the aorta—was not exploited. Due to the weak FSI coupling the immediate feedback of altered flow or changed pressure gradients in the aorta on LV biomechanics was ignored. In our current modeling approach any such feedback must be mediated through changes in the parameterization of the lumped afterload model. However, owing to regulatory mechanism of the circulatory system level this is not directly predictable with the modeling setup used in this study as flow distribution through the four outlets will be influenced by factors which cannot be accounted for in a model comprising only LV, aorta and lumped outflow impedances. In any case, one cannot assume that the computed changes in pressure gradients across a CoA translate directly into a reduction in LV peak pressure. Independently of the modeling approach taken—be it a strongly or weakly coupled FSI model—a lumped model of systemic regulation is likely to be necessary to predict altered LV loading under post-treatment conditions (Arts et al., 2005; Lumens et al., 2009). Compared to a fully coupled FSI model our approach is limited in the sense that CFD simulations do not influence the behavior of the EM model. However, in many clinical settings CFD simulations in the aortic arch and LV with image based kinematics prevail.

Image based kinematic models can only depict the status quo of a patient. With our personalized EM model, based on first principles, we can do simulations altering the motion, simply by changing input paramters. The altered motion is then reflected in the CFD simulation. Examples would include changes in heart beats, infarcts or LBBB conditions.

In this work, the effect of stenting was only accounted for by a geometric change in the computational geometry and an ad hoc adjustement of the lumped model parameters. In future studies, we intend to use a 1-D model of the arterial tree coupled to a 0-D lumped model at the aortic outlets, thus being able to account for the effect of stenting in a more detailed fashion, see for example Quarteroni et al. (2017). As a first step toward our ultimate goal of a fully coupled FSI model, that is based entirely on first principles, we will add the dynamic fluid pressure $\frac{{\rho }_{\text{f}}}{2}{\left|{\mathbf{\text{u}}}_{\text{f}}\right|}^2$ to the pressure of the lumped model (0-D or 1-D). This results in a spatio-temporal pressure inside the LV and the aorta, and to incorporate the dynamic feedback of fluid upon structure we will iterate between a CFD solving step and a EM solving step within each timestep to guarantee a converged solution.

5. Conclusion

Biophysically detailed models of LV EM can be efficiently built and parameterized with clinical data to be considered a viable option for patient-specific simulation. Similar to image-based kinematic models such biophysics-based EM models can be used as a kinematic driver for simulating cardiac and vascular hemodynamics. The cost of model building and execution is comparable between the two approaches. Biophysical EM models offer the significant advantage of being based entirely on first principles and as such, may allow to make predictions of interventions altering pressure and flow patterns onto LV performance. In contrast, image-based kinematics modeling may provide a more accurate representation of blood pool motion, at least under pre-treatment conditions or post-treatment conditions secondary to interventions which do not influence LV kinematics in a significant way.

Author Contributions

EK, GP contributed conception and design of the study; LG, TK acquired and processed clinical data; EK, MG, AN, and CA developed numerical methodology; AP contributed by conceiving modeling workflows and FE meshing; LM and MG developed parameterization of electromechanical model; EK, MG, CA, and GP analyzed and interpreted simulation data; EK, CA, and GP drafted the article; EK, CA, MG, LG, and GP critically revised the article; All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

This research was supported by the grants F3210-N18 and I2760-B30 from the Austrian Science Fund (FWF), the EU grant CardioProof agreement 611232 and a BioTechMed award to GP, and a Marie Skłodowska–Curie fellowship (GA 750835) to CA. We acknowledge PRACE for awarding us access to resource ARCHER based in the UK at EPCC (grant CAMEL) and the Vienna Scientific Cluster VSC-3.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fphys.2018.00538/full#supplementary-material

Footnotes

1. ^http://www.gyrotools.com/products/gt-flow.html

2. ^https://www.mevis.fraunhofer.de/en/solutionpages/mevisflow-non-invasive-interactive-exploration-of-in-vivo-hemodynamics.html

3. ^https://amira.zib.de

4. ^http://www.cgal.org

5. ^http://www.vmtk.org

6. ^http://www.ansys.com/Services/training-center/platform/introduction-to-ansys-icem-cfd-Hexa

References

Andersson, M., Lantz, J., Ebbers, T., and Karlsson, M. (2017). Multidirectional WSS disturbances in stenotic turbulent flows: a pre-and post-intervention study in an aortic coarctation. J. Biomech. 51, 8–16. doi: 10.1016/j.jbiomech.2016.11.064

PubMed Abstract | CrossRef Full Text | Google Scholar

Antiga, L., Piccinelli, M., Botti, L., Ene-Iordache, B., Remuzzi, A., and Steinman, D. A. (2008). An image-based modeling framework for patient-specific computational hemodynamics. Med. Biol. Eng. Comput. 46:1097. doi: 10.1007/s11517-008-0420-1

CrossRef Full Text | Google Scholar

Arts, T., Delhaas, T., Bovendeerd, P., Verbeek, X., and Prinzen, F. (2005). Adaptation to mechanical load determines shape and properties of heart and circulation: the CircAdapt model. Am. J. Physiol. Heart Circ. Physiol. 288, 1943–1954. doi: 10.1152/ajpheart.00444.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Augustin, C. M., Crozier, A., Neic, A., Prassl, A. J., Karabelas, E., Ferreira da Silva, T., et al. (2016a). Patient-specific modeling of left ventricular electromechanics as a driver for haemodynamic analysis. EP Eur. 18(Suppl. 4):iv121–iv129. doi: 10.1093/europace/euw369

PubMed Abstract | CrossRef Full Text | Google Scholar

Augustin, C. M., Holzapfel, G. A., and Steinbach, O. (2014). Classical and all-floating FETI methods for the simulation of arterial tissues. Int. J. Numer. Methods Eng. 99, 290–312. doi: 10.1002/nme.4674

PubMed Abstract | CrossRef Full Text | Google Scholar

Augustin, C. M., Neic, A., Liebmann, M., Prassl, A. J., Niederer, S. A., Haase, G., et al. (2016b). Anatomically accurate high resolution modeling of human whole heart electromechanics: a strongly scalable algebraic multigrid solver method for nonlinear deformation. J. Comput. Phys. 305, 622–646. doi: 10.1016/j.jcp.2015.10.045

PubMed Abstract | CrossRef Full Text | Google Scholar

Balay, S., Abhyankar, S., Adams, M. F., Brown, J., Brune, P., Buschelman, K., et al. (2016a). PETSc Users Manual. Technical Report ANL-95/11 - Revision 3.7, Argonne National Laboratory.

Balay, S., Abhyankar, S., Adams, M. F., Brown, J., Brune, P., Buschelman, K., et al. (2016b). PETSc Web Page. Available online at: http://www.mcs.anl.gov/petsc

Balay, S., Gropp, W. D., McInnes, L. C., and Smith, B. F. (1997). “Efficient management of parallelism in object oriented numerical software libraries,” in Modern Software Tools in Scientific Computing, eds E. Arge, A. M. Bruaset, and H. P. Langtangen (Boston, MA: Birkhäuser Press), 163–202.

Google Scholar

Bayer, J. D., Blake, R. C., Plank, G., and Trayanova, N. A. (2012). A novel rule-based algorithm for assigning myocardial fiber orientation to computational heart models. Ann. Biomed. Eng. 40, 2243–2254. doi: 10.1007/s10439-012-0593-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Bazilevs, Y., Calo, V., Cottrell, J., Hughes, T., Reali, A., and Scovazzi, G. (2007). Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows. Comput. Methods Appl. Mech. Eng. 197, 173–201. doi: 10.1016/j.cma.2007.07.016

CrossRef Full Text | Google Scholar

Bazilevs, Y., Takizawa, K., and Tezduyar, T. E. (2013). Computational Fluid-Structure Interaction: Methods and Applications. Hoboken, NJ: John Wiley & Sons.

Google Scholar

Bertoglio, C., Caiazzo, A., Bazilevs, Y., Braack, M., Esmaily, M., Gravemeier, V., et al. (2017). Benchmark problems for numerical treatment of backflow at open boundaries. Int. J. Numer. Methods Biomed. Eng. 34:e2918. doi: 10.1002/cnm.2918

PubMed Abstract | CrossRef Full Text | Google Scholar

Bochev, P. B., Dohrmann, C. R., and Gunzburger, M. D. (2006). Stabilization of low-order mixed finite elements for the stokes equations. SIAM J. Numer. Anal. 44, 82–101. doi: 10.1137/S0036142905444482

CrossRef Full Text | Google Scholar

Bols, J., Degroote, J., Trachet, B., Verhegghe, B., Segers, P., and Vierendeels, J. (2013). A computational method to assess the in vivo stresses and unloaded configuration of patient-specific blood vessels. J. Comput. Appl. Math. 246, 10–17. doi: 10.1016/j.cam.2012.10.034

CrossRef Full Text | Google Scholar

Braack, M., Mucha, P. B., and Zajaczkowski, W. (2014). Directional do-nothing condition for the Navier–Stokes equations. J. Comput. Math. 32, 507–521. doi: 10.4208/jcm.1405-m4347

CrossRef Full Text | Google Scholar

Brenner, S., and Scott, R. (2007). The Mathematical Theory of Finite Element Methods, Vol. 15. New York, NY: Springer Science & Business Media.

Google Scholar

Chnafa, C., Mendez, S., and Nicoud, F. (2014). Image-based large-eddy simulation in a realistic left heart. Comput. Fluids 94, 173–187. doi: 10.1016/j.compfluid.2014.01.030

CrossRef Full Text | Google Scholar

Choi, Y. J., Constantino, J., Vedula, V., Trayanova, N., and Mittal, R. (2015). A new MRI-based model of heart function with coupled hemodynamics and application to normal and diseased canine left ventricles. Front. Bioeng. Biotechnol. 3:140. doi: 10.3389/fbioe.2015.00140

PubMed Abstract | CrossRef Full Text | Google Scholar

Crozier, A., Augustin, C. M., Neic, A., Prassl, A. J., Holler, M., Fastl, T. E., et al. (2016a). Image-based personalization of cardiac anatomy for coupled electromechanical modeling. Ann. Biomed. Eng. 44, 58–70. doi: 10.1007/s10439-015-1474-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Crozier, A., Blazevic, B., Lamata, P., Plank, G., Ginks, M., Duckett, S., et al. (2016b). The relative role of patient physiology and device optimisation in cardiac resynchronisation therapy: a computational modelling study. J. Mol. Cell. Cardiol. 96, 93–100. doi: 10.1016/j.yjmcc.2015.10.026

PubMed Abstract | CrossRef Full Text | Google Scholar

de Vecchi, A., Gomez, A., Pushparajah, K., Schaeffter, T., Simpson, J., Razavi, R., et al. (2016). A novel methodology for personalized simulations of ventricular hemodynamics from noninvasive imaging data. Comput. Med. Imaging Graph. 51(Suppl. C):20–31. doi: 10.1016/j.compmedimag.2016.03.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Demiray, H. (1972). A note on the elasticity of soft biological tissues. J. Biomech. 5, 309–311. doi: 10.1016/0021-9290(72)90047-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Doenst, T., Spiegel, K., Reik, M., Markl, M., Hennig, J., Nitzsche, S., et al. (2009). Fluid-dynamic modeling of the human left ventricle: methodology and application to surgical ventricular reconstruction. Ann. Thorac. Surg. 87, 1187–1195. doi: 10.1016/j.athoracsur.2009.01.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Donea, J., and Huerta, A. (2003). Finite Element Methods for Flow Problems. Hoboken, NJ: John Wiley & Sons.

Google Scholar

Douglas, J., and Wang, J. P. (1989). An absolutely stabilized finite element method for the stokes problem. Math. Comput. 52, 495–508. doi: 10.1090/S0025-5718-1989-0958871-X

CrossRef Full Text | Google Scholar

Elman, H., Howle, V. E., Shadid, J., Shuttleworth, R., and Tuminaro, R. (2008). A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier–Stokes equations. J. Comput. Phys. 227, 1790–1808. doi: 10.1016/j.jcp.2007.09.026

CrossRef Full Text | Google Scholar

Esmaily Moghadam, M., Bazilevs, Y., Hsia, T.-Y., Vignon-Clementel, I. E., and Marsden, A. L. (2011). A comparison of outlet boundary treatments for prevention of backflow divergence with relevance to blood flow simulations. Comput. Mech. 48, 277–291. doi: 10.1007/s00466-011-0599-0

CrossRef Full Text | Google Scholar

Förster, C., Wall, W. A., and Ramm, E. (2006). On the geometric conservation law in transient flow calculations on deforming domains. Int. J. Numer. Methods Fluids 50, 1369–1379. doi: 10.1002/fld.1093

CrossRef Full Text | Google Scholar

Forti, D. (2016). Parallel Algorithms for the Solution of Large-Scale Fluid-Structure Interaction Problems in Hemodynamics. Ph.D. thesis, Lausanne.

Fouchet-Incaux, J. (2014). Artificial boundaries and formulations for the incompressible Navier–Stokes equations: applications to air and blood flows. SeMA J. 64, 1–40. doi: 10.1007/s40324-014-0012-y

CrossRef Full Text | Google Scholar

Franca, L. P., and Hughes, T. J. (1988). Two classes of mixed finite element methods. Comput. Methods Appl. Mech. Eng. 69, 89–129. doi: 10.1016/0045-7825(88)90168-5

CrossRef Full Text | Google Scholar

Freitag, L. A., and Knupp, P. M. (2002). Tetrahedral mesh improvement via optimization of the element condition number. Int. J. Numer. Methods Eng. 53, 1377–1391. doi: 10.1002/nme.341

CrossRef Full Text | Google Scholar

Goubergrits, L., Mevert, R., Yevtushenko, P., Schaller, J., Kertzscher, U., Meier, S., et al. (2013). The impact of MRI-based inflow for the hemodynamic evaluation of aortic coarctation. Ann. Biomed. Eng. 41, 2575–2587. doi: 10.1007/s10439-013-0879-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Goubergrits, L., Riesenkampff, E., Yevtushenko, P., Schaller, J., Kertzscher, U., Hennemuth, A., et al. (2015). MRI-based computational fluid dynamics for diagnosis and treatment prediction: clinical validation study in patients with coarctation of aorta. J. Magn. Reson. Imaging 41, 909–916. doi: 10.1002/jmri.24639

PubMed Abstract | CrossRef Full Text | Google Scholar

Gresho, P. M., and Sani, R. L. (1987). On pressure boundary conditions for the incompressible Navier–Stokes equations. Int. J. Numer. Methods Fluids 7, 1111–1145. doi: 10.1002/fld.1650071008

CrossRef Full Text | Google Scholar

Guccione, J. M., Costa, K. D., and McCulloch, A. D. (1995). Finite element stress analysis of left ventricular mechanics in the beating dog heart. J. Biomech. 28, 1167–1177. doi: 10.1016/0021-9290(94)00174-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Gurev, V., Pathmanathan, P., Fattebert, J.-L., Wen, H.-F., Magerlein, J., Gray, R. A., et al. (2015). A high-resolution computational model of the deforming human heart. Biomech. Model. Mechanobiol. 14, 829–849. doi: 10.1007/s10237-014-0639-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Hirschvogel, M., Bassilious, M., Jagschies, L., Wildhirt, S. M., and Gee, M. W. (2017). A monolithic 3D-0D coupled closed-loop model of the heart and the vascular system: experiment-based parameter estimation for patient-specific cardiac mechanics. Int. J. Numer. Methods Biomed. Eng. 33:e2842. doi: 10.1002/cnm.2842

CrossRef Full Text | Google Scholar

Hirt, C., Amsden, A. A., and Cook, J. (1974). An arbitrary Lagrangian–Eulerian computing method for all flow speeds. J. Comput. Phys. 14, 227–253. doi: 10.1016/0021-9991(74)90051-5

CrossRef Full Text | Google Scholar

Hughes, T. J. (1995). Multiscale phenomena: Green's functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods. Comput. Methods Appl. Mech. Eng. 127, 387–401.

Google Scholar

Hughes, T. J., Franca, L. P., and Balestra, M. (1986). A new finite element formulation for computational fluid dynamics: V. Circumventing the Babuška–Brezzi condition: a stable Petrov–Galerkin formulation of the Stokes problem accommodating equal-order interpolations. Comput. Methods Appl. Mech. Eng. 59, 85–99. doi: 10.1016/0045-7825(86)90025-3

CrossRef Full Text | Google Scholar

Jansen, K. E., Whiting, C. H., and Hulbert, G. M. (2000). A generalized-α method for integrating the filtered Navier–Stokes equations with a stabilized finite element method. Comput. Methods Appl. Mech. Eng. 190, 305–319. doi: 10.1016/S0045-7825(00)00203-6

CrossRef Full Text | Google Scholar

Kanchi, H., and Masud, A. (2007). A 3D adaptive mesh moving scheme. Int. J. Numer. Methods Fluids 54, 923–944. doi: 10.1002/fld.1512

CrossRef Full Text | Google Scholar

Kelm, M., Goubergrits, L., Bruening, J., Yevtushenko, P., Fernandes, J., Sündermann, S., et al. (2017). Model-based therapy planning allows prediction of haemodynamic outcome after aortic valve replacement. Sci. Rep. 7:9897. doi: 10.1038/s41598-017-03693-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Klotz, S., Dickstein, M. L., and Burkhoff, D. (2007). A computational method of prediction of the end-diastolic pressure–volume relationship by single beat. Nat. Protoc. 2, 2152–2158. doi: 10.1038/nprot.2007.270

PubMed Abstract | CrossRef Full Text | Google Scholar

Krishnamurthy, A., Villongco, C. T., Chuang, J., Frank, L. R., Nigam, V., Belezzuoli, E., et al. (2013). Patient-specific models of cardiac biomechanics. J. Comput. Phys. 244, 4–21. doi: 10.1016/j.jcp.2012.09.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Krittian, S. B., Lamata, P., Michler, C., Nordsletten, D. A., Bock, J., Bradley, C. P., et al. (2012). A finite-element approach to the direct computation of relative cardiovascular pressure from time-resolved MR velocity data. Med. Image Anal. 16, 1029–1037. doi: 10.1016/j.media.2012.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Land, S., Gurev, V., Arens, S., Augustin, C. M., Baron, L., Blake, R., et al. (2015). Verification of cardiac mechanics software: benchmark problems and solutions for testing active and passive material behaviour. Proc. Math. Phys. Eng. Sci. 471:20150641. doi: 10.1098/rspa.2015.0641

CrossRef Full Text | Google Scholar

Long, C., Marsden, A., and Bazilevs, Y. (2013). Fluid–structure interaction simulation of pulsatile ventricular assist devices. Comput. Mech. 52, 971–981. doi: 10.1007/s00466-013-0858-3

CrossRef Full Text | Google Scholar

Lumens, J., Delhaas, T., Kirn, B., and Arts, T. (2009). Three-wall segment (TriSeg) model describing mechanics and hemodynamics of ventricular interaction. Ann. Biomed. Eng. 37, 2234–2255. doi: 10.1007/s10439-009-9774-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Meier, S., Hennemuth, A., Friman, O., Bock, J., Markl, M., and Preusser, T. (2010). “Non-invasive 4D blood flow and pressure quantification in central blood vessels via PC-MRI,” in 2010 Computing in Cardiology, ed A. Murray (Belfast: IEEE), 903–906.

Google Scholar

Mihalef, V., Ionasec, R. I., Sharma, P., Georgescu, B., Voigt, I., Suehling, M., et al. (2011). Patient-specific modelling of whole heart anatomy, dynamics and haemodynamics from four-dimensional cardiac ct images. Interface Focus 1, 286–296. doi: 10.1098/rsfs.2010.0036

PubMed Abstract | CrossRef Full Text | Google Scholar

Mittal, R., Seo, J. H., Vedula, V., Choi, Y. J., Liu, H., Huang, H. H., et al. (2016). Computational modeling of cardiac hemodynamics: current status and future outlook. J. Comput. Phys. 305, 1065–1082. doi: 10.1016/j.jcp.2015.11.022

CrossRef Full Text | Google Scholar

Murray, C. D. (1926). The physiological principle of minimum work I. The vascular system and the cost of blood volume. Proc. Natl. Acad. Sci. U.S.A. 12, 207–214. doi: 10.1073/pnas.12.3.207

PubMed Abstract | CrossRef Full Text | Google Scholar

Neic, A., Campos, F. O., Prassl, A. J., Niederer, S. A., Bishop, M. J., Vigmond, E. J., et al. (2017). Efficient computation of electrograms and ECGs in human whole heart simulations using a reaction-eikonal model. J. Comput. Phys. 346, 191–211. doi: 10.1016/j.jcp.2017.06.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Neic, A., Liebmann, M., Hoetzl, E., Mitchell, L., Vigmond, E. J., Haase, G., et al. (2012). Accelerating cardiac bidomain simulations using graphics processing units. IEEE Trans. Biomed. Eng. 59, 2281–2290. doi: 10.1109/TBME.2012.2202661

PubMed Abstract | CrossRef Full Text | Google Scholar

Nichols, W., O'Rourke, M., and Vlachopoulos, C. (2011). McDonald's Blood Flow in Arteries: Theoretical, Experimental and Clinical Principles. London: CRC press.

Google Scholar

Niederer, S., Mitchell, L., Smith, N., and Plank, G. (2011). Simulating human cardiac electrophysiology on clinical time-scales. Front. Physiol. 2:14. doi: 10.3389/fphys.2011.00014

PubMed Abstract | CrossRef Full Text | Google Scholar

Niederer, S. A., Plank, G., Chinchapatnam, P., Ginks, M., Lamata, P., Rhode, K. S., et al. (2011). Length-dependent tension in the failing heart and the efficacy of cardiac resynchronization therapy. Cardiovasc. Res. 89:336. doi: 10.1093/cvr/cvq318

PubMed Abstract | CrossRef Full Text | Google Scholar

Nordsletten, D., McCormick, M., Kilner, P., Hunter, P., Kay, D., and Smith, N. (2011). Fluid–solid coupling for the investigation of diastolic and systolic human left ventricular function. Int. J. Numer. Methods Biomed. Eng. 27, 1017–1039. doi: 10.1002/cnm.1405

CrossRef Full Text | Google Scholar

Pathmanathan, P., and Whiteley, J. J. P. (2009). A numerical method for cardiac mechanoelectric simulations. Ann. Biomed. Eng. 37, 860–873. doi: 10.1007/s10439-009-9663-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Pauli, L., and Behr, M. (2017). On stabilized space-time FEM for anisotropic meshes: incompressible Navier–Stokes equations and applications to blood flow in medical devices. Int. J. Numer. Methods Fluids 85, 189–209. doi: 10.1002/fld.4378

CrossRef Full Text | Google Scholar

Quarteroni, A., Lassila, T., Rossi, S., and Ruiz-Baier, R. (2017). Integrated heart – coupling multiscale and multiphysics models for the simulation of the cardiac function. Comput. Methods Appl. Mech. Eng. 314, 345–407. doi: 10.1016/j.cma.2016.05.031

CrossRef Full Text | Google Scholar

Ralovich, K., Itu, L., Vitanovski, D., Sharma, P., Ionasec, R., Mihalef, V., et al. (2015). Noninvasive hemodynamic assessment, treatment outcome prediction and follow-up of aortic coarctation from MR imaging. Med. Phys. 42, 2143–2156. doi: 10.1118/1.4914856

PubMed Abstract | CrossRef Full Text | Google Scholar

Rocha, B. M., Kickinger, F., Prassl, A. J., Haase, G., Vigmond, E. J., dos Santos, R. W., et al. (2011). A macro finite-element formulation for cardiac electrophysiology simulations using hybrid unstructured grids. IEEE Trans. Biomed. Eng. 58, 1055–1065. doi: 10.1109/TBME.2010.2064167

PubMed Abstract | CrossRef Full Text | Google Scholar

Schäfer, M., Turek, S., Durst, F., Krause, E., and Rannacher, R. (1996). Benchmark Computations of Laminar Flow Around a Cylinder. Wiesbaden: Vieweg+Teubner Verlag.

Google Scholar

Schenkel, T., Malve, M., Reik, M., Markl, M., Jung, B., and Oertel, H. (2009). MRI-based CFD analysis of flow in a human left ventricle: methodology and application to a healthy heart. Ann. Biomed. Eng. 37, 503–515. doi: 10.1007/s10439-008-9627-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Sellier, M. (2011). An iterative method for the inverse elasto-static problem. J. Fluids Struct. 27, 1461–1470. doi: 10.1016/j.jfluidstructs.2011.08.002

CrossRef Full Text | Google Scholar

Seo, J. H., and Mittal, R. (2013). Effect of diastolic flow patterns on the function of the left ventricle. Phys. Fluids 25:110801. doi: 10.1063/1.4819067

CrossRef Full Text | Google Scholar

Seo, J. H., Vedula, V., Abraham, T., and Mittal, R. (2013). Multiphysics computational models for cardiac flow and virtual cardiography. Int. J. Numer. Methods Biomed. Eng. 29, 850–869. doi: 10.1002/cnm.2556

PubMed Abstract | CrossRef Full Text | Google Scholar

Silvester, D., Elman, H., Kay, D., and Wathen, A. (2001). Efficient preconditioning of the linearized Navier–Stokes equations for incompressible flow. J. Comput. Appl. Math. 128, 261–279. doi: 10.1016/S0377-0427(00)00515-X

CrossRef Full Text | Google Scholar

Stalling, D., Westerhoff, M., and Hege, H.-C. (2005). “Amira: a highly interactive system for visual data analysis,” in The Visualization Handbook, eds C. D. Hansen and C. R. Johnson (London; Oxford; Boston, MA; New York NY; San Diego, CA: Elsevier), 749–767.

Google Scholar

Steinbach, O. (2007). Numerical Approximation Methods for Elliptic Boundary Value Problems: Finite and Boundary Elements. Springer Science & Business Media.

Google Scholar

Su, B., San Tan, R., Le Tan, J., Guo, K. W. Q., Zhang, J. M., Leng, S., et al. (2016). Cardiac MRI based numerical modeling of left ventricular fluid dynamics with mitral valve incorporated. J. Biomech. 49, 1199–1205. doi: 10.1016/j.jbiomech.2016.03.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Tallec, P. L., and Mouro, J. (2001). Fluid structure interaction with large structural displacements. Comput. Methods Appl. Mech. Eng. 190, 3039–3067. doi: 10.1016/S0045-7825(00)00381-9

CrossRef Full Text | Google Scholar

Tang, D., Yang, C., Geva, T., and del Nido, P. J. (2010). Image-based patient-specific ventricle models with fluid-structure interaction for cardiac function assessment and surgical design optimization. Prog. Pediatr. Cardiol. 30, 51–62. doi: 10.1016/j.ppedcard.2010.09.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, D., Yang, C., Geva, T., and Pedro, J. (2008). Patient-specific MRI-based 3D FSI RV/LV/patch models for pulmonary valve replacement surgery and patch optimization. J. Biomech. Eng. 130:041010. doi: 10.1115/1.2913339

CrossRef Full Text | Google Scholar

ten Tusscher, K. H. W. J., Noble, D., Noble, P. J., and Panfilov, A. V. (2004). A model for human ventricular tissue. Am. J. Physiol. Heart Circul. Physiol. 286, H1573–H1589. doi: 10.1152/ajpheart.00794.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

The CGAL Project (2017). CGAL User and Reference Manual, 4.11 Edn. CGAL Editorial Board.

Van Emden, H., and Yang, U. M. (2002). BoomerAMG: a parallel algebraic multigrid solver and preconditioner. Appl. Numer. Math. 41, 155–177. doi: 10.1016/S0168-9274(01)00115-5

CrossRef Full Text | Google Scholar

van Loon, R., Anderson, P., van de Vosse, F., and Sherwin, S. (2007). Comparison of various fluid structure interaction methods for deformable bodies. Comput. Struct. 85, 833–843. doi: 10.1016/j.compstruc.2007.01.010

CrossRef Full Text | Google Scholar

Vázquez, M., Arís, R., Aguado-Sierra, J., Houzeaux, G., Santiago, A., López, M., et al. (2015). Alya Red CCM: HPC-Based Cardiac Computational Modelling. Cham: Springer International Publishing.

Google Scholar

Vigmond, E., Hughes, M., Plank, G., and Leon, L. (2003). Computational tools for modeling electrical activity in cardiac tissue. J. Electrocardiol. 36, 69–74. doi: 10.1016/j.jelectrocard.2003.09.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Vigmond, E. J., Clements, C., McQueen, D. M., and Peskin, C. S. (2008a). Effect of bundle branch block on cardiac output: a whole heart simulation study. Prog. Biophys. Mol. Biol. 97, 520–542. doi: 10.1016/j.pbiomolbio.2008.02.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Vigmond, E. J., Weber dos Santos, R., Prassl, A. J., Deo, M., and Plank, G. (2008b). Solvers for the cardiac bidomain equations. Prog. Biophys. Mol. Biol. 96, 3–18. doi: 10.1016/j.pbiomolbio.2007.07.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., and Liu, W. K. (2004). Extended immersed boundary method using FEM and RKPM. Comput. Methods Appl. Mech. Eng. 193, 1305–1321. doi: 10.1016/j.cma.2003.12.024

CrossRef Full Text | Google Scholar

Westerhof, N., Elzinga, G., and Sipkema, P. (1971). An artificial arterial system for pumping hearts. J. Appl. Physiol. 31, 776–781. doi: 10.1152/jappl.1971.31.5.776

PubMed Abstract | CrossRef Full Text | Google Scholar