Biomech Model Mechanobiol DOI 10.1007/s10237-017-0940-4

ORIGINAL PAPER

Cardiovascular mechanics in the early stages of pulmonary hypertension: a computational study Sebastián Acosta1 · Charles Puelz2 · Béatrice Rivière2 · Daniel J. Penny1 · Ken M. Brady3 · Craig G. Rusin1

Received: 5 January 2017 / Accepted: 12 July 2017 © Springer-Verlag GmbH Germany 2017

Abstract We formulate and study a new mathematical model of pulmonary hypertension. Based on principles of fluid and elastic dynamics, we introduce a model that quantifies the stiffening of pulmonary vasculature (arteries and arterioles) to reproduce the hemodynamics of the pulmonary system, including physiologically consistent dependence between compliance and resistance. This pulmonary model is embedded in a closed-loop network of the major vessels in the body, approximated as one-dimensional elastic tubes, and zero-dimensional models for the heart and other organs. Increasingly severe pulmonary hypertension is modeled in the context of two extreme scenarios: (1) no cardiac compensation and (2) compensation to achieve constant cardiac output. Simulations from the computational model are used to estimate cardiac workload, as well as pressure and flow traces at several locations. We also quantify the sensitivity of several diagnostic indicators to the progression of pulmonary arterial stiffening. Simulation results indicate that pulmonary pulse pressure, pulmonary vascular compliance, pulmonary RC time, luminal distensibility of the pulmonary The second author was supported by a fellowship from the Keck Center of the Gulf Coast Consortia, on the Training Program in Biomedical Informatics, US National Library of Medicine T15LM007093.

B

Sebastián Acosta [email protected] Charles Puelz [email protected]

1

Department of Pediatrics–Cardiology, Baylor College of Medicine and Texas Children’s Hospital, Houston, TX, USA

2

Department of Computational and Applied Mathematics, Rice University, Houston, TX, USA

3

Department of Anesthesiology, Baylor College of Medicine and Texas Children’s Hospital, Houston, TX, USA

artery, and pulmonary vascular impedance are much better suited to detect the early stages of pulmonary hypertension than mean pulmonary arterial pressure and pulmonary vascular resistance, which are conventionally employed as diagnostic indicators for this disease. Keywords Pulmonary arterial hypertension · Computational hemodynamics · Cardiovascular mechanics · Blood flow

1 Introduction Pulmonary hypertension is a vascular disease that affects 1–2 million people each year in the USA and is characterized by a progressive increase in pulmonary pressure which eventually leads to right heart failure (Gaine and Rubin 1998; Farber and Loscalzo 2004). Pulmonary pressure increases due to stiffening of the pulmonary vessels, which causes pulmonary vascular resistance to increase, forcing the heart to work harder to maintain its cardiac output. Left untreated, it is fatal. Early detection of pulmonary hypertension has been linked to improved outcomes (Humbert et al. 2010a, b; Thenappan et al. 2010). However, the typical delay in the diagnosis of this disease has been estimated to be at least 2 years from the appearance of early symptoms (Gibbs 2007; Humbert et al. 2012). Considering that the historical median survival after diagnosis is only 2.8 years (Gibbs 2007; Rich and Rich 2014; Thenappan et al. 2010; Humbert et al. 2010b), a timely diagnosis is critical for improving outcomes. Currently, confirmation and assessment of the severity of the disease often requires cardiac catheterization. This invasive procedure involves threading a catheter from the vena cava through the right atrium and ventricle of the heart to the pulmonary artery, where pressure measurements can be

123

S. Acosta et al.

made. The development of a sensitive and specific noninvasive physiologic indicator of the early stages of pulmonary hypertension would have a significant positive impact on this population. The focus of this work is the development of a highfidelity mathematical model of the hemodynamics associated with pulmonary hypertension. The model is interrogated to establish hemodynamic parameters most sensitive to initial physiologic changes caused by the stiffening of the pulmonary vasculature. This in turn provides physiologic targets for future clinical studies. Our work is complementary to other studies regarding the modeling and detection of pulmonary hypertension. An excellent overview, with an extensive reference list, concerning pulmonary vascular stiffness and its value as a diagnostic marker was published by Hunter et al. (2011). A review of the computational tools to simulate pulmonary blood flow is offered by Kheyfets et al. (2013). These authors provide a clinical overview of pulmonary hypertension and discuss how numerical models can be used to gain insight into cardiopulmonary pathology. Previous work on mathematical modeling of pulmonary hypertension can be found in the literature (Lankhaar et al. 2006; Punnoose et al. 2012; Rausch et al. 2011; Tang et al. 2012; Qureshi et al. 2014). Some studies utilize 0D Windkessel models in which RC circuit elements are used to represent different parts of the pulmonary circulatory system. Lankhaar et al. (2006) used such a model to study different components contributing to right ventricular afterload, including pulmonary vascular resistance, compliance, and impedance. Their conclusion highlighted the clinical importance of vascular compliance in detecting the early stages of pulmonary hypertension. Punnoose et al. (2012) used a similar approach to examine the effect of a right ventricular assist device in patients with pulmonary hypertension. We also mention the work of Segers et al. (2000) for arterial hypertension. The authors employed a 0D model for the arterial circulation and the heart to elucidate the impact of different parameters on the arterial pressure waveform. Included in their model were different types of heart remodeling (Segers et al. 2000). While 0D equations employed by these researchers are computationally inexpensive and suitable for approximating the entire circulation, they yield models with less anatomical detail which often cannot accurately describe pressure wave propagation and physiologic dynamics occurring during the transition from a normotensive to a hypertensive state (Acosta et al. 2015a). Numerical studies of pulmonary hypertension have also been made with higher-fidelity models. Qureshi et al. (2014) developed a multiscale model with one-dimensional representation of large vessels and hierarchical bifurcating trees for the pulmonary vascular beds. This model allows them to resolve inner patterns of vascular remodeling, pressure drop, and pulse wave propagation within the pulmonary

123

bed. Rausch et al. (2011) developed a continuum mechanics approach for mathematically describing the remodeling process of the heart and applied this technique to analyze the effect of elevated afterload in both ventricles. The work of Lungu et al. (2014) uses 0D and 1D models to derive pulmonary hypertension diagnostic parameters based on noninvasive MRI imaging data. The 0D model shows an increase in pulmonary resistance, and the 1D model shows that there is a corresponding increase in the power of the reflected pressure wave as the disease worsens. The same group has recently published work to improve diagnostic accuracy of pulmonary hypertension using machine learning techniques and markers derived from computational models and MRI data (Lungu et al. 2016). Work of Tang et al. (2012) involved 3D fluid dynamics simulations of the pulmonary arteries impacted by hypertension to determine changes in wall shear stress. A related paper from Yang et al. (2016) also modeling the pulmonary circulation, used 3D simulations to study peripheral pulmonary artery stenosis. This work investigated the implementation of boundary conditions and their impact on the accuracy of the computed flow in the pulmonary vasculature. While 3D models provide greater physiologic detail compared to 1D or 0D models, they require significant computational resources to execute, which generally restricts such models to local simulations of the heart and nearby arteries. Our modeling approach seeks to strike a balance between physiologic accuracy and computational efficiency. The model of pulmonary hypertension is embedded in the framework of 1D equations for the major vessels and 0D equations for the organ beds and heart, allowing for the construction of a whole-body description of circulatory hemodynamics. Our work has been influenced primarily by Hughes and Lubliner ˇ c and Kim (2003), Stergiopulos et al. (1992), (1973), Cani´ Formaggia et al. (2002), Alastruey et al. (2008), Formaggia et al. (2006), Mynard (2011), and Acosta et al. (2015b). Similar 0D-1D closed-loop cardiovascular simulations have been recently developed by Blanco and Feijóo (2013), Mynard and Smolich (2015), Audebert et al. (2017), and Strocchi et al. (2017) among others. Reduced 1D vessel models used in this work may be derived from the Navier–Stokes equations ˇ c and Kim 2003) and have been validated against exper(Cani´ imental data (Olufsen et al. 2000; Reymond et al. 2009). The advantage of this approach is its ability to accurately capture the spatial dependence and nonlinearity of fluid flow while being inexpensive to approximate, as compared to three-dimensional simulations employing the full Navier– Stokes equations (Xiao et al. 2014). As a result, this approach allows for the construction of large, closed-loop vessel networks representing each part of the circulatory system. We adopted the parameters for modeling the network, heart, and organ beds from the works of Mynard (2011) and Mynard and Smolich (2015). The major novelty of our work is the mod-

Cardiovascular mechanics in the early stages of pulmonary hypertension: a computational study

pulmonary veins

left heart

lungs

right pulmonary arteries heart

systemic arteries

organs

systemic veins

Fig. 1 A sketch of the closed-loop model used in our simulations along with the systemic arteries and systemic veins represented as onedimensional vessels. The boxes correspond to the heart, while the circles are “lumped” representations of the systemic organ beds and the lungs. The model also contains one-dimensional representations of the main pulmonary artery with its major right and left branches, the upper and

lower pulmonary veins from each lung, and the portal venous system. The left and right sides of the heart and the organ beds are represented with zero-dimensional or lumped-parameter models. The connectivity of this network illustrated in these figures is taken from the closed-loop model by Mynard (2011)

eling approach for the organ hemodynamics including the lungs. This effort seeks to endow the 0D organ model with the nonlinear dependence of compliance and resistance on blood pressure/volume and vascular composition. Based on Acosta et al. (2015a), the model takes into account the multiscale nature of capillary networks which is mathematically aggregated into effective 0D nonlinear governing equations. The paper is structured as follows. In Sect. 2, we describe our modeling approach and discretization schemes. We end this section with our method for modeling the progression of pulmonary hypertension. Section 3 details the results from our model, and in Sect. 4, we state the limitations of our work. Conclusions are provided in Sect. 5.

of Fig. 1 depicts vessel networks for the systemic arteries in red and the systemic veins in blue, with lines in these figures corresponding to individual 1D vessels. The following sections describe each component of the model along with the associated numerical approximation scheme for the 1D vessels. The systems of differential/algebraic equations modeling the organ beds and heart are approximated by implicit or semi-implicit time discretizations to avoid any time step restrictions. Details regarding the implementation are given at the end of this section.

2 Methods The closed-loop hemodynamic model used in this work contains 158 vessels arranged in networks corresponding to the systemic arterial, systemic venous, pulmonary arterial, pulmonary venous, and portal venous circulations. The connectivity of these networks is taken from the closed-loop model given by Mynard (2011). Each vessel is modeled as a 1D elastic tube. Networks are coupled to 0D models for the left and right sides of the heart and organ beds. A schematic of the closed-loop model is provided in the left panel of Fig. 1. Lines in this panel are vessel networks, and nodes are 0D models of the heart and organ beds. The right panel

2.1 One-dimensional vessel model Let t define the time variable in seconds and x define the space variable in cm along the axial dimension of a cylindrical vessel. To derive equations for one-dimensional blood flow through an cylindrical elastic vessel, one begins with the three-dimensional, axially symmetric Navier–Stokes equations modeling an incompressible Newtonian fluid. It is assumed the characteristic vessel diameter D˜ is much ˜ From an smaller than the characteristic vessel length L. asymptotic reduction argument, by dropping second-order ˜ L, ˜ one obtains the following terms containing the ratio D/ reduced equations for vessel cross-sectional area A and fluid velocity U : ∂ A ∂(AU ) + = 0, ∂t ∂x

(1)

123

S. Acosta et al.

1 ∂U 2 1 ∂p α U ∂U + + = −2π ν . ∂t 2 ∂x ρ ∂x α−1 A

(2)

The variable p is the fluid pressure. The parameter ν = μ/ρ is the kinematic viscosity, defined as the ratio of the dynamic viscosity μ = 3.20 × 10−2 g/cm/s and the blood density ρ = 1.06 g/cm3 . In consistency with the tables in “Appendix,” units of pressure are dynes/cm2 , of velocity are cm/s, and of cross-sectional area are cm2 . However, the units of pressure are converted to mmHg for displaying results from the model in accordance with medical convention. For the model employed in this paper, the choice of an axial velocity profile with unidirectional flow which satisfies the no-slip boundary condition determines the parameter α > 1. In this work, we assume α = 1.1 and remark that the above model also assumes a flat profile in the convective term corresponding to a uniform velocity profile equal to the averaged axial velocity. For details of the model derivation, ˇ c and Kim (2003). see, e.g., the work of Cani´ For some of the models below we refer to the flow Q, in units of cm3 /s, defined as Q = AU.

(3)

The state equation provides a direct relationship between area A and pressure p and is specified as follows: p = p0 + ψ(A; A0 ).

(4)

The form for ψ we use in this work models the vessel wall as a linearly elastic membrane (Mikelic et al. 2007): 

1/2

ψ(A; A0 ) = β(D0 ) A1/2 − A0



,

(5)

where A0 = π(D0 /2)2 . For each network of vessels in the circulation model, we specify the parameter p0 . The parameter A0 is the vessel reference cross-sectional area and determines the vessel stiffness β in units dynes/cm3 via a heuristic scaling law. More precisely, let h (cm) and E (dynes/cm2 ) be the thickness and the circumferential Young’s modulus of the vessel wall. The parameter β may be expressed in the following way (Sherwin et al. 2003): β=

16 1 Eh . √ 3 π D0 D0

(6)

We employ the following scaling law as in Olufsen (1999): Eh = k1 exp (k2 D0 ) + k3 , D0

123

(7)

for some parameters k1 , k3 > 0 and k2 < 0. In our simulations, β is defined by: β(D0 ) =

16 1 √ (k1 exp (k2 D0 ) + k3 ) . 3 π D0

(8)

In this work, p0 is independent of x and t, but we allow A0 to vary with x to model vessel tapering. In the case that A0 depends on x, we evaluate (8) at its average over the vessel length. The vector form for (1)–(2) reads:  

   AU ∂ A ∂ 0 2 + (9) ψ = α U . −2π ν α−1 ∂t U ∂ x U2 + ρ A      :=U

:=F(U)

:=S(U)

The Riemann invariants of (9), denoted V1 and V2 , are functions with gradients parallel to the left eigenvectors of the Jacobian of F. With this definition, they satisfy the following system, ∂ V1 u dc0 ∂ V1 + λ1 = S1 = −8π ν − 4(u + c − c0 ) , ∂t ∂x A dx ∂ V2 u dc0 ∂ V2 + λ2 = S2 = −8π ν + 4(u − c + c0 ) , ∂t ∂x A dx (10)  1/2 where c = c(A) := ρA dψ and c0 = c(A0 ). The symdA bols λ1 and λ2 are the eigenvalues of the Jacobian of F. In the case of (9), one has λ1 = U + c(A), V1 = U + 4 (c(A) − c0 ) ,

(11)

λ2 = U − c(A), V2 = U − 4 (c(A) − c0 ) .

(12)

For physiologic selections of the parameters, the eigenvalues satisfy λ1 > 0 and λ2 < 0. The simplicity of (10) allows us to use a characteristics-based numerical method to simulate the flow of blood (Acosta et al. 2015b). 2.2 Numerical approximation for the one-dimensional vessels We choose to simulate blood flow in the one-dimensional vessels by approximating solutions to (10) using the numerical method of characteristics (NMC). In this section, we briefly review this approach based on the presentation in Acosta et al. (2015b). Similar approaches for these equations have been used in the past (Streeter et al. 1963; Melicher and Gajdošík 2008; Wang and Parker 2004). This method relies on the fact that the functions t Si (x, s)ds, i = 1, 2 (13) Vi − 0

Cardiovascular mechanics in the early stages of pulmonary hypertension: a computational study

are constant on the characteristics γi satisfying dγi = λi (γi (s), s), i = 1, 2. ds

(14)

To concisely describe NMC, assume that Si = 0 (i = 1, 2) in (10) and (13). Details on how to handle nonvanishing source terms Si are found in Acosta et al. (2015b). Let V1 and V2 be the exact solutions to (10) on the interval [a, b] and W1 and W2 to be the corresponding approximate solutions computed by NMC. Let h be the operator that projects a function into its continuous piecewise linear interpolant on a uniform discretization of the interval [a, b]. Given some final time T and time step t, we wish to evolve the solutions for N = T / t time steps. In accordance with the notation in Acosta et al. (2015b), let g˜in (x) be an approximation to the point on the characteristic curve γi at time n t which passes through point x at time (n + 1) t. In this scheme, the characteristic curve for each invariant is linearized at every time step. If the tail of the characteristic curve, g˜in (x), falls in between two points on the discrete spatial grid, piecewise linear interpolation is used to recover the value of the invariant. This process is encoded in the operator h . Pseudocode for NMC is displayed in Algorithm 1. Precise implementation details and numerical analysis of the algorithm may be found in Acosta et al. (2015b). We remark that the method is unconditionally stable and fully explicit in time.

1: 2: 3: 4: 5: 6:

Input initial conditions: V10 , V20 . Initialize numerical solution: W10 = h [V10 ] and W20 = h [V20 ]. for n = 1, 2, . . . N Compute g˜in−1 (x) (i = 1, 2) as in Acosta et al. (2015b). Set Win (x) = h [Win−1 (g˜in−1 (x))]. (i = 1, 2) end

Nout vessels

Ak , Uk

A˜j , U˜j

Fig. 2 A sketch of a single junction within a vessel network. The axial coordinate axis x for each vessel is oriented from the vessel inlet to the vessel outlet. Then, we define the incoming vessels as those containing the junction at their outlets and the outgoing vessels as those containing the junction at their inlets. In the definition of the vessel network, each junction of vessels is determined by specifying the corresponding incoming and outgoing vessels

and V2 . System (15)–(19) is solved at every time step using Newton’s method. 

for k = 1, . . . , Nin V1,k = Uk + 4 c (Ak ) − c A0,k

V2, j

(15)      for j = 1, . . . , Nout = U˜ j − 4 c A˜ j − c A˜ 0, j (16)

Nin 

Algorithm 1 NMC pseudocode

Nin vessels

k=1

A k Uk =

Nout 

A˜ j U˜ j

(17)

j=1

ρ ρ (U1 )2 + p (A1 ) = (Uk )2 + p (Ak ) for k = 2, . . . , Nin 2 2 (18)  2   ρ ρ U˜ j + p A˜ j (U1 )2 + p (A1 ) = 2 2 for j = 1, . . . , Nout . (19)

2.3 Vessel junctions 2.4 Organ bed models Suppose a junction of vessels has Nin incoming vessels and Nout outgoing vessels. In this case, we must determine both A and U at the outlets of the incoming vessels and the inlets of the outgoing vessels yielding 2Nin + 2Nout unknowns. Explicitly, let Ak , Uk , V1,k correspond to the values of these variables at the outlet of the incoming vessels and A˜ j , U˜ j , V2, j as values at the inlet of the outgoing vessels. See Fig. 2 for a sketch of a vessel junction. Junctions of vessels are modeled by enforcing the physical principles of conservation of mass and continuity of total pressure in Eqs. (17) and (18)–(19), respectively. The additional Eqs. (15) and (16) require the differential equation governing blood flow (10) to be satisfied at the junction by employing the analytical form for the Riemann invariants V1

Organ beds are modeled using an electric circuit analogy with resistor and capacitor elements. We employ the model used in Mynard (2011) and Mynard and Smolich (2015) but allow the compliance and resistance to depend on the blood pressure, as described in Acosta et al. (2015a). We provide a brief description as follows. We are given reference values for the resistance Rcap,ref , arteriole compliance Cart,ref , venule compliance Cven,ref , and total volume of blood contained in the organ Vref . Then, we define the reference volume of blood contained in the arteriole and venule compartments as Vart,ref = Cart,ref /(Cart,ref + Cven,ref )Vref and Vven,ref = Cven,ref /(Cart,ref + Cven,ref )Vref , respectively, so that Vart,ref + Vven,ref = Vref .

123

S. Acosta et al.

1D arteries

Rart,k

Qart,k , part,k

Qart , part Cart

Rven,j

Rcap

Qven , pven

Qcap

1D veins Qven,j , pven,j

Cven

Fig. 3 A sketch of the organ bed model. The organ bed resistance is Rcap , the compliance on the arterial side is Cart , and the compliance on the venous side is Cven . The characteristic impedances of the incoming and outgoing vessels are given by Rart,i and Rven, j , respectively. Q art and Q ven are defined in (32)–(33)



2



2 

From the approach in Acosta et al. (2015a), we can derive a state equation for the volume of blood contained in the arteriole and venule compartments of the organ bed,

Rcap = Rcap,ref

(1 + Γart ( part − pext ))2 Vart = Vart,ref 2 , 1 + Γart ( part,ref − pext )

where the volumes Vart and Vven depend on the blood pressure as described by (20)–(21). Lastly, the resistances Rart,k = ρc0,k /A0,k and Rven, j = ρ c˜0, j / A˜ 0, j model the arterioles and venules, respectively (Mynard 2011). We remark that unlike Rcap , Cart , and Cven , these parameters do not depend on the blood pressure. However, they do depend on the vascular properties embedded in c0 and A0 . The equations governing the organ bed model are given below. Equations (27)–(29) formulate the analog of Ohm’s law, while Eqs. (30) and (31) couple flow into the capacitors with conservation of mass. The last two Eqs. (32) and (33) describe mass conservation at the incoming and outgoing vessels. We remark that the model for the liver takes a slightly different form from the circuit given in Fig. 3, since it incorporates a high-pressure hepatic artery and a low-pressure portal vein. We use the same liver model as that given (Mynard 2011) and omit details for brevity.

Vven = Vven,ref

(1 + Γven ( pven − pext ))2 2 , 1 + Γven ( pven,ref − pext )

(20) (21)

where part and pven are the arteriole and venule blood pressure and part,ref and pven,ref are their respective reference values. Also pext denotes the external pressure which is assumed to be vanishing in this paper. The effective arteriole and venule wall distensibility are given by Cart,ref , 2Vart,ref − Cart,ref ( part,ref − pext ) Cven,ref . = 2Vven,ref − Cven,ref ( pven,ref − pext )

Γart =

(22)

Γven

(23)

These are inversely proportional to the stiffness (elastic modulus) of the respective vascular tissue. Hence, they are parameters useful for representing the vascular stiffening exhibited in the progression of pulmonary hypertension and other vascular diseases. From (20)–(21), we obtain the effective compliance associated with the organ bed and their dependence on the blood pressure, Cart =

Cven

dVart (1 + Γart ( part − pext )) = 2Γart Vart,ref 2 , d part 1 + Γart ( part,ref − pext )

dVven = d pven

(24) (1 + Γven ( pven − pext )) = 2Γven Vven,ref 2 . 1 + Γven ( pven,ref − pext ) (25)

Vart,ref Vart

+

Vven,ref Vven

,

(26)

part,k − part = Q art,k Rart,k k = 1, . . . , Nin

(27)

pven − pven, j = Q ven, j Rven, j

(28)

j = 1, . . . , Nout

part − pven = Q cap Rcap d part = Q art − Q cap Cart dt d pven = Q cap − Q ven Cven dt Nin  Q art = Q art,k

(29) (30) (31) (32)

k=1

Q ven =

Nout 

Q ven, j .

(33)

j=1

2.5 Heart model for a closed-loop circulation The dependence of the organ resistance Rcap on the blood pressure was derived in Acosta et al. (2015a). For the configuration depicted in Fig. 3, we have that

123

The equations governing our heart model are given as follows:

Cardiovascular mechanics in the early stages of pulmonary hypertension: a computational study

Ain,k , Qin,k

ventricle

atrium pa , Qa , Va

Aout , Qout

pv , Qv , Vv

outgoing vessel incoming vessels

Qext

Qint

Fig. 4 A sketch of the heart model. Each chamber contains three variables: pressure p, momentum Q, and volume V . The diodes in the sketch are the heart valves, with Q int equal to flow through the valve between the atrium and the ventricle, and Q ext equal to flow leaving the ventricle

 Va − Va0 dVa pa = E a (t) + K a pa Va0 dt 

0 Vv − Vv dVv + K v pv pv = E v (t) 0 Vv dt dQ int + (Bint |Q int | + Rint ) Q int L int dt = H ( pa , pv ; γint ) ( pa − pv ) dQ ext + (Bext |Q ext | + Rext ) Q ext L ext dt = H ( pv , p(Aout ); γext ) ( pv − p(Aout )) dVa = −Q a dt dVv = −Q v dt Nin  Q a = Q int − Q in,k

(34)

where T = heart period,

(35)

(45)

m 1 , τ1 = parameters governing shape of contraction, (46) m 2 , τ2 = parameters governing shape of relaxation,

(36)

(37) (38) (39) (40)

k=1

Q v = Q ext − Q int

(41)

Q out = Q ext .

(42)

The first two Eqs. (34) and (35) (for the atrium and ventricle, respectively) use a specified time-periodic elastance function E(t) to relate chamber volume and chamber pressure. Each chamber contains three variables: pressure p, volume V , and flow Q. The residual chamber volume is given as V 0 . Our model also incorporates a viscoelastic term with parameter K similar to the work of Blanco and Feijóo (2013). A sketch of the model for a single side of the heart is depicted in Fig. 4. We use a simplified form for the elastance function from Mynard and Smolich (2015). The time-dependent profile for each chamber’s elastance is defined as   g1 (t) 1 (E max − E min ) + E min , E(t) = k 1 + g1 (t) 1 + g2 (t) (43)

(47)

tonset = time at which contraction begins,

(48)

E min , E max = minimum and maximum elastances.

(49)

The heart rate is defined as HR = 1/T . We remark that this model neglects chamber interactions. Similar to Blanco and Feijóo (2013) and Mynard (2011), the heart valves are modeled as non-ideal diodes with flow governed by a modified version of Bernoulli’s equation in (36) and (37). The variables and parameters for the interior valve, between the chambers, and the exterior valve, between the ventricle and outgoing vessel, are denoted with the subscript “int” and “ext,” respectively. We choose to multiply the pressure gradient by a smoothed version of the heaviside function depending on a parameter γ > 0: −1  . H ( p1 , p2 ; γ ) := 1 + e−γ ( p1 − p2 )

(50)

This parameter γ controls the contribution of the pressure gradient term in Bernoulli’s equation and in turn governs the shape of the dicrotic notch in the pressure waveforms. Lastly, L describes the valve inductance, R determines the linear resistance through the valve, and B scales the nonlinear resistance term depending on the magnitude of the flow. Equations (38) and (39) relate the time derivative of chamber volume to the fluid momentum exiting the chamber, and the last three Eqs. (40)–(42) enforce conservation of mass. Finally, our heart model contains multiple incoming vessels but only one outgoing vessel. The incoming vessels allow us to retain the systemic and pulmonary venous circulations and close the loop of the model. 2.6 Numerical simulation of closed-loop vessel networks

where

 t − tonset m i gi (t) = , i = 1, 2, and τi   g1 1 , k = max t∈[0,T ] 1 + g1 1 + g2 

(44)

The zero-dimensional systems presented in Sects. 2.3, 2.4, and 2.5 are coupled with the model for one-dimensional blood flow, namely Eq. (9), to construct a closed-loop vessel network. At each time step in the numerical scheme, the variables A and U at the inlets of the outgoing vessels and the

123

S. Acosta et al.

outlets of the incoming vessels are computed by Newton’s method from the systems of equations for the junctions, heart, and organ beds. Boundary conditions for one-dimensional vessels are determined by the Riemann invariants from the computed values of A and U , via Eqs. (11)–(12). In the numerical method of characteristics, the right-to-left moving invariant V2 is prescribed at the outlet of the incoming vessel, and the left-to-right moving invariant V1 is prescribed at the inlet of the outgoing vessel, thereby coupling the one-dimensional and zero-dimensional models. The ordinary differential equations in the organ bed and heart models are discretized semi-implicitly in the sense that nonlinearities in the equations are time-lagged (evaluated at the nth time step). For the organ bed model, we employ backward Euler for the equation involving the capacitance, i.e., for the arterial side: n Cart

 1  n+1 n n+1 part − part = Q n+1 art − Q cap , t

and similarly for the venous side. For example, the discretization for the interior valve reads: L int

  1  n+1 Q int − Q nint + Bint |Q nint | + Rint Q n+1 int t  

= H ( pan , pvn ; γint ) pan+1 − pvn+1 .

Equations for conservation of mass (38) are discretized with Crank–Nicolson. For the atrium, combining Eqs. (34) and (38), one obtains:  pan+1

= E a (t

n+1

)

Van −

t 2

  Q an+1 + Q an − Va0 Va0

− K a pan+1 Q an+1 .

(51)

The equations for the ventricle are discretized in the same way. Software for simulating these models was developed in C/C ++ . The implementation is custom and uses linear solvers from Lapack. 2.7 Model of pulmonary vascular stiffening We assume that pulmonary hypertension is caused by the stiffening of the cardio-arterial tissues in the pulmonary circulation (Wang and Chesler 2013, 2011; Wang et al. 2013; Tan et al. 2014; Acosta et al. 2015a; Vonk-Noordegraaf et al. 2013; Rain et al. 2013; Bogaard et al. 2009; Dragu et al. 2015; Hill et al. 2014). This stiffening affects three separate components in our cardiovascular circulatory model: the major pulmonary arteries, arteriole branching of the lungs, and the right ventricle.

123

Pulmonary vascular stiffening affects the main pulmonary artery and its left and right branches (Wang and Chesler 2011; Tan et al. 2014; Wang et al. 2013; Champion et al. 2009). To account for this behavior, the coefficient of vascular elasticity β described in (5)–(6) was scaled up as β = ηβbase where η ≥ 1 is a unitless coefficient representing the degree of stiffening with respect to baseline (normotensive) conditions attained at η = 1. We refer to this multiplicative factor as the normalized pulmonary vascular stiffness. Pulmonary vascular stiffening also affects the smaller arteries and arterioles of the pulmonary bed in both lungs (Wang and Chesler 2011; Wang et al. 2013; Tan et al. 2014; Acosta et al. 2015a; Dragu et al. 2015). In our model of the pulmonary bed, the effective vascular wall distensibility is given by Γart as described by (22). This coefficient is inversely proportional to vascular wall stiffness and was scaled down as Γart = Γart,base /η for η ≥ 1. Progression of pulmonary hypertension also affects the right ventricular diastolic stiffness which prevents this chamber from properly relaxing in diastole (Rain et al. 2013; Wang et al. 2013; Bogaard et al. 2009; Dragu et al. 2015; Moraes et al. 2000; Champion et al. 2009; Komajda and Lam 2014; Kuehne et al. 2004; Hill et al. 2014). Diastolic function of our heart model (43)–(45) is characterized by the minimum elastance E min . This parameter was scaled up E min = η1/2 E min,base to represent a gradual stiffening of the right ventricle. The scaling η1/2 was chosen because the referenced medical literature suggests that the remodeling process of the right ventricle is a delayed response to pulmonary arterial stiffening. We remark that the stiffening of these three components (pulmonary arteries, pulmonary arterioles, right ventricle) simultaneously is a simple but arbitrary model for this disease. Therefore, in Sect. 3, we explore the effect of stiffening each of these three components separately and also combined. 2.8 Cardiac compensation The effects of pulmonary vascular stiffening on the circulatory system include elevated pulmonary arterial pressure and high pulmonary vascular resistance, both of which can cause a reduction in cardiac output (Badesch et al. 2009; McGoon and Kane 2009; Bogaard et al. 2009; McLaughlin and McGoon 2006; Naeije and Huez 2007). However, to ensure that the metabolic demand of the organs is met, physiologic reflexes attempt to compensate for impaired cardiac performance by increasing heart rate (Provencher et al. 2006; Hildenbrand et al. 2012; Schuuring et al. 2013) and myocardial contractility (Katz and Lorell 2000; Bogaard et al. 2009; Rain et al. 2013). To account for these additional physiologic effects, two extreme mechanical scenarios are considered. The first scenario assumes no cardiac compensa-

Cardiovascular mechanics in the early stages of pulmonary hypertension: a computational study Table 1 Physiologic parameters to be included in the quantification of potential diagnostic indicators for early stage pulmonary hypertension Parameter

Abbrev.

Description and method of calculation

Right cardiac power

RCP

Mechanical power delivered by the right side of the heart given by the area enclosed in the PV loop multiplied by the heart rate

RV systolic pressure

RVPsys

Right ventricular pressure at systole

RV end-systolic volume

RVESV

Minimum volume attained by right ventricle during cardiac cycle

RV end-diastolic volume

RVEDV

Maximum volume attained by right ventricle during cardiac cycle

RV ejection fraction

RVEF

Stroke volume divided by the right ventricular end-diastolic volume

Mean pulmonary arterial pressure

MPAP

Mean blood pressure in pulmonary artery

PA pulse pressure

PPAP

Systolic pressure minus diastolic pressure in pulmonary artery

PA area distensibility

PAAD

Maximum change in cross-sectional area of the pulmonary artery divided by the minimum area during cardiac cycle

Pulmonary vascular resistance

PVR

Resistance to flow in pulmonary system given by the difference between mean pulmonary and mean left atrial pressures, divided by the cardiac output

Pulmonary vascular compliance

PVC

Compliance calculated as the stroke volume divided by the pulse pressure of the pulmonary artery

RC time

RC

The product of the pulmonary vascular resistance and the pulmonary vascular compliance

Pulmonary vascular impedance

PVI

Absolute value of the pulmonary vascular impedance (53) at frequency f = 70 bpm

tion and allows the cardiac output to decrease with increasing pulmonary hypertension. The second scenario assumes full cardiac compensation is present, thereby maintaining constant cardiac output with worsening pulmonary hypertension. This is achieved by increasing heart rate and cardiac contractility. In this second scenario, 50 percent of the required increase in cardiac output is achieved by increasing heart rate (HR), while the other 50 percent is attained by proportionally increasing maximum elastance (E max ) of all heart chambers. These two parameters of the heart model are described in precise terms in Sect. 2.5. The baseline value for the heart rate is 70 beats/min. The baseline values for the maximum elastance are displayed in Table 4, and the compensated values for maximum elastance and heart rate are displayed in Table 5. We acknowledge this model for cardiac compensation (half from increasing heart rate and half from increasing cardiac contractility) is a heuristic choice which we are required to make given the lack of quantitative clinical data regarding this compensatory process. However, as discussed in Sect. 3, this choice leads to qualitative agreement with clinical studies of pulmonary hypertension. 2.9 Relative change in diagnostic indicators Studies suggest several hemodynamic quantities relevant to the diagnosis of pulmonary hypertension (Dragu et al. 2015; Chemla et al. 2002; McGoon and Kane 2009; Sanz et al. 2009, 2012; Badesch et al. 2009). The hemodynamic variables that we explore are listed in Table 1 with their respective abbreviations and descriptions. The potential of these variables to be good diagnostic indicators is directly dependent

on how sensitive they are to the stiffening of the pulmonary vasculature. This sensitivity is quantified in terms of the relative change in a potential diagnostic indicator Y = Y (η) to the relative changes in stiffness parameter η,  Y η η dY (η) ≈ , (52) S(η) = Y (η) dη Y η which is the ratio of relative change in Y over the relative change in η. We concentrate on results corresponding to early stage pulmonary hypertension by evaluating (52) at η = 1.

3 Results and discussion Numerical simulations of the hemodynamic model described in Sect. 2 were carried out to measure the change in physiologic parameters as a function of increasing pulmonary hypertension, quantified by the normalized pulmonary stiffness η = 1, 2, 4, 8, 16. This range of values for the variable η is supported by studies, indicating that the stiffness of vascular tissues in the pulmonary arterial tree can increase more than one order of magnitude through the progression of pulmonary hypertension (Sanz et al. 2009; Wang and Chesler 2011; Tan et al. 2014; Wang et al. 2013; Acosta et al. 2015a; Dragu et al. 2015; Lankhaar et al. 2008). Section 3.1 considers the effect induced by the separate stiffening of three components of the pulmonary circulation. These three components are the major pulmonary arteries (1D-model), the pulmonary arterioles (0D-model), and the right ventricle (0D-model). The model for the stiffening of these three components is described in Sect. 2.7. Then, Sect. 3.2 illustrates

123

60

60

50

50

50

40

30

20

10

Pressure (mmHg)

60

Pressure (mmHg)

Pressure (mmHg)

S. Acosta et al.

40

30

20

10

0 30.2

30.3

30.4

30.5

30.6

30.7

30

20

10

0 30.1

40

0

30.8

30.1

30.2

30.3

30.4

30.5

30.6

30.7

30.8

30.1

30.2

30.3

30.4

30.5

Time (s)

Time (s)

Time (s)

(a)

(b)

(c)

30.6

30.7

30.8

600

600

500

500

500

400 300 200

Flow Rate (ml/s)

600

Flow Rate (ml/s)

Flow Rate (ml/s)

Fig. 5 Progression of the pressure waveform in the main pulmonary artery as three components of the pulmonary system stiffen. a Stiffening of the major pulmonary arteries. b Stiffening of the pulmonary arterioles. c Stiffening of the right ventricle

400 300 200

400 300 200

100

100

100

0

0

0

-100

30.1

30.2

30.3

30.4

30.5

30.6

30.7

30.8

-100

30.1

30.2

30.3

30.4

30.5

30.6

30.7

30.8

-100

30.1

30.2

30.3

30.4

30.5

Time (s)

Time (s)

Time (s)

(a)

(b)

(c)

30.6

30.7

30.8

Fig. 6 Progression of the flow waveform in the main pulmonary artery as three components of the pulmonary system stiffen. a Stiffening of the major pulmonary arteries. b Stiffening of the pulmonary arterioles. c Stiffening of the right ventricle

the effects obtain from stiffening all three components simultaneously. 3.1 Separate effect from the stiffening of three components of the pulmonary circulation In this subsection, we explore the effect of separately stiffening (a) the pulmonary arteries, (b) the pulmonary arterioles, and (c) the right ventricle on pulmonary artery waveforms and on ventricle pressure–volume loops. These cases correspond to plots (a), (b), and (c) in Figs. 5, 6, 7, and 8, respectively, and provide some insight into variations of these model outputs from which many of the diagnostic indicators are derived. Figure 5 displays pressure waveforms in the pulmonary artery. We remark that the systolic part of the waveform varies, while the diastolic part remains close to the baseline results (η = 1) in all three cases. The largest variation occurs in case (a): an increase in pulmonary artery stiffness leads to a large increase in the maximum pressure. Case (b) shows the same trend for an increase in pulmonary arteriole stiffness, but the change in maximum pressure is more mild. In

123

case (c), pressure decreases as stiffness of the right ventricle increases. This may be partially explained by the impact of ventricular stiffness on flow into the chamber during diastole. More precisely, normalized pulmonary stiffness η scales the minimum elastance of the ventricle; as η increases, the ventricle becomes less compliant, hindering its ability to stretch and fill with blood. The ventricle ejects less blood into the pulmonary artery, leading to a decrease in pressure. In Fig. 6, we investigate the effect of increasing stiffness on flow waveforms in the pulmonary artery. Maximum flow in each case decreases as η increases. In cases (a) and (b), this trend may be justified by an increase in afterload felt by the right ventricle, as a result of the increased stiffness in the pulmonary artery or arterioles. In case (c), the stiffer right ventricle ejects less blood, thereby decreasing flow. Pressure–volume loops for the right and left ventricles are given in Figs. 7 and 8, respectively. The trends in the right ventricle loops align well with our observations from the waveforms in the pulmonary artery. In particular, the results in case (a) of Fig. 7 show an increase in pressure as stiffness of the pulmonary artery increases. The loops in case (c) of Fig. 7 depict a substantial decrease in end-diastolic volume,

80

80

70

70

70

60

60

60

50 40 30

Pressure (mmHg)

80

Pressure (mmHg)

Pressure (mmHg)

Cardiovascular mechanics in the early stages of pulmonary hypertension: a computational study

50 40 30

50 40 30

20

20

20

10

10

10

0 40

60

80

100

120

140

0 40

160

60

Volume (ml)

80

100

120

140

0 40

160

60

Volume (ml)

(a)

80

100

120

140

160

Volume (ml)

(b)

(c)

160

160

140

140

140

120

120

120

100 80 60

Pressure (mmHg)

160

Pressure (mmHg)

Pressure (mmHg)

Fig. 7 Progression of the right ventricular PV loop as three components of the pulmonary system stiffen. a Stiffening of the major pulmonary arteries. b Stiffening of the pulmonary arterioles. c Stiffening of the right ventricle

100 80 60

100 80 60

40

40

40

20

20

20

0 40

60

80

100

Volume (ml)

120

140

160

0 40

60

80

100

Volume (ml)

(a)

(b)

120

140

160

0 40

60

80

100

120

140

160

Volume (ml)

(c)

Fig. 8 Progression of the left ventricular PV loop as three components of the pulmonary system stiffen. a Stiffening of the major pulmonary arteries. b Stiffening of the pulmonary arterioles. c Stiffening of the right ventricle

as η increases. This result is consistent with the notion that right ventricular stiffness impedes diastolic filling. The left ventricle loops in case (c) of Fig. 8 show a decrease in pressure and end-diastolic volume for increased stiffness of the right ventricle, but are otherwise unaffected in cases (a) and (b). This is partially explained by examining the right ventricle loops in cases (a) and (b) of Fig. 7. In these cases, the end-diastolic volume of the right ventricle does not vary too much, preserving the stroke volume. In contrast, the stroke volume of both ventricles does decrease in case (c) of Figs. 7 and 8. In summary, stiffening of the pulmonary arteries or the pulmonary arterioles, cases (a) and (b), respectively, leads to an increase in the effective resistance felt by the right ventricle. These changes cause a pressure increase in the pulmonary circulation, but do not greatly impact the pressure–volume dynamics of the left ventricle. An increase in the stiffness of the right ventricle, case (c), impedes diastolic filling, implying a decrease in stroke volume for both the left and right ventricles. Stiffening of the pulmonary arterioles has the smallest affect on the waveforms in the pulmonary artery and the ventricle pressure–volume loops.

3.2 Combined effect from the stiffening of three components of the pulmonary circulation Figure 9 illustrates the evolution of the right ventricular PV loops as a function of increasing pulmonary arterial stiffness, with and without cardiac compensation. In both cases, an increase in the ventricular systolic pressure and a decrease in stroke volume were observed, as seen clinically (Kuehne et al. 2004). Figures 10 and 11 display the heart rate, (normalized) cardiac contractility, and cardiac power needed to maintain the cardiac output constant. These figures illustrate how much work must be performed by the heart in order to maintain the baseline cardiac output. We note that including cardiac compensation in the simulation is important given that higher cardiac workload may eventually cause failure of the right heart (Naeije and Huez 2007; Bogaard et al. 2009; Dragu et al. 2015; Rain et al. 2013; McGoon and Kane 2009; VonkNoordegraaf et al. 2013; Wang et al. 2013; Matthews and McLaughlin 2008; Naeije and Manes 2014; Aronson et al. 2013; Voelkel et al. 2006; Komajda and Lam 2014; Kuehne et al. 2004). However, since the model was not designed

123

80

80

70

70

60

60

Pressure (mmHg)

Pressure (mmHg)

S. Acosta et al.

50 40 30

50 40 30

20

20

10

10

0 40

60

80

100

120

140

0 40

160

60

Volume (ml)

80

100

120

140

160

Volume (ml)

Fig. 9 Changes in the right ventricular PV loops as the pulmonary arterial stiffening increases proportionally to the factor η = 1, 2, 4, 8, 16. Left panel: with cardiac compensation (higher heart rate and higher contractility as described in Sect. 2.8). Right panel: without cardiac compensation 2.5

110

With Cardiac Compensation Without Cardiac Compensation

100 90 80 70 60 10 0

Normalized Contractility

Heart Rate (beats/min)

120

With Cardiac Compensation Without Cardiac Compensation

2

1.5

1

0.5 10 0

10 1

Normalized Pulmonary Vascular Stiffness

10 1

Normalized Pulmonary Vascular Stiffness

Fig. 10 Changes in heart rate and cardiac contractility (normalized by its value at η = 1; see Table 5) as functions of pulmonary arterial stiffening

With Cardiac Compensation Without Cardiac Compensation

5.5

5

4.5

0.9

R. Heart Power (W)

Cardiac Output (L/min)

6 With Cardiac Compensation Without Cardiac Compensation

0.8 0.7 0.6 0.5 0.4 0.3

4 10 0

10

1

Normalized Pulmonary Vascular Stiffness

10 0

10 1

Normalized Pulmonary Vascular Stiffness

Fig. 11 Changes in cardiac output (CO) and right heart cardiac power (RCP) as functions of pulmonary arterial stiffening

to simulate the development of cardiac failure, our results are limited to relatively early stages of pulmonary hypertension. Simulation results are qualitatively consistent with empirical comparisons between control subjects and patients with pulmonary hypertension. In Kuehne et al. (2004), the cardiac output was found to decrease to (75.9 ± 8.3)% and the myocardial contractility increased to (184.0 ± 16.2)% of the respective control values. The increase in cardiac contractil-

123

ity indicates that the patient’s body is trying to compensate for the low level of perfusion. However, the compensation does not fully preserve the cardiac output at baseline value. The strength of the compensation is patient dependent and difficult to estimate a priori. Therefore, we choose to simulate extreme scenarios, namely full and absent cardiac compensation, as described in Sect. 2.8, to illustrate bounds for what a patient’s likely response would be. We note that both scenarios render similar behavior for early stages of

Cardiovascular mechanics in the early stages of pulmonary hypertension: a computational study 70 With Cardiac Compensation Without Cardiac Compensation

80 70 60 50

RV end-sys vol (mL)

RV Sys Pressure (mmHg)

90

With Cardiac Compensation Without Cardiac Compensation

65

60

55

40 30 10 0

50 10 0

10 1

Normalized Pulmonary Vascular Stiffness

10 1

Normalized Pulmonary Vascular Stiffness

Fig. 12 Changes in right ventricular systolic pressure (RVPsys ) and right ventricular end-systolic volume (RVESV) as functions of pulmonary arterial stiffening

RV end-dia vol (mL)

With Cardiac Compensation Without Cardiac Compensation

130

120

RV ejection frac (%)

65 140

With Cardiac Compensation Without Cardiac Compensation

60

55

50

110 10 0

10 1

Normalized Pulmonary Vascular Stiffness

45 10 0

10 1

Normalized Pulmonary Vascular Stiffness

Fig. 13 Changes in right ventricular end-diastolic volume (RVEDV) and right ventricular ejection fraction (RVEF) as functions of pulmonary arterial stiffening

pulmonary hypertension characterized by 1 ≤ η  4. Therefore, the magnitude of cardiac compensation does not affect significantly our conclusions concerning the early stages of pulmonary hypertension. Figures 12 and 13 illustrate right ventricular pressure, volume, and ejection fraction. The left plot in Fig. 12 shows how the right ventricular systolic pressure is higher when cardiac compensation takes place. The trend in the right plot of Fig. 12 may be described in the following way. With or without cardiac compensation, the initial decrease in pulmonary vascular compliance prevents the right ventricle from fully ejecting its content, leading to an increasing RVESV. However, at some point the cardiac compensation exerts significantly higher muscular strength during systole (also illustrated in the right panels of Figs. 10 and 11). At that point, the contractions are strong enough to reduce the chamber’s end-systolic volume when compared to dynamics without cardiac compensation. The two scenarios seem to diverge at η ≈ 4. In Fig. 13, the RVEDV and the ejection fraction decrease as η increases, since the minimum elastance of the ventricle increases. More precisely, the right ventricle becomes less compliant, degrading its ability to expand and fill up with blood during diastole. At relatively mild over-

load, the right ventricular wall remodels toward stiffer and thicker composition. However, chronic and excessive pressure overload eventually triggers irreversible dilation, loss of contractility, and overt heart failure (Bogaard et al. 2009; Voelkel et al. 2006; Vonk-Noordegraaf et al. 2013; Wang et al. 2013). We wish to note again that in these simulations we do not account for heart failure manifested as irreversible dilation of the ventricle at severe stages of pulmonary hypertension. Therefore, this dilatory phenomenon is not reflected in Figs. 12 and 13. Due to the nonlinear nature of the model (24)–(26), the stiffening of the pulmonary arterial network affects indirectly other parameters, including pulmonary vascular resistance, pulmonary vascular compliance, the pulmonary RC time, and pulmonary vascular impedance. Figure 14 illustrates the progression of pulmonary vascular resistance (PVR) and pulmonary vascular compliance (PVC). These two properties exhibit nearly inverse behavior as the vascular stiffening increases. The product of PVR and PVC, known as the RC time, is displayed in Fig. 15 along with the absolute value of the pulmonary vascular impedance (PVI). We note that although PVR and PVC progress in the opposite direction, their product does not stay identically constant. In fact, at the

123

S. Acosta et al. 5 With Cardiac Compensation Without Cardiac Compensation

4

3

2

PVC (mL / mmHg)

PVR (mmHg min / L)

5

With Cardiac Compensation Without Cardiac Compensation

4

3

2

1 1 10 0

10 0

10 1

10 1

Normalized Pulmonary Vascular Stiffness

Normalized Pulmonary Vascular Stiffness

Fig. 14 Changes in pulmonary vascular resistance (PVR) and pulmonary vascular compliance (PVC) as functions of pulmonary arterial stiffening 3 With Cardiac Compensation Without Cardiac Compensation

0.5 0.4 0.3 0.2

PVI (mmHg min / L)

Pulmonary RC time (s)

0.6

With Cardiac Compensation Without Cardiac Compensation

2.5 2 1.5 1 0.5

0.1 10 0

0 10 0

10 1

Normalized Pulmonary Vascular Stiffness

10 1

Normalized Pulmonary Vascular Stiffness

Fig. 15 Changes in pulmonary RC time (RC) and the absolute value of the pulmonary vascular impedance (|Z |) at the cardiac frequency as functions of pulmonary arterial stiffening 80 With Cardiac Compensation Without Cardiac Compensation

25

20

15 10 0

10 1

Normalized Pulmonary Vascular Stiffness

Pulse PA Pressure (mmHg)

Mean PA Pressure (mmHg)

30

70

With Cardiac Compensation Without Cardiac Compensation

60 50 40 30 20 10 10 0

10 1

Normalized Pulmonary Vascular Stiffness

Fig. 16 Changes in mean pulmonary arterial pressure (MPAP) and pulse pulmonary arterial pressure (PPAP)

early stages (η ≈ 1), the RC time experiences a significant variation. A similar finding has been reported in clinical studies of thromboembolic pulmonary hypertension (MacKenzie Ross et al. 2013). However, toward the advanced stages of pulmonary arterial stiffening (η  10), the RC time becomes asymptotically constant. This result agrees with other empirical and theoretical studies (Lankhaar et al. 2008, 2006; Reuben 1971; Saouti et al. 2010a, b; Acosta et al. 2015a). The definition of PVI (given by (53)) and the trend in Fig. 15

123

indicate that the increase in PVR and decrease in RC time compound to a significant increase in PVI. This makes PVI very sensitive to changes in pulmonary vascular stiffness. Finally, we see in Fig. 16 that the mean pulmonary arterial pressure and the pulse pressure in the pulmonary artery increase along with the stiffening. However, the variation is much more pronounced in the pulse pressure, which is directly related to the distensibility of this artery as displayed in Fig. 17.

Cardiovascular mechanics in the early stages of pulmonary hypertension: a computational study

PA Area Distensibility

0.5 With Cardiac Compensation Without Cardiac Compensation

0.4 0.3 0.2 0.1 0 10 0

10 1

Normalized Pulmonary Vascular Stiffness

Fig. 17 Changes in pulmonary arterial area distensibility (PAAD) as functions of pulmonary arterial stiffening

3.3 Discussion on the analysis of diagnostic indicators Figure 18 displays the sensitivity S defined by (52) for various hemodynamic variables described in Table 1. Results indicate that at early stages of hypertension, the relative change in diagnostic indicators is almost independent of cardiac compensation. The parameters RVEDV, RVEF, RVESV, MPAP, PVR, RCP and RVPsys were found to have low sensitivity to the early stage of pulmonary arterial stiffening, making them poor candidates for early recognition of pulmonary hypertension. The analysis suggests that PAAD, RC, PVC, PPAP, and PVI may be more useful for early detection. These results agree with some clinical studies. For instance, when compared with conventional markers (such

as the 6-min walk test, cardiac output, right atrial pressure, RVEF, MPAP, and PVR), pulmonary vascular compliance (also known as capacitance) has been found to be a stronger predictor of mortality in patients with idiopathic pulmonary arterial hypertension (Mahapatra et al. 2006). Similar results for the compliance can also be found in the work of Lankhaar et al.; the authors used a 0D computational model of the pulmonary circulation, with parameters estimated from clinical data (Lankhaar et al. 2006). We remark that other variables not considered in this work, derived from wave intensity analysis, also serve as good indicators for pulmonary hypertension (Quail et al. 2015). Among the variables most sensitive to early pulmonary hypertension, PPAP, PVC, and PVI require invasive measurements of blood pressure to obtain reliable estimations of their values. However, PAAD and RC may be obtained noninvasively, without requiring absolute measurements of blood pressure. The pulmonary artery’s area distensibility (PAAD) can be directly estimated using MRI technology since it is a purely geometrical variable (Sanz et al. 2009). The RC time constant can be estimated using measurements of flow and luminal area. At a given frequency f of pressure oscillations, the relationship between the oscillatory flow Q( f ) and the oscillatory pressure P( f ) is P( f ) = Z ( f )Q( f ). The term Z ( f ), the impedance of the pulmonary vasculature, is given by

Z( f ) =

R . 1 + 2πi f RC

(53)

PVI PPAP PVC RC PAAD RVPsys RCP PVR MPAP RVEDV RVEF With Cardiac Compensation Without Cardiac Compensation

RVESV -0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

Sensitivity

Fig. 18 Relative change in diagnostic indicators to changes in pulmonary arterial stiffness as defined by (52). RVEDV right ventricular end-diastolic volume, RVEF right ventricular ejection fraction, RVESV right ventricular end-systolic volume, MPAP mean pulmonary arterial pressure, PVR pulmonary vascular resistance, RCP right cardiac

power, RVPsys right ventricular pressure at systole, PAAD pulmonary arterial area distensibility, RC: product of PVR and PVC; PVC pulmonary vascular compliance, PPAP pulse pulmonary arterial pressure. PVI Absolute value of the pulmonary vascular impedance at 70 bpm. These hemodynamic variables are described in Table 1

123

S. Acosta et al.

According to the state Eq. (5), the pressure P oscillates in phase with the cross-sectional area A. Therefore, the RC time and the phase shift ϕ( f ) between cross-sectional area A and flow Q satisfies the following relationship RC =

tan ϕ( f ) . 2π f

(54)

We note that the state Eq. (5) does not account for viscoelasticity. The viscoelastic behavior of the vascular tissue induces a phase shift between the pressure P and the cross-sectional area A. This phase shift can be estimated and incorporated into (54), provided that the viscoelastic properties of these tissues are known. As a consequence, the RC time can be estimated using simultaneous waveform measurements of flow Q and cross-sectional area A to retrieve their phase difference. Both of these may be acquired using time-resolved magnetic resonance imaging (Sanz et al. 2009; Krishnamurthy et al. 2014; Dyverfeldt et al. 2015; Stankovic et al. 2014; Markl et al. 2014; Lungu et al. 2014). The heart rate is the natural choice for the frequency f since it is the most dominant frequency of the signals. However, it is conceivable that other frequencies could be employed to overdetermine (54) to render more robust estimations of the RC time.

4 Limitations The simulations presented in this work have two main limitations. First, our results are limited to early stages of pulmonary hypertension where the heart has not yet experienced irreversible remodeling and dilation known to be associated with heart failure (Bogaard et al. 2009; Voelkel et al. 2006; Vonk-Noordegraaf et al. 2013). This transition is still an area of active research (Wang et al. 2013). A failing heart dramatically affects the dynamics of the cardiovascular system, and this phenomenon must be taken into account to realistically simulate the advanced stages of pulmonary hypertension. The second limitation is that this analysis does not address situations where patient-to-patient variability may be significant. The magnitude of such variability may be much larger than the response estimated by our deterministic analysis. The methods associated with stochastic uncertainty quantification are beyond the scope of this paper, but will be addressed in future work. We also acknowledge the need for validation of the computational simulations and the results obtained in this paper. Although the results are qualitatively consistent with other computational and clinical studies, there is a need to quantitatively compare the simulated results with clinical data. There are other areas of improvement such as the incorporation of more realistic models of the heart (Rausch et al. 2011) and more detailed stress–strain relations for

123

the 1D and 0D models of vasculature. Following the work of Pursell et al. (2016), Steelet et al. (2011), Valdez-Jasso et al. (2011), and Valdez-Jasso et al. (2009), the resolution of both circumferential and axial strains and/or the dynamics of viscoelastic vessel walls could be included in the model.

5 Conclusions We have constructed a computational model to analyze the hemodynamics of pulmonary hypertension in a closed-loop network capable of representing the major vessels and organs of the human body. The mechanical impact of progressive stiffening of the pulmonary tree (arteries and arterioles) on the rest of the circulatory system was quantified. Since the response to chronic pulmonary hypertension is patient dependent, two extreme scenarios were simulated: a state with constant cardiac output requiring cardiac compensation and variable cardiac output without any cardiac compensation. The cardiac demand of a patient will likely lie between these extreme bounds. For sufficiently early stages of pulmonary hypertension, these two bounds are similar enough for us to draw meaningful conclusions regardless of the patient’s actual physiologic response. The simulated results qualitatively agree with clinical observations which reveal the potential application of this model to assist in the design of diagnostic tools. The model was interrogated to determine hemodynamic variables most highly impacted by the early stages of pulmonary vascular stiffening. The purpose of this investigation was to provide quantitative estimates of the sensitivity of physiologic features to the early stages of pulmonary hypertension. Early detection of pulmonary hypertension is an elusive goal, and better methods for noninvasive tests are still needed (Humbert et al. 2010a, b, 2012; Thenappan et al. 2010; Gibbs 2007; Rich and Rich 2014; Hunter et al. 2008, 2011; Wang and Chesler 2013; Mahapatra et al. 2006). Conventionally, the diagnosis of pulmonary hypertension has been based on detecting elevated mean pulmonary arterial pressure and high pulmonary vascular resistance (Badesch et al. 2009). As shown in Fig. 18, these two variables exhibit low sensitivity to small changes in pulmonary vascular stiffness. Therefore, they are not good candidates to detect the early development of this disease. Our results suggest that the cardiovascular variables most affected by early stages of pulmonary stiffening are: pulse pressure in the pulmonary artery, pulmonary vascular compliance, the RC time of the pulmonary system, and the area distensibility of the pulmonary artery. Of these variables, it is believed that the RC time and area distensibility can be estimated noninvasively using MRI imaging.

Cardiovascular mechanics in the early stages of pulmonary hypertension: a computational study Acknowledgements C. Puelz was supported by a fellowship from the Keck Center of the Gulf Coast Consortia, on the Training Program in Biomedical Informatics, US National Library of Medicine T15LM007093. Compliance with ethical standards Conflict of interest The authors declare that they have no conflict of interest.

6 Appendix: Parameters for baseline model In the tables below, we display the parameters used in our models. The geometrical parameters for the vessels are from Mynard (2011). Initial guesses for the valve parameters were taken from Blanco and Feijóo (2013). All other parameters were manually tuned, with initial guesses from Mynard (2011).

Table 3 continued Rcap,ref Cart,ref Cven,ref Vref (104 dyne s/cm5 ) (10−4 cm5 /dyne) (10−4 cm5 /dyne) (cm3 ) R. arm

1.50

0.20

2.00

120

L. arm

1.50

0.20

2.00

120 100

Thorax

0.40

1.00

10.00

Spleen

0.50

0.40

4.00

70

Stomach

2.00

0.20

2.00

70

Intestines

0.30

1.00

10.00

250

R. kidney

0.25

1.00

10.00

130 130

L. kidney

0.25

1.00

10.00

Liver

0.55

0.50

5.00

425

Pelvic region

1.50

0.30

3.00

160

R. thigh

1.50

0.30

3.00

170

L. thigh

1.50

0.30

3.00

170

R. lower leg

1.00

0.50

5.00

85

6.1 Vessel network parameters

L. lower leg

1.00

0.50

5.00

85

See Table 2.

R. lung

0.40 × 10−2

10.00

20.00

700

L. lung

0.45 × 10−2

10.00

20.00

700

Table 2 Parameters for the different one-dimensional vessel networks p0

(dyne/cm2 )

k1

(dyne/cm2 )

k2

(cm−1 )

k3

(dyne/cm2 )

Systemic arteries

1.20 × 105

1.41 × 106

−0.9

0.0

Systemic veins

6.67 × 103

3.00 × 105

−4.0

2.0×104

Pulmonary 1.47 × 104 arteries

5.00 × 105

−5.0

5.0×105

Pulmonary 1.33 × 104 veins

2.00 × 105

−5.0

2.0×105

1.33 × 104

3.00 × 105

−2.5

1.5×104

Portal veins

See details in Sect. 2.4 The total volume of blood contained in the simulated cardiovascular system (all vessel networks, organ beds and four cardiac chambers) adds to 5186.99 ml

6.3 Heart parameters See Tables 4, 5 and 6.

Table 4 Baseline parameters for the heart chambers L. ventricle L. atrium R. ventricle R. atrium

6.2 Organ bed parameters See Table 3.

E min (103 dyne/cm2 ) 1.2

1.2

1.2

1.2

E max (103 dyne/cm2 ) 120

12

120

12 3

m1

3

3

3

m2

27

27

27

27

τ1 (T)

0.27

0.11

0.27

0.11

0.45

0.18

0.45

0.18

τ2 (T)

Table 3 Organ bed parameters Rcap,ref Cart,ref Cven,ref Vref (104 dyne s/cm5 ) (10−4 cm5 /dyne) (10−4 cm5 /dyne) (cm3 ) Head

0.50

0.30

3.00

120

Neck

1.00

0.40

4.00

80

R. shoulder

2.00

0.15

1.50

120

L. shoulder

2.00

0.15

1.50

120

(10−4

s/cm3 )

5

5

5

5

tonset (T)

0

0.85

0

0.85

V 0 (cm3 )

24

12

46

23

K

See details in Sect. 2.5 The parameters τ1 , τ2 , and tonset are given in terms of the heart period T , and the baseline value for the heart rate (HR) is 70 beats/min

123

S. Acosta et al. Table 5 Parameters for the heart chambers affected by the cardiac compensation described in Sect. 2.8 η=1

η=2

η=4

η=8

η = 16

LA E max (103 dyne/cm2 ) 12.00

14.64

17.16

20.88

26.16

LV E max

(103

dyne/cm2 )

120.00 146.40 171.60 208.80 261.60

RA E max (103 dyne/cm2 ) 12.00

14.64

17.16

20.88

26.16

RV E max (103 dyne/cm2 ) 120.00 146.40 171.60 208.80 261.60 HR (beat/s)

70

73

80

90

105

All the other parameters listed in Table 4 remain the same Table 6 Parameters for the heart valves L. interior L. exterior R. interior R. exterior γ (10−3 cm2 /dyne)

5

3

5

3

B (10−2 dyne s2 /cm8 ) 3

3

3

3

5

5

5

5

1

1

1

1

L

(10−3

dyne

s2 /cm5 )

R (10−3 dyne s/cm5 ) See details in Sect. 2.5

References Acosta S, Penny D, Rusin C (2015) An effective model of blood flow in capillary beds. Microvasc Res 100:40–47. doi:10.1016/j.mvr. 2015.04.009 Acosta S, Puelz C, Rivière B, Penny D, Rusin C (2015) Numerical method of characteristics for one-dimensional blood flow. J Comput Phys 294:96–109. doi:10.1016/j.jcp.2015.03.045 Alastruey J, Moore SM, Parker KH, David T, Peiró J, Sherwin SJ (2008) Reduced modelling of blood flow in the cerebral circulation: coupling 1-D, 0-D and cerebral auto-regulation models. Int J Numer Methods Fluids 56(8):1061–1067. doi:10.1002/fld Aronson D, Darawsha W, Atamna A, Kaplan M, Makhoul B, Mutlak D, Lessick J, Carasso S, Reisner S, Agmon Y, Dragu R, Azzam Z (2013) Pulmonary hypertension, right ventricular function, and clinical outcome in acute decompensated heart failure. J Card Fail 19(10):665–671. doi:10.1016/j.cardfail.2013.08.007 Audebert C, Bucur P, Bekheit M, Vibert E, Vignon-Clementel I, Gerbeau JF (2017) Kinetic scheme for arterial and venous blood flow, and application to partial hepatectomy modeling. Comput Methods Appl Mech Eng 314:102–125. doi:10.1016/j.cma.2016.07.009 Badesch D, Champion H, Gomez Sanchez M, Hoeper M, Loyd J, Manes A, McGoon M, Naeije R, Olschewski H, Oudiz R, Torbicki A (2009) Diagnosis and assessment of pulmonary arterial hypertension. J Am Coll Cardiol 54(1):S55–S66. doi:10.1016/j.jacc.2009. 04.011 Blanco P, Feijóo R (2013) A dimensionally-heterogeneous closed-loop model for the cardiovascular system and its applications. Med Eng Phys 35(5):652–67. doi:10.1016/j.medengphy.2012.07.011 Bogaard H, Abe K, Noordegmaf A, Voelkel N (2009) The right ventricle under pressure. Cellular and molecular mechanisms of right-heart failure in pulmonary hypertension. Chest 135:794–804. doi:10. 1378/chest.08-0492 ˇ c S, Kim E (2003) Mathematical analysis of the quasilinear effects Cani´ in a hyperbolic model blood flow through compliant axi-symmetric vessels. Math Methods Appl Sci 26(14):1161–1186. doi:10.1002/ mma.407 Champion H, Michelakis E, Hassoun P (2009) Comprehensive invasive and noninvasive approach to the right ventricle-pulmonary

123

circulation unit. Circulation 120(11):992–1007. doi:10.1161/ CIRCULATIONAHA.106.674028 Chemla D, Castelain V, Hervé P, Lecarpentier Y, Brimioulle S (2002) Haemodynamic evaluation of pulmonary hypertension. Eur Respir J 20:1314–1331. doi:10.1183/09031936.02.00068002 Dragu R, Rispler S, Habib M, Sholy H, Hammerman H, Galie N, Aronson D (2015) Pulmonary arterial capacitance in patients with heart failure and reactive pulmonary hypertension. Eur J Heart Fail 17(1):74–80. doi:10.1002/ejhf.192 Dyverfeldt P, Bissell M, Barker A, Bolger A, Carlhäll CJ, Ebbers T, Francios C, Frydrychowicz A, Geiger J, Giese D, Hope M, Kilner P, Kozerke S, Myerson S, Neubauer S, Wieben O, Markl M (2015) 4D flow cardiovascular magnetic resonance consensus statement. J Cardiovasc Magn Reson 17:72. doi:10.1186/s12968-015-0174-5 Farber H, Loscalzo J (2004) Pulmonary arterial hypertension. N Engl J Med 351:1655–1665. doi:10.1056/NEJMra035488 Formaggia L, Nobile F, Quarteroni A (2002) A one dimensional model for blood flow: application to vascular prosthesis. In: Math. Model. Numer. Simul. Contin. Mech. Springer, Berlin, pp 137– 153. doi:10.1007/978-3-642-56288-4_10 Formaggia L, Lamponi D, Tuveri M, Veneziani A (2006) Numerical modeling of 1D arterial networks coupled with a lumped parameters description of the heart. Comput Methods Biomech Biomed Eng 9(5):273–288. doi:10.1080/10255840600857767 Gaine S, Rubin L (1998) Primary pulmonary hypertension. Lancet 9129:719–725. doi:10.1016/S0140-6736(98)02111-4 Gibbs JSR (2007) Making a diagnosis in PAH. Eur Respir Rev 16(102):8–12. doi:10.1183/09059180.00010203 Hildenbrand FF, Fauchère I, Huber LC, Keusch S, Speich R, Ulrich S (2012) A low resting heart rate at diagnosis predicts favourable long-term outcome in pulmonary arterial and chronic thromboembolic pulmonary hypertension. A prospective observational study. Respir Res 13(1):76. doi:10.1186/1465-9921-13-76 Hill M, Simon M, Valdez-Jasso D, Zhang W, Champion H, Sacks M (2014) Structural and mechanical adaptations of right ventricle free wall myocardium to pressure overload. Ann Biomed Eng 42(12):2451–2465. doi:10.1007/s10439-014-1096-3 Hughes T, Lubliner J (1973) On the one-dimensional theory of blood flow in the larger vessels. Math Biosci 18:161–170. doi:10.1016/ 0025-5564(73)90027-8 Humbert M, Sitbon O, Chaouat A, Bertocchi M, Habib G, Gressin V, Yaïci A, Weitzenblum E, Cordier JF, Chabot F, Dromer C, Pison C, Reynaud-Gaubert M, Haloun A, Laurent M, Hachulla E, Cottin V, Degano B, Jaïs X, Montani D, Souza R, Simonneau G (2010a) Survival in patients with idiopathic, familial, and anorexigen-associated pulmonary arterial hypertension in the modern management era. Circulation 122(2):156–163. doi:10. 1161/CIRCULATIONAHA.109.911818 Humbert M, Sitbon O, Yaïci A, Montani D, O’Callaghan D, Jaïs X, Parent F, Savale L, Natali D, Günther S, Chaouat A, Chabot F, Cordier JF, Habib G, Gressin V, Jing ZC, Souza R, Simonneau G (2010b) Survival in incident and prevalent cohorts of patients with pulmonary arterial hypertension. Eur Respir J 36(3):549–555. doi:10.1183/09031936.00057010 Humbert M, Gerry Coghlan J, Khanna D (2012) Early detection and management of pulmonary arterial hypertension. Eur Respir Rev 21(126):306–312. doi:10.1183/09059180.00005112 Hunter K, Gross J, Lanning C, Kirby K, Dyer K, Ivy D, Shandas R (2008) Noninvasive methods for determining pulmonary vascular function in children with pulmonary arterial hypertension: application of a mechanical oscillator model. Congenit Heart Dis 3:106–116 Hunter K, Lammers S, Shandas R (2011) Pulmonary vascular stiffness: measurement, modeling, and implications in normal and hypertensive pulmonary circulations. Compr Physiol 1(3):1413–1435. doi:10.1002/cphy.c100005

Cardiovascular mechanics in the early stages of pulmonary hypertension: a computational study Katz A, Lorell B (2000) Regulation of cardiac contraction and relaxation. Circulation 102:69–74. doi:10.1161/01.CIR.102.suppl_4. IV-69 Kheyfets V, O’Dell W, Smith T, Reilly J, Finol E (2013) Considerations for numerical modeling of the pulmonary circulation—a review with a focus on pulmonary hypertension. J Biomech Eng 135(6):061011–061015. doi:10.1115/1.4024141 Komajda M, Lam C (2014) Heart failure with preserved ejection fraction: a clinical dilemma. Eur Heart J 35(16):1022–1032. doi:10. 1093/eurheartj/ehu067 Krishnamurthy R, Cheong B, Muthupillai R (2014) Tools for cardiovascular magnetic resonance imaging. Cardiovasc Diagn Ther 4(2):104–125. doi:10.3978/j.issn.2223-3652.2014.03.06 Kuehne T, Yilmaz S, Steendijk P, Moore P, Groenink M, Saaed M, Weber O, Higgins C, Ewert P, Fleck E, Nagel E, Schulze-Neick I, Lange P (2004) Magnetic resonance imaging analysis of right ventricular pressure-volume loops. Circulation 110(14):2010–2016. doi:10.1161/01.CIR.0000143138.02493.DD Lankhaar JW, Westerhof N, Faes T, Marques K, Marcus J, Postmus P, Vonk-Noordegraaf A (2006) Quantification of right ventricular afterload in patients with and without pulmonary hypertension. Am J Physiol Heart Circ Physiol 291(4):H1731–7. doi:10.1152/ ajpheart.00336.2006 Lankhaar JW, Westerhof N, Faes T, Gan C, Marques K, Boonstra A, van den Berg F, Postmus P, Vonk-Noordegraaf A (2008) Pulmonary vascular resistance and compliance stay inversely related during treatment of pulmonary hypertension. Eur Heart J 29(13):1688– 1695. doi:10.1093/eurheartj/ehn103 Lungu A, Wild J, Capener D, Kiely D, Swift A, Hose D (2014) MRI model-based non-invasive differential diagnosis in pulmonary hypertension. J Biomech 47(12):2941–2947. doi:10.1016/ j.jbiomech.2014.07.024 Lungu A, Swift A, Capener D, Kiely D, Hose R, Wild J (2016) Diagnosis of pulmonary hypertension from magnetic resonance imagingbased computational models and decision tree analysis. Pulm Circ 6(2):181–190. doi:10.1086/686020 MacKenzie Ross R, Toshner M, Soon E, Naeije R, Pepke-Zaba J (2013) Decreased time constant of the pulmonary circulation in chronic thromboembolic pulmonary hypertension. Am J Physiol Heart Circ Physiol 305:H259–H264. doi:10.1152/ajpheart.00128.2013 Mahapatra S, Nishimura R, Sorajja P, Cha S, McGoon M (2006) Relationship of pulmonary arterial capacitance and mortality in idiopathic pulmonary arterial hypertension. J Am Coll Cardiol 47(4):799–803. doi:10.1016/j.jacc.2005.09.054 Markl M, Schnell S, Barker A (2014) 4D flow imaging: current status to future clinical applications. Curr Cardiol Rep 16:481. doi:10. 1007/s11886-014-0481-8 Matthews JC, McLaughlin V (2008) Acute right ventricular failure in the setting of acute pulmonary embolism or chronic pulmonary hypertension: a detailed review of the pathophysiology, diagnosis, and management. Curr Cardiol Rev 4(1):49–59. doi:10.2174/ 157340308783565384 McGoon M, Kane G (2009) Pulmonary hypertension: diagnosis and management. Mayo Clin Proc 84(2):191–207. doi:10.1016/ S0025-6196(11)60828-8 McLaughlin V, McGoon M (2006) Pulmonary arterial hypertension. Circulation 114(13):1417–1431. doi:10.1161/ CIRCULATIONAHA.104.503540 Melicher V, Gajdošík V (2008) A numerical solution of a onedimensional blood flow model-moving grid approach. J Comput Appl Math 215(2):512–520 Mikelic A, Guidoboni G, Canic S (2007) Fluid-structure interaction in a pre-stressed tube with thick elastic walls I: the stationary Stokes problem. Netw Heterog Media 2(3):397–423. doi:10.3934/nhm. 2007.2.397

Moraes D, Colucci W, Givertz M (2000) Secondary pulmonary hypertension in chronic heart failure: the role of the endothelium in pathophysiology and management. Circulation 102:1718–1723. doi:10.1161/01.CIR.102.14.1718 Mynard J (2011) PhD Thesis: Computer modeling and wave intensity analysis of perinatal cardiovascular function and dysfunction. Ph.D. thesis, University of Melbourne Mynard J, Smolich J (2015) One-dimensional haemodynamic modeling and wave dynamics in the entire adult circulation. Ann Biomed Eng 43(6):1443–1460. doi:10.1007/s10439-015-1313-8 Naeije R, Huez S (2007) Right ventricular function in pulmonary hypertension: physiological concepts. Eur Heart J Suppl 9:H5– H9. doi:10.1093/eurheartj/sum023 Naeije R, Manes A (2014) The right ventricle in pulmonary arterial hypertension. Eur Respir Rev 23:476–487. doi:10.1183/ 09059180.00007414 Olufsen M (1999) Structured tree outflow condition for blood flow in larger systemic arteries. Am J Physiol Heart Circ Physiol 276:H257–H268. doi:10.1017/CBO9781107415324.004 Olufsen MS, Peskin CS, Kim WY, Pedersen EM, Nadim A, Larsen J (2000) Numerical simulation and experimental validation of blood flow in arteries with structured-tree outflow conditions. Ann Biomed Eng 28(11):1281–1299 Provencher S, Chemla D, Hervé P, Sitbon O, Humbert M, Simonneau G (2006) Heart rate responses during the 6-minute walk test in pulmonary arterial hypertension. Eur Respir J 27:114–120. doi:10. 1183/09031936.06.00042705 Punnoose L, Burkhoff D, Rich S, Horn E (2012) Right ventricular assist device in end-stage pulmonary arterial hypertension: insights from a computational model of the cardiovascular system. Prog Cardiovasc Dis 55(2):234–243. doi:10.1016/j.pcad.2012.07.008 Pursell E, Vélez-Rendón D, Valdez-Jasso D (2016) Biaxial properties of the left and right pulmonary arteries in a monocrotaline rat animal model of pulmonary arterial hypertension. J Biomech Eng 138(11):111004. doi:10.1115/1.4034826 Quail MA, Knight DS, Steeden JA, Taelman L, Moledina S, Taylor AM, Segers P, Coghlan GJ, Muthurangu V (2015) Noninvasive pulmonary artery wave intensity analysis in pulmonary hypertension. Am J Physiol Heart Circ Physiol 308(12):H1603–H1611 Qureshi M, Vaughan G, Sainsbury C, Johnson M, Peskin C, Olufsen M, Hill N (2014) Numerical simulation of blood flow and pressure drop in the pulmonary arterial and venous circulation. Biomech Model Mechanobiol 13:1137–1154. doi:10.1007/ s10237-014-0563-y Rain S, Handoko M, Trip P, Gan T, Westerhof N, Stienen G, Paulus W, Ottenheijm C, Marcus J, Dorfmuller P, Guignabert C, Humbert M, MacDonald P, dos Remedios C, Postmus P, Saripalli C, Hidalgo C, Granzier H, Vonk-Noordegraaf A, van der Velden J, de Man F (2013) Right ventricular diastolic impairment in patients with pulmonary arterial hypertension. Circulation 128(18):2016– 25. doi:10.1161/CIRCULATIONAHA.113.001873 Rausch M, Dam A, Göktepe S, Abilez O, Kuhl E (2011) Computational modeling of growth: systemic and pulmonary hypertension in the heart. Biomech Model Mechanobiol 10(6):799–811. doi:10.1007/ s10237-010-0275-x Reuben SR (1971) Compliance of the human pulmonary arterial system in disease. Circ Res 29(1):40–50. doi:10.1161/01.RES.29.1.40 Reymond P, Merenda F, Perren F, Rufenacht D, Stergiopulos N (2009) Validation of a one-dimensional model of the systemic arterial tree. Am J Physiol Heart Circ Physiol 297:H208–H222. doi:10.1152/ ajpheart.00037.2009 Rich J, Rich S (2014) Clinical diagnosis of pulmonary hypertension. Circulation 130(20):1820–1830. doi:10.1161/ CIRCULATIONAHA.114.006971 Sanz J, Kariisa M, Dellegrottaglie S, Prat-González S, Garcia M, Fuster V, Rajagopalan S (2009) Evaluation of pulmonary artery

123

S. Acosta et al. stiffness in pulmonary hypertension with cardiac magnetic resonance. JACC Cardiovasc Imaging 2(3):286–295. doi:10.1016/j. jcmg.2008.08.007 Sanz J, García-Alvarez A, Fernández-Friera L, Nair A, Mirelis J, Sawit S, Pinney S, Fuster V (2012) Right ventriculo-arterial coupling in pulmonary hypertension: a magnetic resonance study. Heart 98(3):238–43. doi:10.1136/heartjnl-2011-300462 Saouti N, Westerhof N, Helderman F, Marcus J, Boonstra A, Postmus P, Vonk-Noordegraaf A (2010a) Right ventricular oscillatory power is a constant fraction of total power irrespective of pulmonary artery pressure. Am J Respir Crit Care Med 182(10):1315–20. doi:10.1164/rccm.200910-1643OC Saouti N, Westerhof N, Postmus PE, Vonk-Noordegraaf A (2010b) The arterial load in pulmonary hypertension. Eur Respir Rev 19(117):197–203. doi:10.1183/09059180.00002210 Schuuring M, van Riel A, Vis J, Duffels M, Berger R, Hoendermis E, Van Dijk A, Vliegen H, Mulder B, Bouma B (2013) Resting heart rate is an important determinant of mortality in patients with pulmonary arterial hypertension due to congenital heart disease. J Am Coll Cardiol 61(10):E1271. doi:10.1016/S0735-1097(13)61271-1 Segers P, Stergiopulos N, Westerhof N (2000) Quantification of the contribution of cardiac and arterial remodeling to hypertension. Hypertension 36(5):760–765 Sherwin S, Formaggia L, Peiro J, Franke V (2003) Computational modelling of 1D blood flow with variable mechanical properties and its application to the simulation of wave propagation in the human arterial system. Int J Numer Methods Fluids 43:673–700 Stankovic Z, Allen B, Garcia J, Jarvis K, Markl M (2014) 4D flow imaging with MRI. Cardiovasc Diagn Ther 4(2):173–192. doi:10. 3978/j.issn.2223-3652.2014.01.02 Steelet B, Valdez-Jasso D, Haider M (2011) Predicting arterial flow and pressure dynamics using a 1D fluid dynamics model with a viscoelastic wall. SIAM J Appl Math 71(4):1123–1143 Stergiopulos N, Young D, Rogge T (1992) Computer simulation of arterial flow with applications to arterial and aortic stenoses. J Biomech 25(12):1477–1488 Streeter V, Keitzer W, Bohr D (1963) Pulsatile pressure and flow through distensible vessels. Circ Res 13(1):3–20 Strocchi M, Contarino C, Zhang Q, Bonmassari R, Toro E (2017) A global mathematical model for the simulation of stenoses and bypass placement in the human arterial system. Appl Math Comput 300:21–39. doi:10.1016/j.amc.2016.11.028 Tan W, Madhavan K, Hunter K, Park D, Stenmark K (2014) Vascular stiffening in pulmonary hypertension: cause or consequence? Pulm Circ 4(4):560–80. doi:10.1086/677370 Tang B, Pickard S, Chan F, Tsao P, Taylor C, Feinstein J (2012) Wall shear stress is decreased in the pulmonary arteries of patients with pulmonary arterial hypertension: an image-based, computational fluid dynamics study. Pulm Circ 2(4):470–6. doi:10.4103/ 2045-8932.105035

123

Thenappan T, Shah SJ, Rich S, Tian L, Archer SL, Gomberg-Maitland M (2010) Survival in pulmonary arterial hypertension: a reappraisal of the nih risk stratification equation. Eur Respir J 35(5):1079– 1087. doi:10.1183/09031936.00072709 Valdez-Jasso D, Haider M, Banks H, Santana D, Germán Y, Armentano R, Olufsen M (2009) Analysis of viscoelastic wall properties in ovine arteries. IEEE Trans Biomed Eng 56(2):210–219. doi:10. 1109/TBME.2008.2003093 Valdez-Jasso D, Bia D, Zócalo Y, Armentano R, Haider M, Olufsen M (2011) Linear and nonlinear viscoelastic modeling of aorta and carotid pressure-area dynamics under in vivo and ex vivo conditions. Ann Biomed Eng 39(5):1438–1456. doi:10.1007/ s10439-010-0236-7 Voelkel N, Quaife R, Leinwand L, Barst R, McGoon M, Meldrum D, Dupuis J, Long C, Rubin L, Smart F, Suzuki Y, Gladwin M, Denholm E, Gail D (2006) Right ventricular function and failure: report of a National Heart, Lung, and Blood Institute working group on cellular and molecular mechanisms of right heart failure. Circulation 114(17):1883–1891. doi:10.1161/CIRCULATIONAHA.106. 632208 Vonk-Noordegraaf A, Haddad F, Chin K, Forfia P, Kawut S, Lumens J, Naeije R, Newman J, Oudiz R, Provencher S, Torbicki A, Voelkel N, Hassoun P (2013) Right heart adaptation to pulmonary arterial hypertension: physiology and pathobiology. J Am Coll Cardiol 62(25):D22–D33. doi:10.1016/j.jacc.2013.10.027 Wang Z, Chesler N (2011) Pulmonary vascular wall stiffness: an important contributor to the increased right ventricular afterload with pulmonary hypertension. Pulm Circ 1(2):212–23. doi:10.4103/ 2045-8932.83453 Wang Z, Chesler NC (2013) Pulmonary vascular mechanics: important contributors to the increased right ventricular afterload of pulmonary hypertension. Exp Physiol 98(8):1267–1273. doi:10. 1113/expphysiol.2012.069096 Wang J, Parker K (2004) Wave propagation in a model of the arterial circulation. J Biomech 37(4):457–70. doi:10.1016/j.jbiomech.2003. 09.007 Wang Z, Schreier D, Hacker T, Chesler N (2013) Progressive right ventricular functional and structural changes in a mouse model of pulmonary arterial hypertension. Physiol Rep 1(7):e00184. doi:10. 1002/phy2.184 Xiao N, Alastruey J, Figueroa CA (2014) A systematic comparison between 1-D and 3-D hemodynamics in compliant arterial models. Int J Numer Methods Biomed Eng 30(2):204–31. doi:10.1002/ cnm.2598 Yang W, Feinstein J, Vignon-Clementel I (2016) Adaptive outflow boundary conditions improve post-operative predictions after repair of peripheral pulmonary artery stenosis. Biomech Model Mechanobiol 15(5):1345–1353. doi:10.1007/s10237-016-0766-5