Successful clinical use of patient-specific models for cardiovascular dynamics depends on the reliability of the model output in the presence of input uncertainties. For 1D fluid dynamics models of arterial networks, input uncertainties associated with the model output are related to the specification of vessel and network geometry, parameters within the fluid and wall equations, and parameters used to specify inlet and outlet boundary conditions. This study investigates how uncertainty in the flow profile applied at the inlet boundary of a 1D model affects area and pressure predictions at the center of a single vessel. More specifically, this study develops an iterative scheme based on the ensemble Kalman filter (EnKF) to estimate the temporal inflow profile from a prior distribution of curves. The EnKF-based inflow estimator provides a measure of uncertainty in the size and shape of the estimated inflow, which is propagated through the model to determine the corresponding uncertainty in model predictions of area and pressure. Model predictions are compared to ex vivo area and blood pressure measurements in the ascending aorta, the carotid artery, and the femoral artery of a healthy male Merino sheep. Results discuss dynamics obtained using a linear and a nonlinear viscoelastic wall model.

## Introduction

Numerous studies have used one-dimensional (1D) fluid dynamics models coupled with constitutive arterial wall models to predict wave propagation in arterial networks [1–8]. One-dimensional models are useful for understanding pulse wave propagation in arterial networks [4,7,9–12] and for coupling with 0D and 3D models to generate multiscale representations of larger parts of the cardiovascular system [13–20]. The relatively simple nature of these models allows them to be used in a patient-specific context [21], improving diagnosis of observed hemodynamic disorders, such as pulmonary arterial hypertension [22,23].

Successful use of such patient-specific models in clinical settings depends on the reliability of the model output. Yet, the accuracy of the model output can be traced back to uncertainties in the input parameters [24–26]. Many sources of uncertainty naturally exist in cardiovascular models due to the complexity and variability in the cardiovascular system itself [24,27]. For 1D cardiovascular models, input uncertainties are associated with vessel geometry, specification of equations predicting the fluid dynamics and the arterial wall model, and parameters used to specify the inflow and outflow boundary conditions [1,9], as described below:

- (1)
*Geometry.*For large vessels, geometric properties such as vessel length, diameter and bending, as well as network connectivity, can be obtained from magnetic resonance imaging (MRI) [2,6,22,28] or computed tomography (CT) images [29–31]. Accuracy in these measurements depends on how still and how long the subject lies in the scanner, the choice of image segmentation technique, and the geometric approximations used in reconstructing the vessel and network connectivity from the imaging data [32,33]. This study compares numerical simulations with experimental data where high precision measurements of vessel length (via caliper) and diameter (via ultrasonic crystals) were made directly on single vessels during experimentation, both in vivo before excision and ex vivo after [34,35]. Vessels in this study are assumed to be nontapered (i.e., effects of bends and changes in vessel diameter along vessels are ignored). - (2)
*Fluid dynamics.*In this study, blood is assumed to be Newtonian with constant viscosity and density, and the flow is assumed to be fully developed. As a result, the velocity profile is assumed parabolic [4]. While density and viscosity can be measured fairly accurately in large arteries, the velocity profile is more flat than parabolic [4,36], which gives rise to error in computations. In addition, in smaller vessels viscosity is no longer constant [37]. For the straight vessel studied here, uncertainties associated with these assumptions are minor, particularly since the data are from ex vivo experimental studies in which measurements are not taken until steady pulsating flow has been obtained. - (3)
*One-dimensional model.*This study uses the 1D Navier–Stokes equations for conservation of momentum and volume to relate flow, pressure, and area. To close the system, these are coupled with a constitutive equation relating vessel area and pressure. The latter equation encodes either elastic [6,9,38–40] or viscoelastic [1–3,10,13,41,42] deformation. This 1D approximation is a simplification that does not account for flow disturbances within the network, particularly in the presence of junctions. Moreover, model parameters defining viscosity, flow velocity, vessel stiffness, viscoelastic relaxation, and amplitude are difficult to measure, therefore providing a significant source of uncertainty, even within the single vessel model studied here. - (4)
*Inflow boundary condition.*A boundary condition is specified at the network inlet. This can be done by applying a time-varying flow profile (which may be obtained from MRI measurements) or by using a model that predicts the pumping of the left heart. The flow at the inlet is complex [43], and even if predicted with a 3D heart model [44], it is associated with a significant level of uncertainty [45]. This study focuses on quantifying uncertainty related to the inflow boundary condition. - (5)
*Outflow boundary conditions.*Outflow boundary conditions are typically assigned using either a fractal tree model [6,46–48] or a three-element Windkessel model [2,4,14,49,50] to represent the impact of flow from vessels not included in the network. The measurements of branching angles and ratios in the small vessels used to specify fractal tree models are typically extracted from studies over several species and therefore are not precise. Numerical methods for estimating Windkessel model parameters (two resistors and a capacitor) exist for elastic models but do not hold for viscoelastic models used in this study [8], giving rise to uncertainty in parameter values. - (6)
*Numerical approximation.*In addition to the physical quantities described above, errors associated with the numerical scheme chosen to solve the fluid equations are another source of uncertainty. This study employs a 1D finite element method with specified spatial and time discretization schemes. Previous work studying the convergence of this solver showed that the time discretization significantly affects the accuracy of the numerical solution, while the spatial discretization does not [8]. The numerical error in the time discretization of the solver is used in this study as a means of assessing trust in the model predictions, as described in Sec. 2.3.

Even if measurement techniques for estimating the quantities addressed above were perfect, complete data needed to estimate all of these quantities do not exist. In particular, the lack of concurrent area, pressure, and flow data makes it difficult to validate 1D models. 1D flow and geometry data for large vessels can be extracted from MRI measurements [51], while blood pressure can be measured noninvasively in superficial arteries using Finapres or tonometry methodologies [52,53] or invasively in deeper vessels using a pressure transducer [54,55]. To obtain sufficient resolution, MRI flow measurements must be averaged over several cycles during breath-hold, while blood pressure measurements are instantaneous. Thus, while “complete” validation of 1D models is difficult, acknowledgment of associated uncertainties could significantly improve the use of model-based predictions.

Typical validation of 1D arterial network models includes comparison of predicted and measured flow obtained in networks with geometry specified from MRI measurements. A majority of these studies have been done with models that predict wall deformation using simple linear elastic wall models [6,8,9,14,23,28]. However, it is known that arterial wall deformation (particularly within the aorta) exhibits both nonlinear and viscoelastic characteristics [42]. Some studies have accounted for linear viscoelastic [2,7,13] or nonlinear elastic [31,56] deformation, yet to our knowledge none have accounted for nonlinear viscoelasticity. Accurate model prediction necessitates the use of nonlinear viscoelastic wall models along with matching inflow and outflow boundary conditions and other previously mentioned inputs [2,57]. All of these factors influence wave propagation as well as the type and nature of the reflected waves, thus impacting the shape of the pressure waveform [58].

Most of the studies discussed above are deterministic, i.e., the proposed models fit given waveforms without addressing how changes in input parameters affect output predictions. However, several recent studies [24–26,59–66] have noted the importance of quantifying uncertainty in the model output given variation in the input, which becomes especially important in developing trustworthy models to aid in clinical decision-making. To our knowledge, only a few of these studies have addressed this topic for 1D arterial network models. Chen et al. [24] conducted stochastic simulations within a 1D network to explore the sensitivity of blood flow and pressure to uncertainty in geometric and hemodynamic parameters, including blood density and viscosity. These parameters were randomized using a log-normal distribution, then model predictions of flow and pressure were computed. Sankaran and Marsden [25] took a similar approach, using stochastic collocation to examine uncertainties relating to vessel geometry, boundary conditions, and flow-split in 3D models for abdominal aortic aneurysm, carotid artery bifurcation, and a patient-specific Fontan surgery. Sankaran et al. [59] studied the impact of uncertainties in vessel geometry on patient-specific simulations of the coronary arteries, coupling stochastic collocation with a machine learning-based surrogate to the 3D model. More recently, Sankaran et al. [60] presented a data-driven approach using adaptive stochastic collocation to quantify uncertainties relating to vessel geometry, boundary conditions, and blood viscosity in a patient-specific model of the coronary arteries. This study analyzed the relative impact of uncertainties in minimum lumen diameter, lesion length, boundary resistance, and blood viscosity on blood flow simulations, highlighting the importance of uncertainty in minimum lumen diameter. Eck et al. [26] employed a generalized polynomial chaos (gPC) method computed by stochastic collocation to study the sensitivity of pressure waves to arterial stiffness, while Xiu and Sherwin [61] used gPC to explore uncertainty in pulse wave propagation in a reduced human arterial network model by randomizing a material parameter related to wave propagation speed. Xiao [62] employed a reduced-order unscented Kalman filter to estimate wall stiffness in the carotid artery, as well as Windkessel model resistance and compliance parameters in an idealized aortic bifurcation in a 3D fluid–solid interaction model. Eck et al. [63] recently presented a six-step procedure outlining uncertainty quantification and sensitivity analysis in cardiovascular models, focusing on Monte Carlo and polynomial chaos techniques, while Schiavazzi et al. [64] presented a framework for propagating uncertainties in clinical data to uncertainties in model parameters and model output, using Bayesian parameter estimation techniques, on a patient-specific 3D model simulating virtual stage II single ventricle palliation surgery. The closest related study is that of Brault et al. [65,66], which used projection-based gPC to investigate the sensitivity of 1D model output to uncertainty in the inlet boundary condition. Brault et al. prescribed a periodic inflow model with known form and randomized parameters such as mean and peak flow rate, treating them as independent, uniformly distributed random variables.

The primary objective of this study is to use ex vivo area and pressure data from three single vessels (the ascending aorta, the carotid artery, and the femoral artery) extracted from a healthy male Merino sheep to predict the inflow waveform. In particular, we aim to quantify uncertainty in the shape and magnitude of the estimated inflow waveform, then use this input uncertainty to quantify the corresponding uncertainty in the model output predictions of area and pressure at the center of the vessel for each of these three single vessels of different sizes taken from different locations in the arterial network. This is shown using both a linear and a nonlinear viscoelastic wall model.

The inflow boundary condition is estimated in this study using a Bayesian inverse problems approach [67,68], which provides a natural measure of uncertainty in the inflow estimation. We show how uncertainty in the inflow profile (i.e., the temporal profile prescribed at the vessel inlet boundary) can be propagated through the model to predict uncertainty in the model output at the center of the vessel. More specifically, we develop an iterative scheme based on the ensemble Kalman filter (EnKF) that estimates the inflow from a prior distribution of curves. The EnKF [69,70] has proven to be a useful tool in estimating both unknown model states and parameters in various applications, including weather prediction [71,72] and reservoir modeling [73]. Numerous variations of the original EnKF for state estimation have been derived; in particular, in this study we adapt the augmented EnKF for combined state and parameter estimation derived in Arnold et al. [74].

To our knowledge, this is the first study to propose a methodology for systematically estimating the inflow boundary condition and one of the first studies to address propagation of uncertainty in arterial network models. In our view, prediction of output uncertainty is essential in successfully rendering physiological models patient-specific for use in clinical decision-making. While the current study focuses on inflow estimation, this is just one example that provides a foundation for how such methodologies can be applied in this setting.

In addition to estimating the inflow boundary condition, this study also compares area and pressure predictions obtained using linear and nonlinear viscoelastic models. A previous study showed that including nonlinear elastic deformation is essential within the aorta, while linear elastic deformation was adequate to estimate wall deformation in smaller vessels including the carotid and femoral arteries [42]. However, that study was done without accounting for the fluid dynamics considered here. Results of this study show that the proposed methodology is able to estimate inflow boundary conditions and predict the corresponding output uncertainties for vessels of different size and that accounting for nonlinear viscoelastic deformation improves the estimation.

## Materials and Methods

### Experimental Data.

Data for this study include blood pressure and cross-sectional area measurements from the ascending aorta (AA), the femoral artery (FA), and the carotid artery (CA) of a healthy male Merino sheep, aged 18–24 months with an approximate weight of 32 kg. The excised vessel segments were mounted in the organ chamber of the mock circulation (shown in Fig. 1) and stretched to their in vivo length, which was measured using a high-precision caliper. The vessels were subjected to physiological hemodynamic conditions induced by a Jarvik heart pump. The external arterial diameter of each vessel was measured using sonomicrometry, and blood pressure was measured using a solid-state microtransducer. More details of the experimental protocol can be found in Armentano et al. [75] and Valdez-Jasso et al. [76]. In this work, we analyze the measured area and pressure data over one cardiac cycle in each of the AA, FA, and CA, treated as single vessels of different sizes taken from different locations in the network. The vessels, along with their corresponding area and pressure data, are shown in their network location in Fig. 2.

### Fluid Dynamics Model.

where the blood density $\rho =1.06$ g ml^{−1} and viscosity $\mu =$ 0.049 g s^{−1} cm^{−1} are assumed constant, with kinematic viscosity $\nu =\mu /\rho $. The area of the vessel is denoted by $A=A(x,t)$ (measured in cm^{2}), the pressure by $p=p(x,t)$ (mmHg), and the flow by $q=q(x,t)$ (ml s^{−1}).

where $Am$ denotes the maximal cross-sectional area, $\alpha $ the characteristic pressure at which the vessel begins to saturate, and $k$ an exponent determining the steepness of the sigmoid curve [42].

with two resistors $R1$ and $R2$ (dynes s cm^{−5}) representing the characteristic impedance and total peripheral resistance, respectively, and a capacitor $C$ (cm^{5} dyne^{−1}) denoting arterial compliance.

### The EnKF-Based Inflow Estimator.

This study develops an EnKF-based iterative scheme based on the work by Arnold et al. [74] that, given discrete observations of area and pressure at the center of the vessel, estimates a time-varying flow profile from a prior distribution of curves; this scheme is hereby referred to as the EnKF-based inflow estimator. A brief review of the idea behind Bayesian filtering and the general augmented EnKF algorithm for state and parameter estimation is given in the Appendix. We adapt this methodology for the current application to perform inflow estimation in a 1D vessel.

and we treat the inflow curve $q=q(0,t)$ as an unknown, time-varying parameter to be estimated at $T+1$ discrete time points $tj$, $j=0,1,\u2026,T$. Using the fluid dynamics model (1)–(2), we take a given inflow curve and compute the corresponding area and pressure curves at $xmid$ by employing a stabilized space-time finite element method based on the discontinuous Galerkin method, which we refer to in this paper as the “fluid solver” [3]. The spatial discretization of the fluid solver uses continuous linear polynomials, while the temporal discretization is defined by piecewise constant functions. Previous work studied convergence of the fluid solver using different wall models, varying the number of elements, the number of time steps per period, and the total number of periods [8]. Based on this work, we use 12 elements and a minimum of 400 time steps per period, ensuring convergence of the solver with both the linear and nonlinear wall models. We run the solver for 12 periods to guarantee that the effects of initial conditions wear off, requiring that the last and next-to-last periods differ by no more than a prescribed tolerance.

where $uj=u(xmid,tj)$ and $ej$ is some observation error. For simplicity, we can always discretize the inflow so that its discretized time points match the time points of the data.

#### Generating a Prior.

Let $q$ be a vector denoting the time-discretized inflow curve that we aim to estimate such that $q=[q(0,t0),q(0,t1),\u2026,q(0,tT)]\u2208\mathbb{R}T+1$. Since the fluid solver requires a value of $q$ at each discretized time point $tj$, $j=0,1,\u2026,T$, we begin our EnKF-based scheme to estimate $q$ by assigning each entry $qj=q(0,tj)$ an initial probability distribution with mean $\mu j$ and variance $\sigma j2$, as illustrated in Fig. 3. In this study, we let $qj$ be normally distributed, although the methodology developed here is not restricted by this choice. For each $j$, the mean $\mu j$ is chosen to be a point $q\u0302j=q\u0302(0,tj)$ from a curve $q\u0302$ estimated *a priori*, and the standard deviation $\sigma j$ is selected based on our level of uncertainty in the shape and magnitude of the curve at that point, as described in the bullets below.

The choice of mean curve $q\u0302$ should encompass any

*a priori*knowledge of the data that may help to get a sense of the possible shape of the inflow curve. The $q\u0302$ used in this work was obtained by the following procedure: Originally, measured blood pressure was prescribed at the inlet, and a corresponding flow was predicted using the fluid solver. However, as discussed by Anliker et al. [78], this flow had large negative components following systole. A small negative flow can be expected physiologically, reflecting a reversed flow into the heart upon valve closure. To obtain a positive flow, the Windkessel model parameters were adjusted, resulting in a cardiac output exceeding that known for sheep. Subsequently, the flow profile was scaled to ensure a mean cardiac output of 4.2 l per minute (common for sheep); for more details, see Battista et al. [8].The standard deviation $\sigma j$ can be assigned differently at each $j$ based on some

*a priori*information regarding the certainty in the shape of the curve or set to be the same $\sigma j=\sigma $ for some $\sigma $. For simplicity, in our simulations we assigned the same $\sigma $ at each $j$, but we used a different $\sigma $ for each of the three single vessels to reflect its respective scale.

which are then concatenated into a matrix $S=[s0,s1,\u2026,sT]\u2208\mathbb{R}N\xd7(T+1)$. The rows of the matrix $S$ can be used to form our initial ensemble of $N$ discretized inflow profiles. Since connecting the sample points randomly may result in nonphysical flow profiles, we sort the columns of $S$ in ascending order to generate flow profiles that are physiologically plausible. The rows of the sorted matrix $S$ therefore comprise our initial ensemble of $N$ discretized inflow profiles $qn$, $n=1,\u2026,N$, with mean near $q\u0302$. To enforce regularity and smoothness in the inflow profiles when solving the partial differential equations, each inflow curve is interpolated using a cubic spline so that the time discretization of the flow matches that needed to ensure convergence of the fluid solver.

#### Implementation.

After forming the initial ensemble of inflow curves, the EnKF-based inflow estimator is applied, comprising the following two-step procedure:

- (1)
*Prediction Step*: For each ensemble member $n$, where $n=1,\u2026,N$:- (a)Use the fluid solver to compute the model states $un=[u0n,u1n,\u2026,uTn]$, wherethat correspond to the inflow profile $qn=[q0n,q1n,\u2026,qTn]\u2208\mathbb{R}T+1$ for two different time discretizations: one coarser (but converged) for prediction, denoted by $(un)coarse$; the other finer and converged for error control, denoted by $(un)fine$.$un=[Anpn]=[A0n,A1n,\u2026,ATnp0n,p1n,\u2026,pTn]\u2208\mathbb{R}2\xd7(T+1)$(10)
- (b)Assuming that the model error is related to the time discretization error in the fluid solver, assign the standard deviation $\gamma jn$ of the model innovation term for each ensemble member to be a factor of the difference between the two solutions computed in (a):where $(ujn)coarse$ and $(ujn)fine$ are the model state vectors at time $tj$ computed with the coarse and fine time discretization, respectively, and $\tau >1$ is a multiplicative factor taken here to be $\tau =1.2$. Hence $\gamma jn$ is proportional to the error in the solver difference at each time $tj$ for each ensemble member $n$, allowing for individualized innovation [74,79].$\gamma jn=\tau |(ujn)coarse\u2212(ujn)fine|,j=0,1,\u2026,T$(11)
- (c)
- (d)Combine prediction states and inflow profiles to form an augmented prediction ensemble ${(zn)pred}n=1N={[(z0n)pred,(z1n)pred,\u2026,(zTn)pred]}n=1N$, whereand compute its mean using ensemble statistics$(zn)pred=[(un)predqn]=[(An)pred(pn)predqn]\u2208\mathbb{R}3\xd7(T+1)$(13)$(z\xaf)pred=1N\u2211n=1N(zn)pred\u2208\mathbb{R}3\xd7(T+1)$(14)

- (2)
*Observation update*: For each time point $tj$, $j=0,1,\u2026,T$, sequentially:- (a)
- (b)
- (c)
- (d)
- (e)

The two steps above comprise a single iteration of the algorithm, which may be run successively to ensure convergence of the filter. The algorithm results in a posterior ensemble of inflow profiles, yielding a mean inflow profile (referred to as the EnKF estimate) as well as an estimate of the variance of the ensemble at each discretized time point $tj$, $j=0,1,\u2026,T$. The variance of the resulting ensemble at each discretized time point provides a measure of uncertainty in the estimated inflow at that time point since the variance indicates how wide or narrow the posterior ensemble distribution is around the mean value of the estimated inflow at that point. Therefore, in this work, we use the $\xb12$ standard deviation curves around the estimated mean inflow profile to represent our uncertainty in the inflow estimation.

The resulting mean inflow profile and $\xb12$ standard deviation uncertainty curves are propagated through the model to generate model output predictions of area and pressure at the center of the vessel with corresponding predictions of uncertainty. The model output predictions of uncertainty obtained by propagating the $\xb12$ standard deviation inflow uncertainty curves through the model are interpreted as representing a measure of uncertainty in the model output predictions of area and pressure obtained using the mean estimated inflow curve, thereby connecting the uncertainty in the estimated inflow with uncertainty in the resulting model output.

## Results

The proposed EnKF-based inflow estimator was applied to estimate the unknown temporal flow profile prescribed at inlet boundary given area and pressure data measured at the center of a straight vessel for each of three single vessels (AA, CA, and FA, shown in Fig. 2). Area and pressure time series data for each vessel were provided at 107 discrete time points equispaced over the length of one cardiac cycle. The initial ensemble of inflow profiles for each simulation was formed as described in the methods section on generating prior curves (depicted in Fig. 3). For these simulations, geometry was assumed known and outflow boundary dynamics were determined using a three-element Windkessel model with known parameters.

The method was analyzed on the experimental and simulated AA data *a priori* to determine the effects of varying the number of ensemble members, the convergence of the inflow estimator after multiple iterations of the filtering algorithm, and the filter’s ability to estimate a flow that was able to predict measured area and pressure at the center of the vessel. Based on these results, for the following experiments an ensemble of size $N=100$ was used and the algorithm was run for one iteration. The EnKF computations were performed in matlab, while the fluid equations were solved in C++. Communication between the two platforms was done using text files. Results were obtained running the algorithm sequentially on a Mac OS X desktop computer with a 3.2 GHz Intel Core i3 processor and 8 GB RAM. Simulations took an average of 4.5 h each to complete.

We first applied the EnKF-based inflow estimator to data from the ascending aorta using the linear Kelvin viscoelastic wall model (3) with known parameters to predict radial deformation. Parameters for the wall model were chosen to fit the data when inputting the *a priori* flow curve $q\u0302$ as the inlet boundary condition; these parameter values are listed in Table 1 along with the Windkessel model parameters for the outflow boundary condition. As described in the methods section, the two-step procedure of the inflow estimation algorithm was performed using the innovation model with standard deviation assigned in Eq. (11), with the coarse time discretization set to 428 time steps and the fine set to 856 time steps. The resulting estimated inflow waveform and $\xb1$ 2 standard deviation uncertainty curves are shown in Fig. 4(a), while the corresponding model output predictions and output uncertainty curves are shown in Figs. 4(b)–4(d). While a majority of the area and pressure data is captured within the predicted uncertainty bounds, there is clear underestimation of the area time series shown in Fig. 4(b) toward the beginning of the cardiac cycle and of the pressure time series in Fig. 4(c) toward the end. The portions of the data not captured by the model predictions are more easily seen in the area versus pressure plot in Fig. 4(d).

Motivated by these results, we attempted to capture more of the data within the predicted uncertainty bounds by increasing the standard deviation of the model innovation term defined in Eq. (11) to include a multiple of the difference between the data and the baseline EnKF mean prediction in both area and pressure. This essentially directs the filter in a systematic manner to trust the data more than it trusts the model, which results in increased uncertainty in the estimated inflow and corresponding model predictions shown in Figs. 4(e)–4(h). The plots in Figs. 4(f) and 4(g) better capture the time series data for both area and pressure within the wider predicted uncertainty bounds.

In both of the aforementioned cases, the nonlinearity in the AA pressure-area data (showing decreased deformation with increased pressure) cannot be fully captured by the linear wall model. The area versus pressure plots in Figs. 4(d) and 4(h) make this especially clear. Inspired by results in Valdez-Jasso et al. [42], we extended the wall model to include nonlinear viscoelasticity. In particular, we used the nonlinear sigmoid viscoelastic wall model (4) with known parameters listed in Table 2, keeping the original implementation of the innovation standard deviation in Eq. (11). Table 2 also lists the Windkessel parameters used for the outflow. Results are shown in Figs. 4(i)–4(l), where the nonlinearity in the AA data is much better accounted for in the model predictions (especially visible in Fig. 4(l)). Tests using the sigmoid model (4) and increasing the innovation standard deviation in Eq. (11) did not show significant improvements. However, comparing the plots in Figs. 4(i)–4(l) to Figs 4(a)–4(d), it can be seen that using the sigmoid wall model results in an inflow estimate and corresponding model output predictions of area and pressure with wider uncertainty bounds about the mean curve than when using the Kelvin model. Simulations to confirm this result using synthetic data generated with a known inflow profile are included in the Appendix.

The aforementioned results for the ascending aorta suggest that, while resulting in wider uncertainty bounds, including nonlinear viscoelasticity in the wall model significantly improves predictions of area and pressure by better capturing nonlinearity present in the data. With this in mind, we applied the EnKF-based inflow estimator using the sigmoid wall model (4) with parameters in Table 2 to estimate the inflow profiles for the carotid and femoral arteries. We chose these single vessels as representatives to demonstrate the effectiveness of the proposed algorithm in estimating the inflow and corresponding uncertainty for vessels of different sizes from different locations within the arterial network. Resulting model predictions and uncertainty estimates for area and pressure in each vessel are shown in Fig. 5, which also includes the AA results from Figs. 4(i)–4(l) for completeness. As an additional measure of overall system performance, we used the estimated inflow profiles and $\xb12$ standard deviation uncertainty bounds to predict the flow at the center of each vessel, also shown in Fig. 5 with corresponding uncertainty estimates.

## Discussion

The success of patient-specific models in clinical decision-making relies heavily on the accuracy of the model predictions. Uncertainty quantification plays a vital role in determining how the model output is affected by input uncertainties. For 1D fluid dynamics arterial models, as noted in the introduction, input uncertainties can be found in vessel geometry, fluid model assumptions, parameters specifying the arterial wall model, and parameters relating to the inflow and outflow boundary conditions.

The methodology developed in this study presents a promising technique for quantifying uncertainty in the inflow boundary condition along with the uncertainty in model predictions of area and pressure. More specifically, this study developed an EnKF-based iterative filtering scheme to systematically estimate the inflow boundary condition for a single, nontapered vessel given measurements of area and pressure at the center of the vessel, providing not only an estimate of the size and shape of the inflow profile but also a natural measure of uncertainty in the output.

### Inflow Estimation.

To our knowledge, this is the first study to systematically estimate the inflow boundary condition. Brault et al. [65] acknowledged inflow uncertainty by assigning a periodic model with known form for the inflow and randomizing parameters including mean and peak flow rate, but that study relied on specification of the inflow by a parameterized equation. Brault et al. also relied on the assumption that randomized parameters could be adequately characterized using independent uniformly distributed random variables, and the study did not compare model output results to experimental data. Our methodology, on the other hand, does not require any explicit inflow model, offering more freedom in estimating the shape best suited to predict given area and pressure data, provided some *a priori* estimate of the inflow profile is available (see Fig. 3). Our use of the EnKF also allows for unique posterior distributions to be determined at each time-discretized point of the estimated inflow so that output predictions do not rely on the choice of an initial distribution.

Moreover, our proposed scheme provides a measure of uncertainty arising naturally in the inflow estimation that can be used to assess uncertainty in the model output predictions of area and pressure. The measure of uncertainty obtained in our inflow estimation is an inherent feature of the EnKF and stems from the use of Bayesian filtering to solve the inverse problem of estimating unknown model inputs. Other methods, such as gPC and stochastic collocation, prescribe uncertainty in the input by assigning distributions (with designated means and variances) to unknown input parameters, and then assess the sensitivity of model output with respect to these input uncertainties—input estimation is not explicitly performed. Chen et al. [24], Sankaran and Marsden [25], and Eck et al. [26] took this type of approach in quantifying uncertainties related to various input parameters, such as vessel geometry and physical parameters including blood density and viscosity. The Bayesian approach, on the other hand, assigns a prior distribution to each unknown input parameter, then estimates the input parameters themselves to find posterior input parameter distributions with mean and variance estimates which are then propagated through the model to assess corresponding output uncertainty.

When generating the prior ensemble of inflow profiles (see Sec. 2.3.1), the samples at each time-discretized point are connected in a structured fashion to produce an ensemble of physically plausible inflow profiles. The effects of correlation structure in the prior generation of inflow profiles on the resulting inflow estimates may be analyzed in future work. Further, Eq. (18) uses correlation between the area, pressure, and inflow at time $j$ to update the inflow at time $j$ but not at other time points along the inflow profile. An alternative approach would be to update all $T+1$ points in the inflow profile at each observation time, by modifying the algorithm to account for correlation between all of the points in the inflow profile and the model predictions of area and pressure at that time. However, this approach would require that the model be rerun after each observation update in order to recompute the area and pressure predictions, which would add significantly to the computing time and computational cost of running the algorithm. Future work may include implementing this type of approach and comparing the inflow estimation results with those obtained in this study.

Results using experimental area and pressure data [76] over one cardiac cycle taken from three excised vessels (AA, CA, and FA) in different locations along the arterial network in a healthy sheep show that the proposed EnKF-based inflow estimator is robust in finding an inflow profile that well matches the data starting from an ensemble of flow profiles with varying magnitudes and shapes. The experimental studies considered in this work were done in single, nontapered vessels, but previous work has used these data to show that measurements from individual vessel segments reflect dynamics observed within a network [8]. Future work may include adapting the EnKF-based inflow estimator developed in this study to accommodate use of tapered vessels as well as multiple vessels in network simulations.

Other future work includes optimizing the code to run more efficiently and at a lower computational cost, perhaps by ridding the matlab interface and implementing the filtering portion of the algorithm in C++. Moreover, the embarrassingly parallel nature of Bayesian filtering algorithms during state propagation (since each ensemble member is independently propagated) suggests a very natural way to parallelize the algorithm. Finally, while the EnKF is a very useful and efficient tool in these types of problems, the methodology developed in this work is not restricted by its use, and other Bayesian nonlinear filtering techniques, such as particle filtering, could be used instead if desired. Particle filtering may be particularly well suited for problems in which the resulting parameter distributions are thought to be non-Gaussian.

### Additional Sources of Uncertainty.

While we focused only on inflow estimation in this study, assuming fixed values for the other typically uncertain model inputs listed in the introduction, our novel use of the EnKF in this setting can be extended to include estimation of additional model inputs, such as vessel geometry and outflow boundary condition parameters—the inflow estimation shown in this work is simply one example of how the proposed methodology can be used to quantify uncertainty in patient-specific cardiovascular models. Kalman-type Bayesian filtering was previously used by Xiao [62] to estimate input parameters including arterial wall stiffness in a 3D fluids model, yet its use in systematically estimating a time-dependent boundary condition for 1D arterial network models is new to this work. Since the model predictions of area and pressure are sensitive to the choice of wall model parameters and outflow boundary Windkessel model parameters listed in Tables 1 and 2, future work aims to optimize these parameters using Bayesian filtering methodology prior to or simultaneously while running the EnKF-based inflow estimator.

Further, the effectiveness of the EnKF-based inflow estimator was demonstrated in this study using patient-specific ex vivo area and pressure measurements over one cardiac cycle obtained from three single vessels (AA, CA, and FA) of different sizes and different arterial network locations in a healthy sheep. However, additional sources of uncertainty lie in the data themselves: there is variation, e.g., among different patients and across cardiac cycles, as well as measurement errors (both device and human). Differences in the ex vivo area and pressure data among the healthy sheep considered in this study can be seen in Battista et al. [8], where plots show the variation among different sheep for different single vessels at various locations in the arterial network. While data from one cardiac cycle from one sheep were analyzed in this study, the patient-specific capabilities of the methodology developed in this work could allow for a systematic comparison of inflow profiles estimated at different vessel locations using data over multiple cardiac cycles from different sheep.

### Modeling Aspects.

Most 1D cardiovascular studies use flow data to discuss effects of viscoelasticity, although its expression relates area and pressure. This makes it difficult to validate predicted pressure-area dynamics. The experimental data [76] used in this study, however, contain measurements of both area and pressure, allowing us to directly assess the relation between the two quantities. Our simulations suggest that the nonlinearity in the pressure-area data is better captured using a nonlinear sigmoid wall model with assumed known parameters than a linear Kelvin model.

A previous study by Valdez-Jasso et al. [42] used the same data to compare the linear Kelvin model with the nonlinear sigmoid model as well as the nonlinear arctangent model proposed by Langewouters et al. [80]. Valdez-Jasso et al. concluded that the sigmoid model improved data prediction over the arctangent model for larger arteries like the aorta, while the Kelvin model was preferable for smaller, stiffer vessels like the carotid artery. In particular, the arctangent model was not able to predict unstressed vessel radius. However, that study neglected to consider the fluid dynamics, which may explain why our results show that use of the nonlinear sigmoid model improved prediction of pressure-area dynamics for all three vessels (AA, CA, and FA) regardless of size and network location. We did not consider the arctangent model here.

Although data from only one sheep were considered in this study to demonstrate the effectiveness of the EnKF-based inflow estimator, the patient-specific benefits of this methodology lie in its ability to analyze and quantify uncertainty for individual patients and then make comparisons. This ideology may be especially useful in analyzing patients with diseases such as hypertension and noting any differences that the estimation process may distinguish between patient types. This kind of analysis could have clinical implications if parameters relating to different patient groups (e.g., healthy control versus hypertensive) were distinguishable enough to potentially act has markers for disease.

## Conclusions

This study developed a new iterative scheme based on the EnKF to systematically estimate the inflow boundary condition for a single, nontapered vessel given measurements of area and pressure at the center of the vessel. This estimation scheme provides not only an estimate of the size and shape of the inflow profile but also a natural measure of uncertainty in the model output, offering a promising technique for quantifying uncertainty in patient-specific 1D arterial network models. Results also suggest the use of nonlinear viscoelasticity in the arterial wall model in order to better capture nonlinearities in the data.

## Acknowledgment

This work was supported by Grant No. NSF/DMS-1246991 (Research Training Group in Mathematical Biology at NC State), NSF/DMS-1122424, NSF/DMS-1022688, and Virtual Physiological Rat Project NIH/NIGMS 5P50GM094503 subaward to NC State.

### Appendix A: Bayesian Filtering and the Augmented EnKF

Assume a differential equation model involving both model states $x=x(t)\u2208\mathbb{R}d$ and parameters $\theta \u2208\mathbb{R}k$, whose values may be uncertain or completely unknown. Further assume discrete, noisy observations $y\u2113\u2208\mathbb{R}m$, $\u2113=0,1,\u2026,T$, of some model states. Given these observations, the inverse problem is to estimate the states and/or unknown model parameters. Inverse problems of this type are generally ill-posed [81–83].

This scheme can be viewed as a predictor–corrector method: the first step, known as the prediction step, uses the state evolution equation to make a prediction of the state propagation using only the current data, while the second step, the observation update, corrects that prediction by taking the new data into account.

The posterior mean vectors and covariance matrices for both the states and parameters are computed using the posterior ensemble statistics.

in terms of the projection matrix $G$ from the extended observation model (A9).

computed using the higher order method error control method [79]. Here, $xj+1n$ and $x\u0302j+1n$ are the model solutions using a lower order solver and a higher order solver, respectively, and $\tau >1$ is a safeguard factor to account for the omission of higher order terms. This method provides a systematic approach for estimating the innovation variance in the filter, offering a way to assign the model innovation in Eq. (A5) that is specific to each ensemble member at each time step, while avoiding unnecessary cost in manually choosing a fixed covariance matrix for the innovation term *a priori*.

In estimating the inflow boundary condition in this study, we set $\theta j=q(tj)$ and treat each time-discretized point along the inflow profile as an unknown model parameter. Each time-discretized point in the inflow profile is then updated using its correlation with the corresponding model states (area and pressure) at that time via a modified version of the augmented EnKF, as described in Sec. 2.3. We apply the method of Arnold et al. [74,79] in estimating the variance of the model innovation term, using the numerical error in the time discretization of the fluid solver to set the standard deviation in Eq. (A15), as detailed in step 1(b) of the algorithm presented in Sec. 2.3.2.

Some works in the literature have cited issues using the augmented EnKF for parameter estimation in high-dimensional systems with strong nonlinearities, e.g., see Refs. [85] and [86], or multiplicative combinations of parameters [87]. A suggested alternative approach is to use a dual EnKF scheme involving two interacting filters: one to update the states, and one to update the parameters. Both heuristic [85] and Bayesian consistent [88] versions of the dual EnKF have been proposed. While the 1D model used in this study is nonlinear, particularly with the inclusion of the sigmoid wall model (4), we have not had such issues estimating the inflow profile parameters using the augmented EnKF-type scheme presented in this work. However, it would be possible to implement a dual state-parameter estimation scheme in this setting, which may prove useful if additional parameters (e.g., wall model parameters) were to be estimated along with the inflow profile.

### Appendix B: Wall Model Comparison Using Synthetic Data

The results in Sec. 3, Fig. 4, show that using the nonlinear sigmoid wall model (4) in the fluid dynamics model (1)–(2) provides more accurate model predictions of the pressure-area data than using the linear Kelvin wall model (3), since the sigmoid model better captures the nonlinearity in the data. However, using the sigmoid wall model results in slightly wider uncertainty bounds about the estimated inflow than using the Kelvin wall model, which translates into wider uncertainty bounds about the corresponding model output predictions of area and pressure at the center of the vessel.

Motivated by this result, we performed simulations with synthetic data generated using a known inflow profile to explore the effects that the choice of wall model has on the EnKF-based inflow estimation and corresponding model output predictions, in terms of both accuracy and uncertainty. Taking the flow curve $q\u0302$ described in Sec. 2.3.1 as the known inflow profile, we generated two sets of synthetic area and pressure data: one using the linear Kelvin wall model (which we refer to as the “Kelvin data”) and parameters listed in Table 1, the other using the nonlinear sigmoid wall model (the “sigmoid data”) and parameters for the ascending aorta listed in Table 2. Each set consists of area and pressure time series data at 201 discrete time points, equispaced over the length of one cardiac cycle, with added Gaussian noise.

We then employed the EnKF-based inflow estimator using the different wall models to estimate the inflow profile relating to both the Kelvin and sigmoid data. Fig. 6 shows the results for each of the four combinations: Kelvin wall model on Kelvin data, sigmoid wall model on Kelvin data, Kelvin wall model on sigmoid data, and sigmoid wall model on sigmoid data. To numerically quantify the accuracy in the EnKF-based mean estimated inflow and corresponding model output predictions of area and pressure compared to the true inflow and corresponding synthetic area and pressure data, we compute the coefficient of determination ($R2$) values using the formula

where $SSres$ denotes the residual sum of squares between the estimated and true inflow (or model predictions and synthetic data) and $SStot$ denotes the total sum of squares proportional to the variance in the true inflow (or synthetic data). The ratio $SSres/SStot$ is known as the fraction of unexplained variance: the smaller this fraction is (i.e., the closer to one the $R2$ value is), the smaller the unexplained variance due to error in the model estimate. To numerically quantify the uncertainty, we compute the area between the upper and lower uncertainty bounds ($\xb12$ standard deviation curves about the mean) on the EnKF-based estimated inflow and the corresponding upper and lower uncertainty bounds on the model output predictions. Table 3 lists the $R2$ values and uncertainty estimates for each of the four cases.

In terms of the inflow estimation, the EnKF-based algorithm is able to estimate the true inflow profile fairly well in all four cases, especially when using the Kelvin wall model on the Kelvin data (shown in Fig. 6(a)) and the sigmoid wall model on the sigmoid data (Fig. 6(m)). This is also evident in the inflow $R2$ values listed in Table 3. Although slightly larger when using the sigmoid wall model, the $\xb12$ standard deviation uncertainty bounds are not significantly different in any of the four cases, as is also demonstrated with the inflow uncertainty estimates in Table 3.

Regarding the model output predictions, Figs. 6(b) and 6(c) show that applying the Kelvin wall model on the Kelvin data results in very accurate mean predictions of both area and pressure with tight uncertainty bounds. Use of the Kelvin wall model on the sigmoid data results in less accurate mean predictions but still has tight uncertainty bounds (Figs. 6(j) and 6(k)), while failing to capture the nonlinearities in the data (Fig. 6(l)). On the other hand, Figs. 6(n) and 6(o) show that applying the sigmoid wall model on the sigmoid data gives accurate mean predictions with wider uncertainty bounds than when using the Kelvin wall model. While use of the sigmoid wall model on the Kelvin data does not provide as accurate of mean predictions (Figs. 6(f) and 6(g)) as using the Kelvin wall model, the wider uncertainty bounds on the predictions are better able to capture the data than the tight bounds produced using the Kelvin model. The uncertainty estimates for area and pressure in Table 3 show that the Kelvin wall model has tighter bounds than the sigmoid model for both data sets.

These simulations demonstrate that while use of the Kelvin wall model results in tighter uncertainty bounds on the model predictions for both data sets, it fails to capture the nonlinearities in the sigmoid data. However, while use of the sigmoid wall model may result in better predictions of nonlinear data, the corresponding uncertainty bounds are wider.