Multiscale Modeling Framework of Ventricular-Arterial Bi-directional Interactions in the Cardiopulmonary Circulation

Ventricular-arterial coupling plays a key role in the physiologic function of the cardiovascular system. We have previously described a hybrid lumped-finite element (FE) modeling framework of the systemic circulation that couples idealized FE models of the aorta and the left ventricle (LV). Here, we describe an extension of the lumped-FE modeling framework that couples patient-specific FE models of the left and right ventricles, aorta and the large pulmonary arteries in both the systemic and pulmonary circulations. Geometries of the FE models were reconstructed from magnetic resonance (MR) images acquired in a pediatric patient diagnosed with pulmonary arterial hypertension (PAH). The modeling framework was calibrated with pressure waveforms acquired in the heart and arteries by catheterization as well as ventricular volume and arterial diameter waveforms measured from MR images. The calibrated model hemodynamic results match well with the clinically-measured waveforms (volume and pressure) in the LV and right ventricle (RV) as well as with the clinically-measured waveforms (pressure and diameter) in the aorta and main pulmonary artery. The calibrated framework was then used to simulate three cases, namely, (1) an increase in collagen in the large pulmonary arteries, (2) a decrease in RV contractility, and (3) an increase in the total pulmonary arterial resistance, all characteristics of progressive PAH. The key finding from these simulations is that hemodynamics of the pulmonary vasculature and RV wall stress are more sensitive to vasoconstriction with a 10% of reduction in the lumen diameter of the distal vessels than a 67% increase in the proximal vessel's collagen mass.

Introduction

Ventricular-arterial coupling plays a vital role in the physiologic function of the cardiopulmonary circulation as well as in the evolution of cardiovascular diseases, such as pulmonary arterial hypertension (PAH) (Borlaug and Kass, 2011; Ky et al., 2013). In physiologic conditions, the arterial compliance (endowed by arterial wall tissue constituents) and the ventricular dynamic stiffness (inherent from the contraction of myocytes) confine the dynamic pressure variation to a physiological range to prevent end organ damage, while providing sufficient blood flow to meet oxygen demand of the body under varying workload (Borlaug and Kass, 2011). In pathological conditions, such as PAH, malfunction of one compartment (e.g., microcirculation) in the cardiopulmonary circulation may affect other compartments (e.g., ventricle) through a positive feedback loop that is driven by the tight coupling of ventricular and arterial systems, ultimately leading to end-stage heart failure. A modeling framework that captures the complex ventricular-arterial coupling would help elucidate the mechanisms governing the progression of PAH.

Existing mathematical modeling frameworks describing ventricular-arterial coupling in the cardiopulmonary circulation can be broadly classified as either a lumped parameter or a multi-scale finite element (FE) modeling framework. In a lumped parameter modeling framework, the ventricular-arterial coupling is described by an electrical analog representation of the cardiovascular system (Ursino, 1998; Smith et al., 2004). While such modeling framework is computationally inexpensive, it cannot directly take into account detailed geometrical and microstructural features associated with pathological conditions in the ventricles and arteries. In a hybrid lumped-FE modeling framework, a FE model describing either ventricular mechanics (Kerckhoffs et al., 2007; Shavik et al., 2017, 2019) or arterial hemodynamics (Lau and Figueroa, 2015; Zambrano et al., 2018) is coupled to lumped-parameter representation of the other compartments to provide a detailed description of the cardiovascular system. To overcome limitations associated with simplified representations of cardiovascular components, we previously introduced a hybrid lumped-FE modeling framework that bidirectionally couples FE models of the aorta and left ventricle (LV) mechanics in a closed-loop circulatory system (Shavik et al., 2018). Based on an idealized geometry of the LV and aorta, the modeling framework is able to reproduce pressure, arterial diameter, and LV volume waveforms found in a healthy individual. The modeling framework, however, considers only the systemic circulation and does not take into account the pulmonary circulation.

Here, we describe the extension of our earlier framework (Shavik et al., 2018) in which image-based FE models of the large pulmonary arteries, aorta, and heart (including both ventricles) are coupled bidirectionally in a closed-loop multi-scale FE modeling framework of the cardiopulmonary circulation. The multi-scale framework was calibrated using in vivo clinical measurements of the anatomy, deformation, and hemodynamics from a PAH pediatric patient. Using the calibrated model, we further investigate how changes associated with the mechanical behavior and microstructure of the microcirculation, large pulmonary arteries, and right ventricle (RV), consequent of PAH progression, affect each other.

Methods

This study was approved by the University of Michigan Board of Review (HUM00117706), and informed consent was obtained from the parents/guardians of the patient.

Patient History

Clinical data was prospectively acquired in a 11-year-old female patient who was diagnosed with PAH. The patient had an elevated mean pulmonary arterial pressure (mPAP) of 59 mmHg with normal pulmonary capillary wedge pressure (PCWP) of 6 mmHg and elevated pulmonary vascular resistance (PVR) of 13.5 WU, falling within the clinical classification of PAH (mPAP ≥ 20 mmHg, PCWP ≤ 15 mmHg, and PVR ≥ 3 WU) (Simonneau et al., 2019). She has family history of chronic obstructive pulmonary disease and PAH.

Data Acquisition

Anatomical and hemodynamic data were obtained using magnetic resonance (MR) imaging and arterial catheterization. Cine MR images of the short- and long-axis views of the ventricles were acquired at 30 time points in the cardiac cycle. Using the cine MR images, left and right ventricular endocardial surfaces were segmented with the medical image analysis software MeVisLab (www.mevislab.de) to acquire ventricular volume waveforms. Cardiac-gated gradient echo MR images of the vascular anatomy were acquired in the diastolic phase. Luminal area waveforms were also acquired with phase-contrast MR images (PC-MRI) at the ascending aorta and main pulmonary artery. Arterial catheterization was performed to acquire pressure waveforms in the LV, RV, main pulmonary artery (MPA), and aorta. The ventricular volume and pressure waveforms were synchronized to reconstruct pressure-volume (PV) loops (Xi et al., 2016; Shavik et al., 2019). Hemodynamic and cardiovascular function metrics of the PAH patient are listed in Table 1.

TABLE 1

Table 1. Hemodynamic measurements of PAH patient.

Biventricular and Vascular Geometries

Anatomical models of the LV, RV, aorta, and pulmonary arteries (PA) (consisting of the main, left, and right pulmonary arteries) were reconstructed from the acquired MR images. The biventricular model was reconstructed from images that correspond to the point in the cardiac cycle where ventricular pressures were lowest during filling (Geuzaine and Remacle, 2009). Furthermore, anatomical models of the aorta and large pulmonary arteries were reconstructed using the blood flow modeling software CRIMSON (www.crimson.software) (Figure 1).

FIGURE 1

Figure 1. Reconstruction of biventricular model (left) and large proximal arteries (right) from cine MR images.

Closed Loop Circulatory System

The biventricular, aorta and pulmonary artery FE models were coupled through a closed loop lumped-parameter circulatory model that describes both systemic and pulmonary circulations (Figure 2). The modeling framework consists of eight compartments with four cardiovascular components (ventricle, atrium, artery, and vein) each in the systemic and pulmonary circulations. Conservation of total blood mass in the circulatory model requires the net change of inflow and outflow rates of each compartment to be related to the rate of change of the volume by the following relations

$\begin{array}{ll}\frac{d{V}_{LV}(t)}{dt}=q_{mv}(t)-q_{av}(t),& (\mathrm{1a})\end{array}$

$\begin{array}{ll}\frac{d{V}_{sa}(t)}{dt}=q_{av}(t)-q_{sa}(t),& (\mathrm{1b})\end{array}$

$\begin{array}{ll}\frac{d{V}_{sv}(t)}{dt}=q_{sa}(t)-q_{sv}(t),& (\mathrm{1c})\end{array}$

$\begin{array}{ll}\frac{d{V}_{RA}(t)}{dt}=q_{sv}(t)-q_{tv}(t),& (\mathrm{1d})\end{array}$

$\begin{array}{ll}\frac{d{V}_{RV}(t)}{dt}=q_{tv}(t)-q_{pvv}(t),& (1e)\end{array}$

$\begin{array}{ll}\frac{d{V}_{pa}(t)}{dt}=q_{pvv}(t)-q_{pa}(t)& (\mathrm{1f})\end{array}$

$\begin{array}{ll}\frac{d{V}_{pv}(t)}{dt}=q_{pa}(t)-q_{pv}(t),& (\mathrm{1g})\end{array}$

$\begin{array}{ll}\frac{d{V}_{LA}(t)}{dt}=q_{pv}(t)-q_{mv}(t).& (\mathrm{1h})\end{array}$

FIGURE 2

Figure 2. Schematic of the finite element ventricular-vascular coupling closed loop circulatory modeling framework.

In Equation (1), V_LV, V_sa, V_sv, V_RA, V_RV, V_pa, V_pv, and V_LA are the volumes of the eight compartments with the subscripts denoting the LV, systemic arteries (sa), systemic veins (sv), right atrium (RA), RV, pulmonary arteries (pa), pulmonary veins (pv), and left atrium (LA), respectively. Flow rates at different segments of the circulatory model are denoted by q_mv, q_av, q_sa, q_sv, q_tv, q_pvv, q_pa, and q_pv.

Systemic and pulmonary arteries and veins were modeled using their electrical analogs based on Ohm's law. At each segment, the flow rate depends on the pressure gradient and resistance to the flow as described in the following equation

$\begin{array}{ll}{q}_{mv}(t)=\{\begin{array}{ll}\frac{P_{LA}(t)-P_{LV}(t)}{R_{mv}}\hfill & when,P_{LA}(t)\ge P_{LV}(t)\hfill \\ 0\hfill & when,P_{LA}(t)<P_{LV}(t)\hfill \end{array}& (\mathrm{2a})\end{array}$

$\begin{array}{ll}{q}_{av}(t)=\{\begin{array}{ll}\frac{P_{LV}(t)-P_{sa}(t)}{R_{av}}\hfill & when,P_{LV}(t)\ge P_{sa}(t)\hfill \\ 0\hfill & when,P_{LV}(t)<P_{sa}(t)\hfill \end{array},& (\mathrm{2b})\end{array}$

$\begin{array}{ll}{q}_{sa}(t)=\frac{P_{sa}(t)-P_{sv}(t)}{R_{sa}},& (\mathrm{2c})\end{array}$

$\begin{array}{ll}{q}_{sv}(t)=\frac{P_{sv}(t)-P_{RA}(t)}{R_{sv}},& (\mathrm{2d})\end{array}$

$\begin{array}{ll}{q}_{tv}(t)=\{\begin{array}{ll}\frac{P_{RA}(t)-P_{RV}(t)}{R_{tv}}\hfill & when,P_{RA}(t)\ge P_{RV}(t)\hfill \\ 0\hfill & when,P_{RA}(t)<P_{RV}(t)\hfill \end{array},& (2e)\end{array}$

$\begin{array}{ll}{q}_{pvv}(t)=\{\begin{array}{ll}\frac{P_{RV}(t)-P_{pa}(t)}{R_{pvv}}\hfill & when,P_{RV}(t)\ge P_{pa}(t)\hfill \\ 0\hfill & when,P_{RV}(t)<P_{pa}(t)\hfill \end{array},& (\mathrm{2f})\end{array}$

$\begin{array}{ll}{q}_{pa}(t)=\frac{P_{pa}(t)-P_{pv}(t)}{R_{pa}},& (\mathrm{2g})\end{array}$

$\begin{array}{ll}{q}_{pv}(t)=\frac{P_{pv}(t)-P_{LA}(t)}{R_{pv}}.& (\mathrm{2h})\end{array}$

In Equation (2), R_mv, R_av, R_tv, and R_pvv are the resistances associated with the mitral, aortic, tricuspid, and pulmonary valves, respectively. The valves are each represented by a diode that only permits one-way flow as in previous studies (Punnoose et al., 2012; Shavik et al., 2019). The vessel resistances are denoted by R_sa, R_sv, R_pa, and R_pv, respectively. To describe the compliance of the systemic and pulmonary vessels, we used the following PV relationships

$\begin{array}{ll}{P}_{sv}(t)=\frac{V_{sv}(t)-V_{sv,0}}{C_{sv}},& (\mathrm{3a})\end{array}$

$\begin{array}{ll}{P}_{pv}(t)=\frac{V_{pv}(t)-V_{pv,0}}{C_{pv}},& (\mathrm{3b})\end{array}$

where V_sv,0 and V_pv,0 are the resting volumes and C_sv and C_pv are the total compliance of the systemic and pulmonary veins, respectively.

Contraction of the LA and RA was modeled using a time varying elastance function that is given by the following PV relations

$\begin{array}{ll}{P}_k(t)=e(t)P_{es,k}(V_k(t))+(1-e(t))P_{ed,k}(V_k(t)),& (\mathrm{4a})\end{array}$

where,

$\begin{array}{ll}{P}_{es,k}(V_k(t))=E_{es,k}(V_k(t)-V_{0,k}),& (\mathrm{4b})\end{array}$

$\begin{array}{ll}{P}_{ed,k}(V_k(t))=A_k(e^{B_k(V_k(t)-V_{0,k})}-1).& (\mathrm{4c})\end{array}$

In Equation (4), the subscript k denotes either LA or RA. The volume, end-systolic elastance, and volume-intercept of the end-systolic pressure-volume relationship (ESPVR) of the corresponding atrium are denoted by V_k, E_es,k, and V_0,k, respectively. The parameters A_k and B_k define the atrium curvilinear end-diastolic pressure volume relationship (EDPVR) and the driving function is defined as

$\begin{array}{ll}e(t)=\{\begin{array}{ll}\frac{1}{2}\left(\sin \left[\left(\frac{\pi }{t_{max}}\right)t-\frac{\pi }{2}\right]+1\right);\hfill & 0<t\le 3t_{max}/2\hfill \\ \frac{1}{2}{e}^{\frac{-(t-\frac{3t_{max}}{2})}{\tau }};\hfill & t>3t_{max}/2,\hfill \end{array},& (5)\end{array}$

where t_max is the point of maximal chamber elastance and τ is the time constant of relaxation. The time-varying elastance model has been shown to be able to describe atrium contraction well (Hoit et al., 1994).

The relationships between pressures and volumes in the biventricular unit (i.e., LV and RV), pulmonary artery and aorta were computed from their corresponding FE models. These relationships can be expressed as non-closed form functions.

$\begin{array}{ll}{P_{RV}(t),P}_{LV}(t)=f^{BV}(V_{LV}(t),V_{RV}(t)),& (\mathrm{6a})\end{array}$

$\begin{array}{ll}{P}_{pa}(t)=f^{PA}(V_{pa}(t)),& (\mathrm{6b})\end{array}$

$\begin{array}{ll}{P}_{sa}(t)=f^{AO}(V_{sa}(t)).& (\mathrm{6c})\end{array}$

Finite Element Formulation of the Biventricular Unit

The weak form associated with the biventricular FE model was derived based on minimization of the following Lagrangian functional

$\begin{array}{l}{{L}}_{BV}\left(u_{BV},p_{BV},P_{LV},P_{RV},c_{1,BV},c_{2,BV}\right)\\ ={\displaystyle {\int }_{{\Omega }_{0,BV}}{W}_{BV}(u_{BV})dV}-{\displaystyle {\int }_{{\Omega }_{0,BV}}{p}_{BV}(J_{BV}-1)dV}\\ -P_{LV}\left(V_{LV,\text{cav}}(u_{BV})-V_{LV}\right)+P_{RV}\left(V_{RV,\text{cav}}(u_{BV})-V_{RV}\right)\\ -c_{1,BV}·{\displaystyle {\int }_{{\Omega }_{0,BV}}{u}_{BV}dV}-c_{2,BV}·{\displaystyle {\int }_{{\Omega }_{0,BV}}{X}_{BV}\times u_{BV}dV},& (7)\end{array}$

where Ω_0,BV is the reference configuration of the biventricular unit, u_BV is the displacement field, P_LV and P_RV are, respectively, the Lagrange multipliers that constrain the LV cavity volume V_{LV, cav}(u_BV) to a prescribed value V_LV and the RV cavity volume V_{RV, cav}(u_BV) to a prescribed value V_RV (Pezzuto and Ambrosi, 2014). We note that V_LV and V_RV are prescribed from the closed-loop circulatory model in Equation (6). The Lagrange multiplier p_BV was used to enforce incompressibility of the tissue (i.e., Jacobian of the deformation gradient tensor J = 1). The vectors c_1,BV and c_2,BV are Lagrange multipliers applied to constrain, respectively, the rigid body translation (i.e., zero mean translation) and rotation (i.e., zero mean rotation) (Pezzuto et al., 2014). In Equation (7), X_BV denotes a material point in Ω_0,BV and W_BV is the strain energy function of the myocardial tissue. The cavity volume of the LV and RV were obtained from the displacement field by using the following functional relationship (k = LV or RV)

$\begin{array}{ll}{V}_{k,cav}(u_{BV})={\displaystyle {\int }_{{\Omega }_{inner,k}}d{V}_k}=-\frac{1}{3}{\displaystyle {\int }_{{\Gamma }_{inner,k}}{x}_{BV}.nd{a}_k},& (8)\end{array}$

where Ω_inner,k is the volume enclosed by the inner surface Γ_inner,k of the LV or RV, and n denotes the outward unit normal vector of those surfaces. Taking the first variation of the Lagrangian functional given in Equation (7) leads to

$\begin{array}{l}\delta {{L}}_{BV}={\displaystyle {\int }_{{\Omega }_{0,BV}}(P_{BV}-p_{BV}{F}_{BV}{}^{-\text{T}}):\nabla \delta u_{BV}dV}\\ -{\displaystyle {\int }_{{\Omega }_{0,BV}}\delta p_{BV}(J-1)dV}-(P_{LV}\\ +P_{RV})-{\displaystyle {\int }_{{\Omega }_{0,BV}}cof(F_{BV}):\nabla \delta u_{BV}dV}-\delta P_{LV}(V_{LV,\text{ cav}}(u_{BV})\\ -V_{LV})-\delta P_{RV}(V_{RV,\text{ cav}}(u_{BV})-V_{RV})\\ -\delta c_{1,BV}·{\displaystyle {\int }_{{\Omega }_{0,BV}}{u}_{BV}dV}-\delta c_{2,BV}·{\displaystyle {\int }_{{\Omega }_{0,BV}}{X}_{BV}\times u_{BV}dV}\\ -c_{1,BV}·{\displaystyle {\int }_{{\Omega }_{0,BV}}\delta u_{BV}dV}-c_{2,BV}·{\displaystyle {\int }_{{\Omega }_{0,BV}}{X}_{BV}\times \delta u_{BV}dV}.& (9)\end{array}$

In Equation (9), P_BV is the first Piola Kirchhoff stress tensor and F_BV is the deformation gradient tensor. The variations of the displacement field, Lagrange multiplier for enforcing incompressibility and volume constraint, zero mean translation, and rotation are denoted by δu_BV, δp_BV, δP_{LV, cav}, δP_RV,cav, δc_1,BV, and δc_2,BV, respectively. Together with the constraint that the basal deformation at z = 0 is in-plane in the biventricular unit, the solution of the Euler-Lagrange problem was obtained by finding $u_{BV}{\in H}^1({\Omega }_0)$, $p_{BV}\in L^2({\Omega }_0)$, $P_{LV,\mathrm{\text{cav}}}\in ℝ,P_{RV,\mathrm{\text{cav}}}\in ℝ,c_{1,BV}\in {ℝ}^3,c_{2,BV}\in {ℝ}^3$ that satisfies

$\begin{array}{ll}{\delta {L}}_{BV}=0,& (\mathrm{10a})\end{array}$

$\begin{array}{ll}{u}_{BV}\left(x,y,0\right).\mathrm{n}{|}_{base}=0,& (\mathrm{10b})\end{array}$

for all $\delta u_{BV}{\in H}^1({\Omega }_0),\delta p_{BV}\in L^2({\Omega }_0)$, δP_{LV, cav}∈ ℝ, δP_{RV, cav}∈ ℝ, ${\delta c}_{1,BV}\in {ℝ}^3$, $\delta c_{2,BV}\in {ℝ}^3.$ The solution of Equation (10) gives the relationship between P_RV, P_LV, V_RV, V_LV in Equation (6).

Mechanical Behavior of the Cardiac Tissue

Mechanical behavior of the myocardial tissue was described by an active stress formulation in which the first Piola-Kirchhoff stress tensor P_BV in Equation (9) was additively decomposed into a passive and an active component, i.e.,

$\begin{array}{ll}{P}_{BV}=P_{BV,p}+P_{BV,a}{e}_f\otimes e_{f_0}.& (11)\end{array}$

In Equation (11), P_{BV, p} is the passive stress tensor, P_{BV, a} is the magnitude of the active stress, whereas e_f and e_f0 are the local basis vectors that define the cardiac muscle fiber directions in the current and reference configuration, respectively. The passive stress tensor P_{BV, p} is related to the strain energy function W_{BV, p} and deformation gradient tensor F_BV by

$\begin{array}{ll}{P}_{BV,p}=\frac{d{W}_{BV,p}}{d{F}_{BV}}.& (12)\end{array}$

A Fung-type transversely-isotropic hyperelastic strain energy function (Guccione et al., 1991)

$\begin{array}{ll}{W}_{BV,p}=\frac{1}{2}{C}_{BV}(e^Q-1),& (\mathrm{13a})\end{array}$

with

$\begin{array}{l}Q=b_{ff}{E}_{ff}^2+b_{xx}(E_{ss}^2+E_{nn}^2+E_{sn}^2+E_{ns}^2)\\ +b_{fx}\left(E_{fn}^2+E_{nf}^2+E_{fs}^2+E_{sf}^2\right)& (\mathrm{13b})\end{array}$

was prescribed. In Equation (13b), E_ij with (i, j) ϵ (f, s, n) denote the components of the Green-Lagrange strain tensor $\mathit{E}=\frac{1}{2}(F_{BV}{}^T{F}_{BV}-I)$ with f, s, n denoting the myofiber, sheet and sheet normal directions, respectively. Material parameters of the Fung-type constitutive model are C_BV, b_ff, b_xx, and b_fx.

To describe the active stress behavior, a previously developed active contraction model (Kerckhoffs et al., 2003) was used. The magnitude of the active stress P_{BV, a} was described by

$\begin{array}{ll}{P}_{BV,a}=\frac{l_s}{l_{s0}}{f}^{iso}\left(l_c\right)f^{twitch}\left({t,l}_s\right)\left(l_s-l_c\right)E_a,& (14)\end{array}$

where l_s is the sarcomere length, l_c is the length of the contractile element, l_s0 is the sarcomere length in a prescribed reference state (relaxed sarcomere length), and E_a is the stiffness of the serial elastic element. The function $f^{iso}(l_c)$ denotes the dependency of the isometrically developed active stress on l_c and is given by

$\begin{array}{ll}{f}^{iso}\left(l_c\right)=\{\begin{array}{ll}{T}_0{\text{tanh}}^2[{\text{a}}_6({\text{l}}_{\text{c}}-{\text{a}}_7)]\hfill & when,l_c<a_7\hfill \\ \hfill 0\hfill & when,l_c>a_7\hfill \end{array},& (15)\end{array}$

where T₀ is a model parameter that scales the active tension. Both a₆ and a₇ are model parameters. The time course of the active tension development is controlled by

$\begin{array}{ll}{f}^{twitch}\left(t,l_s\right)=\{\begin{array}{ll}0\hfill & when,t<0\hfill \\ \text{tan}{\text{h}}^2(\frac{t}{t_r}){\text{tanh}}^2(\frac{t_{max-}t}{t_d})\hfill & when,0<t<t_{max}\hfill \\ 0\hfill & when,t>0,\hfill \end{array}& (\mathrm{16a})\end{array}$

$\begin{array}{ll}{t}_{max}=b(l_s-l_d).& (\mathrm{16b})\end{array}$

In Equation (16), t_r is the activation rise time constant, t_d is the activation decay time constant, b relates activation duration t_max to the sarcomere length l_s, and l_d is the sarcomere length at the start of the activation time, i.e., when t_max = 0. The time course of the contractile element l_c was expressed by an ordinary differential equation

$\begin{array}{ll}\frac{\partial l_c}{\partial t}=\left[E_a\left(l_s-l_c\right)-1\right]v_0,& (17)\end{array}$

where v₀ is the unloaded shortening velocity. The sarcomere length l_s was calculated from the myofiber stretch λ and the relaxed sarcomere length l_s0 by

$\begin{array}{ll}\lambda =\sqrt{e_{f_0}·F_{BV}{}^T{F}_{BV}{e}_{f_0}},& (\mathrm{18a})\end{array}$

$\begin{array}{ll}{l}_s=\lambda l_{s0}.& (\mathrm{18b})\end{array}$

Finite Element Formulation of the Arteries

The pulmonary artery and aorta were modeled as 3D membranes. In the formulation that follows, the subscript k = AO denotes the aorta and k = PA denotes the pulmonary artery. Similar to that of the biventricular unit, the finite element formulation of these two arteries can be generalized from the minimization of the following Lagrangian functional, described in the following equation

$\begin{array}{l}{{L}}_k(u_k,P_{k,\text{cav}},c_{1,k},c_{2,k})\\ ={\displaystyle {\int }_{{\Omega }_{0,k}}{W}_k(u_k)dV}-P_{k,\text{cav}}(V_{k,\text{cav}}(u_k)-V_k)\\ -c_{1,k}·{\displaystyle {\int }_{{\Omega }_{0,k}}{u}_{k}dV}-c_{2,k}·{\displaystyle {\int }_{{\Omega }_{0,k}}{X}_k\times u_{k}dV},& (19)\end{array}$

where Ω_0,k is the reference configuration of the arteries, u_k is the displacement field and P_k,cav is the Lagrange multiplier that constrains the arterial cavity volume V_k,cav(u_k) to a prescribed value V_k. The vectors c_1,k and c_2,k are Lagrange multipliers applied to constrain rigid body motions. The inlet and outlets of the arteries were constrained to move only in-plane. Therefore, the solution of the Euler-Lagrange problem was obtained by finding $u_k{\in H}^1({\Omega }_0),P_{k,\mathrm{\text{cav}}}\in ℝ,c_{1,k}\in {ℝ}^3,c_{2,k}\in {ℝ}^3$that satisfies

$\begin{array}{ll}{\delta {L}}_k=0,& (\mathrm{20a})\end{array}$

$\begin{array}{ll}{u}_k\left(x,y,0\right).\mathrm{n}{|}_{inlet,outlets}=0,& (\mathrm{20b})\end{array}$

for all $\delta u_k{\in H}^1({\Omega }_0)$, $\delta P_{k,\mathrm{\text{cav}}}\in ℝ,\delta c_{1,k}\in {ℝ}^3,\delta c_{2,k}\in {ℝ}^3.$ The solution above gives the relationships between P_pa, V_pa, and P_sa, V_sa in Equations (6b) and (6c), respectively.

Mechanical Behavior of the Vascular Tissue

The mechanical behavior of the arteries were described by the strain energy function W_k in Equation (21), which is given as the sum of the key tissue constituents, namely, elastin-dominated matrix W_k,e, collagen fiber families W_k,c,i and vascular smooth muscle cells (SMC) W_k,m (Baek et al., 2007; Zeinali-Davarani et al., 2011), i.e.,

$\begin{array}{ll}{W}_k=W_{k,e}+{\displaystyle \sum _{i=1}^4}{W}_{k,c,i}+W_{k,m}.& (21)\end{array}$

Strain energy function of the elastin-dominated amorphous matrix in the arteries is given by

$\begin{array}{ll}{W}_{k,e}=M_{k,e}\left(\frac{C_{k,1}}{2}\right)(\text{tr}(C_k)-3),& (22)\end{array}$

where M_k,e is the mass per unit volume of the elastin in the tissue, C_k,1 is a stiffness parameter and, ${\mathbf{\text{C}}}_{\mathit{\text{k}}}={\mathbf{\text{F}}}_{\mathit{\text{k}}}^{\mathbf{\text{T}}}{\mathbf{\text{F}}}_{\mathit{\text{k}}}$ is the right Cauchy-Green deformation tensor associated with the arteries.

In the membrane models, four collagen fiber families were considered. The first and second families of collagen fibers (i = 1 and 2) were oriented in the longitudinal and circumferential directions, whereas the third and fourth families of collagen fibers (i = 3 and 4) were oriented, respectively, at an angle α = 45° and −45° with respect to the longitudinal axis based on a previous study (Zeinali-Davarani et al., 2011). We assumed that the same strain energy function for all the families of collagen fibers is given by

$\begin{array}{ll}{W}_{k,c,i}=M_{k,i}\frac{c_{k,2}}{{4c}_{k,3}}\left\{exp\left[c_{k,3}{({{\lambda }_{k,i}}^2-1)}^2\right]-1\right\}.& (23)\end{array}$

In Equation (23), M_k,i is the mass per unit volume of ith family of collagen fibers, λ_k,i is the corresponding stretch of those fibers, and both c_k,2 and c_k,3 are the material parameters that govern the collagen stiffness. The stretch in the ith family of collagen fibers was defined by ${\lambda }_{k,i}=\sqrt{e_{k,i0}·C_k{e}_{k,i0}}$ where e_k,i0 is the local unit vector that defines the corresponding fiber orientation.

Strain energy function of the SMC W_k,m is given by

$\begin{array}{ll}{W}_{k,m}=M_{k,m}\frac{c_{k,4}}{{4c}_{k,5}}\left\{exp\left[c_{k,5}{({{\text{λ}}_{k,m}}^2-1)}^2\right]-1\right\}.& (24)\end{array}$

Here, M_k,m is the mass per unit volume of the SMC in the tissue, λ_k,m is the stretch of the SMC, whereas c_k,4 and c_k,5 are the stiffness parameters. The SMC was assumed to be aligned in the circumferential direction. Mass per unit volume for the different constituents were calculated using following relations

$\begin{array}{ll}{M}_{k,e}={\varphi }_{k,e}\rho ,& (\mathrm{25a})\end{array}$

$\begin{array}{ll}{M}_{k,m}={\varphi }_{k,m}\rho ,& (\mathrm{25b})\end{array}$

$\begin{array}{ll}{M}_{k,i}={\varphi }_{k,i}(1-{\varphi }_{k,e}-{\varphi }_{k,m})\rho ,& (\mathrm{25c})\end{array}$

where ϕ_k,e, ϕ_k,m, ϕ_k,i denote the mass fraction for elastin, SMC and ith family of collagen fibers, respectively. Twenty percent of the total collagen mass is assumed to be equally distributed in the longitudinal and circumferential fiber families and the remaining 80% was distributed equally to α = 45° and −45° fiber directions. Constitutive parameters, mass fraction of each constituent and other parameters of the pulmonary artery and aorta membrane models are listed in Table 2.

TABLE 2

Table 2. Model parameters for FE models for the baseline case.

Solution Algorithm

An explicit time integration scheme was used to solve the ODEs in Equation (1). Specifically, compartment volumes (V_LV, V_sa, V_sv, V_RA, V_RV, V_pa, V_pv, V_LA) at each time t_i were determined from their respective values and the segmental flow rates (q_mv, q_av, q_sa, q_sv, q_tv, q_pvv, q_pa, q_pv) at previous time t_i−1 in Equation (2). The computed compartment volumes at t_i were used to update the corresponding pressures (P_LA, P_RA, P_LV, P_RV, P_sa, P_pa, P_sv, P_pv). Pressures in the atrium (P_LA, P_RA) and veins (P_sv, P_pv) were computed from Equations (4) and (3), respectively. On the other hand, pressures in the LV (P_LV), RV (P_RV), were computed from the FE solutions of Equation (10) for the biventricular unit with the volumes (V_LV, V_RV) at time t_i as input. Similarly, pressures in the aorta (P_sa) and pulmonary artery (P_pa) were computed from the FE solutions of Equation (20) with their corresponding volumes (V_sa, V_pa) at time t_i. We note here that (P_LV, P_RV, P_sa, P_pa) are scalar Lagrange multipliers in the FE formulation for constraining the cavity volumes to the prescribed values (V_LV, V_RV, V_sa, V_pa). The computed pressures at time t_i were then used to update the segmental flow rates in Equation (2) that will be used to compute the compartment volumes at time t_i+1 in the next iteration.

Model Parameterization and Simulation

The biventricular FE model was divided into three material regions, namely the LV free wall (LVFW), the septum, and the RV free wall (RVFW). Similar to a previous study (Finsberg et al., 2018), passive stiffness C and contractility T₀ were prescribed to be the same values in the LVFW and septum (denoted as C_LV and T_0,LV) and had different values in the RVFW (denoted as C_RV and T_0,RV). In the baseline case, model parameters were adjusted to fit the clinically measured LV and RV PV loops, volume and pressure waveforms throughout the cardiac cycle. Specifically, the LV and RV end diastolic pressures were matched by adjusting the passive stiffness parameters C_LV and C_RV. Stroke volume (SV) of the LV and RV were matched by adjusting the regional contractility parameters (i.e., T_0,LV, T_0,RV). While other model parameters can also affect the SV (e.g., peripheral resistances R_sa and R_pa of the systemic and pulmonary circulations as well as preload), the parameters T_0,LV and T_0,RV, which scale the active stress generated by the myofiber, have a larger effect on the LV and RV SV, respectively. On the other hand, the contraction model parameters t_r, t_d and b were adjusted to match the time course of the volume and pressure waveforms measured in the LV and RV. Parameters t_r and t_d were adjusted to match the time to peak tension and b was adjusted to achieve the desirable relaxation of the myofibers. Circulatory model parameters (resistances and compliances) were also adjusted to match the systolic pressure (afterload), preload and systemic and pulmonary vein pressures. Aortic and PA peripheral resistances (R_sa, R_pa) were calibrated to match the systolic pressures of LV and RV. The parameters related to LA and RA time-varying elastance models were prescribed based on a previous study (Shavik et al., 2019). Parameters related to the aorta and PA constitutive models (that alter the vessel's compliance) were adjusted to match the measured pressure waveforms, and the diameters estimated from the PC-MRI. All the model parameters for the biventricular, aorta and PA FE models are listed in Table 2.

The multiscale modeling framework was implemented using FEniCS (Alnæs et al., 2015). The biventricular unit was meshed with ~7,700 tetrahedral elements based on a previous study (Finsberg et al., 2018) showing that local fiber stress and global features related to cardiac contraction are not sensitive to mesh resolution beyond ~4,000 elements. Furthermore, the aorta and pulmonary arteries FE models contain ~8,000 triangular elements based on previous study (Zeinali-Davarani et al., 2011) that used ~1,500 elements. Steady state PV loop was established by running the simulation over several cardiac cycles until cycle-to-cycle periodicity was achieved. The prescribed cardiac cycle time (690 ms) was derived from the heart rate (87 bpm) measured during the PC-MRI acquisition.

Since it is known that key features of the progression of PAH include stiffening of main PA, reduced RV contraction, and increased distal resistance of PA (Fan et al., 1997; Shimoda and Laurie, 2013), we used our calibrated model to investigate how these changes affect the cardiopulmonary circulation. Specifically, a sensitivity analysis on the parameters associated with PAH progression was performed by simulating the following cases: (1) a 67% increase in PA collagen mass fraction ϕ_PA,c, (2) a 50% decrease in RV contractility T_0,RV, and (3) a 50% increase in the pulmonary arterial resistance R_pa.

Results

Comparison Between Simulated Results and Clinical Measurements

Model predictions of the LV and RV PV loops, volume waveforms, and pressure waveforms in the baseline case matched reasonably well with the clinically measured PAH patient data described in Section Data Acquisition (Figure 3). Good overall fitting was obtained for the volume and pressure in both the LV and RV with the coefficients of determination R² value of 0.901 and 0.903, respectively (Figure 4). Pressure waveforms in the pulmonary and systemic circulations predicted by the model also agree, in general, reasonably well with the measurements, except for the diastolic pressure. The model predicted smaller diastolic pressure in the aorta (by ~17 mmHg) and PA (by ~15 mmHg) when compared to the measurements (Figure 3B). The simulated ascending AO and PA diameter waveforms compared well with the clinical measurements of the dynamic cross-sectional area from the PC-MRI (Figure 3C). Specifically, the simulated and clinically measured diameter waveforms in the ascending AO are in good agreement (max. abs difference ~10%) while the model predicted a larger change of the diameter compared to the measurements for the MPA (max. abs difference ~28%).

FIGURE 3

Figure 3. Measurements and model predictions for the baseline case. (A) LV and RV PV loops; (B) pressure waveforms of pulmonary circulation; (C) pressure waveforms of systemic circulations; (D) LV and RV volume waveforms; (E) MPA and AO diameter waveforms.

FIGURE 4

Figure 4. Scatter plot of the simulated vs. measured volume (left) and pressure (right) at all cardiac time points of the baseline case with a linear fit showing the zero-error reference.

Effects of the Changes in Vascular Microstructure on Cardiac Function

Changing the mass fractions of the constituents in the PA wall led to changes in its function, which in turn affects the RV function. Specifically, increasing the mass fraction ϕ_PA,c of the collagen of PA wall by 67% (from 0.42 to 0.70) with a corresponding decrease in the mass fraction of the elastin (from 0.35 to 0.15) and SMC (from 0.23 to 0.15) produced an increase in the PA pressure of 10% (from 71 to 78 mmHg). The RV systolic pressure also increased by 11% (from 68 to 76 mmHg) correspondingly (Figure 5A). Because of the more exponential mechanical response of the PA with higher collagen fraction, the PA pressure also decayed more rapidly during the diastolic phase resulting in an increased pulse pressure (from 45 mmHg baseline to 55 mmHg) (Figure 5C). The LV and RV SV and EF remained relatively unchanged (Figures 5A,B). In the aorta, systolic, diastolic, and pulse pressures did not change significantly from the baseline case (Figure 5D). The change in PA diameter was slightly reduced when compared to baseline (Figure 5F) as the vessel becomes stiffer with higher collagen mass fraction. Spatially averaged RV fiber stress did not change when compared to the baseline case. Maximum arterial wall stress located at the bifurcation increased (~7.4%) but the spatially averaged wall stress did not change significantly from baseline (Figure 6).

FIGURE 5

Figure 5. Hemodynamic comparison between the different simulations. (A) RV and (B) LV PV loops; (C) Pulmonary and (D) systemic circulation pressure waveforms; (E) LV and RV volume waveforms; and (F) MPA and AO diameter waveforms.

FIGURE 6

Figure 6. Comparison of wall stresses in the different simulations. (A) Von-mises stress map of the pulmonary artery FE model. (B) Average von-mises wall stress waveforms in the pulmonary artery FE model. (C) Average fiber stress waveforms in the biventricular FE model and fiber stress map of baseline model.

Effects of the Change in RV Contractility on Vasculature

Decreasing the RV contractility T_0,RV by 50% (from 1,800 kPa baseline to 900 kPa) reduced the RV EF by 5% (from 58 to 53%) (Figure 5A). Due to less contractile force being generated by the RV, both RV and PA peak systolic pressure decreased by about 9% (RV: 71 to 65 mmHg; PA: 68 to 62 mmHg) (Figure 5C). In addition, the LV EF as well as peak systolic pressure in both the LV and aorta were slightly decreased compared to the baseline (Figures 5B,D,E). Because of the reduced pressure, PA diameter was slightly reduced during systole when compared to baseline (Figure 5F). Average RV fiber stress also decreased by 37% (from 195 to 124 kPa) compared to baseline. Both maximum arterial and spatially averaged RV wall stress were reduced by about 9% (Figure 6).

Effects of the Change in PA Resistance

Increasing the pulmonary arterial resistance R_pa by 50% led to an increase in PA pulse pressure by 36% (from 45 to 61 mmHg), which was also accompanied by an increase in PA systolic and diastolic pressure (Figure 5C). The RV peak systolic pressure increased by 34% (from 71 to 95 mmHg) and the RV EF decreased by 2% (from 58 to 56%) (Figure 5A). Due to the higher pressure, the PA diameter waveform shifted upwards and became higher than the baseline throughout the cardiac cycle. Similar to the case with reduced RV contractility, LV EF as well as peak systolic pressure in both the LV and aorta were slightly decreased compared to the baseline (Figures 5B,D,E). A 7% (195 to 208 kPa) increase in average RV fiber stress as well as a 41% increase in maximum arterial wall stress were also found in the PA (Figure 6).

Discussion

In order to characterize the intricate progression of PAH, we developed the first closed-loop multiscale modeling framework (consisting of image-based FE models of the left and right ventricles, large pulmonary arteries, and aorta) that captures detailed bi-directional ventricular-arterial interactions. We have shown that our proposed model describes the cardiopulmonary circulation reasonably well by reproducing patient-specific measurements of (1) LV and RV PV loops, (2) LV and RV volume and pressure waveforms, and (3) aorta and PA pressure and diameter waveforms of a PAH patient.

This framework extends our previously developed hybrid lumped-FE model of the systemic circulation (Shavik et al., 2018) by including the RV, large pulmonary arteries and the pulmonary micro-circulation (represented with a lumped model). Previous modeling frameworks have coupled a FE biventricular model with a lumped representation of the pulmonary circulation (Kerckhoffs et al., 2007; Xi et al., 2016) but not with FE model of the large pulmonary arteries. The ability to couple a FE model of the large arteries and both ventricles in this framework enables us to investigate PAH progression reflected in the large pulmonary arteries and the RV. Specifically, the framework allows us to alter the microstructural, geometrical and mechanical behaviors of the pulmonary arteries and characterize how these changes affect the RV, and vice versa. Implementing 3D FE models of the arteries in the framework also allow us to capture non-homogeneous stress distribution in the vessels (e.g., high stress concentration at the bifurcation of the pulmonary artery in Figure 6) which would not be possible using lumped-parameter models. Using the calibrated framework, we have created three cases to simulate progressive pathological changes associated with PAH in the (1) large pulmonary arteries (increase in collagen mass and degradation of elastin) (Wang et al., 2013), (2) RV (decrease in contractility due to right ventricular failure) (Naeije and Manes, 2014), and (3) pulmonary microcirculation (increase in resistance due to remodeling) (Kobs et al., 2005).

Increasing the collagen mass in the elastic proximal pulmonary arteries increased PA pulse pressure from baseline. This behavior is due to the stiffening of the PA, which results from a more exponential stress-strain behavior associated with the higher concentration of collagen fibers. This result is consistent with animal experiments where an increase in PA pulse pressure has been associated with an increase in collagen mass in PAH (Wang et al., 2017). Furthermore, the connection between pulse pressure and changes in collagen can also be found in the aorta during aging, where a loss of elastin (which results in a more collagen-dominated extracellular matrix) produces an increase in systemic pulse pressure (Safar et al., 2003). A decrease in PA compliance that is caused by an increase in collagen mass produced an increase in RV afterload as reflected by an increase in RV systolic pressure in our model, consistent with previous studies (Mahapatra et al., 2006; Gan et al., 2007). Consistent with our previous study (Shavik et al., 2018), the more pulsatile PA waveform can also be observed in the ejection phase of the RV PV loop, where the pressure-volume curve became steeper toward end-of-systole (Figure 5A). Our model did not predict a significant reduction in the SV, which could be attributed to a high RV end-systolic elastance in the model. We note that a high RV end-systolic elastance has also been associated with PAH (Vélez-Rendón et al., 2018), especially during the compensatory phase.

Decreasing RV contractility (by 50%) in the model, which reflects the transition to decompensated heart failure, produced an expected decrease in EF and peak systolic pressure that results in a substantial decrease in myofiber stress in the RV. Reducing the RV contractility also reduces the PA peak and pulse pressures, only decreasing the arterial wall stress in the PA slightly. Based on consensus that arterial wall stress is the driver for vascular remodeling (Humphrey, 2008), this result suggests that remodeling in the large pulmonary arteries may attenuate the transition to the decompensated phase. This result also suggests that negative inotropic agent targeted at the RV may help attenuate remodeling in the PA vasculature.

Lastly, increasing the distal pulmonary arterial resistance, which reflects remodeling of the distal vessels, increased pressures in the proximal PA and RV. A 50% increase in the distal pulmonary arterial resistance (equivalent to a ~10% reduction of the vascular lumen diameter based on Poiseuille's law) causes ejection to start at a higher pressure and the EF to be slightly reduced in the RV. These results are broadly consistent with the effects on the RV measured in patients under acute hypoxia (Akgül et al., 2007), which shows an increase in both end-systolic and end-diastolic volume and a slight (but not significant) decrease in EF. The same increase in resistance also produced a significantly higher increase in the systolic PA pressure than the simulation with a 67% increase in collagen mass in the proximal pulmonary arteries. These results suggest that remodeling in the microcirculation contributes more to changes in the pulmonary pressure than remodeling in the proximal pulmonary arteries, suggesting that PAH is primarily driven by distal arterial remodeling. In summary, we have shown that isolated changes in both the arteries and ventricles as predicted by our modeling framework lead to expected effects in the cardiopulmonary circulation. This confirms that the modeling framework can capture bi-directional ventricular-arterial interactions, which can be used to further our understanding of PAH progression.

Model Limitations

Though our modeling framework is able to predict behaviors that are consistent with the measurements there are, however, some limitations associated with it. First, the local myofiber orientation was varied transmurally from 60° in the endocardium to −60° at the epicardium using a “rule based” method. Thus, we did not take into account any changes in myofiber orientation during RV remodeling (Hill et al., 2014) that may occur in PAH. Second, we have assumed a uniform wall thickness and homogeneous material properties for both aorta and PA in our model. We believe that this assumption contributes to the mismatch in the MPA diameter waveforms. Third, we have assumed that FE models of the pulmonary arteries and aorta account for the compliance of the entire pulmonary and systemic arterial system, respectively. This is a limitation because the FE models are associated with only a segment of their corresponding arterial systems. We show in a preliminary study (see Appendix) that the addition of a lumped-parameter compliance to the modeling framework can be used to provide a better match of pulmonary artery pressure and diameter waveforms, as well as the pressure-volume loops. Fourth, we have neglected the dynamic behavior of the fluid and its interaction with the vessel walls and the spatial variation of pressure waveform along the aortic and pulmonary tree and shear stress on the luminal surface of the vessels. We note, however, that the computed shear stress (~Pa) is several order of magnitude smaller than the normal traction force (pressure) on the surface of the vessel (~kPa) and variation of peak pressure within the vessel is <10%. For these reasons, the omission of shear traction should not affect the computed arterial stresses. Last, the modeling framework was calibrated using data acquired from one PAH patient. Caution must be exercised in extrapolating results to the general population of pediatric PAH patients.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author.

Ethics Statement

This study was approved by the University of Michigan Board of Review (HUM00117706). Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.

Author Contributions

SS, SB, and LL developed the theoretical formulation and computational framework of the model. SS and CT-B carried out the simulations for different cases and prepared the results. CT-B and CF acquired the clinical data. LL, SB, and CF planned and supervised the work. All authors helped in interpretation of the results and contributed to the final manuscript.

Funding

This work was supported by American Heart Association (AHA) 17SDG33370110, NIH R01HL134841, and NIH U01HL135842 grants.

Conflict of Interest

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

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2020.00002/full#supplementary-material

Abbreviations

AO, Aorta; EDPVR, End-diastolic Pressure Volume Relationship; EF, Ejection Fraction; ESPVR, End-systolic Pressure Volume Relationship; FE, Finite Element; HR, Heart Rate; LA, Left Atrium; LPA, Left Pulmonary Artery; LV, Left Ventricle; LVEDV, Left Ventricular End-diastolic Volume; LVEF, Left Ventricular Ejection Fraction; LVESV, Left Ventricular End-systolic Volume; LVFW, Left Ventricular Free Wall; MAP, Mean Aortic Pressure; MPA, Main Pulmonary Artery; mPAP, Mean Pulmonary Arterial Pressure; MR, Magnetic Resonance; PA, Pulmonary Artery; PAH, Pulmonary Arterial Hypertension; PC-MRI, Phase-Contrast Magnetic Resonance Image; PCWP, Pulmonary Capillary Wedge Pressure; PV, Pressure-Volume; pv, Pulmonary Veins; PVR, Pulmonary Vascular Resistance; RA, Right Atrium; RPA, Right Pulmonary Artery; RV, Right Ventricle; RVEDV, Right Ventricular End-diastolic Volume; RVEF, Right Ventricular Ejection Fraction; RVESV, Right Ventricular End-systolic Volume; RVFW, Right Ventricular Free Wall; sa, Systemic Arteries; SMC, Smooth Muscle Cells; sv, Systemic Veins; WU, Wood Unit.

References

Akgül, F., Batyraliev, T., Karben, Z., and Pershukov, I. (2007). Effects of acute hypoxia on left and right ventricular contractility in chronic obstructive pulmonary disease. Int. J. Chron. Obstruct. Pulmon. Dis. 2:77. doi: 10.2147/copd.2007.2.1.77

PubMed Abstract | CrossRef Full Text | Google Scholar

Alnæs, M., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., et al. (2015). The FEniCS project version 1.5. Arch. Numer. Softw. 3, 9–23. doi: 10.11588/ans.2015.100.20553

CrossRef Full Text | Google Scholar

Baek, S., Valentín, A., and Humphrey, J. D. (2007). Biochemomechanics of cerebral vasospasm and its resolution: II. Constitutive relations and model simulations. Ann. Biomed. Eng. 35, 1498–1509. doi: 10.1007/s10439-007-9322-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Borlaug, B. A., and Kass, D. A. (2011). Ventricular-vascular interaction in heart failure. Cardiol. Clin. 29, 447–459. doi: 10.1016/j.ccl.2011.06.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, D., Wannenburg, T., and De Tombe, P. P. (1997). Decreased myocyte tension development and calcium responsiveness in rat right ventricular pressure overload. Circulation 95, 2312–2317. doi: 10.1161/01.CIR.95.9.2312

PubMed Abstract | CrossRef Full Text | Google Scholar

Finsberg, H., Xi, C., Tan, J. L., Zhong, L., Genet, M., Sundnes, J., et al. (2018). Efficient estimation of personalized biventricular mechanical function employing gradient-based optimization. Int. J. Numer. Method. Biomed. Eng. 34:e2982. doi: 10.1002/cnm.2982

PubMed Abstract | CrossRef Full Text | Google Scholar

Gan, C. T. J., Lankhaar, J. W., Westerhof, N., Marcus, J. T., Becker, A., Twisk, J. W. R., et al. (2007). Noninvasively assessed pulmonary artery stiffness predicts mortality in pulmonary arterial hypertension. Chest. 132, 1906–1912. doi: 10.1378/chest.07-1246

PubMed Abstract | CrossRef Full Text | Google Scholar

Geuzaine, C., and Remacle, J. F. (2009). Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities. Int. J. Numer. Methods Eng. 79, 1309–1331. doi: 10.1002/nme.2579

CrossRef Full Text | Google Scholar

Guccione, J. M., McCulloch, A. D., and Waldman, L. K. (1991). Passive material properties of intact ventricular myocardium determined from a cylindrical model. J. Biomech. Eng. 113, 42–55. doi: 10.1115/1.2894084

PubMed Abstract | CrossRef Full Text | Google Scholar

Hill, M. R., Simon, M. A., Valdez-Jasso, D., Zhang, W., Champion, H. C., and Sacks, M. S. (2014). Structural and mechanical adaptations of right ventricle free wall myocardium to pressure overload. Ann. Biomed. Eng. 42, 2451–2465. doi: 10.1007/s10439-014-1096-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoit, B. D., Shao, Y., Gabel, M., and Walsh, R. A. (1994). In vivo assessment of left atrial contractile performance in normal and pathological conditions using a time-varying elastance model. Circulation. 89, 1829–1838. doi: 10.1161/01.CIR.89.4.1829

PubMed Abstract | CrossRef Full Text | Google Scholar

Humphrey, J. D. (2008). Mechanisms of arterial remodeling in hypertension coupled roles of wall shear and intramural stress. Hypertension. 52, 195–200. doi: 10.1161/HYPERTENSIONAHA.107.103440

PubMed Abstract | CrossRef Full Text | Google Scholar

Kerckhoffs, R. C. P., Bovendeerd, P. H. M., Kotte, J. C. S., Prinzen, F. W., Smits, K., and Arts, T. (2003). Homogeneity of cardiac contraction despite physiological asynchrony of depolarization: a model study. Ann. Biomed. Eng. 31, 536–547. doi: 10.1114/1.1566447

PubMed Abstract | CrossRef Full Text | Google Scholar

Kerckhoffs, R. C. P., Neal, M. L., Gu, Q., Bassingthwaighte, J. B., Omens, J. H., and McCulloch, A. D. (2007). Coupling of a 3D finite element model of cardiac ventricular mechanics to lumped systems models of the systemic and pulmonic circulation. Ann. Biomed. Eng. 35, 1–18. doi: 10.1007/s10439-006-9212-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Kobs, R. W., Muvarak, N. E., Eickhoff, J. C., and Chesler, N. C. (2005). Linked mechanical and biological aspects of remodeling in mouse pulmonary arteries with hypoxia-induced hypertension. Am. J. Physiol. Heart Circ. Physiol. 288, 1209–1217. doi: 10.1152/ajpheart.01129.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Ky, B., French, B., May Khan, A., Plappert, T., Wang, A., Chirinos, J. A., et al. (2013). Ventricular-arterial coupling, remodeling, and prognosis in chronic heart failure. J. Am. Coll. Cardiol. 62, 1165–1172. doi: 10.1016/j.jacc.2013.03.085

PubMed Abstract | CrossRef Full Text | Google Scholar

Lau, K. D., and Figueroa, C. A. (2015). Simulation of short-term pressure regulation during the tilt test in a coupled 3D−0D closed-loop model of the circulation. Biomech. Model. Mechanobiol. 14, 915–929. doi: 10.1007/s10237-014-0645-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Mahapatra, S., Nishimura, R. A., Sorajja, P., Cha, S., and McGoon, M. D. (2006). Relationship of pulmonary arterial capacitance and mortality in idiopathic pulmonary arterial hypertension. J. Am. Coll. Cardiol. 48, 850–851. doi: 10.1016/j.jacc.2005.09.054

CrossRef Full Text | Google Scholar

Naeije, R., and Manes, A. (2014). The right ventricle in pulmonary arterial hypertension. Eur. Respir. Rev. 23, 476–487. doi: 10.1183/09059180.00007414

PubMed Abstract | CrossRef Full Text | Google Scholar

Pezzuto, S., and Ambrosi, D. (2014). Active contraction of the cardiac ventricle and distortion of the microstructural architecture. Int. J. Numer. Method. Biomed. Eng. 30, 1578–1596. doi: 10.1002/cnm.2690

PubMed Abstract | CrossRef Full Text | Google Scholar

Pezzuto, S., Ambrosi, D., and Quarteroni, A. (2014). An orthotropic active-strain model for the myocardium mechanics and its numerical approximation. Eur. J. Mech. A/Solids. 48, 83–96. doi: 10.1016/j.euromechsol.2014.03.006

CrossRef Full Text | Google Scholar

Punnoose, L., Burkhoff, D., Rich, S., and Horn, E. M. (2012). Right ventricular assist device in end-stage pulmonary arterial hypertension: insights from a computational model of the cardiovascular system. Prog. Cardiovasc. Dis. 55, 234–243.e2. doi: 10.1016/j.pcad.2012.07.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Safar, M. E., Levy, B. I., and Struijker-Boudier, H. (2003). Current perspectives on arterial stiffness and pulse pressure in hypertension and cardiovascular diseases. Circulation. 107, 2864–2869. doi: 10.1161/01.CIR.0000069826.36125.B4

PubMed Abstract | CrossRef Full Text | Google Scholar

Shavik, S. M., Jiang, Z., Baek, S., and Lee, L. C. (2018). High spatial resolution multi-organ finite element modeling of ventricular-arterial coupling. Front. Physiol. 9:119. doi: 10.3389/fphys.2018.00119

PubMed Abstract | CrossRef Full Text | Google Scholar

Shavik, S. M., Wall, S. T., Sundnes, J., Burkhoff, D., and Lee, L. C. (2017). Organ-level validation of a cross-bridge cycling descriptor in a left ventricular finite element model: effects of ventricular loading on myocardial strains. Physiol. Rep. 5:e13392. doi: 10.14814/phy2.13392

PubMed Abstract | CrossRef Full Text | Google Scholar

Shavik, S. M., Zhong, L., Zhao, X., and Lee, L. C. (2019). In-silico assessment of the effects of right ventricular assist device on pulmonary arterial hypertension using an image based biventricular modeling framework. Mech. Res. Commun. 97, 101–111. doi: 10.1016/j.mechrescom.2019.04.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Shimoda, L. A., and Laurie, S. S. (2013). Vascular remodeling in pulmonary hypertension. J. Mol. Med. 91, 297–309. doi: 10.1007/s00109-013-0998-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Simonneau, G., Montani, D., Celermajer, D. S., Denton, C. P., Gatzoulis, M. A., Krowka, M., et al. (2019). Haemodynamic definitions and updated clinical classification of pulmonary hypertension. Eur. Respir. J. 53:1801913. doi: 10.1183/13993003.01913-2018

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, B. W., Chase, J. G., Nokes, R. I., Shaw, G. M., and Wake, G. (2004). Minimal haemodynamic system model including ventricular interaction and valve dynamics. Med. Eng. Phys. 26, 131–139. doi: 10.1016/j.medengphy.2003.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Ursino, M. (1998). Interaction between carotid baroregulation and the pulsating heart: a mathematical model. Am. J. Physiol. 275(5 Pt 2), H1733–H1747. doi: 10.1152/ajpheart.1998.275.5.H1733

PubMed Abstract | CrossRef Full Text | Google Scholar

Vélez-Rendón, D., Zhang, X., Gerringer, J., and Valdez-Jasso, D. (2018). Compensated right ventricular function of the onset of pulmonary hypertension in a rat model depends on chamber remodeling and contractile augmentation. Pulm. Circ. 8:2045894018800439. doi: 10.1177/2045894018800439

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Z., Lakes, R. S., Eickhoff, J. C., and Chesler, N. C. (2013). Effects of collagen deposition on passive and active mechanical properties of large pulmonary arteries in hypoxic pulmonary hypertension. Biomech. Model. Mechanobiol. 12, 1115–1125. doi: 10.1007/s10237-012-0467-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Z., Schreier, D. A., Abid, H., Hacker, T. A., and Chesler, N. C. (2017). Pulmonary vascular collagen content, not cross-linking, contributes to right ventricular pulsatile afterload and overload in early pulmonary hypertension. J. Appl. Physiol. 122, 253–263. doi: 10.1152/japplphysiol.00325.2016

CrossRef Full Text | Google Scholar

Xi, C., Latnie, C., Zhao, X., Tan, J. L., Wall, S. T., Genet, M., et al. (2016). Patient-specific computational analysis of ventricular mechanics in pulmonary arterial hypertension. J. Biomech. Eng. 138:111001. doi: 10.1115/1.4034559

PubMed Abstract | CrossRef Full Text | Google Scholar

Zambrano, B. A., McLean, N. A., Zhao, X., Tan, J.-L., Zhong, L., Figueroa, C. A., et al. (2018). Image-based computational assessment of vascular wall mechanics and hemodynamics in pulmonary arterial hypertension patients. J. Biomech. 68, 84–92. doi: 10.1016/j.jbiomech.2017.12.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeinali-Davarani, S., Sheidaei, A., and Baek, S. (2011). A finite element model of stress-mediated vascular adaptation: application to abdominal aortic aneurysms. Comput. Methods Biomech. Biomed. Eng. 14, 803–817. doi: 10.1080/10255842.2010.495344

PubMed Abstract | CrossRef Full Text | Google Scholar