Parameter inference in a computational model of hemodynamics in pulmonary hypertension
Mitchel J. Colebank, Amanda L. Colunga, Macie King, Christopher Schell, Matt Sheldon, Mariam Kharbat, Robert Sternquist, Mette S. Olufsen
PParameter inference in a computational model of hemodynamics in pulmonary hypertension
Mitchel J. Colebank ∗ , Amanda L. Colunga ∗ , Macie King, Christopher Schell, Matt Shelden, Mariam Kharbat,
Robert Sternquist, Mette S. Olufsen
Department of Mathematics, North Carolina State University, Raleigh, NC
Running title: Parameter inference in pulmonary hypertension * Mitchel Colebank and Amanda Colunga are joint first authors.
Corresponding author
Mette Olufsen SAS Hall 3216 Department of Mathematics, Box 8205 North Carolina State University Raleigh, NC 27695 Email: [email protected] Phone: 919-515-2678 bstract
Pulmonary hypertension (PH), defined by a mean pulmonary arterial pressure (mPAP) mmHg higher than 20 mmHg, is caused by vascular remodeling, increasing pulmonary vascular resistance and decreasing pulmonary compliance. The disease is progressive with a low 5-year survival rate. There are few measurable biomarkers of PH progression, and conclusive diagnosis of the disease requires invasive right heart catheterization (RHC). This study develops a systems-level model calibrated to RHC data, providing a noninvasive tool for identifying disease progression indicators. The model is formulated using an electrical circuit analog, including all four heart chambers and the pulmonary and systemic circulations. We calibrate the model considering two quantities of interest: systolic, diastolic, and mean RHC pressure measurements and a combination of static RHC data and time-varying pressure waveforms. Local and global sensitivity analyses are applied to determine influential, identifiable parameter subsets estimated to fit the model to data from five PH patients, three with chronic thromboembolic PH (CTEPH) and two with pulmonary arterial hypertension (PAH). From the calibrated model, we compute relevant outcomes, including cardiac power, resistance and compliance ratios, and a pulsatility index, and compare predictions between PH patients and normotensive controls. Results show that both CTEPH and PAH patients have elevated pulmonary vascular resistance and decreased compliance relative to normotensive pa-tients. Both the pulsatility index and right ventricular cardiac power increase in PH. Lastly, we simulate treatment strategies for all five PH patients. Consistent with clinical knowledge, treatment predictions show that CTEPH is curable, whereas PAH is not.
Key Points • PH is a progressive disease that leads to elevated PVR and right ventricular dysfunction. • An eight-compartment systems-level cardiovascular model is calibrated using a combination of static and dynamic data from PH patients, estimating parameters indicative of disease pro-gression. • Sensitivity analysis and model selection are used to identify the most informative parameter subset for parameter inference. • Model outcomes are compared between normotension and PH, revealing an increase in pul-monary vascular resistance and a decrease in pulmonary vascular compliance.
Simulations show that CTEPH treatment reduces the main pulmonary artery pressure below the hypertensive limit (20 mmHg), while PAH treatment is less effective.
Keywords
Pulmonary hypertension; computational model; parameter inference; cardiovascular modeling
Abbreviations
0D Zero dimensional (a model that only depends on time) 1D One dimensional (a model that depends on time and one spatial component) AIC Akaike information criteria CO Cardiac output CTEPH Chronic thromboembolic pulmonary hypertension ECG Electrocardiogram ET-1 Endothelin-1 iid Independent and identically distributed IRB Institutional review board MPA Main pulmonary artery mPAP Mean pulmonary arterial pressure mSAP Mean systemic arterial pressure NO Nitric oxide ODE Ordinary differential equation PDE5 Phosphodiesterase type 5 PAH Pulmonary arterial hypertension PAWP Pulmonary arterial wedge pressure PH Pulmonary hypertension PV Pressure-volume PVR Pulmonary vascular resistance RA Right atrium RHC Right heart catheterization RV Right ventricle VVI Ventricular-ventricular interaction
NTRODUCTION
Patients with a resting mean pulmonary arterial blood pressure (mPAP) higher than 20 mmHg are diagnosed with pulmonary hypertension (PH) (Simonneau et al. , 2019). This disease has no cure and, if left untreated, it progresses rapidly, causing thickening and stiffening the pulmonary arteries and veins, vascular-ventricular decoupling, and right ventricular (RV) dysfunction (Fayyaz et al. , 2018; Hoeper et al. , 2017). There are five main etiologies of PH, including (1) pulmonary arterial hypertension (PAH), (2) PH due to left heart disease, (3) PH due to lung disease and/or hypoxia, (4) chronic thromboembolic PH (CTEPH), and (5) PH with unclear multifactorial mechanisms (Foshat & Boroumand, 2017). Of these, only patients in groups (1, PAH) and (4, CTEPH) have PH as their primary disease. In contrast, PH is a comorbidity in groups (2) and (3) patients, and patients in group (5) have multiple symptoms and therefore do not constitute a coherent group. This study analyzes data from PAH and CTEPH patients, as these patient-groups have PH as their primary disease. PAH and CTEPH patients typically experience symptoms found in numerous diseases, includ-ing shortness of breath, dizziness and fainting, fatigue, and swelling of the legs and abdomen (Dunlap & Weyer, 2016), making the disease difficult to detect. These patients undergo several tests, and a definite diagnosis requires right heart cardiac catheterization (RHC), providing inva-sive measurements of pulmonary arterial blood pressure (Montani et al. , 2013). PH symptoms do not appear until 1-2 years after disease onset (Hoeper et al. , 2017); at this stage, patients typically have significant vascular remodeling. For this reason, an understanding of how disease markers (e.g., pulmonary vascular resistance (PVR) and compliance) are modulated with the disease can assist in early detection and developing better therapeutic interventions. To address this, we con-struct a systems-level model calibrated using static (diastolic and systolic pressure values) and dynamic pressure waveforms. We employ local and global sensitivity analyses to determine which parameters change with disease and study the effects of standard treatments on model predictions. Mathematical modeling is becoming increasingly useful as a tool for monitoring and under-standing cardiovascular disease progression (Ellwein et al. , 2008; Kung et al. , 2014; Marquis et al. , 2018; Colebank et al. , 2019; Colunga et al. , 2020). For example, Colunga et al. (2020) developed a zero-dimensional (0D) systems-level model, predicting pressure-volume (PV) loops and the power of the left ventricle, to study the recovery of heart transplant patients. Kung et al. (2014) developed a 0D model to simulate exercise in Fontan patients, an essential indicator of ontan patient survival rate. Both studies used models to predict patient outcomes; however, as noted by Colunga et al. (2020), reliable results require that only influential and identifiable model parameters are estimated given the model and available data. Parameter influence on model predictions can be quantified using local (Ellwein et al. , 2008; Olufsen & Ottesen, 2013) or global (Sumner et al. , 2012; Eck et al. , 2016) sensitivity analyses combined with subset selection. The latter technique identifies parameters that can be inferred uniquely given the model and available data. Schiavazzi et al. (2017) used sensitivity analysis in combination with optimization to estimate model parameters, comparing model predictions to data from single-ventricle patients with a Norwood physiology. Our group (Colunga et al. , 2020) used similar methods to analyze data from heart-transplant patients showing that model predictions align with static RHC data from single patient recordings and longitudinal patient data. The aforementioned studies used noninvasive and static data, but only a few have studied the importance of calibrating the model to the complete pressure waveforms. Examples include the study by Marquis et al. (2018), who calibrated a systems-level model to simultaneously recorded pressure and volume waveforms measured in the rat left ventricles. Other studies utilize one-dimensional (1D) fluid dynamics models (Colebank et al. , 2019; Chambers et al. , 2020) to match predictions of pulmonary pressure to time-series waveforms measured in-vivo in normotension and PH. Finally, Gerringer et al. (2018) matched model predictions from 3 and 4 element Wind-kessel models to main pulmonary artery (MPA) pressure waveforms in control and PAH mice. To compare the two models, they used the Akaike information criteria (AIC). They found that the simpler 3 element Windkessel model was preferred. These studies matched model simulations to data but only used a single MPA pressure waveform and were limited to predictions in the pulmo-nary vasculature, i.e., they could not study the effect on the heart or the systemic circulation. The only PH type that can be cured is CTEPH, but for patients with other disease types (e.g., PAH), further progression can be alleviated by medical intervention. PAH phenotypes include small artery and arteriolar muscularization and narrowing (Rol et al. , 2017). For these patients, disease management utilizes vasodilators, decreasing distal vasoconstriction and vascular wall thickening. PAH treatments traditionally target endothelial driven processes by (1) reducing the release of endothelin-1 (ET-1), a vasoconstrictor, via endothelin receptor antagonists; (2) increas-ing circulating levels of nitric oxide (NO), a vasodilator, using phosphodiesterase type 5 (PDE-5) inhibitors; and (3) increasing the release of prostacyclin, a vasodilator, by using prostanoids Humbert, 2010). More recently, soluble guanylate cyclase stimulators, such as riociguat, have successfully alleviated PAH by increasing the available cyclic guanosine monophosphate in the vascular wall inducing vasorelaxation (Ghofrani et al. , 2013 b ). CTEPH is initiated by vascular lesions and obstructions due to unresolved thromboembolic material and, similar to PAH, leads to small vessel disease (Gerges et al. , 2020). In contrast to PAH, CTEPH is considered nearly curable by surgical intervention (Siennicka et al. , 2019). The gold-standard treatment for CTEPH is pulmonary endarterectomy (PEA), which surgically removes embolic material. However, patients with distal lesions far from the MPA are considered inoperable. These patients receive pharmaceutical treatment, e.g., with riociguat (Ghofrani et al. , 2013 a ) and/or balloon pulmonary angioplasty (BPA). The ultimate goal of all available treatment protocols is to decrease distal PVR and reestablish high pulmonary vascular compliance to lower RV pressure. To study the effect of treatments, we adjust PVR and pulmonary compliance parameters, noting how changing these parameters change predictions for individual patients. In summary, this study uses a lumped parameter model, including the systemic and pulmonary circulations, the atria, and the ventricles, calibrated to static (systolic and diastolic) and dynamic (waveform) pressure data measured using RHC. We use local and global sensitivity analyses and subset selection to determine which parameters can be estimated given the model and available data. To quantify the benefits of including waveform data, we examine two residual vectors: com-paring model predictions to static data (systolic, diastolic, and mean pressures and cardiac output (CO)) and using a combination of static and dynamic data (RHC time-series). Results show that including dynamic waveform data in the residual allows us to estimate cardiac timing parameters, necessary for accurate prediction of pressure-volume dynamics in the heart. By integrating math-ematical modeling, patient-specific data, and physiological intuition, we categorize each patient's functional state, and by running simulations with estimated parameters, we calculate patient-spe-cific physiological biomarkers. We compare parameters across five PH patients to a normotensive control subject and simulate the effect of changing parameters in response to treatment strategies. METHODS Ethical Approval atient-specific data was obtained from two hospitals, adhering to their respective institutional review boards (IRB) guidelines. Deidentified RHC patient data were obtained from the Scottish Pulmonary Vascular Unit at Golden Jubilee National Hospital and from the Center from Pulmo-nary Vascular Disease at Duke University Medical Center.
Blood Pressure Data
This study utilizes RHC data from five patients with confirmed PH: two from group 1 (PAH) and three patients from group 4 (CTEPH). Two CTEPH datasets are from Duke University, and one CTEPH and the two PAH datasets are from the Scottish Pulmonary Vascular Unit at Golden Jubi-lee Hospital in Glasgow, Scotland. Static data include demographics (height, weight, sex, age), heart rate, systolic, diastolic, and mean systemic blood pressure measured with a sphygmomanom-eter. Patients underwent RHC, during which a catheter was advanced from the right atrium (RA) to the right ventricle (RV), and finally, the MPA, recording dynamic pressure waveforms at each location. The pulmonary arterial wedge pressure (PAWP, mmHg), an estimate of left atrial pres-sure, was captured by inflating a catheter balloon in a distal pulmonary artery and measuring the pressure downstream. CO (L/min) was measured during RHC by thermodilution. All pressure readings were obtained over 7-8 heartbeats. Static data are listed in Table 1, and dynamic pressure waveforms are shown in Fig. 1.
TABLE 1 FIGURE 1 Data Extraction
Time-series data are extracted from RHC reports using GraphClick version 3.0.3 for Mac OS. Beat-to-beat hemodynamic profiles for each patient are obtained by aligning each pressure wave-form to the electrocardiogram (ECG) signals measured simultaneously with the RHC. The time-series profiles are separated using the ECG R-R intervals and stored as separate data files. The analysis uses representative RA, RV, and MPA waveforms for each patient (shown in Fig. 1). The RHC data are not measured simultaneously. To ensure consistency, the representative waveforms are selected during expiration and assigned a cardiac cycle length equal to the averaged RA, RV, and MPA cycle for each patient. To align waveforms within the cardiac cycle, we shift the RA and PA signals such that RA contraction occurs before the start of RV isovolumic contraction and the peak RV pressure occurs immediately before peak MPA pressure. Magnitudes of the RA, RV, and MPA pressure signals are shifted slightly to ensure physiological valve dynamics. In addition to signal alignment, time-series data are shifted forward or backward in time to align with the computational model. We choose an optimal shift by minimizing the misfit between the model and data after selecting an identifiable, influential, and physiological parameter set. We simultaneously shift the RV and MPA data since their anatomical proximity supports nearly syn-chronous hemodynamics. After choosing the RV-MPA optimal shift, we carry out a similar mini-mization procedure using the RA time-series data and model predictions. We repeat this procedure for each patient and use this data when reporting optimal model fits and parameter estimates. Lastly, we construct a normotensive control subject using pressure and volume values from literature (Krohn Therkelsen et al. , 2006; Boron & Boulpaep, 2017); these pressure values are displayed in Table 2. Normotensive parameters and model predictions are compared to those ob-tained using PH data.
TABLE 2 Mathematical Model
This study develops a systems-level ordinary differential equations (ODE) compartment model (shown in Fig. 2), predicting dynamic pressure 𝑝 (mmHg), flow 𝑞 (mL/s), and volume 𝑉 (mL) waveforms in the systemic and pulmonary vasculature. The model consists of 8 compartments: four heart chambers (the left and right atria and ventricles), the systemic arteries and veins, and the pulmonary arteries and veins. The model is formulated using an electrical circuit analogy, with pressure analogous to voltage, flow to current, volume to charge, and compliance to capacitance. Equivalent to an RC-circuit, equations relating the three dependent quantities are given by 𝑑𝑉 !, 𝑑𝑡 = 𝑞 − 𝑞 (1) 𝑞 = 𝑝 − 𝑝 𝑅 (2) 𝑉 !, = 𝑉 −𝑉 ’(, = 𝐶 (𝑝 − 𝑝 ), (3) where the subscripts 𝑖 − 1, 𝑖, 𝑖 + 1 refer to the prior, current, and next compartments in the sys-tem. 𝑉 !, (mL) denotes the stressed volume (the circulating volume), 𝑉 ’(, (mL) the unstressed olume (the non-circulating volume, assumed constant), 𝑅 (mmHg s / mL) the resistance between two compartments, and 𝐶 (ml / mmHg) the compartment compliance. Equation (1) ensures con-servation of volume, equation (2) is the analog of Ohm’s law, and equation (3) relates volume and pressure. Flow in and out of heart compartments are controlled by valves modeled as diodes. The valves are either open or closed depending on the pressure gradient between compartments, i.e., 𝑞 = 1 𝑝 − 𝑝 𝑅 , valve open0, valve closed. The model includes four heart valves: two semilunar valves (the tricuspid and mitral valves) and two atrioventricular valves (the pulmonary and aortic valves). Heart dynamics are modeled using a time-varying function (Ellwein et al. , 2008; Marquis et al. , 2018), relating the pressure in each chamber to a time elastance 𝐸 (𝑡) (mmHg / mL) and volume by 𝑝 (𝑡) = 𝐸 (𝑡̃)𝑉 !, , (4) where 𝑖 = 𝑟𝑎, 𝑙𝑎, 𝑟𝑣, 𝑙𝑣 denote the left (𝑙) and right ( 𝑟) atria ( 𝑎 ) and ventricle ( 𝑣 ) (Note: elastance is the inverse of compliance, i.e., 𝐸 = 1/𝐶) . The elastance function is defined over one cardiac cycle, i.e., time 𝑡̃ = mod (𝑡, 𝑇), where
𝑇 (s) is the length of the cardiac cycle. The elastance func-tion 𝐸 ) (𝑡̃) determining ventricular contraction over a cardiac cycle ( ) is described by the piecewise continuous function (Ellwein et al. , 2008) 𝐸 ) (𝑡) = ⎩⎪⎨⎪⎧ (𝐸 ),* − 𝐸 ),+ )2 M1 − cos M𝜋𝑡̃ 𝑇 ),, OO + 𝐸 ),+ , 0 ≤ 𝑡̃ ≤ 𝑇 ),, (𝐸 ),* − 𝐸 ),+ )2 M1 + cos M 𝜋P𝑡̃ − 𝑇 ),, Q(𝑇 ),- − 𝑇 ),, )OO + 𝐸 ),+ , 𝑇 ),, ≤ 𝑡̃ ≤ 𝑇 ),- 𝐸 ),+ , 𝑇 ),- ≤ 𝑡̃ ≤ 𝑇 (5) where 𝐸 ),+ and 𝐸 ),* (mmHg / mL) are the minimal and maximal ventricular elastances, and 𝑇 ),, (s) and 𝑇 ),- (s) are the durations of ventricular contraction and relaxation, respectively. The atrial elastance function is described in a similar fashion (Liang et al. , 2009) . (𝑡̃) = ⎩⎪⎪⎪⎨⎪⎪⎪⎧𝐸 .,* − 𝐸 .,+ .,- Q(𝑇 − 𝑇 .,, + 𝑇 .,-
OO + 𝐸 .,+ , 0 ≤ 𝑡̃ ≤ 𝑇 .,- 𝐸 .,+ , 𝑇 .,- ≤ 𝑡̃ ≤ 𝜏 .,, 𝐸 .,* − 𝐸 .,+ .,, Q𝑇 .,, − 𝜏 .,, OO + 𝐸 .,+ , 𝜏 .,, ≤ 𝑡̃ ≤ 𝑇 .,, 𝐸 .,* − 𝐸 .,+ .,, Q𝑇 − 𝑇 .,, + 𝑇 .,-
OO + 𝐸 .,+ 𝑇 .,, ≤ 𝑡̃ ≤ 𝑇. (6) Here, 𝐸 .,+ and 𝐸 .,* (mmHg / mL) denote the minimum and maximum elastances of the atria, and 𝑇 .,- , 𝜏 .,, and 𝑇 .,, (s) denote the start of atrial relaxation ( 𝐸(𝑡̃) = 𝐸 .,+ ), the start of atrial contrac-tion, and the point of maximum atrial contraction, parameterized such that .,- ≤ 𝜏 .,, ≤𝑇 .,, ≤ 𝑇 . For each cardiac cycle ( , ventricular elastance 𝐸 ) (𝑡̃) is initialized at the be-ginning of the isovolumic contraction. In contrast, atrial elastance 𝐸 . (𝑡̃) decreases to 𝐸 .,+ during isovolumic contraction and then increases before ventricular contraction, consistent with cardiac physiology (Boron & Boulpaep, 2017). The time course of atrial and ventricle elastance is shown in Fig. 3. FIGURES 2 and 3 Model Outcomes
In addition to hemodynamic predictions of pressure, flow, and volume, we compute five physio-logical indices consistent with clinical understanding of PH progression. Several of these indices can be measured in the clinic, but we propose new indices that can only be determines using a computational modeling approach. i.
Cardiac power per cycle (CPW) is defined as the time-averaged integral of the PV loop,
𝐶𝑃𝑊 = %/ ∫ 𝑝(𝑡)𝑑𝑉′ , calculated in each heart chamber (Rimehaug et al. , 2013; Colunga et al. , 2020). ii. Resistance ratio of pulmonary and systemic resistance is defined as 𝑅 /𝑅 ! (Yang et al. , 2018). iii. Compliance ratio of pulmonary and systemic compliance is defined as 𝐶 /𝐶 !. . iv. Ventricular elastance ratio of minimal elastance in the right and left ventricles is defined as 𝐸 -),+ /𝐸 . . Pulsatility index (PI) is computed as the ratio of pulmonary arterial pulse pressure to average right atrial pressure, P𝑝 − 𝑝 Q/𝑝̅ -. (Mazimba et al. , 2019). Parameter Values and Initial Conditions
The combination of sparse data and a large number of parameters makes it imperative that nominal parameter values and initial conditions are physiological and patient-specific. Following the ap-proach described in (Marquis et al. , 2018; Colunga et al. , 2020), we use a combination of patient-specific data (where available) and values reported in the literature. Table A1 in the Appendix lists the nominal parameter values and their calculation.
Compartment volumes and cardiac output . Using Hidalgo’s formula (Hidalgo et al. , 1962), for each patient, the total blood volume ( 𝑉 , mL) is calculated as a function of height ( 𝐻, cm), weight ( 𝑊 , kg), and sex (Williams et al. , 2019) as 𝑉 = Y (3.47 ⋅ 𝐵𝑆𝐴 − 1.954) ⋅ 1000, 𝑖𝑓 Female(3.29 ⋅ 𝐵𝑆𝐴 − 1.229) ⋅ 1000, 𝑖𝑓 Male, (7) where 𝐵𝑆𝐴 = g𝑊 ⋅ 𝐻/3600 ( m ) is the body surface area (DeCherney & Berkowitz, 1982). For each compartment, the initial stressed volume (initial conditions for the states) is calculated as a proportion of the total volume using previously published data (Beneken & DeWit, 1966). The systemic arterial and venous volumes are approximately 13% and 65% of 𝑉 . The corre-sponding pulmonary arterial and venous volumes are 3% and 11% of 𝑉 . Of these, approximately 27% and 7.5% of the systemic arterial and venous volume are stressed, while 60% and 11% of the pulmonary arterial and venous volume are stressed. These distributions cannot be measured but are taken from previous studies (Ellwein et al. , 2008; Marquis et al. , 2018). For the heart, we assume that the atrial and ventricular volumes are 1.5% and 2.5% of the total blood volume and that the unstressed volume is 5 and 10 ml, respectively. CO is calculated, assuming that the entire blood volume circulates in one minute (Ellwein et al. , 2008; Boron & Boulpaep, 2017). An ejection fraction of 60% is assumed in the ventricles, while the atrial ejection fraction is 47% (Lin et al. , 2008). Pressure.
Nominal values for the pulmonary circulation pressure are extracted from the RHC data, and nominal systolic and diastolic pressures in the systemic arteries are taken from cuff measure-ments. These values are listed in Table 2. Nominal pressure values for compartments from which e do not have measurements (the left atrium, left ventricle, and systemic veins) are calculated by scaling pressures in their adjacent compartments where data is available (Williams et al. , 2019). For the left atrium, we assume that the pulse pressure is approximately 5 mmHg, consistent with previous studies (Pironet et al. , 2013). These considerations can be summarized as 𝑝 !) = 1.01 𝑝 -),+ , (8) 𝑝 = 0.95 𝑝 , (9) 𝑝 = 𝑝 + 5, (10) 𝑝 = 0.97 𝑝 , (11) 𝑝 = 1.01 𝑝 !.,* . (12) Resistance.
Each compartment is separated by a resistance to flow. Utilizing Ohm’s law, the nom-inal vascular resistance is calculated as 𝑅 = 𝑝 − 𝑝 CO , (13) where the resistance in compartment 𝑖 depends on the pressure in the current and previous com-partment and the CO. The heart valves’ resistances are calculated in a similar fashion; the aortic and pulmonary valve resistances are calculated as 𝑅 .). = 𝑝 − 𝑝 !.,* CO and 𝑅 = 𝑝 -),* − 𝑝 CO (14) whereas the mitral valve resistance is calculated as 𝑅 +). = 𝑝 − 𝑝 CO (15) to ensure the left atrium drains into the left ventricle during diastole. In the case of PH, RA pressure is elevated (Alenezi et al. , 2020), i.e., equation (15) overestimates valve resistance. To circumvent this, we set 𝑅 = 0.055 for all five PH patients. Compliance is defined as the relative change in volume for a given change in pressure (Wang et al. , 2013). It quantifies the ability of the vasculature to distend under load. In this study, nominal compliance estimates are 𝐶 = 𝑉 𝑝 and 𝐶 = 𝑉 𝑝̅ (16) for the arterial and venous compartments, respectively. eart parameters include elastance and timing parameters. Noting that compliance is the inverse of elastance and that the compliance in the heart is minimal during end-systole (computed at the maximum pressure and minimal volume) (Marquis et al. , 2018), we calculate the maximum and minimum elastances as 𝐸 = 𝑝 𝑉 and 𝐸 = 𝑝 𝑉 , (17) respectively, where 𝑖 = 𝑙𝑎, 𝑟𝑎, 𝑙𝑣, 𝑟𝑣 . Nominal timing parameters for the RA and RV elastance functions are calculated using the time-series data. Specifically, the maximum and minimum RV elastance occur at peak systole and the beginning of diastole, at time 𝑡 = 𝑇 ),, and 𝑡 = 𝑇 ),- , respectively. For the atrium, the end of atrial systole, the start of atrial contraction, and peak atrial contraction are used to calculate the timing parameters 𝜏 .,- , 𝑇 .,, and 𝑇 .,- . Since dynamic data is not available for the left atrium and ventricle, we define the left-heart chamber timing parameters relative to the right-heart timing parameters (Boron & Boulpaep, 2017). Summary
The model is a system of eight ODE’s (one per compartment) with 25 parameters. The rate of change for the eight stressed volumes, 𝑉 !, , for each compartment are computed as 𝒚 = 𝑔(𝑡, 𝒙; 𝜽), 𝑑𝒙𝑑𝑡 = 𝑓(𝑡, 𝒙; 𝜽), 𝒙 = p𝑉 , 𝑉 , 𝑉 !. , 𝑉 !) , 𝑉 -. , 𝑉 -) , 𝑉 , 𝑉 q 𝜽 = {𝑅 ! , 𝑅 , 𝑅 .). , 𝑅 +). , 𝑅 , 𝑅 , 𝑅 , 𝑅 !) , 𝐶 !. , 𝐶 !) , 𝐶 , 𝐶 , 𝑇 .,- , 𝜏 .,, , 𝑇 .,, , 𝑇 ),, , 𝑇 ),- 𝐸 , 𝐸 , 𝐸 -.,* , 𝐸 -.,+ , 𝐸 , 𝐸 , 𝐸 -),* , 𝐸 -),+ } (18) where 𝒙 are the state variables ( 𝑉 !, in compartment 𝑖 ; we drop the subscript 𝑠 for simplicity). The functions 𝑓(𝑡, 𝑥; 𝜽) denote the states’ evolution (see equation (1)), and 𝜽 denote the parameters. The vector 𝒚 includes pressure and CO predictions used for parameter inference; both pressures 𝑝 and flows 𝑞 are computed from volumes 𝑉 via equations (2, 3). arameter Inference To determine which biomarkers change during PH progression, we estimate model parameters by minimizing the relative least-squares error between model predictions and data. This study uses the Levenberg-Marquardt method to solve the generalized least-squares problem (Kelley, 1999). In general, the observed data 𝒚 (static or time-series) is assumed to be of the form 𝒚 = 𝑔(𝑡, 𝒙; 𝜽) + 𝜺, (19) where 𝑔(𝑡, 𝒙; 𝜽) is the model predictions (here, pressure and CO), and 𝜺 is the measurement error, assumed to be independent and identically distributed (iid) white Gaussian noise, i.e., 𝜺 ~ 𝒩(0, 𝜎 𝐈) . Using this framework, we minimize the relative sum of squared errors, 𝐽 = 𝒓 / 𝒓 , where 𝒓 is the residual vector, encompassing the relative differences between the measured data 𝒚 and the model predictions 𝒚 = 𝑔(𝑡, 𝒙; 𝜽) . To understand how different data improve the inference procedure, we compute residual vec-tors with combinations of static and dynamic RHC data. The static residual is defined as 𝒓 ! = 1g𝑁 ! 𝒚 − 𝒚 𝒚 , (20) 𝒚 = ~𝑝 -.,* , 𝑝 -.,+ , 𝑝 -),* , 𝑝 -),+ , 𝑝 , 𝑝 , 𝑝 !.,* , 𝑝 !.,+ , 𝑝 , CO(cid:127) where the vector 𝒚 contains maximum and minimum pressure predictions, as well as CO, 𝒚 the corresponding data, and 𝑁 ! is the number of points in the static residual. The three dynamic resid-uals accounting for dynamic variation of the cardiac cycle are given by 𝒓 -. = 1g𝑁 -. 𝒑 -. (𝑡) − 𝒑 -.7 (𝑡)𝒑 -.7 (𝑡) (21) 𝒓 -) = 1g𝑁 -) 𝒑 -) (𝑡) − 𝒑 -)7 (𝑡)𝒑 -)7 (𝑡) (22) 𝒓 = 1g𝑁 𝒑 (𝑡) − 𝒑 (𝑡)𝒑 (𝑡) (23) where 𝒑 (𝑡) , 𝒑 (𝑡), and 𝑁 are the pressure predictions, pressure data, and number of residual points, respectively, for the RA ( 𝑟𝑎 ), RV ( 𝑟𝑣 ), and pulmonary arteries ( 𝑝𝑎) . Utilizing the above, we consider two combined residuals as our quantity of interests: 𝒓 = 𝒓 ! , 𝒓 = ~𝒓 ! , 𝒓 -. , 𝒓 -) , 𝒓 (cid:127). or the atrial residual 𝒓 -. , we include data before and after peak contraction at 𝑡 = 𝑇 .,, , as well as the first and last time point. We do this to ensure that parameter estimation is not biased by the model’s inability to capture the entire atrial data. As seen in previous studies (Ellwein et al. , 2008; Liang et al. , 2009), elastance models for the atria typically cannot capture all the atrial dynamics seen in-vivo (Pironet et al. , 2013). Similar to the approach in (Marquis et al. , 2018), each residual is computed over the last 30 cycles of the model predictions compared to the data. Sensitivity Analysis
To determine parameters that affect model predictions, we conduct sensitivity analyses with re-spect to the residual vectors 𝒓 and 𝒓 . We utilize both local, derivative-based, and global, vari-ance-based sensitivity analyses. The former methods are valid at the nominal parameter values and quantify the gradient of the residual vectors with respect to the parameters. The latter measure model sensitivity throughout the physiological parameter space and vary multiple factors at a time. The local sensitivities of the residuals, 𝝌 , are used to rank parameters based on their influence (Marquis et al. , 2018; Colunga et al. , 2020). Global sensitivity analysis (GSA) provides first ( 𝑆 ) and total order ( 𝑆 / ) Sobol’ indices (Sobol, 2001); the former index measures the parameters’ con-tribution to the total output variance of the residual, whereas the latter quantifies all higher-order interactions between the parameters on the residual variance. The total order Sobol’ indices 𝑆 / are used to order parameters from most to least influential. Detailed derivations of both local and global methods can be found in Appendices 1 and 2. Parameter Subset Selection
Given the limited data and the large number of parameters in the ODE model summarized in equa-tion (18), not all parameters are identifiable, i.e., the parameters that best fit the data cannot be estimated uniquely. Moreover, the model is formulated from electrical circuit theory. In such mod-els, resistors and capacitors in series and parallel cannot be estimated uniquely unless data is avail-able adjacent to the circuit component. A subset of the most influential parameters is identified using the sensitivity analyses, insights from the underlying circuit theory, and physiological intuition. The latter is vital to estimate pa-tient-specific biomarkers (Pope et al. , 2008; Colunga et al. , 2020). We take several steps to ensure that the parameter subsets are identifiable and influential. ocal sensitivity-based subset selection . Local sensitivity analysis provides information about which parameters (at their nominal values) are most influential with respect to the residuals. As-ymptotic analyses (Cintrón-Arias et al. , 2009) show that these sensitivities construct a linear ap-proximation of the Fisher information matrix , 𝑭 = 𝝌 / 𝝌 , from which the correlation matrix of the parameters 𝒄 can be calculated as 𝑐 = 𝐶 g𝐶 𝐶 ;; , 𝑪 = 𝜎 𝑭 $% , (24) where −1 ≤ 𝑐 ≤ 1 for 𝑖, 𝑗 = 1. . . , ℳ , and 𝜎 <6 is the noise variance as before. The correlation matrix 𝑐 provides information about pairwise relationships between parameters; a high correla-tion ( |𝑐 | → 1 ) implies that parameters 𝑖 and 𝑗 cannot be estimated simultaneously (Ellwein et al. , 2008; Pope et al. , 2008; Marquis et al. , 2018). In this study, we denote parameters for which (cid:138)𝑐 (cid:138) ≥ 0.95 correlated. For each correlated parameter pair, we fix the least influential parameter at its nominal value. Global sensitivity-based subset selection.
The first-order effects 𝑆 from the GSA describe a single parameter’s influence on the residuals’ variance. Since the Sobol’ indices are defined in terms of variance, several properties hold. By definition , 𝛴 𝑆 ≤ 1 , hence the difference 𝑆 de-scribes the variance attributed to higher-order effects. Second, the difference between the total and first-order effects, 𝑆 / ! − 𝑆 , quantifies the proportion of higher-order (interaction) effects at-tributed to the parameter, with 𝑆 = 𝑆 / ! suggesting negligible higher-order effects. Lastly, if 𝑆 / ! =0, then by definition parameter 𝑖 does not affect the variance of the quantity of interest and can be fixed before parameter inference (Sumner et al. , 2012; Eck et al. , 2016) . Model Selection
As noted in our previous section, it is possible to determine several identifiable subsets. For each subset, the ℳ model parameters are split into two groups 𝜽 = p𝜽 > , 𝜽 ℳ$> q, where parameters in 𝜽 > are identifiable and parameters in 𝜽 ℳ$> are not identifiable and kept fixed at their nominal value, i.e., each subset gives rise to a separate model biased by the
ℳ − 𝜌 fixed parameters. Sim-ilar to previous studies (Gerringer et al. , 2018; Qureshi et al. , 2019), we compute two information criteria, balancing the goodness-of-fit and the model complexity, to select the best parameter ubset: the corrected Akaike information criteria (AICc) and the Bayesian information criteria (BIC). These are given by
𝐴𝐼𝐶𝑐 = 2 log(𝐽) + 2𝜌 + 2 𝜌(𝜌 + 1)𝑁 − 𝜌 − 1 , (25)
𝐵𝐼𝐶 = 2 log(𝐽) + 𝜌 log(𝑁). (26) where 𝜌 ≤ ℳ is the number of parameters in the subset, 𝑁 is the number of data points in the residual, and 𝐽 is the least-squares cost. Confidence and Prediction Intervals
Both model parameters and predictions come with some level of uncertainty, which can be quan-tified using an asymptotic analysis (Colebank et al. , 2019). Under the assumption of iid errors 𝜺 , we construct the variance estimator 𝜎(cid:144) <6 and parameter covariance estimator 𝑽(cid:146) using asymptotic analysis for weighted nonlinear least-squares. This derivation is included in Appendix 3 . For the estimated model parameters 𝜽(cid:146) > , the parameter confidence intervals for each 𝜃(cid:149) are computed as ~𝜃(cid:149) , 𝜃(cid:149) (cid:127) = 𝜃(cid:149) ± t A$>B.DEF (cid:152)𝑽(cid:146) , (27) where t A$ℳ " %$G/6 is a two-sided t-statistic with a confidence level, and (cid:152)𝑽(cid:146) repre-sents the standard error for the 𝑖′ th parameter estimator. Similarly, the confidence and prediction intervals for the optimal model output 𝑦(cid:144) ; at time 𝑡 ; are given by ~𝑦(cid:144) ;?@$ , 𝑦(cid:144) ;?@& (cid:127) = 𝑦(cid:144) ; ± t A$>B.DEF (cid:152)𝝌 ;I 𝑽(cid:146) 𝝌 ; (28) ~𝑦(cid:144) ;J@$ , 𝑦(cid:144) ;J@& (cid:127) = 𝑦(cid:144) ; ± t A$>B.DEF (cid:152)𝜎 <6 + 𝝌 ;I 𝑽(cid:146) 𝝌 ; , (29) where 𝝌 ;I is the sensitivity vector at 𝑡 ; evaluated at 𝜽(cid:146) = p𝜽(cid:146) 𝝆 , 𝜽 ℳ$𝝆 q. Note that the prediction in-tervals account for the variance in both the model output and the data; hence they are wider. The small sample size for 𝒓 % limits the analysis of this residual. Therefore, we perform our uncertainty quantification using 𝒓 as the quantity of interest. reatments RV dysfunction in PH patients directly affects mPAP and the elastic properties of the arterial tree (Tabima et al. , 2017) and is regulated by both PVR and pulmonary arterial compliance. To simu-late different treatment procedures for PH patients, we vary the PVR and pulmonary arterial com-pliance. These are the ultimate targets for PH treatment for PAH and CTEPH patients. In particu-lar, we analyze the effects of two drug classes, Phosphodiesterase-5 (PDE-5) inhibitors, and solu-ble guanylate cyclase stimulators, as well as two surgical interventions, PEA and BPA, used for CTEPH patients. Since the effects of these treatments vary with patient and disease severity, we simulate several treatment intensities. Table 3 provides a summary of the effects of each treatment.
PAH specific treatments . PDE-5 Inhibitors reduce pulmonary and systemic resistance by inducing vasodilation (Lindman et al. , 2012), increasing pulmonary vascular compliance to its normoten-sive value. Systemic and pulmonic vasodilation can be simulated by reducing parameters associ-ated with the resistance in the systemic and pulmonary arteries, 𝑅 ! and 𝑅 , while simultaneously increasing compliance, 𝐶 !. and 𝐶 . CTEPH specific treatments . PEA involves surgical removal of lesions in the proximal pulmonary arteries and is deemed the gold standard treatment for operable CTEPH (Siennicka et al. , 2019). For inoperable CTEPH patients, BPA is an alternative treatment. BPA disrupts pulmonary artery lesions by inflating a balloon catheter and leads to a reduced obstruction to blood flow. Both of these procedures affect the pulmonary arteries by reducing the impedance to flow. We simulate these treatments by reducing the resistance in the pulmonary arteries, 𝑅 . Long-term effects of clot removal, including normalized pulmonary vascular compliance, is simulated simultaneously by increasing the compliance in the pulmonary arteries 𝐶 . PAH and CTEPH treatment . Riociguat, a soluble guanylate cyclase stimulator, is the only FDA-approved therapeutic drug for PAH and CTEPH (Ghofrani et al. , 2013 a , 2013 b ). Taken orally, riociguat targets pulmonary vascular smooth muscle cells leading to vasodilation and increased lumenal area, reducing both mPAP and the mean systemic arterial pressure (mSAP) (Ghofrani et al. , 2013 b ). Like PDE-5 inhibitors, the relaxing and widening of the pulmonary arteries and re-duction of mean systemic artery pressure can be simulated by reducing 𝑅 ! and 𝑅 , while increasing 𝐶 !. and 𝐶 . ABLE 3 Simulations
To study the impact of PH, we compare the simulations and outcomes for normotensive controls and PH patients and analyze the effects of PH interventions.
Control:
Simulations for a normotensive control subject are conducted using pressure values listed in Table 2. Hemodynamic predictions of pressure, flow, and volumes are compared to those ob-tained for PH patients.
Static:
Similar to Colunga et al. (2020), we simulate and match model predictions utilizing only static pressure data for each PH patient. For these simulations, we infer parameters specific to PH patients using residual 𝒓 % . Simulation results are used as a benchmark to determine whether the addition of dynamic waveforms results in better model predictions. Dynamic waveforms:
Our model is also calibrated using static and dynamic waveform data. RHC recordings are not measured simultaneously. Therefore, waveforms are shifted, both forward and backward, ensuring that predictions are in phase with time-series profiles. Predictions of systolic, diastolic, and mean pressure are computed in combination with dynamic RA, RV, and pulmonary artery predictions. These predictions are matched to shifted data that best aligns to the model, utilizing the residual 𝒓 . Treatments:
After deducing the best parameter subsets and residual vector, we simulate improve-ments in hemodynamics in response to PH treatments. Inferred resistance and compliance param-eters are altered following treatment interventions summarized in Table 3. Treatment predictions are computed for all five patients and subsequently compared to predictions for a normotensive control subject.
RESULTS
Local and global sensitivity analyses are applied to both residuals 𝒓 and 𝒓 distinguish influential and noninfluential parameters. The noninfluential parameters are fixed at their nominal values, and correlation analysis is used to construct subsets of influential, identifiable parameters. The best subset for each residual is determined using information criteria and subsequently used for param-eter inference. Waveforms are shifted to obtain the model predictions best fitting the measured ata. Finally, we contrast model predictions between normotensive and PH patients and simulate improvements in pulmonary and systemic arterial pressure in response to treatments. Sensitivity Analysis
Figure 4a-b shows the patient-specific local sensitivity parameter ranking for 𝒓 % (static values, panel (a)) and 𝒓 (static and time-series data, panel (b)). The two-norm of the residual sensitivities is divided by the largest sensitivity. The compliance, resistance, and elastance parameters near the right heart are the most influential for all five patients. Results show that by accounting for dynamic predictions encoded in residual 𝒓 , the timing parameters 𝑇 ,,-) , 𝑇 -,-) and 𝜏 ,,-. become influential. Four of the five patients display consistent pa-rameter ranking for both residual vectors. The exception is patient 3, for which 𝑇 ,,-) and 𝑇 -,-) are less influential. Overall, the parameters 𝜽 AL = ~𝑅 .). , 𝑅 +). , 𝑅 , 𝑅 , 𝑅 !) (cid:127) are noninfluential for both residuals, suggesting that these parameters can be kept fixed at their nominal values, i.e., they need not be included in parameter subset and model selection. For the GSA, 𝓃 = 10 M samples are generated for each parameter using a Sobol’ sequence, ensuring adequate parameter space coverage. First-order and total effects are shown in Fig. 4c-d for the cost functional 𝐽 (𝜽), 𝑖 = 1,2 calculated using residual 𝒓 % and 𝒓 , respectively. The total Sobol’ indices 𝑆 / ! are almost zero for parameter subset 𝜽 AL , consistent with the local sensitivity results. The most influential parameters for 𝒓 % are 𝑅 ! , 𝐶 !) , and 𝑅 . The next group includes 𝑅 , 𝐶 !. , 𝐶 , and the ventricular elastance parameters. In contrast, for 𝒓 , the timing parameters 𝑇 ),, and 𝑇 ),- are the most influential, and the pulmonary resistance 𝑅 and systemic vein compli-ance 𝐶 !) are the next most influential parameters. FIGURE 4
Parameter Subsets, Model Selection, and Parameter Inference
To identify parameters that can be inferred, we construct the correlation matrix defined in equation (30). The noninfluential parameters 𝜽 AL are not considered for this analysis. Parameter pairs with a correlation greater than 𝛾 = 0.95 are denoted correlated. Following Olufsen (2013), the least nfluential parameter satisfying this condition is fixed at its nominal value and removed from the subset. This procedure is repeated until the subset contains no correlated parameter pairs. Since the subset selection procedure leads to several sets of influential, identifiable parameters, we need additional criteria to deduce the best parameter subset. To do so, we infer parameters for one patient (patient 5) using all subsets, recording the minimal least-squares cost for each. We subsequently compute AICc and BIC values for each subset. Table 4 lists the parameter subsets along with the cost, AICc, and BIC scores. In general, the best model has the smallest AICc or BIC; hence, we select the subset that gives the smallest AICc and BIC values, 𝜽 %- and 𝜽 M- $ , for 𝒓 % and 𝒓 , respectively. After selecting optimal subsets, we determine optimal shifts to align the model predictions with the dynamic data. Figure 5a shows changes in the least-squares error for different data shifts for patient 5; the optimal shift is the minimum in the curve. This is done for the RV and PA first, and then successively for the RA. Figure 5b shows the data shifts and the model predictions. The optimal shifts are used when inferring the parameter subsets 𝜽 %- and 𝜽 M- $ for each patient. Table 5 shows the estimated parameters in 𝜽(cid:146) M- $ for each patient, along with 95% parameter confidence intervals calculated using eq. (34) (expressed as 𝜃(cid:149) ± 2𝜎 N ! ). Optimal parameter estimates 𝜽(cid:146) %- are listed in Table A2. TABLES 4 AND 5 FIGURE 5
Hemodynamic Predictions
To understand how parameters change with PH, we display the inferred patient parameters in Fig. 6 as the relative change from the normotensive parameters in box-and-whisker plots. Note that parameters shared between 𝜽 %- and 𝜽 M- $ are nearly identical even with additional parameters in 𝜽 M- $ . Post-optimization predictions of pressure and CO using either 𝒓 % or 𝒓 are depicted in Fig. 7 along with the measured data from patient 5. Predictions for all PH patients are shown in Figure A1. Inference using 𝒓 minimizes the mismatch between the model predictions and static and time-series data. RA and RV dynamics improve significantly when including time-series data. In con-trast, PA predictions improve only marginally; the diastolic decay in PA pressure occurs more uickly than the PA data. CO predictions are not as well-matched when using 𝒓 instead of 𝒓 % , but maximum and minimum pressure values still match the data well. A benefit of computational models is that essential but unmeasurable outcomes can be pre-dicted, such as pressure-volume (PV) loops. We contrast PV loops from all four heart chambers for the normotensive subject and the five PH patients (using inferred parameters from 𝒓 ) in Fig. 8. Except for Patient 5, all PH patients have reduced left atrial pressure, while the RA pressure is increased in 2 of 5 PH patients. Moreover, the maximum RA volume is significantly higher in PH patients. In the normotensive simulations, both ventricles display a rapid increase in pressure dur-ing systole with a steady decrease in volume. Eventually, the pressure decreases before the onset of diastole. In contrast, ventricular PV loops from the PH patients have a constant increase in pressure throughout systole. Also, the RV PV loops are similar in shape to the left ventricular PV loops in PH, while the RV PV loop in normotensive conditions has a more district shape compared to the left ventricle. These PV loops allow us to calculate CPW (the integral of the PV loop) for all four heart chambers listed in Table 6. We also contrast other model outcomes, including the resistance, compliance, and minimum elastance ratios, 𝑅 /𝑅 ! , 𝐶 /𝐶 !. , 𝐸 -),+ /𝐸 , and the pul-satility index PI . LA CPW is lower in PH for all but patient 5, while the RA CPW is lower in all CTEPH patients (3, 4, and 5). Left ventricle CPW varies across PH patients, but RV CPW is in-creased in all five PH patients. Finally, as expected 𝑅 /𝑅 ! , 𝐸 -),+ /𝐸 , and PI are increased, while 𝐶 /𝐶 !. is decreased in PH. The optimal parameters 𝜽(cid:146) M- $ are used to predict and contrast pressure and CO predictions be-tween the normotensive control and patient 5 in Fig. 9a. In addition, optimal parameters are used to construct both parameter confidence intervals (listed in Table 5) and model confidence and prediction intervals, shown in Fig. 9b. The confidence and prediction intervals show uncertainty in mean pulmonary venous pressure (matched to PAWP data), CO, maximum and systemic artery pressures, and uncertainty intervals for the time series predictions in the RA, RV, and PA. The uncertainty of systemic artery pressure and CO predictions is smaller than that of pulmonary ve-nous pressure predictions; the prediction intervals for the former two vary only ±5% from the predicted value, while pulmonary venous pressure prediction intervals have an uncertainty of ±20 %. The RA and RV data nearly all fall within the 95% prediction intervals; however, the diastolic portion of the PA time-series data does not fall within the 95% prediction intervals. Note hat RA uncertainty is only calculated during the systolic phase, which corresponds to which data are used in our parameter inference. TABLES 5 AND 6 FIGURES 6-9
Treatment Interventions
Patients with PH have increased PVR and decreased arterial compliance. Clinical studies indicate that specific treatment strategies decrease the resistance of the pulmonary and systemic vasculature and increase pulmonary and systemic artery compliance. To study the effects of administering treatment, the parameters 𝑅 , 𝐶 , 𝑅 ! and 𝐶 !. are varied one at a time for each of the five patients. 𝑅 and 𝑅 ! are gradually decreased in 10% intervals from 0% to 80%, while 𝐶 and 𝐶 !. are in-creased. Lastly, we combine increases in resistance and decreases in compliance. Results in Fig. A2 show that 𝑅 and 𝑅 ! are more influential on mPAP and mSAP than 𝐶 and 𝐶 !. . We conduct a more thorough investigation into the effects of specific PH treatments, analyzing the effect of vasodilators and surgical intervention with or without the combination of a vasodilator. Specific treatments based on each patient’s etiology are summarized in Table 3; patients 1-5 (PAH and CTEPH) receive vasodilators, and patients 3-5 (CTEPH) receive surgical intervention with and without drug therapy. Figure 10 illustrates each simulated treatments’ effect on the mPAP, systolic, and diastolic systemic artery pressures. All treatment simulations induce a de-crease in mPAP in all five patients. In particular, treatments T5 and T8, related to surgical inter-ventions for CTEPH, result in a reduction of mPAP below the PH threshold of 20 mmHg for patients three and four. Additionally, SA systolic pressures decreased in all patients for T1-4 and T8-10 while increasing in T5-7; diastolic SA pressure follows a similar pattern. FIGURE 10
DISCUSSION
This study investigates improvements in parameter inference when combining static and time-series data in a 0D cardiovascular model. Using a combination of sensitivity analyses, correlation analysis, and information criteria, we determine parameter subsets that best fit the data without ncluding non-informative parameters. We shift time-series data and fit the model to time-depend-ent hemodynamic waveforms, revealing changes in both model predictions and parameters due to PH. Model outcomes, including
𝐶𝑃𝑊, 𝑃𝐼 , compliance ratios, and resistance ratios, change in PH, consistent with physiological understanding of PH. We simulate PH treatment in each patient, noting changes in mPAP and mSAP. Our results show that surgical interventions for CTEPH pa-tients can lower mPAP below 20 mmHg, while vasodilator treatment alone for all five PH patients cannot.
Sensitivity Analysis
Sensitivity analysis is crucial for determining which parameters are influential on the model out-puts or outcomes. The model in this study has 25 parameters, yet limited data makes inferring all these parameters infeasible. Both the local analysis and GSA reveal a noninfluential set of param-eters 𝜽 AL that can be fixed. Within this set are the mitral and aortic valve resistances, which agree with our previous study (Colunga et al. , 2020). We speculate that if data are measured in the left ventricle, these parameters likely become influential. This could be important when studying group 2 PH attributed to left-heart failure (Philip et al. , 2019). Local sensitivity analysis results depend on the nominal parameter values used. As shown in Fig. 4a-b, ventricular timing parameters 𝑇 ),, and 𝑇 ),- are less influential for patient 3 than the other PH patients. Marquis et al. (2018) reported similar results for their 0D model. These discrepancies are circumvented when using GSA, which quantifies the model sensitivity throughout the param-eter space. Results in Fig. 4c-d reveal that the ventricular timing parameters are most influential for 𝒓 , supporting the local analysis results from patients 1, 2, 4, and 5. GSA results suggest that 𝑅 ! , 𝐶 !) , 𝑅 , 𝑅 , and 𝐸 +,3) are the most influential with static predic-tions used in 𝒓 % . This result agrees with findings by Colunga et al. (2020), who showed that elas-tance of the systemic veins (here, 𝐶 !) ) and 𝑅 ! are influential for predicting systolic and diastolic pressure values in the systemic and pulmonary arteries and veins, as well as the RV. These two parameters, along with 𝐸 +,3) , determine systemic circulation pressure ranges in the model. The parameters 𝑅 and 𝑅 determine the interactions between the right heart and MPA. The tricuspid valve resistance is determined from RA-RV interactions, while elevated pulmonary resistance, 𝑅 , increases RV afterload. Results using 𝒓 suggest that parameters in the RV and PA are most influ-ential. Physiologically, previous studies show that RV function degrades in PH, captured here by hanging 𝑇 ),, and 𝑇 ),- . Increased PVR and decreased PA compliance are common phenotypes of PH (Wang & Chesler, 2011; Kheyfets et al. , 2013) and are controlled by 𝑅 and 𝐶 . GSA ranking in Fig. 4c-d supports this and shows that these two parameters are more influential than a majority of other parameters. Deficiencies in RA reservoir function and active contraction dynamics are strong predictors of mortality in PH (Alenezi et al. , 2018). The RA filling is dictated by systemic vein dynamics and tricuspid valve integrity, supporting the fact that 𝐶 !) and 𝑅 are relatively more influential than other parameters. In the model, RA contraction is dictated by minimum elastance 𝐸 -.,+ , another influential parameter based on the GSA results. Interestingly, our results suggest that pulmonary valve resistance ( 𝑅 ) is noninfluential, though this valve is directly associated with the coupling between the RV and PA. This parameter may be more influential if a different valve model, see, e.g., (Mynard et al. , 2012), is used. Parameter Inference and Model Selection
We determine parameter subsets using a combination of correlation analysis and information cri-teria (AICc and BIC). While previous studies utilize sensitivity based subset selection (Olufsen & Ottesen, 2013; Olsen et al. , 2018) or information criteria (Gerringer et al. , 2018; Guan et al. , 2019) for parameter reduction, a combination of the two methods to obtain the most informative, identi-fiable subset is innovative. Of the subsets analyzed in Table 4, subsets 𝜽 %- and 𝜽 M- $ have the lowest AICc and BIC scores. The subset 𝜽 M- $ contains all the parameters identified for 𝒓 % , suggesting that base parameters in 𝜽 %- are the most informative to the model even with additional time-series data. This consistency strengthens the use of this model type in PH assessment, as the number, not the type, of parameters inferred will increase with more data. Including time-series data in the residual allows us to estimate RA and RV timing parameters. However, the time-series data was first shifted to align the model and data, as shown in Fig. 5. Simulations begin during ventricular isovolumic contraction; hence MPA waveforms align with the upstroke of RV pressure. Shifting the RA data ensures that RA systole occurs before RV con-traction. Atrial dynamics consist of three distinct phases, where the atria function as a reservoir, conduit, and contractile chambers (Alenezi et al. , 2018). Our model captured atrial contraction and motivated us to use only the contractile phase of the data during parameter inference. Previous modeling approaches that couple 0D models with higher fidelity models (1D or three dimensional) an predict all three phases (Liang et al. , 2009; Mynard et al. , 2012). We acknowledge and utilize the model’s limitations when selecting portions of the data for parameter inference. This discrep-ancy between the model and true physiological process is an essential source of structural uncer-tainty that can bias parameter estimates (Paun et al. , 2020). Hemodynamic Simulations
As shown in Fig. 6, the most significant biomarker, PVR (represented by 𝑅 ), increases by a factor of four in PH compared to normotension. This agrees with findings reported in the literature for PAH and CTEPH patients (Humbert, 2010; Vonk Noordegraaf et al. , 2017). In addition, cardiac contractility ( 𝐸 * − 𝐸 + ) of the right heart is reduced since both 𝐸 +,-. and 𝐸 +,-) are larger in PH than normotension. Changes in these parameters are associated with both RA and RV function, which is altered in PH; PAH patients have decreased RA longitudinal strain (Sakata et al. , 2016), corresponding to an increased 𝐸 +,-. , which has been linked to RV failure and overload (Tello et al. , 2019). Moreover, increased end-diastolic elastance, 𝐸 +,-) , is negatively correlated with RA a reservoir, passive, and active strain (Tello et al. , 2019), suggesting that RA and RV function dete-riorate simultaneously during PH progression. Our simulations show that pulmonary venous compliance decreases in all but one PH patient. This agrees with results from prior studies that show intimal thickening of the pulmonary veins in CTEPH, decreasing 𝐶 (Gerges et al. , 2020). The additional estimation of 𝑇 .,- , 𝑇 .,, , and 𝑇 ),, does not drastically alter the optimal value of parameters estimated with 𝒓 % . This indicates that addi-tional time-series data do not change the primary conclusion reached using static data. Sensitivity analyses show that 𝑇 ),, is most influential with respect to the residual 𝒓 . This parameter increases slightly in PH patients as shown in Fig. 6b. As shown in Figs. 7 and A1, predictions of RV pressure during systole take up more of the cardiac cycle length in PH (systole is 59 ±
7% of the cardiac cycle in PH vs. 50% in normotension). Physiologically, RV function degrades in PH, leading to an elongated ventricular systole (Driessen et al. , 2018) also seen in our PH simulations. Several authors have used 0D systems-level models to predict hemodynamic function in the pulmonary circulation (Kheyfets et al. , 2016; Colunga et al. , 2020; Tang et al. , 2020). The study by Tang et al. (2020) constructs an expansive 0D circuit model with 25 compartments and simu-lates four PH cases, including distal pulmonary artery stenosis. Their results show model predic-tions are comparable to PH hemodynamics when changing 0D model parameters indicative of the isease but did not estimate patient-specific parameters. Kheyfets et al. (2016) combine an elas-tance RV model with a 3-element Windkessel model of the MPA and verify that their RV-MPA model matches systolic and diastolic MPA pressure across 115 pediatric patients. This model did not consider the interaction between the systemic and pulmonary circulations, accounted for in our study here. To the authors’ knowledge, no investigations have studied improvements in parameter inference when using different RHC data modalities. As shown in Fig. 7, model predictions of systolic, diastolic, and mean pressure are matched to static RHC data, i.e., using 𝒓 % , as done pre-viously (Kheyfets et al. , 2016; Colunga et al. , 2020). Our study is the first to match maximum and minimum RA pressure to model predictions under PH conditions. RA dysfunction is common in PH and is an independent predictor of mortality and hospitalization (Alenezi et al. , 2018). The study by Gerringer et al. (2018) matches dynamic MPA predictions from a Windkessel model to pressure data during PAH progression in rats. This study showed an increased PVR and decreased pulmonary compliance with increasing PAH severity but did not account for systems-level dy-namics or right heart function. In contrast, we estimate RA and RV timing parameters to fit model predictions to dynamic RHC waveforms. Our results provide a -0.75 Pearson correlation coeffi-cient between mean MPA pressure and 𝐶 , supporting the idea that PH severity increases with decreasing compliance. Besides, both 𝐶 and 𝑅 have positive correlations with systolic pressure in the MPA (-0.80 and 0.80, respectively) and RV (-0.79 and 0.80, respectively), again illustrating that pulmonary artery pressure increases due to changes in pulmonary artery compliance and PVR. We fit time-series data using a closed-loop compartmental model, including the right heart and pulmonary arteries, and show that PVR increases in PH. The addition of time series data allows for patient-specific RA and RV predictions which describe the RA-RV coupling. Alenezi et al. (2020) noted that optimal RA-RV coupling is important for efficient right heart performance. Our RV predictions in Fig. 7 fit both systolic and diastolic phases of the data when using 𝒓 ; this former section of data is matched by estimating 𝑇 ),, , whereas improved fits during diastole are due to a calibrated RA model. We conclude that the inclusion of time-series data is crucial for understand-ing the synergistic relationship between the right heart chambers in PH. Our results also show that 𝑅 /𝑅 ! , 𝐸 -),+ /𝐸 , and PI increase with PH, while 𝐶 /𝐶 ! de-creases. The study by Yang et al. (2018) shows that an increased 𝑅 /𝑅 ! correlates with unstable pediatric PAH progression. The pulsatility index, PI , is strongly associated with PH mortality (Mazimba et al. , 2019) and reflects worsening of pulmonary vascular bed stiffness. These two ndices agree with our findings that these ratios increase in PH. We introduce two new indices: the elastance ratio and the compliance ratio. The former increases in PH, reflecting a relative increase in the resting contractile state of the RV, while the latter decreases due to PH-induced pulmonary vascular remodeling (Fayyaz et al. , 2018). We suspect that these indices may differentiate PH types and severity if computed from a larger patient data cohort. PV loops provide useful information regarding heart function, yet these typically cannot be measured in-vivo . Predictions in Fig. 8 show that RV PV loops are similar to those of the left ventricle in PH, whereas normotensive RV simulations have a distinct triangular shape (Tabima et al. , 2017). This is also shown in the computational study by Tang et al. (2020), finding that the PV loop shape is similar between the ventricles in PH. Calculations of CPW (see Table 6) in the left ventricle are similar for normotension and PH, but RV CPW is increased in PH. This agrees with Yang et al. (2018), who use a lumped parameter model to predict CPW in pediatric PAH. Their results show an increase in RV CPW for patients with worsening PH, suggesting a negative cor-relation between RV CPW and disease stability. RA PV loops are not typically used in PH assess-ment, yet these may illustrate how RA function leads to increased RV dysfunction (Alenezi et al. , 2018). PV loops for the RA have an elongated shape in PH and left atrial PV loops have a unique “left-tail” that corresponds to the non-contractile phase of atrial dynamics. The atrial PV loops produced by Tang et al. (2020) are similar in shape to our results, but their model is more complex with 25 compartments. In contrast ours model can predict similar (and patient-specific) dynamics with only six compartments. It is unclear if RA PV loops can expound right heart function and should be investigated further. Pressure predictions in all eight compartments of the model change with PH. Figure 9 (and Fig. A1) shows that pressures in the RA, RV, and pulmonary arteries change significantly in PH; however, systemic arterial and venous pressures are also changed for PH patients. For instance, the elevated systemic arterial pressure in patient 5 requires that both left atrial and ventricular pressures increase, illustrating systems-level changes in cardiovascular function due to PH. RV deterioration ultimately affects the systemic circulation via impaired left ventricle function (Philip et al. , 2019), supporting the use of systems-level modeling for PH. The uncertainty in these predictions, due to uncertainties in the measured data, are also quan-tified. Uncertainty quantification is a crucial step in the model analysis and reveals the volatility of model predictions (Mirams et al. , 2016). Previous studies (Marquis et al. , 2018; Colunga et al. , 020) predicted uncertainty in 0D model predictions using more advanced, Bayesian techniques, requiring numerous simulations. In this study, we utilize a frequentist framework, reducing com-putational complexity. Uncertainty in RA predictions are only available during RA contraction (where we match the model to data), and results in Fig. 9b show that uncertainty in the RA is the largest. Confidence and prediction intervals depend on both model sensitivity and the residual magnitude; hence, the high sensitivity of RA pressure and a more pronounced discrepancy between predictions and time-series data lead to larger uncertainty in our predictions. As detailed in the review by Rosenkranz & Preston (2015), there are several pitfalls in RHC measurements. Meas-urements of PAWP can fluctuate with the respiratory cycle and can be overestimated if catheter balloons are over-inflated. Fig. 9b shows that the relative error in pulmonary venous pressure, matched to PAWP, has a higher relative uncertainty (25% coefficient of variation) than the other pressures; hence, the prediction intervals likely capture the true pulmonary venous pressure even in the presence of measurement noise. Moreover, Rosenkranz & Preston (2015) conclude that RHC measurement errors increase the probability of misdiagnosis. This again suggests that accounting for measurement error in model predictions is crucial in identifying PH severity. Treatment Interventions
Patient-specific PH treatments include three subgroups: 1) vasodilators, 2) surgical intervention, and 3) surgical intervention with vasodilators; note that all patients are treated in subgroup one, whereas subgroups two and three are simulated in patients 3-5 as those treatments are CTEPH-specific. Results show that the treatments associated with surgical intervention reduce mPAP in patients 3 and 4 below the PH threshold. In general, surgical intervention with and without drug therapy yields a larger relative change in mPAP than drug therapy alone; 41.9 ± ± et al. , 2020). The average change in mPAP in PAH patients, one and two, is 28.2 ± et al. , 2019; Hien et al. , 2020). These results are consistent with the notion that PAH is not curable, but disease progression is manageable through therapeutic drugs. ystemic hypertension is defined by the European Society of Cardiology as a systolic/diastolic systemic arterial blood pressure above 140/90 mmHg, and hypotension is defined as pressures below 90/60 mmHg (Sharma et al. , 2020). Results in Fig. 10, show that patients 1 and 4 have a baseline systemic systolic arterial blood pressure above 140 mmHg while no patients exhibit base-line diastolic pressures near the hypotensive regime. Systemic artery pressures decrease in all pa-tients for treatments T1-T4 and T8-T9, consistent with pharmaceutical knowledge that vasodilators reduce blood pressure. Prior studies show that PH vasodilators decrease systemic resistance (Shanmugam et al. , 2015), and adverse effects include potential systemic hypotension. Our model depicts system-wide changes that do affect the systemic arteries, and treatment T1 shows a de-crease is diastolic systemic arterial blood pressure below the 60 mmHg threshold for patient 4, while patient 2 is right at the cutoff. Alternatively, the treatments associated with surgical inter-vention alone show an increased systolic and diastolic systemic arterial blood pressure with no indication of hypotension. This may be associated with immediate physiological effects from re-moving obstructions in the PA, leading to an immediate increase in blood flow through the left heart and systemic arteries. Limitations
There are several shortcomings of this study, in both the data acquisition and the modeling. We do not have simultaneous data in the RA, RV, and PA compartments. To remedy this, we aligned data to ECG and shifted data to align with the model as part of the parameter inference. Alternatively, we could have adjusted initial conditions in the model. This strategy is essential, as it is not possible with RHC to measure PH pressures simultaneously. RHC measurement of PAWP is used as a surrogate for pulmonary venous pressure, yet PAWP may overestimate the true pressure in these vessels. Moreover, in the absence of left-heart dynamic data, we cannot infer parameters specific to the left atrium and ventricle timing. Adding these measurements are essential when studying group 2 PH due to left-sided heart failure and should be considered in future studies. We selected a single RHC waveform for each compartment, yet data included multiple heartbeats of data. Fit-ting the model to multiple heartbeats may reveal which parameters adapt with the respiratory cycle. Our model describes valves as simple diodes, yet valve dynamics are complex in both normo-tensive and PH conditions. Also, our atrial model does not predict the distinctive “v loop” in sim-ulated PV loops (Pironet et al. , 2013). Prior modeling studies incorporated more sophisticated alve dynamics (Mynard et al. , 2012) circumvented both of these issues and may provide better fits to dynamic RA data. Lastly, we did not consider the effects of the ventricular septum, which play a role in normotensive and PH RV dynamics. While ventricular-ventricular interaction (VVI) is always present when a patient develops PH, the interaction between the ventricles has a more significant influence on cardiovascular function through the thickening of the heart muscles and the deflection of the septal wall towards the left ventricle. Accounting for VVI in this model type leads to a system of differential algebraic equations, requiring simultaneous solutions of an implicit nonlinear algebraic equation and a system of differential equations. We suspect that including VVI in the model will allow us to predict indices indicative of RV dysfunction and will pursue this in future studies.
CONCLUSION
This study uses a 0D, systems-level hemodynamics model to predict changes in cardiovascular parameters due to PH. We utilize sensitivity analyses, subset selection techniques, and information criteria to deduce the best parameter subsets for two residuals: one with static data and one with static data and dynamic RA, RV, and PA pressure tracings. Our results show that time-series data allows for timing parameters in the right heart to be estimated. Overall, model outcomes are con-sistent with the physiological understanding of the disease. PH increases PVR ( 𝑅 ), the RV elas-tance ratio, 𝑃𝐼 , and RV 𝐶𝑃𝑊 and decreases pulmonary arterial compliance and the value of 𝐶 /𝐶 !. . By simulating treatment strategies, we predict an effective reduction in mPAP pressure when CTEPH patients undergo surgical treatment and that PAH patients reduce their mPAP to a “mild” PH level. ADDITIONAL INFORMAITON Computer Code
Computer code with waveforms and static data for each subject as well as for treatment simula-tions is available on the GitHub repository https://github.com/mjcolebank/CDG_NCSU in the “Cardiovascular-Systems-Model-PH” folder.
Competing Interests
None of the authors has any conflicts of interests on the submission.
Funding
This work is funded by the National Science Foundation, Division of Mathematical Sciences Award
Acknowledgements
We thank Drs. Geeshath Jayasekera and Martin Johnson from the Scottish Pulmonary Vascular Unit at Golden Jubilee Hospital in Glasgow, UK, for providing RHC data from PH patients 1,2 and 3. We also thank Dr. Sudarshan Rajagopal at Duke University Medical Center for providing RHC data for patients 4 and 5.
Author Contributions
All persons designated as authors qualify for authorship, and all those who qualify for authorship are listed.
Mitchel J Colebank:
Conception or design of the work; Acquisition, analysis or interpretation of data for the work; Drafting the work or revising it critically for important intellectual content; Approval of the final version of the manuscript; Agreement to be accountable for all aspects of the work.
Amanda L Colunga:
Conception or design of the work; Acquisition, analysis or interpretation of data for the work; Drafting the work or revising it critically for important intellectual content; Approval of the final version of the manuscript; Agreement to be accountable for all aspects of the work.
Macie King:
Acquisition, analysis or interpretation of data for the work; Drafting the work or revising it critically for important intellectual content; Approval of the final version of the manu-script; Agreement to be accountable for all aspects of the work. hristopher Schell:
Acquisition, analysis or interpretation of data for the work; Drafting the work or revising it critically for important intellectual content; Approval of the final version of the man-uscript; Agreement to be accountable for all aspects of the work.
Matt Shelden:
Acquisition, analysis or interpretation of data for the work; Drafting the work or revising it critically for important intellectual content; Approval of the final version of the manu-script; Agreement to be accountable for all aspects of the work.
Mariam Kharbat:
Acquisition, analysis or interpretation of data for the work; Drafting the work or revising it critically for important intellectual content; Approval of the final version of the man-uscript; Agreement to be accountable for all aspects of the work.
Robert Sternquist:
Acquisition, analysis or interpretation of data for the work; Drafting the work or revising it critically for important intellectual content; Approval of the final version of the manuscript; Agreement to be accountable for all aspects of the work.
Mette Olufsen:
Conception or design of the work; Acquisition, analysis or interpretation of data for the work; Drafting the work or revising it critically for important intellectual content; Approval of the final version of the manuscript; Agreement to be accountable for all aspects of the work.
APPENDICIES Appendix 1: Local sensitivity analysis
Using the residual vectors 𝒓 𝒋 , the log-scaled sensitivity with respect to 𝜃 is computed by 𝜕𝒓 𝒋 (𝑡, 𝒙; 𝜽)𝜕 log(𝜃 ) = 𝜕𝒓 𝒋 (𝑡, 𝒙; 𝜽)𝜕𝜃 𝜃 , 𝑖 = 1, 2, … , ℳ, 𝑗 = 1, 2 (A1) where 𝑡 (s) denotes time, 𝒙(𝑡; 𝜽) the dependent state variables, 𝜽 the model parameters, and ℳ the number of parameters. Log-scaling the sensitivities puts parameters of different magnitude on a similar scale (Marquis et al. , 2018; Colebank et al. , 2019; Colunga et al. , 2020). We approximate the sensitivity equations via a forward finite difference, where the sensitivity of the residual with respect to 𝜃 is the column vector 𝝌 (𝑡) . Lastly, we use the 2-norm of the sensitivity results, i.e., (cid:138)|𝝌 |(cid:138) , 𝑖 = 1,2, … , ℳ, to rank the parameters from most to least influential. ppendix 2: Global Sensitivity Analysis GSA quantifies model sensitivity by sampling over the plausible parameter space. While GSA methods are more computationally expensive than local methods, they can expose undiscovered relationships between parameters by varying multiple parameters at a time (Eck et al. , 2016).
In this study, we use variance-based Sobol’ indices (Sobol, 2001; Saltelli et al. , 2010) to quan-tify parameters’ influence on the variance of the residual vectors. We begin by constructing phys-iological parameter regimes, enforcing a parameter space of ±30% from the nominal parameter estimates; these regimes were first analyzed to ensure that model outputs are physiological. Consider the quantity of interest 𝒓 𝒋 (𝑡, 𝒙; 𝜽) as before, with parameters 𝜽 = [𝜃 % ,··· , 𝜃 ℳ ] (note that we drop the dependence of 𝒓 on 𝒙 and 𝑡 for clarity hereon). Each parameter 𝑖 lies within the physiologically admissible parameter space 𝛤 , i.e., the entire parameter domain is ⋃ Γ ℳ = Ω ℳ ⊂ 𝐑 ℳ . The expectation and variance of 𝒓 𝒋 , 𝑗 = 1, 2, are defined as 𝐸 “𝒓 𝒋 (𝜽)« = ‹ 𝒓 𝒋 (𝜽) 𝑑𝜽 P ℳ , (A2) 𝑉 “𝒓 ; (𝜽)« = ‹ ›𝒓 𝒋 (𝜽) − 𝐸 “𝒓 𝒋 (𝜽)«fi 𝑑𝜽 P ℳ = 𝐸P𝒓 𝒋 (𝜽) Q − 𝐸 “𝒓 𝒋 (𝜽)« . (A3) We are interested in computing the conditional expectation and variance of residuals when a single parameter 𝜃 is known; hence, we can define the operators 𝐸 N~ “𝒓 𝒋 (𝜽|𝜃 )« = ‹ 𝒓 𝒋 (𝜽|𝜃 )𝑑𝜃 P ℳ& , Ω
ℳ$% = Ω ℳ \Γ 𝑗 = 1, 2, (A4) which does not include 𝜃 , and the partial variance 𝑉 N ! = 𝑉 “𝐸 N~ P𝒓 𝒋 |𝜃 Q« , 𝑗 = 1, 2. (A5) Equation (A4) measures the variance of the expected value of 𝒓 ; , conditioned on the fixed, known parameter 𝜃 ; i.e., it measures the variance of the mean not attributed to 𝜃 . From this we define the first-order sensitivity measure 𝑆 and the total order effects 𝑆 / ! 𝑆 = 𝑉 N ! 𝑉P𝒓 𝒋 Q, 𝑆 / ! = 𝐸 N~ “𝑉 N ! P𝒓 𝒋 (cid:138)𝜃 ~ Q«𝑉P𝒓 𝒋 Q , 𝑗 = 1 ,2 . (A6) The former measures the relative contribution of 𝜃 to the output variance, while 𝑆 / ! measures all higher-order interactions with 𝜃 that contribute to the variance. Note that “𝐸 N ! P𝒓 𝒋 |𝜃 ~ Q«𝑉P𝒓 𝒋 Q = 1 ⇒ 𝑉 N ~! “𝐸 N ! P𝒓 𝒋 |𝜃 ~ Q« = 𝑉P𝒓 𝒋 Q, 𝑗 = 1, 2 (A7) hence 𝑆 / ! ≈ 0 implies that higher-order effects are negligible. To compute the Sobol’ indices, we use the Saltelli algorithm (Saltelli et al. , 2006) as described in Algorithm 1. Algorithm 1:
Sobol’ Indices (1) Generate two ( 𝓃 × ℳ) sample matrices 𝐴 and 𝐵 , with pseudorandom entries 𝜃 and 𝜃(cid:149) drawn from respective densities 𝐴 = ‡ 𝜃 %% ⋯ 𝜃 ℳ% ⋮ ⋱ ⋮𝜃 %𝓃 ⋯ 𝜃 ℳ𝓃 • 𝐵 = ‡ 𝜃(cid:149) %% ⋯ 𝜃(cid:149) ℳ% ⋮ ⋱ ⋮𝜃(cid:149) %𝓃 ⋯ 𝜃(cid:149) ℳ𝓃 • (2) Generate ℳ matrices 𝐴 S which is equal to the matrix 𝐴 except the 𝑖 column is the 𝑖 column from 𝐵 . Similarly, create 𝐵 U . (3) To estimate the total variance, generate a matrix appending 𝐵 to 𝐴 𝐶 = ‚𝐴𝐵„ (4) Evaluate the model at each row of the matrices 𝐴 and 𝐵, with outputs 𝑓(𝐴) ; and 𝑓(𝐵) ; for 𝑗 = 1, ⋯ , 𝓃 . This requires model evaluations (5) Evaluate the model at each row of the matrices 𝐴 S and 𝐵 U , with outputs 𝑓P𝐴 S Q ; and 𝑓P𝐵 U Q ; for 𝑗 = 1, ⋯ , ℳ . This requires model evaluations (6) Estimate the first-order Sobol’ indices approximated by 𝑆 ≈ 1𝓃 ∑ »𝑓(𝐴) ; 𝑓P𝐵 U Q ; − 𝑓(𝐴) ; 𝑓(𝐵) ; … 𝓃;V%
12𝓃 ∑ 𝑓(𝐶) ; 𝑓(𝐶) ; − 𝐸 [𝑓(𝐶)] (7) Estimate the total Sobol’ indices approximated by 𝑆 / ! ≈ 12𝓃 ∑ »𝑓(𝐴) ; − 𝑓P𝐴 S Q ; … 𝓃;V%
12𝓃 ∑ 𝑓(𝐶) ; 𝑓(𝐶) ; − 𝐸 [𝑓(𝐶)] TABLE A1-A2 FIGURE A1 Appendix 3: Uncertainty Quantification
Confidence and prediction intervals are computed following methods from nonlinear regression (Seber & Wild, 2003). In this study, we compute the relative residual by scaling the difference between the measured data and model prediction by the data itself. The weighted sum of squares error (WSSE) is then
WSSE = P𝒚 − 𝒚 𝒅 Q I 𝑾P𝒚 − 𝒚 𝒅 Q, 𝑾 = diag “P𝑦 %7 Q $6 , P𝑦 Q $6 , … « (A8) where 𝑾 is the weight matrix with the squared inverse of the data along the diagonal. Applying an asymptotic analysis and using the model sensitivity 𝜒 ≈ 𝜕𝑔/𝜕𝜃 gives the parameter estimator 𝛉(cid:146) XYYZ = 𝜽 + (𝝌 I 𝑾𝝌) $% 𝝌 I 𝑾𝝐 (A9) where 𝜽 are the true, unknown parameters and 𝝐 is a normally distributed, iid, random variable. Note that E~𝜽(cid:146)
XYYZ (cid:127) = 𝜽 and the covariance
Var~𝛉(cid:146)
XYYZ (cid:127) = 𝑽
XYYZ is 𝑽 XYYZ = 𝜎 <6 (𝝌 I 𝑾𝝌) $% 𝝌 I 𝑾𝑾𝝌(𝝌 I 𝑾𝝌) $𝟏 (A10) where 𝜎 < is the error variance, and its estimator 𝜎(cid:144) <6 can be derived using the weight matrix and model sensitivity (see Meermeyer (2015) for details). This definition of 𝑉 XYYZ is used in the con-fidence and prediction interval construction in equations (27) - (30). eferences
Alenezi F, Mandawat A, Il’Giovine ZJ, Shaw LK, Siddiqui I, Tapson VF, Arges K, Rivera D, Romano MMD, Velazquez EJ, Douglas PS, Samad Z & Rajagopal S (2018). Clinical Utility and Prognostic Value of Right Atrial Function in Pulmonary Hypertension.
Circ Cardiovasc Imaging e006984. Alenezi F, Rajagopal S & Kutty S (2020). Assessing right atrial function in pulmonary hypertension: Window to the soul of the right heart?
Am J Physiol Heart Circ Physiol
H154–H155. Beneken JEW & DeWit B (1966). A physical approach to hemodynamic aspects of the human cardiovascular system. In
Physical Bases of Circulatory Transport: Regulation and Exchange. , ed. Guyton A & Reeve E, pp. 1–45. Saunders, Philadelphia. Boron WF & Boulpaep EL (2017).
Medical physiology . Philadelphia, PA : Elsevier. Chambers MJ, Colebank MJ, Qureshi MU, Clipp R & Olufsen MS (2020). Structural and hemodynamic properties of murine pulmonary arterial networks under hypoxia-induced pulmonary hypertension.
Proc Inst Mech Eng Part H J Eng Med , J Inverse Ill-Posed Probl
Int J Numer Method Biomed Eng , e3242. Colunga AL, Kim KG, Woodall NP, Dardas TF, Gennari JH, Olufsen MS & Carlson BE (2020). Deep phenotyping of cardiac function in heart transplant patients using cardiovascular system models. J Physiol
N Engl J Med
PLoS One
Int J Numer Method Biomed Eng
Cardiovasc Eng Circulation
Arch Pathol Lab Med
Circulation , Physiol Rep a ). Riociguat for the treatment of chronic thromboembolic pulmonary hypertension. N Engl J Med b ). Riociguat for the treatment of pulmonary arterial hypertension. N Engl J Med
Biomech Model Mechanobiol , 1213–1232. Hidalgo JU, Nadler SB & Bloch T (1962). The use of the electronic digital computer to determine best fit of blood volume formulas. J Nucl Med Prog Pediatr Cardiol doi:10.1016/j.ppedcard.2020.101230. Hoeper MM, Ghofrani H-A, Grünig E, Klose H, Olschewski H & Rosenkranz S (2017). Pulmonary hypertension.
Dtsch Arztebl Int
Eur Respir Rev
Sci Rep
Iterative Methods for Optimization . Frontiers in Appplied Mathematics. Society for Industrial and Applied Mathematics. Philadelpha, PA. Kheyfets VO, Dunning J, Truong U, Ivy D, Hunter K & Shandas R (2016). A zero-dimensional model and protocol for simulating patient-specific pulmonary hemodynamics from limited clinical data.
J Biomech Eng
J Biomech Eng
Am J Cardiol
J Biomech Eng
Med Biol Eng Comput
JACC Cardiovasc Imaging Circulation
Math Biosci
Heart Lung Circ
Comput Stat
J Physiol
Orphanet J Rare Dis Int J Numer Method Biomed Eng
Biol Cybern , 121–138. Olufsen MS & Ottesen JT (2013). A practical approach to parameter estimation applied to model predicting heart rate regulation. J Math Biol
J R Soc Interface
Am J Physiol Heart Circ Physiol, , H1167-H1177. Pironet A, Dauby PC, Paeme S, Kosta S, Chase JG & Desaive T (2013). Simulation of left atrial function using a multi-scale model of the cardiovascular system.
PLoS One , e65146. Pope S, Ellwein LM, Zapata C, Novak V, Kelley C & Olufsen MS (2008). Estimation and iden-tification of parameters in a lumped cerebrovascular model. Math Biosci Eng Biomech Model Mechanobiol
Int J Cardiovasc Imaging
Physiol Rep Physiol Rep Eur Respir Rev
J Echocardiogr
Comput Phys Commun
Reliab Eng Syst Saf
Int J Numer Method Biomed Eng
Nonlinear Regression . Wiley-Interscience, Hoboken, NJ, USA. Shanmugam E, Jena A & George M (2015). Riociguat: Something new in pulmonary hypertension therapeutics?
J Pharmacol Pharmacother StatPearls . Siennicka A et al. (2019). Treatment of chronic thromboembolic pulmonary hypertension in a multidisciplinary team.
Ther Adv Respir Dis , 1753466619891529. Simonneau G, Montani D, Celermajer DS, Denton CP, Gatzoulis MA, Krowka M, Williams PG & Souza R (2019). Haemodynamic definitions and updated clinical classification of pulmonary hypertension. Eur Respir J , 1801913. Sobol IM (2001). Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Math Comput Simul
J R Soc Interface Physiology
J Cardiovasc Transl Res
13, 826–852. Tello K, Dalmer A, Vanderpool R, Ghofrani HA, Naeije R, Roller F, Seeger W, Wiegand M, Gall H & Richter MJ (2019). Right ventricular function correlates of right atrial strain in pulmonary hypertension: A combined cardiac magnetic resonance and conductance catheter study.
Am J Physiol Heart Circ Physiol
J Am Coll Cardiol
Pulm Circ PLoS One , e78569. Williams ND, Brady R, Gilmore S, Gremaud P, Tran HT, Ottesen JT, Mehlsen J & Olufsen MS 2019). Cardiovascular dynamics during head-up tilt assessed via pulsatile and non-pulsatile models. J Math Biol
Pulm Circ TABLES
Table 1.
Patient demographics for this study; group 1: pulmonary arterial hypertension (PAH); group 4: chronic thromboembolic pulmonary hypertension (CTEPH).
Patient PH group Age Sex Height (cm) Weight (kg) CO (L/min)
1 1 64 Male 164.0 72.6 4.0 2 1 58 Male 161.0 70.0 4.3 3 4 27 Female 151.0 81.1 2.6 4 4 71 Female 167.6 93.3 6.1 5 4 51 Male 179.1 117.2 3.6 CO: cardiac output; PH: pulmonary hypertension able 2.
Static values obtained from patient data are used in nominal parameter calculation. Mean and standard deviation values are calculated for PH data only. † Normotensive control values obtained from (Boron & Boulpaep, 2017). ‡ left atrial diastolic value used in place of PAWP.
Data Control Patient 1 Patient 2 Patient 3 Patient 4 Patient 5 Mean± SD 𝑝 "
15 10 28 9 16 33 19 ± 11 𝑝 "
3 6 15 2 8 25 11 ± 9 𝑝 "
26 87 91 93 69 81 84±10 𝑝 "
2 3 5 3 1 18 6±7 𝑝 ")$,&
17 85 90 92 68 80 83±10 𝑝 ")$,’
9 32 38 34 23 30 32±5 𝑝̅ ")$
12 48 55 54 43 50 50±5 𝑝 "* 𝑝 "+$,&
119 158 112 127 148 118 133±20 𝑝 "+$,’
76 85 76 90 78 77 81±6 𝑝̅ "+$
88 109 88 102 101 91 98±9 able 3.
Parameter changes used to simulate treatment strategies in PAH and CTEPH.
PAH & CTEPH Parameter 𝑅 ) 𝑅 + 𝐶 +$ 𝐶 )$ Drug therapy Treatment 1 -40% -30% +30% +40% Treatment 2 -40% -20% +20% +40% Treatment 3 -40% -10% +10% +40% Treatment 4 -20% -10% +10% +20%
CTEPH Parameter 𝑅 ) 𝑅 + 𝐶 +$ 𝐶 )$ Surgical intervention Treatment 5 -80% - - +80% Treatment 6 -60% - - +60% Treatment 7 -40% - - +40% Surgical intervention & drug therapy Treatment 8 -85% -20% +20% +85% Treatment 9 -65% -20% +20% +65% Treatment 10 -45% -20% +20% +45% able 4.
Parameters included in residual vectors 𝒓 , and 𝒓 and the AICc and BIC scores for each param-eter subset. Columns in red are selected as the best parameter subsets. 𝒓 𝒓 Parameter 𝜽 , ! 𝜽 / ! 𝜽 ! 𝜽 ! 𝜽 ! 𝜽 , " 𝜽 / " 𝜽 " 𝜽 " 𝜽 " 𝑅 + X X X X X X X X X X 𝑅 ) X X X X X X X X X X 𝑅 X X X X X X X X X X 𝐶 +$ X X X X X X X X X X 𝐶 )$ X X X X X X X X X X 𝐶 )( X X X X X X X X X X 𝐸 X X X X X X X X X X 𝐸 X X X 𝐸 X X X X X X X X 𝐸 X X X X X X X 𝑇 $, X X 𝜏 $,5 X X X X X 𝑇 $,5 X X X 𝑇 (,5 X X X X X X X X AICc 18.5 18.8 18.7 18.7 18.8 23.8 23.8 25.8 23.5 23.8 BIC 90.6 90.9 90.9 90.9 91.0 200.5 200.5 217.2 158.3 158.7 able 5.
Estimated parameter values using 𝒓 / . Results are presented as the optimal value ± 𝜽 Patient 1 Patient 2 Patient 3 Patient 4 Patient 5 𝑅 + 𝑅 ) 𝑅 𝐶 +$ 𝐶 )$ 𝐶 )( 𝐸 𝐸 𝐸 𝑇 $, 𝑇 $,5 𝑇 (,5 able 6. Model outcomes from normotensive and PH simulations.
Patient LA CPW LV CPW RA CPW RV CPW 𝑹 𝒑 /𝑹 𝒔 𝑪 𝒑𝒂 /𝑪 𝒔𝒂 𝑬 𝒓𝒗,𝒎 /𝑬 𝒍𝒗,𝒎 𝐏𝐈 Normotensive
PAH
PAH
CTEPH
CTEPH
CTEPH
Indices include cardiac power (CPW, Watts) in all four heart chambers, resistance ratios (dimensionless), compliance ratios (dimensionless), ventricular elastance ratios (dimensionless), and pulsatility index (PI, dimensionless) calculated after estimating parameters using 𝑟 / . LA – left atrium, LV – left ventricle, RA – right atrium, RV – right ventricle. able A1.
Model parameters and how they are calculated.
Parameter Units Equation Reference Heart Valves 𝑅 $($ mmHg smL 𝑝 − 𝑝 +$,& 𝑞 Ohm’s Law 𝑅 ’($ mmHg smL 𝑝 )( − 𝑝 𝑞 Ohm’s Law 𝑅 )($ mmHg smL 𝑝 − 𝑝 )$,& 𝑞 Ohm’s Law 𝑅 mmHg smL 𝑅 +( mmHg smL 𝑝̅ +( − 𝑝 𝑞 Ohm’s Law 𝑅 )( mmHg smL 𝑝̅ )( − 𝑝 𝑞 Ohm’s Law
Systemic Vasculature 𝑅 + mmHg smL 𝑝̅ +$ − 𝑝̅ +( 𝑞 Ohm’s Law 𝐶 +$ mLmmHg 𝑉 +$,& 𝑝 +$,& (Marquis et al. , 2018) 𝐶 +( mLmmHg 𝑉 +(,& 𝑝 +( (Marquis et al. , 2018) Pulmonary Vasculature 𝑅 ) mmHg smL 𝑝̅ )$ − 𝑝̅ )( 𝑞 Ohm’s Law 𝐶 )$ mLmmHg 𝑉 )$,& 𝑝 )$,& (Marquis et al. , 2018) 𝐶 )( mLmmHg 𝑉 )(,& 𝑝 )(,& (Marquis et al. , 2018) Heart Elastance 𝐸 mmHgmL 𝑝 𝑉 (Marquis et al. , 2018) 𝐸 mmHgmL 𝑝 𝑉 (Marquis et al. , 2018) 𝐸 mmHgmL 𝑝 𝑉 (Marquis et al. , 2018) mmHgmL 𝑝 𝑉 (Marquis et al. , 2018) 𝐸 mmHgmL 𝑝 𝑉 (Marquis et al. , 2018) 𝐸 mmHgmL 𝑝 𝑉 (Marquis et al. , 2018) 𝐸 mmHgmL 𝑝 𝑉 (Marquis et al. , 2018) 𝐸 mmHgmL 𝑝 , 𝑚𝑉 , 𝑀 (Marquis et al. , 2018) Heart Timing 𝜏 𝑠 data - 𝑇 𝑠 data - 𝑇 𝑠 data - 𝑇 𝑠 data - 𝑇 𝑠 data - 𝜏 𝑠 (Boron & Boulpaep, 2017) 𝑇 𝑠 (Boron & Boulpaep, 2017) 𝑇 𝑠 𝑇 (Boron & Boulpaep, 2017) 𝑇 𝑠 (Boron & Boulpaep, 2017) 𝑇 𝑠 𝑇 (Boron & Boulpaep, 2017) able A2. Optimized Parameters using 𝒓 , . 𝜽 Patient 1 Patient 2 Patient 3 Patient 4 Patient 5 𝑅 + 𝑅 ) 𝑅 𝐶 +$ 𝐶 )$ 𝐶 )( 𝐸 𝐸 𝐸 Figure 1. Data processing
Dynamic data from the right atrium (RA), right ventricle (RV), and main pulmonary artery (MPA) for each patient are digitized from right heart catheterization recordings and used for model calibration. Circles denote the maximum and minimum values used for static data.
Figure 2. Model schematic
Schematic of the compartmental model following an electrical RC-circuit analog. There are eight com-partments, including systemic and pulmonary arteries and veins, and four heart chambers with diode-type valves. Each compartment is connected via a resistance. The right atrium, right ventricle, and pulmonary arteries (boxes in red) use both dynamic and static data for parameter inference. The pulmonary vein and systemic arteries also contain static data used for parameter inference. RHC: right heart catheterization; CO: cardiac output.
Dynamic RHC data- Static data - CO data
Figure 3. Heart chamber elastance function
Representative elastance function for the atrial (red) and ventricular (blue) heart chamber. Timing param-eters are shown above their respective phases of the cardiac cycle. Note that ventricular isovolumic con-traction occurs while the atrium is still relaxing.
Figure 4. Local and global sensitivity analysis results
The two-norm of the local sensitivities using 𝒓 , (a) and 𝒓 / (b) are compared for each PH patient. The sen-sitivity values are scaled by the maximum parameter influence for each patient. First-order ( 𝑆 ? ) and total order ( 𝑆 @ ) Sobol’ indices are computed over the plausible parameter space (c). Parameter ranking based on these indices are shown (d), with the parameter order based on 𝑆 @ magnitude. R p v E M , l a E m , l a SS T SS T a , c (cid:111) R s R p R a v a R m v a R p v a R t v a R sv C s a C sv C pa C p v E M , r a E m , r a E M , r v E m , r v E M , l v E m , l v T a , r T a , c T v , c T v , r -4 -2 S T S T (c)(d) R p v E M , l a E m , l aa , c (cid:111) R s R p R a v a R m v a R p v a R t v a R sv C s a C sv C pa C p v E M , r a E m , r a E M , r v E m , r v E M , l v E m , l v T a , r T a , c T v , c T v , r -3 -2 -1 Lo c a l s en s i t i v i t i e s : -4 -2 R p v E M , l a E m , l aa , c (cid:111) R s R p R a v a R m v a R p v a R t v a R sv C s a C sv C pa C p v E M , r a E m , r a E M , r v E m , r v E M , l v E m , l v T a , r T a , c T v , c T v , r R p v E M , l a E m , l aa , c (cid:111) R s R p R a v a R m v a R p v a R t v a R sv C s a C sv C pa C p v E M , r a E m , r a E M , r v E m , r v E M , l v E m , l v T a , r T a , c T v , c T v , r r (a)(b) Parameters
Pat 1Pat 2Pat 3Pat 4Pat 5 Lo c a l s en s i t i v i t i e s : r (cid:54) (cid:82)(cid:69)(cid:82) (cid:79)(cid:183) (cid:3) (cid:76) (cid:81)(cid:71) (cid:76) (cid:70) (cid:72) (cid:86) G l oba l s en s i t i v t y r an k i ng Figure 5. Data alignment
Representative data (from Patient 5) shifted forward and backward in time to best align with model pre-dictions. (a) Least-squares cost as a function of the shift in the time series data. This minimizer of the least-squares is shown in red. (b) Original data (black) along with the shifted data (red) that minimizes the least-squares cost using the model predictions (blue). All possible shifts are also plotted (grey). -0.15 -0.05 0.05 0.15Shift00.40.81.2
Lea s t s qua r e s c o s t RV -0.05 0.05 0.15Shift1.21.622.42.8 10 -4 RA -2 (a)(b) RV RAPA P r e ss u r e ( mm H g ) Time (s)0 T 0 T 0 TT/2T/2 T/2
Figure 6. Changes in parameters due to PH
Box and whisker plots with quantiles and outliers for the estimated parameters. Each PH patients’ param-eters are estimated using static data only ( 𝒓 , , (a)) or a combination of static and dynamic data ( 𝒓 / , (b)). Results are shown as the relative difference from the normotensive parameters. The lower and upper lim-its of the box represent the 25 th and 75 th quantile, and the red line indicates the 50 th quantile (i.e., the me-dian). -2-1012345-2-1012345 (a)(b) r r R e l a t i v e d i ff e r en c e be t w een P H and no r m o t en s i v e pa r a m e t e r s R s R p R t v a C s a C pa C p v E m , r a E m , r v E m , l v T a , r T a , c T v , c R s R p R t v a C s a C pa C p v E m , r a E m , r v E m , l v Figure 7. Comparison of PH data and model predictions
Predictions of both cardiac output and pressure in all eight compartments using the optimal parameter es-timates from Patient 5. Predictions from the model using 𝒓 , (dashed, red line) and 𝒓 / (solid, red line) are contrasted to the data (solid, black line). Note that time-series data is better matched when using 𝒓 / , yet some static predictions are not as accurate compared to predictions with 𝒓 , . C O ( m L / s e c ) ppa ( mm H g ) p sv ( mm H g ) p l a ( mm H g ) p l v ( mm H g ) p r a ( mm H g ) pp v ( mm H g ) p s a ( mm H g ) p r v ( mm H g ) _ Data r r
22 20.2 20.4 20.6 20.2 20.4 20.620.2 20.4 20.620.2 20.4 20.620.2 20.4 20.620.2 20.4 20.6
Figure 8. Predicted pressure-volume loops in PH and normotensive conditions
Post-optimization model predictions of pressure and volume are integrated to provide pressure-volume (PV) loops. Results shown are from using 𝒓 / in the objective function. Note that simulated PV loops in both the left and right atrium (LA and RA, respectively) do not exhibit the typical “figure-8” dynamic seen in-vivo . Left and right ventricular volume loops (LV and RV) are also shown. Note that the shape of the RV PV loop in normotensive conditions has a more rounded shape, whereas the RV PV loop in PH is similar in shape to the LV PV loop. P r e ss u r e ( mm H g ) Volume (mL) P r e ss u r e ( mm H g ) Volume (mL)Normotensive Patient 1 Patient 2 Patient 3 Patient 4 Patient 50
LA RALV RV
Figure 9. Comparison of normotensive and PH hemodynamics and predicted uncertainty
Parameters are used to generate a normotensive control subject (a) contrasted to a representative PH pa-tient (patient 5) (b). Using the estimated parameters and available data, 95% confidence and prediction intervals can be produced for both systolic and diastolic values and for dynamic pressure in the right atrium (RA), right ventricle (RV), and pulmonary arteries (PA) (b). Note that the uncertainty in RA pres-sure predictions can only be calculated for atrial systole, since that is the only component used in the re-sidual. LA: left atrium; LV: left ventricle; SA: systemic arteries; SV: systemic veins; PV: pulmonary veins; CO: cardiac output; CI: confidence interval; PI: prediction interval.
Figure 10. Simulated treatment strategies
Treatment (T) by either vasodilators (a) or a combination of surgical intervention and vasodilators (b) are simulated by reducing resistance with and without an increase in compliance (see Table 3 for details). Mean pulmonary artery pressure (mPAP) and systolic and diastolic systemic artery (SA, unshaded and shaded bars, respectively) predictions. Note that panel (a) is applied to all five patients, whereas surgical intervention in panel (b) only applies to CTEPH patients. Changes in resistance and compliance in the vasodilator only scenarios do not decrease pressure below 20 mmHg, whereas treatments 5 and 8 do re-duce mean PA pressure below 20 mmHg, suggesting that these patients are “cured” after treatment. Cut-offs for PH (20 mmHg), systolic systemic hypertension (140 mmHg), and diastolic systemic hypotension (60 mmHg) are shown with dashed lines.
T10T10
PH cutoff
Baseline T1 T2 T3 T4204060 m PAP ( mm H g ) Baseline T1 T2 T3 T4
Treatment SA P r e ss u r e ( mm H g ) Baseline T5 T6 T7 T8 T9Baseline T5 T6 T7 T8 T9
Treatment (a) Vasodilator only (PAH & CTEPH) (b) Surgical intervention and vasodilator (CTEPH)
Patient 1 Patient 2 Patient 3 Patient 4 Patient 5 SA systolichypertensionSA diastolichypotension
Figure A1. Model predictions for all five patients
Normotensive predictions (a) are contrasted to optimal PH predictions in patients 1-5 (panels (b)-(f), re-spectively). PH data is shown along with optimal model predictions using either 𝒓 , or 𝒓 / .
29 29.5 3005 29 29.5 300100200 29 29.5 3010015029 29.5 30678 29 29.5 306810 29 29.5 3005010029 29.5 30406080 29 29.5 3005 788029 29.5 30 C O ( m L / s e c ) p ( mm H g ) pa p ( mm H g ) sv p ( mm H g ) l a p ( mm H g ) l v p ( mm H g ) r a p ( mm H g ) s a p ( mm H g ) r v p ( mm H g ) p v _ C O ( m L / s e c ) p ( mm H g ) pa p ( mm H g ) sv p ( mm H g ) l a p ( mm H g ) l v p ( mm H g ) r a p ( mm H g ) s a p ( mm H g ) r v p ( mm H g ) p v _
50 5151015 50 51050100 50 519011013050 511.52 50 510510 50 5105010050 5150100 50 517.588.5 47 48 49607080 C O ( m L / s e c ) p ( mm H g ) pa p ( mm H g ) sv p ( mm H g ) l a p ( mm H g ) l v p ( mm H g ) r a p ( mm H g ) s a p ( mm H g ) r v p ( mm H g ) p v _ Time (s)
Time (s)
Time (s) C O ( m L / s e c ) p ( mm H g ) pa p ( mm H g ) sv p ( mm H g ) l a p ( mm H g ) l v p ( mm H g ) r a p ( mm H g ) s a p ( mm H g ) r v p ( mm H g ) p v _ Time (s)
Time (s)
Time (s) C O ( m L / s e c ) p ( mm H g ) pa p ( mm H g ) sv p ( mm H g ) l a p ( mm H g ) l v p ( mm H g ) r a p ( mm H g ) s a p ( mm H g ) r v p ( mm H g ) p v _ (b)(c) (d)(e) (f) Data r r
29 29.5 3001020 29 29.5 30050100 29 29.5 308010029 29.5 3055.2 29 29.5 300510 29 29.5 300102029 29.5 301015 29 29.5 302.42.62.8 29 29.5 30585960 (a) C O ( m L / s e c ) p ( mm H g ) pa p ( mm H g ) sv p ( mm H g ) l a p ( mm H g ) l v p ( mm H g ) r a p ( mm H g ) p v _ p ( mm H g ) s a p ( mm H g ) r v Figure A2. Incremental changes in resistance and compliance
Incremental decreases in resistance and increases in compliance and their effects on mean pulmonary ar-tery pressure (mPAP, panel (a)) and mean systemic arterial pressure (mSAP, panel (b)). Note that each patient responds slightly differently to the changes in parameter values. m SAP ( mm H g )
0% 20% 40% 60% 80% R p R p & R s & m PAP ( mm H g ) R p C pa R s C sa R p C pa & R s C sa & (a)(b) R s
0% 20% 40% 60% 80% 0% 20% 40% 60% 80%0% 20% 40% 60% 80% 0% 20% 40% 60% 80% 0% 20% 40% 60% 80% C pa C sa C pa C sasa