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

TL;DR

This paper develops a comprehensive electro-mechano-fluidic model of the left ventricle and aorta to better understand the impact of aortic coarctations on ventricular load, integrating detailed physics and clinical data.

Contribution

It introduces a novel first-principles-based EMF model coupling electromechanics with CFD, enabling mechanistic analysis of LV load changes due to aortic coarctations.

Findings

Model demonstrates numerical stability with large deformations.

Workflow cost is comparable to existing image-based models.

Simulation of pre- and post-treatment conditions shows potential for clinical application.

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.

Tables

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

	Sex	Age	EDV	ESV	SV	EF	HR	CO	$P_{\mathrm{ao}/\mathrm{cuff}}^{\mathrm{dia}}$	$P_{\mathrm{ao}/\mathrm{cuff}}^{\mathrm{sys}}$	MAP	$P_{\mathrm{open}}$
			$\mathrm{ml}\text{/}$	$\mathrm{ml}\text{/}$	$\mathrm{ml}\text{/}$	$\%\text{/}$	$\mathrm{bpm}$	$\mathrm{ml}\text{/}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	$71.1/62$	$122.7/138$	$88.3/87.3$	$71.33$
44-Pre	M	12	91.7	31.6	60.09	65.5	76	76.31	$74.6/120$	$125.2/154$	$91.5/131.3$	$74.78$

Table 2. Table 2 : Labels used for defining the subdomains of Ω s , total 0 superscript subscript Ω 𝑠 total 0 \Omega_{s,\mathrm{total}}^{0} .

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 3. Table 3 : Fitted parameters for EM Model.

	EM Fitting
	$S_{\mathrm{peak}}$	$t_{\mathrm{dur}}$	${\tau }_{\mathrm{c0}}$	${\tau }_r$	$t_{\mathrm{emd}}$	$C_{\mathrm{Guc}}$	$R$	$Z$	$C$
	$\mathrm{kPa}\text{/}$	$\mathrm{ms}\text{/}$	$-$	$-$	$\mathrm{ms}\text{/}$	$\mathrm{kPa}\text{/}$	$\mathrm{kPa}\mathrm{ms}\text{/}\mathrm{ml}$	$\mathrm{kPa}\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 4. Table 4 : Physical parameters for the Navier–Stokes equations with their respective ranges and units. 𝜺 ​ ( 𝐮 f ) 𝜺 subscript 𝐮 f \boldsymbol{\mathsf{\varepsilon}}(\mathbf{u}_{\mathrm{f}}) denotes the symmetric gradient defined as 1 2 ​ ( ∇ 𝐱 𝐮 f + ∇ 𝐱 𝐮 f ⊤ ) 1 2 subscript ∇ 𝐱 subscript 𝐮 f subscript ∇ 𝐱 superscript subscript 𝐮 f top \frac{1}{2}\left(\nabla_{\mathbf{x}}\mathbf{u}_{\mathrm{f}}+\nabla_{\mathbf{x}}\mathbf{u}_{\mathrm{f}}^{\top}\right)

Name	Variable	Units	Expression / Value
Velocity	${𝐮}_f$	$m\text{/}s$	$-$
Pressure	$p_f$	$\mathrm{Pa}\text{/}$	$-$
Fluid stress tensor	${𝝈}_f({𝐮}_f,p_f)$	$\mathrm{Pa}\text{/}$	$-p_f𝗜+2{\mu }_f𝜺\left({𝐮}_f\right)$
Fluid density	${\rho }_f$	$\mathrm{kg}\text{/}{m}^3$	$1060$
Dynamic viscosity	${\mu }_f$	$\mathrm{Pa}s\text{/}$	$0.004$

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

	EM Comparison
	${\mathrm{EDV}}_{\mathrm{cl},\mathrm{sim}}$	${\mathrm{ESV}}_{\mathrm{cl},\mathrm{sim}}$	${\mathrm{SV}}_{\mathrm{cl},\mathrm{sim}}$	${\mathrm{EF}}_{\mathrm{cl},\mathrm{sim}}$	${\mathrm{CO}}_{\mathrm{cl},\mathrm{sim}}$	$P_{\mathrm{cl},\mathrm{sim}}^{\mathrm{sys}}$
	$\mathrm{ml}\text{/}$	$\mathrm{ml}\text{/}$	$\mathrm{ml}\text{/}$	$\%\text{/}$	$\mathrm{ml}\text{/}s$	$\mathrm{mmHg}\text{/}$
28-Pre	$88.16/87.47$	$30.62/31.02$	$57.54/57.14$	$65.27/64.81$	$87.46/86.85$	$146.037/139.362$
44-Pre	$91.68/91.67$	$31.59/30.95$	$60.10/60.72$	$65.54/66.24$	$76.31/76.32$	$158.413/135.236$
rel. error $\left[\%\text{/}\right]$	$0.78/0.01$	$1.3/2.0$	$0.69/1.03$	$0.70/1.07$	$0.69/0.013$	$4.57/14.63$

Table 6. Table 6 : Discretization details for the studied cases. Shown are the number of elemens (NE), number of vertices (NV), average edge length h ℎ h in µ ​ m / divide micrometer absent \mathrm{\SIUnitSymbolMicro m}\text{/} , degrees of freedom for displacement (DOF), degrees of freedom for velocity (DOFU), degrees of freedom for pressure (DOFP).

	Electromechanics Model				CFD Model
	NE	NV	$h$ [$\mu m\text{/}$]	DOF	NE	NV	$h$ [$\mu 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 7. Table 7 : Windkessel parameters for case 28-Pre .

	Outlet
	DCA	BCA	RSC	LSC
$R\left[\mathrm{kPa}\mathrm{ms}\text{/}\mathrm{ml}\right]$	$590.46$	$264.6$	$1058.24$	$1058.24$
$Z\left[\mathrm{kPa}\mathrm{ms}\text{/}\mathrm{ml}\right]$	$29.52$	$13.23$	$52.91$	$52.91$
$C\left[\mathrm{ml}\text{/}\mathrm{kPa}\right]$	$1.69$	$3.78$	$0.944$	$0.944$

Table 8. Table 8 : Windkessel parameters for case 44-Pre .

	Outlet
	DCA	BCA	RSC	LSC
$R\left[\mathrm{kPa}\mathrm{ms}\text{/}\mathrm{ml}\right]$	$2480.07$	$276.01$	$466.23$	$7440.2$
$Z\left[\mathrm{kPa}\mathrm{ms}\text{/}\mathrm{ml}\right]$	$124.003$	$13.9$	$23.31$	$372.01$
$C\left[\mathrm{ml}\text{/}\mathrm{kPa}\right]$	$0.403$	$3.62$	$2.14$	$0.134$

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

Flux Comparison
		28-Pre		44-Pre
	Unit	${\epsilon }_{\mathrm{DCA}}$	${\epsilon }_{\mathrm{ASC}}$	${\epsilon }_{\mathrm{ASC}}$	${\epsilon }_{\mathrm{BCA}}$	${\epsilon }_{\mathrm{LSC}}$
$Q_{\mathrm{peak},\mathrm{sim}}$	$\mathrm{ml}\text{/}s$	$85.5073$	$286.056$	$316.713$	$160.493$	$132.540$
$Q_{\mathrm{peak},\mathrm{cl}}$	$\mathrm{ml}\text{/}s$	$70.3071$	$290.719$	$352.114$	$171.571$	$109.290$
rel. error	$\%\text{/}$	$21.62$	$1.604$	$10.054$	$6.46$	$21.27$

Equations

$${\Omega }_{\mathrm{s},\mathrm{total}}^0=\bigcup _{i\in I}{\Omega }_{\mathrm{s},i}^0,$$

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

$${\Omega }_{\mathrm{s}}^0$$

$${\Omega }_{\mathrm{s},\mathrm{bp}}^0={\tilde{\Omega }}_f^0$$

$${\Gamma }_{\mathrm{s},\mathrm{N}}^0$$

$${\Gamma }_{\mathrm{s},\mathrm{H}}^0$$

$${\Gamma }_{\mathrm{s},\mathrm{bp}}^0$$

$$\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.$$

$$\beta C_{\mathrm{m}}\frac{\partial V_{\mathrm{m}}}{\partial t}={\nabla }_{\mathbf{X}}\cdot {\sigma }_{\mathrm{i}}{\nabla }_{\mathbf{X}}{V}_{\mathrm{m}}+I_{\mathrm{foot}}-\beta I_{\mathrm{ion}},$$

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

$${\mathbf{d}}_{\mathrm{s}}=0\text{on}{\Gamma }_{\mathrm{s},\mathrm{D}}^0,$$

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

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

$$\mathbf{S}={\mathbf{S}}_{\mathrm{pas}}+{\mathbf{S}}_{\mathrm{act}},$$

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

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

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

$$\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\}.$$

$$S_{\mathrm{a}}(t,\lambda )=S_{\mathrm{peak}}\varphi (\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),\text{for }0<t_{\mathrm{s}}<t_{\mathrm{dur}},$$

$$\varphi =\tanh (\mathrm{ld}(\lambda -{\lambda }_0)),{\tau }_{\mathrm{c}}={\tau }_{{\mathrm{c}}_0}+{\mathrm{ld}}_{\mathrm{up}}(1-\varphi ),t_{\mathrm{s}}=t-t_{\mathrm{a}}-t_{\mathrm{emd}}$$

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

$$\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}},$$

$$C_i\frac{\mathrm{d}{p}_{\mathrm{d},i}}{\mathrm{d}t}+\frac{p_{\mathrm{d},i}}{R_i}$$

$$p_{\mathrm{wk},i}$$

$${\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 {\sigma }_{\mathrm{f}}({\mathbf{u}}_{\mathrm{f}},p_{\mathrm{f}})$$

$${\nabla }_{\mathbf{x}}\cdot {\mathbf{u}}_{\mathrm{f}}$$

$${\mathbf{u}}_{\mathrm{f}}$$

$${\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\mathbf{u}_{\mathrm{f}}\big{|}_{t=0}$$

$${\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|{\mathbf{u}}_{\mathrm{f}}\cdot {\mathbf{n}}_{\mathrm{f}}\right|\right)$$

$${\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\}.$$

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

$${\Gamma }_{\mathrm{f},\mathrm{mov}}^t:=\partial {\Omega }_{\text{f}}^t\backslash \bigcup _{i=1}^{n_{\text{outlets}}}{\Gamma }_{\mathrm{f},\mathrm{outflow},i}^t,$$

$$\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}})$$

Full text

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 [math]; 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 [math]; 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

[TABLE]

with the index set

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

homogeneous Neumann boundary conditions

[TABLE]

and inhomogeneous Neumann boundary conditions

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

with

[TABLE]

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

[TABLE]

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

[TABLE]

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$

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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],

[TABLE]

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}})$

[TABLE]

with the classical bilinear form of the Navier–Stokes equations

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where stress and strain tensor are

[TABLE]

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

[TABLE]

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 [math] 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](#S2.SS5.SSS5), are summarized in Tables [7](#Sx4.T7) and [8](#Sx4.T8). 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](#Sx5.F8)) 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.

Tables

Figures