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

Elias Karabelas  Institute of Biophysics, Medical University of Graz, Graz,
Austria  Matthias A. F. Gsell  Institute of Biophysics, Medical University of
Graz, Graz, Austria  Christoph M. Augustin  Institute of Biophysics, Medical
University of Graz, Graz, Austria  Shadden Research Group, Department of
Mechanical Engineering, University of California, Berkeley, CA, USA  Laura
Marx  Institute of Biophysics, Medical University of Graz, Graz, Austria
Aurel Neic  Institute of Biophysics, Medical University of Graz, Graz, Austria
Anton J. Prassl  Institute of Biophysics, Medical University of Graz, Graz,
Austria  Leonid Goubergrits  Department of Congenital Heart Disease/Pediatric
Cardiology, German Heart Institute Berlin, Berlin, Germany  Institute for
Imaging Science and Computational Modelling in Cardiovascular Medicine,
Charité-Universitätsmedizin Berlin, Berlin, Germany  Titus Kuehne  Department
of Congenital Heart Disease/Pediatric Cardiology, German Heart Institute
Berlin, Berlin, Germany  Institute for Imaging Science and Computational
Modelling in Cardiovascular Medicine, Charité-Universitätsmedizin Berlin,
Berlin, Germany  Gernot Plank

###### Abstract

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 towards 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.

Keywords: cardiac mechanics, computational fluid dynamics, finite element
model, arbitrary Lagrangian-Eulerian formulation, patient-specific modeling,
translational cardiac modeling, total heart function.

##  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 [1, 34, 66], or by image-based
kinematic models [24, 18, 69, 72, 54, 76] 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 [46] or stenting of a coarctation [35], 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 [20] 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 [22, 62, 79, 78, 83] and immersed boundary (IB) methods [19, 71,
86]. While ALE formulations often rely on severe simplifications or automatic
remeshing strategies [51], 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 [82].

In this study, we report on the progress made towards 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 [5] 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 to
1536 compute cores. The overall cost of the entire workflow for building,
fitting and execution of EMF simulations is $\approx$ 48 hours which is
comparable to plain image-based kinematic driver models [55].

##  2 Methods

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

  1. 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.

  2. 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).

  3. iii)

The solution procedure and software implementation details are outlined in
Section 2.4.

  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\text{\,}\mathrm{mmHg}\text{/}$$
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.

####  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\times 2.5\times$2.5\text{\,}\mathrm{mm}\text{/}$$;
reconstructed voxel $1.7\times 1.7\times$2.5\text{\,}\mathrm{mm}\text{/}$$;
repetition time $3.5\text{\,}\mathrm{ms}\text{/}$; echo time
$2.2\text{\,}\mathrm{ms}\text{/}$; flip angle $None$; 25 reconstructed cardiac
phases; number of signal averages 1; High velocity encoding
($36\text{\,}\mathrm{m}\text{/}\mathrm{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
software111http://www.gyrotools.com/products/gt-flow.html (Gyrotools, Zurich,
Switzerland).

The 3DWH exemplary sequence parameters were: acquired voxel $0.66\times
0.66\times$3.2\text{\,}\mathrm{mm}\text{/}$$; reconstructed voxel $0.66\times
0.66\times$1.6\text{\,}\mathrm{mm}\text{/}$$; repetition time
$4.0\text{\,}\mathrm{ms}\text{/}$; echo time
$2.0\text{\,}\mathrm{ms}\text{/}$; flip angle $None$; 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\times 0.86\times 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 towards 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
MevisFlow222https://www.mevis.fraunhofer.de/en/solutionpages/mevisflow-non-
invasive-interactive-exploration-of-in-vivo-hemodynamics.html. Briefly, the
PPE can be derived from the Navier–Stokes equations by taking the divergence
of the momentum equation (25), see [36, 49] for more details. The processing
and analysis pipeline of the pressure mapping consists of the following four
steps.

  1. i)

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

  2. ii)

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

  3. 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.

  4. iv)

Solving the PPE at each time step having 4D VEC-MRI data as input.
Furthermore, a $5\text{\,}\mathrm{\char 37\relax}\text{/}$ mask size reduction
is applied in order to avoid numerical inconsistencies close to the vessel
wall as suggested earlier [53].

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.10.2\text{\,}\mathrm{mg}\text{/}\mathrm{kg}$, max.
$5\text{\,}\mathrm{mg}\text{/}$), followed by a bolus of propofol
($12\text{\,}\mathrm{mg}\text{/}\mathrm{kg}$, as needed) and continuous
infusion of propofol (approximately
$4\text{\,}\mathrm{m}\mathrm{g}\mathrm{/}\mathrm{k}\mathrm{g}\mathrm{/}\mathrm{h}$,
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
software333https://amira.zib.de [74]. The segmentations were smoothed and
upsampled to a $0.1\text{\,}\mathrm{mm}\text{/}$ isotropic resolution using a
variational smoothing method [20]. The resulting high resolution multi-label
segmentation was meshed using CGAL444http://www.cgal.org [81], giving a global
mesh $\Omega_{\mathrm{s},\mathrm{total}}^{0}$ consisting of tetrahedral
elements. Here, $(\bullet)^{0}$ denotes the mechanical reference configuration
at end-diastolic pressure. The mesh was subdivided into various subdomains
corresponding to predefined labels which are summarized in Table 2. We write

$\Omega_{\mathrm{s},\mathrm{total}}^{0}=\bigcup_{i\in
I}\Omega_{\mathrm{s},i}^{0},$ (1)

with the index set

$I:=\\{\mathrm{lv},\mathrm{ao},\mathrm{cushion},\mathrm{av},\mathrm{mv},\mathrm{lvbp},\mathrm{aobp}\\},$
(2)

see Figure 1 (E–G) for illustration. With this, we define the following
submeshes

$\displaystyle\Omega_{\mathrm{s}}^{0}$
$\displaystyle:=\Omega_{\mathrm{s,total}}^{0}\backslash\left(\Omega_{\mathrm{s,lbvp}}^{0}\cup\Omega_{\mathrm{s,aobp}}^{0}\right),$
(3) $\displaystyle\Omega_{\mathrm{s,bp}}^{0}=\widetilde{\Omega}_{f}^{0}$
$\displaystyle:=\Omega_{\mathrm{s,av}}^{0}\cup\Omega_{\mathrm{s,lvbp}}^{0}\cup\Omega^{0}_{\mathrm{s,aobp}},$
(4)

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

$\displaystyle\Gamma_{\mathrm{s,N}}^{0}$
$\displaystyle:=\partial\left(\left(\Omega_{s,\mathrm{lv}}^{0}\cup\Omega_{\mathrm{s,av}}^{0}\cup\Omega_{\mathrm{s,mv}}^{0}\right)\cap\Omega_{\mathrm{s,lvbp}}^{0}\right),$
(5) $\displaystyle\Gamma_{\mathrm{s,H}}^{0}$
$\displaystyle:=\partial\Omega_{\mathrm{s}}^{0}\backslash\left(\Gamma_{\mathrm{s,N}}^{0}\cup\Gamma_{\mathrm{s,D}}^{0}\right),$
(6) $\displaystyle\Gamma_{\mathrm{s,bp}}^{0}$
$\displaystyle:=\partial\Omega_{\mathrm{s,bp}}^{0}\backslash\Gamma_{\mathrm{s,D}}^{0},$
(7)

where $\Gamma_{\mathrm{s,D}}^{0}$ denote the cutoff faces as indicated by blue
lines in Figure 1; $\Gamma_{\mathrm{s,N}}^{0}$ are surfaces subject to
pressure; and $\Gamma_{\mathrm{s,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_{\mathrm{s,bp}}^{0}$ was extracted and smoothed using the VMTK
toolbox555http://www.vmtk.org [2]. The smoothed surface,
$\Gamma_{\mathrm{f,wall}}^{0}$, was used to define the boundary of the fluid
domain reference configuration, $\Omega_{\mathrm{f}}^{0}$, for volumetric FE
meshing using ANSYS ICEM CFD666http://www.ansys.com/Services/training-
center/platform/introduction-to-ansys-icem-cfd-Hexa. Refined boundary layers
were included in this process to better resolve sharp gradients in the
vicinity of $\Gamma_{\mathrm{f,wall}}^{0}$ occurring during simulation of
hemodynamics. The various processing stages for building EM and CFD models are
illustrated in Figures 1 and 4, respectively.

###  2.2 Electromechanical Model

####  2.2.1 Electrophysiology of the LV

A recently developed reaction-eikonal (R-E) model [58] 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

$\left\\{\begin{array}[]{rcll}\sqrt{\nabla_{\mathbf{X}}t_{\mathrm{a}}^{\top}\,\boldsymbol{\mathbf{V}}\,\nabla_{\mathbf{X}}t_{\mathrm{a}}}&=&1&\text{in
}\Omega_{\mathrm{s,lv}}^{0},\\\ t_{\mathrm{a}}&=&t_{0}&\text{on
}\Gamma_{\mathrm{s},\ast}^{0},\end{array}\right.$ (8)

where $(\nabla_{\mathbf{X}})$ is the gradient with respect to the end-
diastolic reference configuration $\Omega_{\mathrm{s,lv}}^{\mathrm{0}}$;
$t_{\mathrm{a}}$ is a positive function describing the wavefront arrival time
at location $\mathbf{X}\in\Omega_{\mathrm{s,lv}}^{\mathrm{0}}$; and $t_{0}$
are initial activations at locations
$\Gamma_{\mathrm{s},\ast}^{0}\subseteq\Gamma_{\mathrm{s,N}}^{\mathrm{0}}$. The
symmetric positive definite $3\times 3$ tensor
$\boldsymbol{\mathbf{V}}(\mathbf{X})$ holds the squared velocities
$\left(v_{\mathrm{f}}(\mathbf{X}),v_{\mathrm{s}}(\mathbf{X}),v_{\mathrm{n}}(\mathbf{X})\right)$
associated to the tissue’s eigenaxes, referred to as fiber, $\mathbf{f}_{0}$,
sheet, $\mathbf{s}_{0}$, and sheet normal, $\mathbf{n}_{0}$, orientations. The
arrival time function $t_{\mathrm{a}}(\mathbf{X})$ was subsequently used in a
modified monodomain R-D model given as

$\beta C_{\mathrm{m}}\frac{\partial V_{\mathrm{m}}}{\partial
t}=\nabla_{\mathbf{X}}\cdot\boldsymbol{\mathbf{\sigma}}_{\mathrm{i}}\nabla_{\mathbf{X}}V_{\mathrm{m}}+I_{\mathrm{foot}}-\beta
I_{\mathrm{ion}},$ (9)

where an arrival time dependent foot current,
$I_{\mathrm{foot}}(t_{\mathrm{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 [80]. As indicated in Equations (8, 9),
activation sequences and electrical source distribution in the LV were
computed in its end-diastolic configuration
$\Omega_{\mathrm{s,lv}}^{\mathrm{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

$-\nabla_{\mathbf{X}}\cdot\boldsymbol{\mathbf{F}}\boldsymbol{\mathbf{S}}(\mathbf{d}_{\mathrm{s}},t)=\mathbf{0}\quad\mbox{in
}{\Omega}_{\mathrm{s}}^{0},$ (10)

for $t\in[0,T]$, where $\mathbf{d}_{\mathrm{s}}$ is the unknown displacement;
$\boldsymbol{\mathbf{F}}$ is the deformation gradient;
$\boldsymbol{\mathbf{S}}$ is the second Piola–Kirchhoff stress tensor; and
$(\nabla_{\mathbf{X}}\;\cdot)$ denotes the divergence operator in the Lagrange
reference configuration. Homogeneous Dirichlet boundary conditions

$\mathbf{d}_{\mathrm{s}}=\mathbf{0}\quad\text{on}\quad\Gamma_{\mathrm{s,D}}^{0},$
(11)

homogeneous Neumann boundary conditions

$\boldsymbol{\mathbf{F}}\boldsymbol{\mathbf{S}}(\mathbf{d}_{\mathrm{s}},t)\,\mathbf{n}_{\mathrm{s,0}}=\mathbf{n}_{\mathrm{s,0}}\quad\text{on}\quad{\Gamma}_{\mathrm{s,H}}^{0},$
(12)

and inhomogeneous Neumann boundary conditions

$\boldsymbol{\mathbf{F}}\boldsymbol{\mathbf{S}}(\mathbf{d}_{\mathrm{s}},t)\,\mathbf{n}_{\mathrm{s,0}}=p(t)J\,\boldsymbol{\mathbf{F}}^{-\top}(\mathbf{d}_{\mathrm{s}},t)\,\mathbf{n}_{\mathrm{s,0}}\quad\text{on}\quad{\Gamma}_{\mathrm{s,N}}^{0}$
(13)

were imposed, where $\mathbf{n}_{\mathrm{s,0}}$ is the outward unit normal
vector; $p(t)$ is the pressure; and $J=\det\boldsymbol{\mathbf{F}}$. For sake
of clarity, boundary conditions are illustrated in Figure 1 (C).

The total stress $\boldsymbol{\mathbf{S}}$ was additively decomposed according
to

$\boldsymbol{\mathbf{S}}=\boldsymbol{\mathbf{S}}_{\mathrm{pas}}+\boldsymbol{\mathbf{S}}_{\mathrm{act}},$
(14)

where $\boldsymbol{\mathbf{S}}_{\mathrm{pas}}$ and
$\boldsymbol{\mathbf{S}}_{\mathrm{act}}$ refer to the passive and active
stresses, respectively. Passive stresses were modeled based on the
constitutive equation

$\boldsymbol{\mathbf{S}}_{\mathrm{pas}}=2\frac{\partial\Psi(\boldsymbol{\mathbf{C}})}{\partial\boldsymbol{\mathbf{C}}}$
(15)

given a hyper-elastic strain-energy function $\Psi$ and the right Cauchy–Green
strain tensor
$\boldsymbol{\mathbf{C}}=\boldsymbol{\mathbf{F}}^{\top}\boldsymbol{\mathbf{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_{\mathrm{s,lv}}^{0}$ and fiber orientations
$(\mathbf{f}_{0},\mathbf{s}_{0},\mathbf{n}_{0})$ are the same as for the EP
model, Section 2.2.1, the transversely isotropic constitutive relation

$\Psi_{\mathrm{Guc}}(\boldsymbol{\mathbf{C}})=\frac{\kappa}{2}\left(\log\,J\right)^{2}+\frac{C_{\mathrm{Guc}}}{2}\left[\exp(\mathcal{Q})-1\right].$
(16)

by Guccione et al. [38] was employed. Here, the term in the exponent is

$\mathcal{Q}=b_{\mathrm{f}}(\mathbf{f}_{0}\cdot\overline{\boldsymbol{\mathbf{E}}}\mathbf{f}_{0})^{2}+b_{\mathrm{t}}\left[(\mathbf{s}_{0}\cdot\overline{\boldsymbol{\mathbf{E}}}\mathbf{s}_{0})^{2}+(\mathbf{n}_{0}\cdot\overline{\boldsymbol{\mathbf{E}}}\mathbf{n}_{0})^{2}+2(\mathbf{s}_{0}\cdot\overline{\boldsymbol{\mathbf{E}}}\mathbf{n}_{0})^{2}\right]+2b_{\mathrm{fs}}\left[(\mathbf{f}_{0}\cdot\overline{\boldsymbol{\mathbf{E}}}\mathbf{s}_{0})^{2}+(\mathbf{f}_{0}\cdot\overline{\boldsymbol{\mathbf{E}}}\mathbf{n}_{0})^{2}\right]$
(17)

and
$\overline{\boldsymbol{\mathbf{E}}}=\frac{1}{2}(\overline{\boldsymbol{\mathbf{C}}}-\boldsymbol{\mathbf{I}})$
is the modified isochoric Green–Lagrange strain tensor, where
$\overline{\boldsymbol{\mathbf{C}}}:=J^{-\nicefrac{{2}}{{3}}}\boldsymbol{\mathbf{C}}$.
Default values of $b_{\mathrm{f}}=18.48$, $b_{\mathrm{t}}=3.58$, and
$b_{\mathrm{fs}}=1.627$ were used. The parameter $C_{\mathrm{Guc}}$ was varied
for the different cases, see Table 3. In the aorta
$\Omega_{\mathrm{s,ao}}^{0}$, unlike in previous studies [4], 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 [23] was used

$\Psi_{\mathrm{Dem}}(\boldsymbol{\mathbf{C}}):=\frac{\kappa}{2}\left(\log\,J\right)^{2}+\frac{a}{2\,b}\left\\{\exp\left[b\,\big{(}\operatorname{tr}(\overline{\boldsymbol{\mathbf{C}}})-3\big{)}\right]-1\right\\}.$
(18)

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

A simplified phenomenological contractile model was used to represent active
stress generation [61]. Owing to its small number of parameters and its direct
relation to clinically measurable quantities such as peak pressure,
$p_{\mathrm{lv}}$, and the maximum rate of rise of pressure,
$\,\mathrm{d}p_{\mathrm{lv}}/\,\mathrm{d}t_{\mathrm{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

$S_{\mathrm{a}}(t,\lambda)=S_{\mathrm{peak}}\,\phi(\lambda)\,\tanh^{2}\left(\frac{t_{\mathrm{s}}}{\tau_{\mathrm{c}}}\right)\,\tanh^{2}\left(\frac{t_{\mathrm{dur}}-t_{\mathrm{s}}}{\tau_{\mathrm{r}}}\right),\qquad\text{for
}0<t_{\mathrm{s}}<t_{\mathrm{dur}},$ (19)

with

$\phi=\tanh(\mathrm{ld}(\lambda-\lambda_{0})),\quad\tau_{\mathrm{c}}=\tau_{\mathrm{c_{0}}}+\mathrm{ld}_{\mathrm{up}}(1-\phi),\quad
t_{\mathrm{s}}=t-t_{\mathrm{a}}-t_{\mathrm{emd}}$ (20)

and $t_{\mathrm{s}}$ is the onset of contraction; $\phi(\lambda)$ is a non-
linear length-dependent function in which $\lambda$ is the fiber stretch and
$\lambda_{0}$ is the lower limit of fiber stretch below which no further
active tension is generated; $t_{\mathrm{a}}$ is the local activation time
from Eq. (8); $t_{\mathrm{emd}}$ is the EM delay between the onsets of
electrical depolarization and active stress generation; $S_{\mathrm{peak}}$ is
the peak isometric tension; $t_{\mathrm{dur}}$ is the duration of active
stress transient; $\tau_{\mathrm{c}}$ is time constant of contraction;
$\tau_{\mathrm{c_{0}}}$ is the baseline time constant of contraction;
$\mathrm{ld}_{\mathrm{up}}$ is the length-dependence of $\tau_{\mathrm{c}}$;
$\tau_{\mathrm{r}}$ is the time constant of relaxation; and $\mathrm{ld}$ is
the degree of length dependence. Thus, active stresses in this simplified
model are only length-dependent, but dependence on fiber velocity,
$\dot{\lambda}$, is ignored. Unlinke in previous studies [61] we set the
nonlinear length-dependent function $\phi(\lambda)=1$ for the whole
simulation. The active stress tensor in the reference configuration
$\Omega_{\mathrm{s,lv}}^{0}$ induced in fiber direction $\mathbf{f}_{0}$ is
defined as

$\boldsymbol{\mathbf{S}}_{\mathrm{a}}=S_{\mathrm{a}}\left(\mathbf{f}_{0}\cdot\boldsymbol{\mathbf{C}}\mathbf{f}_{0}\right)^{-1}\mathbf{f}_{0}\otimes\mathbf{f}_{0},$
(21)

with $S_{a}$ defined in Equation (19). This active stress involves a scaling
by $\lambda^{2}=\mathbf{f}_{0}\cdot\boldsymbol{\mathbf{C}}\mathbf{f}_{0}$, see
[63] for details.

####  2.2.3 Mechanical and Hemodynamic Afterload Models

Hydrostatic pressures in the LV, $p_{\mathrm{lv}}$, and the proximal aorta,
$p_{\mathrm{ao}}$, were modeled using a 3-element Windkessel model [88], and
the system of PDEs (10) 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_{\mathrm{lv}}$, exceeds the pressure in the
proximal aorta, $p_{\mathrm{ao}}$, and closes at the end of ejection when
$p_{\mathrm{lv}}$ drops below $p_{\mathrm{ao}}$ and the flow $q_{\mathrm{lv}}$
starts to reverse. In its open state the aortic valve was modeled as a linear
resistor, $R_{\mathrm{av}}$, in series with the characteristic impedance of
the aorta, $Z_{\mathrm{c}}$. During ejection, the pressure in the LV was then
computed by the Windkessel equation

$\frac{\,\mathrm{d}p_{\mathrm{lv}}}{\,\mathrm{d}t}=\frac{1}{C}\left(1+\frac{Z_{\mathrm{c}}+R_{\mathrm{av}}}{R}\right)q_{\mathrm{lv}}+(Z_{\mathrm{c}}+R_{\mathrm{av}})\frac{\,\mathrm{d}q_{\mathrm{lv}}}{\,\mathrm{d}t}-\frac{1}{RC}\,p_{\mathrm{lv}},$
(22)

which predicts the rate of change of pressure in the LV as a function of flow
$q_{\mathrm{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.

A similar form of Equation (22) was also used to estimate the pressure in the
aorta, $p_{\mathrm{ao}}$. In this case, there is no additional resistance due
to an outlet valve and hence $R_{\mathrm{av}}$ is omitted. Balancing of the
PDE (10) and the ODE (22) was achieved by recasting Equation (10) as a saddle
point problem, see [39, 40].

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

$\displaystyle
C_{i}\frac{\mathrm{d}p_{\mathrm{d},i}}{\mathrm{d}t}+\frac{p_{\mathrm{d},i}}{R_{i}}$
$\displaystyle=q_{i},$ (23) $\displaystyle p_{\mathrm{wk},i}$
$\displaystyle=Z_{i}q_{i}+p_{\mathrm{d},i},$ (24)

see [31, 13] for more details. During ejection the Windkessel pressure
$p_{\mathrm{wk}}$ at an outlet was then applied as an outflow boundary
condition for the fluid flow model, see Section 2.5.5. In Equations (23, 24),
$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 [59]. Let $\Omega_{\mathrm{f}}\subsetneq\mathbb{R}^{3}$ denote the
fluid domain, then the evolution of flow is governed by the incompressible
Navier–Stokes equations

$\displaystyle\rho_{\mathrm{f}}\left(\frac{\partial{}}{\partial{t}}\mathbf{u}_{\mathrm{f}}+\mathbf{u}_{\mathrm{f}}\cdot\nabla_{\mathbf{x}}\mathbf{u}_{\mathrm{f}}\right)-\nabla_{\mathbf{x}}\cdot\boldsymbol{\mathsf{\sigma}}_{\mathrm{f}}(\mathbf{u}_{\mathrm{f}},p_{\mathrm{f}})$
$\displaystyle=\mathbf{0}$ $\displaystyle\text{in }\Omega_{\mathrm{f}},$ (25)
$\displaystyle\nabla_{\mathbf{x}}\cdot\mathbf{u}_{\mathrm{f}}$
$\displaystyle=0$ $\displaystyle\text{in }\Omega_{\mathrm{f}},$ (26)
$\displaystyle\mathbf{u}_{\mathrm{f}}$ $\displaystyle=\mathbf{0}$
$\displaystyle\text{on }\Gamma_{\text{noslip}},$ (27)
$\displaystyle\mathbf{u}_{\mathrm{f}}$ $\displaystyle=\mathbf{g}_{\mathrm{f}}$
$\displaystyle\text{on }\Gamma_{\text{inflow}},$ (28)
$\displaystyle\boldsymbol{\mathsf{\sigma}}_{\mathrm{f}}\mathbf{n}_{\mathrm{f}}-\rho_{\mathrm{f}}\beta\left(\mathbf{u}_{\mathrm{f}}\cdot\mathbf{n}_{\mathrm{f}}\right)_{-}\mathbf{u}_{\mathrm{f}}$
$\displaystyle=-p_{\mathrm{wk}}\mathbf{n}_{\mathrm{f}}$ $\displaystyle\text{on
}\Gamma_{\text{outflow}},$ (29)
$\displaystyle\mathbf{u}_{\mathrm{f}}\big{|}_{t=0}$
$\displaystyle=\mathbf{u}_{0},$ (30)

where $\mathbf{u}_{\mathrm{f}}$ denotes fluid velocity; $p_{\mathrm{f}}$ is
fluid pressure; $\rho_{\mathrm{f}}$ is the density of blood;
$\boldsymbol{\mathsf{\sigma}}_{\mathrm{f}}$ is the fluid stress tensor;
$\mathbf{g}_{\mathrm{f}}$ is a velocity inlet; $p_{\mathrm{wk}}$ is the
Windkessel pressure solution to Equations (23, 24); $\mathbf{u}_{0}$ refers to
the initial condition; $\mathbf{n}_{\mathrm{f}}$ is the outward normal of the
fluid domain; and $(\nabla_{\mathbf{x}})$ is the gradient and
$(\nabla_{\mathbf{x}}\cdot)$ is the divergence operator in the fluid domain
$\Omega_{\mathrm{f}}$. The sets $\Gamma_{\text{noslip}}$,
$\Gamma_{\text{inflow}}$, and $\Gamma_{\text{outflow}}$ denote the
complementary subsets of $\Gamma_{\text{f}}:=\partial\Omega_{\text{f}}$ and we
assume that $\left\lvert{\Gamma_{\text{outflow}}}\right\rvert>0$. Note that
Equation (28) 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_{\mathrm{wk}}\equiv 0$, boundary condition (29) is referred to as
_directional do-nothing boundary condition_ [28, 16] and the term

$\left(\mathbf{u}_{\mathrm{f}}\cdot\mathbf{n}_{\mathrm{f}}\right)_{-}:=\frac{1}{2}\left(\mathbf{u}_{\mathrm{f}}\cdot\mathbf{n}_{\mathrm{f}}-\left\lvert{\mathbf{u}_{\mathrm{f}}\cdot\mathbf{n}_{\mathrm{f}}}\right\rvert\right)$
(31)

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\leq\frac{1}{2}$ were also used without causing
numerical issues, see [28]. All physical parameters in Equations (25)–(30) are
summarized in Table 4. In presence of multiple outlets outflow boundary
conditions as given in Equation (29) were prescribed at each of the outlets.

####  2.3.1 Extension to Moving Geometries

For time-dependent fluid domains, i.e.,
$\Omega_{\text{f}}=\Omega_{\text{f}}^{t}$, Equations (25)–(30) 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 [12,
Section 1.3] and [41]. Similar to structural mechanics, a reference fluid
configuration $\Omega_{\text{f}}^{0}\subsetneq\mathbb{R}^{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 $\mathbf{x}$
and the reference coordinate system is denoted by $\mathbf{X}$. Their relation
is given by the ALE mapping
$\mathbf{x}=\mathbf{X}+\mathbf{d}_{\mathrm{f}}(t,\mathbf{X})$. Here,
$\mathbf{d}_{\mathrm{f}}(t,\mathbf{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

$\Omega_{\text{f}}^{t}:=\left\\{\mathbf{x}:\mathbf{x}=\mathbf{X}+\mathbf{d}_{\mathrm{f}}(t,\mathbf{X}),\,\forall\mathbf{X}\in\Omega_{\text{f}}^{0}\right\\}.$
(32)

Further, we define the fluid domain velocity $\mathbf{w}_{\mathrm{f}}$ as

$\mathbf{w}_{\mathrm{f}}:=\frac{\partial{}}{\partial{t}}\mathbf{d}_{\mathrm{f}}\big{|}_{\mathbf{X}},$
(33)

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

$\Gamma_{\mathrm{f,mov}}^{t}:=\partial\Omega_{\text{f}}^{t}\
\backslash\bigcup_{i=1}^{n_{\text{outlets}}}\Gamma_{\mathrm{f,outflow},i}^{t},$
(34)

where $\Gamma_{\mathrm{f,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., [12, 29],

$\displaystyle\rho_{\mathrm{f}}\left(\frac{\partial{}}{\partial{t}}\mathbf{u}_{\mathrm{f}}\big{|}_{\mathbf{X}}+\left(\mathbf{u}_{\mathrm{f}}-\mathbf{w}_{\mathrm{f}}\right)\cdot\nabla_{\mathbf{x}}\mathbf{u}_{\mathrm{f}}\right)-\nabla_{\mathbf{x}}\cdot\boldsymbol{\mathsf{\sigma}}_{\mathrm{f}}(\mathbf{u}_{\mathrm{f}},p_{\mathrm{f}})$
$\displaystyle=\mathbf{0}$ $\displaystyle\text{on }\Omega_{\mathrm{f}}^{t},$
(35) $\displaystyle\nabla_{\mathbf{x}}\cdot\mathbf{u}_{\mathrm{f}}$
$\displaystyle=0$ $\displaystyle\text{on }\Omega_{\mathrm{f}}^{t},$ (36)
$\displaystyle\mathbf{u}_{\mathrm{f}}$
$\displaystyle=\mathbf{g}_{\mathrm{mov}}$ $\displaystyle\text{on
}\Gamma_{\mathrm{f,mov}}^{t},$ (37)
$\displaystyle\boldsymbol{\mathbf{\sigma}}_{\mathrm{f}}(\mathbf{u}_{\mathrm{f}},p_{\mathrm{f}})\mathbf{n}_{\mathrm{f}}-\rho_{\mathrm{f}}\beta((\mathbf{u}_{\mathrm{f}}-\mathbf{w}_{\mathrm{f}})\cdot\mathbf{n}_{\mathrm{f}})_{-}\mathbf{u}_{\mathrm{f}}$
$\displaystyle=-p_{\mathrm{wk,i}}\mathbf{n}_{\mathrm{f}}$
$\displaystyle\text{on each }\Gamma_{\mathrm{f,outflow},i}^{t},$ (38)
$\displaystyle\mathbf{u}_{\mathrm{f}}\big{|}_{t=0}$
$\displaystyle=\mathbf{u}_{0}$ $\displaystyle\text{in
}\Omega_{\mathrm{f}}^{0}.$ (39)

Along $\Gamma_{\mathrm{f,mov}}^{t}$ we imposed equality between fluid velocity
and the velocity of the moving surfaces. Boundary condition (38) is the ALE
equivalent of the outflow stabilization in Equation (29), see [12, 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 [11, 12, 64], the discrete variational formulation of the ALE
equations (35)–(39) can be stated in the following abstract form: find
$\mathbf{u}^{\mathrm{h}}_{\mathrm{f}}\in[\mathcal{S}_{\mathrm{h},\mathbf{g}}^{1}(\mathcal{T}_{\mathrm{N}})]^{3},p^{\mathrm{h}}_{\mathrm{f}}\in\mathcal{S}_{\mathrm{h}}^{1}(\mathcal{T}_{\mathrm{N}})$
such that for all
$\mathbf{v}^{\mathrm{h}}\in[\mathcal{S}_{h,\mathbf{0}}^{1}(\mathcal{T}_{\mathrm{N}})]^{3}$
and for all
$q^{\mathrm{h}}\in\mathcal{S}_{\mathrm{h}}^{1}(\mathcal{T}_{\mathrm{N}})$

$A_{\mathrm{NS}}(\mathbf{v}^{\mathrm{h}},q^{\mathrm{h}};\mathbf{u}^{\mathrm{h}}_{\mathrm{f}},p^{\mathrm{h}}_{\mathrm{f}})+S_{\mathrm{VMS}}(\mathbf{v}^{\mathrm{h}},q^{\mathrm{h}};\mathbf{u}^{\mathrm{h}}_{\mathrm{f}},p^{\mathrm{h}}_{\mathrm{f}})=F_{\mathrm{NS}}(\mathbf{v}^{\mathrm{h}}),$
(40)

with the classical bilinear form of the Navier–Stokes equations

$\displaystyle
A_{\mathrm{NS}}(\mathbf{v}^{\mathrm{h}},q^{\mathrm{h}};\mathbf{u}^{\mathrm{h}}_{\mathrm{f}},p^{\mathrm{h}}_{\mathrm{f}}):=$
$\displaystyle\rho_{\mathrm{f}}\int\limits_{\Omega_{\mathrm{f}}^{t}}\mathbf{v}^{\mathrm{h}}\cdot\left(\frac{\partial{}}{\partial{t}}\mathbf{u}^{\mathrm{h}}_{\mathrm{f}}+\left(\mathbf{u}_{\mathrm{f}}^{\mathrm{h}}-\mathbf{w}_{\mathrm{f}}^{\mathrm{h}}\right)\cdot\nabla_{\mathbf{x}}\mathbf{u}_{\mathrm{f}}^{\mathrm{h}}\right)\,\mathrm{d}\mathbf{x}+\int\limits_{\Omega_{\mathrm{f}}^{t}}\boldsymbol{\mathsf{\varepsilon}}(\mathbf{v}^{\mathrm{h}}):\boldsymbol{\mathbf{\sigma}}_{\mathrm{f}}(\mathbf{u}^{\mathrm{h}}_{\mathrm{f}},p^{\mathrm{h}}_{\mathrm{f}})\,\mathrm{d}\mathbf{x}$
$\displaystyle+\int\limits_{\Omega_{\mathrm{f}}^{t}}q^{\mathrm{h}}\nabla_{\mathbf{x}}\cdot\mathbf{u}^{\mathrm{h}}_{\mathrm{f}}\,\mathrm{d}\mathbf{x}-\rho_{\mathrm{f}}\beta\sum_{i=1}^{n_{\text{outlets}}}\int\limits_{\Gamma_{\mathrm{f,outflow},i}^{t}}((\mathbf{u}^{\mathrm{h}}_{\mathrm{f}}-\mathbf{w}^{\mathrm{h}}_{\mathrm{f}})\cdot\mathbf{n}_{\mathrm{f}})_{-}\mathbf{v}^{\mathrm{h}}\cdot\mathbf{u}^{\mathrm{h}}_{\mathrm{f}}\,\mathrm{d}s_{\mathbf{x}},$
(41)

the bilinear form $S_{\mathrm{VMS}}$, which is explained later in Equation
(44), and the right-hand side contribution

$F_{\text{NS}}(\mathbf{v}^{\mathrm{h}}):=-\sum_{i=1}^{n_{\text{outlets}}}p_{\mathrm{wk,i}}\int\limits_{\Gamma_{\mathrm{f,outflow},i}^{t}}\mathbf{v}^{\mathrm{h}}\cdot\mathbf{n}_{\mathrm{f}}\,\mathrm{d}s_{\mathbf{x}}.$
(42)

In Equation (41), $\boldsymbol{\mathbf{\varepsilon}}$ is the strain-rate
tensor and $\mathbf{w}^{\mathrm{h}}_{\mathrm{f}}$ is the discrete counterpart
of the fluid domain velocity $\mathbf{w}_{\mathrm{f}}$, i.e.,

$\mathbf{w}^{\mathrm{h}}_{\mathrm{f}}(t^{n+1},\mathbf{X})=\frac{\mathbf{d}_{\mathrm{f}}(t^{n+1},\mathbf{X})-\mathbf{d}_{\mathrm{f}}(t^{n},\mathbf{X})}{\Delta
t}.$ (43)

The FE function space
$\mathcal{S}_{\mathrm{h},*}^{1}(\mathcal{T}_{\mathrm{N}})$ is the conformal
trial space of piecewise linear, globally continuous basis functions over a
decomposition $\mathcal{T}_{\mathrm{N}}$ of $\Omega_{\mathrm{f}}^{t}$ into $N$
simplicial elements constrained by $\mathbf{v}^{\mathrm{h}}=\ast$ on essential
boundaries. The FE function space
$\mathcal{S}_{\mathrm{h}}^{1}(\mathcal{T}_{\mathrm{N}})$ denotes the same
space without constraints. For further details we refer to [17, 75].

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 [25]. 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
[43, 32, 26, 14]. 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 [42, 11, 12, 64] 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 [11].
The term $S_{\mathrm{VMS}}$ in Equation (40) denotes the bilinear form of the
RBVMS formulation and reads as

$\displaystyle
S_{\mathrm{VMS}}(\mathbf{v}^{\mathrm{h}},q^{\mathrm{h}};\mathbf{u}^{\mathrm{h}}_{\mathrm{f}},p^{\mathrm{h}}_{\mathrm{f}}):=$
$\displaystyle\frac{1}{\rho_{\mathrm{f}}}\sum\limits_{l=1}^{n_{\text{el}}}\int\limits_{\tau_{\ell}}\tau_{\text{MOM}}\left(\rho_{\mathrm{f}}\left(\mathbf{u}^{\mathrm{h}}_{\mathrm{f}}-\mathbf{w}^{\mathrm{h}}_{\mathrm{f}}\right)\cdot\nabla_{\mathbf{x}}\mathbf{v}^{\mathrm{h}}+q^{\mathrm{h}}\right)\cdot\mathbf{r}_{\text{MOM}}(\mathbf{u}^{\mathrm{h}}_{\mathrm{f}},p^{\mathrm{h}}_{\mathrm{f}})\,\mathrm{d}\mathbf{x}$
$\displaystyle+\sum\limits_{l=1}^{n_{\text{el}}}\int\limits_{\tau_{\ell}}\tau_{\text{CONT}}\nabla_{\mathbf{x}}\cdot\mathbf{v}^{\mathrm{h}}\nabla_{\mathbf{x}}\cdot\mathbf{u}^{\mathrm{h}}_{\mathrm{f}}\,\mathrm{d}\mathbf{x}$
$\displaystyle-\sum\limits_{l=1}^{n_{\text{el}}}\int\limits_{\tau_{\ell}}\tau_{\text{MOM}}\mathbf{v}^{\mathrm{h}}\cdot\left(\nabla_{\mathbf{x}}\mathbf{u}^{\mathrm{h}}_{\mathrm{f}}\mathbf{r}_{\text{MOM}}(\mathbf{u}^{\mathrm{h}}_{\mathrm{f}},p^{\mathrm{h}}_{\mathrm{f}})\right)\,\mathrm{d}\mathbf{x}$
$\displaystyle-\frac{1}{\rho_{\mathrm{f}}}\sum\limits_{l=1}^{n_{\text{el}}}\int\limits_{\tau_{\ell}}\tau_{\text{MOM}}^{2}\boldsymbol{\mathsf{\varepsilon}}(\mathbf{v}^{\mathrm{h}}):(\mathbf{r}_{\text{MOM}}(\mathbf{u}^{\mathrm{h}}_{\mathrm{f}},p^{\mathrm{h}}_{\mathrm{f}})\otimes\mathbf{r}_{\text{MOM}}(\mathbf{u}^{\mathrm{h}}_{\mathrm{f}},p^{\mathrm{h}}_{\mathrm{f}}))\,\mathrm{d}\mathbf{x},$
(44)

where the vector $\mathbf{r}_{\text{MOM}}$ is defined as

$\mathbf{r}_{\text{MOM}}(\mathbf{u}^{\mathrm{h}}_{\mathrm{f}},p^{\mathrm{h}}_{\mathrm{f}}):=\rho_{\mathrm{f}}\left(\frac{\partial{}}{\partial{t}}\mathbf{u}^{\mathrm{h}}_{\mathrm{f}}+\left(\mathbf{u}^{\mathrm{h}}_{\mathrm{f}}-\mathbf{w}^{\mathrm{h}}_{\mathrm{f}}\right)\cdot\nabla_{\mathbf{x}}\mathbf{u}^{\mathrm{h}}_{\mathrm{f}}\right)-\nabla_{\mathbf{x}}\cdot\boldsymbol{\mathsf{\sigma}}_{\mathrm{f}}(\mathbf{u}^{\mathrm{h}}_{\mathrm{f}},p^{\mathrm{h}}_{\mathrm{f}}).$
(45)

The definition of the parameters $\tau_{\text{MOM}},\tau_{\text{CONT}}$
according to [64] is given by

$\tau_{\text{MOM}}:=\min\left\\{\left(\frac{4}{\Delta
t^{2}}+(\mathbf{u}^{\mathrm{h}}_{\mathrm{f}}-\mathbf{w}^{\mathrm{h}}_{\mathrm{f}})\cdot\boldsymbol{\mathsf{G}}(\mathbf{u}^{\mathrm{h}}_{\mathrm{f}}-\mathbf{w}^{\mathrm{h}}_{\mathrm{f}})\right)^{-\frac{1}{2}},\frac{\rho_{\mathrm{f}}C_{\mathrm{M}}}{\mu_{\mathrm{f}}\sqrt{\boldsymbol{\mathsf{G}}:\boldsymbol{\mathsf{G}}}}\right\\},$
(46)

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

$\boldsymbol{\mathsf{K}}:=\frac{1}{2\sqrt[3]{2}}\begin{pmatrix}3&-1&-1\\\
-1&3&-1\\\ -1&-1&3\end{pmatrix}$ (47)

and the constant $C_{\mathrm{M}}=0.0285$. Further, the stabilization parameter
$\tau_{\text{CONT}}$ is defined as

$\displaystyle\tau_{\text{CONT}}$
$\displaystyle:=\frac{1}{\tau_{\text{MOM}}\mathbf{g}_{\mathrm{f}}\cdot\mathbf{g}_{\mathrm{f}}},$
(48) $\displaystyle g_{\mathrm{f},i}$
$\displaystyle:=\sum\limits_{j=1}^{3}\left(\frac{\partial{\mathbf{\xi}}}{\partial{\mathbf{x}}}\right)_{ji}.$
(49)

####  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{g}_{\mathrm{mov}}=\frac{\partial{}}{\partial{t}}\mathbf{d}_{\mathrm{s}}$
in Equation (37). Since the surface of the reference CFD blood pool mesh,
$\partial\Omega_{\mathrm{f}}^{0}$, is not conformal with the surface of the
reference EM blood pool mesh, $\Omega_{\mathrm{s,bp}}^{0}$, and the overlap of
the two surfaces is imperfect due to smoothing of
$\partial\Omega_{\mathrm{f}}^{0}$ and remeshing of $\Omega_{\mathrm{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{d}}_{\mathrm{s}}$ that form
the endocardial interface with the blood pool, $\Gamma_{\mathrm{s,bp}}^{0}$,
were extracted from the solution $\mathbf{d}_{\mathrm{s}}$ defined at
$\Omega_{\mathrm{s}}^{0}$. Since the mesh interface between
$\Omega_{\mathrm{s}}^{0}$ and $\Omega_{\mathrm{s,bp}}^{0}$ is conformal the
extracted displacements can be applied as inhomogeneous time-varying Dirichlet
boundary conditions to the blood pool mesh $\Omega_{\mathrm{s,bp}}^{0}$ to
solve a linear elastic problem given as

$\displaystyle-\nabla_{\mathbf{X}}\cdot\boldsymbol{\mathsf{\sigma}}(\mathbf{d}_{\mathrm{s}}(t))$
$\displaystyle=\mathbf{0}$ $\displaystyle\text{in
}\Omega_{s,\mathrm{bp}}^{0},$ (50) $\displaystyle\mathbf{d}_{\mathrm{s}}(t)$
$\displaystyle=\widetilde{\mathbf{d}}_{\mathrm{s}}(t)$ $\displaystyle\text{on
}\partial\Omega_{s,\mathrm{bp}}^{0},$ (51)

where stress and strain tensor are

$\displaystyle\boldsymbol{\mathsf{\sigma}}(\mathbf{d}_{\mathrm{s}})$
$\displaystyle:=\frac{E}{1+\nu}\left(\frac{\nu}{1-2\nu}\nabla_{\mathbf{X}}\cdot\mathbf{d}_{\mathrm{s}}\boldsymbol{\mathsf{I}}+\boldsymbol{\mathbf{\varepsilon}}(\mathbf{d}_{\mathrm{s}})\right),$
(52) $\displaystyle\boldsymbol{\mathbf{\varepsilon}}(\mathbf{d}_{\mathrm{s}})$
$\displaystyle:=\frac{1}{2}\left(\nabla_{\mathbf{X}}\mathbf{d}_{\mathrm{s}}+(\nabla_{\mathbf{X}}\mathbf{d}_{\mathrm{s}})^{\top}\right),$
(53)

the constant $E$ is Young’s modulus in $\mathrm{kPa}\text{/}$ and the constant
$\nu$ is Poisson’s ratio which is dimensionless in the range of $[-1,0.5)$.
Combining the solutions $\mathbf{d}_{\mathrm{s}}$ computed for
$\Omega_{\mathrm{s}}^{0}$ and $\Omega_{\mathrm{s,bp}}^{0}$ yields
displacements $\mathbf{d}_{\mathrm{s}}$ for $\Omega_{\mathrm{s,total}}^{0}$.
Since $\partial\Omega_{\mathrm{f}}^{0}$ is fully embedded in this domain,
$\Omega_{\mathrm{s,total}}^{0}$ $\Omega_{\mathrm{s,total}}^{0}$ can be used as
a hanging background mesh for interpolating displacements onto the blood pool
mesh, $\Omega_{\mathrm{f}}^{0}$, used for CFD simulations. However, for
reasons of mesh quality, interpolation is solely applied on the boundary
$\Omega_{\mathrm{f}}^{0}$ itself, and to find the interior displacement field
the exact same linear elastic problem (50)–(53), is solved for
$\mathbf{d}_{\mathrm{f}}$ instead of $\mathbf{d}_{\mathrm{s}}$. In both
patient cases studied, ejection fractions were large leading to a substantial
deformation of the blood pool mesh $\Omega_{\mathrm{f}}^{t}$. To maintain mesh
quality under such large deformations the parameters $E$ and $\nu$ governing
stiffness and incompressibility of the material were altered accordingly.
Initially, a fixed $E_{0}$ and $\nu_{0}$ was chosen while the subsequent
modification of $E$ and $\nu$ was guided by a combination of the two following
strategies.

  1. i)

_Quality based stiffening:_ For each element $\tau_{\ell}$ in the fluid mesh a
tetrahedral quality indicator $\kappa(\tau_{\ell})$ based on the movement from
the previous time step was calculated, see [33, 45], and rescaled such that
for elements of good quality $\kappa$ is close to 1, while for elements with
poor quality $\kappa$ tends towards infinity. Eventually, the parameter $E$
was multiplied by $\kappa$ within each element.

  2. ii)

_$\nu$ -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 (50)–(53), the current element volumes were tracked and a
volume ratio relative to an undeformed reference element was computed as
$\frac{\left\lvert{\tau_{\ell}}\right\rvert}{\left\lvert{\hat{\tau}_{\ell}}\right\rvert}$.
For ratios below a predefined critical value the parameter $\nu$ 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. [84]. Numerical details on FE discretization [67]
and solution of EP [85, 57, 58] and EM [6] 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 [30, Section 1.4.2], showed similar results to the more
expensive fully-implicit generalized-$\alpha$ method [44]. Hence, all results
in Section 3 were obtained using the semi-implicit scheme; to advance from
time step $t^{n}$ to $t^{n+1}$, 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 [8, 9, 7]. We used a right
preconditoned flexible GMRES method with PETSc fieldsplit preconditioning [73,
27] which in turn uses BoomerAMG [89] to approximate sub-block inverses. While
the time step size for mechanics and CFD was the same, $\Delta
t_{\mathrm{mech}}=\Delta t_{\mathrm{CFD}}=$0.5\text{\,}\mathrm{ms}\text{/}$$,
it was significantly smaller for EP, where $\Delta
t_{\mathrm{EP}}=$25\text{\,}\mathrm{\SIUnitSymbolMicro s}\text{/}$$.

The implementation of the CFD solvers has been subjected to various validation
procedures against standard CFD benchmarks [68]. 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 [5],
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
[38]. 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 [5]. The passive biomechanical model
governed by the strain-energy function given in Equation (16) 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 [47] by adjusting the isotropic scaling
parameter $C_{\mathrm{Guc}}$ in Equation (16). 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 [70, 15, 48]. The unloading procedure was repeated with
varying trial material parameters, $C_{\mathrm{Guc}}$, until the difference
between the unstressed LV volume of the model and the prediction of the Klotz
relation was less than $5\text{\,}\mathrm{\char 37\relax}\text{/}$.

####  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 [39] and the parameters of
the active stress given in Equation (19) rate of contraction,
$\tau_{\mathrm{c}}$, and peak active stress, $S_{\mathrm{peak}}$, were
manually adjusted to fit the maximum rate of rise of pressure,
$(\mathrm{d}P/\mathrm{d}t)_{\mathrm{max}}$, and peak pressure,
$p_{\mathrm{lv}}$.

####  2.5.4 Afterload

When the LV pressure $p_{\mathrm{lv}}$ exceeded the aortic pressure,
$p_{\mathrm{ao}}$, ejection was initiated by connecting the LV model with the
lumped 3-element Windkessel model [88]. Volume traces recorded from a given
patient during ejection were used as input to compute aortic pressure traces
by solving Equation (22). 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_{\mathrm{lv}}(t)$ with the instant of opening of the aortic valve in the
pressure trace $p_{\mathrm{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_{\mathrm{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_{\mathrm{c}}$, as well as
resistance, $R$, and compliance, $C$, of the arterial system was sampled using
a recently developed stochastic sampling approach [21].

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
$\pm 20\text{\,}\mathrm{\char 37\relax}\text{/}$ around the mean. Due to high
frequency errors introduced by the pressure transducer we refrained from
computing norms $||p_{\mathrm{ao,meas}}-p_{\mathrm{ao,fit}}||$ to quantify the
deviations of fitted from measured pressure and opted for manual selection
using three criteria, aortic peak pressure, $p_{\mathrm{ao}}$, closing
pressure of aortic valve and exponential decay of $p_{\mathrm{ao}}$ during
diastole. For the sake of fitting $Z_{\mathrm{c}}$ we assumed
$p_{\mathrm{ao}}\approx p_{\mathrm{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,
$\mathbf{d}_{\mathrm{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 $\mathbf{d}_{\mathrm{f}}(t,\mathbf{x})$ defined on the whole CFD
mesh. Figure 4 (G) shows a summary of the boundary conditions. On the boundary
$\Gamma_{\mathrm{f},\mathrm{mov}}^{t}$ a Dirichlet boundary condition
enforcing the mesh velocity $\mathbf{w}^{\mathrm{h}}_{\mathrm{f}}$ is applied.
On each aortic outlet $\Gamma_{\mathrm{f,outflow},i}(t)$ a 3-Element
Windkessel model as described in Section 2.2.3 is attached. Further, the
stabilization parameter $\beta$ in Equation (38) 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 $\alpha_{i}$ of total CO running through
the outlet the resistance $R_{i}$ was estimated as

$R_{i}\approx\frac{\mathrm{MAP}}{\alpha_{i}\mathrm{CO}}.$ (54)

The percentages $\alpha_{i}$ were obtained either by measurement or by
applying Murray’s law [56]. The impedances $Z_{i}$ were chosen as
$5\text{\,}\mathrm{\char 37\relax}\text{/}$ of $R_{i}$, and the compliances
$C_{i}$ were chosen such that
$R_{i}C_{i}\approx$1000\text{\,}\mathrm{ms}\text{/}$$. 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.

##  3 Results

###  3.1 Building electromechanical kinematic driver models

Using a previously developed automated workflow [20], 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.

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 3. 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 5.

###  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
$\mathbf{d}_{\mathrm{s}},\mathbf{d}_{\mathrm{f}}$ in Panel (II) of Figure 6.
Due to the large EF of about $65\text{\,}\mathrm{\char 37\relax}\text{/}$ for
both _28-Pre_ and _44-Pre_ , the blood pool underwent a significant
deformation. However, using a combination of element quality and $\nu$-Volume
based stiffening with an initial Young’s Modulus
$E_{0}=$100\text{\,}\mathrm{kPa}\text{/}$$ and Poisson’s ratio $\nu_{0}=0.3$,
sufficient element quality was preserved throughout the entire ejection phase
and numerical instabilities could be avoided. Panel (I) of Figure 6 shows the
$80^{\mathrm{th}}$-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 [33, 45] 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.

###  3.3 Numerical CFD benchmarks

The implementation of the Navier–Stokes solver was verified by solving a set
of standardized benchmark problems, see [68]. 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 $1536$. Details on computational complexity and
costs are summarized in Table 6. For temporal discretization a time step of
$\Delta t=$0.5\text{\,}\mathrm{ms}\text{/}$$ was used to simulate the ejection
phase lasting for $208\text{\,}\mathrm{ms}\text{/}$. 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\text{\,}\mathrm{\char 37\relax}\text{/}$
at $96$ cores down to $55\text{\,}\mathrm{\char 37\relax}\text{/}$ at $768$
cores. Scalability stalled when doubling the core count to $1536$ which
reduced the degrees of freedom per parallel partition down to $3386$. Parallel
efficiency dropped to $27\text{\,}\mathrm{\char 37\relax}\text{/}$ which is
attributed due to the unfavorable ratio between local compute work and
communication.

###  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 $\varepsilon_{\mathrm{DSC}}$,
$\varepsilon_{\mathrm{BCA}}$, $\varepsilon_{\mathrm{LCA}}$ and
$\varepsilon_{\mathrm{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
$\alpha_{i}$ expressed in percent of the total peak flow through the aorta
ascendens as determined over the plane $\varepsilon_{\mathrm{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 Subfigures (A), and (E) of Figure 8.
MAP and computed mean flow through each outlet orifice were used to determine
the parameters of the coupled Windkessel models of afterload in Equations (23,
24), see Tables 7 and 8. In the _28-Pre_ case this resulted in flow splits of
$\alpha_{i}\approx$23\text{\,}\mathrm{\char
37\relax}\text{/}$,$51.3\text{\,}\mathrm{\char
37\relax}\text{/}$,$12.83\text{\,}\mathrm{\char 37\relax}\text{/}$$ and
$12.83\text{\,}\mathrm{\char 37\relax}\text{/}$ whereas in the _44-Pre_ case
the flow split ratios were $\alpha_{i}\approx$5.68\text{\,}\mathrm{\char
37\relax}\text{/}$$, $57.45\text{\,}\mathrm{\char 37\relax}\text{/}$ and
$34.01\text{\,}\mathrm{\char 37\relax}\text{/}$ for
$\varepsilon_{\mathrm{DSC}}$, $\varepsilon_{\mathrm{BCA}}$,
$\varepsilon_{\mathrm{LCA}}$ and $\varepsilon_{\mathrm{LSCA}}$, respectively.

For the CFD analysis a time step of $\Delta
t=$0.5\text{\,}\mathrm{ms}\text{/}$$ was used. The ejection phases of the EM
simulations were chosen as time horizons for the CFD simulation which lasted
from $t=$90\text{\,}\mathrm{ms}\text{/}$$ to
$t=$302\text{\,}\mathrm{ms}\text{/}$$ in the _28-Pre_ case and from
$t=$70\text{\,}\mathrm{ms}\text{/}$$ to $t=$329\text{\,}\mathrm{ms}\text{/}$$
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 Tables 7 and 8. Pressure $p_{\mathrm{f}}$ along the
centerline $s_{\rm c}$ and fluxes through the planes
$\varepsilon_{\mathrm{DSC}}$, $\varepsilon_{\mathrm{LSC}}$,
$\varepsilon_{\mathrm{BCA}}$, and $\varepsilon_{\mathrm{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
Subfigure (D) of Figure 8) and 4D VEC MRI fluxes. For case _28-Pre_ pressure
drops were calculated from the pressure values on the intersection of the
centerline and $\varepsilon_{\mathrm{DSC}}$, $\varepsilon_{\mathrm{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 $\approx$29.2\text{\,}\mathrm{mmHg}\text{/}$$ which agreed
well with the clinically estimated pressure drop of
$\approx$30\text{\,}\mathrm{mmHg}\text{/}$$. 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 9.
Subfigures (C), (G), and (H) of Figure 8 show velocity profiles at peak flow
condtions.

###  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 3. Only $S_{\mathrm{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 Subsection 3.4 for both scenarios. For
_28-Pre_ we calculated a pressure drop of
$\approx$29.2\text{\,}\mathrm{mmHg}\text{/}$$ while for _28-Post_ a pressure
drop of $\approx$14.15\text{\,}\mathrm{mmHg}\text{/}$$ was calculated.

##  4 Discussion

In this study, we report on the progress made towards 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 [37]. 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
[55], suggesting that the proposed methodology may be, in principle,
compatible with clinical time scales.

###  4.1 Biomechanical modeling versus 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 $\times$
$1.5\text{\,}\mathrm{mm}\text{/}$, but more limited through-plane resolution
(typically about $8\text{\,}\mathrm{mm}\text{/}$) while CT is capable of
isotropic spatial resolution on the sub millimeter scale
($\approx$0.5\text{\,}\mathrm{mm}\text{/}$$) and clear delineation of
trabeculae and lumen boundaries. CMR has the advantage of higher temporal
resolution ($3050\text{\,}\mathrm{ms}\text{/}$) while temporal resolution in
CT depends on the scanning system ($50200\text{\,}\mathrm{ms}\text{/}$). This
is orders-of-magnitude lower than the temporal resolution required for the
flow simulation ($\approx 1000$ 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, $\mathrm{O}_{2}$
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\text{\,}\mathrm{\char 37\relax}\text{/}$$, that is, EF was
$\approx$65\text{\,}\mathrm{\char 37\relax}\text{/}$$ 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 [19, 71, 86] are known to
be more convenient to cope with the large deformation of the CFD blood pool
[65]. 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 [87]
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 [82, 77]. However,
it is often assumed that unstructured FE approaches, as implemented in this
paper, critically depend on automatic remeshing strategies [51] to keep mesh
quality within acceptable bounds [55]. 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\text{\,}\mathrm{\char
37\relax}\text{/}$$ (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 2506987 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 $\approx 17$ iterations during early
ejection up to $\approx 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 [60] and mechanics [6]. 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 [58]. In this study, this recent reaction-eikonal
approach was used for simulating EM using the same FE grid with an average
resolution of $\approx 1$ mm for both EP and mechanics. Such lower resolutions
suffice for solving for mechanics with sufficient accuracy [50]. 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
$\approx$180\text{\,}\mathrm{min}\text{/}$$ 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 $350000/1500000$ nodes/degrees of freedom for
case _28-Pre_ and $400000/1700000$ 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 $1300000/5000000$
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 1536
cores. Local compute load with 1536 was $900/2600$ nodes/degrees of freedom
per core. The patient simulations were performed using 384 cores, resulting in
a load per core of about $900/2700$ nodes/dofs, respectively. At these
resolutions CFD simulations were executed in
$\approx$40\text{\,}\mathrm{min}\text{/}$$, 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
[4]. 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 $\Delta V$ of its lumen. In the CFD simulations $\Delta
V\approx 0$ translates into a stiff aorta of low compliance which may cause a
bias towards 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
[3, 52]. 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 [65]. As a first step towards 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_{\mathrm{f}}}{2}\left\lvert{\mathbf{u}_{\mathrm{f}}}\right\rvert^{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.

## 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.

## 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;
AJP 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.

## References

  * Andersson et al. [2017] Magnus Andersson, Jonas Lantz, Tino Ebbers, and Matts Karlsson.  Multidirectional WSS disturbances in stenotic turbulent flows: A pre-and post-intervention study in an aortic coarctation.  _Journal of Biomechanics_ , 51:8–16, 2017.
  * Antiga et al. [2008] Luca Antiga, Marina Piccinelli, Lorenzo Botti, Bogdan Ene-Iordache, Andrea Remuzzi, and David A Steinman.  An image-based modeling framework for patient-specific computational hemodynamics.  _Medical & Biological Engineering & Computing_, 46(11):1097, 2008.
  * Arts et al. [2005] Theo Arts, Tammo Delhaas, Peter Bovendeerd, Xander Verbeek, and Frits Prinzen.  Adaptation to mechanical load determines shape and properties of heart and circulation: the CircAdapt model.  _Am. J. Physiol. Heart Circ. Physiol._ , 288:1943–1954, 2005.
  * Augustin et al. [2014] Christoph M. Augustin, Gerhard A. Holzapfel, and Olaf Steinbach.  Classical and all-floating FETI methods for the simulation of arterial tissues.  _International Journal for Numerical Methods in Engineering_ , 99(4):290–312, 2014.
  * Augustin et al. [2016a] Christoph M. Augustin, Andrew Crozier, Aurel Neic, Anton J. Prassl, Elias Karabelas, Tiago Ferreira da Silva, Joao F. Fernandes, Fernando Campos, Titus Kuehne, and Gernot Plank.  Patient-specific modeling of left ventricular electromechanics as a driver for haemodynamic analysis.  _EP Europace_ , 18(Supplement 4):iv121–iv129, 2016a.
  * Augustin et al. [2016b] Christoph M Augustin, Aurel Neic, Manfred Liebmann, Anton J Prassl, Steven A Niederer, Gundolf Haase, and Gernot Plank.  Anatomically accurate high resolution modeling of human whole heart electromechanics: A strongly scalable algebraic multigrid solver method for nonlinear deformation.  _Journal of Computational Physics_ , 305:622–646, 2016b.
  * Balay et al. [1997] Satish Balay, William D. Gropp, Lois Curfman McInnes, and Barry F. Smith.  Efficient management of parallelism in object oriented numerical software libraries.  In E. Arge, A. M. Bruaset, and H. P. Langtangen, editors, _Modern Software Tools in Scientific Computing_ , pages 163–202. Birkhäuser Press, 1997.
  * Balay et al. [2016a] Satish Balay, Shrirang Abhyankar, Mark F. Adams, Jed Brown, Peter Brune, Kris Buschelman, Lisandro Dalcin, Victor Eijkhout, William D. Gropp, Dinesh Kaushik, Matthew G. Knepley, Lois Curfman McInnes, Karl Rupp, Barry F. Smith, Stefano Zampini, Hong Zhang, and Hong Zhang.  PETSc Web page.  http://www.mcs.anl.gov/petsc, 2016a.
  * Balay et al. [2016b] Satish Balay, Shrirang Abhyankar, Mark F. Adams, Jed Brown, Peter Brune, Kris Buschelman, Lisandro Dalcin, Victor Eijkhout, William D. Gropp, Dinesh Kaushik, Matthew G. Knepley, Lois Curfman McInnes, Karl Rupp, Barry F. Smith, Stefano Zampini, Hong Zhang, and Hong Zhang.  PETSc users manual.  Technical Report ANL-95/11 - Revision 3.7, Argonne National Laboratory, 2016b.
  * Bayer et al. [2012] J. D. Bayer, R. C. Blake, G. Plank, and N. A. Trayanova.  A novel rule-based algorithm for assigning myocardial fiber orientation to computational heart models.  _Annals of Biomedical Engineering_ , 40(10):2243–2254, 2012.
  * Bazilevs et al. [2007] Y Bazilevs, VM Calo, JA Cottrell, TJR Hughes, A Reali, and G Scovazzi.  Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows.  _Computer Methods in Applied Mechanics and Engineering_ , 197(1):173–201, 2007.
  * Bazilevs et al. [2013] Yuri Bazilevs, Kenji Takizawa, and Tayfun E Tezduyar.  _Computational fluid-structure interaction: methods and applications_.  John Wiley & Sons, 2013.
  * Bertoglio et al. [2017] Cristóbal Bertoglio, Alfonso Caiazzo, Yuri Bazilevs, Malte Braack, Mahdi Esmaily, Volker Gravemeier, Alison L. Marsden, Olivier Pironneau, Irene E. Vignon-Clementel, and Wolfgang A. Wall.  Benchmark problems for numerical treatment of backflow at open boundaries.  _International Journal for Numerical Methods in Biomedical Engineering_ , page e2918, 2017.
  * Bochev et al. [2006] Pavel B Bochev, Clark R Dohrmann, and Max D Gunzburger.  Stabilization of low-order mixed finite elements for the stokes equations.  _SIAM Journal on Numerical Analysis_ , 44(1):82–101, 2006.
  * Bols et al. [2013] Joris Bols, Joris Degroote, Bram Trachet, Benedict Verhegghe, Patrick Segers, and Jan Vierendeels.  A computational method to assess the in vivo stresses and unloaded configuration of patient-specific blood vessels.  _Journal of Computational and Applied Mathematics_ , 246:10–17, 2013.
  * Braack et al. [2014] Malte Braack, Piotr Boguslaw Mucha, and Wojciech Zajaczkowski.  Directional do-nothing condition for the Navier–Stokes equations.  _Journal of Computational Mathematics_ , 32(5):507–521, 2014.
  * Brenner and Scott [2007] Susanne Brenner and Ridgway Scott.  _The mathematical theory of finite element methods_ , volume 15.  Springer Science & Business Media, 2007.
  * Chnafa et al. [2014] C. Chnafa, S. Mendez, and F. Nicoud.  Image-based large-eddy simulation in a realistic left heart.  _Comput. Fluids_ , 94:173–187, 2014.  ISSN 00457930.
  * Choi et al. [2015] Young Joon Choi, Jason Constantino, Vijay Vedula, Natalia Trayanova, and Rajat Mittal.  A new MRI-based model of heart function with coupled hemodynamics and application to normal and diseased canine left ventricles.  _Frontiers in Bioengineering and Biotechnology_ , 3:140, 2015.
  * Crozier et al. [2016a] A Crozier, C M Augustin, A Neic, A J Prassl, M Holler, T E Fastl, A Hennemuth, K Bredies, T Kuehne, M J Bishop, S A Niederer, and G Plank.  Image-based personalization of cardiac anatomy for coupled electromechanical modeling.  _Annals of Biomedical Engineering_ , 44:58–70, 2016a.
  * Crozier et al. [2016b] Andrew Crozier, Bojan Blazevic, Pablo Lamata, Gernot Plank, Matthew Ginks, Simon Duckett, Manav Sohal, Anoop Shetty, Christopher A Rinaldi, Reza Razavi, et al.  The relative role of patient physiology and device optimisation in cardiac resynchronisation therapy: A computational modelling study.  _Journal of Molecular and Cellular Cardiology_ , 96:93–100, 2016b.
  * de Vecchi et al. [2016] A. de Vecchi, A. Gomez, K. Pushparajah, T. Schaeffter, J.M. Simpson, R. Razavi, G.P. Penney, N.P. Smith, and D.A. Nordsletten.  A novel methodology for personalized simulations of ventricular hemodynamics from noninvasive imaging data.  _Computerized Medical Imaging and Graphics_ , 51(Supplement C):20 – 31, 2016.
  * Demiray [1972] Hilmi Demiray.  A note on the elasticity of soft biological tissues.  _Journal of Biomechanics_ , 5(3):309–311, 1972\.
  * Doenst et al. [2009] Torsten Doenst, Kathrin Spiegel, Michael Reik, Michael Markl, Jürgen Hennig, Stefan Nitzsche, Friedhelm Beyersdorf, and Herbert Oertel.  Fluid-dynamic modeling of the human left ventricle: methodology and application to surgical ventricular reconstruction.  _The Annals of Thoracic Surgery_ , 87(4):1187–1195, 2009.
  * Donea and Huerta [2003] Jean Donea and Antonio Huerta.  _Finite element methods for flow problems_.  John Wiley & Sons, 2003.
  * Douglas and Wang [1989] Jim Douglas and Jun Ping Wang.  An absolutely stabilized finite element method for the stokes problem.  _Mathematics of Computation_ , 52(186):495–508, 1989.
  * Elman et al. [2008] Howard Elman, Victoria E Howle, John Shadid, Robert Shuttleworth, and Ray Tuminaro.  A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier–Stokes equations.  _Journal of Computational Physics_ , 227(3):1790–1808, 2008.
  * Esmaily Moghadam et al. [2011] Mahdi Esmaily Moghadam, Yuri Bazilevs, Tain-Yen Hsia, Irene E. Vignon-Clementel, and Alison L. Marsden.  A comparison of outlet boundary treatments for prevention of backflow divergence with relevance to blood flow simulations.  _Computational Mechanics_ , 48(3):277–291, 2011\.
  * Förster et al. [2006] Ch Förster, Wolfgang A Wall, and Ekkehard Ramm.  On the geometric conservation law in transient flow calculations on deforming domains.  _International Journal for Numerical Methods in Fluids_ , 50(12):1369–1379, 2006.
  * Forti [2016] Davide Forti.  _Parallel Algorithms for the Solution of Large-Scale Fluid-Structure Interaction Problems in Hemodynamics_.  PhD thesis, SB, Lausanne, 2016.
  * Fouchet-Incaux [2014] Justine Fouchet-Incaux.  Artificial boundaries and formulations for the incompressible Navier–Stokes equations: applications to air and blood flows.  _SeMA Journal_ , 64(1):1–40, 2014.
  * Franca and Hughes [1988] Leopoldo P Franca and Thomas JR Hughes.  Two classes of mixed finite element methods.  _Computer Methods in Applied Mechanics and Engineering_ , 69(1):89–129, 1988.
  * Freitag and Knupp [2002] Lori A Freitag and Patrick M Knupp.  Tetrahedral mesh improvement via optimization of the element condition number.  _International Journal for Numerical Methods in Engineering_ , 53(6):1377–1391, 2002.
  * Goubergrits et al. [2013] L Goubergrits, R Mevert, P Yevtushenko, J Schaller, U Kertzscher, S Meier, S Schubert, E Riesenkampff, and T Kuehne.  The impact of MRI-based inflow for the hemodynamic evaluation of aortic coarctation.  _Annals of Biomedical Engineering_ , 41(12):2575–2587, 2013.
  * Goubergrits et al. [2015] Leonid Goubergrits, Eugenie Riesenkampff, Pavlo Yevtushenko, Jens Schaller, Ulrich Kertzscher, Anja Hennemuth, Felix Berger, Stephan Schubert, and Titus Kuehne.  MRI-based computational fluid dynamics for diagnosis and treatment prediction: Clinical validation study in patients with coarctation of aorta.  _Journal of Magnetic Resonance Imaging_ , 41(4):909–916, 2015.
  * Gresho and Sani [1987] Philip M Gresho and Robert L Sani.  On pressure boundary conditions for the incompressible Navier–Stokes equations.  _International Journal for Numerical Methods in Fluids_ , 7(10):1111–1145, 1987.
  * Gsell et al. [2018] M A F Gsell, C M Augustin, A J Prassl, J F Fernandes, M Kelm, L Goubergrits, T Kuehne, and G Plank.  Assessment of wall stresses and internal mechanical heart power in the left ventricle: Finite element modeling versus Laplace analysis.  _International Journal for Numerical Methods in Biomedical Engineering_ , 2018.
  * Guccione et al. [1995] Julius M. Guccione, Kevin D. Costa, and Andrew D. McCulloch.  Finite element stress analysis of left ventricular mechanics in the beating dog heart.  _Journal of Biomechanics_ , 28(10):1167 – 1177, 1995.
  * Gurev et al. [2015] Viatcheslav Gurev, Pras Pathmanathan, Jean-Luc Fattebert, Hui-Fang Wen, John Magerlein, Richard A Gray, David F Richards, and J Jeremy Rice.  A high-resolution computational model of the deforming human heart.  _Biomechanics and Modeling in Mechanobiology_ , 14:829–849, 2015.
  * Hirschvogel et al. [2017] Marc Hirschvogel, Marina Bassilious, Lasse Jagschies, Stephen M. Wildhirt, and Michael W. Gee.  A monolithic 3D-0D coupled closed-loop model of the heart and the vascular system: Experiment-based parameter estimation for patient-specific cardiac mechanics.  _International Journal for Numerical Methods in Biomedical Engineering_ , 33(8):e2842, 2017.
  * Hirt et al. [1974] CW Hirt, Anthony A Amsden, and JL Cook.  An arbitrary Lagrangian–Eulerian computing method for all flow speeds.  _Journal of Computational Physics_ , 14(3):227–253, 1974.
  * Hughes [1995] Thomas JR Hughes.  Multiscale phenomena: Green’s functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods.  _Computer Methods in Applied Mechanics and Engineering_ , 127(1-4):387–401, 1995.
  * Hughes et al. [1986] Thomas JR Hughes, Leopoldo P Franca, and Marc Balestra.  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.  _Computer Methods in Applied Mechanics and Engineering_ , 59(1):85–99, 1986.
  * Jansen et al. [2000] Kenneth E Jansen, Christian H Whiting, and Gregory M Hulbert.  A generalized-$\alpha$ method for integrating the filtered Navier–Stokes equations with a stabilized finite element method.  _Computer Methods in Applied Mechanics and Engineering_ , 190(3):305–319, 2000.
  * Kanchi and Masud [2007] Harish Kanchi and Arif Masud.  A 3D adaptive mesh moving scheme.  _International Journal for Numerical Methods in Fluids_ , 54(6-8):923–944, 2007.
  * Kelm et al. [2017] M Kelm, L Goubergrits, J Bruening, P Yevtushenko, JF Fernandes, SH Sündermann, F Berger, T Kuehne, S Nordmeyer, et al.  Model-based therapy planning allows prediction of haemodynamic outcome after aortic valve replacement.  _Scientific Reports_ , 7:9897, 2017.
  * Klotz et al. [2007] Stefan Klotz, Marc L Dickstein, and Daniel Burkhoff.  A computational method of prediction of the end-diastolic pressure–volume relationship by single beat.  _Nature Protocols_ , 2(9):2152–2158, 2007.
  * Krishnamurthy et al. [2013] Adarsh Krishnamurthy, Christopher T Villongco, Joyce Chuang, Lawrence R Frank, Vishal Nigam, Ernest Belezzuoli, Paul Stark, David E Krummen, Sanjiv Narayan, Jeffrey H Omens, et al.  Patient-specific models of cardiac biomechanics.  _Journal of Computational Physics_ , 244:4–21, 2013.
  * Krittian et al. [2012] Sebastian B.S. Krittian, Pablo Lamata, Christian Michler, David A. Nordsletten, Jelena Bock, Chris P. Bradley, Alex Pitcher, Philip J. Kilner, Michael Markl, and Nic P. Smith.  A finite-element approach to the direct computation of relative cardiovascular pressure from time-resolved MR velocity data.  _Medical Image Analysis_ , 16(5):1029 – 1037, 2012.
  * Land et al. [2015] Sander Land, Viatcheslav Gurev, Sander Arens, Christoph M Augustin, Lukas Baron, Robert Blake, Chris Bradley, Sebastian Castro, Andrew Crozier, Marco Favino, Thomas E Fastl, Thomas Fritz, Hao Gao, Alessio Gizzi, Boyce E Griffith, Daniel E Hurtado, Rolf Krause, Xiaoyu Luo, Martyn P Nash, Simone Pezzuto, Gernot Plank, Simone Rossi, Daniel Ruprecht, Gunnar Seemann, Nicolas P Smith, Joakim Sundnes, J Jeremy Rice, Natalia Trayanova, Dafang Wang, Zhinuo Jenny Wang, and Steven A Niederer.  Verification of cardiac mechanics software: benchmark problems and solutions for testing active and passive material behaviour.  _Proceedings. Mathematical, Physical, and Engineering Sciences_ , 471:20150641, 2015.
  * Long et al. [2013] CC Long, AL Marsden, and Y Bazilevs.  Fluid–structure interaction simulation of pulsatile ventricular assist devices.  _Computational Mechanics_ , 52(5):971–981, 2013\.
  * Lumens et al. [2009] Joost Lumens, Tammo Delhaas, Borut Kirn, and Theo Arts.  Three-wall segment (TriSeg) model describing mechanics and hemodynamics of ventricular interaction.  _Ann. Biomed. Eng._ , 37(11):2234–55, 2009.
  * Meier et al. [2010] S. Meier, A. Hennemuth, O. Friman, J. Bock, M. Markl, and T. Preusser.  Non-invasive 4D blood flow and pressure quantification in central blood vessels via PC-MRI.  In A. Murray, editor, _2010 Computing in Cardiology_ , pages 903–906, 2010.
  * Mihalef et al. [2011] Viorel Mihalef, Razvan Ioan Ionasec, Puneet Sharma, Bogdan Georgescu, Ingmar Voigt, Michael Suehling, and Dorin Comaniciu.  Patient-specific modelling of whole heart anatomy, dynamics and haemodynamics from four-dimensional cardiac ct images.  _Interface Focus_ , 1(3):286–296, 2011.
  * Mittal et al. [2016] Rajat Mittal, Jung Hee Seo, Vijay Vedula, Young J Choi, Hang Liu, H Howie Huang, Saurabh Jain, Laurent Younes, Theodore Abraham, and Richard T George.  Computational modeling of cardiac hemodynamics: current status and future outlook.  _Journal of Computational Physics_ , 305:1065–1082, 2016\.
  * Murray [1926] Cecil D Murray.  The physiological principle of minimum work I. The vascular system and the cost of blood volume.  _Proceedings of the National Academy of Sciences_ , 12(3):207–214, 1926.
  * Neic et al. [2012] A Neic, M Liebmann, E Hoetzl, L Mitchell, E J Vigmond, G Haase, and G Plank.  Accelerating cardiac bidomain simulations using graphics processing units.  _IEEE Transactions on Bio-medical Engineering_ , 59:2281–2290, 2012.
  * Neic et al. [2017] Aurel Neic, Fernando O. Campos, Anton J. Prassl, Steven A. Niederer, Martin J. Bishop, Edward J. Vigmond, and Gernot Plank.  Efficient computation of electrograms and ECGs in human whole heart simulations using a reaction-eikonal model.  _Journal of Computational Physics_ , 346:191–211, 2017\.
  * Nichols et al. [2011] Wilmer Nichols, Michael O’Rourke, and Charalambos Vlachopoulos.  _McDonald’s blood flow in arteries: theoretical, experimental and clinical principles_.  CRC press, 2011.
  * Niederer et al. [2011a] Steven Niederer, Lawrence Mitchell, Nicolas Smith, and Gernot Plank.  Simulating human cardiac electrophysiology on clinical time-scales.  _Frontiers in Physiology_ , 2:14, 2011a.
  * Niederer et al. [2011b] Steven A. Niederer, Gernot Plank, Phani Chinchapatnam, Matthew Ginks, Pablo Lamata, Kawal S. Rhode, Christopher A. Rinaldi, Reza Razavi, and Nicolas P. Smith.  Length-dependent tension in the failing heart and the efficacy of cardiac resynchronization therapy.  _Cardiovascular Research_ , 89(2):336, 2011b.
  * Nordsletten et al. [2011] D Nordsletten, M McCormick, PJ Kilner, Peter Hunter, D Kay, and NP Smith.  Fluid–solid coupling for the investigation of diastolic and systolic human left ventricular function.  _International Journal for Numerical Methods in Biomedical Engineering_ , 27(7):1017–1039, 2011.
  * Pathmanathan and Whiteley [2009] Pras Pathmanathan and JP Jonathan P. Whiteley.  A numerical method for cardiac mechanoelectric simulations.  _Ann. Biomed. Eng._ , 37(5):860–73, 2009.
  * Pauli and Behr [2017] L. Pauli and M. Behr.  On stabilized space-time FEM for anisotropic meshes: Incompressible Navier–Stokes equations and applications to blood flow in medical devices.  _International Journal for Numerical Methods in Fluids_ , pages 189–209, 2017.
  * Quarteroni et al. [2017] Alfio Quarteroni, Toni Lassila, Simone Rossi, and Ricardo Ruiz-Baier.  Integrated Heart – Coupling multiscale and multiphysics models for the simulation of the cardiac function.  _Comput. Methods Appl. Mech. Eng._ , 314(314):345–407, 2017.
  * Ralovich et al. [2015] Kristóf Ralovich, Lucian Itu, Dime Vitanovski, Puneet Sharma, Razvan Ionasec, Viorel Mihalef, Waldemar Krawtschuk, Yefeng Zheng, Allen Everett, Giacomo Pongiglione, et al.  Noninvasive hemodynamic assessment, treatment outcome prediction and follow-up of aortic coarctation from MR imaging.  _Medical Physics_ , 42(5):2143–2156, 2015.
  * Rocha et al. [2011] Bernardo M. Rocha, Ferdinand Kickinger, Anton J. Prassl, Gundolf Haase, Edward J. Vigmond, Rodrigo Weber dos Santos, Sabine Zaglmayr, and Gernot Plank.  A macro finite-element formulation for cardiac electrophysiology simulations using hybrid unstructured grids.  _IEEE Trans Biomed Eng_ , 58(4):1055–1065, 2011\.
  * Schäfer et al. [1996] M. Schäfer, S. Turek, F. Durst, E. Krause, and R. Rannacher.  _Benchmark Computations of Laminar Flow Around a Cylinder_ , pages 547–566.  Vieweg+Teubner Verlag, Wiesbaden, 1996.
  * Schenkel et al. [2009] Torsten Schenkel, Mauro Malve, Michael Reik, Michael Markl, Bernd Jung, and Herbert Oertel.  MRI-based CFD analysis of flow in a human left ventricle: methodology and application to a healthy heart.  _Annals of Biomedical Engineering_ , 37(3):503–515, 2009.
  * Sellier [2011] M. Sellier.  An iterative method for the inverse elasto-static problem.  _Journal of Fluids and Structures_ , 27(8):1461 – 1470, 2011.
  * Seo and Mittal [2013] Jung Hee Seo and Rajat Mittal.  Effect of diastolic flow patterns on the function of the left ventricle.  _Physics of Fluids_ , 25(11):110801, 2013.
  * Seo et al. [2013] Jung Hee Seo, Vijay Vedula, Theodore Abraham, and Rajat Mittal.  Multiphysics computational models for cardiac flow and virtual cardiography.  _International Journal for Numerical Methods in Biomedical Engineering_ , 29(8):850–869, 2013.
  * Silvester et al. [2001] David Silvester, Howard Elman, David Kay, and Andrew Wathen.  Efficient preconditioning of the linearized Navier–Stokes equations for incompressible flow.  _Journal of Computational and Applied Mathematics_ , 128(1):261–279, 2001.
  * Stalling et al. [2005] Detlev Stalling, Malte Westerhoff, and Hans-Christian Hege.  Amira: a highly interactive system for visual data analysis.  In Charles D. Hansen and Chris R. Johnson, editors, _The Visualization Handbook_ , pages 749–767. Elsevier, 2005.
  * Steinbach [2007] Olaf Steinbach.  _Numerical approximation methods for elliptic boundary value problems: finite and boundary elements_.  Springer Science & Business Media, 2007.
  * Su et al. [2016] Boyang Su, Ru San Tan, Ju Le Tan, Kenneth Wei Qiang Guo, Jun Mei Zhang, Shuang Leng, Xiaodan Zhao, John Carson Allen, and Liang Zhong.  Cardiac MRI based numerical modeling of left ventricular fluid dynamics with mitral valve incorporated.  _Journal of Biomechanics_ , 49(7):1199–1205, 2016\.
  * Tallec and Mouro [2001] P. Le Tallec and J. Mouro.  Fluid structure interaction with large structural displacements.  _Computer Methods in Applied Mechanics and Engineering_ , 190(24):3039 – 3067, 2001.
  * Tang et al. [2008] Dalin Tang, Chun Yang, Tal Geva, and J Pedro.  Patient-specific MRI-based 3D FSI RV/LV/patch models for pulmonary valve replacement surgery and patch optimization.  _Journal of Biomechanical Engineering_ , 130(4):041010, 2008.
  * Tang et al. [2010] Dalin Tang, Chun Yang, Tal Geva, and Pedro J. del Nido.  Image-based patient-specific ventricle models with fluid–structure interaction for cardiac function assessment and surgical design optimization.  _Prog. Pediatr. Cardiol._ , 30(1-2):51–62, 2010\.
  * ten Tusscher et al. [2004] K H W J ten Tusscher, D Noble, P J Noble, and A V Panfilov.  A model for human ventricular tissue.  _American Journal of Physiology. Heart and Circulatory Physiology_ , 286:H1573–H1589, 2004.
  * The CGAL Project [2017] The CGAL Project.  _CGAL User and Reference Manual_.  CGAL Editorial Board, 4.11 edition, 2017.  URL http://doc.cgal.org/4.11/Manual/packages.html.
  * van Loon et al. [2007] R. van Loon, P.D. Anderson, F.N. van de Vosse, and S.J. Sherwin.  Comparison of various fluid–structure interaction methods for deformable bodies.  _Computers & Structures_, 85(11):833 – 843, 2007.
  * Vázquez et al. [2015] M. Vázquez, R. Arís, J. Aguado-Sierra, G. Houzeaux, A. Santiago, M. López, P. Córdoba, M. Rivero, and J. C. Cajas.  _Alya Red CCM: HPC-Based Cardiac Computational Modelling_ , pages 189–207.  Springer International Publishing, Cham, 2015.
  * Vigmond et al. [2003] E. Vigmond, M. Hughes, G. Plank, and L. Leon.  Computational tools for modeling electrical activity in cardiac tissue.  _J Electrocardiol_ , 36:69–74, 2003.
  * Vigmond et al. [2008a] E J Vigmond, R Weber dos Santos, A J Prassl, M Deo, and G Plank.  Solvers for the cardiac bidomain equations.  _Progress in Biophysics and Molecular Biology_ , 96:3–18, 2008a.
  * Vigmond et al. [2008b] Edward J. Vigmond, Clyde Clements, David M. McQueen, and Charles S. Peskin.  Effect of bundle branch block on cardiac output: A whole heart simulation study.  _Progress in Biophysics and Molecular Biology_ , 97(2):520 – 542, 2008b.
  * Wang and Liu [2004] Xiaodong Wang and Wing Kam Liu.  Extended immersed boundary method using FEM and RKPM.  _Computer Methods in Applied Mechanics and Engineering_ , 193(12):1305 – 1321, 2004.
  * Westerhof et al. [1971] Nicolaas Westerhof, GIJS Elzinga, and P Sipkema.  An artificial arterial system for pumping hearts.  _Journal of Applied Physiology_ , 31(5):776–781, 1971.
  * Yang et al. [2002] Ulrike Meier Yang et al.  BoomerAMG: a parallel algebraic multigrid solver and preconditioner.  _Applied Numerical Mathematics_ , 41(1):155–177, 2002.

## Tables

| Sex | Age | EDV | ESV | SV | EF | HR | CO | $P_{\mathrm{ao/cuff}}^{\mathrm{dia}}$ | $P_{\mathrm{ao/cuff}}^{\mathrm{sys}}$ | MAP | $P_{\mathrm{open}}$
---|---|---|---|---|---|---|---|---|---|---|---|---
|  |  | $\mathrm{ml}\text{/}$ | $\mathrm{ml}\text{/}$ | $\mathrm{ml}\text{/}$ | $\mathrm{\char 37\relax}\text{/}$ | $\mathrm{bpm}$ | $\mathrm{ml}\text{/}\mathrm{s}$ | $\mathrm{mmHg}\text{/}$ | $\mathrm{mmHg}\text{/}$ | $\mathrm{mmHg}\text{/}$ | $\mathrm{mmHg}\text{/}$
_28-Pre_ | F | 9 | 88.2 | 30.6 | 57.54 | 65.3 | 91 | 87.46 | $\nicefrac{{71.1}}{{62}}$ | $\nicefrac{{122.7}}{{138}}$ | $\nicefrac{{88.3}}{{87.3}}$ | $71.33$
_44-Pre_ | M | 12 | 91.7 | 31.6 | 60.09 | 65.5 | 76 | 76.31 | $\nicefrac{{74.6}}{{120}}$ | $\nicefrac{{125.2}}{{154}}$ | $\nicefrac{{91.5}}{{131.3}}$ | $74.78$
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_{\mathrm{ao/cuff,dia}}$ and $P_{\mathrm{ao/cuff,sys}}$), mean arterial pressure (MAP) computed from pressure recorded invasively in the aorta or estimated from $P_{\mathrm{cuff,dia}}$ and $P_{\mathrm{cuff,sys}}$, and aortic valve open pressure $P_{\mathrm{open}}$ determined from invasive pressure recordings.  Tag | Label
---|---
$\mathrm{lv}$ | Myocardium
$\mathrm{ao}$ | Aortic wall
$\mathrm{cushion}$ | Elastic cushion
$\mathrm{av}$ | Aortic valve
$\mathrm{mv}$ | Mitral valve
$\mathrm{lvbp}$ | Left ventricular bloodpool
$\mathrm{aobp}$ | Aortic bloodpool
Table 2: Labels used for defining the subdomains of $\Omega_{s,\mathrm{total}}^{0}$. | EM Fitting
---|---
| $S_{\mathrm{peak}}$ | $t_{\mathrm{dur}}$ | $\tau_{\mathrm{c}0}$ | $\tau_{\mathrm{r}}$ | $t_{\mathrm{emd}}$ | $C_{\mathrm{Guc}}$ | $R$ | $Z$ | $C$
| $\mathrm{kPa}\text{/}$ | $\mathrm{ms}\text{/}$ | $-$ | $-$ | $\mathrm{ms}\text{/}$ | $\mathrm{kPa}\text{/}$ | $\mathrm{kPa}\text{\,}\mathrm{ms}\text{/}\mathrm{ml}$ | $\mathrm{kPa}\text{\,}\mathrm{ms}\text{/}\mathrm{ml}$ | $\mathrm{ml}\text{/}\mathrm{kPa}$
_28-Pre_ | 60.0 | 380 | 30.0 | 30.0 | 15.0 | 0.48 | 170.65 | 12.00 | 6.75
_28-Post_ | 55.0 | 380 | 30.0 | 30.0 | 15.0 | 0.48 | 170.65 | 12.00 | 6.75
_44-Pre_ | 90.0 | 400 | 50.0 | 50.0 | 15.0 | 0.48 | 166.65 | 13.33 | 7.42
Table 3: Fitted parameters for EM Model. Name | Variable | Units |  Expression / Value
---|---|---|---
Velocity | $\mathbf{u}_{\mathrm{f}}$ | $\mathrm{m}\text{/}\mathrm{s}$ | $-$
Pressure | $p_{\mathrm{f}}$ | $\mathrm{Pa}\text{/}$ | $-$
Fluid stress tensor | $\boldsymbol{\mathsf{\sigma}}_{\mathrm{f}}(\mathbf{u}_{\mathrm{f}},p_{\mathrm{f}})$ | $\mathrm{Pa}\text{/}$ | $-p_{\mathrm{f}}\boldsymbol{\mathsf{I}}+2\mu_{\mathrm{f}}\boldsymbol{\mathsf{\varepsilon}}(\mathbf{u}_{\mathrm{f}})$
Fluid density | $\rho_{\mathrm{f}}$ | $\mathrm{kg}\text{/}{\mathrm{m}}^{3}$ | $1060$
Dynamic viscosity | $\mu_{\mathrm{f}}$ | $\mathrm{Pa}\text{\,}\mathrm{s}\text{/}$ | $0.004$
Table 4: Physical parameters for the Navier–Stokes equations with their respective ranges and units. $\boldsymbol{\mathsf{\varepsilon}}(\mathbf{u}_{\mathrm{f}})$ denotes the symmetric gradient defined as $\frac{1}{2}\left(\nabla_{\mathbf{x}}\mathbf{u}_{\mathrm{f}}+\nabla_{\mathbf{x}}\mathbf{u}_{\mathrm{f}}^{\top}\right)$ | EM Comparison
---|---
| $\mathrm{EDV}_{\mathrm{cl},\mathrm{sim}}$ | $\mathrm{ESV}_{\mathrm{cl},\mathrm{sim}}$ | $\mathrm{SV}_{\mathrm{cl,sim}}$ | $\mathrm{EF}_{\mathrm{cl,sim}}$ | $\mathrm{CO}_{\mathrm{cl,sim}}$ | $\mathrm{P}_{\mathrm{cl,sim}}^{\mathrm{sys}}$
| $\mathrm{ml}\text{/}$ | $\mathrm{ml}\text{/}$ | $\mathrm{ml}\text{/}$ | $\mathrm{\char 37\relax}\text{/}$ | $\mathrm{ml}\text{/}\mathrm{s}$ | $\mathrm{mmHg}\text{/}$
_28-Pre_ | $\nicefrac{{88.16}}{{87.47}}$ | $\nicefrac{{30.62}}{{31.02}}$ | $\nicefrac{{57.54}}{{57.14}}$ | $\nicefrac{{65.27}}{{64.81}}$ | $\nicefrac{{87.46}}{{86.85}}$ | $\nicefrac{{146.037}}{{139.362}}$
_44-Pre_ | $\nicefrac{{91.68}}{{91.67}}$ | $\nicefrac{{31.59}}{{30.95}}$ | $\nicefrac{{60.10}}{{60.72}}$ | $\nicefrac{{65.54}}{{66.24}}$ | $\nicefrac{{76.31}}{{76.32}}$ | $\nicefrac{{158.413}}{{135.236}}$
rel. error $[$\mathrm{\char 37\relax}\text{/}$]$ | $\nicefrac{{0.78}}{{0.01}}$ | $\nicefrac{{1.3}}{{2.0}}$ | $\nicefrac{{0.69}}{{1.03}}$ | $\nicefrac{{0.70}}{{1.07}}$ | $\nicefrac{{0.69}}{{0.013}}$ | $\nicefrac{{4.57}}{{14.63}}$
Table 5: Comparison of clinical indicators and indicators computed from simulation for the EM models. | Electromechanics Model | CFD Model
---|---|---
| NE | NV |  $h$ [$\mathrm{\SIUnitSymbolMicro m}\text{/}$] | DOF | NE | NV |  $h$ [$\mathrm{\SIUnitSymbolMicro m}\text{/}$] | DOFU | DOFP
_28-Pre_ | 747266 | 167509 | 897 | 502527 | 1943060 | 352006 | 746.5 | 1056018 | 352006
_28-Post_ | 632635 | 149174 | 954 | 447522 | 7405128 | 1294264 | 531.6 | 3882792 | 1294264
_44-Pre_ | 727194 | 168804 | 997 | 506412 | 2285005 | 412728 | 717 | 1238184 | 412728
Table 6: Discretization details for the studied cases. Shown are the number of elemens (NE), number of vertices (NV), average edge length $h$ in $\mathrm{\SIUnitSymbolMicro m}\text{/}$, degrees of freedom for displacement (DOF), degrees of freedom for velocity (DOFU), degrees of freedom for pressure (DOFP). | Outlet
---|---
| DCA | BCA | RSC | LSC
$R\leavevmode\nobreak\ [$\mathrm{kPa}\text{\,}\mathrm{ms}\text{/}\mathrm{ml}$]$ | $590.46$ | $264.6$ | $1058.24$ | $1058.24$
$Z\leavevmode\nobreak\ [$\mathrm{kPa}\text{\,}\mathrm{ms}\text{/}\mathrm{ml}$]$ | $29.52$ | $13.23$ | $52.91$ | $52.91$
$C\leavevmode\nobreak\ [$\mathrm{ml}\text{/}\mathrm{kPa}$]$ | $1.69$ | $3.78$ | $0.944$ | $0.944$
Table 7: Windkessel parameters for case _28-Pre_. | Outlet
---|---
| DCA | BCA | RSC | LSC
$R\leavevmode\nobreak\ [$\mathrm{kPa}\text{\,}\mathrm{ms}\text{/}\mathrm{ml}$]$ | $2480.07$ | $276.01$ | $466.23$ | $7440.2$
$Z\leavevmode\nobreak\ [$\mathrm{kPa}\text{\,}\mathrm{ms}\text{/}\mathrm{ml}$]$ | $124.003$ | $13.9$ | $23.31$ | $372.01$
$C\leavevmode\nobreak\ [$\mathrm{ml}\text{/}\mathrm{kPa}$]$ | $0.403$ | $3.62$ | $2.14$ | $0.134$
Table 8: Windkessel parameters for case _44-Pre_. Flux Comparison
---
|  | _28-Pre_ | _44-Pre_
| Unit | $\varepsilon_{\mathrm{DCA}}$ | $\varepsilon_{\mathrm{ASC}}$ | $\varepsilon_{\mathrm{ASC}}$ | $\varepsilon_{\mathrm{BCA}}$ | $\varepsilon_{\mathrm{LSC}}$
$Q_{\mathrm{peak,sim}}$ | $\mathrm{ml}\text{/}\mathrm{s}$ | $85.5073$ | $286.056$ | $316.713$ | $160.493$ | $132.540$
$Q_{\mathrm{peak,cl}}$ | $\mathrm{ml}\text{/}\mathrm{s}$ | $70.3071$ | $290.719$ | $352.114$ | $171.571$ | $109.290$
rel. error | $\mathrm{\char 37\relax}\text{/}$ | $21.62$ | $1.604$ | $10.054$ | $6.46$ | $21.27$
Table 9: Comparison of clincal estimated flow rates and simulated flow rates
through the various planes for cases _28-Pre_ and _44-Pre_.

## Figures

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 [10]
(D). A labeled FE geometry $\Omega_{\mathrm{s,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 Table
2. Boundaries $\Gamma_{\mathrm{s,D}}^{0}$ used for prescribing homogeneous
Dirichlet boundary conditions are sketched as blue curves. From this mesh the
EM submesh $\Omega_{\mathrm{s}}^{0}$ (E) and the unsmoothed blood pool (F)
were extracted. Boundary $\Gamma_{\mathrm{s,N}}^{0}$ was used to prescribe
pressure boundary conditions inside the LV and $\Gamma_{\mathrm{s,bp}}^{0}$ is
the surface of the blood pool.  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_{\mathrm{av}}$) to a 3-element Windkessel model
representing aortic impedance, $Z_{\mathrm{c}}$, and peripheral arterial
compliance, $C$, and resistance, $R$, during ejection, and through the mitral
valve ($R_{\mathrm{mv}}$) to a constant pressure $p_{\mathrm{la}}$ in the left
atrium during filling. Negative flows $-q_{\mathrm{la}}$ and
$-q_{\mathrm{lv}}$ mean the respective cavity is ejecting, while positive flow
means cavity is being filled.  Figure 3: (A) Invasive clinical recordings from
cases _28-Pre_ and _44-Pre_. Top: Recorded aortic pressure
$\mathrm{P}_{\mathrm{ao}}$ (black curve) and recorded LV pressure
$\mathrm{P}_{\mathrm{LV}}$ (blue curve). Marked with dashed lines are Systolic
pressure $\mathrm{P}_{\mathrm{sys}}$, mean arterial pressure $\mathrm{MAP}$,
and diastolic pressure $\mathrm{P}_{\mathrm{dia}}$; Center: Volume change in
the LV, $\mathrm{V}_{\mathrm{LV}}$, in red ranging from end-diastolic volume
$\mathrm{EDV}$ to end-systolic volume $\mathrm{ESV}$. Bottom: LV flow
$\mathrm{Q}_{\mathrm{LV}}$ in orange with marked peak flow
$\mathrm{Q}_{\mathrm{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 $\mathrm{V}_{0}$. Figure 4: Processing
workflow used for generating blood pool FE models: (A) and (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) and
(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_{\mathrm{f,mov}}^{t}$ colored in orange, outlet boundaries
$\Gamma_{\mathrm{f,outflow},i}^{t}$ colored in blue with attached illustration
of the 3-element Windkessel models.  Figure 5: CoA anatomy of case 28 before
and after virtual stenting procedure. CoA location is indicated with a red
circle. Figure 6: Panel (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). Panel (II) shows the processing stages of
kinematic transfer for the _28-Pre_ case at maximum displacement. (A)
Displacement $\mathbf{d}_{\mathrm{s}}$ on EM mesh $\Omega_{\mathrm{s}}^{0}$.
(B) Displacement $\mathbf{d}_{\mathrm{s}}$ extended to conformal EM blood pool
mesh $\Omega_{\mathrm{s,total}}^{0}$ which serves as hanging background mesh
for the kinematic transfer onto the CFD blood pool mesh
$\Omega_{\mathrm{f}}^{0}$. (C) Displacement $\mathbf{d}_{\mathrm{s}}$ on
$\Omega_{\mathrm{s}}^{0}$ superimposed with fluid mesh displacement
$\mathbf{d}_{\mathrm{f}}$ on $\Omega_{\mathrm{f}}^{0}$.  Figure 7: Results of
strong scaling benchmark based on case _28-Post_ with $5.2$ million overall
degrees of freedom. $T_{\mathrm{AvgSolv}}$ is the total solving time divided
by the total amount of linear iterations per simulation run. Figure 8: CFD
results. Subfigures (A), (E) show the given clinical measurements for flow
through different planes. The planes are depicted in Subfigures (B) and (F).
Subfigures (B) and (F) also depict the pressure along the centerlines at peak
flow conditions at $t=$167\text{\,}\mathrm{ms}\text{/}$$ and
$t=$142\text{\,}\mathrm{ms}\text{/}$$ respectively. Subfigure (C) shows
velocity streamlines at peak flow. Subfigure (D) shows the relative pressure
map from the Pressure–Poisson mapping used for validating the pressure drop in
our simulations. Subfigures (G), (H) show velocity streamlines at peak flow
and $t=$200\text{\,}\mathrm{ms}\text{/}$$ for case _44-Pre_. 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.