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

Hemodynamic measurements of PAH patient.

Quantities	Values
HR, bpm	75
LVEDV, ml	72
LVESV, ml	25
LVEF, %	65
MAP, mmHg	68
RVEDV, ml	77
RVESV, ml	30
RVEF, %	61
RVEDV/LVEDV	1.07
mPAP, mmHg	59
PCWP, mmHg	6

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

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

Figure 2

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

In Equation (1), V, V, V, V, V, V, V, and V 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, q, q, q, q, q, q, and q.

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

In Equation (2), R, R, R, and R 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, R, R, and R, respectively. To describe the compliance of the systemic and pulmonary vessels, we used the following PV relationships

where V and V are the resting volumes and C and C 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

where,

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, E, and V0,, respectively. The parameters A and B define the atrium curvilinear end-diastolic pressure volume relationship (EDPVR) and the driving function is defined as

where t 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.

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

where Ω0, is the reference configuration of the biventricular unit, is the displacement field, P and P are, respectively, the Lagrange multipliers that constrain the LV cavity volume V() to a prescribed value V and the RV cavity volume V() to a prescribed value V (Pezzuto and Ambrosi, 2014). We note that V and V are prescribed from the closed-loop circulatory model in Equation (6). The Lagrange multiplier p was used to enforce incompressibility of the tissue (i.e., Jacobian of the deformation gradient tensor J = 1). The vectors 1, and 2, 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), denotes a material point in Ω0, and W 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)

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

In Equation (9), is the first Piola Kirchhoff stress tensor and 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 δ, δp, δP, δP, δ1,, and δ2,, 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 , , that satisfies

for all , δP∈ ℝ, δP∈ ℝ, , The solution of Equation (10) gives the relationship between P, P, V, V 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 in Equation (9) was additively decomposed into a passive and an active component, i.e.,

In Equation (11), is the passive stress tensor, P is the magnitude of the active stress, whereas and are the local basis vectors that define the cardiac muscle fiber directions in the current and reference configuration, respectively. The passive stress tensor is related to the strain energy function W and deformation gradient tensor by

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

with

was prescribed. In Equation (13b), E with (i, j) ϵ (f, s, n) denote the components of the Green-Lagrange strain tensor with f, s, n denoting the myofiber, sheet and sheet normal directions, respectively. Material parameters of the Fung-type constitutive model are C, b, b, and b.

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 was described by

where l is the sarcomere length, l is the length of the contractile element, l is the sarcomere length in a prescribed reference state (relaxed sarcomere length), and E is the stiffness of the serial elastic element. The function denotes the dependency of the isometrically developed active stress on l and is given by

where T0 is a model parameter that scales the active tension. Both a6 and a7 are model parameters. The time course of the active tension development is controlled by

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

where v0 is the unloaded shortening velocity. The sarcomere length l was calculated from the myofiber stretch λ and the relaxed sarcomere length l by

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

where Ω0, is the reference configuration of the arteries, is the displacement field and P is the Lagrange multiplier that constrains the arterial cavity volume V() to a prescribed value V. The vectors 1, and 2, 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 that satisfies

for all , The solution above gives the relationships between P, V, and P, V 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 in Equation (21), which is given as the sum of the key tissue constituents, namely, elastin-dominated matrix W, collagen fiber families W and vascular smooth muscle cells (SMC) W (Baek et al., 2007; Zeinali-Davarani et al., 2011), i.e.,

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

where M is the mass per unit volume of the elastin in the tissue, C is a stiffness parameter and, 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

In Equation (23), M is the mass per unit volume of ith family of collagen fibers, λ is the corresponding stretch of those fibers, and both c and c are the material parameters that govern the collagen stiffness. The stretch in the ith family of collagen fibers was defined by where is the local unit vector that defines the corresponding fiber orientation.

Strain energy function of the SMC W is given by

Here, M is the mass per unit volume of the SMC in the tissue, λ is the stretch of the SMC, whereas c and c 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

where ϕ, ϕ, ϕ 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

Model parameters for FE models for the baseline case.

Biventricular FE model
Passive material model	CLV = 280 Pa, CRV = 170 Pa
Active contraction model	T0,LV = 2000 kPa, T0,RV = 1800 kPa, tr= 280 ms, td = 80 ms, b = 0.17 ms.μm−1
Circulatory model	Csv = 0.02 Pa•ml, Cpv = 0.09 Pa•ml, Rsa = 125 kPa•ms•ml−1, Rpa = 75 kPa•ms•ml−1, Rsv = Rpv = 2 kPa•ms•ml−1, Rav = 3.2 kPa•ms•ml−1, Rmv = 0.9 kPa•ms•ml−1, Rtv = 0.4 kPa•ms•ml−1, Rpvv = 2 kPa•ms•ml−1, Vsv,0 = 3570 ml, Vpv,0 = 485 ml
Time varying elastance model of LA and RA	Ees = 60 Pa/ml, V0 = 10 ml, tmax = 135 ms, τ = 50 ms, A = 58.67 Pa, B = 0.049 ml−1
Aorta FE model
Elastin	cAO,1 = 120 kPa, ϕAO,e = 0.35
Collagen families	cAO,2 = 0.2 kPa, cAO,3 = 8.0, ϕAO,c = 0.20 (ϕAO,1 = 0.1ϕAO,c, ϕAO,2 = 0.1ϕAO,c, ϕAO,3 = 0.4ϕAO,c, ϕAO,4 = 0.4ϕAO,c)
SMC	cAO,4 = 0.08 kPa, cAO,5 = 3.5, ϕAO,m = 0.45
Pulmonary artery FE model
Elastin	cPA,1 = 45 kPa, ϕPA,e = 0.35
Collagen families	cPA,2 = 100.0 kPa, cPA,3 = 3.0, ϕPA,c = 0.42 (ϕPA,1 = 0.1ϕPA,c, ϕPA,2 = 0.1ϕPA,c, ϕPA,3 = 0.4ϕPA,c, ϕPA,4 = 0.4ϕPA,c)
SMC	cPA,4 = 5 kPa, cPA,5 = 3.5, ϕPA,m = 0.23

Solution Algorithm

An explicit time integration scheme was used to solve the ODEs in Equation (1). Specifically, compartment volumes (V, V, V, V, V, V, V, V) at each time t were determined from their respective values and the segmental flow rates (q, q, q, q, q, q, q, q) at previous time t in Equation (2). The computed compartment volumes at t were used to update the corresponding pressures (P, P, P, P, P, P, P, P). Pressures in the atrium (P, P) and veins (P, P) were computed from Equations (4) and (3), respectively. On the other hand, pressures in the LV (P), RV (P), were computed from the FE solutions of Equation (10) for the biventricular unit with the volumes (V, V) at time t as input. Similarly, pressures in the aorta (P) and pulmonary artery (P) were computed from the FE solutions of Equation (20) with their corresponding volumes (V, V) at time t. We note here that (P, P, P, P) are scalar Lagrange multipliers in the FE formulation for constraining the cavity volumes to the prescribed values (V, V, V, V). The computed pressures at time t were then used to update the segmental flow rates in Equation (2) that will be used to compute the compartment volumes at time t 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 T0 were prescribed to be the same values in the LVFW and septum (denoted as C and T0,) and had different values in the RVFW (denoted as C and T0,). 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 and C. Stroke volume (SV) of the LV and RV were matched by adjusting the regional contractility parameters (i.e., T0,, T0,). While other model parameters can also affect the SV (e.g., peripheral resistances R and R of the systemic and pulmonary circulations as well as preload), the parameters T0, and T0,, 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, t and b were adjusted to match the time course of the volume and pressure waveforms measured in the LV and RV. Parameters t and t 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, R) 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 ϕ, (2) a 50% decrease in RV contractility T0,, and (3) a 50% increase in the pulmonary arterial resistance R.

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 R2 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

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

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 ϕ 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

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

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

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.