Assessment of wall stresses and mechanical heart power in the left ventricle: Finite element modeling versus Laplace analysis

TL;DR

This study compares finite element modeling and Laplace analysis for assessing wall stresses and mechanical power in the left ventricle of patients with aortic stenosis, highlighting the superior accuracy of FE models.

Contribution

It demonstrates that patient-specific finite element models provide more accurate assessments of LV stresses and power than Laplace-based methods in clinical settings.

Findings

Laplace estimates approximate global mean stress but miss spatial heterogeneity.

FE models accurately capture stress and power distribution in the LV wall.

Laplace methods are less accurate for assessing mechanical power.

Abstract

Introduction: Stenotic aortic valve disease (AS) causes pressure overload of the left ventricle (LV) that may trigger adverse remodeling and precipitate progression towards heart failure (HF). As myocardial energetics can be impaired during AS, LV wall stresses and biomechanical power provide a complementary view of LV performance that may aide in better assessing the state of disease. Objectives: Using a high-resolution electro-mechanical (EM) in silico model of the LV as a reference, we evaluated clinically feasible Laplace-based methods for assessing global LV wall stresses and biomechanical power. Methods: We used N = 4 in silico finite element (FE) EM models of LV and aorta of patients suffering from AS. All models were personalized with clinical data under pre-treatment conditions. LV wall stresses and biomechanical power were computed accurately from FE kinematic data and compared to Laplace-based estimation methods which were applied to the same FE model data. Results and Conclusion: Laplace estimates of LV wall stress are able to provide a rough approximation of global mean stress in the circumferential-longitudinal plane of the LV. However, according to FE results spatial heterogeneity of stresses in the LV wall is significant, leading to major discrepancies between local stresses and global mean stress. Assessment of mechanical power with Laplace methods is feasible, but these are inferior in accuracy compared to FE models. The accurate assessment of stress and power density distribution in the LV wall is only feasible based on patient-specific FE modeling.

Figures

Tables

Table 1. Table 1: Pre-treatment AS patient characteristics from MRI and non-invasive cuff pressure recordings including end-diastolic volume (EDV), end-systolic volume (ESV), stroke volume (SV), ejection fraction (EF), heart rate (HR), diastolic and systolic cuff pressures ( p dia subscript 𝑝 dia p_{\rm{dia}} and p sys subscript 𝑝 sys p_{\rm{sys}} ), mean arterial pressure (MAP), wall thickness at the LV equator measured in septum / lateral free wall ( h ℎ h ), pressure drop across aortic valve ( Δ ​ p Δ 𝑝 \Delta p ), presence of hypertension (HT) and mitral valve regurgitation (MVR).

	Sex	Age	EDV	ESV	SV	EF	HR	$p_{\mathrm{dia}}$	$p_{\mathrm{sys}}$	MAP	$h$	$\Delta p$	HT	MVR
		[years]	[ml]	[ml]	[ml]	[%]	[min${}_{}{}^{-1}]$	[mmHg]	[mmHg]	[mmHg]	[mm]	[mmHg]
A	F	$63$	$112.0$	$46.0$	$66.00$	$58.93$	$53$	$74$	$126$	$91.33$	$12.0/12.5$	$95$	No	No
B	M	$73$	$121.0$	$54.7$	$66.32$	$54.81$	$81$	$75$	$134$	$94.67$	$11.2/13.8$	$62$	No	Mild
C	M	$54$	$118.2$	$42.2$	$76.14$	$64.42$	$75$	$71$	$141$	$94.33$	$16.0/18.2$	$79$	Yes	Mild
D	M	$85$	$172.0$	$103.0$	$69.00$	$40.12$	$68$	$79$	$144$	$100.67$	$14.0/15.2$	$59$	Yes	No

Table 2. Table 2: Geometric parameters inner radius, r 𝑟 r , outer radius, R 𝑅 R , wall thickness, h ℎ h and wall thickness to radius of curvature ratio, h / r ℎ 𝑟 h/r of spherical shell models 𝐒𝐩𝐡 5 subscript 𝐒𝐩𝐡 5 \mathrm{\mathbf{Sph}_{5}} , 𝐒𝐩𝐡 25 subscript 𝐒𝐩𝐡 25 \mathrm{\mathbf{Sph}_{25}} and 𝐒𝐩𝐡 150 subscript 𝐒𝐩𝐡 150 \mathrm{\mathbf{Sph}_{150}} and of image-based anatomical LV models in the stress-free reference configuration.

	$r$ [mm]	$R$ [mm]	$h$ [mm]	$h/r$
${\mathrm{𝐒𝐩𝐡}}_5/{\mathrm{𝐒𝐩𝐡}}_{25}/{\mathrm{𝐒𝐩𝐡}}_{150}$	$15/15/15$	$15.5/17.5/30.0$	$0.5/2.5/15.0$	$0.033/0.166/1.0$
${\mathrm{𝐋𝐕}}_A$/ ${\mathrm{𝐋𝐕}}_B$/ ${\mathrm{𝐋𝐕}}_C$/ ${\mathrm{𝐋𝐕}}_D$	$16.9/19.5/17.1/22.1$	$30.6/33.6/36.9/38.5$	$13.7/14.1/19.8/16.4$	$0.82/0.72/1.16/0.74$

Table 3. Table 3: Fitted parameters of circulatory, active stress and passive mechanical model components.

	afterload				active stress		passive mechanical model
	$Z_v\left[\frac{\text{kPa.ml}}{\text{ms}}\right]$	$Z_a\left[\frac{\text{kPa.ml}}{\text{ms}}\right]$	$R\left[\frac{\text{kPa.ml}}{\text{ms}}\right]$	$C\left[\frac{\text{ml}}{\text{kPa}}\right]$	$S_{\mathrm{peak}}$ [kPa]	${\tau }_c$ [ms]	$a$ [kPa]
${\mathrm{𝐋𝐕}}_A$	$35.82$	$26.00$	$187.74$	$15.23$	$69$	$80$	$0.5$
${\mathrm{𝐋𝐕}}_B$	$16.08$	$11.03$	$72.65$	$26.62$	$85$	$35$	$0.65$
${\mathrm{𝐋𝐕}}_C$	$15.93$	$12.78$	$77.34$	$25.50$	$63$	$40$	$0.5$
${\mathrm{𝐋𝐕}}_D$	$22.34$	$11.09$	$62.73$	$30.93$	$98$	$58$	$0.5$

Table 4. Table 4: Comparison of FE-based mean wall stresses σ ¯ r ​ r subscript ¯ 𝜎 𝑟 𝑟 {{\overline{\sigma}}_{rr}} , σ ¯ θ ​ θ subscript ¯ 𝜎 𝜃 𝜃 {{\overline{\sigma}}_{\theta\theta}} and σ ¯ φ ​ φ subscript ¯ 𝜎 𝜑 𝜑 {{\overline{\sigma}}_{\varphi\varphi}} in radial, azimuthal and meridional direction, respectively, with the Laplace-based wall stress estimates σ L,h subscript 𝜎 L,h {{\sigma}_{\text{L,h}}} , σ L,H subscript 𝜎 L,H {{\sigma}_{\text{L,H}}} and σ L,V subscript 𝜎 L,V {{\sigma}_{\text{L,V}}} . All stresses refer to the maximum applied pressure of p = 4 𝑝 4 p=4~{} kPa.

Setup	${\overline{\sigma }}_{rr}$ [kPa]	${\overline{\sigma }}_{\theta \theta }$ [kPa]	${\overline{\sigma }}_{\phi \phi }$ [kPa]	${\sigma }_{\text{L,h}}$ [kPa]	${\sigma }_{\text{L,H}}$ [kPa]	${\sigma }_{\text{L,V}}$ [kPa]	#elements	$\overline{dx}$ [mm]
${\mathrm{𝐒𝐩𝐡}}_5$	$-1.77$	$99.46$	$99.47$	$100.47$	$99.48$		$83825$	$0.65$
${\mathrm{𝐒𝐩𝐡}}_{25}$	$-1.16$	$15.63$	$15.63$	$16.93$	$15.99$		$40974$	$1.23$
${\mathrm{𝐒𝐩𝐡}}_{150}$	$-0.63$	$1.44$	$1.44$	$2.35$	$1.65$		$213742$	$1.70$
${\mathrm{𝐋𝐕}}_{A-\mathrm{Gu}}$	$-0.90$	$7.05$	$3.33$	$5.24$	$4.40$	$6.83$	$420704$	$1.52$
${\mathrm{𝐋𝐕}}_{B-\mathrm{Gu}}$	$-0.85$	$8.71$	$4.29$	$5.94$	$5.08$	$8.19$	$332221$	$1.74$
${\mathrm{𝐋𝐕}}_{C-\mathrm{Gu}}$	$-0.67$	$4.47$	$2.12$	$3.54$	$2.76$	$4.44$	$456553$	$1.84$
${\mathrm{𝐋𝐕}}_{D-\mathrm{Gu}}$	$-0.83$	$6.97$	$3.93$	$5.73$	$4.88$	$7.07$	$394808$	$1.86$

Table 5. Table 5: Comparison of FE-based peak mechanical power P ^ int subscript ^ 𝑃 int {{\hat{P}}_{\text{int}}} with different estimates that were all evaluated at the instant of peak pressure, t p ^ subscript 𝑡 ^ 𝑝 t_{\hat{p}} : P int ​ ( t p ^ ) subscript 𝑃 int subscript 𝑡 ^ 𝑝 {{P}_{\text{int}}}(t_{\hat{p}}) , P ext ​ ( t p ^ ) subscript 𝑃 ext subscript 𝑡 ^ 𝑝 {{P}_{\text{ext}}}(t_{\hat{p}}) , P int,h ​ ( t p ^ ) subscript 𝑃 int,h subscript 𝑡 ^ 𝑝 {{P}_{\text{int,h}}}(t_{\hat{p}}) , P int,H ​ ( t p ^ ) subscript 𝑃 int,H subscript 𝑡 ^ 𝑝 {{P}_{\text{int,H}}}(t_{\hat{p}}) , P int,V ​ ( t p ^ ) subscript 𝑃 int,V subscript 𝑡 ^ 𝑝 {{P}_{\text{int,V}}}(t_{\hat{p}}) ; peak hydrodynamic power P ^ ext subscript ^ 𝑃 ext {{\hat{P}}_{\text{ext}}} ; product of peak pressure and peak flow in the LV p ^ LV ⋅ q ^ LV ⋅ subscript ^ 𝑝 LV subscript ^ 𝑞 LV \hat{p}_{\mathrm{LV}}\cdot\hat{q}_{\mathrm{LV}} ; relative internal heart power marker IHP and external heart power EHP ; cardiac power efficiency P eff subscript 𝑃 eff {{P}_{\text{eff}}} and its IHP -based approximation P eff,clin subscript 𝑃 eff,clin {{P}_{\text{eff,clin}}} .

	${\hat{P}}_{\text{int}}$	$P_{\text{int}}\left(t_{\hat{p}}\right)$	${\hat{P}}_{\text{ext}}$	$P_{\text{ext}}\left(t_{\hat{p}}\right)$	${\hat{p}}_{\mathrm{LV}}\cdot {\hat{q}}_{\mathrm{LV}}$	$P_{\text{int,h}}\left(t_{\hat{p}}\right)$	$P_{\text{int,H}}\left(t_{\hat{p}}\right)$	$P_{\text{int,V}}\left(t_{\hat{p}}\right)$	IHP	EHP	$P_{\text{eff}}$	$P_{\text{eff,clin}}$
	[W]	[W]	[W]	[W]	[W]	[W]	[W]	[W]	[W]	[W]	[1]	[1]
${\mathrm{𝐋𝐕}}_A$	$-5.13$	$-5.02$	$-5.50$	$-5.42$	$-5.54$	$-3.93$	$-2.62$	$-7.94$	$-6.50$	$-2.56$	$0.49$	$0.39$
${\mathrm{𝐋𝐕}}_B$	$-13.64$	$-13.52$	$-13.56$	$-13.48$	$-13.59$	$-9.24$	$-7.30$	$-12.94$	$-32.76$	$-6.00$	$0.43$	$0.18$
${\mathrm{𝐋𝐕}}_C$	$-17.94$	$-17.75$	$-17.44$	$-17.39$	$-17.48$	$-14.62$	$-8.56$	$-27.09$	$-22.27$	$-6.30$	$0.35$	$0.28$
${\mathrm{𝐋𝐕}}_D$	$-22.01$	$-21.92$	$-21.80$	$-21.74$	$-21.84$	$-20.32$	$-14.82$	$-30.20$	$-29.59$	$-8.02$	$0.36$	$0.27$

Table 6. Table 6: Comparison of FE-based mean wall stresses σ ¯ r ​ r subscript ¯ 𝜎 𝑟 𝑟 {{\overline{\sigma}}_{rr}} , σ ¯ θ ​ θ subscript ¯ 𝜎 𝜃 𝜃 {{\overline{\sigma}}_{\theta\theta}} and σ ¯ φ ​ φ subscript ¯ 𝜎 𝜑 𝜑 {{\overline{\sigma}}_{\varphi\varphi}} in radial, azimuthal and meridional direction, respectively, with the Laplace-based wall stress estimates σ L,h subscript 𝜎 L,h {{\sigma}_{\text{L,h}}} , σ L,H subscript 𝜎 L,H {{\sigma}_{\text{L,H}}} and σ L,V subscript 𝜎 L,V {{\sigma}_{\text{L,V}}} . All stresses refer to the maximum applied pressure of p = 4 𝑝 4 p=4~{} kPa.

Setup	${\overline{\sigma }}_{rr}$ [kPa]	${\overline{\sigma }}_{\theta \theta }$ [kPa]	${\overline{\sigma }}_{\phi \phi }$ [kPa]	${\sigma }_{\text{L,h}}$ [kPa]	${\sigma }_{\text{L,H}}$ [kPa]	${\sigma }_{\text{L,V}}$ [kPa]	#elements	$\overline{dx}$ [mm]
${\mathrm{𝐋𝐕}}_{A-\mathrm{Dem}}$	$-0.39$	$6.56$	$3.77$	$5.16$	$4.32$	$6.81$	$420704$	$1.52$
${\mathrm{𝐋𝐕}}_{B-\mathrm{Dem}}$	$-0.38$	$8.26$	$4.78$	$5.92$	$5.07$	$8.24$	$332221$	$1.74$
${\mathrm{𝐋𝐕}}_{C-\mathrm{Dem}}$	$-0.21$	$4.07$	$2.26$	$3.43$	$2.65$	$4.41$	$456553$	$1.84$
${\mathrm{𝐋𝐕}}_{D-\mathrm{Dem}}$	$-0.31$	$6.53$	$4.28$	$5.81$	$4.69$	$7.06$	$394808$	$1.86$

Table 7. Table 7: Comparison of FE-based biomechanical and hemodynamic work, W int subscript 𝑊 int {{W}_{\text{int}}} and W ext subscript 𝑊 ext {{W}_{\text{ext}}} , with the Laplace-based work estimates W int,h subscript 𝑊 int,h {{W}_{\text{int,h}}} , W int,H subscript 𝑊 int,H {{W}_{\text{int,H}}} and W int,V subscript 𝑊 int,V {{W}_{\text{int,V}}} for passive inflation with a pressure of p = 4 𝑝 4 p=4 kPa.

Setup	$W_{\text{ext}}$ [mJ]	$W_{\text{int}}$ [mJ]	$W_{\text{int,h}}$ [mJ]	$W_{\text{int,H}}$ [mJ]	$W_{\text{int,V}}$ [mJ]	#elements	$\overline{dx}$ [mm]
${\mathrm{𝐒𝐩𝐡}}_5$	$16.82$	$16.92$	$19.86$	$19.65$		$83825$	$0.65$
${\mathrm{𝐒𝐩𝐡}}_{25}$	$10.94$	$10.99$	$11.58$	$10.87$		$40974$	$1.23$
${\mathrm{𝐒𝐩𝐡}}_{150}$	$8.30$	$8.15$	$10.89$	$7.51$		$54449$	$2.31$
${\mathrm{𝐋𝐕}}_{A-\mathrm{Dem}}$	$67.03$	$65.48$	$67.94$	$55.67$	$89.30$	$420704$	$1.52$
${\mathrm{𝐋𝐕}}_{B-\mathrm{Dem}}$	$111.67$	$110.01$	$124.41$	$104.31$	$183.99$	$332221$	$1.74$
${\mathrm{𝐋𝐕}}_{C-\mathrm{Dem}}$	$96.19$	$94.90$	$100.48$	$75.44$	$128.44$	$456553$	$1.84$
${\mathrm{𝐋𝐕}}_{D-\mathrm{Dem}}$	$129.19$	$126.98$	$140.30$	$117.23$	$171.12$	$394808$	$1.86$
${\mathrm{𝐋𝐕}}_{A-\mathrm{Gu}}$	$93.19$	$89.34$	$89.11$	$72.65$	$116.26$	$420704$	$1.52$
${\mathrm{𝐋𝐕}}_{B-\mathrm{Gu}}$	$149.55$	$145.27$	$167.23$	$131.08$	$219.04$	$332221$	$1.74$
${\mathrm{𝐋𝐕}}_{C-\mathrm{Gu}}$	$136.84$	$133.42$	$146.17$	$109.42$	$183.75$	$456553$	$1.84$
${\mathrm{𝐋𝐕}}_{D-\mathrm{Gu}}$	$182.77$	$177.51$	$173.92$	$143.91$	$216.12$	$394808$	$1.86$

Equations

$$-\nabla \cdot \mathbf{\sigma }(\mathbf{u},t)$$

$$\mathbf{\sigma }(\mathbf{u},t)\cdot \mathbf{n}$$

$$\mathbf{u}$$

$${\mathbf{\sigma }}_{\text{pas}}$$

$${\mathbf{\sigma }}_{\text{act}}$$

$${\Psi }_{\mathrm{Gu}}(\mathbf{C})=\frac{\kappa }{2}{\left(\log J\right)}^2+\frac{a}{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}}(\mathbf{C})=\frac{\kappa}{2}\left(\operatorname{log}\,J\right)^{2}+\frac{a}{2\,b}\left\{\operatorname{exp}\left[b\,\big{(}\operatorname{tr}(\overline{\mathbf{C}})-3\big{)}\right]-1\right\},$$

$$S_{\mathrm{a}}=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_s<t_{\mathrm{dur}},$$

$$\boldsymbol{\sigma}=\left(\begin{array}[]{ccc}\sigma_{xx}&\sigma_{xy}&\sigma_{xz}\\ \sigma_{yx}&\sigma_{yy}&\sigma_{yz}\\ \sigma_{zx}&\sigma_{zy}&\sigma_{zz}\end{array}\right)=\mathbf{P}\left(\begin{array}[]{ccc}\sigma_{rr}&\sigma_{r\varphi}&\sigma_{r\theta}\\ \sigma_{\varphi r}&\sigma_{\varphi\varphi}&\sigma_{\varphi\theta}\\ \sigma_{\theta r}&\sigma_{\theta\varphi}&\sigma_{\theta\theta}\\ \end{array}\right)\mathbf{P}^{\top}=\mathbf{P}\,{{\boldsymbol{\sigma}}_{\text{sph}}}\,\mathbf{P}^{\top}$$

$${\overline{\sigma }}_{\star \star }(\mathbf{u},t)=\frac{1}{\left|\Omega \right|}\underset{\Omega }{\int }{\sigma }_{\star \star }(\mathbf{u},t)\mathrm{d}\mathbf{x}$$

$${{\boldsymbol{\sigma}}_{\text{sph}}}\approx\left(\begin{array}[]{ccc}0&0&0\\ 0&{{\sigma}_{\text{circ}}}&0\\ 0&0&{{\sigma}_{\text{circ}}}\\ \end{array}\right).$$

$${\sigma }_{\text{circ}}=\frac{pr}{2h\left(1+\frac{h}{2r}\right)}={\sigma }_{\text{L,H}}$$

$${\sigma }_{\text{circ}}=\frac{pr}{2h}={\sigma }_{\text{L,h}}$$

$${\sigma }_{\text{circ}}=\frac{p}{{\left(\frac{V_{\text{cav}}+V_{\text{myo}}}{V_{\text{cav}}}\right)}^{2/3}-1}={\sigma }_{\text{L,V}},$$

$$p_{\mathrm{int}}(\mathbf{x},t)=\mathbf{\sigma }(\mathbf{u},t):\dot{\mathbf{\epsilon }}(\mathbf{u},t),$$

$$P_{\text{int}}(t)=\underset{\Omega }{\int }\mathbf{\sigma }(\mathbf{u},t):\dot{\mathbf{\epsilon }}(\mathbf{u},t)\mathrm{d}\mathbf{x}$$

$$W_{\text{int}}=\underset{t_0}{\overset{t}{\int }}{P}_{\text{int}}(\tau )d\tau .$$

$$P_{\text{int,⋆}}(t)=V_{\text{myo}}(t){\sigma }_{\text{L,⋆}}(t)\underset{\approx {\dot{\mathbf{\epsilon }}}_{\mathrm{circ}}}{\underbrace{\left(\frac{\dot{r}(t)}{r_0}+\frac{\dot{R}(t)}{R_0}\right)}}$$

$$W_{\text{int,⋆}}(t)=\underset{t_0}{\overset{t}{\int }}{P}_{\text{int,⋆}}(\tau )\mathrm{d}\tau .$$

$$\text{IHW}=V_{\text{myo}}{\sigma }_{\text{L,⋆}}.$$

$$\text{IHP}=\frac{\text{IHW}}{T_{\text{sct}}}.$$

$$P_{\text{ext}}=pq=p\frac{\mathrm{d}{V}_{\text{cav}}}{\mathrm{d}t}$$

$$W_{\text{ext}}=\underset{t_0}{\overset{t}{\int }}p(\tau )q(\tau )\mathrm{d}\tau ,$$

$$W_{\text{ext}}=\underset{{V_{\text{cav}}}_0}{\overset{{V_{\text{cav}}}_1}{\int }}p(V_{\text{cav}})\mathrm{d}{V}_{\text{cav}}.$$

$$\text{EHP}={\bar{P}}_{\text{ext,ao}}=\frac{1}{T_{\mathrm{sys}}}\underset{t_{\mathrm{ED}}}{\overset{t_{\mathrm{ES}}}{\int }}{p}_{\mathrm{ao}}q\mathrm{d}t\approx \text{MAP}\cdot \text{CO},$$

$$P_{\text{eff,clin}}=\frac{\text{EHP}}{\text{IHP}}.$$

$$P_{\text{eff}}=\frac{{\bar{P}}_{\text{ext,ao}}}{{\hat{P}}_{\text{int}}}\approx \frac{{\bar{P}}_{\text{ext,ao}}}{{\hat{P}}_{\text{ext}}}.$$

$${\bar{r}}_i(t)=\frac{1}{N}\sum _{j=1}^N{r}_{i,j}(t)\text{and}{\bar{h}}_i(t)=\frac{1}{N}\sum _{j=1}^N{h}_{i,j}(t).$$

$$\bar{r}(t)=\frac{1}{M}\sum _{i=1}^M{\bar{r}}_i(t)\text{and}\bar{h}(t)=\frac{1}{M}\sum _{i=1}^M{\bar{h}}_i(t).$$

$$p{r}^2\pi \stackrel{!}{=}(R^2-r^2){\sigma }_{\text{circ}}\pi .$$

Full text

Keywords: Aortic stenosis, transvalvular pressure gradient, heart failure

00footnotetext: Abbreviations: AVD, aortic valve disease; LV, left ventricle; HF, heart failure; PV, pressure-volume; EM, electro-mechanical; IHP, internal mechanical heart power; FE, finite element; IVC, isovolumetric contraction; IVR, isovolumetric relaxation.

1 Introduction

In AS elevated pressure gradients impose a higher load upon the LV. Under such conditions, the pressure produced by the LV must increase in order to achieve an adequate cardiac output that meets the metabolic demands. This requires the LV wall to generate higher active forces, which can be achieved either by an increase in wall stresses or a change in ventricular shape and mass. Such pressure overload conditions, if persistent for long enough, trigger adverse remodeling processes, eventually precipitating progression towards HF [17]. Treatments aim at alleviating pressure overload by reducing transvalvular pressure gradients closer to normal levels by surgical or catheter based aortic valve replacement [3]. However, re-stenosis frequently occurs and despite a successful reduction of transvalvular pressure gradients, a majority of patients remains hypertensive, consequently showing increased risk for irreversible course of HF and higher morbidity and mortality [14]. Thus, a successful reduction of pathologically elevated pressure gradients alone cannot be considered a reliable prognostic marker of long-term post-treatment outcomes in these patient cohorts.

As a consequence, alternative biomarkers beyond pressure gradients are sought to that provide a complementary view of cardiac function and, potentially, offer a higher predictive power with regard to outcomes. In a recent study, the use of end-diastolic or end-systolic wall stresses as assessed by a wall stress index has been proposed as a novel diagnostic criterion of HF [4]. This is physiologically motivated as elevated wall stress levels are assumed to impair the balance between metabolic supply and demand [36] by hindering perfusion and, thus, contribute towards adverse remodeling [1]. Wall stresses are directly linked to the mechanical power generated by the myocardial muscle and the work performed by it and as such can be considered a metabolic marker. Different approaches have been proposed to assess work and the energy expenditure of the myocardium. As a direct measurement of energy metabolism Positron Emission Tomography (PET) was used [19, 20], however the method is limited due to its complexity including the need for tracers involving ionizing radiation. More recently, a concept of biomechanical internal myocardial heart power (IHP), necessary to maintain adequate cardiac output (EHP, external heart power), has been introduced in patients with aortic coarctation [15]. Findings in this cohort suggest that the ratio $\text{EHP}/\text{IHP}$, referred to as power efficiency, improved mostly in those cases with elevated IHP. While potential marker qualities of such concepts need to be further evaluated, it has remained a yearned-for goal to lower energy expenditure and increase efficiency of the myocardium likewise in any treatment procedure, including those for stenotic valvular and vascular disease.

Another method used to determine the work performed by the muscle is pressure-volume (PV) relations. These are usually measured using conductance catheter techniques. However, these procedures are invasive, time-consuming and expensive. Alternatively, PV-loops are measured with 3D echo or MRI [2], but even these methods are complex and pressure and volume traces are not recorded simultaneously.

Despite the diagnostic potential of markers based on wall stresses and expended mechanical power, this assessment has not evolved towards a routinely used diagnostic tool in the clinic mainly due to methodological limitations. Attempts to address this relied upon different variants of Laplace’s law, which require the acquisition of only a small number of measures representing LV cavity volume, wall width and pressure [4, 15]. However, these approaches are based on simplifications with regard to LV geometry, tissue structure and biomechanical properties as they assume the LV a thin-walled mechanically isotropic spherical shell. The accuracy and validity of these simplifications have not been firmly established, thus casting doubt on the reliability and fidelity of any metrics based on them [27, 41]. Experimental validation based on direct measurements of stresses in vivo is challenging and not feasible yet with currently available technologies. However, an indirect inference is viable using computational tools such as FE modeling where 3D wall stresses can be computed from a set of reliable strains – either measured in vivo [5, 23] or computed in silico [7] – using constitutive material models which are derived from ex vivo measurements [18] and material parameters fitted to clinical data [7, 40].

To evaluate the accuracy of Laplace analysis for estimating global wall stresses and mechanical power in the LV we employed four FE based EM LV models which have been previously fitted and validated with clinical data [7] under pre-treatment conditions. These models provide reliable strain data at a high spatio-temporal resolution from which wall stresses, biomechanical power and work in the LV can be determined at the best possible accuracy. Laplace analysis was applied to these in silico models to estimate hoop stress and mechanical power over a cardiac cycle and compared to the global ground truth data based on FE analysis.

2 Methods

2.1 Patient Data

Data from four AS patients with clinical indication for aortic valve treatment, all preceding a cardiac magnetic resonance study, were used (Tab. 1). AS treatment indicators included valve area and/or systolic pressure drop across the valve. 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.

2.2 Biomechanical FE model

The ventricular myocardium was modeled as a non-linear, hyperelastic, nearly incompressible and anisotropic material with a layered organization of myocytes and fibres that is characterized by a right-handed orthonormal set of basis vectors [18, 22]. These basis vectors consist of the fiber axis $\mathbf{f}_{0}(\mathbf{x})$, which coincides with the prevailing orientation of the myocytes at location $\mathbf{x}$, the sheet axis $\mathbf{s}_{0}(\mathbf{x})$, and the sheet-normal axis $\mathbf{n}_{0}(\mathbf{x})$. The mechanical deformation of the tissue is described by Cauchy’s equation of motion under stationary equilibrium assumptions leading to a quasi-static boundary value problem: for a given pressure $p(t)$, find the unknown displacement $\mathbf{u}$ such that

[TABLE]

holds for $t\in[0,T]$. By $\Omega$ we denote the deformed geometry and by $\Gamma=\partial\Omega$ we define its boundary with $\Gamma=\overline{\Gamma_{D}}\cup\overline{\Gamma_{H}}\cup\overline{\Gamma_{N}}$ and $\left|{\Gamma_{D}}\right|>0$. The normal outward vector of $\Gamma$ is denoted by $\mathbf{n}$. The total Cauchy stress tensor $\boldsymbol{\sigma}$ refers to the sum of a passive stress tensor ${{\boldsymbol{\sigma}}_{\text{pas}}}$ and an active stress tensor ${{\boldsymbol{\sigma}}_{\text{act}}}$. That is, $\boldsymbol{\sigma}={{\boldsymbol{\sigma}}_{\text{pas}}}+{{\boldsymbol{\sigma}}_{\text{act}}}$ with

[TABLE]

where $\mathbf{F}$ is the deformation gradient, $\Psi$ is the strain energy function, $\mathbf{f}_{0}$ is fiber orientation in the reference configuration, $J=\operatorname{det}\mathbf{F}$ is the Jacobian, $\mathbf{C}=\mathbf{F}^{\top}\mathbf{F}$ is the right Cauchy–Green strain tensor and $S_{a}$ is the scalar active contractile stress generated by the myocytes acting along $\mathbf{f}_{0}$.

The passive behavior of myocardial tissue was modeled using two material models, either the transversely-isotropic Guccione et al. model [18],

[TABLE]

where

[TABLE]

and $\overline{\mathbf{E}}=\frac{1}{2}(J^{-\frac{2}{3}}\mathbf{C}-\mathbf{I})$ is the modified isochoric Green–Lagrange strain tensor, or the isotropic Demiray model [12]

[TABLE]

with $\overline{\mathbf{C}}=J^{-\frac{2}{3}}\mathbf{C}$ the modified isochoric right Cauchy–Green tensor. In both models, Eqs. (4) and (6), the bulk modulus $\kappa$, which serves as a penalty parameter to enforce near incompressibility, was chosen as $\kappa=650$ kPa.

A simplified phenomenological contractile model [31] was used to represent active stress generation:

[TABLE]

where $S_{\mathrm{peak}}$ is the peak isometric tension, $\phi(\lambda)$ is a non-linear function dependent on fiber stretch $\lambda=\left|\mathbf{F}\mathbf{f}_{0}\right|$ describing the length-dependence of active stress generation, $t_{\mathrm{s}}$ is the onset of contraction, $\tau_{\mathrm{c}}$ is the upstroke time constant, $t_{\mathrm{dur}}$ is the active stress transient duration and $\tau_{\mathrm{r}}$ is the downstroke time constant. This simplified model allows efficient fitting to patient data as the parameters for peak stress, $S_{\rm peak}$, and time constant of contraction, $\tau_{\mathrm{c}}$, and twitch duration, $t_{\mathrm{dur}}$, are related to the two clinical key parameters of interest, peak pressure and maximum rate of pressure increase, in an intuitive manner.

Solving these equations under given mechanical boundary conditions using the FE method at a sufficiently high spatio-temporal discretization provides an accurate description of tissue kinematics. Computed displacement $\mathbf{u}$ serve as input then in a post-processing procedure to evaluate wall stresses $\boldsymbol{\sigma}(\mathbf{x},t)$ and to compute internal power expended by the LV (see Sec. 2.6).

A Newton scheme was applied in each time step to linearize the nonlinear boundary value problem (1) yielding a non-symmetric FE system. The linear FE system was solved by a parallel GMRES algorithm with an algebraic multigrid preconditioner. For the Newton scheme a relative tolerance of $1.0e^{-5}$ and an absolute tolerance of $1.0e^{-8}$ was used as stopping criterion.

2.3 Verification of finite element model

To verify the FE-based calculation of stress-derived metrics, a geometrically simple and well-studied benchmark problem was chosen for which circumferential hoop stresses can be found from Laplace’s law under the following assumptions:

• (A1)

The wall material is isotropic.

• (A2)

The shape is a symmetric spherical shell with inner radius, $r$, and outer radius, $R$.

• (A3)

The thickness of the wall, $h=R-r$, is sufficiently small, that is, the wall thickness to radius of curvature ratio is small, $h/r\ll 1$.

Since all these assumptions are violated in the LV which is orthotropic (A1), non-spherical in shape (A2), and thick-walled (A3) with $h/r\approx 1$, differences between Laplace analysis and FE computation are to be expected. Three configurations were considered, an ideal thin-walled spherical shell, $\mathrm{\mathbf{Sph}_{5}}$, which complies with all assumptions (A1)–(A3) and thus can serve as a reference for FE validation, and two thicker-walled spheres, $\mathrm{\mathbf{Sph}_{25}}$ and $\mathrm{\mathbf{Sph}_{150}}$, where assumption (A3) is increasingly violated. Geometries and mechanical boundary conditions are illustrated in Fig. LABEL:sub@fig:sph:geom–LABEL:sub@fig:sph:bcs.

The inner radius was chosen as $r=15.0$ mm in all models with $h$ varying from $h=0.5$ mm to $h=15.0$ mm in $\mathrm{\mathbf{Sph}_{5}}$ and $\mathrm{\mathbf{Sph}_{150}}$, respectively. The choices for $\mathrm{\mathbf{Sph}_{150}}$ are representative of the $h/r$ ratios found in the LVs of patients in this study (Tab. 2). In line with assumption (A1), the nonlinear isotropic material law stated in Eq. (6) was employed with $a=10$ kPa and $b=8$. Passive inflation experiments were performed by solving (1) with ${{\boldsymbol{\sigma}}_{\text{act}}}\equiv 0$ and applying a pressure $p$ in the range from [math] to $4$ kPa to the endocardial surface, $\Gamma_{\mathrm{endo}}$, which covers the range of pressures observed in vivo during diastole. Pressure at the epicardial surface, $\Gamma_{\mathrm{epi}}$, was assumed to be zero. To render the solution of this pure Neumann problem unique, displacement boundary conditions were enforced at the intersections of the Cartesian axes with the epicardial surface by restricting displacements to the respective intersecting axes (see Fig LABEL:sub@fig:sph:bcs). Unstructured tetrahedral FE meshes were generated for the $\mathrm{\mathbf{Sph}_{5}}$, $\mathrm{\mathbf{Sph}_{25}}$ and $\mathrm{\mathbf{Sph}_{150}}$ geometries where the mean spatial resolution, $\bar{\;\mathrm{d}x}$, was increased until solutions were deemed converged.

2.4 LV model

2.4.1 Anatomical Modeling

FE meshes of the LV anatomy and aortic root were generated from 3D whole heart MRI acquired at end diastole (ED) with $1.458\times 1.548\times 2$ mm resolution at the German Heart Center Berlin. Multi-label segmentation of the LV myocardium, LV cavity and aortic lumen was done using the ZIB Amira software (https://amira.zib.de/). Segmentations were smoothed and upsampled to a $0.1$ mm isotropic resolution [11]. The wall of the aorta was automatically generated by dilation of the aortic lumen with a thickness of $1.2$ mm, and the aorta was clipped before the branch of the brachiocephalic artery. Due to limited resolution valves were not segmented, but were included in the FE model as a thin layer of tissue for applying pressure boundary conditions and computation of cavity volume. The multi-label segmentations were meshed using CGAL (http://www.cgal.org/) with a target discretization of $1.25$ mm in the LV myocardium and $1$ mm in the aortic wall. For the transversely-isotropic Guccione material model, see Eq. (4), we equipped all models with a rule-based fiber architecture [9], where fibers rotated linearly from $-75^{\circ}$ at the epicardium to $+75^{\circ}$ at the endocardium (Fig. 2.A). All anatomical models built are shown in Fig. 2.

2.4.2 Model fitting

To remove rigid body motion and provide physiological boundary conditions that allow a vertical movement of the LV base, as observed in-vivo, mechanical boundary conditions were applied by fixing the terminal rim of the clipped aorta (Fig. 2.A) and resting the apex of the LV on an elastic cushion which was rigidly anchored at its base. Constitutive relations were represented by Eq. (4). Using the ED geometry, default material parameters and an estimated ED pressure (EDP), an initial guess of the stress-free reference configuration was computed by unloading the model using a backward displacement method [35, 10]. Since clinically recorded data of the ED PV-relation (EDPVR) are often limited, the Klotz relation [24] providing an empiric description of EDPVR, $p({{V}_{\text{cav}}})$, was used as target to steer the fitting of constitutive parameters. In absence of accurate measurements of EDP we refrained from fitting all material parameters to $p({{V}_{\text{cav}}})$. Rather, default values for the parameters $b_{f}=18.48$, $b_{t}=3.58$ and $b_{fs}=1.627$ were used as reported in the literature [18], and only the scaling parameter $a$ was adjusted individually for each patient. With a given data point (EDV,EDP) $a$ was fitted to minimize the difference in stress-free residual volume, $V_{\rm{0,dia}}$, between model and Klotz curve. This yielded values for $a$ of $0.5$, $0.65$, $0.5$ and $0.5$ for the cases $\mathrm{\mathbf{LV}_{A}}$, $\mathrm{\mathbf{LV}_{B}}$, $\mathrm{\mathbf{LV}_{C}}$ and $\mathrm{\mathbf{LV}_{D}}$, respectively.

A three-element Windkessel model of LV afterload was used to provide the pressure-flow relationship during ejection [39] (see Fig. 3). LV models were parameterized to match clinically recorded PV-data using LV cavity volume traces, ${{V}_{\text{cav}}}(t)$, determined from Cine-MRI with a temporal resolution of $45.28$, $29.63$, $32.00$ and $35.29$ ms for $\mathrm{\mathbf{LV}_{A}}$, $\mathrm{\mathbf{LV}_{B}}$, $\mathrm{\mathbf{LV}_{C}}$ and $\mathrm{\mathbf{LV}_{D}}$, respectively. Continuously monitored invasive pressure recordings were not available as catheterization was not indicated. Peak pressure in the LV was determined by estimating peak pressure in the aortic root from cuff pressure measurements and by determining the pressure drop at peak flow across the aortic valve from ultrasound flow measurements using Bernoulli’s law [13]. Windkessel parameters representing the aortic input impedance, $Z$, comprising the flow resistance of aortic valve, $Z_{v}$, and the characteristic input impedance of the aorta, $Z_{\rm c}$, as well as resistance $R$ and compliance $C$ of the arterial system were fit to reproduce estimated LV peak pressure using measured volume traces ${{V}_{\text{cav}}}(t)$ as input.

In a final step, active mechanical properties were fit to the same hemodynamic data used for fitting the afterload model. A reaction-eikonal model was used to generate activation sequences and simulate action potential propagation in the LV [29]. Active stress generation was triggered with a prescribed electromechanical delay when the upstroke of the action potential crossed the $-40$ mV threshold. Parameters peak stress, $S_{\rm peak}$, time constant of contraction, $\tau_{\mathrm{c}}$, and twitch duration, $t_{\mathrm{dur}}$, were adjusted manually to fit peak pressure, $\hat{p}_{\mathrm{LV}}$, duration of pressure pulse and flow. Due to the intuitive link of the active stress model given in Eq. (7) with the fitting targets, a satisfactory fit was achieved within $\leq$ 5 simulation runs. The goodness of fit was deemed sufficiently accurate when the clinically measured metrics EF, SV, MAP and peak LV pressure, $\hat{p}_{\mathrm{LV}}$, were matched within a margin of error of $\pm 5\%$. Clinical input data and fitted model parameters are summarized in Tab. 3.

2.5 Myocardial wall stresses

Stresses can be computed from deformations $\mathbf{u}$ using constitutive material models based on ex vivo experimental data which link stresses to strains. In the FE model stress tensors $\boldsymbol{\sigma}(\mathbf{x},t)$ are computed by evaluating Eqs. (2) and (3), which yields a $3\times 3$ tensor where only 6 components are independent for symmetry reasons. For the models $\mathrm{\mathbf{Sph}_{5}}$, $\mathrm{\mathbf{Sph}_{25}}$ and $\mathrm{\mathbf{Sph}_{150}}$ the stress tensor simplifies. Due to the assumption of isotropy in (A1) and symmetry in (A2), any solution, if expressed in a spherical coordinate system, must also be symmetric. Quantities computed in the FE Cartesian coordinate system are recast in spherical coordinates as defined in Fig. LABEL:sub@fig:sph:csys using a projection matrix, $\mathbf{P}$. For the total Cauchy stress tensor $\boldsymbol{\sigma}$ we obtain

[TABLE]

with the projection matrix $\mathbf{P}=(\mathbf{e}_{r},\mathbf{e}_{\varphi},\mathbf{e}_{\theta})$ and $r\in\mathbb{R}^{+},\theta\in[0,\pi],\varphi\in[0,2\,\pi)$.

While all quantities in the spherical models must be perfectly symmetric this is not necessarily the case in the FE solutions. Depending on spatial resolution and boundary conditions, a minor numerical jitter around mean values will inevitably occur. For comparing FE with Laplace analysis averaged mean quantities were therefore computed over the entire domain by

[TABLE]

with $\mathbf{u}$ being the FE solution at time $t$ and $\star\in\{r,\theta,\varphi\}$.

In a thin-walled spherical shell stresses in azimuthal and meridional direction must be equal due to symmetry, that is, ${{\sigma}_{\text{circ}}}=\sigma_{\varphi\varphi}=\sigma_{\theta\theta}$, and with $h\ll r$ radial stresses can be assumed to be negligible relative to circumferential hoop stresses, that is, $0\approx\sigma_{rr}\ll{{\sigma}_{\text{circ}}}$. The stress tensor in spherical coordinates simplifies therefore to

[TABLE]

As a reference for verifying FE-based stresses different variants of Laplace’s law were used. In particular, we use an extension of Laplace’s law which takes into account the finite thickness of the wall

[TABLE]

and refer to stress estimates based on this formula as Laplace stress in thick-walled spherical shells, ${{\sigma}_{\text{L,H}}}$. Exploiting assumption (A3), that is, $\frac{h}{r}\ll 1$ and thus $\left(1+\frac{h}{2\,r}\right)\approx 1$, allows a further simplification yielding

[TABLE]

which we refer to as Laplace stress in thin-walled spherical shells, ${{\sigma}_{\text{L,h}}}$. Finally, we consider a volume-based estimation of ${{\sigma}_{\text{circ}}}$ [26] defined as

[TABLE]

which has been used previously in clinical studies [4], We note that Eq. (13) is equivalent to Eq. (11) for a spherical shell. However, when applied to a non-spherical structure such as the LV this is not the case. Eq. (13) may offer advantages as the determination of ${{V}_{\text{cav}}}$ and ${{V}_{\text{myo}}}$ may be less ambiguous than the determination of a representative inner radius $r$ and wall thickness $h$ (see Sec. 2.8.1).

2.6 Myocardial power and work

For a given displacement $\mathbf{u}$ at time $t\in[0,T]$, where $T$ refers to the duration of a cardiac cycle and $t_{0}=t_{\mathrm{ED}}=0$ marks the end of diastole, the biomechanical power density, $p_{\mathrm{int}}$, generated or consumed at location $\mathbf{x}$ within the LV wall can be computed by evaluating

[TABLE]

where $\dot{\boldsymbol{\varepsilon}}$ is the strain rate tensor and $\mathbf{A}:\mathbf{B}=\operatorname{tr}(\mathbf{A}^{\top}\,\mathbf{B})$ denotes the double contraction of two tensors, see e.g. [21] for further details. Integration of Eq. (14) over the entire myocardial wall yields the global biomechanical power, ${{P}_{\text{int}}}$,

[TABLE]

and integration of Eq. (15) over time yields an expression of biomechanical work, ${{W}_{\text{int}}}$, performed

[TABLE]

Based on Laplace’s law biomechanical power can be estimated using

[TABLE]

where $\star$ denotes which formula was used for estimating the circumferential wall stress that is, $\star\in\left\{\text{h},\text{H},\text{V}\right\}$ and $\dot{\boldsymbol{\varepsilon}}_{\mathrm{circ}}$ approximates circumferential strains. For a derivation of Eq. (17) see supplementary material. Laplace-based mechanical work is estimated analogously to Eq. (16) yielding

[TABLE]

In addition, a recently introduced also Laplace-based relative power indicator, IHP, was evaluated which attempts to estimate the power generated by the LV around the instant of peak pressure, $t_{\hat{p}}$. Based on [32], the mechanical work expended, IHW, during contraction time, $T_{\text{sct}}$, defined as the time elapsed between the onset of isovolumetric contraction at $t_{\mathrm{ED}}$ and the instant of peak stress in the LV at $t_{\hat{p}}$, is approximated by

[TABLE]

IHW is interpreted as a measure of the mechanical potential energy stored in the LV from which a measure of the peak biomechanical power generated by the LV between $t_{\mathrm{ED}}$ and $t_{\hat{p}}$ is derived then by

[TABLE]

Note that Eq. (20), in contrast to Eq. (17), does not include any measure of $\dot{\boldsymbol{\varepsilon}}$. Thus, while consistent in terms of physical units, IHP must be considered a relative indicator and not a physical measure of power.

2.7 Hydrodynamic power and work

Hydrodynamic power, ${{P}_{\text{ext}}}$, is given by

[TABLE]

where $p$ is the hydrostatic pressure acting at endocardial surface, $\Gamma_{\mathrm{endo}}$ and $q$ represents blood flow out of the LV cavity during ejection. Hydrodynamic work, ${{W}_{\text{ext}}}$, is then the work expended by the LV myocardium when changing the volume of its cavity, ${{V}_{\text{cav}}}$, given by

[TABLE]

or, equivalently, expressed as PV work by

[TABLE]

In absence of active stresses, i.e. ${{\boldsymbol{\sigma}}_{\text{act}}}=0$, and isovolumetric contstraints imposed by valves upon ${{V}_{\text{cav}}}$, ${{P}_{\text{int}}}\equiv{{P}_{\text{ext}}}$ must hold. Under such conditions, external work can therefore serve as a reference for validating the FE-based computation of internal power and work. This is not necessarily the case during the isovolumetric phases of a heartbeat where internal work may be expended which does not necessarily manifest as external work. During these phases, changes in $\boldsymbol{\sigma}$ occur which may entail shape changes of the LV myocardium and thus induce a non-zero strain rate tensor $\dot{\boldsymbol{\varepsilon}}$. However, due to the isovlumetric constraints imposed by the incompressibility of the blood pool and the closed state of all valves, no global change in cavity volume can occur, i.e. $\mathrm{d}{{V}_{\text{cav}}}=0$.

Under healthy conditions, hydrodynamic power in the LV cavity equals the power delivered to the arterial system as transvalvular power losses are small. However, in AS cases where transvalvular pressure gradients, $\Delta p$, are significant, the effective hydrodynamic power externally delivered to the arterial system is reduced. Following [15], we define external hydrodynamic heart power, EHP, as

[TABLE]

where MAP and CO are mean arterial pressure and cardiac output, respectively, and $p_{\mathrm{ao}}=p_{\mathrm{LV}}-\Delta p$ is the pressure in the aorta ascendens. Power efficiency, ${{P}_{\text{eff}}}$, has been defined previously in [15] as the ratio

[TABLE]

Since ${{P}_{\text{eff}}}$ essentially relates the mean hydrodynamic power delivered to the arterial system to the peak biomechanical power generated by the LV myocardium during systole, ${{P}_{\text{eff}}}$ can be expressed as

[TABLE]

where ${{\hat{P}}_{\text{int}}}$, ${{\hat{P}}_{\text{ext}}}$ and ${{\bar{P}}_{\text{ext,ao}}}$ are determined based on Eqs. (15), (21) and (24), respectively. That is, ${{P}_{\text{eff}}}$ can be estimated from hemodynamic data using ${{\hat{P}}_{\text{ext}}}$ or from LV deformation analysis using ${{\hat{P}}_{\text{int}}}$. Unlike ${{P}_{\text{eff,clin}}}$ in Eq. (25), which compares an absolute measure of external hydrodynamic power to a relative indicator of internal biomechanical peak power, Eq. (26) provides a physically consistent comparison.

2.8 Evaluation of Laplace-based assessment of wall stress and power

Human EM LV models which were validated against clinical data (see Sec. 2.4.2), provided accurate ground truth data on strains $\boldsymbol{\varepsilon}(\mathbf{u},t)$ and stresses $\boldsymbol{\sigma}(\mathbf{u},t)$ in the LV wall. Using these as reference, Laplace analysis was applied to the in silico models to assess its accuracy and validity.

2.8.1 Determination of clinical input data for Laplace-based analysis

Geometric input parameters $r$ and $h$ required for Laplace analysis must be determined from clinical imaging datasets. As LV shape deviates markedly from a spherical shell representative mean parameters of $r$ and $h$ must be found. Since there is no unique best solution to establish a geometric correspondence between LV shape and a spherical shell, various methods have been used in clinical applications. Typically, transverse slices from short-axis Cine-MRI scans were analyzed to measure, either manually or semi-automatically, $r$ and $h$, where $h$ is measured in the postero-lateral wall, the septum or an average is taken. The analysis is either carried out in one representative mid-cavity LV short axis slice, or a number of slices is selected to capture representative basal, mid-cavity and apical LV cross-sections.

Similar issues arise when applying Laplace analysis to in silico datasets. In order to extract $r$ and $h$ as objectively as possible without operator bias, automated processing workflows were implemented (see Fig. 4). Analogous to the $z$-slice selection in MRI protocols, the unstructured FE meshes of the LV models were decomposed into slices of $\approx 8$ mm resolution, comparable to the MRI out-of-plane resolution. Decomposition was achieved by first determining the long axis, $\mathbf{z}$, of the LV using principal component analysis (see Fig. LABEL:sub@fig:lv:slicing), which yielded, depending on the spatial extent of the LV long-axis of a given model, between $10$ and $14$ slices, each slice $i$ is centered around $z_{i}$. A mean $z$ coordinate of the LV in its current configuration, $\bar{z}_{c}(t)$, was computed to define the center slice plane using the long-axis unit vector, $\mathbf{e_{\mathrm{z}}}$, and the center $z_{i}(t)$ of individual slices was shifted, keeping slice width and distance to the LV center $\bar{z}_{c}(t)$ constant. Within a selected plane, radial vectors, $\mathbf{r}_{i,j}(t)$, were computed that emanated from $z_{i}$ and were oriented in polar angles $\varphi_{j}$ ranging from [math] to $360^{\circ}$ with an angular sampling of $\Delta\varphi=9^{\circ}$. For each vector $\mathbf{r}_{i,j}(t)$, the intersection with surfaces, $\Gamma_{\mathrm{endo}}(t)$ and $\Gamma_{\mathrm{epi}}(t)$, was determined, yielding $N=2\pi/\Delta\varphi$ inner radii, $r_{i,j}(t)$, outer radii, $R_{i,j}(t)$ and wall widths, $h_{i,j}(t)=R_{i,j}(t)-r_{i,j}(t)$ (Fig. LABEL:sub@fig:lv:slice_avg). Mean radius, $\bar{r}_{i}(t)$, and wall width, $\bar{h}_{i}(t)$, were determined as the arithmetic average

[TABLE]

Finally, multi-slice mean $\bar{r}$ and $\bar{h}$ were computed by averaging over $M$ slices

[TABLE]

The time course of the mean values $\bar{r}(t)$ and $\bar{h}(t)$ (see Fig. LABEL:sub@fig:lv:mean_r_h_of_t) was plugged then into the respective Laplace equations to compute stress in Eqs. (11), (12), and (13), power in Eq.  (17) and work in Eq. (18).

2.8.2 Simulation protocols and data analysis

To evaluate the influence of violating the assumptions (A1) and (A2), passive inflation experiments were performed with LV models and the anisotropic material given in Eq. (4) following the same protocol as applied before to the spherical shell models $\mathrm{\mathbf{Sph}_{5}}$, $\mathrm{\mathbf{Sph}_{25}}$ and $\mathrm{\mathbf{Sph}_{150}}$ in Sec. 2.3. Laplace-based stress estimates ${{\sigma}_{\text{L,h}}}$, ${{\sigma}_{\text{L,H}}}$ and ${{\sigma}_{\text{L,V}}}$ were compared to the mean stresses obtained from the FE solution. Stresses were evaluated with respect to an ellipsoidal coordinate system to facilitate a comparison with stresses in the spherical shell models (see Fig. LABEL:sub@fig:sph:csys). An ellipsoidal coordinate system was constructed for the LV models by assigning fiber and sheet orientations using a rule-based method with a constant fiber angle of $0^{\circ}$ [9]. Stress components ${{\sigma}_{rr}}(\mathbf{x})$, ${{\sigma}_{\varphi\varphi}}(\mathbf{x})$ and ${{\sigma}_{\theta\theta}}(\mathbf{x})$ were averaged according to Eq. (9) yielding ${{\overline{\sigma}}_{rr}}$, ${{\overline{\sigma}}_{\varphi\varphi}}$ and ${{\overline{\sigma}}_{\theta\theta}}$, respectively. Laplace-based estimation of power, {{P}_{\text{int,\star}}}, was compared to FE-based power, ${{P}_{\text{int}}}$, and to external hydrodynamic power in the LV cavity, ${{P}_{\text{ext}}}$.

Laplace analysis was applied to clinically fitted EM LV models $\mathrm{\mathbf{LV}_{A}}$–$\mathrm{\mathbf{LV}_{D}}$ to compare LV stress {{\sigma}_{\text{L,\star}}}, power {{P}_{\text{int,\star}}} and IHP over an entire systolic cycle to the FE-based stresses ${{\overline{\sigma}}_{rr}}$, ${{\overline{\sigma}}_{\varphi\varphi}}$, ${{\overline{\sigma}}_{\theta\theta}}$ and ${{\overline{\sigma}}_{\text{mean}}}$, and biomechanical power ${{P}_{\text{int}}}$. Further, biomechanical power due to LV deformation, ${{P}_{\text{int}}}$, and hydrodynamic power, ${{P}_{\text{ext}}}$, derived from PV data were also compared to assess differences during isovolumetric phases.

2.9 Numerical Solution

Discretization of all PDEs and the solution of the arising systems of equations relied upon the Cardiac Arrhythmia Research Package framework [37]. Details on FE discretization [33] as well as numerical solution of electrophysiology [38, 28, 29] and electro-mechanics [8] equations have been described in detail previously. Both electrophysiolgy and mechanics FE solvers were validated previously in N-version benchmark studies [30, 25].

3 Results

3.1 Verification of FE model

The FE implementation was verified by performing passive inflation experiments with spherical shell models for which Laplace laws are known to be almost exact ($\mathrm{\mathbf{Sph}_{5}}$) or, at least, sufficiently accurate ($\mathrm{\mathbf{Sph}_{25}}$ and $\mathrm{\mathbf{Sph}_{150}}$). The resulting EDPVRs and principal components of the Cauchy stress ${{\boldsymbol{\sigma}}_{\text{sph}}}$, evaluated in spherical coordinates and globally averaged to yield mean stresses ${{\overline{\sigma}}_{\text{circ}}}$ and ${{\overline{\sigma}}_{rr}}$, are shown in Figs. LABEL:sub@fig:_edpvr_sph and LABEL:sub@fig:_edpvr_stress. A numerical comparison of stresses and work at the maximum pressure of $p=4$ kPa is provided in Tab. 4 and Tab. 7 in the supplement.

Agreement of FE-based mean circumferential stress ${{\overline{\sigma}}_{\text{circ}}}$ with Laplace laws was very close, that is, ${{\overline{\sigma}}_{\text{circ}}}=\frac{1}{2}({{\overline{\sigma}}_{\theta\theta}}+{{\overline{\sigma}}_{\varphi\varphi}})\approx{{\overline{\sigma}}_{\theta\theta}}\approx{{\overline{\sigma}}_{\varphi\varphi}}\approx{{\sigma}_{\text{L,H}}}\approx{{\sigma}_{\text{L,h}}}$ held. With increasing $h$ ${{\sigma}_{\text{L,H}}}$ provided estimates that were closer to the FE-based stress ${{\overline{\sigma}}_{\text{circ}}}$ than ${{\sigma}_{\text{L,h}}}$ (Tab. 4). The simple Laplace overestimated ${{\overline{\sigma}}_{\text{circ}}}$ in $\mathrm{\mathbf{Sph}_{5}}$, $\mathrm{\mathbf{Sph}_{25}}$ and $\mathrm{\mathbf{Sph}_{150}}$ by $1.02\%$, $8.32\%$ and $63.19\%$ whereas with ${{\sigma}_{\text{L,H}}}$ deviations were much smaller with $0.02\%$, $2.3\%$ and $14.58\%$. Radial stresses were negligible in the thinner-walled models $\mathrm{\mathbf{Sph}_{5}}$ and $\mathrm{\mathbf{Sph}_{25}}$, that is, ${{\overline{\sigma}}_{rr}}\ll{{\overline{\sigma}}_{\text{circ}}}$, but not in the thick-walled model $\mathrm{\mathbf{Sph}_{150}}$ where ${{\overline{\sigma}}_{rr}}$ amounted to $\approx 43.75\%$ of ${{\overline{\sigma}}_{\text{circ}}}$. A comparison of FE-based work ${{W}_{\text{int}}}$ to Laplace-based ${{W}_{\text{int,h}}}$ and ${{W}_{\text{int,H}}}$ is given in Tab. 7 in the supplement.

3.2 Evaluation of Laplace-based assessment of wall stresses and power

After verification with spherical shell models FE analysis was applied to a validated in silico EM LV model to compute stresses, power and work during both diastolic and systolic phases. Since all assumptions underlying Laplace laws are violated in LV models the FE-based results were considered the ground truth and, thus, could be used to gauge the accuracy of Laplace-based assessment of LV mechanics.

3.2.1 Passive inflation of LV models

LV models $\mathrm{\mathbf{LV}_{A}}$–$\mathrm{\mathbf{LV}_{D}}$ were inflated following the same protocol as in Sec. 2.3 (see Fig. 6.A). The temporal evolution of FE- and Laplace-based stresses, power and work are shown in Fig. 6 for model $\mathrm{\mathbf{LV}_{D}}$. Minor quantitative differences to other models $\mathrm{\mathbf{LV}_{A}}$–$\mathrm{\mathbf{LV}_{C}}$ were observed, but qualitatively the overall behavior was identical. Stresses at $p=4$ kPa are summarized in Tab. 4, incurred work is given in Tab. 7 in the supplement.

3.2.2 Analysis of LV cycle experiments

Using Cine-MRI-based LV volume traces and estimated $\hat{p}_{\mathrm{LV}}$ as inputs the models $\mathrm{\mathbf{LV}_{A}}$–$\mathrm{\mathbf{LV}_{D}}$ were fitted over the cycle phases isovolumetric contraction (IVC), ejection and isovolumetric relaxation (IVR) (Fig. 7). All models replicated the clinical metrics of interest such as SV, EF or peak aortic pressure $\hat{p}_{\mathrm{ao}}$ with sufficient accuracy ($<5\%$).

Fig. 8 compares the time course of the averaged FE-based quantities azimuthal, meridionial, radial and circumferential mean stresses, ${{\overline{\sigma}}_{\varphi\varphi}}$, ${{\overline{\sigma}}_{\theta\theta}}$, ${{\overline{\sigma}}_{rr}}$ and ${{\overline{\sigma}}_{\text{circ}}}=\frac{1}{2}({{\overline{\sigma}}_{\varphi\varphi}}+{{\overline{\sigma}}_{\theta\theta}})$, respectively and power ${{P}_{\text{int}}}$ to the Laplace-based estimation of stresses {{\sigma}_{\text{L,\star}}} and power {{P}_{\text{int,\star}}}. In all cases, the Laplace-based stresses ${{\sigma}_{\text{L,h}}}$ and ${{\sigma}_{\text{L,H}}}$ tended to underestimate the FE-based mean circumferential stress ${{\overline{\sigma}}_{\text{circ}}}$, being closer to the azimuthal stress ${{\overline{\sigma}}_{\varphi\varphi}}$, whereas ${{\sigma}_{\text{L,V}}}$ overestimated ${{\overline{\sigma}}_{\text{circ}}}$ and was closer to ${{\overline{\sigma}}_{\theta\theta}}$. Further, both Laplace stresses or globally averaged mean stresses deviate noticeably from the true local stresses acting at a given location (Fig.  9).

Laplace-based power estimates ${{P}_{\text{int,h}}}$, ${{P}_{\text{int,H}}}$ and ${{P}_{\text{int,V}}}$ were qualitatively comparable to the exact FE-based ${{P}_{\text{int}}}$, but quantitatively marked discrepancies were observed. The time course of Laplace-based power showed both a faster onset and decay with an early peak in power. Quantitative differences between the Laplace estimates were also significant with ${{P}_{\text{int,V}}}>{{P}_{\text{int,h}}}>{{P}_{\text{int,H}}}$. Around the instant $t_{\hat{p}}$ deviations were in the range of $-2.41/+2.92$, $-6.22/-0.58$, $-9.19/+9.34$ and $-7.10/+8.28$ W for $\mathrm{\mathbf{LV}_{A}}$, $\mathrm{\mathbf{LV}_{B}}$, $\mathrm{\mathbf{LV}_{C}}$ and $\mathrm{\mathbf{LV}_{D}}$, respectively. The relative IHP marker led to large deviations (see Fig. 13 in the supplement) from the true FE-based mechanical power, even around the instant of $\hat{p}_{\mathrm{LV}}$ the IHP marker was intended for.

A numerical comparison between markers of LV peak power is given in Tab. 5. Differences between true mechanical peak power ${{\hat{P}}_{\text{int}}}$ and ${{P}_{\text{int}}}$ at the instant of peak pressure, $t_{\hat{p}}$, were minor, with the maximum difference being $\left|{{{\hat{P}}_{\text{int}}}-{{P}_{\text{int}}}(t_{\hat{p}})}\right|<0.19$ W or $2.15\%$. Laplace estimation of ${{P}_{\text{int}}}$, evaluated at $t_{\hat{p}}$, misestimated ${{\hat{P}}_{\text{int}}}$ by $5.1\%$ - $54.7\%$. Interestingly, ${{P}_{\text{int,h}}}$ performed better than ${{P}_{\text{int,H}}}$ in all cases, but with non-negligible maximum relative errors of $23.37\%$, $32.23\%$, $18.51\%$ and $7.68\%$ for $\mathrm{\mathbf{LV}_{A}}$, $\mathrm{\mathbf{LV}_{B}}$, $\mathrm{\mathbf{LV}_{C}}$ and $\mathrm{\mathbf{LV}_{D}}$, respectively. IHP overestimated ${{\hat{P}}_{\text{int}}}$ significantly, in the range between $24.10\%$ to $140.28\%$. Since the mechanical power developed during isovolumetric phases was marginal (see Fig. 14 in the supplement), the most accurate estimate of ${{P}_{\text{int}}}$ is obtained by computing ${{P}_{\text{ext}}}$ from standard hemodynamic PV data. In terms of peak mechanical power the difference $\left|{{{\hat{P}}_{\text{int}}}-{{\hat{P}}_{\text{ext}}}}\right|$ was less than $0.5\;W$ or $7.18\%$ in all cases where ${{\hat{P}}_{\text{ext}}}$ can be estimated with high accuracy by taking the product of peak power and flow, $\hat{p}_{\mathrm{LV}}\cdot\hat{q}_{\mathrm{LV}}$.

4 Discussion

Wall stress and mechanical power generated by the LV are considered important biomarkers that promise potential clinical utility for diagnosis and as a predictor of post-treatment LV remodeling after interventions [4, 15]. Moreover, the modelling of stresses and power would allow to gain an improved understanding of mechanisms that contribute to adverse remodelling. Laplace analysis would have the charme that inputs such as $p$, $r$, $h$, ${{V}_{\text{myo}}}$ and ${{V}_{\text{cav}}}$ are accessible within routine clinical procedures. However, Laplace analysis is based on a global force balance calculation and relies upon simplifying assumptions on LV shape, tissue structure and biomechanical behavior. This study attempts to establish validity, accuracy and potential limitations of Laplace analysis of stresses and mechanical power generated by the LV by comparing against a FE model for which these quantities can be determined with high accuracy.

4.1 FE verification

FE computation of stresses and mechanical power was verified by performing passive inflation experiments with geometrically well defined spherical shell models of varying wall width for which Laplace laws hold with sufficient accuracy. FE computed circumferential stresses ${{\sigma}_{\text{circ}}}$ in all models agreed closely with the Laplace stresses {{\sigma}_{\text{L,\star}}} (see Fig. 5 and Tab. 4). As expected, with increasing $h$, deviations became more pronounced and the thick-walled Laplace stresses ${{\sigma}_{\text{L,H}}}$ agreed closer with FE stresses than the standard Laplace stress ${{\sigma}_{\text{L,h}}}$. In terms of work expended, more noticeable discrepancies were observed between ${{W}_{\text{int}}}$ and Laplace-based {{W}_{\text{int,\star}}} (see Tab. 7 in the supplement). However, since the agreement between FE-computed internal work ${{W}_{\text{int}}}$ and external work ${{W}_{\text{ext}}}$ was essentially perfect, as expected on grounds of conservation of energy, we concluded that our FE implementation for evaluating stresses, power and work is correct and that the observed deviations are rather attributable to inherent inaccuracies in the Laplace approximations. In particular, we consider the mean strain rate approximation in Eq. (17) and the omission of radial stresses likely candidate causes.

4.2 Laplace versus FE-based stress and power analysis

The validated high resolution in silico model served as a reference for evaluating the accuracy of the Laplace-based approximation of $\boldsymbol{\sigma}$, ${{P}_{\text{int}}}$ and ${{W}_{\text{int}}}$. While the FE models which were built from, fitted to and validated against clinical data, may deviate from clinical data within the limits of clinical data uncertainty, for assessing Laplace analysis, the FE model represents the ground truth as it provides accurate data on stresses $\boldsymbol{\sigma}(\mathbf{x},t)$ and strains $\boldsymbol{\varepsilon}(\mathbf{x},t)$ which can serve to compute local power and work densities ${{p}_{\text{int}}}(\mathbf{x},t)$ and ${{w}_{\text{int}}}(\mathbf{x},t)$, respectively, as well as global ${{P}_{\text{int}}}$, ${{W}_{\text{int}}}$ and ${{W}_{\text{ext}}}$ with highest possible accuracy. All input parameters needed for Laplace analysis can be derived from the FE model with higher accuracy than what is achievable clinically. In this regard, the application of Laplace analysis to the in silico model can be considered a best case scenario.

4.2.1 Wall stress in the LV

Wall stress $\boldsymbol{\sigma}(\mathbf{x})$ in the LV is a tensorial quantity that varies in space (see Fig. 9). The tensor comprises six independent components whereas Laplace stresses {{\sigma}_{\text{L,\star}}} provide only one scalar stress value representing a global circumferential or hoop stress, ${{\sigma}_{\text{circ}}}$. While ${{\sigma}_{\text{circ}}}$ is equivalent to ${{\sigma}_{\varphi\varphi}}$ and ${{\sigma}_{\theta\theta}}$ in a thin-walled spherical shell such as $\mathrm{\mathbf{Sph}_{5}}$ (see Tab. 4) this is not the case in the LV as there is no direct equivalence to any component of $\boldsymbol{\sigma}$. As shown in a FE modeling study by Zhang et al. [41], the correlation of Laplace stresses to fiber and cross-fiber stresses is poor. Conceptually, the force balance consideration used in the derivation of Laplace laws suggests that Laplace stresses are most likely representative of the mean stresses in the longitudinal-circumferential plane, ${{\overline{\sigma}}_{\text{circ}}}=\frac{1}{2}\left({{\overline{\sigma}}_{\varphi\varphi}}+{{\overline{\sigma}}_{\theta\theta}}\right)$. Indeed, a fair qualitative agreement was observed between ${{\overline{\sigma}}_{\text{circ}}}$ and {{\sigma}_{\text{L,\star}}} during passive LV inflation as illustrated in Fig. 6. During ejection the time course of {{\sigma}_{\text{L,\star}}}(t) followed a similar trend as ${{\overline{\sigma}}_{\text{circ}}}$ although waveforms deviated to different degrees owing to the marked differences in the LV anatomies. However, quantitatively discrepancies were significant during both passive inflation and over a LV cycle as evident in Fig. 6.B and in the stress panels of Fig. 8 with substantial differences in stress magnitudes between the various Laplace laws and the global circumferential mean stress with ${{\sigma}_{\text{L,V}}}>{{\overline{\sigma}}_{\text{circ}}}>{{\sigma}_{\text{L,h}}}>{{\sigma}_{\text{L,H}}}$.

Besides the fundamental problem of stress heterogeneity and tensorial properties of LV wall stress, Laplace calculations are afflicted with significant uncertainties. The meaning of geometric parameters $r$ and $h$ required for the evaluation of Eqs. (11) or (12) is ambiguous when applied to the LV which deviates in shape markedly from a spherical shell. Therefore, $r$ and $h$ must be determined from averaging over a number of short axis Cine MRI scans to find representative values. Due to longitudinal shortening additional averaging occurs as different slices of the heart are being imaged during ejection. Thus, the determination of parameters $r$ and $h$ cannot be unique as the particular method employed for averaging, such as the one described in Eq. (28), influences, to some extent, the results. Using Eq. (13) seems to circumvent this problem since ${{V}_{\text{cav}}}$ and ${{V}_{\text{myo}}}$ are used as inputs which may be determined uniquely for the LV. However, in our simulations ${{\sigma}_{\text{L,V}}}$ led to larger misestimations than ${{\sigma}_{\text{L,h}}}$ and ${{\sigma}_{\text{L,H}}}$.

It is well known that Laplace-based calculation of stresses is afflicted with various inaccuracies [27]. Nonetheless, Laplace-based calculation of LV wall stresses has been used in clinical studies as a diagnostic criterion [4]. However, according to observations in this study based on an in silico model and in line with other studies [41], the scope for clinical applications appears narrow. Laplace stresses may provide information of diagnostic value, but, if so, rather as an empirical than a mechanistic marker. As a biomarker representing LV wall stresses in a physical sense Laplace-based calculations suffer from severe fundamental limitations.

4.2.2 Mechanical heart power and power efficiency

Mechanical heart power ${{P}_{\text{int}}}$ and cardiac power efficiency ${{P}_{\text{eff}}}$ defined as the ratio between peak mechanical power expended by the LV, ${{\hat{P}}_{\text{int}}}$, and the mean hydrodynamic power delivered to the arterial system, EHP, have been proposed recently as a diagnostic marker [15]. On grounds of conservation of energy the global mechanical power ${{P}_{\text{int}}}$ expended by the LV and the hydrodynamic power transferred to the LV blood pool, ${{P}_{\text{ext}}}$, must be equal. Discrepancies may occur due to isovolumetric phases during which hydrodynamic power is close to zero, but mechanical power is expended by the LV to some extent as conformational changes of the LV myocardium and the shape of the LV cavity occur. However, in all models studied ${{P}_{\text{int}}}$ during isovolumetric phases was negligible (see Fig. 14 in the supplement). This does not conflict with experimental studies providing evidence of heterogeneous circumferential strains, longitudinal shortening and wall thickening during IVC [6]. Qualitatively similar behavior is observed in our in silico models, but magnitude and velocity of strain development is much smaller during IVC than during ejection. Thus, the strain rate tensors $\dot{\boldsymbol{\varepsilon}}$ remained small during IVC and mechanical power expenditure was minor. In all LV models under study global mechanical power ${{P}_{\text{int}}}$ and the hydrodynamic power in the LV cavity ${{P}_{\text{ext}}}$ were virtually identical (see Tab. 5 and Fig. 14 in the supplement). Hence, mechanical heart power can be determined either by analyzing the deformation of the LV myocardium or from PV relations in the LV.

The estimation of ${{P}_{\text{int}}}$ is feasible directly from LV deformation either by using FE models or, as suggested in [15], based on Laplace’s law where the latter approach is more readily applicable in the clinic. However, when global LV power is of interest, Laplace-based approaches do not seem to offer any additional benefits over more standard approaches relying on hemodynamic data for a number of reasons.

First of all, the evaluation of {{P}_{\text{int,\star}}} based on Eq. (17) introduces a systematic error which leads to a misestimation of the actual ${{P}_{\text{int}}}$, even in the spherical shell models, since any work expended in the radial direction is ignored. Laplace’s law takes into account only circumferential stresses and neglects any radial stresses. As shown in Fig. 5, this simplification is only well justified in thin-walled structures such as $\mathrm{\mathbf{Sph}_{5}}$, but introduces pronounced discrepancies for increased $h$ (see passive inflation experiments in Tab. 4 and Tab. 6 in the supplement as well as ${{\overline{\sigma}}_{rr}}$ traces in Fig. 8). Secondly, in addition to the parameters needed for wall stress estimation which are afflicted with substantial uncertainties as discussed above, the parameters $V_{\mathrm{myo}}$ and $\dot{\boldsymbol{\varepsilon}}$ are required. Using the approximation given by Eq. (37) in the supplement the estimation of $\dot{\boldsymbol{\varepsilon}}_{\rm{circ}}$ requires that both inner and outer radii $r$ and $R$ of the LV can be tracked with sufficient temporal resolution and accuracy. However, as evidenced in Fig. 8, even when evaluated in an in silico model where tracking of these quantities is feasible with the highest possible accuracy, the overall accuracy of the method is rather poor with significant under- or overestimation of the true ${{P}_{\text{int}}}$, depending on whether ${{P}_{\text{int,h}}}$, ${{P}_{\text{int,H}}}$ or ${{P}_{\text{int,V}}}$ is used and whether an early or late phase of ejection is considered (see Fig. 8).

The evaluation of cardiac power efficiency ${{P}_{\text{eff}}}$ or ${{P}_{\text{eff,clin}}}$ requires only point estimates of peak mechanical power.

Following [15], this is feasible by assuming that ${{\hat{P}}_{\text{int}}}$ occurs at the instant of peak pressure, $t_{\hat{p}}$. Consistent with expectations based on Laplace’s law this was not the case in any of our LV models. As $p\propto\sigma/\left(r/h\right)$ peak pressure $\hat{p}_{\mathrm{LV}}$ and peak stress would only coincide under isometric conditions. In the contracting LV during ejection, the ratio $\left(r/h\right)$ decreases, thus facilitating a further increase of $p$ beyond the instant of peak pressure (see pressure and stress panels in Fig. 8). Nonetheless, the instants of peak power and peak pressure fell sufficiently close together with $\left|{t_{\hat{p}}-t_{\rm{pow,peak}}}\right|$ of $32$, $12$, $11$ and $9$ ms for $\mathrm{\mathbf{LV}_{A}}$, $\mathrm{\mathbf{LV}_{B}}$, $\mathrm{\mathbf{LV}_{C}}$ and $\mathrm{\mathbf{LV}_{D}}$, respectively. Indeed, inspection of Tab. 5 and the power panels in Fig. 8 suggests that the Laplace-based estimation of ${{\hat{P}}_{\text{int}}}$ seems feasible by evaluating power at the instant of peak pressure (compare ${{\hat{P}}_{\text{int}}}$, ${{P}_{\text{int}}}(t_{\hat{p}})$, ${{P}_{\text{int,h}}}(t_{\hat{p}})$,${{P}_{\text{int,H}}}(t_{\hat{p}})$ and ${{P}_{\text{int,V}}}(t_{\hat{p}})$ in Tab. 5), albeit with inferior accuracy compared to estimations based on hemodynamic PV data.

Alternatively, the simpler ${{P}_{\text{eff,clin}}}$ marker can be used as in [15] which relies on IHP and does not require an estimation of $\dot{\boldsymbol{\varepsilon}}$. While simpler, its use brings about a number of drawbacks. Since $\dot{\boldsymbol{\varepsilon}}$ is ignored IHP is only a relative marker that is non-linearly related to ${{P}_{\text{int}}}$. Therefore, IHP provided highly inconsistent relative estimates of ${{P}_{\text{int}}}$ with errors varying in the range from $24.1\%$ to $140.28\%$ (see Tab. 5). Thus, IHP as an indicator of ${{P}_{\text{int}}}$ appears to be of insufficient accuracy even for clinical applications of modest accuracy demands. Overall, the scope for Laplace-based power estimation as proposed in [15] seems limited as standard methods based on hemodynamic data are afflicted with less uncertainty, offer higher accuracy and are easier to evaluate. As shown in Tab. 5, ${{\hat{P}}_{\text{int}}}$ is straight forwardly approximated – with higher accuracy than any Laplace-based method – as the product of peak pressure and flow, $\hat{p}\cdot\hat{q}$.

The mechanical power generated by the LV is an indicator of metabolic demands. Local wall stresses and power densities governing energetic demand and supply ratios in the LV myocardium are known to play important roles as drivers of remodeling in the pressure-overloaded LV of AS patients. However, analogous to the stresses shown in Fig. 8 the distribution of power density ${{p}_{\text{int}}}(\mathbf{x},t)$ in the LV wall is highly heterogeneous as well with significant regional variability around the global mean power density. In this view Laplace-based global markers derived from mechanical deformation such as ${{\overline{\sigma}}_{\text{circ}}}$ or {{P}_{\text{int,\star}}} are not representative of local stresses and power within the LV myocardium and appear to offer limited insight and predictive power beyond standard PV analysis.

An accurate representation of local mechanical stresses $\boldsymbol{\sigma}(\mathbf{x},t)$ and power ${{p}_{\text{int}}}(\mathbf{x},t)$ over a cardiac cycle depends on reliable sets of strains $\boldsymbol{\varepsilon}(\mathbf{x},t)$. While techniques for measuring strains in 3D throughout the LV myocardium are available [23], such recordings are not part of clinical routine, their analysis requires expensive non-trivial post-processing and spatio-temporal resolution and accuracy are limited. A carefully fitted and validated FE-based EM LV model which replicates a patients physiology in terms of PV relations as well as LV kinematics, provides accurate data on strains $\boldsymbol{\varepsilon}(\mathbf{x},t)$ at a high spatio-temporal resolution. Using an appropriate parameterized patient-specific constitutive model such as given in Eq. (4), $\boldsymbol{\varepsilon}(\mathbf{x},t)$ can be used to compute LV wall stresses $\boldsymbol{\sigma}(\mathbf{x},t)$ and ${{P}_{\text{int}}}$ or any other stress-related biomarker efficiently with high accuracy. Such models are able to provide either global power ${{P}_{\text{int}}}(t)$, but also fine-grained distributed power density ${{p}_{\text{int}}}(\mathbf{x},t)$. A spatio-temporal view on $\boldsymbol{\varepsilon}(\mathbf{x},t)$, $\boldsymbol{\sigma}(\mathbf{x},t)$ and ${{p}_{\text{int}}}(\mathbf{x},t)$ in the LV may provide additional insights as regions of elevated strain, stress or power are assumed to be implicated in the mechanisms driving remodeling in the pressure overloaded LV [16, 34].

5 Conclusions

Laplace estimates of LV wall stress are able to provide a rough approximation of global mean stress in the circumferential-longitudinal plane of the LV. However, according to FE results spatial heterogeneity of stresses in the LV wall is significant, leading to major discrepancies between local stresses and global mean stress. Assessment of mechanical power with Laplace methods is feasible, but these are inferior in accuracy compared to FE models and do not offer any benefits compared to standard methods based on hemodynamic data. In this view, the scope for Laplace-based analysis in clinical applications seems narrow. The accurate assessment of stress and power density distribution in the LV wall is only feasible based on patient-specific FE modeling.

6 Supplementary Material

6.1 Derivation of Laplace’s law

To derive the law of Laplace we consider a spherical shell of inner radius, $r$ and thin walls of a given thickness, $h$. Inflating the shell by applying a pressure, $p$, within the shell’s cavity, induces deformation which causes the buildup of stresses within the wall. If we consider one half of the sphere, the total force $\mathbf{F}_{p}$ acting on the inner surface must be balanced with the total force acting over the cut surface (see Fig. LABEL:sub@supp:fig:sph:bof). Due to spherical symmetry, the circumferential stresses ${{\sigma}_{\text{circ}}}(\rho)$ at any radius $r\leq\rho\leq R$ must be the same and the shear stress is zero. Integrating circumferential stresses over the cut surface yields the total force balancing the force due to the applied pressure. That is, we have

[TABLE]

Assuming that $h\ll r$, radial stresses are small compared to circumferential stresses, ${{\sigma}_{rr}}\ll{{\sigma}_{\text{circ}}}$, and the total stress tensor is approximated by

[TABLE]

with respect to the spherical coordinate system and the projection matrix $\mathbf{P}=(\mathbf{e}_{r},\mathbf{e}_{\varphi},\mathbf{e}_{\theta})^{\top}$, see Fig. LABEL:sub@supp:fig:sph:csys. Note that $\boldsymbol{\sigma}$ in a spherical shell differs from a stress tensor in the LV in various ways. Unlike in the LV, stresses in circumferential and meridonial/longitudinal direction are equal whereas in the LV longitudinal stresses tend to be larger than circumferential stresses. Further, the assumption $r\ll h$ is not justified, rather $r\approx h$ holds. Thus radial stresses in the LV are non-negligible, that is, ${{\sigma}_{rr}}$ is at an order of magnitude comparable to ${{\sigma}_{\text{circ}}}$.

6.1.1 Laplace’s law for a thick-walled sphere

From Eq. (29) the circumferential stress in a thick-walled sphere is found as

[TABLE]

which we denote as ${{\sigma}_{\text{L,H}}}$.

6.1.2 Laplace’s law for a thin-walled sphere

Using assumption (A3), i.e. $h/r\ll 1$, we have $\left(1+\frac{h}{2\,r}\right)\approx 1$ which yields the simple law of Laplace for a thin-walled sphere given by

[TABLE]

which we denote as ${{\sigma}_{\text{L,h}}}$.

6.1.3 Volume-based stress

Since radius $r$ and wall thickness $h$ are two quantities which are not always available or are hard to determine for a general geometry like a LV, we rewrite Eq. (31) in terms of the cavity volume ${{V}_{\text{cav}}}$ and the myocardial volume ${{V}_{\text{myo}}}$. For the spherical shell geometry there holds ${{V}_{\text{cav}}}=\frac{4}{3}\pi r^{3}$ and ${{V}_{\text{cav}}}=\frac{4}{3}\pi\left(R^{3}-r^{3}\right)$ which entails $r=\left(\frac{3\,{{V}_{\text{cav}}}}{4\pi}\right)^{1/3}$ and $(r+h)=\left(\frac{3\,({{V}_{\text{cav}}}+{{V}_{\text{myo}}})}{4\pi}\right)^{1/3}$. Using the volume-based representations of $r$ and $r+h$, we can rewrite Eq. (31)

[TABLE]

which we denote as ${{\sigma}_{\text{L,V}}}$.

6.2 Computation of power and work

The approximations ${{\sigma}_{\text{L,h}}}$, ${{\sigma}_{\text{L,H}}}$ and ${{\sigma}_{\text{L,V}}}$ for the circumferential stress ${{\sigma}_{\text{circ}}}$ can be used to derive an estimator for the internal power

[TABLE]

and internal work

[TABLE]

For this sake, we consider Eq. (34) and the simplified representation of the total stress tensor Eq. 30. In Eq. (34) an approximation of the strain rate $\dot{\boldsymbol{\varepsilon}}$ is required. Rewriting the strain tensor $\boldsymbol{\varepsilon}$ in spherical coordinates, as done for the stress tensor $\boldsymbol{\sigma}$, we obtain

[TABLE]

where $\mathbf{P}$ is the projection matrix introduced in Sec. 6.1. Similarly, the strain rate $\dot{\boldsymbol{\varepsilon}}$ is expressed as

[TABLE]

Using Eq. (30), an approximation of the internal power density, ${{p}_{\text{int}}}$ can be derived as $(\boldsymbol{\sigma}:\dot{\boldsymbol{\varepsilon}})\approx{{\sigma}_{\text{circ}}}\left(\dot{\varepsilon}_{\varphi\varphi}+\dot{\varepsilon}_{\theta\theta}\right)$. Due to the assumption of symmetry (A2), strains in circumferential direction do not vary with space, i,e. $\varepsilon_{\varphi\varphi}=\varepsilon_{\theta\theta}={{\varepsilon}_{\text{circ}}}$, and the approximation for the internal power density simplifies to

[TABLE]

An approximation of the circumferential strain ${{\varepsilon}_{\text{circ}}}$ can be found based the Cauchy strain and considerations illustrated in Fig. 11. Accordingly, for a given radius $r$ circumferential strain can be approximated as

[TABLE]

and for the circumferential strain rate we obtain ${{\dot{\varepsilon}}_{\text{circ}}}\approx\frac{\dot{r}}{r_{0}}$.

To approximate the circumferential strain rate ${{\dot{\varepsilon}}_{\text{circ}}}(\mathbf{u},t)$ of a spherical shell of thickness $h=R-r$, we take the arithmetic mean of the strain rate at inner radius $r$ and outer radius $R$, that is

[TABLE]

where $r_{0}$ is the initial inner radius and $R_{0}$ is the initial outer radius of the spherical shell at its stress free configuration, i.e. $p=0$.

Using Eqs. (36) and (37), the internal power can be estimated by

[TABLE]

where ${{V}_{\text{myo}}}(t)$ is the volume of the shell’s wall at time $t$.

Using an estimate for ${{\sigma}_{\text{circ}}}$ we obtain

[TABLE]

with $\star\in\left\{\text{h},\text{H},\text{V}\right\}$ and by integrating over time, we get an estimate for the internal work

[TABLE]

6.3 Methods

6.3.1 Model fitting

To delineate anisotropy from pure geometry effects, passive inflation experiments were also performed with LV models using the Demiray model and passive mechanical behavior was compared to the spherical shell models $\mathrm{\mathbf{Sph}_{5}}$, $\mathrm{\mathbf{Sph}_{25}}$and $\mathrm{\mathbf{Sph}_{150}}$. In these cases, parameters were set to $b=7$ and $a$ was chosen in a patient-specific manner to obtain the same volume at maximum inflation pressure as with the Guccione model. Note that in none of the simulations of a full cardiac cycle the Demiray model was considered as the resulting kinematics was in stark contrast to the clinical data.

6.3.2 Analysis of LV inflation experiments

To evaluate the influence of violating the assumption on the geometry (A2), passive inflation experiments were performed with LV models and the isotropic material due to Demiray, Eq. (6), following the same protocol as applied to the spherical shell models $\mathrm{\mathbf{Sph}_{5}}$, $\mathrm{\mathbf{Sph}_{25}}$ and $\mathrm{\mathbf{Sph}_{150}}$. The parameter $a$ in the Demiray model was set to $0.45$, $0.63$, $0.41$ and $0.38$ kPa for the models $\mathrm{\mathbf{LV}_{A}}$, $\mathrm{\mathbf{LV}_{B}}$, $\mathrm{\mathbf{LV}_{C}}$ and $\mathrm{\mathbf{LV}_{D}}$, respectively. The Laplace-based stress estimates ${{\sigma}_{\text{L,h}}}$, ${{\sigma}_{\text{L,H}}}$ and ${{\sigma}_{\text{L,V}}}$ were compared to the mean stresses obtained from the FE solution. Stresses were evaluated with respect to an ellipsoidal coordinate system to facilitate a comparison with stresses computed in the spherical shell models $\mathrm{\mathbf{Sph}_{5}}$, $\mathrm{\mathbf{Sph}_{25}}$ and $\mathrm{\mathbf{Sph}_{150}}$ where spherical coordinates were used for stress analysis. The ellipsoidal coordinate system for the LV models was constructed by assigning fiber and sheet orientations using a rule-based method with a constant fiber angle of $0^{\circ}$. Stress components ${{\sigma}_{rr}}(\mathbf{x})$, ${{\sigma}_{\varphi\varphi}}(\mathbf{x})$ and ${{\sigma}_{\theta\theta}}(\mathbf{x})$ were averaged yielding ${{\overline{\sigma}}_{rr}}$, ${{\overline{\sigma}}_{\varphi\varphi}}$ and ${{\overline{\sigma}}_{\theta\theta}}$, respectively. Note that all models except $\mathrm{\mathbf{Sph}_{5}}$ showed marked spatial stress variations. Thus, the reported mean stresses $\bar{\sigma}$ may deviate considerably from the true local stresses $\sigma(\mathbf{x})$. Laplace-based estimations of power {{P}_{\text{int,\star}}} and work {{W}_{\text{int,\star}}}, were compared to those obtained by FE simulation, ${{P}_{\text{int}}}$ and ${{W}_{\text{int}}}$ and to external hydrodynamic power and work in the LV cavity, ${{P}_{\text{ext}}}$ and ${{W}_{\text{ext}}}$.

6.4 Results

6.4.1 Verification of the FE model

Similarly, with increasing $h$ the accuracy of the thick-walled Laplace estimate ${{W}_{\text{int,H}}}$ performed better than the simpler thin-walled Laplace estimate ${{W}_{\text{int,h}}}$. As expected on grounds of conservation of energy, the agreement between biomechanical work ${{W}_{\text{int}}}$ and hemodynamic work ${{W}_{\text{ext}}}$ was essentially perfect with differences $<2\%$ for all models.

6.4.2 Passive inflation of LV models

The LV models $\mathrm{\mathbf{LV}_{A}}$–$\mathrm{\mathbf{LV}_{D}}$ were inflated following the loading protocol in Fig. 12.A. Passive material behavior was represented compliant with (A1) by the isotropic Demiray model. The temporal evolution of FE- and Laplace-based stresses, power and work are shown in Fig. 12.B for model $\mathrm{\mathbf{LV}_{D}}$. Minor quantitative differences to other models $\mathrm{\mathbf{LV}_{A}}$–$\mathrm{\mathbf{LV}_{C}}$ were observed, but qualitatively the overall behavior was identical. Stresses at $p=4$ kPa and the amount of work incurred during inflation up to this pressure are summarized in Tabs. 6 and  7.

6.4.3 Analysis of LV cycle experiments