Lumped‐parameter models of the pulmonary vasculature during the progression of pulmonary arterial hypertension

Abstract A longitudinal study of monocrotaline‐induced pulmonary arterial hypertension (PAH) was carried out in Sprague‐Dawley rats to investigate the changes in impedance (comprising resistance and compliance) produced by elevated blood pressure. Using invasively measured blood flow as an input, blood pressure was predicted using 3‐ and 4‐element Windkessel (3WK, 4WK) type lumped‐parameter models. Resistance, compliance, and inductance model parameters were obtained for the five different treatment groups via least‐squares errors. The treated animals reached levels of hypertension, where blood pressure increased two folds from control to chronic stage of PAH (mean pressure went from 24 ± 5 to 44 ± 6 mmHg, P < 0.0001) but blood flow remained overall unaffected. Like blood pressure, the wave‐reflection coefficient significantly increased at the advanced stage of PAH (0.26 ± 0.09 to 0.52 ± 0.09, P < 0.0002). Our modeling efforts revealed that resistances and compliance changed during the disease progression, where changes in compliance occur before the changes in resistance. However, resistance and compliance are not directly inversely related. As PAH develops, resistances increase nonlinearly (R d exponentially and R at a slower rate) while compliance linearly decreases. And while 3WK and 4WK models capture the pressure‐flow relation in the pulmonary vasculature during PAH, results from Akaike Information Criterion and sensitivity analysis allow us to conclude that the 3WK is the most robust and accurate model for this system. Ninety‐five percent confidence intervals of the predicted model parameters are included for the population studied. This work establishes insight into the complex remodeling process occurring in PAH.


Introduction
Pulmonary arterial hypertension (PAH) is characterized by elevated blood pressure in the pulmonary circulation. During the development of PAH, the pulmonary vasculature undergoes structural remodeling that compromises its normal physiological function. In particular, the pulmonary blood vessels undergo currently untreatable changes and develop lesions that alter their biomechanics and mechanobiology (Stenmark et al. 2009). PAH also subjects the right ventricle of the heart to pressure overload. While cardiac remodeling initially compensates for the increase in pulmonary vascular stiffening, eventually this response becomes maladaptive and, if left untreated, PAH eventually results in cardiac decompensation and heart failure (Vonk-Noordegraaf et al. 2013). The disease is very aggressive with a 3-year survival rate of only 48% (Barst 2008) and a prevalence among women is twice as high as in men (Rich et al. 1987). To date, transplantation of the lung remains the only curative treatment option.
Although it has been established that in PAH the pulmonary vasculature stiffens, the causes initiating this remodeling process remain unknown. Clinical studies have shown that pulmonary vascular impedance, a combination of resistance and compliance, is an independent risk factor (Weinberg et al. 2004;Dyer et al. 2006;Stenmark et al. 2009). Thus, it is critical to study the alteration in pulmonary vascular resistance and capacitance to understand the pathophysiology of PAH. Monocrotalineinduced PAH is a commonly used animal model for the study of this disease. After 4 weeks of being injected with monocrotaline (MCT), the animals reach chronic PAH pressure levels and display vascular lesions (Gomez-Arroyo et al. 2012;Maarman et al. 2013) such as those found in human PAH patients. Since no longitudinal study of the disease has been carried out, our goal was to characterize the time-course of changes in vascular resistance and capacitance during the development of the disease in an MCT animal model of PAH. Blood pressure and flow were measured invasively at four stages of the disease. The evolving properties of the pulmonary vasculature were studied by applying the commonly used three-and four-element Windkessel models. Blood flow measurements were used to initiate the Windkessel models, and the resulting blood pressure predictions were compared to the in-vivo measured data. Model parameters were then used to determine alterations of the resistive and compliant properties of the pulmonary system in the time-course of PAH. The results showed distinct time-courses of altered resistance and compliance and will be useful for formulating and testing new computational models.

Methods & Materials Data acquisition
All experimental procedures performed were in accordance with the University of Illinois at Chicago Animal Care and Use Committee. Male Sprague-Dawley rats (Charles River Laboratories, Chicago, IL) with a mean body weight of 272 AE 30.1 g were included in the study. At 8 weeks of age, the animals were divided into two groups and injected subcutaneously on the posterior aspect of the neck. A group of 27 rats was injected with a dose of 60 mg/kg of monocrotaline (Sigma Aldrich, St Louis, MO, USA) dissolved in HCl and diluted in di-H 2 O and NaOH at a concentration of 0.025 g/mL. A control group (PL) that accounted for aging effects was injected with saline solution at a volume calculated from the monocrotaline dosage. When possible, a treated (PAH) animal and control (PL) animal were housed together for up to 4 weeks, and underwent a terminal surgical procedure at the end of their specified stage.
For each week postinjections, paired PL (N = 17) and PAH rats underwent an open-chest surgery to obtain invivo hemodynamic measurements. First, the animals were anesthetized with oxygen at 4% isoflurane, intubated, and placed on a ventilator (SAR-1000 Ventilator, CWE Inc., Ardmore, PA, USA). The breathing rate was set to 50 breaths/minute and a tidal volume that was based on the animal weight at time of surgery (Pacher et al. 2008). Following a thoracotomy procedure, the heart and large blood vessels were exposed. A 1.6F dual pressure sensor catheter (Transonic Scisense, Ontario, Canada) was inserted into the right ventricle and advanced to the main pulmonary artery (MPA) via the pulmonic valve. An ultrasonic flow probe (Transonic Scisense, Ontario, Canada) was wrapped around the MPA, while continuously measuring blood pressure in the right ventricle and MPA. An illustration of the set-up for in vivo hemodynamic data measurements is shown in Figure 1A. Time series of pulmonary arterial and right ventricular pressures and blood flow (Fig. 1B) were collected in LabChart software (ADInstruments Inc. Colorado Springs, CO) and further analyzed in a custom-written MATLAB (version R2016a, The Mathworks, Natick, MA) code for the calculation of hemodynamic parameters such as heart rate, end-systolic, and end-diastolic pressures.

Windkessel models
To identify the individual changes and contributions of resistance and compliance in the vasculature, Windkessel models were used to relate blood flow Q to blood pressure P. These circuit element-based models shown in Figure 2 are used as analogies to known physiological properties of the vasculature. In particular, the buffering of pulsatile blood flow performed by the large, proximal vessels can be described through a combination of capacitive, or total vascular compliance (C), and proximal resistive (R) elements. The function of the distal vessels to convert pulsatile blood flow to steady flow, and thus preventing capillary rupture, can be represented through a highly resistive element (R d ).
To describe the inertial aspect of the pulmonary system, an inductor (L) circuit element is included to account for the inertial effects of blood. The number of electrical elements used leads to the three-or fourelement Windkessel models, with L being absent in the three-element model.
gis the set of model parameters. In the absence of inductance, the four-element model reduces to the three-element Windkessel model (Equation 2), viz.
with the set of model parameters h 3WK ¼ R; R d ; C f g .

Numerical implementation
One cardiac cycle of blood flow was used to predict blood pressure for each Windkessel model. The model parameters h model were estimated via solution of the inverse  problem of the predicted and measured values of blood pressure. Starting from a set of initial values h 0 , the inverse algorithm iteratively estimated the set of parametersĥ that minimize the normalized root mean square error between the experimental p j and model predicted P (t j ,h) measured at time points t j . The objective function was defined using a two-norm error function as shown in Equation 3.
Here m is the number of data points. To numerically implement this objective function, a constrained Nelder-Mead simplex method (fminsearchbnd function in MATLAB) was used (D'Errico 2012). This bounded function constrains the parameters to be positive to ensure physical properties (D'Errico 2012). At each iteration of the minimization method, the algorithm solved for pressure P(t j ,h) using a fourth-order Runge-Kutta method (ode45 function in MATLAB). Once a tolerance of 10 À6 normalized error was met, fminsearchbnd function exited, output an optimal set of model parametersĥ, and convergence of the algorithm was confirmed.

Initial conditions
The model parameters related to characteristic and total vascular resistance were initialized via input impedance analysis. By mapping three consecutive blood pressure and flow waveforms to the frequency domain using Fourier series, input impedance spectrum (Z in ) for each individual animal was computed (O'Rourke 1982). Briefly, the magnitude of these Fourier coefficients of the pressure and flow signals were then expressed as function of the harmonic term P(n) and Q(n), respectively. The ratio between P(n) and Q(n) amplitudes was taken to determine Z in as shown in Equation 4. As used in Segers et al., the reflection coefficient (Γ) was computed by taking the following ratio (Z in -Z c )/(Z in +Z c ) and the wave-reflection intensity was determined by evaluating the reflection coefficient (Γ 1 ) at the first harmonic or Z in (n = 1) .
Total vascular resistance (Z 0 ) and characteristic impedance (Z c ) are metrics based on the impedance magnitude at the 0 th order harmonic, and from the average of the higher frequencies from 3rd to 10th harmonics, respectively . Using these values from each individual animal, the model parameters R d and R were initialized using the relations R 0 (Coogan et al. 2011). Inductance was initially set to 0.002 AE 0.001 mmHgÁmin 2 ÁmL À1 based on work in mice aorta by Segers et al. (2005). Compliance C, not being a parameter obtainable from impedance analysis, was determined using the pulse pressure method (PPM) as described by Stergiopulos et al. (1994Stergiopulos et al. ( , 1999. Here, initial compliance is estimated by adjusting the capacitance value of an electrical body (resistor and capacitor in parallel). Blood flow and pressure are used as input and the resistor value is fixed at Z 0 . Once the pulse pressure from the model and data match, the value of the capacitor is assigned to the initial compliance C 0 . The initial set of parameters h 0 is shown in Table 1.

Sensitivity analysis
To identify the contribution of each model parameterĥ i in the prediction of pressure P(t j ,h) for each model, sensitivity analysis was carried out. The sensitivity functions Table 1. Summary statistics of input impedance Z 0 , characteristics impedance Z c , wave-reflection intensity |Γ 1 |, and arterial compliance based on pressure and flow measurements for each group. .
where n p is the number of model parameters, and m the number of data points in each time series.

Model selection
The two models were compared using the Akaike Information Criterion (AIC). This criterion allows for model comparison when models contain different number of parameters. The number of parameters used in a model is weighted into the error criterion E as follows:(Valdez-Jasso 2010)

Statistical analysis
Summary statistics were computed for the initial condi- , hemodynamics and the set of optimized parameters from the Windkessel models (R d , R, C, L). An analysis of variance (ANOVA) was performed using a Standard Least Squares statistical model to identify differences in means of the parameters according to the disease stage (PL, PAH1, PAH2, PAH3, PAH4). For the univariate recordings of R d , R, and C parameters effects of disease stage (PL, PAH1, PAH2, PAH3, PAH4) and estimation method (3WK, 4WK) were investigated. All the statistical analyses were carried out in JMP Pro Statistical Software (Version 13.0.0, © SAS Institute Inc., Cary, NC). Multiple comparisons were performed using Tukey HSD (Honest Significant Difference) test to determine which differences in average parameter recording among disease stages and estimation methods were significant. For all statistical tests, significance level a was set at 0.05. All values are presented as mean AE standard deviation. Based on the optimized model parameters a 95% confidence interval is included representing each treatment group.
Goodness-of-fit was measured via coefficient of determination as shown in Equation 8.
the mean of the observed data p j , P(t j , h) the predicted data, SS tot the total sum of squares and SS res the sum of squares of residuals.

Results
Blood pressure rose in the animals while blood flow remained overall constant. The changes in resistive and compliant properties of the vessels are manifested in the parameter estimation of the 3WK and 4WK model parameters. Blood pressure predictions of the Windkessel models resembled closely the measured pressure. The metrics of model robustness, sensitivity, and accuracy quantified these differences.

Hemodynamics
Representative MPA blood pressure (Fig. 1B-top) and right ventricular pressure (Fig. 1B-bottom) measurements indicate how the pressure ranges steadily increases as time (weeks) progresses. However, blood flow remains overall constant ( Fig. 1B-middle). The normotensive group had a systolic pressure of 33.6 AE 4.9 mmHg with a mean pulmonary arterial pressure (mPAP) of 24.1 AE 4.9 mmHg, fitting within a normotensive group. In early disease state (PAH1, N = 7) the waveform is almost identical to PL (N = 17), with a systolic pressure of 33.8 AE 1.5 mmHg and mPAP of 24.7 AE 2.1 mmHg. Between mild and advanced PAH (PAH2, N = 7 and PAH3, N = 7) the waveforms displayed a sharper increase during systole and a more rapid decrease during diastole coupled with overall increase in mPAP. At the advanced stages of PAH (PAH4, N = 5), the systolic pressure increased to 72.5 AE 8 mmHg, with a mPAP of 44.6 AE 5.7 mmHg, an exaggerated phenotype seen in PAH2 (N = 8) and PAH3 (N = 7). The statistically significant increases in mPAP throughout the different weeks confirmed that MCT induced different stages of PAH. However, the mean flow (or cardiac output) did not show a significant difference among the groups, indicating that the right ventricular performance had not been compromised. The summary statistics of these hemodynamic parameters are presented in Figure 3.

Model predictions
All the Windkessel models investigated closely predicted the in-vivo measured blood pressures. Qualitatively, the three-and four-element models predicted blood pressure in the same manner and were able to mimic the dicrotic notch at all the PAH stages, with more accuracy at the early stages of PAH. The measured pressure (solid black lines) and the predicted pressures (dashed red lines) can be visualized in Figure 4. Quantitatively, values of the coefficient of determination R 2 were close to 1.0, and the values in the normalized root mean square error were nearly 0 (Table 2). R 2 values had a decreasing trend with 0.97 in placebos and 0.93 by advanced PAH. It is important to note the R 2 values were identical for the two models. Summary statistics of the optimized set of parameters are shown in Table 2.

Parameter estimation
The prediction of the parameters R d , R, and C using the three-and four-element models was sensitive to the stage of PAH. Both distal and proximal resistances R d and R nonlinearly increased as the disease progressed in the 3WK and 4WK models ( Fig. 5A and C, respectively). Specifically, R d exponential increased while R nonlinearly rose at a slower rate. On the other hand, the total arterial compliance C decreased in a linear fashion (Fig. 5B). Thus, the rate of change in resistances and compliance are not symmetric or proportional to each other. More interestingly, the changes in compliance occur by week 2, earlier than any major changes in resistance (R and R d ). This calls to further investigate the mechanisms underlying the remodeling process of the pulmonary vasculature and identify the relation between the resistive and compliant properties of this nonlinear system. R d was statistically larger in the later stages of PAH compared to the normotensive state. In the models, C significantly decreased from normotensive to the chronic stage of PAH. The parameter R displayed a pattern like R d , but leveled off at weeks 3 and 4 with a significant increase from PL (N = 17) and PAH2 (N = 8) to PAH4 (N = 5). Lastly, the inductor L from the 4WK model varied across the PL (N = 17) and treated MCT groups with PAH4 (N = 5) ranging close to zero. The parameter L was transformed using the natural logarithm, ln(L), in an effort to detect statistical differences, however, there was The Physiological Society and the American Physiological Society not a significant change in the parameter with PAH progression (P < 0.3834). Method predictions were not statistically significant different between 3WK and 4WK for none of the parameters (R d : P < 0.9698, R: P < 1.0000, C: P < 0.8775). P-values for treatment and method effects are summarized in Table 3, where the values less than 0.05 were considered statistically significant.

Sensitivity analysis
The ranking of the optimized R d , R, and C parameters were on the same order of magnitude. R d had the largest magnitude value, almost twofold higher than the other parameters indicating that R d is the biggest contributor in the model. This was followed by the compliance parameter C and R. Specifically, in the 3WK model R d rank ranged between 5.4 AE 0.8 and 8.3 AE 0.8, whereas C from 2.6 AE 0.2 to 3.7 AE 0.3, and R from 1.6 AE 6.1 to 2.0 AE 6.7, respectively ( Table 4). The R d and C parameters had the highest numerical ranking and similar values across both Windkessel models and their relative contribution to the models remained unaffected by the disease. For the four-element Windkessel models, L was two orders of magnitude lower than the rest of the parameters. Due to this large difference in magnitude and thus low ranking, this indicated the contribution of L to be significantly smaller than the other parameters. In this case, the contribution of the inductor L to the model prediction can be neglected.

Model selection
Normalized root mean square error E and the AIC values were used to determine the model closest to the data while accounting for the number of parameters used. As shown in Table 2, E and AIC values for 4WK and 3WK are within the same ranges for each stage of PAH. The cost of adding a fourth parameter was, therefore, incremental. Furthermore, confidence interval results show the parameters, excluding L, initially have a relative narrow range (Table 5). As the disease progresses, these intervals become larger, which may be attributed to the less accurate fitting of the severe PAH conditions. L was shown not to statistically change over the progression of PAH and the confidence intervals further reflect the low contribution of this parameter to the model predictions.

Discussion
In this study, pulmonary vascular remodeling properties were investigated in an animal model of PAH. Using lumped-parameter Windkessel models, in-vivo measurements from the main pulmonary artery were used to identify changes in vascular resistance, characteristic impedance, vascular compliance, and inertial effects at each stage of the disease. As blood pressure increased in the pulmonary system, the resistive (R and R d ) and compliant (C) properties of the vessels changed, compliance changing earlier in the time-course. Inductance contributions were not significant in the 4WK model. Outcomes from AIC and sensitivity analysis indicate that the 3WK is sufficient to represent the impedance characteristics of the pulmonary system during the progression of PAH in the MCT-animal model.

Parameter estimation
The computed initial conditions followed the expected behavior and provided with a good starting point for the optimization process. Several studies have shown a significant increase in total vascular resistance (Z 0 ) with increased afterload in both systemic and pulmonary systems (Murgo et al. 1980;Segers et al. 2005;Hunter et al. 2008). Our findings showed this trend from normotensive to chronic disease state, as well as the significant increases during the disease progression. The average of the higher harmonics, represented by the characteristic impedance Z c , also tended to increase, albeit less significantly. As expected, our initial C 0 was lower in value compared to the optimized parameter C, as has been seen previously Table 3. P-values for treatment and method effects in the Windkessel parameters (R d , R, C, and L).

Parameter/effect Treatment Method
Statistical significance is considered when P < 0.05. Treatment had significant additive effect on the R d , R, and C. Method did not show any difference between the 3-and 4-element models. by Segers et al. (1999). This was attributed to the pulsepressure method focusing on lower harmonics, which is where compliance is mostly affected, therefore possibly giving a more accurate representation ).
Overall, decreases in C from week to week in each model indicate a stiffening of the pulmonary vasculature (Thenappan et al. 2016), which has also been observed in the systemic circulation (Shirwany and Zou 2010;Wu et al. 2015). The impedance values Z 0 -Z c and Z c differed and were underestimated from those based on the parameter estimation routine. This can be explained as the model parameter C is not identified via input impedance analysis, overestimating R d and R. Finally, the wave-reflection intensity increased progressively with the disease. This parameter is important in confirming the change in overall stiffness in the pulmonary vasculature as the disease progressed. As PAH becomes more severe, the wave-reflection intensity, and therefore stiffness in the vasculature increased (Nichols et al. 2008). As PAH progressed, the optimized parameters followed the same trend as the initial conditions. The differences, however, were that the optimized parameters were able to distinguish the unique contribution of resistance and compliance of the vasculature. More importantly, the optimized model parameters offered a mechanistic understanding of the system allowing for a more physiologically accurate representation. It was observed that the more severe PAH is, the stiffer the vasculature becomes, primarily because of decreased C (Thenappan et al. 2016). More specifically, the rate of increase in the resistance is not proportional to the rate of decrease in compliance. Also, changes in compliance occur before major changes in resistance. This indicates an inverse nonlinear relationship between the two parameters and further modeling efforts are needed to address the underlying mechanism of pulmonary vascular remodeling during PAH. To the best of our knowledge, the outcomes presented in this study are the first to indicate how the resistance and compliant properties are changing as a function of the PAH progression in the pulmonary vasculature. Lankhaar et al. (2006) conducted a study with human patients with different severities of PAH (chronic thromboembolic and idiopathic), however, this does not display how the severity changes as a function of disease in a specific individual.
In a week-to-week comparison, R d follows an exponential increase as the disease progresses, whereas C follows a linear decrease. In fact, sooner changes in compliance can be attributed to an early remodeling process occurring in the larger vessels instead of the resistive vessels. Due to the models being a nondimensional model, the exact interpretation of these results is limited (Westerhof et al. 2009). Future work will involve the development of a one-dimensional fluid model in order to incorporate wall shear stress and add spatial components to determine the locations of where resistance and compliance change in the vasculature. This would reveal whether the behavior of the proximal vessels during disease progression imposes the behavior on the distal vessels, or vice versa.

Model selection
Both models described the physiological system by accurately predicting the dicrotic notch in the pressure waveform, as well as providing physiologically realistic values. However, given that the 3WK model offers the same evaluation than the 4WK does, including inductance for the modeling of PAH impedance is not necessary. In work done by Lambermont et al. (1997) the inclusion of L was important to improve their data fitting but they did not evaluate the contribution of the inductor through a sensitivity analysis. Our findings show, however, that the Table 5. Confidence intervals for the estimated model parameters for the three-and four-element Windkessel models for a male Sprague-Dawley rat population of MCT-induced PAH.

Model
Group R d (mmHgÁminÁmL À1 ) R (mmHg min mL À1 ) C [mL/mmHg] ln(L) (mmHgÁmin 2 ÁmL À1 ) contribution of the inductor in the four-element model is negligible. In addition, the large confidence interval of L further confirms the uncertainty of this parameter to model the pulmonary vascular hemodynamics. We suggest that the three-element Windkessel model is sufficient to accurately and robustly represent the pulmonary system dynamics. While both models predicted the blood pressure, the model-to-data discrepancy increased as PAH progresses. We speculate that the pressure waveforms may become more complex as the disease progresses and the Windkessel models may be limited to capture these changes in hemodynamics. This in turn motivates the development of a 1D fluid model to more accurately capture pulmonary vascular remodeling in PAH.
To summarize, we found that the three-element Windkessel model organized in series provided accurate and physiological predictions for resistive and compliant properties in the pulmonary vasculature. It is evident that with the progression of PAH, total vascular resistance is increasing exponentially, while the overall compliance has a steady, linear decrease. Furthermore, more pronounced changes are occurring earlier in compliance. The interaction of the resistant and compliant properties may give insight into where and how the pulmonary vasculature is remodeling as a function of disease state and be useful in future modeling (i.e., 1D fluid model) of PAH. More specifically, based on where the linear compliant and nonlinear resistant changes are occurring, targeted treatment methods could be devised to mitigate the issues in those regions.

Limitations
The main limitations of this model were the lack of spatial representation of the pulmonary system as well as lack of wave transmission information. In the future, with proper data acquisition, this model can be incorporated as a boundary condition for a one-dimensional fluid model to give understanding of the complex properties of the pulmonary vasculature as it undergoes PAH progression. In addition, the study was only limited to male rats. Even though the incidence of PAH in women is almost twice as in men (Foderaro and Ventetuolo 2016), a study by Ventetuolo et al. (2014) reported that female patients have a better survival than male patients. Furthermore, while PAH incidence is higher in females than in males, rodent models of PAH have suggested that estrogen play an important role in the pulmonary vasculature, producing beneficial effects in response to PAH (Pugh and Hemnes 2010). Future work will study how compliance and resistance change over the disease progression in female rats.