STUDY OF UNCERTAINTIES ON A 0D MODEL OF THE SYSTEMIC CIRCULATION

: Cardiovascular system is intensely researched to understand the intricate nature of the heart and blood circulation. Nowadays we have well evolved computational models which are useful in many ways for the understanding and analysis of physiological and pathophysiological conditions of the heart. However, the practical use of these models and their results for clinical decision making in speciﬁc patients is not straightforward. In this context, models predictions must be accurate and reliable, which can be assessed by quantiﬁcation of uncertainties in the predictions and sensitivity analysis of the input parameters. Lumped parameter models for the cardiovascular physiology can provide useful data for clinical patient-speciﬁc applications. However, the accurate estimation of all parameters of these models is a difﬁcult task, and therefore the determination of the most sensitive parameters is an important step towards the calibration of these models. We perform uncertainty quantiﬁcation and sensitivity analysis based on generalized polynomial chaos expansion in a lumped parameter model for the systemic circulation. The objective of this work is to verify the effect of uncertainties from input parameters on the predictions of the models and to identify parameters that contribute signiﬁcantly to relevant quantities of interest. Numerical experiments are performed and results indicate a set of the most relevant parameters in the context of these models. is the volume of blood in the ventricle at the end of ﬁlling. The stroke volume (SV) is the volume of blood pumped from the left ventricle per beat and is given by SV = EDV − ESV . The ejection fraction (EF) is the portion of total blood volume ejected from the left ventricle and is calculated by EF = SV/EDV . Lastly, the ﬁnal QOIs calculated are the maximum and minimum arterial pressures (Max AP and Min AP, respectively). We only computed metrics for the upper body arterial pressure, since that is the QoI commonly useful in clinical scenarios.


INTRODUCTION
The evolution of methods for obtaining experimental data, invasive or non-invasive, generates a large amount of physiological information that is increasingly more accurate and detailed. Such comprehensive data enables the construction of more complex and robust models, capable of performing simulations of cardiovascular system, cardiac electrophysiology, or dynamics related to blood flow (BESSONOV et al., 2016), whose predictions bring up relevant questions and answers for a better understanding of the cardiovascular system both in healthy and pathological conditions. Nowadays, there is an understanding that the predictions provided by these models will be more and more used in medical research and clinical applications in a personalized manner.
In particular, information related to cardiovascular hemodynamics, are very important in clinical settings, since they can be used to obtain essential physiological measures of great relevance for treatment and therapy design. A simple way of obtaining such hemodynamic measures is through the use of lumped parameter (LP) cardiovascular models. In fact, through a representation analogous to electrical circuits, LP models are designed to represent the cardiovascular dynamics and can be as simple as a three-compartment windkessel model of the arterial system, and as complex as a model for the global hemodynamics with over 153 parameters (LIANG;LIU, 2005). Thanks to their capability of representing the cardiovascular dynamics under various physiological and pathological conditions, the simulations of LP models find many clinical and patient specific uses (DUANMU et al., 2018;WILLIAMS et al., 2014), although these models, especially the most complex ones, are commonly developed for research purposes (SHI; LAWFORD; HOSE, 2011).
In a clinical context, where patient-specific complexities are relevant, a model that accurately describes the patient's individual reality is sought. One way to capture patient specific information to serve as model inputs is by identifying model parameters aided by measured observations. However, most of the times, clinical data are inevitably hampered with measurement uncertainties. A better understanding of the model behavior is necessary to increase its robustness and accuracy. Great insights into the model functioning can be achieved through uncertainty quantification (UQ) of model predictions that is inherently associated to the uncertainties in the inputs. In addition, sensitivity analysis (SA) is useful to find out which uncertain inputs are responsible for the more variations in the outputs.
In this work, we apply the generalized polynomial chaos method to perform global UQ and SA on a five-compartment LP model of the systemic circulation in order to study parametric uncertainty propagation and to determine the most influential parameters with respect to quantities of interest usually employed in clinical applications. The information acquired with Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 2, p. 221-01, 221-13, 2020 DOI: 10.21575/25254782rmetg2020vol5n21152 this study can help grasp the overall performance of the model and provide useful insights into model calibration, simplification, and robustness enhancement.

CARDIOVASCULAR MODEL
The model described here was presented in the work of (WILLIAMS et al., 2014) to examine cardiovascular regulation during the head-up tilt (HUT) protocol and more recently in a study to predict patient-specific blood flow and pressure during the HUT protocol (WILLIAMS et al., 2019). As shown in Figure 1, it models the dynamics of systemic circulation by estimating the volume (V ) in ml and pressure (p) in mmHg in the upper and lower systemic arteries (subscripts au and al, respectively) and veins (subscripts vu and vl, respectively), as well as the left heart (subscript lh). It also predicts the blood flow (q) in ml/s between upper and lower vessels, in the upper and lower peripheral arteries (subscripts aup and alp, respectively) and in the mitral and aortic valves (subscripts mv and ao, respectively).  Figure 1: Systemic circulation representation of the 5-compartment model.
Pressure-volume relation for the systemic vessels can be defined as where V i and V i,us denote total and unstressed volumes, respectively, in compartment i = {au, al, vu, vl}. Here, C i (given in ml/mmHg) is the compartment compliance and p i is the instantaneous blood pressure. For the left heart, the pressure-volume equation is given by where E lh (t) is the time-varying elastance function that provides the pumping mechanisms of Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 2, p. 221-01, 221-13, 2020 DOI: 10.21575/25254782rmetg2020vol5n21152 the heart chambers, given by wheret is the time within a cardiac cycle of total duration T . Maximum and minimum elastances are represented by the parameters E M and E m , respectively, while the times at which they occur are given by T s and T r , respectively. Flow between vessels follows where p in and p out are the pressures in the upstream and downstream compartments, and R (mmHg s/ml) is the resistance of the vessel. The opening and closing mechanisms of the mitral and aortic valves are modeled by assuming the flow to zero when the valve is closed and when it is opened the flow follows Equation (4). Finally, the rate of volume change in compartments is denoted by Once the volumes are calculated on Equations (5), and using the appropriate flows q in and q out , pressure and flow for each compartment are obtained using Equations (1), (2) and (4).

Parameters, Initial Conditions and Numerical Solution of the Model
The parameters and initial conditions considered for the simulation of the fivecompartment model are taken from the work of (WILLIAMS et al., 2019), which were obtained from experimental data and model fitting. Table 1 shows all the parameter values and initial conditions used in the simulations.
The opening and closing mechanism of the valves introduces a discontinuity, which characterizes the system of differential equations as stiff (MARQUIS et al., 2018). It is essential for stiff differential equations to be solved with appropriate methods, since these problems can introduce several numerical complications. Also, when solved with explicit methods they perform poorly and require extremely small time steps. An efficient approach to solve stiff problems is the backward differentiation formula (BDF) methods (WANNER; HAIRER, 1996). Hence, to solve the LP equations, we used the BDF method implemented in the SciPy library of the Python programming language. Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 2, p. 221-01, 221-13, 2020 DOI: 10.21575/25254782rmetg2020vol5n21152

ANALYSIS
We assume that the inputs are subject to uncertainties due to measurement errors and as a consequence the outputs are also affected by these uncertainties. The functional representation of the model is given by: where Y is the stochastic output and Z = [Z 1 , Z 2 , . . . , Z D ] is the stochastic input vector which is a random multivariate variable, where D is the number parameters. The uncertain parameters Z i are assumed independent and its joint probability distribution is denoted by ρ Z (z). The output Y is also a random variable with an associated probability density function ρ Y , which depends on Z and is governed by the model f . To describe the uncertainty in the output Y the following statistical moments are used: the expected value µ, the standard deviation σ and the coefficient of variation CoV , that is: where V [Y ] is variance of Y , given by The SA is performed using Sobol's main and total indices (SOBOL, 2001) to quantify the contributions of uncertain inputs Z i to the output Y and their interactions. The first-order or main sensitivity index gives the direct effect of the uncertain input Z i by quantifying the portion that Z i contributes (without interactions) to the total variance of V [Y ]. It is given by: where E [Y |Z i ] represents the expected value of the output Y for a fixed value of Z i , with 1 ≤ i ≤ D. The main sensitivity index is useful for input prioritization, since it gives the expected reduction in the total variance that could be achieved if the unknown true value of Z i could be measured. The total sensitivity index is the sum of all first and higher order effects (interactions between inputs) where the uncertain input Z i is involved. It is given by where Z −i is the set of uncertain inputs except Z i . The total sensitivity index represents the variance linked to the effect of all interactions of Z i as well as the direct effect of Z i .

Polynomial Chaos
To obtain the measures of uncertainty and sensitivity presented in Equations (6)-(11) we used the generalized polynomial chaos method (gPC) (ECK et al., 2016;XIU, 2010). In this method the output Y is expanded into a series of orthogonal polynomials, which are functions of the uncertain inputs Z. The orthogonal polynomials form a basis that span the output space of Y , which is expressed as: where c e are the expansion coefficients, Φ e are the orthogonal polynomials of order p and N p = (D+p)! D!p! is the number of basis functions. The polynomial basis is chosen according to the input parameter distribution (XIU, 2010). To determine the polynomial for Y the coefficients c e must be computed. The probabilistic collocation method (PCM) (TATANG et al., 1997) is a non-intrusive approach to compute these coefficients and was employed here. After the Revista Mundi, Engenharia e Gestão, Paranaguá, PR, v. 5, n. 2, p. 221-01, 221-13, 2020 DOI: 10.21575/25254782rmetg2020vol5n21152 computation of the coefficients, the UQ and SA metrics can be calculated by directly evaluating Equations (6), (7), (10) and (11) with the gPC expansion (12). In this work we used the ChaosPy (FEINBERG; LANGTANGEN, 2015) library implementation of the gPC method.
The accuracy of the uncertain and sensitivity measures produced by the gPC method depends on the polynomial degree p and the number of samples N s . Normally choosing the number of samples as N s = dN p , where d > 1 yields better accuracy (ECK et al., 2016) and requires the solution of a least squares problem. To check the accuracy of the surrogate model Y ≈ f gPC (Z) we used the leave one out cross validation (LOO) approach (BLATMAN; SUDRET, 2010). The LOO cross validation approach takes one sample out of the N s samples available and then builds the gPC with N s − 1 samples. This is performed for i = 1, . . . , N s , leaving at each time one different sample out of the construction of the polynomial. LetŶ i be the model output computed by the i-th gPC approximation, whereas Y i is the forward model output for i-th sample left out. The predicted residual is given by ∆ i = Y i −Ŷ i . The approximation error (Err LOO ) is the mean value of the squared predicted values and the determination coefficient Q 2 to evaluate the accuracy of the gPC approximations are given by: The determination coefficient Q 2 takes a value between 0 and 1, and the closer it is to 1, the more accurate is the surrogate model.

Quantities of Interest
We computed UQ and SA measures for the left heart pressure and volume and considered six scalar QoIs for the study of the model. The end systolic volume (ESV) is the volume of blood in the ventricle at the end of contraction. The end diastolic volume (EDV) is the volume of blood in the ventricle at the end of filling. The stroke volume (SV) is the volume of blood pumped from the left ventricle per beat and is given by SV = EDV − ESV . The ejection fraction (EF) is the portion of total blood volume ejected from the left ventricle and is calculated by EF = SV /EDV . Lastly, the final QOIs calculated are the maximum and minimum arterial pressures (Max AP and Min AP, respectively). We only computed metrics for the upper body arterial pressure, since that is the QoI commonly useful in clinical scenarios.

RESULTS
The simulation results for the LP model are presented next. First, we present the predictions obtained using the baseline values from Table 1. Then, we show the results of the UQ and SA analysis performed for the model.
The simulation performed considered 10 cardiac cycles of 0.98 s. Results shown in Figure 2, with baseline values for the parameters, depict the blood pressure in all the compartments, together with the blood volume in the left heart and the blood flow through the aortic and mitral valves. These results are in accordance with those obtained by (WILLIAMS et al., 2019). To carry out UQ and SA analysis, the following set of parameters were considered as uncertain inputs: Z = {R aup , R alp , R al , R vl , R mv , R av , C au , C al , C vu , C vl , E M , E m , T s , T r }.
Since we have no previous information about the probability distribution of these inputs, we set all parameters with a normal distribution with a coefficient of variation of 20% around their baseline values (Table 1).
Considering the computational demand of using a large polynomial degree p and number of samples N s , we started with a low polynomial degree and a small number of samples, until we found a combination of p and N s that yielded a desired accuracy for the gPC expansion. Thus, we performed the LOO cross-validation for a polynomial with order p = 2. Since the set of uncertain parameter Z is comprised of 14 parameters, the gPC expression has N p = 120 polynomial base functions. Using N s = 2N p samples, the LOO test yielded a Q 2 value next to 1 for all QoIs (see Table 2), showing that the gPC expression obtained yields a desirable convergence with this configuration, which was used in this work.
The UQ results are shown in Table 2. The EF presented a variation of 6%, ESV presented a 23% of variation, while all the other QoIs had from 14% to 18% of variation.  whereas ESV, SV, EF and Min AP are asymmetrical. Mean values for a healthy male young adult (between 20 and 29 years old) are shown for all QoIs, together with the 95% confidence interval of the same type of subject for the ESV, EDV, SV and EF. These data are from the study of (MACEIRA et al., 2006), which did not provide the 95% confidence interval for the arterial pressures.

CONCLUDING REMARKS
In this work, a LP cardiovascular model of the systemic circulation which describes the circulation through the left heart, veins and arteries of the upper and lower parts of the body was implemented. UQ and global SA techniques based on the gPC method were applied to the model, to verify the effects of the uncertainties on the model's input parameters on the quantities computed by it. The results showed that the ESV was the most affected QoI by the uncertainties in the model inputs. Still, the distributions of ESV, EDV, SV and EF are in accordance with the physiological ranges of a young adult subject. In addition, the sensitivity indices calculated for each uncertain input showed that the maximum and minimum elastances (E M and E m ), as well as the compliance of the upper vessels (C au ) and the upper peripheral arterial resistance (R aup ) are essential parameters for the QoIs studied here. This means that small uncertainties in these model input parameters can cause great variations in the predictions obtained by the model. It is interesting to observe the result of the sensitivity analysis, since almost all the QoIs considered are directly related to the functioning of the left ventricle and the most sensitive parameters are all influential in the dynamics of the left ventricle.
The methodology for quantifying uncertainties and global sensitivity analysis presented in this work was not applied to the proposed five-compartment model, in the work of (WILLIAMS et al., 2019). Thus, it is expected that the results obtained in the present study can increase knowledge about the model in order to make it robust for use in clinical contexts with data from specific patients.