Hemodynamic assessment of pulmonary hypertension in mice: A model based analysis of the disease mechanism
Abstract
This study uses a one dimensional fluid dynamics arterial network model to infer changes in hemodynamic quantities associated with pulmonary hypertension in mice. Data for this study include blood flow and pressure measurements from the main pulmonary artery for 7 control mice with normal pulmonary function and 5 hypertensive mice with hypoxia induced pulmonary hypertension. Arterial dimensions for a 21 vessel network are extracted from microCT images of lungs from a representative control and hypertensive mouse. Each vessel is represented by its length and radius. Fluid dynamic computations are done assuming that the flow is Newtonian, viscous, laminar, and has no swirl. The system of equations is closed by a constitutive equation relating pressure and area, using a linear model derived from stressstrain deformation in the circumferential direction assuming that the arterial walls are thin, and also an empirical nonlinear model. For each dataset, an inflow waveform is extracted from the data, and nominal parameters specifying the outflow boundary conditions are computed from mean values and characteristic time scales extracted from the data. The model is calibrated for each mouse by estimating parameters that minimize the least squares error between measured and computed waveforms. Optimized parameters are compared across the control and the hypertensive groups to characterize vascular remodeling with disease. Results show that pulmonary hypertension is associated with stiffer and less compliant proximal and distal vasculature with augmented wave reflections, and that elastic nonlinearities are insignificant in the hypertensive animal.
1 Introduction
Pulmonary hypertension (PH) is defined as an invasively measured mean pulmonary arterial blood pressure (mPAP) greater than 25 mmHg (Simonneau et al., 2013). It is associated with vascular remodeling, which adversely affects the properties of the cardiopulmonary system, including pulmonary vascular resistance (PVR), proximal and distal arterial stiffness, compliance, and amplitude of wave reflections (Nichols et al., 2011). The mPAP and PVR are conventionally used diagnostic markers for PH but are not good indicators of disease severity (Wang & Chesler, 2011). Here we use proximal arterial stiffness and the wave reflection amplitude as physiomarkers for detecting disease progression (Castelain et al., 2001; Hunter et al., 2011). In particular, the proximal arterial stiffness is an excellent predictor of mortality in patients with pulmonary arterial hypertension (Gan et al., 2007). Quantifying relative distributions of proximal and distal arterial stiffness (or compliance) and wave reflections in elevating the mPAP and PVR is vital for understanding disease mechanisms.
In this study, we setup and calibrate a mathematical model predicting wave propagation in the pulmonary vasculature in C57BL6/J male mice with normal pulmonary function (control group (CTL), n = 7) and in mice with hypoxiainduced pulmonary hypertension (hypertensive group (HPH), n = 5) (Tabima et al., 2012; Vanderpool et al., 2011). The novelty of this study is the integration of high fidelity morphometric and hemodynamic data from multiple mice with a one dimensional (1D) model of large pulmonary arteries coupled with a zero dimensional (0D) model of the vascular beds. This is achieved by incorporating available data at each stage of the modeling including network extraction, parameter estimation and model validation. The outcome is used to infer disease progression by quantifying relative changes in PVR, proximal and distal arterial stiffness, compliance, and amplitudes of wave reflections, across the two groups (CTL and HPH). Moreover, we investigate the influence of presumed elastic nonlinearities in the wall model on the parameter inference. This approach allows us to study disease impact on the distal vasculature by predicting waveforms from multiple locations in the simulated network, which are difficult to obtain experimentally.
Experimental studies.
To understand the relation between hemodynamics and vascular remodeling, it is essential to analyze morphometric and hemodynamic timeseries data during disease progression. Morphometric data can be obtained by non invasive procedures, like magnetic resonance imaging (MRI) or computed tomography (CT) scans (Meaney & Beddy, 2012), but dynamic pulmonary blood pressure can only be obtained from right heart catheterization (Rich et al., 2011). Moreover, in humans, the disease takes years to develop making it difficult to study its progression. An alternative means to gain understanding is to study disease progression in mouse models. An advantage is that mice have a relatively short lifespan and it is feasible to generate specific disease groups (e.g. mice with HPH) within a short timespan (1 month). Experimental studies in mice are typically done within a specific genetic strain, in order to limit variation among individuals. In most studies (e.g. Pursell et al. (2016); Tabima et al. (2012); Vanderpool et al. (2011)), hypoxia is used to induce pulmonary hypertension. This type of PH (Group III) is clinically manifested in human patients with hypoxia and lung disease (Simonneau et al., 2013), believed to be initiated by remodeling of the vascular beds (e.g microvascular vasoconstriction and rarefaction (Wang & Chesler, 2011)) followed by progressive remodeling of the large arteries (Tuder et al., 2007). Therefore, investigation of pulmonary hypertension in mice with hypoxia may provide vital understanding of disease mechanisms in humans with similar pathology.
Modeling studies.
Examples of previous modeling efforts include lumped 0D (Lankhaar et al., 2006; Lumens et al., 2009), distributed 1D (Acosta et al., 2017; Lungu et al., 2014; Qureshi et al., 2014) and locally complex three dimensional (3D) (Tang et al., 2012; Yang et al., 2016) models. Lankhaar et al. (2006) combined a 0D 3element Windkessel model with hemodynamic data from PH patients to quantify right ventricular afterload. Lumens et al. (2009) developed a geometric heart model coupled with a closed loop 0D model of the entire circulation to predict ventricular hypertrophy in patients with PH. Tang et al. (2012) and Yang et al. (2016) constructed patient specific 3D models of pulmonary arteries to study shear stress (Tang et al., 2012), and PVR in pre and post operative situations (Yang et al., 2016). To our knowledge, only Lungu et al. (2014) used a coupled 1D0D model of the main pulmonary artery (MPA) to study diagnostic parameters including PVR, mPAP, arterial stiffness and compliance in patients with and without PH. Studies by Acosta et al. (2017) and Qureshi et al. (2014) used a 1D framework to develop distributed models of the pulmonary arteries and veins to study cardiopulmonary coupling (Acosta et al., 2017) and microvascular remodeling (Qureshi et al., 2014) during PH. Yet, none of these studies used actual pressure data for parameter estimation and model validation. Other notable studies using 1D models, but not investigating pulmonary hypertension, include Blanco et al. (2014); Mynard & Smolich (2015); Olufsen et al. (2000); Reymond et al. (2009); Willemet & Alastruey (2015). See Boileau et al. (2015); Safaei et al. (2016); van de Vosse & Stergiopulos (2011) for an overview of modeling approaches and computational tools, and Hunter et al. (2011); Kheyfets et al. (2013); Tawhai et al. (2011) for focused reviews on modeling pulmonary hypertension.
Most of the above studies consider application to humans. Only a few of 1D modeling studies have investigated wave propagation in systemic (Aslanidou et al., 2016) and pulmonary (Lee et al., 2016; Qureshi et al., 2017) arteries in mouse networks. To our knowledge, no studies combine disease specific imaging and hemodynamics data to predict changes in disease. In this study, we expand upon prior results from Qureshi et al. (2017) by developing a 1D fluid dynamics network model of pulse wave propagation in the large pulmonary arteries of control and hypertensive mice. We do so by extracting detailed network information and by optimizing hemodynamic predictions. To set up the 1D domain, we extract the arterial networks from microCT images of a representative CTL mouse and a mouse with HPH. We combine these networks with hemodynamic data, measured in the MPA of each mouse, to predict pressure and flow waveforms. For each vessel in the network, we solve (numerically) fluid dynamics equations derived from the Navier–Stokes equations coupled with a constitutive equation relating pressure and vessel area (i.e. describing the vessel wall mechanics). We use the measured flow waveform from each mouse as the inflow boundary condition and use 0D Windkessel models (Westerhof et al., 2009) to characterize the impedance of the vascular beds.
As the disease progresses, vascular remodeling changes the composition of wall constituents (Humphrey, 2008). This affects the local stressstrain response, which is important for inferring arterial wall stiffness and compliance (Hunter et al., 2011). We focus on understanding how wall deformation changes under hypertensive conditions. Specifically, we aim to study if presumed elastic nonlinearities have a significant impact on disease progression. To this end we use a well established linear wall model (Willemet & Alastruey, 2015; Safaei et al., 2016) as well as an empirical nonlinear wall model.
The advantage of the linear wall model is that it is easy to derive from first principles and it has been shown to be successful within physiological pressure/area values for systemic arteries (Olufsen et al., 2000). Yet, it does not account for the fact that arteries stiffen with pressure. On the other hand, detailed, structural hyperelastic wall models can be derived from first principles, e.g. Gerhard & Ray (2010), but they are difficult to analyze due to the large number of parameters. Inspired by Langewouters et al. (1985), we introduce a simple empirical nonlinear model with three key properties: (a) it predicts vessel stiffening with pressure so that the lumen area approaches a finite limit as pressure increases, (b) it incorporates an undeformed area (or radius) corresponding to zero transmural pressure, and (c) it reduces to the linear wall model under a small strain assumption providing a basis for model comparison and nominal parameter estimation.
The process of modeling requires a priori specification of parameters, including proximal arterial stiffness, total vascular resistance, and peripheral compliance, which are known to vary across individuals. This creates the need for methods to estimate parameters that predict observed hemodynamics across individual mice. A few recent studies have performed parameter estimation in blood flow models, including the study by Eck et al. (2017), which used polynomial chaos expansion to analyze a stochastic model of pressure waves in the large systemic arteries, the study by Arnold et al. (2017), which used an ensemble Kalman filter (EnKF) to estimate the unknown inflow to a single vessel, representing the ovine aorta, and the study by Tran et al. (2017), which used a Bayesian Markov chain Monte Carlo (MCMC) approach to estimate parameters in a multiscale three dimensional model of coronary arteries. Similarly, in the recent study by Paun et al. (2018), we used MCMC to estimate parameters for a 1D model of mouse pulmonary arteries. It should be noted that MCMC algorithms come at a substantial computational cost, making them infeasible for multisubject studies. Other studies estimating pressure dynamics using optimization algorithms were done using either 0D (Valdez et al., 2011; Williams et al., 2014) or 1D (Lungu et al., 2014) models. Yet, none of these studies investigated hemodynamic variation in the pulmonary arterial network.
In this study, we estimate global network parameters that allow prediction of observed dynamics in both CTL and HPH groups. We first determine a priori parameter values for the wall models and boundary conditions. This is done by combining available data and existing results in the literature (Alastruey et al., 2016; Reymond et al., 2009). Second, we solve the overall model using linear and nonlinear wall models and conduct constrained nonlinear optimization to estimate parameters predicting the observed dynamics. A similar approach combining apriori nominal parameters with iterative tuning was used by Alastruey et al. (2016) in a model of the human aorta. In this study blood flow data were available for all terminal vessels, making it easier to compute downstream resistance. To our knowledge, our approach estimating parameters for a 1D network model of the pulmonary circuation using morphometric and dynamic pressure and flow data is novel relative to these studies.
Similar to previous studies (Lee et al., 2016; Qureshi et al., 2014, 2017), we compare the hemodynamic signatures in the time and frequency domains for the control and hypertensive animals. We follow Acosta et al. (2017); Lankhaar et al. (2006); Lungu et al. (2014) to analyze the estimated parameters, inferring HPH progression, and to investigate, using a model selection criterion (Burnham & Anderson, 2001; Schwarz, 1978), the extent to which the nonlinear wall model enables accurate prediction of the observed dynamics.
The manuscript is organized as follows: Sec. 2 presents experimental and mathematical methods including data extraction procedures, governing equations, parameter estimation and numerical simulations. In Sec. 3 we present results comparing CTL and HPH hemodynamics, analyzing: waveforms predictions along the arterial network, estimated parameters and their sensitivities to hypertensive conditions, and impedance and wave reflection analysis. Key findings are discussed in Sec. 4, followed by limitations in Sec. 4.4. Finally, we state key conclusions in Sec. 5.
2 Methods
2.1 Experimental Methods
This study uses existing hemodynamic data and micro computed tomography (microCT) images from control and hypertensive mice. Detailed experimental protocols for extracting the hemodynamic and image data can be found in Tabima et al. (2012) and Vanderpool et al. (2011), respectively. Both procedures were approved by the University of Wisconsin Institutional Animal Care and Use Committee. Below we summarize these protocols and highlight the data analyzed herein.
Hemodynamic data.
The hemodynamic data include dynamic pressure and flow waveforms from male C57BL6/J mice, average age 1213 weeks and average body weight of 24 g. The mice were divided into CTL (n = 7) and HPH groups (n = 5). The mice in the HPH group were exposed to 21 days of chronic hypoxia (10% O partial pressure) and both groups were exposed to a 12 hour lightdark cycle. Mice were instrumented to obtain dynamic pressure and flow waveforms in the main pulmonary artery as described previously (Tabima et al., 2012). In brief, mice were anesthetized with an intraperitoneal injection of urethane solution (2mg/g body weight) and placed on a heating pad to maintain physiological heart rate. After intubation for ventilation, the chest wall was removed to expose the right ventricle. A 1.0F pressuretip catheter (Millar Instruments, Houston, TX) was inserted into the apex of right ventricle and advanced to the MPA. Flow was measured with ultrasound (Visualsonics, Toronto, Ontario, Canada) using a 40 MHz probe during catheterization and recorded synchronously with pressure and ECG (Cardiovascular Engineering, Norwood, MA). Pressure and flow waveforms were signalaveraged using the ECG as a fiducial point. Twenty consecutive cardiac cycles were averaged to produce average pressure and average flow waveforms. Representative hemodynamic data and associated frequency domain signatures are shown in Figure 1 and essential cardiovascular parameters are summarized in Table 1.
Healthy (n = 7)  Hypoxia (n = 5)  
HR (beats/min)  53327  55921 
CO (ml/s)  0.180.03  0.150.01 
mPAP (mmHg)  13.43.1  22.12.3 
sPAP (mmHg)  20.13.5  32.33.1 
dPAP (mmHg)  8.03.3  14.02.3 
pPAP (mmHg)  12.21.4  18.32.8 
PVR (mmHg.s/ml)  77.022.7  146.123.6 
Abbreviations: Heart rate (HR), cardiac output (CO), mean pressure (mPAP), systolic pressure (sPAP), diastolic pressure (dPAP), pulse pressure (pPAP) in the main pulmonary artery, total pulmonary vascular resistance (PVR).
Imaging data.
Stacked planar Xray microCT images of pulmonary arterial trees were obtained from male C57BL6/J mice, average age 1012 weeks, under control (healthy mice) and 10 days of hypobaric hypoxia at 10% O partial pressure protocols. Detailed descriptions of animal handling and lung preparation can be found in Vanderpool et al. (2011), whereas details of the microCT image acquisition are described in Karau et al. (2001). In brief, mice were anesthetized with intraperitoneal injection of pentobarbital sodium (52 mg/kg body weight) and then euthanized by exsanguination. The trachea and the main pulmonary artery were cannulated and the heart was dissected away. Pulmonary artery cannula (PE90 tubing, 1.27mm external and 0.86mm internal diameter) was positioned well above the first bifurcation. Excised lungs were treated with Rho kinase inhibitor and ventilated and perfused with perfluorooctyl bromide (PFOB), a vascular contrast agent, and placed in the imaging chamber. The arterial trees were then imaged under a static filling pressure of 6.3 mmHg, while rotating the lungs in the Xray beam at 1 increments to obtain 360 planar images. Each planar image was averaged over seven frames to minimize noise and maximize vascular contrast. Fo each lung, isometric 3D volumetric dataset (497497497 pixels) was obtained, by reconstructing the 360 planer images using the Feldkamp conebeam algorithm (Feldkamp et al., 1984), and converted into Dicom 3.0.
For this study, two representative networks with 21 vessels were extracted from images of the control and hypertensive mice. The 21 vessel network was chosen since it was the most expansive network that could be identified with a onetoone vessel map in both control and hypertensive animals. Network dimensions and connectivities were obtained using the segmentation protocol described by Ellwein et al. (2015). This protocol uses ITKSNAP (Yushkevich et al., 2006) to create a 3D geometry from Dicom 3.0 files, using semiautomated “snake evolution” in the regions of interest (the 21 vessels). The image pixel threshold was set at 45 to reduce artifact detection. Paraview (Kitware; Clifton Park, NY) was used to convert file types to vtk polygonal data (.vtp) allowing us to compute centerlines and connectivity using the Vascular Modeling ToolKit (VMTK (Antiga et al., 2008)). The output from VMTK is a matrix representing each vessel by a unique set of coordinates and the associated radius value, , computed from the maximally inscribed sphere within the 3D vessel. Known internal diameter of the contrast filled cannula (PE90 tubing) was used to compute a scaling factor to convert voxels into cm. For the MPA, the radius was computed as a mean of all slices between the cannulated region and the first bifurcation, whereas the radii of other vessels were computed as the mean of all slices along the vessels but away from the junctions. For each vessel, the length, , was calculated as the sum of the shortest distances () between successive points, i.e.
(1) 
and is the total number of samples for the vessel. The network, in its load free state (i.e. at zero transmural pressure) is then constructed by assigning and to each vessel. The shared coordinates information is used to generate a connectivity map of the 3D structure (see Figure 2(c)) using the “digraph” function in MATLAB (version 16a). Both control and hypertensive networks have the same connectivity illustrated in Figure 2(c), but the individual vessel radii and length vary as shown in Table 2.
Control  Hypertensive  
vessel  connectivity  
index  (daughters)  (cm)  (cm)  (cm)  (cm) 
1  (2,3)  0.47  4.10  0.51  3.58 
2  (4,5)  0.26  4.45  0.26  4.03 
3  (6,7)  0.37  3.72  0.37  3.08 
4  (8,9)  0.24  2.41  0.25  2.92 
5  –  0.13  0.52  0.17  0.65 
6  (14,15)  0.32  2.02  0.28  1.60 
7  –  0.17  2.12  0.19  0.93 
8  (10,11)  0.23  3.11  0.24  2.06 
9  –  0.17  1.77  0.17  0.51 
10  (12,13)  0.20  2.62  0.22  2.37 
11  –  0.16  0.69  0.17  0.88 
12  –  0.15  1.40  0.19  1.27 
13  –  0.14  0.62  0.15  0.51 
14  (16,17)  0.26  0.81  0.27  1.20 
15  –  0.19  1.84  0.19  1.55 
16  (18,19)  0.25  0.83  0.26  0.71 
17  –  0.15  3.02  0.18  1.68 
18  (20,21)  0.24  4.69  0.24  3.55 
19  –  0.15  1.77  0.18  1.86 
20  –  0.22  1.78  0.23  2.24 
21  –  0.18  0.55  0.19  1.07 
* root vessel. Connectivity , denotes the left daughter and the right daughter. Vessels indicated by – are terminal.
2.2 Fluid Dynamics Model
Assuming that blood is incompressible, flow is Newtonian, laminar and axisymmetric, and has no swirl, conservation of mass and momentum (Olufsen et al., 2000) are given by
(2) 
where and are the axial and temporal coordinates, (mmHg) is the transmural blood pressure, and are the pressures acting in and outside the arterial wall, respectively, (ml/s) is the volumetric flow rate, (cm) is the crosssectional area and (cm) is the vessel radius. The blood density (g/ml) and the kinematic viscosity (cm/s) are assumed constant. The momentum equation is derived under the noslip condition assuming that the wall is impermeable, and that the velocity of the fluid at the wall equals the velocity of the wall. To satisfy this condition, we assume that the velocity profile over the lumen area is flat, decreasing linearly within the boundary layer with thickness (van de Vosse & Stergiopulos, 2011).
2.3 Wall Model
To close the system of equations, a constitutive equation relating pressure and crosssectional area is needed. In this study, we compare two models. A linear elastic model (Safaei et al., 2016) derived from balancing circumferential stress and strain, and an empirical nonlinear wall model inspired by Langewouters et al. (1985).
The linear wall model
is derived under the assumptions that the vessels are cylindrical and purely elastic, that the walls are thin (), incompressible and homogeneous, that the loading and deformation are axisymmetric, and that the vessels are tethered in the longitudinal direction. Under these conditions, the external forces can be reduced to stresses in the circumferential direction, and Laplace’s law (Nichols et al., 2011) yields the linear stressstrain relation
(3) 
is the stiffness parameter, defined in terms of Young’s modulus in the circumferential direction. The associated Poisson ratio , the wall thickness , and the undeformed radius at zero transmural pressure (). (cm) denotes the undeformed cross sectional area. Similar to previous studies (e.g. Olufsen et al. (2000); Safaei et al. (2016)), we use .
The nonlinear wall model
relates pressure and area as
(4) 
where (mmHg) is a material parameter that describes the halfwidth pressure, and is a scaling parameter that determines the maximal lumen area as , giving
(5) 
This model is formulated to ensure that at .
2.4 Boundary Conditions
Since the system of equations is hyperbolic, boundary conditions must be specified at the inlet and outlet of each vessel, i.e. the network needs an inflow condition, junction conditions, and outflow conditions.
At the network inlet we specify a flow waveform extracted from hemodynamic data (see Figure 1). At junctions (all bifurcations in the network studied), we impose pressure continuity and conservation of flow, i.e.
(6) 
where, the subscripts and () denote the parent and daughter vessels, respectively.
At the terminal vessels, a Windkessel model (represented by an circuit) is used to prescribe the outflow boundary condition (see Fig. 2c) by computing input impedance as
(7) 
where and are the pressure and flow in the frequency domain, is the angular frequency and is the length of the cardiac cycle. (mmHg s/ml) denote the two resistances, and (ml/mmHg) the capacitance. Moreover, the total peripheral resistance is , where represents the resistance of the proximal vasculature and is the resistance of the distal vasculature, whereas denotes the total peripheral compliance of the vascular region in question. Similar to Qureshi et al. (2017), relates the pressure and flow at the outlet of each terminal vessel via a convolution integral over .
2.5 Parameter Values
The model parameters are divided into three groups: hemodynamics , vessel wall stiffness for the linear wall model and for the nonlinear wall model, Windkessel for the vascular beds.
Hemodynamic parameters
are assumed constant. For each mouse (control and hypertensive), the length of the cardiac cycle 1/HR (s) is extracted from the data (mean HR values for each group are given in Table 1). The blood density g/ml (Riches et al., 1973) and the kinematic viscosity cm/s, measured at a shear rate of 94 s (Windberger et al., 2003). These values represent average values for mice. As discussed earlier the boundary layer thickness is .
The nominal vessel wall parameters
for the linear and nonlinear wall models are approximated by combining an analytical approach, involving linearization of the fluid dynamics model, with available pressure and flow data. Since the hemodynamic data is available from one location only, we assume a constant throughout the network.
Linearization of (2) and (3) about the reference state provides an expression for the characteristic impedance (Nichols et al., 2011) as
(8) 
where is the wave speed at (see Appendix B for the definition of for the linear () and the nonlinear () wall models). For each mouse, is esitmated from the slope of the pressureflow loop during early ejection, using the ‘upslope method’ (see Figure 3a) (Dujardin, & Stone, 1981) including 95% of the flow during ejection phase.
For the linear wall model, substituting in (8) with gives
(9) 
Similarly, for the nonlinear wall model, substituting in (8) with and using (9), gives
(10) 
where is given by (9) and is the ‘stiffness’ of the nonlinear wall model. To fully specify the nominal (initial) values for the parameter inference (see Sec. 2.6 and Table 4), we set and . These values give cm. Note that (10) is only valid if both the linear and nonlinear models incorporate the same value of at .
Nominal parameters for outflow boundary conditions
must be computed for each terminal vessel. A priori values for these parameters can be obtained from distributing the total peripheral resistance computed in the MPA (, the ratio of the mean pressure gradient to mean flow) to each terminal vessel as
where is the mean flow in vessel . This expression is valid, under the assumption that capillary pressure drops to zero and that mean pressure remains constant across the arterial network. The flow to vessel is estimated by applying Poiseuille’s law recursively at each junction, giving
(11) 
is the vessel conductance. Here denotes the mean flow to vessel . Similar to previous studies, e.g. Reymond et al. (2009), the total resistance is distributed as and , where the a priori value of (McDonald & Attinger, ).
Finally, as suggested by Stergiopulos et al. (1995), the total vascular compliance (defined in Appendix A) can be calculated by estimating the timeconstant by fitting data to an exponential function of the form
(12) 
where the diastolic pressure decays in time (see Figure 3b). Assuming that is constant throughout the network, the compliance for each Windkessel can be computed as
(13) 
2.6 Numerical Simulations
Similar to previous studies (Olufsen et al., 2000, 2012), we use the twostep Lax Wendroff method to solve the model equations in Sections 2.2 and 2.3, for the arterial networks presented in Table 2. We combine data from Sec. 2.1 with methods described in Sec. 2.5 to estimate nominal parameters, specifying arterial wall stiffness and boundary conditions (Sec. 2.4). Convergence to steady state is ensured by running simulations until the least square error between consecutive cycles of pressure is less than , taking about 6 cycles on average. For this study, we fixed the number of cycles to 7 for all mice.
To match the model to data, parameters inferred using optimization include for the linear model and for the nonlinear model. We define the objective function using sum of squared errors as
(14) 
where is the length of the time vector spanning one cardiac cycle, and are the measured pressure and flow, and are the computed pressure and flow from the inlet of the MPA, and are the parameters to be optimized using the linear or the nonlinear wall models. We use the function fmincon in MATLAB under the Sequential Quadratic Programming (SQP) gradientbased method (Boggs & Tolle, 2000) to solve the associated leastsquares estimation problem
Optimization was run in parallel on an iMac (3.1 GHz Intel Core i5, 16GB RAM, OS 10.12.6) initializing the parameters with four initial values including the nominal and nearby values. Nominal values for individual mice are given in Table 4 (Appendix C), whereas upper and lower bounds as well as averaged optimal parameters are given in Table 5 and 6 in Appendix D. The algorithm was iterated until the convergence criterion was satisfied with a tolerance .
2.7 Model Analysis
Similar to Qureshi & Hill (2015) and Qureshi et al. (2017), we analyze wave intensity and impedance spectra for further insight. Wave intensity analysis allows us to quantify the type and nature of the reflected waves in the time domain whereas the impedance analysis provides a frequency domain signature, which as shown in Figure 1 differs between the two groups.
Wave intensity analysis
allows us to separate the simulated waveforms into their incident () and reflected () components assuming negligible frictional losses. By setting , where is the fluid velocity, the incident and reflected waves can be approximated by
(15)  
and is the PWV computed by (31) (Appendix B). Timenormalized wave intensity (W/cms) is defined as
(16) 
along with or characterize the incident waves as compressive or decompressive, while and or characterize the reflected waves as compressive or decompressive, respectively. Finally, we compute the wave reflection coefficient as the ratio of the amplitudes of the reflected to the incident pressure waves, as
(17) 
Impedance analysis (IA).
Under the assumptions of periodicity and linearity, the pulsatile pressure and flow waveforms can be approximated by a Fourier series of the form
(18) 
where is the Fourier series approximation of the original waveform , is the time vector for a given sampling rate Fs, is the length of the signal , ) are the angular frequencies, is the mean of , and and (rad) are the moduli and phase spectra, associated with each harmonic , and is the smallest resolution of harmonics required for the impedance analysis. Both, and , are defined in terms of and , the coefficients of basic trigonometric Fourier series, i.e.
Setting as and in (18), the impedance spectrum Z, can be computed as ratios of harmonics of pressure to flow by
(19)  
where and are the moduli and and the phase angles of the pressure and flow harmonics, respectively. s are the impedance moduli and the corresponding phases at a given frequency. Note if then the pressure harmonic lags the flow harmonic, and vice versa. The zeroth harmonic is the total pulmonary vascular resistance (PVR).
2.8 Statistical analysis
We implement statistical analysis methods to study parameter interference and devise a model selection criteria to determine the ability of the linear and nonlinear wall models to predict hemodynamics.
Parameter inference.
Optimized parameters and hemodynamic quantities are compared to assess changes with hypertension. For this analysis we compare predictions of total vascular resistance (), total vascular compliance (), the compliance ratio (), the resistance ratio (), the wave reflection coefficient (), and characteristic time scales (). Impact of the disease on a given quantity , averaged across the two groups, is inferred by computing an importance index as a relative change in due to HPH as
(20) 
Model selection criterion.
The 1D fluid dynamics model is coupled with a linear and a nonlinear wall model, leading to 4 dimensional (4D) and 5 dimensional (5D) parameter spaces, respectively. To identify the model that is more consistent with the data, we employ a statistical criterion that trades off goodness of fit versus model complexity (i.e. the number of estimated parameters). To this end, we use the corrected Akaike Information Criterion (AICc) (Burnham & Anderson, 2001) and the Bayesian Information Criterion (BIC) (Schwarz, 1978) defined in (2.8) and (22). The model with the lower AICc and BIC score is preferred.
AICc  
BIC  (22) 
Here is the maximum loglikelihood, is the number of parameters in the model and is the total number of measurements. If the noise is white, i.e. independent and identically Gaussian distributed, then is just a rescaled version of the sumofsquares function , defined in (14). However, a plot of the residuals (see Figure 10 in the Results section) indicates that the independence assumption is violated and that the correlation structure of the noise needs to be taken into account. Under the assumption of multivariate normal noise, the log likelihood is given by
(23) 
where is the covariance matrix and is the vector of residuals. For the estimation of , we tried different approaches. We fitted an autoregressive moving average model, ARMA (p,q), to the time series of residuals, and identified the optimal parameters p and q by minimising the BIC score. We then used the standard procedure proposed by Box & Jenkins (Box & Jenkins, 1970) to estimate the covariance matrix. Alternatively, we fitted a Gaussian process (GP) to the time series, using a variety of standard kernels: squared exponential, Matérn 5/2, Matérn 3/2, neural network kernel and a periodic kernel; see Rasmussen & Williams (2006) for details.
For the BIC and AICc scores, we need the maximum likelihood configuration of the parameters. This would require an iterative optimization scheme, where for each parameter adaptation, the covariance matrix would have to be recomputed. As this would lead to a substantial increase in the computational costs, we approximated the maximum likelihood parameters by the parameters that minimize the residual sumofsquares error of equation (14).
Obtaining requires inverting the covariance matrix, see (23). To avoid numerical instabilities in the inversion of the covariance matrix, we introduce a regularization scheme: . This adds jitter to the diagonal elements for improving the conditioning number of the matrix while leaving the marginal variance (i.e. the diagonal elements of ) invariant. To minimize the bias, the interpolation parameter is chosen to be the minimum value subject to the constraint that the inversion of is numerically stable, as measured by a conditioning number below . We discarded covariance matrices that required values of above 0.1, due to the fact that in these cases the numerical stabilization incurred too high a bias.
3 Results
In this section, we present results of numerical simulations predicting pressure and flow waveforms along the pulmonary arterial network for a representative control and hypertensive mouse followed by a comparison of estimated parameter values.
3.1 Hemodynamics
Figure 4 shows optimized pressure and flow waveforms for a representative control (left) and hypertensive (right) animal using the linear (dashed) and nonlinear (solid) wall models. Results were obtained estimating the least squares error between measured and computed waveforms in the MPA. The box and whisker plot (Figure 4A) summarizes the least squares errors (reported in Table 6 in Appendix D) comparing the measured and computed waveforms in the MPA (vessel 1). Results show that both wall models (linear and nonlinear) are able to fit the data well, and the least squares error is smaller for the hypertensive animals.
Even though both models provide similar fits in the MPA, the two wall models lead to different predictions in the downstream vasculature. For the control mouse, the nonlinear wall model leads to higher pressure in the downstream vasculature for the control mouse. With the linear wall model, the mean pressure drop () mmHg from the MPA to the terminal vessels, whereas for the nonlinear model mmHg. Moreover, the nonlinear wall model predicts a notch in the pressure and flow waveforms during the ejection phase, not observed with the linear wall model (displayed in panels 2, 8 and 13). Finally it should be noted that these differences cannot be observed in predictions from the hypertensive animal, likely since the wall is almost rigid.
Figure 5, depicting the pressurearea relationship in the MPA, gives further insight into the behavior of the linear and nonlinear wall mechanics. As expected, results show that for both groups the nonlinear wall model yields area predictions that are concave down implying increased stiffening at higher pressures. Higher compliance in control animals lead to more significant deformation despite the observation that the unstressed radius is larger in the hypoxic animals ( vs. cm, from Table 2).
3.2 Parameter estimates
Figures 6 summarizes the variation in estimated parameters for the two groups. Parameters, used to to generate Figures 6, are given in Table 6 in the Appendix.
The arterial wall stiffness
is larger in mice with HPH (a–c). For the linear model, this claim is supported by larger estimates of the stiffness parameter (a). For the nonlinear wall model increased stiffness is predicted by a statistically significant increase in the parameter representing the half width pressure (b). This is not surprising since and have the same units and are proportional (see eq. 10). In addition, the constant (estimated using the nonlinear wall model) is smaller leading to a smaller maximal area (see eq. (5)).
Compliance distributions
between the network () and the vascular beds () are shown in Figure 6 (e–f). Here the total compliance () for the entire vasculature is computed as
where is determined by recursively adding local compliance estimates (eqs. (25) and (26) in Appendix A), computed as functions of the diastolic blood pressure and area for each vessel , while denote the weighted Windkessel capacitance (eq. (7)) associated with each terminal vessel . For simplicity, we omit the symbol in the rest of the manuscript.
Results show that the total vascular compliance decreases with hypertension (Figure 6e). Given the reciprocal relation between compliance and stiffness (eqs. (25) and (26)), this behavior is anticipated, shown in Figure 6(a)–(c). The nonlinear model predicts a smaller compliance for both groups due to an increased stiffening with pressures. Finally, Figure 6f shows the compliance distribution via the ratio . Although the majority of the compliance is located in the vascular beds (i.e. )(f), it decreases along the vasculature (from larger to smaller vessels). Finally, results with the linear wall model reveal that reduces to 95% of under hypertensive conditions, indicating that large vessels have been remodeled by the disease.
The total vascular resistance
is computed as
where is the terminal vessel index. Its distribution within the vascular beds is depicted in Figure 6 (d,h). In this study, we estimated the proximal and distal components and as described in Sec. 2.6. Results show that for both wall models increases under hypertensive conditions (d). This increase is dominated by an increase in , evidenced by the increase in the resistance ratio (h).
The Characteristic time scale
, decreases under hypertensive conditions (Figure 6g). However, the nonlinear model predicts a smaller for all animals.
3.3 Wave intensity analysis.
Results of the wave intensity analysis (WIA) allow us to investigate the behavior of the incident and reflected waves. Since WIA requires measurement of dynamic crosssectional area and PWV, this analysis is done on the simulated waveforms.
Figure 7 shows separation of the pressure waveforms (a,b) into their incident and reflected components and the associated wave intensities (eh). The group averaged reflection coefficients , quantifying the magnitude of the reflections, are shown in panels (c) and (d). Results show that the amplitudes of the reflected waves are bigger in the hypertensive mice (b) leading to a bigger wave reflection coefficient, shown in panel (d). The effects of elastic nonlinearities on the amplitudes of the incident and reflected waves are minor (the increase in is negligible) (c,d).
Figure 7(e)(h) present the wave intensity profiles associated with the separated pressure and velocity waves. The inner panels show zoomed intensities around peak systole. Again, the amplitude of the incident wave is higher for hypertensive mice (compare panels (e,g) with (f,h)). For both groups the nonlinear wall model gives rise to a smaller amplitude incident wave (compare panels (e,f) to (g,h)). Similar patterns are observed for the reflected waves. One notable difference, for both groups, is that the nonlinear model gives rise to more oscillations of the reflected wave.
3.4 Impedance analysis
Figure 8 depicts the impedance moduli and phase spectra computed from the measured and simulated waveforms. Dashed lines show simulation results and solid black lines show data. Panels (a,c,e,g) show results from a control mouse and panels (b,d,f,h) show results from a hypertensive mouse. The impedance spectra were generated using (19) and plotted for the first 14 harmonics including the mean component (zeroth harmonic).
While timevarying simulations fit the data well, Figure 4 shows characteristic differences in frequency domain signatures between the linear and nonlinear wall models. First, the zeroth frequency components do not vary between models. Comparison of moduli spectra (ad) show that the linear wall model better captures the frequency response of the original system. However, both models miss the spike in the impedance moduli at the 3rd harmonic (about 30 Hz) observed for the control mouse (a,c). This should be contrasted with results from the hypertensive animal, where both models predict the low frequency behavior well (b,f). At higher frequencies, particularly after the 9th harmonic, the nonlinear wall model deviates from the measured impedance. In addition, the associated phase () dips below zero indicating persistence of pressure harmonics, which recede the flow harmonic. Again, the linear wall model deviates less at higher frequencies and its phase oscillates about zero.
3.5 Statistical analysis
In this section we compare estimated hemodynamics quantities pertinent to analysis of disease progression in HPH mice. To do so we calculate an importance index () computed from (20) using quantities averaged across the CTL and HPH groups.
Figure 9(a) shows predictions of for the essential cardiovascular quantities (summarized in Table 1), whereas Figure 9(b) does the same for the model parameters. Positive or negative values indicate an increase or a decrease in the quantity due to HPH. Moreover, a value of 1 (or 1) denote a 100% change in the quantity.
Results show that peripheral vascular resistance () and compliance () significantly contribute to differentiating between control and disease for the hypertensive animals (), whereas the resistance ratio is important for predictions with the linear model (), but not for predictions with the nonlinear model . A similar observation was made for the nonlinear stiffness parameters, that is more important than the linear stiffness parameter . Though both of these contributed significantly to distinguishing control vs. hypertension. Finally, the wave reflection coefficient and the time constant also differ between control and disease but at a smaller scale ().
Finally, we use statistical model selection criteria to determine which wall model is more consistent with the available hemodynamic data. Since the simple least squares error shows no difference between the two wall models for the hypertensive mouse, the model selection was done for the representative control mouse only. Specifically, we employ two model selection criteria (AICc and BIC), and model the residual correlation using the statistical ARMA model, and the GP models with different kernels (squared exponential, Matérn 3/2, Matérn 5/2, neural network and periodic kernel), as described in Sec. 2.8. We estimated the hyperparameters of the GP covariance matrices from the time series of residuals plotted in Figure 10 by maximum likelihood. This was done by either using standard optimization algorithms involving maximizing the combined likelihood for the linear and the nonlinear models (Rasmussen & Williams, 2006), or by separately finding the hyperparameters for the linear and the nonlinear models and then averaging the covariance matrices. Since the correlation structure of the noise depends on the experimental protocol and is independent of the model assumptions, the later approach can be implemented under the constraint that the covariance matrices are the same for both wall models. We found that both approaches yield similar results to optimizing the GP hyperparameters for the residuals of both wall models. For simplicity, we applied the latter procedure to the ARMA model as well.
Score type  Wall model  GP Matérn 3/2  GP Matérn 5/2  GP neural network 

AICc  Linear  621  693  652 
AICc  Nonlinear  607  681  641 
BIC  Linear  604  676  635 
BIC  Nonlinear  586  660  620 
We found that for three covariance matrices (GPs with squared exponential and periodic kernels and the ARMA model), the amount of regularization required for numerical stabilization was high (). Therefore, we discarded these results due to the high expected bias (although the results obtained with the ARMA model were consistent with our findings). We obtained the AICc and BIC scores for the covariance matrices that only required a low amount of numerical stabilization, (GP with Matérn 5/2, Matérn 3/2 and neural network kernel), and found that they were consistent with each other indicating lower AICc and BIC scores for the linear model compared to the nonlinear model, see Table 3.
4 Discussion
In this section we discuss our findings in physiological (CTL vs. HPH) and modeling (Linear vs. Nonlinear wall mechanics) contexts.
4.1 Inference of Disease Progression: CTL vs. HPH
In this study we used microCT images to set up anatomical networks and combined them with hemodynamic predictions from a 1D model. We calibrated the model for seven control and five hypertensive mice, estimating key parameters (large vessel stiffness, peripheral vascular resistance, peripheral vascular compliance, wave reflection coefficient) that characterize vascular remodeling due to pulmonary hypertension.
Stiffness, compliance, and resistance
Comparison between the two groups (independent of the wall model) shows that hypertensive mice exhibit stiffer vessels both proximally (within the network), and in the microcirculation (represented by the Windkessel models) reflecting a rather advanced stage of HPH. In this study, we attribute changes in total resistance and compliance to both functional (vasoconstriction) and structural (rarefaction) microvascular remodeling, known to elevate the mPAP (Qureshi et al., 2014; Wang & Chesler, 2011). Changes in the large vessel parameters, and for the nonlinear model and for the linear model as well as the unstressed vessel radius , contribute to the stiffness of the proximal vessels.
The results presented here show that the parameter with the largest impact is . This is expected physiologically, as the disease progresses from the microvasculature to the large vessels. Almost as significant are changes in vessel stiffness ( and ) reflecting that the large vessels also stiffen, but with a relatively small change in indicating that the large vessels have not decreased in size. These results also reveal a significant decrease in compliance with hypertension.
In summary, these findings suggest that the hypertensive animals analyzed here exhibit remodeling of the entire vasculature, with distal vascular remodeling playing a dominant role in elevating the larger arterial blood pressure.
Wave reflection.
Analysis of computed waveforms show that the wave reflection coefficient is significantly higher in hypertension (Figure 7 (c,d) and 9). This observation can be explained by hypothesized proximal and distal remodeling. As shown in Figure 7 (a,c,e,g) the system has some degree of compressive wave reflections under control conditions. Our results show that the coefficient is amplified in hypertension (Figure 7 d) as the remodeled (stiffened) proximal arteries require a higher ejection pressure from the heart to maintain the cardiac output, which gets reflected as an augmented pressure wave.
The advantage in using is that it provides insight into systemic effects rather than focusing on either small or large vessels, and as such could be an important indicator of disease progression. However, its usefulness cannot be fully analyzed in the current model since it employs a fixed inflow, simulating the cardiac compensation scenario during which the cardiac output is maintained with hypertension (Acosta et al., 2017). A model including a heart component (e.g. Acosta et al. (2017); Mynard & Smolich (2015)) and 1D vascular beds (e.g Chen et al. (2016); Olufsen et al. (2012); Qureshi et al. (2014)) would be ideal to analyze the complex role of wave reflections in disease progression.
Impedance analysis
Overall, the characteristic features of the impedance spectra resemble those reported by Nichols et al. (2011) (Ch. 16). We note that the impedance moduli are higher in hypertension, whereas the phase show that it starts out negative and then crosses zero at higher frequencies. These results supplement observations from Vanderpool et al. (2011), who studied the impedance spectra using in vitro pulsatile hemodynamics in isolated lungs of the same type of mice.
The results of low frequency components are more subtle. One observation is that the low frequency response is more dynamic in the control animals, reflecting that large vessels are more compliant. A characteristic feature that we were not able to reproduce was a pronounced third harmonic observed in control animals. To our knowledge, this feature has not been reported elsewhere, and could be a consequence of vessel tapering or wave propagation from the vascular beds (which are not included in this study), or it could be an artifact from cycle averaged data. Unfortunately we do not have access to the raw data to confirm or deny this characteristic.
In summary, pulmonary hypertension is characterized by a more resistive and less compliant vasculature with augmented wave reflections, leading to high blood pressure in the pulmonary arteries. This is in line with observations by Lankhaar et al. (2006) and Lungu et al. (2014) for patients with and without pulmonary hypertension. Lankhaar et al. (2006) used a 0D Windkessel model neglecting the effects of wave propagation and arterial stiffness, whereas Lungu et al. (2014) used a coupled 0D1D model representing the pulmonary vasculature by a single vessel.
4.2 Inference of Disease Progression: Linear vs. nonlinear wall mechanics
It is well known that arterial deformation acts nonlinearly imposing increased stiffening with increased pressure (Valdez et al., 2011). Although, most studies confirming this behavior are done in systemic arteries (Eck et al., 2017; Langewouters et al., 1985; Valdez et al., 2011), pulmonary arteries are composed of the same type of tissue and therefore should exhibit similar behavior under dynamic loading (Lee et al., 2016). To understand the effects of nonlinearities on hemodynamic predictions and parameter estimation, we implemented a linear mechanistic and a nonlinear empirical wall model. A qualitatively reasonable outcome, i.e. increased stiffening with pressure, in the form of nonlinear pressure–area curve is evident from Figure 5. However, in the absence of actual data for area deformation, the area predictions using the linear and nonlinear models cannot be validated.
The nonlinear wall model.
Previous empirical nonlinear stress strain models are formulated using sigmoidal functions, relating pressure and area, saturating at both high and low pressures (Langewouters et al., 1985; Valdez et al., 2011). This type of model is characterized by an inflection point determined by a parameter representing the halfsaturation (or maximum compliance) pressure. Valdez et al. (2011) validated these models against dynamic loading data, from ovine thoracic descending aorta (TDA) and carotid artery (CA) under in vivo and ex vivo conditions, showing that pressurearea dynamics lie on the upper (concave down) part of the sigmoidal curves. While these models were able to fit data well, their estimates of the halfsaturation parameter lie outside the known physiological range.
The advantage of the model used here, is that it provides similar estimates, but uses one less parameter, eliminating the need to estimate a parameter that likely is unidentifiable. Moreover, our model provides a basis for theoretical comparison with the linear model at a reference state , since under the small strain assumption, the nonlinear wall model (4) approximates to
(24) 
which is the linear model with (see (10)). This property justifies the interchangeable use of and in (8), used to estimate nominal values for and .
In summary, simulations with the linear wall model are predominantly governed by the 0D boundary conditions, i.e. vascular beds, whereas the nonlinear model modulates both the 1D, i.e. large vessels and 0D domain to predict the hypertensive hemodynamics. This suggests that the linear wall model predicts greater remodeling of the vascular beds due to HPH. However, the general inferences about control and hypertensive hemodynamics, decreased compliance, increased stiffness, resistance and amplitudes of wave reflections, remain the same.
The results shown here indicate that measurements in the MPA can be predicted with either wall model, but predictions in the small vessels differ. Without more data confirming behavior in small vessels it is not possible to conclude which model is better. Moreover, we showed that in hypertensive animals both models provide comparable predictions, likely a result of increased vessel stiffness.
4.3 Model selection
To our knowledge no previous 1D wave propagation studies have implemented statistical criterion for model selection. In this study we have carried out statistical model selection to discriminate between the linear and the nonlinear wall model given available data in the main pulmonary artery. Using AICc and BIC scores. While we have made a parametric assumption about the measurement noise, we have taken its correlation structure into consideration by fitting a set of standard time series models to the residuals. We have focused on the control mice, for which the difference between the linear and the nonlinear model was significant. Our results suggest that the linear model is preferred. One study by Valdez (2010), examining wall properties used the Akaike Information Criterion (AIC) for selecting the wall model. Their results showed that for control animals (they studied sheep) the nonlinear wall model performs better, whereas for stiffer vessels, the linear model performs better. Their former conclusion is contradictory to our model selection analysis of control mouse but the later is consistent with our findings, which show no difference in the predictions using the linear and the nonlinear wall models. However, Valdez (2010)only examined the stressstrain relation in single vessels in absence of fluid dynamics. It should be noted that results presented in this study were done using classical AICc and BIC selection criteria, which have an asymptotic justification. Better approximations that are less reliant on asymptotics are the WatanabeAkaike information criterion (WAIC) (Watanabe, 2010) or the WatanabeBayesian information criterion (WBIC) (Watanabe, 2013). However, these more accurate model selection methods, which we have explored in a recent proofofconcept study (Paun et al., 2018), are based on Monte Carlo simulations and are thus computationally considerably more expensive.
4.4 Limitations
In this study, we did not account for uncertainties in network dimensions and topology associated with image segmentation. We also ignored the uncertainty in the hemodynamic data resulting from ensemble encoding of pressure and flow waveforms over multiple cardiac cycles. Although uncertainty quantification is beyond the scope of this study, it is an important aspect that may have significant impact on the parameter inference. Moreover, it is not clear from image studies if the individual vessels taper, thus the effects of vessel tapering were not considered in this study. Yet, it is known that tapering introduces significant augmentation of pressure along each vessel and throughout the network, which makes it a sensitive model parameter, which could be investigated. Other modeling limitations of this work include the use of a fixed inflow and lumped 0D model for the vascular bed. Finally, the model selection outcome is only valid for the current hemodynamic data, which is available from one location in the main pulmonary artery. The AICc and BIC scores may vary significantly if the waveforms are available from multiple location within the network.
5 Conclusions
We found that the hypertensive mice display significant disease progression associated with remodeling within both large and small vessels. Microvascular remodeling includes reduced compliance and increased resistance, augmenting wave reflections, that combined with increased stiffening of the large vessels lead to significant increase in blood pressure. We also conclude that both linear and nonlinear models can be used to predict the control and hypertensive hemodynamics in the MPA with high accuracy, yet the prediction in the smaller vessels differ. Without more data, it is not possible to select which model better reflects wave propagation along the entire network. These differences were only displayed for control mice, with more compliant vessels. For the hypertensive mice, both large and small vessels are almost rigid and the two models predict the same behavior. Although, the model selection criteria pick the linear wall model for the control muse, these results should not be considered in the statistical context alone as the availability of more physiological data for optimization may alter the present outcome. Finally, analysis of network hemodynamics, wave intensities and impedance moduli indicates an increased presence of wave reflections using the nonlinear model. For this reason, parameter inference and characterization of normal and remodeled vasculature should be regarded as qualitative while using our nonlinear model.
Acknowledgements
Funding This study was supported by the National Science Foundation (NSF) awards NSFDMS # 1615820, NSFDMS # 1246991 and Engineering and Physical Sciences Research Council (EPSRC) of the UK, grant reference number EP/N014642/1.
Conflict of interest The authors declare that they have no conflict of interest.
Appendix
Appendix A Total Vascular Compliance
The volumetric compliance, defined as (ml/mmHg) for a cylindrical vessel with volume , is computed from the linear () and nonlinear () models. For a longitudinally tethered vessel in the network
(25) 
where is the fixed length of the vessel and is computed from (3) and (4), giving
(26) 
where denotes the reference compliance at , given by
(27) 
Computation of the total compliance in a bifurcating network requires compliances to be added in series and in parallel. The total compliance of two vessels and connected in parallel and in series are given by
(28) 
As an example, the total compliance in a single bifurcation network shown in Figure 3c, is given by
(29) 
The total vascular compliance in a vascular region consisting of single bifurcation network and two vascular beds is given by
(30) 
where is computed for the linear or nonlinear wall models using (25)–(29), while is the weighted compliance, as suggested by Alastruey et al. (2016), i.e. where is estimated from the Windkessel models attached to the terminal daughter vessels in the example network.
Appendix B Pulse wave velocity
The pulse wave velocity (PWV), (cm/s), is computed from the eigenvalues of the hyperbolic system of equations (2), from. where
(31) 
Setting and in (31) gives the squared PWV computed for the linear and nonlinear wall models, respectively
(32) 
where is the square of the reference PWV at , given by
(33) 
We use PWV in wave intensity analysis, described in Sec. 2.7, for separating the incident and reflected waves.
Appendix C Nominal parameter values
Control  (mmHg)  (mmHg s/ml)  (s) 
1  37.5  108  0.15 
2  44.5  56  0.06 
3  37.7  47  0.06 
4  40.7  70  0.10 
5  15.6  69  0.13 
6  26.0  87  0.14 
7  31.7  101  0.13 
Hypertensive  (mmHg)  (mmHg s/ml)  (s) 
1  150.6  164  0.09 
2  56.8  107  0.11 
3  123.8  163  0.15 
4  100.7  154  0.13 
5  78.5  143  0.11 
Appendix D Optimized parameter values
For all cases, we optimized and for the wall models, and the global scaling parameters for the Windkessel model, such that
where indicate the nominal quantity. Upper and lower bounds for the optimization intervals are given in Table 5.
Parameter  ()  
Lower bound  0.05  
Upper bound  2.5 
Linear  Nonlinear  
Control  (%)  (%)  
1  67.5  90.7  2.52  29.6  97.1  2.9  199  28.6  2.6  94.1  1.67  31.7  94.0  8.6  147 
2  61.9  28.0  3.15  30.4  97.9  4.5  640  41.0  2.2  34.7  2.33  31.8  98.5  12.4  654 
3  36.1  22.4  2.43  40.3  95.2  8.4  587  15.3  2.5  34.4  1.96  33.2  95.9  6.0  376 
4  56.3  46.9  3.22  29.8  97.7  3.0  369  23.9  3.0  58.9  1.88  33.8  97.7  11.0  348 
5  49.9  52.3  3.42  31.7  97.4  4.1  606  19.8  4.7  65.1  2.24  33.4  96.0  13.5  702 
6  52.1  65.6  2.22  36.0  96.2  5.4  125  22.1  4.0  78.7  1.44  38.8  94.8  12.3  99 
7  70.4  102.9  2.54  32.3  97.6  5.8  737  47.8  3.5  99.8  2.32  32.5  98.2  11.6  698 
HPH  (%)  (%)  
1  151.3  146.2  4.52  56.5  93.6  15.4  49  102.7  1.6  147.5  4.79  58.4  97.5  12.9  32 
2  86.3  94.3  0.13  36.3  95.9  10.2  593  36.6  2.6  100.9  0.12  39.3  96.6  13.0  450 
3  121.0  155.6  9.41  33.3  96.0  7.0  407  82.2  2.1  158.7  8.79  35.7  98.0  6.9  467 
4  117.4  142.2  8.61  42.9  95.5  14.4  63  49.8  2.2  146.4  7.36  44.3  96.5  13.4  41 
5  108.2  130.9  7.81  48.0  94.5  15.0  113  73.5  2.3  133.5  6.82  49.2  95.8  13.0  42 
(mmHg), (mmHg s/ml), (ml/mmHg), (dimensionless), (mmHg), (dimensionless), : least square error.
References
 Acosta et al. (2017) Acosta S, Puelz C, Riviï¿½re B, Penny DJ, Brady KM, Rusin CR (2017) Cardiovascular mechanics in the early stages of pulmonary hypertension: a computational study. Biomech Model Mechanobiol 16: 2093–2112
 Alastruey et al. (2016) Alastruey J, Xiao N, Fok H, Schaeffter T, Figueroa CA (2016) On the impact of modelling assumptions in multiscale, subjectspecific models of aortic haemodynamics. J Royal Society, Interface 13(119) https://doi:10.1098/rsif.2016.0073
 Arnold et al. (2017) Arnold A, Battista C, Bia D, German YZ, Armentano RL, Tran HT, Olufsen MS (2017) Uncertainty quantification in a patientspecific onedimensional arterial network model: EnKFbased inflow estimator. ASME J Verif Valid Uncert 2(1):14
 Antiga et al. (2008) Antiga L, Piccinelli M, Botti L, EneIordache B, Remuzzi A, Steinman DA (2008) An imagebased modeling framework for patientspecific computational hemodynamics. (http://www.vmtk.org). Med Biol Eng Comput, 46:10971112
 Aslanidou et al. (2016) Aslanidou L, Trachet B, Reymond P, Fragasilva RA. Segers P, Stergiopulos N (2016) A 1D Model of the Arterial Circulation in Mice. ALTEX 33:13–28
 Blanco et al. (2014) Blanco PJ, Watanabe SM, Dari EA, Passos MARF, Feijï¿½o RA (2014) Blood flow distribution in an anatomically detailed arterial network. Biomech Model Mechanobiol,13(6):1303–1330
 Boileau et al. (2015) Boileau E, Nithiarasu P, Blanco PJ, Mï¿½ller LO, Fossan FE, Hellevik LR, Donders WP, Huberts W, Willemet M, Alastruey J (2015) A benchmark study of numerical schemes for onedimensional arterial blood flow modelling. Int J Numer Method Biomed Eng, 31(10) https://doi:10.1002/cnm.2732
 Castelain et al. (2001) Castelain V, Hervï¿½ P, Lecarpentier Y, Duroux P, Simonneau G, Chemla D (2001) Pulmonary artery pulse pressure and wave reflection in chronic pulmonary thromboembolism and primary pulmonary hypertension. Journal of the American College of Cardiology, 37(4):1085–1092
 Chen et al. (2016) Chen WW, Gao HG, Luo XY, Hill NA (2016) Study of cardiovascular function using a coupled left ventricle and systemic circulation model. J Biomech, 49(12):2445–2454
 Boggs & Tolle (2000) Boggs P and Tolle J (2000) Sequential quadratic programming for largescale nonlinear optimization Journal of Computational and Applied Mathematics, 124:123–137ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½
 Box & Jenkins (1970) Box GEP, Jenkins GM (1970) Time series analysis: forecasting and control, HoldenDay
 Burnham & Anderson (2001) Burnham KP, Anderson DR (2002) Model Selection and Multimodel Inference: A practical informationtheoretic approach, 2nd ed. SpringerVerlag
 Dujardin, & Stone (1981) Dujardin JP, Stone DN (1981) Characteristic impedance of the proximal aorta determined in the time and frequency domain: a comparison. Med Biol Eng Comput 19:565–568
 Ellwein et al. (2015) Ellwein LM, Marks DS, Migrino RQ, Foley WD, Sherman S, LaDisa JF (2016) ImageBased Quantification of 3D Morphology for Bifurcations in the Left Coronary Artery: Application to Stent Design. Catheter Cardiovasc Interv 87:1244–1255
 Eck et al. (2017) Eck VG, Sturdy J, Hellevik LR (2017) Effects of arterial wall models and measurement uncertainties on cardiovascular model predictions. J Biomech 50:188–194
 Feldkamp et al. (1984) Feldkamp LA, Davis LC, Kress JW (1984) Practical conebeam algorithm. Journal of the Optical Society of America A 1:612–619
 Gan et al. (2007) Gan CT, Lankhaar JW, Westerhof N, Marcus JT, Becker A, Twisk JW, Boonstra A, Postmus PE, VonkNoordegraaf A (2007) Noninvasively assessed pulmonary artery stiffness predicts mortality in pulmonary arterial hypertension. Chest 132(6):1906–1912
 Gerhard & Ray (2010) Gerhard AH, Ray WO (2010) Constitutive modelling of arteries. Proc R Soc A 466:1551–1597
 Humphrey (2008) Humphrey JD (2008) Mechanisms of Arterial Remodeling in Hypertension: Coupled Roles of Wall Shear and Intramural Stress. Hypertension 52(2): 195200
 Hunter et al. (2011) Hunter KS, Lammers SR, Shandas S (2011) Pulmonary Vascular Stiffness: Measurement, Modeling, and Implications in Normal and Hypertensive Pulmonary Circulations. Comp Physiol 1: 1413–1435
 Karau et al. (2001) Karau K, Johnson R, Molthen R, Dhyani A, Haworth S, Hanger C, Roerig D, Dawson C (2011) Microfocal Xray CT imaging and pulmonary arterial distensibility in excised rat lungs. Am J Physiol Heart Circ Physiol 281:H1447?H1457
 Kheyfets et al. (2013) Kheyfets VO, O’Dell W, Smith T, Reilly JJ, Finol EA (2013) Considerations for Numerical Modeling of the Pulmonary Circulation–A Review With a Focus on Pulmonary Hypertension. J Biomed Engin 135:061011–2
 Langewouters et al. (1985) Langewouters GJ, Wesseling KH, Goedhard WJ (1985) The pressure dependent dynamic elasticity of 35 thoracic and 16 abdominal human aortas in vitro described by a five component model. J Biomech 18:613–620
 Lee et al. (2016) Lee P, Carlson BE, Chesler N, Olufsen MS, Qureshi MU, Smith NP, Sochi T, Beard DA (2016) Heterogeneous mechanics of the mouse pulmonary arterial network. Biomech Model Mechanobiol 15:1245–1261
 Lankhaar et al. (2006) Lankhaar JW, Westerhof N, Faes T, Marques K, Marcus J, Postmus P, VonkNoordegraaf A (2006) Quantification of right ventricular afterload in patients with and without pulmonary hypertension. Am J Physiol Heart Circ Physiol 29(4):H1731–173
 Li et al. (2016) Li Y, Parker KH, Khir AW (2016) Using wave intensity analysis to determine local reflection coefficient in flexible tubes. J Biomech 49:2709–2717
 Lungu et al. (2014) Lungu A, Wild JM, Capener D, Kiely DG, Swift AJ, Hose DR (2014) MRI modelbased noninvasive differential diagnosis in pulmonary hypertension. J Biomech 47: 2941–2947
 Lumens et al. (2009) Lumens J, Delhaas T, Kirn B, Arts T (2009) ThreeWall Segment (TriSeg) Model Describing Mechanics and Hemodynamics of Ventricular Interaction. Ann Biomed Eng 37(11): 2234–2255
 (29) McDonald and Attinger The characteristics of arterial pulse wave propagation in the dog. Inform Exch Gp. No 3, Sci Mem 7
 Meaney & Beddy (2012) Meaney JFM, Beddy P (2012) Pulmonary MRA. In: Carr J., Carroll T. (eds) Magnetic Resonance Angiography. Springer, New York, NY
 Mynard et al. (2008) Mynard J, Penny DJ, Smolich JJ (2008) Wave intensity amplification and attenuation in nonlinearflow: Implications for the calculation of local reflection coefficients. J of Biomech 41:3314–3321
 Mynard & Smolich (2015) Mynard JP, Smolich JJ (2015) Onedimensional haemodynamic modeling and wave dynamics in the entire adult circulation. Ann Biomed Eng 43: 144–1460
 Nichols et al. (2011) Nichols WW, O’Rourke MF, Vlachopoulos C (2011) MCDonald’s Blood Flow in Arteries: Theoretical, experimental and clinical principles, 6th edn. London, England: Hodder Arnold.
 Olufsen et al. (2000) Olufsen MS, Peskin CS, Kim WY, Pedersen EM, Nadim A, Larsen J (2000) Numerical simulation and experimental validation of blood flow in arteries with structuredtree outflow conditions. Ann Biomed Eng 28:1281–1299
 Olufsen et al. (2012) Olufsen MS, Hill NA, Vaughan GD, Sainsbury C, Johnson M (2012) Rarefaction and blood pressure in systemic and pulmonary arteries. J Fluid Mech 705:280–305
 Paun et al. (2018) Paun LM, Qureshi MU, Colebank M, Hill NA, Olufsen MS, Haider MA, Husmeier D (2018) MCMC methods for inference in a mathematical model of pulmonary circulation. Statistica Neerlandica (accepted) arXiv:1801.07742.
 Pursell et al. (2016) Pursell ER, Vï¿½lezRendï¿½n D, ValdezJasso D (2016) Biaxial Properties of the Left and Right Pulmonary Arteries in a Monocrotaline Rat Animal Model of Pulmonary Arterial Hypertension. ASME J Biomech Eng 138:11100411100411
 Qureshi et al. (2014) Qureshi MU, Vaughan GD, Sainsbury C, Johnson M, Peskin CS, Olufsen MS, Hill NA (2014) Numerical simulation of blood flow and pressure drop in the pulmonary arterial and venous circulation. Biomech Model Mechanobiol 13(5): 1137–1154
 Qureshi & Hill (2015) Qureshi MU, Hill NA (2015) A computational study of pressure wave reflections in the pulmonary arteries. J Math Biol 71:1525–1549.
 Qureshi et al. (2017) Qureshi MU, Haider MA, Chesler NC, Olufsen MS (2017) Simulating the effects of hypoxia on pulmonary haemodynamics in mice. In: Proc CMBE vol. 1:271–274
 Qureshi et al. (2018) Qureshi MU, Colebank MJ, Schreier DA, Tabima DM, Haider MA, Chesler NC, Olufsen MS (2018) Characteristic Impedance: Frequency or time domain approach? Physiol Meas 39(1):014004 https://doi.org/10.1088/13616579/aa9d60
 Rasmussen & Williams (2006) Rasmussen CE, Williams CKI (2006) A computational study of pressure wave reflections in the pulmonary arteries. J Math Biol 71:1525–1549.
 Reymond et al. (2009) Reymond P, Merenda F, Perren F, Rufenacht D, Stergiopulos N (2009) Validation of a onedimensional model of the systemic arterial tree. Am J Physiol Heart Circ Physiol 297: H208–H222
 Riches et al. (1973) Riches AC, Sharp JG, Thomas DB, Smith SV (1973) Blood volume determination in mouse. J Physiol 228(2):279–284
 Rich et al. (2011) Rich JD, Shah SJ, Swamy RS, Kamp A, Rich S (2011) Inaccuracy of Doppler echocardiographic estimates of pulmonary artery pressures in patients with pulmonary hypertension: implications for clinical practice. Chest 139:988–993
 Safaei et al. (2016) Safaei S, Bradley CP, Suresh V, Mithraratne K, Muller A, Ho H, Ladd D, Hellevik L, Omholt SW, Chase JG, Mï¿½ller LO, Watanabe SM, Blanco PJ, de Bono B, Hunter PJ (2016) Roadmap for cardiovascular circulation model. J Physiol 594(23):6909–6928
 Schwarz (1978) Schwarz G (1978) Estimating the dimension of a model. Ann Statist 6:461–464
 Stergiopulos et al. (1995) Stergiopulos N, Meister JJ, Westerhof N (1995) Evaluation of methods for estimation of total arterial compliance. Am J Physiol 268:H1540–1548
 Simonneau et al. (2013) Simonneau G, Gatzoulis MA, Adatia I, Celermajer D, Denton C, Ghofrani A, Sanchez MA, Kumar RK, Landzberg M, Machado RF, Olschewski H, Robbins IM, Souza R (2013) Updated clinical classification of pulmonary hypertension. J Am Coll Cardiol 62:D34D41
 Tabima et al. (2012) Tabima DM, RoldanAlzate A, Wang Z, Hacker TA, Molthen RC. Chesler NC (2012) Persistent vascular collagen accumulation alters hemodynamic recovery from chronic hypoxia. J Biomech 45:799–804
 Tang et al. (2012) 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 imagebased, computational fluid dynamics study. Pulm Circ 2(4): 470–476
 Tran et al. (2017) Tran JS, Schiavazzi DE, Ramachandra AB, Kahnb AM, Marsden AL (2017) Automated tuning for parameter identification and uncertainty quantification in multiscale coronary simulations. Computers and Fluids 142:128–138
 Tuder et al. (2007) Tuder RM, Marecki JC, Richter A, Fijalkowska I, Flores S (2007) Pathology of Pulmonary Hypertension. Clinics in Chest Medicine 28(1):23–vii
 Tawhai et al. (2011) Tawhai MH, Clark AR, Burrowes KS (2011) Computational models of the pulmonary circulation: Insights and the move towards clinically directed studies. Pulm Circ 1(2):224–238
 Vanderpool et al. (2011) Vanderpool RR, Kim AR, Chesler NC (2011) Effects of acute Rho kinase inhibition on chronic hypoxiainduced changes in proximal and distal pulmonary arterial structure and function. J Appl Physiol 110: 188–198
 Valdez et al. (2011) ValdezJasso D, Bia D, Zï¿½calo Y, Armentano RL, Haider MA, Olufsen MS (2011) Linear and nonlinear viscoelastic modeling of aorta and carotid pressurearea dynamics under in vivo and ex vivo conditions. Ann Biomed Eng 39:1438–1456
 Valdez (2010) ValdezJasso D (2010) Modeling and Identification of Vascular Biomechanical Properties in Large Arteries. PhD Thesis, North Carolina State University, Raleigh, NC
 van de Vosse & Stergiopulos (2011) van de Vosse FN, Stergiopulos N (2011) Pulse wave propagation in the arterial tree. Annu Rev Fluid Mech 43:467–499
 Yushkevich et al. (2006) Yushkevich PA, Piven J, Hazlett HC, Smith RG, Ho S, Gee JC, Gerig G (2006) Userguided 3D active contour segmentation of anatomical structures: Significantly improved efficiency and reliability. (www.itksnap.org) Neuroimage 31:1116–1128
 Watanabe (2010) Watanabe S (2010) Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. J Machi Learn Res 11:3571–3594
 Watanabe (2013) Watanabe S (2013) A Widely Applicable Bayesian Information Criterion. JJ Machi Learn Res 14:867–897
 Wang & Chesler (2011) Wang Z, Chesler NC (2011) Pulmonary vascular wall stiffness: An important contributor to the increased right ventricular afterload with pulmonary hypertension. Pulm Circ 1(2): 212–223, 2011
 Westerhof et al. (2009) Westerhof N, Lankhaar J, Westerhof B (2009) The arterial windkessel. Med Biol Eng Comput 47:131–141
 Williams et al. (2014) Williams ND, WindWillassen O, Wright AA, Program REU, Mehlsen J, Ottesen JT, Olufsen MS (2014) Patient specific modeling of headup tilt. Math Med Biol 31:365–392
 Willemet & Alastruey (2015) Willemet M, Alastruey J (2015) Arterial Pressure and Flow Wave Analysis Using TimeDomain 1D Hemodynamics. Ann Biomed Eng 43:190–206
 Windberger et al. (2003) Windberger U, Bartholovitsch A, Plasenzotti R, Korak KJ, Heinze G (2003) Whole blood viscosity, plasma viscosity and erythrocyte aggregation in nine mammalian species: reference values and comparison of data. Exp Physiol 88:431–440

Yang et al. (2016)
Yang W, Feinstein J, VignonClementel I (2016)
Adaptive outflow boundary conditions improve postoperative predictions after repair of peripheral pulmonary artery stenosis.
Biomech Model Mechanobiol 15(5):1345–1353
M. Umar Qureshi (Corresponding author), Department of Mathematics, North Carolina State University, Raleigh, North Carolina 27695, USAEmail address: muquresh@ncsu.edu Mitchel Colebank, Mansoor A. Haider, Mette S. Olufsen Department of Mathematics, North Carolina State University, Raleigh, North Carolina 27695, USALaura Ellwein Department of Mathematics and Applied Mathematics, Virginia Commonwealth University, Richmond, VA 23284, USA Naomi Chesler Department of Biomedical Engineering and Madison, University of Wisconsin at Madison, WI, 537061609, USA Mihaela Paun, Dirk Husmeier, Nicholas A. Hilll School of Mathematics and Statistics, University of Glasgow, Glasgow G12 8QQ, UK