Predictive Modeling of Secondary Pulmonary Hypertension in Left Ventricular Diastolic Dysfunction

Diastolic dysfunction is a common pathology occurring in about one third of patients affected by heart failure. This condition may not be associated with a marked decrease in cardiac output or systemic pressure and therefore is more difficult to diagnose than its systolic counterpart. Compromised relaxation or increased stiffness of the left ventricle induces an increase in the upstream pulmonary pressures, and is classified as secondary or group II pulmonary hypertension (2018 Nice classification). This may result in an increase in the right ventricular afterload leading to right ventricular failure. Elevated pulmonary pressures are therefore an important clinical indicator of diastolic heart failure (sometimes referred to as heart failure with preserved ejection fraction, HFpEF), showing significant correlation with associated mortality. However, accurate measurements of this quantity are typically obtained through invasive catheterization and after the onset of symptoms. In this study, we use the hemodynamic consistency of a differential-algebraic circulation model to predict pulmonary pressures in adult patients from other, possibly non-invasive, clinical data. We investigate several aspects of the problem, including the ability of model outputs to represent a sufficiently wide pathologic spectrum, the identifiability of the model's parameters, and the accuracy of the predicted pulmonary pressures. We also find that a classifier using the assimilated model parameters as features is free from the problem of missing data and is able to detect pulmonary hypertension with sufficiently high accuracy. For a cohort of 82 patients suffering from various degrees of heart failure severity, we show that systolic, diastolic, and wedge pulmonary pressures can be estimated on average within 8, 6, and 6 mmHg, respectively. We also show that, in general, increased data availability leads to improved predictions.

Diastolic dysfunction is a common pathology occurring in about one third of patients affected by heart failure. This condition may not be associated with a marked decrease in cardiac output or systemic pressure and therefore is more difficult to diagnose than its systolic counterpart. Compromised relaxation or increased stiffness of the left ventricle induces an increase in the upstream pulmonary pressures, and is classified as secondary or group II pulmonary hypertension (2018 Nice classification). This may result in an increase in the right ventricular afterload leading to right ventricular failure. Elevated pulmonary pressures are therefore an important clinical indicator of diastolic heart failure (sometimes referred to as heart failure with preserved ejection fraction, HFpEF), showing significant correlation with associated mortality. However, accurate measurements of this quantity are typically obtained through invasive catheterization and after the onset of symptoms. In this study, we use the hemodynamic consistency of a differential-algebraic circulation model to predict pulmonary pressures in adult patients from other, possibly non-invasive, clinical data. We investigate several aspects of the problem, including the ability of model outputs to represent a sufficiently wide pathologic spectrum, the identifiability of the model's parameters, and the accuracy of the predicted pulmonary pressures. We also find that a classifier using the assimilated model parameters as features is free from the problem of missing data and is able to detect pulmonary hypertension with sufficiently high accuracy. For a cohort of 82 patients suffering from various degrees of heart failure severity, we show that systolic, diastolic, and wedge pulmonary pressures can be estimated on average within 8, 6, and 6 mmHg, respectively. We also show that, in general, increased data availability leads to improved predictions.

INTRODUCTION
Diastolic heart failure (sometimes referred to as heart failure with preserved ejection fraction or HFpEF, see Table 1) is a serious, often fatal, cardiovascular pathology. Recent reviews (Obokata et al., 2020) report that this pathology is often associated with several comorbidities, making the selection of homogeneous treatment groups difficult. It is, however, commonly characterized by elevated left ventricular filling pressures and normal systemic circulatory indicators such as left ventricular ejection fraction, cardiac output, and mean arterial pressure (Bonow and Udelson, 1992). In 2013, heart failure was mentioned in one of every nine death certificates in the United States, and was the underlying condition in roughly 20% of these cases. The number of deaths attributable to heart failure was approximately as high in 1995 as it was in 2013, with hospital discharges remaining stable from 2000 to 2010 (Mozaffarian et al., 2016). It is also estimated that about one-third of the patients with congestive heart failure (CHF) have a normal left ventricular ejection fraction (LVEF, see e.g., Zile and Brutsaert, 2002).
It has been observed that pulmonary hypertension (PH) is highly prevalent and often severe in HFpEF and that both pulmonary venous and arterial hypertension contribute to the severity of HFpEF with a marked correlation between systolic pulmonary arterial pressure (sPAP) and mortality (Lam et al., 2009;Obokata et al., 2020). A recent paper reviewing the current understanding of etiology and treatment of HFpEF, reports that multiple non-diastolic abnormalities contribute to the syndrome of HFpEF, including LV systolic dysfunction, pulmonary hypertension, RV and LA dysfunction, vascular stiffening, ventricular interdependence, and coronary microcirculation dysfunction (Obokata et al., 2020). Additionally, homogeneous treatment of HFpEF is made difficult by the multitude of phenotypes induced by obesity, ischemia, and cardiometabolic abnormalities. This is consistent with the findings in Shah et al. (2015) where a machine learning-based approach is used to cluster patients with HFpEF in three phenogroups with common characteristics, providing a way to go beyond a homogeneous treatment of HFpEF that has so far produced unsatisfactory results. Even though these studies highlight the current challenges in the treatment of HFpEF, they both seem to agree on the strong association between PH and HFpEF. In this regard, Obokata et al. (2020) mentions how "PH is extremely common in HFpEF, seen in roughly 80% of patients, and mortality is increased in this cohort, " whereas in Shah et al. (2015), elevated PH was one of the main criteria for patient recruitment.
While non-invasive echocardiography and machine learning may be useful for phenotyping classification and treatment selection (Obokata et al., 2020), early diagnosis of HFpEF relies on invasive pressure acquisition through right-heart catheterization, often performed following the manifestation of symptoms (Fisher et al., 2009;Galiè et al., 2009). Therefore, it is evident that methods enabling accurate indirect estimation of pulmonary pressures using minimally invasive clinical data would be extremely beneficial for early diagnosis of HFpEF in a way that could trigger lifestyle changes that will, in turn, prevent other co-morbidities from developing. Thus, in this study, we investigate how the physics-based consistency of a lumped parameter hemodynamic model containing three compartments, i.e., a four-chamber heart, systemic and pulmonary circulation compartments, may be used to monitor PH in patients from non-invasive and uncertain clinical measurements. The development of computer models to study hemodynamics in humans started in the 1960s and 1980s (Snyder and Rideout, 1969;Avanzolini et al., 1985Avanzolini et al., , 1988Ursino, 1998), with application in pediatrics developed in the 2000s for single-ventricle congenital heart disease (Pennati et al., 1997), Norwood physiology (Migliavacca et al., 2001), and systemicto-pulmonary artery shunts (Pennati et al., 2010). Approaches for automatic parameter estimation date back to the late 1970s (Deswysen, 1977;Deswysen et al., 1980), ranging from two-stage Prony-Marquardt optimization (Clark et al., 1980) and adaptive control systems for left ventricular bypass assist devices (McInnis et al., 1985;Shimooka et al., 1991) to Kalman filters (Yu et al., 1998(Yu et al., , 2001 and recursive least squares (Ruchti et al., 1993). An iterative, proportional gain-based identification method is presented in Revie et al. (2013) and an application to coronary artery disease is discussed in Sughimoto et al. (2013). Other studies include estimation of three-element Windkessel boundary conditions (Spilker and Taylor, 2010) and left ventricular viscoelasticity (Cappello et al., 1987;Avanzolini et al., 1992). More recently, examples of automatic parameter tuning in lumped circulatory models have included the physiology of children with congenital heart disease undergoing the first stage (i.e., Norwood) of single ventricle palliation surgery (Schiavazzi et al., 2016), construction of optimally trained patient-specific models for coronary artery disease (Tran et al., 2017) and predicting time evolution of ventricular dilation and thickening (Witzenburg and Holmes, 2018). A study using lumped parameter models in diastolic heart failure is finally discussed in Luo et al. (2011), while parameter identification for a mice model with chronic hypoxia and drug-induced pulmonary hypertension is proposed in Tewari et al. (2013).
Circuit models in hemodynamics typically contain a large number of parameters which need to be trained from clinical records collected at multiple visits, including a variable but typically sparse number of clinical measurements. This aggravates the ill-conditioning of the inverse problem, where model outputs do not change in response to perturbations along a number of unidentifiable linear combinations of parameters. In these circumstances, optimization may not be successful in identifying global optima and sequential Monte Carlo techniques (Del Moral et al., 2006) may underperform in practice, as data typically represent extremes (maxima/minima) or mean values of clinical indicators over one heart cycle.
Two technological trends make the present contribution particularly timely. On the one hand, there is increasing importance attributed to the availability of large training datasets which is at the base of the current revolution in AI and deep learning (see, e.g., LeCun et al., 2015). This includes anonymous electronic health records (EHRs) for specific sub-populations affected by a common clinical condition. On the other hand, there is an increasing availability of computational resources on the cloud, creating a perfect infrastructure for distributed computing with lightweight models. For these reasons, we envision an increased adoption of numerical models as regularizers to determine physics-informed predictive distributions for missing data in EHRs, going beyond currently adopted, physics-agnostic multiple imputation methods (Schafer, 1999). This study aims to be the first step in this direction and provides the following new contributions: • We propose a systematic approach to train patient-specific circulatory models with clinical data uncertainty, and demonstrate the results obtainable on a modest cohort of 82 patients. • Explore optimal parameter training as a possible approach to increase the feature space in order to facilitate classification of cardiovascular anomalies.
In section 2.2, we discuss the differential formulations of a compartmental circulation model for human adults, including circuit elements, a generic heart model, and the governing equations for the aortic, systemic, and pulmonary compartments. This is followed by an analysis of two datasets in section 2.3, the first used for validation, while the second contains EHRs for 82 patients. Our numerical investigation is articulated through answers to the questions formulated in section 2.4, using the numerical algorithms and tools briefly introduced in sections 2.5 and 2.6. Section 3 highlights the results of both the model analysis and the predictive results of the model. This section addresses the physiological admissibility of the selected model, the ability of the model to capture dysfunction mechanisms, and the sensitivity and identifiability of input parameters. On the predictive side, this section addresses the predictive performance of pulmonary pressures, non-pulmonary target ranking, and viability of modelbased PH classifiers. Conclusions are discussed in section 4. Finally, for convenience, a list of acronyms is provided in Table 1.

Compliance With Ethical Standards
This study was classified as research not involving human subjects and was approved on June 13th, 2019, by the Office of Research Compliance and Institutional Review Board at the University of Notre Dame under IRB#19-05-5371. This work utilizes data resulting from external studies that involved human participants. The medical procedures that occurred in these studies were performed in accordance with both the ethical standards of the institutional and national research committee and with the 1964 Helsinki declaration and its later amendments. As this study is a retrospective study, formal consent from the involved human participants was not required.

Modeling Approach
Blood circulation in adults can be simulated through a lumped parameter model (sometimes also referred to as zero-dimensional or 0D model). The foundation beneath a lumped parameter model lies in the equations that govern the voltage and current in an electrical circuit. The equations that underlie the voltage and current of an electrical circuit stem from the principles of the conservation of energy. Utilizing the hydrodynamic analogy (see, e.g., Rideout and Dick, 1967), these equations can also be used to govern the pressure and flow of blood where voltage corresponds to pressure and current corresponds to flow. In this study, we consider patient-specific 0D representations containing seven compartments: the four chambers in a biventricular heart, a compliant aortic arch, and the pulmonary and systemic circulations. The pulmonary compartment is represented through an RC circuit, as shown in Figure 1.

Generic Heart Model Compartment
The four heart chambers are represented by a series arrangement of a pressure-volume generator, an inductor, an ideal unidirectional valve, and a resistor, following prior work. The atrial and ventricular pressure-volume generators are formulated through a combination of an activation function and active/passive pressure curves (Avanzolini et al., 1985; FIGURE 1 | Lumped parameter hemodynamic model with RC pulmonary circuit. Migliavacca et al., 2001), where the index i refers to either the right-left or atrialventricular chamber, i ∈ {ra, rv, la, lv}. The active pressure curve is assumed linear and characterized by an active elastance E i,act and an unstressed chamber volume V i,0 . Additionally, passive volumes and pressures are related through an exponential relation, characterized by two elastance coefficients K i,pas,1 and K i,pas,2 (Avanzolini et al., 1985;Migliavacca et al., 2001). This is the simplest formulation compatible with the finite strains experienced by the ventricle during filling (see, e.g., Mirsky, 1976). The pressure-volume generator is completed by an atrial and ventricular activation functions A a,i , A v,i of the form where t sa = t c · t sas and t sv = t c · t svs are the atrial and ventricular relative activation times, respectively, and t c = 60.0/HR is the heart cycle duration in seconds. Additionally, t mv = t − ⌊t/t c ⌋ t c is the time from the beginning of systole (start of the cardiac cycle), t pw = t c /t pws is the atrial relative activation time delay and t ma = (t + t sa − t pw ) − (t + t sa − t pw )/t c t c . The model inputs are t sas , t svs and t pws . The chamber volume is determined through the equations where φ i , i ∈ {T, P, M, A} are valve activation functions for the tricuspid (T), pulmonary (P), mitral (M), and aortic (A) valve, respectively. These are equal to one for a negative pressure gradient through the valve and zero otherwise. Moreover, valves are modeled as perfect, without accounting for possible leakage or regurgitation. For models including valve prolapse and consequent regurgitation the reader is referred to, e.g., Pant et al. (2016). An inductance element located downstream with respect to the pressure-volume generator simulates the inertia of the blood in the chamber, according to the differential equation where Q i is the volumetric blood flow going through the i-th chamber, P i and P dn the pressure in the i-th and downstream chamber, respectively, R i the viscous resistance located between chamber i and chamber dn, and L i an inductance parameter.
Additionally, the selected model is fully capable of representing the physiologic consequences of a stenotic valve, as an increase in the resistance of the associated compartment would accentuate the pressure drop across the valve and lead to an increase of the upstream pressure. Finally, we remark that systolic and diastolic functions are separately represented in this model, using three parameters for each chamber, i.e., E i,act , K i,pas,1 , and K i,pas,2 . Identification of these parameters from clinical health records would therefore be informative of systolic and diastolic chamber function.
Note that the RL parameters of each cardiac chamber remain fixed in this study (see Supplementary Table 1). In other words, we assume that valvular stenosis has been excluded as a possible cause of pulmonary hypertension, for example, exclusion through a non-invasive echocardiographic assessment.

Aortic and Systemic Compartment
An aortic compartment consisting of an isolated capacitor is positioned downstream of the left ventricular outflow, modeled through an equation of the form where Q up and Q dn are the volumetric flow rate from the left ventricle and abdominal aorta, respectively, P ao is the aortic pressure and C ao the aortic compliance. We compare the value of P ao computed by this model with the clinically acquired brachial pressure. An RCR circuit simulates the systemic circulation, with C sys used to represent the overall systemic compliance, while two resistors simulate the viscous resistance in arteries and veins R sys,a and R sys,v , respectively. The algebraic-differential equations for the systemic compartments are therefore:

Pulmonary Compartment
The pulmonary circulation is represented through a RC circuit with equations where the pulmonary, left atrial and right ventricular pressures are denoted by P pa , P la , and P rv , and pulmonary capacitance and resistance are C pa and R pa , respectively. The pulmonary flow rate is denoted as Q pa , while Q rv,pa indicates the flow across the pulmonary valve, having activation equal to,

Initial Conditions
Initial conditions are specified for all state variables in the system of ODE. These include ventricular and atrial chamber volumes (V rv,ini , V lv,ini , V ra,ini , V la,ini ). Additionally, the model contains one state variable for every capacitor and inductor. Inductors are located at each valve, so initial conditions must be specified within the model for the initial flow across the tricuspid (Q ra,rv ), pulmonary (Q rv,pa ), mitral (Q la,lv ), and aortic (Q lv,ao ) valve. Finally, an initial pressure is specified for capacitors located in the aorta (P ao ), pulmonary (P pa ) and systemic (P sys ) circulation.

Available Data Sets
Two data sets are used throughout this study. Synthetic patientagnostic clinical measurements representing increasing severity of diastolic left ventricular dysfunction are used initially, while anonymized patient-specific electronic health records (EHR) for a cohort of 82 patients are utilized in the second part of the study.

Validation Data Set
The validation data set ( Table 2) contains the mean and standard deviation of thirteen different clinical targets for three different heart failure groups: healthy patients, mild heart failure patients, and patients with severe heart failure. Normal physiologic clinical targets were determined from the literature (Edwards Lifesciences Corporation, 2009), while targets associated with severe left ventricular diastolic dysfunction were assigned with the supervision of a clinician. The values for mild heart failure patients were obtained through a linear interpolation between severe dysfunction and healthy conditions. The selected targets for severe HF conditions are characterized by a normal SBP and DBP accompanied by a slight reduction in CO.

EHR Data Set
Completely anonymized patient-specific clinical measurements for 82 adult patients were provided in the context of a research project funded by Google through its ATAP initiative, focusing on Modeling Non-invasive Measurements of Cardiovascular Dynamics. There are 26 clinical data targets, which are listed below in Table 3. Missing data were present with the pattern highlighted in Figure 2, with patients having zero to nineteen of the clinical targets. Two of the 84 patients in the dataset did not have any of the relevant clinical targets and were therefore excluded from the study (see the two patients with zero available targets in Figure 2A). The remaining 82 patients had between one and nineteen clinical targets. All patients but one had measurements for three clinical targets: heart rate, diastolic blood pressure, and systolic blood pressure. No single patient had all 26 measurements. Finally, the standard deviations for each target are also shown in Table 3. Their value was determined through a preliminary literature review (see,  e.g., Gordon et al., 1983;Maceira et al., 2006;Yared et al., 2011). A histogram of clinical data occurrences is illustrated in Figure 2A and displays three main modes. Most patients have either four, eight, or 17 available clinical measurements, likely due to the data aggregation produced by screening protocols. Additionally, a heat map of the EHRs is illustrated in Figure 2B. Each row of the heat map was normalized to a zero to one range to highlight the relative magnitude of the clinical target. Note how the size of the selected cohort (84 patients) is modest but sufficient to investigate the effectiveness of Bayesian inference in the context of a simple physiologic model. Acquisition of larger data sets is possible but made non-trivial by the need to automatically extract large volumes of clinical targets from text reports written following echocardiographic and catheter lab examinations. In such cases, using natural language processing tools is key to make larger EHR data sets available for research.
Finally, this data set focuses on cases of secondary PH, where a reversible increase of PVR follows an increase in left ventricular filling pressures. Therefore, this study does not consider primary PH and does not make any claim of differentiating primary from secondary PH.

Methodological Approach
This study is articulated through a number of logically consequential questions driving our numerical experiments. These questions are: 1. Physiological admissibility of 0D representations under normal and heart failure conditions-Are the model outputs able to reproduce sets of clinical targets ranging from healthy to pathological conditions? In other words, is the identification problem well-posed, in the sense that model outputs are able to represent a wide spectrum of conditions from health to disease? We answer this question in section 3.1.
2. Ability to model distinct diastolic/systolic dysfunction mechanisms-Is the selected model formulation able to separately represent the systolic and diastolic functions of the heart muscle? And does the alteration of these properties produce expected modifications in the physiology represented through model outputs? We answer this question in section 3.2. 3. Parameter sensitivity and identifiability-Once a set of quantities whose prediction is of interest (e.g., pulmonary arterial pressure) has been identified, do they show nonnegligible sensitivity with respect to changes in the parameters associated with physiologically relevant mechanisms affecting these quantities? Moreover, are these parameters identifiable so that it is possible to uniquely estimate their distribution from the available clinical data? Is this estimate robust (i.e., characterized by a limited uncertainty)? We answer this question in section 3.3. 4. Predictive ability of optimally trained models-Are models trained from clinical data other than the pulmonary pressures able to predict such pressures, and what is their accuracy? We answer this question in section 3.4.

Relative importance of non-pulmonary clinical targets-
Which minimal set of clinical targets should be collected to guarantee accurate model predictions for systolic, diastolic pulmonary pressure and pulmonary venous wedge pressure?
In other words, we would like to rank the non-pulmonary targets, starting with those producing maximally accurate predictions on the pulmonary arterial pressure. This allows for identification of a minimal set of maximally informative clinical quantities for predicting specific model outputs. We answer this question in section 3.5. 6. Detecting pulmonary arterial hypertension from assimilated circulation models-Group II pulmonary arterial hypertension is typically detected by a mean pulmonary pressure higher than 25 mmHg or a systolic pulmonary pressure higher than 35 mmHg (Simonneau et al., 2013). Instead of direct characterization based on clinical data, would it be possible to detect pulmonary arterial hypertension by classification from the parameters of a model trained with non-invasive measurements? Once assimilated, a model can be used to generate a large number of features, leading to a higher dimensional space with possibly improved separability (Cover, 1965). We answer this question in section 3.6.

Inference
Consider a set of m measured clinical targets represented through the random vector d ∈ R m with each component In other words, each quantity d i , i = 1, . . . , m measured in the clinic has marginal density ρ i (d i ), i = 1, . . . , m, and we assume all these measurements to be independent, i.e., their joint probability factors. We design a physiological 0D circuit model with n parameters y ∈ R n , so its outputs match the observed targets or, in other words, we introduce a statistical model of the form and assume each noise component ǫ i ∼ N (0, σ 2 i ) to follow a zero-mean Gaussian distribution. Note that the ith realization from the parameter vector y is denoted by y (i) and that the i-th model output is denoted by o i = G i (y), while the vector o ∈ R m contains the complete set of model outputs.
Each model is trained using two different approaches, i.e., by determining a maximum a posteriori estimate of the parameters y using repeated Nelder-Mead optimization (Nelder and Mead, 1965), and by solving an inverse problem through adaptive Markov chain Monte Carlo sampling, specifically, through the differential evolution adaptive Metropolis algorithm (Vrugt et al., 2009;Vrugt, 2016) and assessing convergence through the Gelman-Rubin diagnostic (Gelman and Rubin, 1992). In both cases, the posterior distribution P(y|d) is obtained by combining a uniform prior P(y) (see the admissible parameter ranges in Supplementary Table 1) with a Gaussian likelihood, P(y|d) ∝ P(d|y) · P(y),

Computational Tools
The TULIP (Tools for Uncertainty quantification, Lumped modeling and Identification of Parameters) software framework was developed to answer the above research questions. Tulip is a OOP C++ code designed to simplify the task of estimating parameters of lumped models for human circulation and contains abstractions for computational models, operations performed on these models (e.g., optimization, Bayesian estimation, local and global sensitivity analysis, etc.) and data sources used to store the available clinical targets. For an overview of the procedures for statistical data assimilation used in this study, the interested reader is referred to Schiavazzi et al. (2016) and Akintunde et al. (2019).

Physiological Admissibility
The circulation model discussed in section 2.2 was first trained on the validation dataset in section 2.3.1. We then generated a collection of model outputs using a subset y (i) , i = 1, . . . , 5, 000 parameter realizations from the converged MCMC samples, and compared the resulting distributions (after post-processing with Gaussian kernel density estimation) with the distributions assumed for the targets d. Specifically, we assumed that each clinical target follows a normal distribution with the mean and standard deviations listed in Table 2. The Kullback-Leibler (KL) divergence (Kullback and Leibler, 1951) was used to determine the agreement between the model-based predictions and measurements. The KL divergence measures how much information is lost when one uses the predicted, instead of the assumed, target distribution and a small KL divergence suggests physiological admissibility. As shown in Figure 3, the KL divergence is negligible for all targets, and the relative percent difference between mean target value and MAP model outputs is also generally small (see Table 4).

Inverse Assessment of Dysfunction Mechanisms
We now focus on the ability of our model to distinguish between diastolic and systolic dysfunction mechanisms when trained based on conditions reported in Table 2. Specifically, we aim to determine whether the selected model parameterization is able to separately represent the relaxation or contraction (i.e., diastolic or systolic) function of the heart muscle. Recall that diastolic function, and therefore left ventricular stiffness, during relaxation relate to the linear and exponential passive curve parameters K i,pas,1 and K i,pas,2 in equation (1). These two parameters directly relate to the condition that we aim to assess. Figure 4 shows the mean diastolic and systolic chamber function parameters and associated 10-90% confidence intervals, grouped by anatomical relevance (i.e., left ventricle, right ventricle) and plotted for healthy patients and patients with mild and severe heart failure, respectively. Comparison between heart failure conditions allows one to assess how model parameters change due to disease progression. This change is shown in Figure 4, where many of the parameters related to pulmonary FIGURE 3 | KL divergence between predicted and assumed clinical targets for varying heart failure severity.  hypertension change as patients progress from healthy to mild, and mild to severe diastolic dysfunction. The results show that the passive left ventricular curve slope K lv,pas,1 and the pulmonary resistance increase from healthy to mild and mild to severe PH. Additionally, changes in the left ventricular active elastance parameter E max,lv remain limited, i.e., the model correctly predicts a practically unaltered systolic function. This confirms the ability of the physiological model used in this study to link PH with diastolic dysfunction and increased pulmonary resistance. Since Figure 4 combines quantities with different units, we have also looked at the left ventricular active and passive pressure-volume slope at V 0 . The active curve slope at V 0 is equal to E max,lv , whose MAP estimates are 3.01, 2.23, 2.02 Barye/mL for healthy, mild, and severe PH, respectively. The passive curve slope at V 0 is instead K pas,lv,1 * K pas,lv,2 (see Equation 1) equal to 0.06, 0.12, 0.15 Barye/mL. Therefore, from healthy to mild PH the active slope changes by −25.8%, while the passive slope changes by +113.1%. From mild PH to severe PH, the active slope drops by −9.4%, while the passive slope increases by +24.3%. This perspective supports our conclusion that our model consistently captures increased severity in underlying PH through more relevant variations in the diastolic rather than systolic left ventricular properties. Finally, the significant variability in the unstressed ventricular volumes is explained by the presence of multiple local peaks in the posterior distribution, i.e., reasonable PV loops compatible with the clinical targets produced by wildly different unstressed volumes.

Sensitivity and Identifiability of Circulation Model Parameters
Results from the above sections confirm the well-posedness of the selected model formulation for the full spectrum of diastolic dysfunction, from mild to severe. We now focus on determining the most relevant model parameters that significantly alter the main quantities of interest, particularly the systolic, diastolic, and pulmonary wedge pressures. In addition, we study both the structural and practical identifiability in an effort to determine unimportant parameters and their nonidentifiable combinations.

Average Local Sensitivities
Non-dimensional local sensitivities are computed for all outputs as so that we consider the relative change in model outputs that correspond to a 1% variation in each parameter. The maximum a posteriori parameter vector y map and the corresponding model outputs o map = G(y map ) are computed from MCMC for each of the 82 patients in the cohort, and used to compute the sensitivities in Equation (11). The resulting sensitivities are then averaged across all patients. Figure 5A illustrates the average sensitivities obtained by training our model with the complete list of clinical targets, including systolic, diastolic, and venous wedge pulmonary pressures and pulmonary vascular resistance. We note that large sensitivities are apparent for the heart rate across all outputs while at the same time, accurate measurements of heart rate are easy to obtain non-invasively. Additionally, to check how the sensitivities were affected by the availability of pulmonary pressure targets (i.e., the very same quantities we would like to predict), we re-computed sensitivities using parameter estimates y map obtained by excluding the pulmonary pressure targets during training. The average sensitivities appear to be minimally affected by the selective exclusion of pulmonary pressure targets. Figure 5B suggests the following important parameters. Changes in the heart rate or t svs directly affect the amount of blood flow ejected by the ventricles, altering the mean pulmonary pressures under a constant PVR and initial condition (i.e., P pa,ini ). As already discussed in the above sections, the left ventricular diastolic pressure/volume ratio K lv,pas,1 and the associated exponential factor K lv,pas,2 govern the diastolic properties of the left ventricle, while the E max,lv is instead responsible for the systolic function. The left atrioventricular and aortic valve resistance R la,lv and R lv,ao , respectively, govern the pressure drop from the left atrium to the left ventricle and from the left ventricle to the aorta. These two parameters therefore affect the left atrial and ventricular pressures and, in turn, the upstream pulmonary pressures. Mitral or aortic valve stenosis are typical examples of this mechanism (see, e.g., Tracy et al., 1990). The pulmonary resistance and capacitance parameters, R pa and C pa , clearly affect the mean pulmonary pressures and their range. In contrast, changes in systemic vascular resistance, R sys,a and R sys,v , affect the left ventricular afterload and, in turn the pulmonary pressures.

Structural Identifiability
Structural identifiability analysis is performed to gain an understanding of our ability to recover a given set of model parameters through the solution of an inverse problem from idealized noiseless clinical targets that belong to the model range. This contrasts with the analysis in the next section of the practical identifiability, where real clinical data are used instead. We solve the model for the default parameter combination (see Supplementary Table 1) and regard the outputs as data, from which a MAP estimate of the default parameter values is recomputed by MCMC followed by NM optimization. We would like to point out that the parameters in Supplementary Table 1 correspond to a healthy patient. Therefore no attempt is made in this section to represent hypertensive physiological conditions. Additionally, histograms are generated using 5,000 samples from the MCMC parameter traces and compared to the default (true) parameter set. The results for right ventricular, left ventricular, and resistor parameters can be seen in Figures 6A-C, respectively. The true parameters are always found within the parameter distributions from MCMC and often correspond to the mode of the histogram. Some deviations may be observed in Figure 6C regarding the arterial and venous systemic resistance. In such a case, since SVR is the sum of such resistances,  the increase in the first is compensated by a reduction in the second.
The coefficients of variation for the model parameters marginal posteriors and associated learning factors are shown in Figure 7. While parameters with a higher coefficient of variation have a greater spread relative to their mean, the learning factor quantifies how much the marginal variance is reduced by conditioning the model output to the available observations or, in other words, how much the marginal variance is reduced from the prior to the posterior (Schiavazzi et al., 2016). Figure 7 shows that the parameters with the largest coefficient of variation are also the parameters with the smallest learning factor, as expected. Heart timing parameters, parameters for systemic and pulmonary resistance and compliance, the aortic compliance parameter, and the active ventricular parameters are generally well-learned.
The analysis is completed by a comparison between the true and optimal physiology, calculated using 5,000 parameter combinations from the MCMC traces, and analyzed over a single heart cycle. The results for the aortic, pulmonary, and left ventricular pressures and flows, and the left and right ventricular pressure-volume loop are reported in Figure 8. Pressures, flows, and volumes agree well with those generated from the default parameter set.

Practical Identifiability
We also perform local identifiability analysis through the Fisher information matrix rank (Rothenberg, 1971) to determine  the presence of non-identifiable parameter combinations and unimportant parameters. To do so, we compute the matrix ∂ G(y map )/∂ y of local derivatives for our output quantities o map = G(y map ) with respect to the parameters y. The Fisher Information matrix I(y map ) can be computed as where the precision matrix B contains the inverse target variances, i.e., B i,i = 1/σ 2 i and B i,j = 0 for i = j. Rank deficiency total radar plots of all parameters whose eigenvalues are less than the selected cut-off. (C) Example of unimportant initial condition, i.e., whose perturbation has no effect on the model results. (D) Example of non-identifiable parameter combinations where no dominant parameter can be identified and where the combination significantly changes across patients.
in I(y map ) reveals the presence of non-identifiable parameter combinations as illustrated in Figure 9A by plotting the Fisher information matrix (FIM) eigenvalues ordered by magnitude for all patients. Small eigenvalues are observed across all analyzed patients, confirming a lack of local identifiability.
Eigenvectors for all patients whose corresponding eigenvalues were less than a selected cut-off were superimposed on the same radar plot to search for situations characterized by dominant components, representative of unimportant parameters. This process is visualized in Figure 9A, illustrating all eigenvalues colored uniquely per patient and Figure 9B, showing the effect of changing the cut-off value on the number of selected patients. Specifically, it shows that no significant changes result by adopting cut-offs in the range [1 × 10 −12 , 1 × 10 −16 ]. Figures 9C,D shows an example of two of the 17 radar plots generated for this study. The plot on the left shows an example of an unimportant parameter with dominant eigenvalue associated with the parameter V la,ini , i.e., the initial left atrial volume. Instead, the radar plot on the right does not show any clear pattern involving parameter combinations that significantly change across patients. This local identifiability analysis confirms how initial conditions for pressures, flows, and volumes (except P pa,ini , P ao,ini , and P sys,ini as per the results of the previous sensitivity analysis) are generally unimportant.

Prediction of Pulmonary Pressures
We now focus on the problem of using the physiological consistency of our compartmental model to predict pulmonary pressures from other possibly non-invasive clinical targets. To do so, we have trained our models without including the pulmonary targets (i.e., systolic, diastolic, wedge pressures, and pulmonary vascular resistance) and propagated forward the estimated parameters to quantify the marginal distributions of these targets. We then evaluated the error between predicted and true pulmonary pressures together with their variability. To do so, we first show four Bland-Altman plots for sPAP, dPAP, mPAP, and PVR, respectively (Figures 10A,B). Minimal bias is generated when training the model with the pulmonary pressures and PVR. Predictions generated by pulmonary pressure-blind training FIGURE 10 | (A,B) Bland-Altman plots of sPAP, dPAP, mPAP, and PVR for clinical targets, and as predicted from models trained with and without pulmonary targets, respectively. Solid horizontal lines represent the mean of all differences, while dashed lines are drawn at 1.96 times the standard deviation. Predictions associated with gray markers in the Bland-Altman plots result from using less than 15 clinical targets. (C) Predictive performance for pulmonary pressure. Absence of PAP targets in average prediction errors is represented using dashed lines. The shaded region represents the area bounded by the 5th and 95th percentile from 5,000 random subsamples from MCMC. (D) Zoom on average errors which correspond to patients with more than 15 available clinical targets.
does not seem to generate proportional bias except for the predicted PVR.
In Figure 10C, we also plot the average absolute pressure error, with associated uncertainty, across 31 patients (those for which the pulmonary pressure was available and characterized by having more than six REDCap entries) versus the minimum number of prescribed clinical targets. Finally, as illustrated in the closeup shown in Figure 10D, which focuses on patients with at least 15 REDCap entries, the average errors on the predicted pressures is around 8 mmHg for systolic PAP, and 6 mmHg for Diastolic PAP and PCW.

Relative Importance of Non-pulmonary Targets
In this section, we investigate which clinical targets are the most important to include during training, in order to minimize errors in pulmonary pressure predictions. In other words, we rank clinical targets starting from those having a more beneficial impact on the accuracy in predicting pulmonary hypertension. We achieve this goal through a sequence of optimization steps. We begin by performing training using optimization with a single target at a time. The target found to minimize the average  combined prediction error for pulmonary pressures is ranked first and permanently added to the list of targets included in all successive optimization steps. The percent error was computed for each clinical pulmonary target that was available for each patient. Percentage error expresses the percentage difference between the model output and the clinical target. A cumulative percentage error was calculated by taking the sum of the percentage errors for each pulmonary target clinically available for each patient. To avoid bias toward patients with fewer available pulmonary targets, the cumulative percentage error is transformed to an average combined prediction error through division by the number of available pulmonary targets. Let p ∈ {1, 2, 3, 4} be the number of pulmonary targets available in the data of a certain patient. The percent error of clinical target s i,target , i = 1, . . . , p is computed as e i = 100 · |s i,target − s i,output |/s i,target , i = 1, . . . , p, while the average combined prediction error is p i=1 , e i /p. All remaining targets are iteratively tested, and those producing the minimum average combined PAP errors (shown in Figure 11A) are progressively ranked until all targets have been considered. The process above is repeated for each patient. Figure 11B shows the resulting average ranking and associated occurrences, i.e., the number of patients where a specific target was collected.
The list of targets ordered by average rank and filtered by occurrence is also shown in Table 5. Most of the quantities (except mPAP, PCW, PVR, and SVR) can be estimated noninvasively. From RAP (or CVP, RVEDP), SBP-DBP, HR, and SVR it is possible to estimate the cardiac output. From CO, HR, and LVEF, it is possible to estimate LVEDV, which is correlated with PCW. This might explain why PCW is always estimated better than sPAP and dPAP. The mitral valve deceleration time and velocity E/A ratio are indicators of diastolic LV function and therefore correlated with PAP. If PVR, CO, and PCW are known, it is possible to estimate mPAP, which is correlated to sPAP and dPAP. Finally, the parameter ranking is summarized for all patients in Figure 12.

PH Classifiers From Assimilated Circulation Models
This section explores the use of trained lumped parameter models for the automatic detection of abnormal pulmonary pressures from minimally invasive clinical measurements. We employ a naïve Bayes classifiers (see, e.g., Lewis, 1998) to detect pulmonary hypertension from our dataset. The first step of this classification process was to generate a ground truth variable for hypertension that could be used to analyze the accuracy of our hypertensive classification. Using the patient's clinical data, a binary hypertensive variable was defined using the criteria from Simonneau et al. (2013). If a patient had mPAP > 25 mmHG or sPAP > 35 mmHG, then the patient was classified as hypertensive. Of the 82 patients, 65 of the patients contained sufficient clinical data to define a ground truth binary hypertensive variable. The remaining 17 patients did not contain enough clinical data to determine a ground truth value for hypertension, so these patients were excluded from the testing and training datasets.
Before a naïve Bayes classification can be performed, the problem of missing data must first be addressed. Our dataset representing the cohort of 82 patients contains a non-negligible ratio of missing data, with patients missing between 1 and 19 of the 24 total clinical targets. To overcome this issue and to verify how classification results depend on the strategy selected for missing data imputation, five different missing data imputation approaches were tested. The first method tested was complete-case analysis, considering only the 4 variables that were available for all patients, i.e., heart rate, systolic blood pressure, diastolic blood pressure, and cardiac output. The remaining four missing data methods consisted of replacing the missing values with (1) zeros, (2) the max value of the data set, (3) a value far outside the range of the data set (10 times the max), or (4) the median. In addition to the training of five separate classifiers using five different missing data methods, a sixth naïve Bayes classifier was trained on the MAP parameters of the model. The data for the MAP parameters from the model is a complete set, so no missing data method was required. Figure 13 shows how the above imputation approaches affect classification accuracy. At first sight, the high accuracy produced by a multiple imputation strategy using the maximum inter-patient clinical value (or 10 times its value) may seem surprising. This can be explained by observing how, in such a case, the probability of the feature associated with the missing data is practically zero, thus essentially reducing the number of features and increasing the resulting accuracy. We then grouped patients according to a training and a testing set, consisting of 4/5 and 1/5 of the available data, respectively. Training and testing is performed either using the raw clinical data or model parameters assimilated though optimization and MCMC, resulting in ∼88% accuracy for the classifier trained with the 24 raw clinical data points and about 76.4% accuracy for the classifier trained with the 45 assimilated parameter values (see Table 6).
Of the 82 patients, 65 contained enough information for ground truth classification of whether the patient should be identified as hypertensive. Based on these data, we were able to classify 22 of the 65 patients, or about one-third of the patients, as being affected by pulmonary hypertension and the remaining two-thirds with normal PAP pressures. This imbalance is known to introduce bias (Sun et al., 2009), with the classifier more likely to label patients as non-hypertensive (as 2/3 of the training data are non-hypertensive). As a remedy, six different approaches were applied to deal with this unbalance, which represents a mixture of over-sampling methods, under-sampling methods, and a combination of both. Figure 14 shows the principal component decomposition and contingency table for the unbalanced data and the data balanced with centroid clustering. Note how the confusion matrices are evaluated based on test data, hence 13 patients (1/5 of the 65 patients) for the unbalanced dataset. Undersampling in centroid clustering reduces the training/testing data from 65 to 33: hence one-fifth of the 33 samples becomes ∼7, as FIGURE 13 | Accuracy of a naïve Bayes classifiers for pulmonary hypertension using different approaches for multiple imputation. Frontiers in Physiology | www.frontiersin.org FIGURE 15 | Accuracy of a Naïve Bayes classifiers for pulmonary hypertension, using identified lumped parameters as features, and various methods to handle unbalanced data. The methods used were Random Over-sampling (ros), Random Under-sampling (rus), an over-sampling method: SMOTE (Synthetic Minority Oversampling TEchnique) (sm), two under-sampling methods: Tomek Links (tl) and Cluster Centroids (cc), and a hybrid approach that performs over-sampling followed by under-sampling using SMOTE followed by Tomek Links (smt). shown in the confusion matrix for the centroid clustering undersampling method. Figure 15 shows the effect of each unbalanced data method on the overall area under the receiver operating characteristic curve (ROC) for the classifier trained on the model parameters. For such a classifier, centroid clustering (Lin et al., 2017) increases accuracy from 76.4% up to 100%. In addition, data balance using Synthetic Minority Oversampling TEchnique (SMOTE) (Chawla et al., 2002) increases the accuracy of the model trained on the data from 88 to 100%, as shown in Table 6. We again remark that the number of patients in this study is small, and future studies will investigate the generalization of this approach to larger datasets.

CONCLUSION
This study demonstrates that a relatively simple lumped parameter compartmental model can represent a wide range of physiologies, spanning healthy patients to patients affected by severe diastolic left ventricular dysfunction. In this context, an activation formulation for the heart chambers has proven important to separately account for systolic and diastolic pressure-volume behavior. When trained using ideal clinical data from subjects with diastolic left ventricular dysfunction, the parameters associated with the governing physiological mechanisms (i.e., diastolic ventricle relaxation) change in a way that is consistent with the underlying physiology of the dysfunction. The average parameter sensitivities are determined after training with real data from 82 patients, confirming intuition for the most important parameters. Pulmonary pressures were found to be primarily sensitive to the following parameters: heart rate and contraction timing parameters, diastolic and systolic pressure/volume ratio parameters (K pas,1,lv , K pas,2,lv , and E max,lv ), left atrioventricular and aortic valve resistance (R la,lv , R lv,ao ), pulmonary resistance and capacitance (R pa , C pa ) and systemic resistance (R sys ). Additionally, the model was found to be locally unidentifiable, with initial conditions generally unimportant. The structural identifiability analysis showed that inference of model parameters is feasible under perfect, noiseless conditions and with a sufficiently large number of available clinical targets. Target ranking based on sequential optimization reveals the most important non-invasively acquired clinical targets: heart rate, systemic pressures, peak pressure gradient across aortic and pulmonary valves, pulmonary valve acceleration time, mitral valve deceleration time, and mitral valve E/A peak ratio. The clinical targets that positively affect the prediction of pulmonary pressures, but require an invasive practice for their measurement, are central venous pressure, right ventricular systolic and diastolic pressure, systemic vascular resistance, cardiac output, and mean right atrial pressure. After investigating parameter identifiability/sensitivity, the pulmonary pressures of the 82 patients with various heart failure severities were predicted from a lumped hemodynamic model, which was trained based on the remaining clinical targets. The average absolute pressure error on the 11 patients characterized by at least 11 distinct clinical entries was found ≃8 and 6 mmHg for systolic and diastolic/wedge pulmonary pressures, respectively. While these errors may seem large to detect cases of mild hypertension, they may be more than reasonable depending on how high are the PA pressures we are trying to predict/detect. In other words, while an 8 mmHg pressure error may seem relevant with an actual pulmonary pressure of 20 mmHg, for cases of patients with the real disease (i.e., PA pressures of 40, 50, 60 s, etc.), then the difference between (say) 50 and 58 would not be nearly as large. However, more importantly, our approach would still identify those patients as patients that would require further evaluation.
Finally, we have shown that MAP estimates of circulation model parameters can be used to detect elevated pulmonary pressures, and that a simple classifier provides high accuracy on balanced data, even when these parameters are identified without measuring systolic, diastolic, wedge pulmonary pressures, and pulmonary vascular resistance. Future work will be devoted to the systematic demonstration of the proposed approach on larger patient data sets, opening new avenues for translational applications of model-based diagnostics, improving model formulations that target specific diseases, and developing improved and more efficient estimation approaches.

DATA AVAILABILITY STATEMENT
The datasets analyzed for this study plus a Python/Cython implementation of the adult circulation model are available at https://github.com/desResLab/supplMatHarrod20.

ETHICS STATEMENT
The study was classified as research not involving human subjects and approved on June 13th, 2019 by the Office of Research Compliance and Institutional Review Board at the University of Notre Dame under IRB#19-05-5371. All procedures performed in studies involving human participants were in accordance with the ethical standards of the institutional and/or national research committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.