Practical Identifiability and Uncertainty Quantification of a Pulsatile Cardiovascular Model
Abstract
Mathematical models are essential tools to study how the cardiovascular system maintains homeostasis. The utility of such models is limited by the accuracy of their predictions, which can be determined by uncertainty quantification (UQ). A challenge associated with the use of UQ is that many published methods assume that the underlying model is identifiable (e.g. that a onetoone mapping exists from the parameter space to the model output). In this study we present a novel methodology that is used here to calibrate a lumpedparameter model to left ventricular pressure and volume time series data sets. Key steps include using (1) literature and available data to determine nominal parameter values; (2) sensitivity analysis and subset selection to determine a set of identifiable parameters; (3) optimization to find a point estimate for identifiable parameters; and (4) frequentist and Bayesian UQ calculations to assess the predictive capability of the model. Our results show that it is possible to determine 5 identifiable model parameters that can be estimated to our experimental data from three rats, and that computed UQ intervals capture the measurement and model error.
keywords:
Cardiovascular dynamics; modeling; parameter estimation; uncertainty quantification; patientspecific modeling.1 Introduction
Precision medicine is a growing model of healthcare that proposes to customize of care, medical decisions, practices, and products to each individual patient. This approach is important, as pathologies such as cancer, autoimmune disorders, and cardiovascular diseases are unique to a given individual making it challenging to develop diagnostic and treatment protocols. One approach, to studying patientspecific complexities, is to use mathematical modeling to estimate function and predict features that are difficult to measure, thus providing a more comprehensive set of information to distinguish between individual patients.
A rich history of cardiovascular modeling exists in the literature, typically written either from a fluid dynamics perspective (resulting in systems of PDEs) FormaggiaCVBook (); QuarteroniCVBook (); van2011pulse (), to study local flow properties, or from a compartment perspective (resulting in systems ODEs) to study systems level dynamics OttesenOlufsenCVBook (); ottesen2013development (). This study focuses on analysis of compartment models predicting systems level propagation of flow and pressure through the cardiovascular system Blanco10 (); Kokalari13 (); Shi11 (). In this model type, compartments represent groups of vessels (e.g. large or small arteries or veins, capillaries, or vessels supplying specific tissues or organs) coupled to a heart compartment that act as a pump to drive the system. Some models include both pulmonary and systemic circulations neal2007subject (), while others analyze one of the two systems zinemanas1994relating (). This model can be used to extract vascular properties such as vascular resistance, cardiac contractility, or compliance by fitting models to pressure and flow data from noninvasive imaging studies Olufsen05 (); Pope09 (); Williams14 () and/or from invasive catheterization Revie13 (); pacher2008measurement (); mackenzie1979effects () studies.
One of the biggest challenges in calibrating compartment models to data is obtaining accurate parameter estimates. Even in its basic form, where the model is formulated using systems of linear differential equations, forced by a contracting heart, it is typically not possible to uniquely estimate all model parameters. To overcome this, we propose to use sensitivity analysis and subset selection for apriori study of the model structure followed by parameter estimation and uncertainty quantification. In general, parameters that are unidentifiable as a result of model structure are referred to as structurally unidentifiable mahdi2014structural (), whereas parameters that are unidentifiable as a result of practical restrictions, such as availability and quality of data, are referred to as practically unidentifiable miao2011identifiability (). Theoretically, structural identifiability is a prerequisite for practical identifiability. However, in practice it can be difficult to establish the former, since analysis is restricted to models for which it is possible to define a unique inputoutput relation mahdi2014structural ().
Only a few studies have addressed structural identifiability in cardiovascular models. Kirk et al. Kirk13 (), studying Windkessel models, showed that three of four parameters are identifiable, and Pironet et al. Pironet16 () demonstrated that every parameter in a linear sixcompartment model including a left and right heart, systemic and pulmonary arteries and veins are structurally identifiable if outputs contain both pressure (in all arteries and veins) and left/right ventricular stroke volume, while models relying on either pressure or volume alone are structurally unidentifiable. Other studies have employed sensitivity and practical as opposed to structural identifiability analysis in models predicting arterial blood pressure and cardiac output Ellwein08 (); Pope09 (); Williams14 (); Gul16 (). Several recent studies have addressed uncertainty quantification, mostly for analysis of 1D fluid dynamics models, however to our knowledge, these methodologies have not previously been applied to analysis of compartment models. The study by Eck et al. Eck16 () develops a guide to uncertainty quantification in cardiovascular models presenting a number of methodologies. Several studies have predicted uncertainties in specific onedimensional fluid mechanics models Cheng13 (); Arnold16 (); Eck15 (); eck2015stochastic (); eck2016guide (); Paun2018 (). Of these, three studies accounted for uncertainty using Kalman filtering Cheng13 (); Arnold16 (); Eck15 (), two used polynomial chaos expansion, and one Paun2018 () used an MCMC approach based on the Delayed Rejection Adaptation Metropolis (DRAM) algorithm haarioDRAM (). To our knowledge, no study has combined these techniques into an organized workflow for the determination of model parameters in compartmental CV models given a specific data set.
In this study, we present a potentially general multistage methodology to establish reliable parameter estimation approaches and uncertainty quantification (UQ) for lumpedparameter cardiovascular models which is used here to characterize left ventricular volume and blood pressure data from three Sprague Dawley rats. The key steps in our methodology include: (1) the use of literature and available data to compute nominal parameter values specific to each rat; (2) sensitivity analysis and subset selection to determine a set of identifiable parameters; (3) optimization to compute point estimates for the identifiable parameters; and (4) statistical techniques to quantify uncertainty of the model output.
2 Methods
2.1 Experimental Data
Data analyzed here are extracted from experiments performed on 3 SpragueDawley (SD) rats (2 male, 1 female). The average weight of these animals was g. Rats were anesthetized with sodium pentobarbital (50 mg/kg, ip), and catheters were placed in a femoral vein and artery for administration of anesthetics and monitoring of systemic blood pressure respectively. A pressurevolume conduction catheter (Millar SPR869, 2F tip with four electrodes and 6mm spacing) was inserted through the right carotid artery into the left ventricle to simultaneously obtain pressure and volume measurements. For each rat basic physiological measures (sex, weight, heart rate, average stroke volume and cardiac output, Table 1) were recorded along with continuous measurements of left ventricular volume and pressure. For this study, approximately 20second pressurevolume timeseries, measured at rest, were selected for model identification and the final 0.5 seconds of each data set was used to calibrate the model, shown in Figure 1.
Rat  Sex  Weight  Heart rate  Stroke volume  Cardiac output 

(g)  (beats/min)  (l)  (ml/min)  
Rat 1  Male  339  2403  3081  740.2 
Rat 2  Male  350  2403  2161  520.2 
Rat 3  Female  342  4203  1431  600.2 
VolumeConductance Calibration
Conduction catheters measure blood volume indirectly by measuring the conductance through the changing volume of blood in the left ventricle during the experiment baan1981continuous (). This method traditionally requires the infusion of a hypertonic saline solution and a cuvette calibration at the end of the experiment to convert the raw voltage time series data into a left ventricular volume time series. However, in these experiments, perturbations performed between the saline calibration and cuvette measures caused this calibration approach to be inaccurate. So to obtain volumes with a standard physiological range, the cuvette measures were scaled against literature estimates pooling end diastolic and stroke volume data from 20 previously published studies using conduction catheters, ultrasound, and magnetic resonance imaging (MRI) of the left ventricle in SD, WistarKyoto, and Lewis rat strains (from 191 animals total). These data are shown as a function of animal body weight in Figure 2. Since MRI based volume measurement techniques have been adopted as the goldstandard for ventricular volume studies kjaergaard2006evaluation (), one can observe that measurements from nonMRI methods tend to under estimate both the end diastolic and stroke volume. To obtain realistic volumes, we fit 13 MRIbased end diastolic and stroke volume data sets to body weight using exponential functions
(1) 
where BW is the body weight, EDV is the end diastolic volume, SV is the stroke volume, and are estimated parameters. Predicted values of EDV and SV from (1) are calculated for each animal in our study. The raw conduction (volume) signal is recorded in volts, therefor the predicted EDV is used to convert the peak voltage for each cardiac cycle to the maximum volume and the predicted SV is used to scale the amplitude of the voltage signal.
2.2 Model
Similar to previous studies Williams14 (), we use a fivecompartment model to predict pressure, flow, and volume in the left ventricle, the large and small arteries and veins; see Figure 3. Using an electrical circuit analogy, blood flow () is analogous to current, pressure () to voltage, volume () to charge, vessel resistance () to electric resistance, and vessel elastance () to the reciprocal of capacitance. For each compartment, dynamics are predicted from three equations relating pressure, flow and volume. The flow in and out of each compartment is proportional to pressure via Ohm’s law
(2) 
the pressure and volume in each compartment is related to elastance via
(3) 
where is the external tissue pressure, is the unstressed blood volume (both assumed constant), and is the stressed blood volume; and conservation of volume gives
(4) 
For the circuit shown in Figure 3, we derive a system of five differential equations for the stressed volume of the form (4), detailed in the Appendix.
The beating of the heart (Figure 4) is modeled by a periodic timevarying elastance function Ellwein08 () defined over one cardiac cycle of length as
(5) 
where and denote the minimum and maximum elastance, respectively, of the left ventricle. denotes time for end systole and the time at which the heart has relaxed to its diastolic value.
We model the two heart valves using its electrical equivalent, a diode, i.e.
(6) 
where , representing the aortic (av) and mitral (mv) valves, respectively.
In summary, the model takes the form
where denotes the model states, the model parameters, and the model output.
2.3 Nominal parameters and initial values
Nominal parameters and initial conditions for the ODEs can be obtained from analysis of data and known physiological features extracted from literature. In this study, parameters include resistance, elastance, and timing of the cardiac cycle, which are determined as a function of resting blood volumes, pressures, cardiac output, and heart rate. In the following sections we discuss how to calculate a priori values for all model parameters, listed in Table 2.
2.3.1 Blood volume
Total blood volume for healthy adult Wistar rats is 57l/g body weight trippodo () (listed in Table 1). Measurements analyzed in this study are from 2 male and 1 female healthy Sprague Dawley rats (Charles River Laboratories, USA). Following suggestions by Young young () and Gelman gelman (), the total volume for the purpose of nominal value estimation is distributed with 2.5% in the large arteries, 7.5% in the large veins, 20% in the small arteries and 70% in the small veins to 70%, i.e. for each compartment, the volume is given by
(7) 
where is the total blood volume, refers to the th compartment volume, and denotes the percentage of the total blood volume to compartment .
As noted above, the model is formulated in stressed volume, which is a small percentage of the total blood volume. Literature estimates for the total stressed volume vary significantly, from about 10% to 40% beneken (); trippodo (); gelman (); young (); Gotwals00 (). We found no resources describing distribution of stressed and unstressed volume across organs or between arteries and veins in rats, therefore we used the same stressed volume fraction (30%) in all compartments, i.e.
(8) 
where denotes ”left heart” (lh), ”large arteries” (la), ”small arteries” (sa), ”small veins” (sv), or ”large veins” (lv).
2.3.2 Pressure
Measurements of left ventricular pressure is used to approximate pressures in all vascular compartments. Starting at the large arteries, the aortic valve open when the pressure in the left ventricle exceeds the aortic pressure. The maximum pressure in the systemic arteries follow the one in the left ventricle. Assuming a constant pressure drop from left ventricle through the small arteries we let
where denotes the left ventricular pressure data, denotes maximum larger artery pressure, and denotes maximum small artery pressure. Without pressure measurements in the systemic arteries, more assumptions are needed to approximate the mean arterial pressure. Assuming a pulse pressure of mmHg, large and small artery diastolic pressure can be estimated as , and from which we can estimate mean arterial blood pressure as
using the common clinical approximation KlabundeWeb ().
This model does not represent the right ventricle or pulmonary vasculature so the pressure in the systemic veins becomes the filling pressure for the left ventricle. Thus similar to the arterial side, we assume a 10% pressure drop between large and small veins, so we let
2.3.3 Cardiac output
For each rat, the average stroke volume (Table 1) is extracted from measurements of maximum and minimum left ventricular volume. Cardiac output can subsequently be calculated as the product of stroke volume and heart rate (i.e. ).
2.3.4 Vascular resistance
Using the baseline cardiac output and pressure, the vascular resistances can be computed from Ohm’s law (2). For example, arterial resistance is predicted as
(9) 
Instead of using the mean pressures to calculate valve resistance we consider the maximum and minimum pressure to estimate the resistance across the aortic and mitral which can then are expressed as
(10) 
where is the aortic valve resistance and is the mitral valve resistance. Since the venous pressure do not vary significantly, the average pressure of the large veins is used to predict .
2.3.5 Vascular compliance
Each vascular compartment is associated with an elastance (constant). Assuming that the tissue pressure is negligible (e.g. ), elastance can be estimated as
(11) 
following the pressurevolume relation (3), set up for the large arteries as a representative example.
2.3.6 Heart parameters
The elastance function (5) has four parameters, including timing parameters denoting the length of the cardiac contraction () and relaxation (), as well as the minimum () and maximum () elastance. Rat unstressed ventricular volume has been estimated at 37l london2012 (). For each rat, the timing parameters and can be extracted from data, as shown in Figure 4. denotes the time at which the left ventricular volume reaches it maximum (at ), and is the time at which reaches its baseline after contraction. The minimum elastance is associated with enddiastole, where the left ventricular pressure is minimal and ventricular volume is maximal, and the maximum elastance is associate with systole, where the ventricular pressure is maximal and ventricular volume is minimal. These considerations let us set
(12) 
2.3.7 Initial conditions
Assuming that the model simulation begins at the end of systolic contraction (beginning of diastolic filling), we set the initial value of the volume in each compartment to their stressed volume (). This implies that
(13) 
where is the initial volume of the th compartment.
Quantity  Equation  Reference  Quantity  Equation  Reference 
young (); gelman ()  data  
KlabundeWeb ()  
young (); gelman ()  data  
KlabundeWeb ()  
young (); gelman ()  data  
young (); gelman ()  data  
(10)  (11)  
(10)  (11)  
(9)  (11)  
(9)  (11)  
(9)  
(12)  
(12) 
3 Model Analysis
The model described in Section 2.2 is linear with respect to the states but nonlinear with respect to the parameters. This gives rise to multiple parameter interactions that inevitably complicate the parameter estimation process. While a model with many nonlinear parameter interactions can be structurally identifiable, unidentifiable parameter relationships often result in an illconditioned optimization problem ipsen2011rank ().
In this section, we outline a systematic methodology for parameter estimation, comprising identification of

Nominal parameters from literature and available data followed by a baseline simulation to ensure an appropriate model response. For the model analyzed here this step requires two parts: i) using data to compute subject specific nominal parameter values, ii) determine a shift in data to ensure that model predictions and data are in phase.

Local sensitivities used to study how the parameters influence the model output.

Structured correlations used to determine a subset of parameters with minimal parameter interactions.

Parameter estimates obtained using the LevenbergMarquardt optimization method, estimating the subset of identifiable parameters minimizing the least squares error between the model output and available data.

Frequentist prediction and confidence intervals used to quantify uncertainty in the model solutions.

Parameter distributions and credible intervals obtained using the DRAM (Delayed Rejection Adaptive Metropolis) algorithm.
While the analysis is devised for a relatively simple cardiovascular model with a specific set of output data (left ventricular pressure and volume), the approach introduced here is applicable to any predictive model fitted to timevarying data. Obtaining the nominal parameter of step one is model specific however, steps 26 are more generic and the approach presented here can be more generally applied
3.1 Local Sensitivity Analysis
Sensitivity analysis quantifies how the model output changes in response to changes in parameter values Ellwein08 (). In this study we use derivativebased sensitivities to quantify the local influence of the model output on each parameter. Similar to the study by Pope et al. Pope09 (), we computed sensitivities from partial derivatives of the model output residual with respect to each model parameter, i.e. we define the sensitivity matrix as
(14) 
where is the th parameter and is the th time step. When estimating parameters in a logtransformed space (14) can be expressed as
The matrix can be calculated analytically for simple models; however, numerical approximation of is more practical for complex models. Here we employ finite differences to approximate by
where is chosen to reflect the precision of the model output and is the unit vector in the th component direction. If the error in the model evaluation (ODE solver error tolerance) is on the order of , the step should be to get an error of the same magnitude in the sensitivities Pope09 (). We tested the stability of our finite difference approximation by reducing the ODE solver tolerance and observing that the results of the sensitivity analysis converged the same values (results not shown). To get a rough approximation of how a parameter influences model behavior we use ranked sensitivities , defined as the twonorm over each column of
3.2 Subset Selection: Structured Correlation Analysis
The model considered here has inherent parameter interactions that necessitate the need for selecting parameter subsets with minimal unidentifiable parameter interactions. While various subset selection algorithms exist (e.g. Miao11 ()), we employ the structured correlation method by Ottesen et al. OlufsenCorr () to construct a set of identifiable model parameters. According to this method, a pair of parameters with a large correlation (and strongly coupled uncertainty) cannot both be uniquely estimated. This method systematically removes the least sensitive correlated parameters parameters until an identifiable set remains. The input to this method is the sensitivity matrix  which means that any correlation is only valid within some neighborhood of the parameter values used to construct the sensitivity matrix.
Using the sensitivity matrix in (14), the covariance matrix is given by the inverse of the Fisher Information Matrix (FIM)
(15) 
where . Note that can only be inverted if has full rank. Linearly dependent columns of are a result of parameters being perfectly correlated, meaning that a parameter can be algebraically expressed in terms of other parameters. Additionally, columns of insensitive parameters () can lead to having a large condition number and thus should be removed a priori from the sensitivity matrix input. The entries of are used to compute the correlation matrix with entries given by
is a symmetric matrix with ones along the diagonal. Parameter pairs with a correlation value greater than (user defined threshold) are denoted correlated. Correlations between parameters show how the parameter values depend on each other when fitting experimental data from the same system.
3.3 Nonlinear Least Squares Optimization: LevenbergMarquardt
The goal of nonlinear least squares optimization is to estimate the set of parameters that minimizes the difference between the model output and the data , assumed to be some function of the model output corrupted with additive noise; e.g.
(16) 
where denotes the th data point, the model response at the th time point, the observation function mapping the model variables to the measured states, and the normally distributed observation error. In this study, extracts the left ventricular volume and pressure timeseries.
To fit our model to the data, we use generalized nonlinear least squares to determine a parameter set , which minimizes the sum of squares cost function
(17) 
The residual vectors are defined as
(18) 
The superscripts “m” and “d” denote the model and data, respectively. Dividing the residuals by the amplitude of the data ensures that the optimization procedure gives equal weight to both model outputs.
To estimate we employ the LevenbergMarquardt optimization routine by Kelley kelley (). To obtain a well scaled problem, we estimate the natural log of the true parameter values. As part of the optimization routine we discard parameter values that are 20 times larger or smaller than the nominal estimate. Like other gradient based optimization routines, the choice of the initial parameter vector is important. If the cost function has multiple minima the optimization routine may end in an undesirable or unrealistic local minimum. Using the physiologicallyjustified nominal parameter values described in Section 2.3 helps to place our closer to a physiologically realistic minimum. To verify that the optimization converged we randomly perturbed the nominal values by 10% and observed that the parameters all converged to the same values.
3.4 Uncertainty Quantification
Uncertainty quantification (UQ) is the process of determining uncertainties in estimated model parameters given uncertainties in the model formulation and experimental measurements (the inverse problem), as well as establishing how uncertainties in model inputs (such as parameters) affect the model output (forward propagation of uncertainty). In this study, we utilize UQ procedures from both frequentist and Bayesian statistical frameworks. We calculate confidence intervals and Bayesian credible intervals to measure the precision of the model in predicting the mean response. These approaches are outlined below; for more details, see SmithUQ ().
3.4.1 Frequentist Approach
Frequentist uncertainty propagation methods are computationally inexpensive compared to their Bayesian analog. Most frequentist statistics are derived from asymptotic assumptions assuming that the uncertainty distributions take a Gaussian shape, whereas Bayesian methods make no assumption about the shape of the underlying distributions. One of the main benefits of the asymptotic assumptions is that uncertainty distributions can be expressed as explicit formulas in the frequentist perspective. Frequentist confidence intervals can be calculated from
(19) 
where is the student tdistribution with degrees of freedom ( is the number of data points and is the number of parameters), is the estimate of the model standard deviation , is the sensitivity matrix, and is the th row of stacked as a column vector. Frequentist prediction intervals can be calculated in a similar manner by
(20) 
3.4.2 Bayesian Approach: Delayed Rejection Adaptive Metropolis (DRAM)
While frequentist methodology is fundamentally rooted in quantifying uncertainty in terms of repeating the data generating procedure, Bayesian inference is conditioned on a single data set; this allows for uncertainty about parameters to be expressed by probability distributions. In the Bayesian framework, represents a vector of random variables. Given observations , Bayes’ formula
(21) 
describes how the posterior density relates to the prior density , encompassing any a priori information known about the parameters, and the likelihood of observing the data for the model given . The marginal density of the data typically functions as a normalization factor and can be determined by
(22) 
Under the hypothesis (16), the likelihood function is given by
(23) 
where denotes the least square cost defined by (17), is the number of data points, and is the model variance  the index denotes the variance for pressure () or volume (). Both and are parameters represented by probability distributions that are also estimated. With a known likelihood and prior density , it is possible to estimate the posterior density if the integral (22) in the normalizing constant can be estimated. While this route is possible, the evaluation of highdimensional integrals is a difficult and expensive, and is currently an active area of research; see, e.g., sparse grid methods smolyak1963quadrature (); bungartz2004sparse () and quasiMonte Carlo methods halton1960efficiency (); joe2003remark ().
An alternative is to use Monte Carlo integration to randomly sample from the density . Many suitable Markov chain Monte Carlo (MCMC) methods exist in the literature (see Andrieu2008 () for an overview). This study uses the DRAM algorithm haarioDRAM (), which combines two methods for improving efficiency of MetropolisHastings type MCMC algorithms: delayed rejection (DR) MiraDR2001 () and adaptive Metropolis (AM) HaarioAM2001 (). These Metropolistype methods are acceptancerejection algorithms that accept new parameter samples only if the likelihood of the new candidate is higher than the current sample. DR allows for additional proposals per step if the initially proposed step is not accepted, thereby increasing the acceptance rate and wellmixing of the sample. AM allows for updating of the covariance matrix based on the history of the sample, thereby helping the algorithm to make better proposals and move toward the correct posterior distribution faster, reducing the burnin period.
We use (23) as our likelihood distribution, uninformed diffuse Gaussian distributions centered around the optimized point estimate from Section 3.3 as our prior distribution, and a Gaussian proposal distribution. To ensure a well scaled problem we have DRAM estimate the natural log of the true parameter densities like we do Section 3.3. Samples are taken from the DRAMestimated parameter probability distributions to compute Bayesian credible and prediction intervals; for more details, see haarioDRAM (); SmithUQ (). In this study, we utilize the MCMC toolbox by Haario et al. (2006) at http://helios.fmi.fi/lainema/mcmc/, which includes code for running DRAM as well as for computing Bayesian credible and prediction intervals.
4 Results
In this section, we use the stepbystep procedure described in Section 3 to systematically estimate an identifiable subset of parameters for the cardiovascular compartment model derived in Section 2.2 that fit left ventricular volume and blood pressure data described in Section 2.1.
4.1 Nominal parameter values
Given that sensitivity and subset selection analysis are determined at an estimate of the local minimum, the first step involves using available data and literature to determine the best possible nominal parameter values, as described in Section 2.3. Second, we manually shift the starting point of the data so the model predictions and data both begin at the same point of the cardiac cycle, as shown in Figure 5(a) and described in Section 3. This estimate of the data shift is refined in a later model optimization step.
4.2 Sensitivity Analysis and Subset Selection
Next, we employ local sensitivity analysis and structured correlation, along with physiological knowledge, to construct a subset of identifiable model parameters.
Figure 6 depicts normalized ranked parameter sensitivities and the subset of identifiable parameters (24) for all three animals. Ranking the parameters by sensitivity provides initial insight into what parameters are most influential in determining model behavior. Note that the ranked sensitivities of parameters and are orders of magnitude smaller than the other parameters, thus these parameters were removed a priori from the set of parameters analyzed. We used the structured correlation algorithm from Section 3.2 to obtain an identifiable parameter set. For all three rats we used a correlation threshold as an upper bound on the pairwise correlations between parameters and found the identifiable subset
(24) 
4.3 Optimization
Using the nominal parameters listed in Table 3, we optimize the parameter subset (24) using a LevenbergMarquardt optimization scheme; we tested convergence of the optimizer by varying the initial guess around the nominal values. Results were verified by randomly perturbing the nominal values by 10% and optimizing these perturbed parameters. For ten unique perturbations, all of the parameters converged to the same values (results not included). To ensure that the solution has settled to steady state, computations were done over entire 20 seconds of data, with the least squares cost evaluated over the last 0.5 seconds of the data shown in Figure 1. As part of the optimization procedure we also optimize the data shift which was manually determined in step 1. This was done by repeating optimizations for a large number of shifts. Results showed (5(b)) that there are multiple parabolic minimums for the cost reflecting the pulsatile nature of cardiovascular mechanics. We chose the minima that provided the smallest cost (17). Results depicted in Figure 5 shows the cost (17) and gradient plotted as a function of relative data shifts for Rat 1 (similar results were obtained for rats 2 and 3).
Optimization results for all three rats, with optimal shifts, are reported in Table 3 listing the nominal and optimized parameter values, and Figure 8, showing the model fit with optimized parameter values compared to the data. Note that we report the true parameter values and not the natural log of the parameters used by the LevenbergMarquardt optimizer.
4.4 Uncertainty Quantification
The pointestimates obtained using the LevenbergMarquart routine were used to initialize both frequentist and Bayesian UQ methods, to construct confidence and credibility intervals along with prediction intervals for the model output predictions. This choice was motivated by wanting to minimize the already high computational cost of solving stiff ODEs. To demonstrate our methodology is reproducible across data sets we repeated computations for all three data sets shown in Figure 1.
To perform UQ, we considered both the frequentist formulas stated in Section 3.4.1 and the Bayesian inference using DRAM as described in Section 3.4.2. Using the optimized values, we applied equations (19) and (20) to compute frequentist confidence and prediction intervals, respectively. The optimized point estimates were used to construct the prior distribution of parameters as diffuse Gaussian distributions (uninformed prior centered around the optimized point estimate) for the DRAM algorithm. This choice of a prior is justified by the fact there is no information we could use to construct a more informative prior distribution. Chains of 100,000 sample points were generated using DRAM. Figure 7 shows the resulting DRAMestimated parameter chains, parameter densities, and pairwise parameter correlations for rat 1, similar results were obtained for the other animals. Frequentist UQ intervals were computed from equations (19) and (20), and Bayesian UQ intervals were constructed by solving the model over randomly sampled values from the posterior parameter chains. Figure 8 compares the resulting frequentist and Bayesian UQ for all three animals.
Parameter  Units  Nominal  LevMar  DRAM Mean 

s mmHg/l  0.011, 0.014, 0.012      
s mmHg/l  0.942, 1.155, 1.003  1.600, 1.630, 2.196  1.588, 1.631, 2.187  
s mmHg/l  0.002, 0.001, 0.002      
l/mmHg  9.641, 8.446, 8.535      
l/mmHg  1.193, 1.045, 1.056      
l/mmHg  0.006, 0.003, 0.006      
l/mmHg  0.050, 0.027, 0.049      
s  2.500, 2.500, 2.500  4.544, 4.161, 4.166  4.518, 4.163, 4.213  
s  6.000, 6.000, 6.000  7.309, 6.409, 8.149  7.403, 6.409, 8.161  
mmHg/l  0.043, 0.026, 0.042  0.050, 0.025, 0.039  0.049, 0.025, 0.039  
mmHg/l  5.476, 3.864, 2.595  5.965, 4.366, 3.252  5.904, 4.364, 3.254 
5 Discussion
It is well known that blood vessels dilate and constrict to maintain homeostatic levels of blood flow and pressure. However it is challenging to quantify the functionality of the heart and vasculature. Cardiovascular quantities including volumes, flows, and pressures (e.g. ventricular volume and pressure  studied here) can be measured experimentally, yet the system properties that underlies specific outputs can typically not be measured directly, e.g. peripheral vascular resistance, vessel compliance, or cardiac contractility. The simple cardiovascular model presented here can provide experimental biologists a framework to estimate these quantities young (); gelman (). Using the rigorous approach presented here, we are able to assess model uncertainty that would be essential for clinical application.
In summary, our stepbystep approach showed that by removing parameter interactions by means of local sensitivity and structured correlation analysis, we can obtain the same subset of practically identifiable model parameters for all three rats. We used nonlinear least squares optimization to estimate this subset of parameters to provide a high quality fit to the data. Finally, we used both frequentist and Bayesian UQ methods to predict confidence, credibility and prediction intervals around the model output. The approached discussed here was demonstrated on a simple cardiovascular model, yet it is applicable to any model with data that can be formulated as a system of ODEs.
Model Formulation and nominal parameter values
The model type used in this study is not new, and even though previous studies, e.g. Smith04 (); Lu01 (); Puelz17 (); Guyton72 () involve more complex interactions (e.g. nonlinear venous compliance, systemic and pulmonary subcircuits, inductors, etc.), the important conclusion is that to obtain reliable parameter estimates it is essential to compute nominal parameters that are physiologically reasonable. This allows one to determine what parameters are identifiable and to assess model uncertainty. The specific model discussed here includes measurements of left ventricular volume and pressure. Other commonly available data include measurements of arterial pressure, central venous pressure, and cardiac output. For most models, including the one presented here, assumptions must be made since not all information is available, e.g. we used literature estimate to determine the unstressed volume fraction. This is a quantity that is important, but as discussed by Spiegel Spiegel16 (), no good experimental method exists for determining the stressedtounstressed ratio. The model presented here was fitted to left ventricular pressure and volume, the advantage being that both arterial and venous pressures and stroke volume can be determined, while aortic pulse pressure had to be estimated. Many other studies, e.g. the study by Williams et al. Williams14 () rely on noninvasive measurements typically only available at the arterial side. To identify reliable parameters this model can only be done by incorporating assumptions about venous pressure and cardiac output. As discussed in the study by Pironet et al. Pironet16 () reliable parameter estimates require measurements of both pressures and left ventricular stroke volume. In general, we show how to determine a subset of identifiable parameters that can be estimated by fitting model output to data.
The model and data analyzed here are pulsatile. To accurately fit model to data in the case of this periodic cardiovascular systems model it is essential to align the phase of the periodic driving function with that of the data. The phase of model predictions is determined by initial conditions. Since for most compartments only either mean and/or minimum or maximum predictions are available it is difficult to set up initial conditions that provide precise alignment of model predictions to data. Mathematically, we want to initialize the model at “steady state,” so that with constant periodic forcing the solution oscillates around a known mean. Depending on timescale of the model, typically dictated by compartment elastance, computations should be done over a long enough time to allow the system to reach a steady state. Once steady state is achieved, we show how to systematically determine the optimal shift of the data best aligning the model predictions with the data. To our knowledge this element has not been addressed in previous studies with respect to cardiovascular models.
Finally, the pulsatile nature of this model is facilitated by periodic forcing of the system combined with valves. This model uses diodes to predict flow into and out of the valves, while this is the most typical formulation Pironet16 (); sun1995mathematical (), the discrete opening and closing of valves make the system of differential equations stiff. While the diodes approach is easier to formulate, it introduces a discontinuity in the system of equations not characteristic of the pressuredependent resistances Williams14 () or valves accounting for inertia Shi11 (); smith2004minimal (); Mynard12 (). Even when accounting for inertia, the system of ODEs arising in pulsatile cardiovascular models is stiff and requires careful attention when solving numerically; for more details on solving stiff differential equations, see, e.g. LeVeque (); Iserles ().
Below we discuss results of each component of our stepbystep approach engaged to estimate identifiable parameters and determine model uncertainty.
Sensitivity analysis and subset selection
Most cardiovascular models, as the one formulated here, are overparameterized compared to available data. While data available for parameter estimation often is dense in time, only a few states are measured, e.g. the model presented here are based on predictions of left ventricular flow and pressure. Therefore, care must be taken both in the model building phase, by not including more compartments than can be justified, and in the parameter estimation phase, only estimating identifiable parameters. The first step in the analysis phase involves sensitivity analysis. There are many methodologies to consider, and care must be taken since the system of equations often is nonlinear and stiff. For this study we elected to calculate sensitivities using a forward difference methodology, that we validated by reducing our ODE solver tolerance and observing that the results converged to stable values. While this approach is easy, more effort could be put into study the system globally, e.g. using Sobol indices or Morris screening ChristianPhD (); Eck16 (). Alternatively, local results could be obtained by solving sensitivity equations in the complex plane martins2003complex ()  however this method requires that model be analytic and the diode valve formulation we employ violates this condition.
Parameters were ranked from most to least sensitive using a twonorm averaging over the part of the data used for model validation (i.e. after steady oscillations have been achieved, when the effects of initial conditions disappear). Given the periodic nature of the model studied here this approach makes sense, yet for other model types, generalized sensitivity Banks10 () may be more appropriate as they also provide information about what part of the data specific parameters influence.
Results of our sensitivity analysis revealed that and were insensitive compared to the other parameters, as a result we chose to keep those fixed at their nominal values. Another use of sensitivities is to detect if specific model components can be eliminated or should be described in more detail. It should be noted that the sensitives we use reflect how an individual parameter influences the proximity of the model solution to the data rather than providing a metric of how it influences the model behavior.
Parameters being sensitive do not equate to being identifiable Ottesen13 (); Miao11 (). Thus it is important to combine sensitivity analysis with identifiability analysis, or subset selection (the term used here) to identify relationships among parameters. However, determining if a specific parameter interaction is identifiable is nontrivial. Pironet et al. Pironet16 () illustrated that elastance and resistance parameters form a structurally identifiable relationship given left ventricular stroke volume and pressure data for every compartment in their model. Their model uses a simplified elastance driving function with no parameters and assumed that stroke volume and all compartment pressures were available. Consequently results from the Pironet model cannot directly be translated to our model. Mahdi et al. mahdi2014structural () provided some guidelines for constructing arbitrarily complex structurally identifiable springdashpot networks, however their results were constrained to linear viscoelastic formulations. Our methodology is motivated by analyzing the practically identifiable components of cardiovascular models to make predictive inference through the use of uncertainty quantification.
Intuitively one would expect that parameter interactions hinder the identifiability of a model. However it is crucial to understand, even though intimately related, that correlation and identifiability are distinct concepts. Correlation refers to the precise structure and relationship between interacting parameters, whereas identifiability refers to mapping between the parameter space and the model output, which ideally would be onetoone. So while minimizing parameter interactions and correlations can be an effective strategy to find an identifiable subset of parameters, it is not an exhaustive or universally applicable strategy for any arbitrary model. We found that approaching the parameter estimation problem from a Bayesian perspective was particularly useful in this regard as it allows one to visualize parameter interactions. It should also be noted that the ideal subset of “identifiable” parameters is highly dependent on the type of data available (e.g. pressure in large systemic arteries rather vs. the left ventricle) and may vary across individuals.
For this study we chose to use structural correlation analysis Ottesen13 (); Miao11 (), which is local in nature and it only provides a first order approximation of the parameter correlations, i.e. it is not able to capture nonlinear parameter interactions. However, with good initial parameter estimates, this methodology has been shown to work ChristianPhD (); Ottesen13 (); Miao11 (). In this study, we leverage this limitation using DRAM to estimate the parameters in a Bayesian framework. By representing the parameters as random variables, we can trace the exact shape of the joint densities and can visualize the parameter interactions, as shown in Figure 7. Unidentifiable parameter subsets typically take more MCMC iterations to converge and often have more correlation in their pairwise density plots. We have illustrated that in Figure 9 (using a slightly altered subset including , a subset that structured correlation analysis method determined correlated). Results demonstrate very mild correlations between and despite the subset (24) being identifiable. Our recommendation is to combine Bayesian parameter estimation techniques together with asymptotic subset selection methods in an iterative process to refine ones understanding of different model parameterizations and potentially find a practically identifiable subset of parameters.
Optimization
Similar to numerous other studies, e.g. Pope09 (), we used a nonlinear least squares methodology to obtain pointestimates for the identifiable parameters. This approach is commonly used for ODE models and data describing the kinetics and/or dynamics of a system hovorka2001parameter (). There are many different styles of optimization algorithms with corresponding benefits and limitations. Gradientbased optimization methods such as LevenbergMarquart are limited in that they do not search the entire the parameter space while minimizing the cost function. They are designed to find the nearest local minimum relative to the initial parameter estimate, rather than the global minimum. Other optimization routines that do explore the entirety of the parameter space could have been employed (e.g. NelderMead, SPQ algorithm, or interior point optimization), but for well conditioned problems LevenbergMarquart is computationally very efficient. As a result we put significant effort into a priori parameterization of the model to yield a well conditioned optimization problem, a topic we have not seen addressed in many cardiovascular models.
Uncertainty Quantification
Having a highquality model fit is rarely the end of one’s modeling efforts. In particular if the model is to be adapted for analysis of clinical data it is essential to assess the limitations of a model’s predictive power. This study used asymptotic and Bayesian methods to determine confidence, prediction and credible intervals. The field of uncertainty quantification is an active area of research schiavazzi2016uncertainty (); chen2013simulation (); sankaran2016uncertainty (); pathmanathan2015uncertainty (); Paun2018 (); Eck16 ()  these references are but a few examples of recent studies applying UQ to models of cardiovascular mechanics. For the utility that UQ promises, the obvious problem that has been neglected with regard effective implementation is that most methods assume that the mapping from the parameter space to the model output is onetoone. This is a nontrivial condition that few models satisfy. Furthermore, very few studies have reported results of identifiability analysis prior to the implementation of the UQ procedure. While one can forcibly propagate the uncertainties of unidentifiable model inputs (parameters and/or data) through a model to obtain uncertainty bounds, the reproducibility and stability of such results are suspect.
Our analysis, utilizing the estimated subset of identifiable parameters (24) has allowed us to use both frequentist and Bayesian uncertainty propagation methods to achieve interpretable and tight uncertainty bounds around our model solutions. As a result of using a practically identifiable subset of parameters, our UQ bounds are reproducible and stable across distinct data sets. Results of analysis showed that asymptotic prediction intervals and the Bayesian prediction intervals are almost identical for pressure data, whereas prediction intervals are slightly wider for the volume data, the latter is likely a result of a higher noise level in the volume data than the pressure data.
Due to its flexibility in representing parameters as random variables, an advantage of the Bayesian framework is that UQ is an intrinsic feature. Assuming the model is identifiable and that the parameter chains has converged, credible and prediction intervals can be obtained by extracting model evaluations from the parameter chains. These uncertainty bounds are not subject to any potentially limiting assumption about their distribution, as their frequentist analogs are. The largest limitation to MCMCtype Bayesian approaches is the sheer computational cost. MCMC methods are robust, yet they require thousands of model evaluations to construct a converged posterior parameter distribution. For sufficiently complex models, MCMC may be an impractical route. As an alternative, sequential Bayesian methods such as particle filtering Kaipio2005 (); LiuWest2001 (); Pitt1999 (); Arnold2013 () or ensemble Kalman filtering Evensen1994 (); Burgers1998 (); Evensen2009 (); Arnold2014 () may reduce computational time by evaluating the model from one data point to the next, as opposed to integrating over the entire data set at once. Another option is to use model emulation within the DRAM evaluation rasmussen2006gaussian (). For the cardiovascular model presented here, we can see in Figure 8 that there is very little appreciable difference between the employed statistical approaches. Thus for future studies using this model and corresponding data sets, frequentist UQ alone may be sufficient.
In conclusion, we developed a model and estimated subject specific nominal parameters using available data and literature. We used local sensitivity analysis and subset selection to determine a set of identifiable parameters that were optimized against left ventricular pressure and volume data from three rats. Subsequently we used asymptotic and Bayesian analysis to estimate prediction and credible intervals, which since done in a close neighborhood of the optimal values were comparable. Given the computational cost associated with Bayesian methods, we propose to identify possible subsets and conduct estimations using both local and Bayesian approaches for a few test animals and then carry out larger population studies using only local methods.
Appendix
Model Equations
The complete system of differential equations describing the rates of change of the compartments of the model analyzed in this study are given as follows:
where
and
Unidentifiable subset
Representing parameters as probability distributions can provide information about the identifiability of a multidimensional parameter distribution. Plots of parameter chains and densities such as those seen in Figures 7 are the fundamental diagnostic tools to asses if the MCMC parameter optimizer has converged. The posterior parameter chain plot should ideally be white noise of the distribution, and the posterior parameter density should have a single clearly defined mean. Unidentifiable parameter distributions often have multimodal distributions that can be observed from plots of the parameter chain and the density. One can interpret a multimodal distribution as a parameter chain that has not converged, or an unidentifiable parameter. Should the first scenario be the case be, one can run the MCMC algorithm for more iterations until the distribution converges to the true posterior. However, if the converged parameter distribution is still multimodal, it indicates that multiple values of a parameter can be used to produce the same model output. By definition, this means that the parameter is unidentifiable.
To illustrate how an unidentifiable subset of parameters can result in multimodal posterior parameter distribution we ran DRAM using a subset which included . Results from this simulation, shown in Figure 9, revealed that it takes approximately 12,000 iterations for the parameter chains to burnin (compared with the identifiable subset that start near the converged distribution), even with this long burnin the chains are not well mixed. Poorly mixed parameter chains are indicative of an illconditioned optimization problem cowles1996markov (). Further note that the parameter densities have lumpy shapes, e.g. parameters and exhibit bimodality. Furthermore, we observe that and have a highly correlated joint density. One may be able to run the MCMC routine for more iterations to get better mixed chains and smooth out the densities, but for an arbitrarily complex model it may not be worth the effort given the raw computational time and resources MCMC routines require. We note that while the modality of posterior parameter distributions can be used to assess the identifiability of a subset of parameters, an individual parameter exhibiting multimodal behavior is not itself unidentifiable. From these results alone, we could not explicitly determine if or are unidentifiable. Additionally, there are situations where a third parameter with a relatively smooth density may in fact be the unidentifiable culprit. The posterior parameter density is a jointdistribution of all the parameters, so any multimodal behavior must be considered in the context of the other parameters in the distribution.
References
References
 (1) L. Formaggia, A. M. Quarteroni, A. Veneziani (Eds.), Cardiovascular mathematics: modeling and simulation of the circulatory system, Springer, 2009.
 (2) A. M. Quarteroni (Ed.), Modeling the Heart and the Circulatory System, Springer, 2015.
 (3) F. N. Van de Vosse, N. Stergiopulos, Pulse wave propagation in the arterial tree, Ann Rev Fluid Mech 43 (2011) 467–499.
 (4) J. T. Ottesen, M. S. Olufsen, J. K. Larsen, Applied Mathematical Models in Human Physiology, SIAM, 2004.
 (5) J. T. Ottesen, V. Novak, M. S. Olufsen, Development of patient specific cardiovascular models predicting dynamics in response to orthostatic stress challenges, in: Mathematical Modeling and Validation in Physiology, Springer, 2013, pp. 177–213.
 (6) P. Blanco, F. RA, A 3d1d0d computational model for the entire cardiovascular system, Mecánica Computacional 24 (2010) 5887–5911.
 (7) I. Kokalari, T. Karaja, M. Guerrisi, Review on lumped parameter method for modeling the blood flow in systemic arteries, J Biomed Sci Eng 6 (2013) 92–99.
 (8) S. Yubing, P. Lawford, R. Hose, Review of zerod and 1d models of blood flow in the cardiovascular system, Biomed Eng Online 10.
 (9) M. L. Neal, J. B. Bassingthwaighte, Subjectspecific model estimation of cardiac output and blood volume during hemorrhage, Cardiovasc Eng 7 (2007) 97–120.
 (10) D. Zinemanas, R. Beyar, S. Sideman, Relating mechanics, blood flow and mass transport in the cardiac muscle, Int J Heat Mass Transf 37 (1994) 191–205.
 (11) M. Olufsen, J. Ottesen, H. Tran, L. Ellwein, L. Lipsitz, N. V, Blood pressure and blood flow variation during postural change from sitting to standing: model development and validation, J Appl Physiol 99 (2005) 1523–1537.
 (12) S. Pope, L. Ellwein, C. Zapata, V. Novak, C. Kelley, M. Olufsen, Estimation and identification of parameters in a lumped cerebrovascular model, Math Biosci Eng 6 (2009) 93–115.
 (13) N. Williams, O. WindWillassen, R. Program, J. Mehlsen, J. Ottesen, M. Olufsen, Patientspecific modelling of headup tilt, Math Med Biol 31 (2014) 365–392.
 (14) J. Revie, D. Stevenson, J. Chase, C. Hann, B. Lambermont, A. Ghuysen, P. Kolh, G. Shaw, S. Heldmann, T. Desaive, Validation of subjectspecific cardiovascular system models from porcine measurements, Comp Meth Prog Biomed 109 (2013) 197–210.
 (15) P. Pacher, T. Nagayama, P. Mukhopadhyay, S. Bátkai, D. A. Kass, Measurement of cardiac function using pressure–volume conductance catheter technique in mice and rats, Nature protocols 3 (2008) 1422.
 (16) E. T. Mackenzie, J. K. Farrar, W. Fitch, D. I. Graham, P. C. Gregory, A. M. Harper, Effects of hemorrhagic hypotension on the cerebral circulation. i. cerebral blood flow and pial arteriolar caliber, Stroke 10 (1979) 711–718.
 (17) A. Mahdi, N. Meshkat, S. Sullivant, Structural identifiability of viscoelastic mechanical systems, PloS one 9 (2014) e86411.
 (18) H. Miao, X. Xia, A. S. Perelson, H. Wu, On identifiability of nonlinear ode models and applications in viral dynamics, SIAM rev 53 (2011) 3–39.
 (19) J. Kirk, M. Saccomani, S. Shroff, A priori identifiability analysis of cardiovascular models, Cadiovasc Eng Technol 4 (2013) 500–512.
 (20) A. Pironet, P. Dauby, J. Chase, P. Docherty, J. Revie, T. Desaive, Structural identifiability analysis of a cardiovascular system model, Med Eng Physics 38 (2016) 433–441.
 (21) L. Ellwein, H. Tran, C. Zapata, V. Novak, M. Olufsen, Sensitivity analysis and model assessment: mathematical models for arterial blood flow and blood pressure, Cardiovasc Eng 8 (2008) 94–108.
 (22) R. Gul, Mathematical modeling and sensitivity analysis of lumpedparameter model of the human cardiovascular system, Ph.D. thesis, The Free University, Berlin, Germany (2016).
 (23) V. Eck, W. Donders, J. Sturdy, J. Feinberg, T. Delhaas, L. Hellevik, W. Huberts, A guide to uncertainty quantification and sensitivity analysis for cardiovascular applications, Num Meth biomed Eng doi: 10.1002/cnm.2755.
 (24) P. Chen, A. Quarteroni, G. Rozza, Simulationbased uncertainty quantification of human arterial network hemodynamics, Int J Numer Meth Biomed Eng 29 (2013) 698–721.
 (25) A. Arnold, C. Battista, D. Bia, Y. Zocalo, R. Armentano, H. Tran, M. Olufsen, Uncertainty quantification in a patientspecific onedimensional arterial network model: EnKFbased inflow estimator, ASME J Ver Val Uncer Quant 2 (2017) 011002.
 (26) V. Eck, J. Feinberg, H. Langtangen, L. Hellevik, Stochastic sensitivity analysis for timing and amplitude of pressure waves in the arterial system, Int J Numer Meth Biomed Eng DOI: 10.1002/cnm.2711.
 (27) V. Eck, J. Feinberg, H. Langtangen, L. Hellevik, Stochastic sensitivity analysis for timing and amplitude of pressure waves in the arterial system, Int J Num Meth Biomed Eng 31.
 (28) V. G. Eck, W. P. Donders, J. Sturdy, J. Feinberg, T. Delhaas, L. R. Hellevik, W. Huberts, A guide to uncertainty quantification and sensitivity analysis for cardiovascular applications, Int J Num Meth Biomed Eng 32.
 (29) M. Paun, Q. MU, M. Colebank, N. Hill, D. Husmeier, MCMC methods for inference in a mathematical model of pulmonary circulation, Statistica Neerlandica, Accepted for publication.
 (30) H. Haario, M. Laine, A. Mira, E. Saksman, Dram: efficient adaptive mcmc, Stat Comput 16 (2006) 339–354.
 (31) J. Baan, T. T. A. Jong, P. L. Kerkhof, R. J. Moene, A. D. Van Dijk, E. T. van Der Velde, J. Koops, Continuous stroke volume and cardiac output from intraventricular dimensions obtained with impedance catheter, Cardiovascular Research 15 (1981) 328–334.
 (32) J. Kjaergaard, C. L. Petersen, A. Kjaer, B. K. Schaadt, J. K. Oh, C. Hassager, Evaluation of right ventricular volume and function by 2d and 3d echocardiography compared to mri, European journal of echocardiography 7 (2006) 430–438.
 (33) H. S. Hwang, M. O. Boluyt, K. Converso, M. W. Russell, B. E. Bleske, Effects of hawthorn on the progression of heart failure in a rat model of aortic constriction, Pharmacotherapy 29 (2009) 639–648.
 (34) B. Cosyns, S. Droogmans, C. Weytjens, T. Lahoutte, G. Van Camp, D. Schoors, P. R. Franken, Effect of streptozotocininduced diabetes on left ventricular function in adult rats: an in vivo pinhole gated spect study, Cardiovasc Diabetol 6 (2007) 30.
 (35) A. I. AlShafei, R. G. Wise, G. Gresham, T. Carpenter, L. Hall, C. L.H. Huang, Magnetic resonance imaging analysis of cardiac cycle events in diabetic rats: the effect of angiotensinconverting enzyme inhibition, J Physiol 538 (2002) 555–572.
 (36) J. Holt, E. Rhode, H. Kines, Ventricular volumes and body weight in mammals, Am J Physiol 215 (1968) 704–715.
 (37) P. Nordbeck, L. Bönhof, K.H. Hiller, S. Voll, P. AriasLoza, L. Seidlmayer, T. Williams, Y.X. Ye, D. Gensler, T. Pelzer, et al., Impact of thoracic surgery on cardiac morphology and function in small animal models of heart disease: a cardiac mri study in rats, PloS one 8 (2013) e68275.
 (38) S. E. Litwin, T. E. Raya, P. G. Anderson, C. M. Litwin, R. Bressler, S. Goldman, Induction of myocardial hypertrophy after coronary ligation in rats decreases ventricular dilatation and improves systolic function., Circulation 84 (1991) 1819–1827.
 (39) S. K. Engle, P. F. Solter, K. M. Credille, C. M. Bull, S. Adams, M. J. Berna, A. E. Schultze, E. C. Rothstein, M. D. Cockman, M. L. Pritt, et al., Detection of left ventricular hypertrophy in rats administered a peroxisome proliferator–activated receptor / dual agonist using natriuretic peptides and imaging, Toxicol Sci 114 (2009) 183–192.
 (40) R. Wise, C.H. Huang, G. Gresham, A. AlShafei, T. Carpenter, L. Hall, Magnetic resonance imaging analysis of left ventricular function in normal and spontaneously hypertensive rats, J Physiol 513 (1998) 873–887.
 (41) M. Nahrendorf, F. Wiesmann, K.H. Hiller, K. Hu, C. Waller, J. Ruff, T. E. Lanz, S. Neubauer, A. Haase, G. Ertl, et al., Serial cinemagnetic resonance imaging of left ventricular remodeling after myocardial infarction in rats, J Magn Res Imag 14 (2001) 547–555.
 (42) C. Vanhove, T. Lahoutte, M. Defrise, A. Bossuyt, P. R. Franken, Reproducibility of left ventricular volume and ejection fraction measurements in rat using pinhole gated spect, Eur J Nucl Med Mol Imag 32 (2005) 211–220.
 (43) M. P. Bal, W. Vries, F. Leij, M. Oosterhout, J. Baan, E. Wall, F. Bel, P. Steendijk, Left ventricular pressure–volume relationships during normal growth and development in the adult rat–studies in 8and 50weekold male wistar rats, Acta Physiol 185 (2005) 181–191.
 (44) M. Nahrendorf, K.H. Hiller, A. Greiser, S. Köhler, T. Neuberger, K. Hu, C. Waller, M. Albrecht, S. Neubauer, A. Haase, et al., Chronic coronary artery stenosis induces impaired function of remote myocardium: Mri and spectroscopy study in rat, Am J Physiol 285 (2003) H2712–H2721.
 (45) T. Radovits, A. Oláh, Á. Lux, B. T. Németh, L. Hidi, E. Birtalan, D. Kellermayer, C. Mátyás, G. Szabó, B. Merkely, Rat model of exerciseinduced cardiac hypertrophy: hemodynamic characterization using left ventricular pressurevolume analysis, Am J Physiol 305 (2013) H124–H134.
 (46) S. KorkmazIcöz, A. Lehner, S. Li, A. Vater, T. Radovits, M. Brune, M. Ruppert, X. Sun, P. Brlecic, M. Zorn, et al., Left ventricular pressurevolume measurements and myocardial gene expression profile in type2 diabetic gotokakizaki rats, Am J Physiol (2016) ajpheart–00956.
 (47) A. Todica, G. Böning, S. Lehner, E. Weidl, P. Cumming, C. Wängler, S. G. Nekolla, M. Schwaiger, P. Bartenstein, R. Schirrmacher, et al., Positron emission tomography in the assessment of left ventricular function in healthy rats: A comparison of four imaging methods, J Nucl Cardiol 20 (2013) 262–274.
 (48) C. A. Carr, D. J. Stuckey, L. Tatton, D. J. Tyler, S. J. Hale, D. Sweeney, J. E. Schneider, E. MartinRendon, G. K. Radda, S. E. Harding, et al., Bone marrowderived stromal cells home to and remain in the infarcted rat heart but fail to improve function: an in vivo cinemri study, Am J Physiol 295 (2008) H533–H542.
 (49) J. R. Jones, J. F. Mata, Z. Yang, B. A. French, J. N. Oshinski, Left ventricular remodeling subsequent to reperfused myocardial infarction: evaluation of a rat model using cardiac magnetic resonance imaging, J Cardiovasc Magn Res 4 (2002) 317–326.
 (50) D. J. Stuckey, C. A. Carr, D. J. Tyler, E. Aasum, K. Clarke, Novel mri method to detect altered left ventricular ejection and filling patterns in rodent models of disease, Magn Res Med 60 (2008) 582–587.
 (51) J.L. Daire, J.P. Jacob, J.N. Hyacinthe, P. Croisille, K. MontetAbou, S. Richter, D. Botsikas, M. LepetitCoiffé, D. Morel, J.P. Vallée, Cine and tagged cardiovascular magnetic resonance imaging in normal rat at 1.5 t: a rest and stress study, J Cardiovasc Magn Res 10 (2008) 48.
 (52) M. Ruppert, S. KorkmazIcöz, S. Li, B. T. Németh, P. Hegedűs, P. Brlecic, C. Mátyás, M. Zorn, B. Merkely, M. Karck, et al., Myocardial reverse remodeling after pressure unloading is associated with maintained cardiac mechanoenergetics in a rat model of left ventricular hypertrophy, Am J Physiol 311 (2016) H592–H603.
 (53) N. Trippodo, Total circulatory capacity in the rat. effects of epinephrine and vasopressin on compliance and unstressed volume., Circ Res 49 (1981) 923–931.
 (54) D. B. Young, Control of cardiac output, Integrated Systems Physiology: From Molecule to Function 2 (2010) 1–97.
 (55) S. Gelman, Venous function and central venous pressurea physiologic story, J Am Soc Anesthesiol 108 (2008) 735–748.
 (56) J. Beneken, B. DeWit, A physical approach to hemodynamic aspects of the human cardiovascular system, Physical bases of circulatory transport: Regulation and exchange (1967) 1–45.
 (57) R. J. Gotwals. Cardiovascular physiology: The windkessel model [online] (20002003).
 (58) R. Klabunde. Cardiovascular physiology concepts [online] (2014).
 (59) G. London, A. C. Simon, Y. Weiss, M. E. Safar, Arterial and venous systems in essential hypertension, Vol. 63, Springer Science & Business Media, 2012.
 (60) I. C. Ipsen, C. Kelley, S. Pope, Rankdeficient nonlinear least squares problems and subset selection, SIAM Journal on Numerical Analysis 49 (2011) 1244–1266.
 (61) H. Miao, X. Xia, A. Perelson, H. Wu, On identifiability of nonlinear ode models and applications in viral dynamics, SIAM Rev 53 (2011) 3–39.
 (62) J. T. Ottesen, J. Mehlsen, M. S. Olufsen, Structural correlation method for model reduction and practical estimation of patient specific parameters illustrated on heart rate regulation, Math Biosci 257 (2014) 50–59.
 (63) C. T. Kelley, Iterative methods for optimization, Vol. 18, Siam, 1999.
 (64) R. C. Smith, Uncertainty quantification: theory, implementation, and applications, Vol. 12, SIAM, 2013.
 (65) S. Smolyak, Quadrature and interpolation formulas for tensor products of certain classes of functions, Sov Phys Dokl 148 (1963) 1042–1045.
 (66) H.J. Bungartz, M. Griebel, Sparse grids, Acta Numerica 13 (2004) 147–269.
 (67) J. Halton, On the efficiency of certain quasirandom sequences of points in evaluating multidimensional integrals, Numerische Mathematik 2 (1960) 84–90.
 (68) S. Joe, F. Kuo, Remark on algorithm 659: Implementing sobol’s quasirandom sequence generator, ACM Trans Mathe Software (TOMS) 29 (2003) 49–57.
 (69) C. Andrieu, J. Thoms, A tutorial on adaptive mcmc, Stat Comput 18 (2008) 343–373.
 (70) A. Mira, On metropolishastings algorithm with delayed rejection, Metron LIX (2001) 231–241.
 (71) H. Haario, E. Saksman, J. Tamminen, An adaptive metropolis algorithm, Bernoulli 7 (2001) 223–242.
 (72) B. Smith, J. Chase, R. Nokes, G. Shaw, G. Wake, Minimal haemodynamic system model including ventricular interaction and valve dynamics, Med Eng Phy 26 (2004) 131–139.
 (73) K. Lu, J. Clark, F. Ghorbel, D. Ware, A. Bidani, A human cardiopulmonary system model applied to the analysis of the valsalva maneuver, Am J Physiol 281 (2001) H2661ÐH2679.
 (74) C. Puelz, S. Acosta, B. Riviere, D. Penny, B. KM, C. Rusin, A computational study of the fontan circulation with fenestration or hepatic vein exclusion, Comp Biol Med 89 (2017) 405–418.
 (75) A. Guyton, T. Coleman, H. Granger, Circulation: Overall regulation, Annu Rev Physiol 34 (1972) 13–46.
 (76) R. Spiegel, Stressed vs. unstressed volume and its relevance to critical care practitioners, Clin Exp Emerg Med 3 (2016) 52–54.
 (77) Y. Sun, B. Sjoberg, P. Ask, D. Loyd, B. Wranne, Mathematical model that characterizes transmitral and pulmonary venous flow velocity patterns, American Journal of PhysiologyHeart and Circulatory Physiology 268 (1995) H476–H489.
 (78) B. W. Smith, J. G. Chase, R. I. Nokes, G. M. Shaw, G. Wake, Minimal haemodynamic system model including ventricular interaction and valve dynamics, Med Eng Phy 26 (2004) 131–139.
 (79) J. Mynard, M. Davidson, D. Penny, J. Smolich, A simple, versatile valve model for use in lumped parameter and onedimensional cardiovascular models, Int J Numer Meth Biomed Eng 28 (2012) 626–641.
 (80) R. J. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations: SteadyState and TimeDependent Problems, SIAM, 2007.
 (81) A. Iserles, A First Course in the Numerical Analysis of Differential Equations, Cambridge University Press, 2008.
 (82) C. Olsen, Modeling heart rate regulation by the baroreflex, Ph.D. thesis, NC State University, Raleigh, NC (2014).
 (83) J. R. Martins, P. Sturdza, J. J. Alonso, The complexstep derivative approximation, ACM Transactions on Mathematical Software (TOMS) 29 (2003) 245–262.
 (84) H. Banks, S. Dediu, S. Ernstberger, F. Kappel, Generalized sensitivities and optimal experimental design, J Inv Illposed Problems 18 (2010) 25–83.
 (85) J. Ottesen, M. Olufsen, A practical approach to parameter estimation applied to model predicting heart rate regulation, J Math Biol 67 (2013) 39–68.
 (86) R. Hovorka, P. Vicini, Parameter estimation, Mod Meth Phys Med (2001) 107–151.
 (87) D. Schiavazzi, G. Arbia, C. Baker, A. M. Hlavacek, T.Y. Hsia, A. Marsden, I. VignonClementel, et al., Uncertainty quantification in virtual surgery hemodynamics predictions for single ventricle palliation, International journal for numerical methods in biomedical engineering 32.
 (88) P. Chen, A. Quarteroni, G. Rozza, Simulationbased uncertainty quantification of human arterial network hemodynamics, International journal for numerical methods in biomedical engineering 29 (2013) 698–721.
 (89) S. Sankaran, H. J. Kim, G. Choi, C. A. Taylor, Uncertainty quantification in coronary blood flow simulations: impact of geometry, boundary conditions and blood viscosity, Journal of biomechanics 49 (2016) 2540–2547.
 (90) P. Pathmanathan, M. S. Shotwell, D. J. Gavaghan, J. M. Cordeiro, R. A. Gray, Uncertainty quantification of fast sodium current steadystate inactivation for multiscale models of cardiac electrophysiology, Progress in biophysics and molecular biology 117 (2015) 4–18.
 (91) J. P. Kaipio, E. Somersalo, Statistical and Computational Inverse Problems, Springer, New York, 2005.
 (92) J. Liu, M. West, Combined parameter and state estimation in simulationbased filtering, in: A. Doucet, N. de Freitas, N. Gordon (Eds.), Sequential Monte Carlo Methods in Practice, Springer, New York, 2001, pp. 197–223.
 (93) M. Pitt, N. Shephard, Filtering via simulation: auxiliary particle filters, J Amer Statist Assoc 94 (1999) 590–599.
 (94) A. Arnold, D. Calvetti, E. Somersalo, Linear multistep methods, particle filtering and sequential Monte Carlo, Inverse Problems 29 (2013) 085007.
 (95) G. Evensen, Sequential data assimilation with a nonlinear quasigeostrophic model using Monte Carlo methods to forecast error statistics, J Geophys Res 99 (1994) 10143–10162.
 (96) G. Burgers, P. J. Leeuwen, G. Evensen, Analysis scheme in the ensemble Kalman filter, Mon Weather Rev 126 (1998) 1719–1724.
 (97) G. Evensen, The ensemble Kalman filter for combined state and parameter estimation, IEEE Control Syst Mag 29 (2009) 83–104.
 (98) A. Arnold, D. Calvetti, E. Somersalo, Parameter estimation for stiff deterministic dynamical systems via ensemble Kalman filter, Inverse Problems 30 (2014) 105008.
 (99) C. E. Rasmussen, C. K. Williams, Gaussian processes for machine learning, Vol. 1, MIT press Cambridge, 2006.
 (100) M. K. Cowles, B. P. Carlin, Markov chain monte carlo convergence diagnostics: a comparative review, Journal of the American Statistical Association 91 (1996) 883–904.