D E Schiavazzi

G Arbia

C Baker

A M Hlavacek

T Y Hsia

MAE Alison L Marsden
email: amarsden@eng.ucsd.edu

I E Vignon-Clementel

Uncertainty quantification in virtual surgery hemodynamics predictions for single ventricle palliation

Keywords: uncertainty quantification, Fontan palliation surgery, single-ventricle congenital heart disease, sparse grids, inverse Bayesian parameter estimation, virtual surgery

The adoption of simulation tools to predict surgical outcomes is increasingly leading to questions about the variability of these predictions in the presence of uncertainty associated with the input clinical data. In the present study we propose a methodology for full propagation of uncertainty from clinical data to model results that, unlike deterministic simulation, enables estimation of the confidence associated with model predictions. We illustrate this problem in a virtual Stage II single-ventricle palliation surgery example. First, probability density functions (PDFs) of right pulmonary artery (PA) flow split ratio and average pulmonary pressure are determined from clinical measurements, complemented by literature data. Starting from a 0D semi-empirical approximation, Bayesian parameter estimation is used to find the distributions of boundary conditions that produce the expected PA flow split and average pressure PDFs as pre-operative model results. To reduce computational cost, this inverse problem is solved using a Kriging approximant. Second, uncertainties in the boundary conditions are propagated to simulation predictions. Sparse grid stochastic collocation is employed to statistically characterize model predictions of post-operative hemodynamics in models with and without PA stenosis. The results quantify the statistical variability in virtual surgery predictions, allowing for placement of confidence intervals on simulation outputs.

INTRODUCTION

The maturity of cardiovascular simulation tools has significantly increased in recent years, offering new insight into disease progression and new capabilities for prediction of physiologic response to to numerical integration on multi-variate grids, computation of the pseudo-spectral stochastic coefficients using sparsity promoting greedy heuristics has been proposed in [START_REF] Doostan | A non-adapted sparse approximation of PDEs with stochastic inputs[END_REF]; extensions to multiresolution representations are discussed in [START_REF] Schiavazzi | Sparse multiresolution regression for uncertainty propagation[END_REF].

A few prior studies have begun to account for clinical data uncertainties in cardiovascular simulation. Sparse grid SC was used in [START_REF] Xiu | Parametric uncertainty analysis of pulse wave propagation in a model of a human arterial network[END_REF] to solve the stochastic differential equations governing the propagation of blood pressure in a realistic one-dimensional CFD model of the human circulation. A robust optimization algorithm, combining the surrogate management framework (SMF [START_REF] Marsden | A computational framework for derivative-free optimization of cardiovascular geometries[END_REF]) with SC on adaptive sparse grids, was applied in [START_REF] Sankaran | The impact of uncertainty on shape optimization of idealized bypass graft models in unsteady flow[END_REF] to parameterized models of bypass graft anastomoses, and in [START_REF] Sankaran | A stochastic collocation method for uncertainty quantification and propagation in cardiovascular simulations[END_REF] to several cardiovascular simulation examples. A similar approach was also applied [START_REF] Sankaran | An efficient framework for optimization and parameter sensitivity analysis in arterial growth and remodeling computations[END_REF] to growth and remodeling predictions in arteries. A systematic application of SC is discussed in [START_REF] Chen | Simulation-based uncertainty quantification of human arterial network hemodynamics[END_REF] where numerical experiments of increasing complexity were performed to quantify the impact of log-normally distributed input random inflow, boundary distal resistances, and vessel cross sectional area on blood pressures and flow rates at locations of interest.

In the above studies, the PDFs of the stochastic inputs were selected a priori, i.e., they were only partially inferred by the process of acquiring data in clinical practice. This should instead, in our opinion, be the starting point of any quantification of confidence in the results of patient-specific cardiovascular simulations. For example, input distributions of resistance or compliance elements in 0D circulation networks should be inferred from distributions in the clinical data determined from the patient's pathology and physiology. This will enable improved quantification of PDFs and correlations of input random boundary conditions, in contrast to the typical assumption of independence and identical distribution. Finally, we note that strategies combining clinical data assimilation with physics-informed Navier-stokes constraint have been proposed in [START_REF] D'elia | A variational data assimilation procedure for the incompressible navier-stokes equations in hemodynamics[END_REF]. Another example is provided in [START_REF] Moireau | Sequential identification of boundary support parameters in a fluid-structure vascular model using patient image data[END_REF] where visco-elastic support constants are inferred from clinical image data for a coupled fluid-structural model of the human aorta.

While applicable to other problems and disciplines, the present study is motivated by the pressing need for uncertainty quantification in multiscale models of single ventricle palliation [START_REF] Corsini | An integrated approach to patient-specific predictive modeling for single ventricle heart palliation[END_REF][START_REF] Kung | Predictive modeling of the virtual hemi-Fontan operation for second stage single ventricle palliation: two patient-specific cases[END_REF]. Complex cardiac malformations, such as Hypoplastic Left Heart Syndrome or Tricuspid Atresia, are characterized by the existence of only one functional ventricle and treated with a staged surgical transition to the Fontan circulation [START_REF] Fontan | Surgical repair of tricuspid atresia[END_REF]. As a direct transition to the Fontan circulation is contraindicated in the neonatal period due to high pulmonary vascular resistance, a three staged palliation is usually adopted, which allows the heart and lungs to progressively adapt to the new physiological paradigm: Norwood or BT shunt procedure (stage I), superior cavopulmonary connection (SCPC stage II, bi-directional Glenn or hemi-Fontan surgeries) and Fontan completion (stage III, lateral tunnel or extracardiac conduit). Numerical simulations were previously performed for a cohort of patients to understand the hemodynamic consequences of performing Stage II surgery with and without combined left pulmonary artery LPA arterioplasty on the resulting RPA/PA flow split ratio, average pulmonary pressure and mean pressure drop across the stenosis [START_REF] Schiavazzi | Hemodynamic effects of left pulmonary artery stenosis following superior cavopulmonary connection: a patient-specic multiscale modeling study[END_REF]. The boundary conditions for these simulations were determined using preoperative clinical data affected by uncertainty. In this study, we use the single ventricle example to present a framework for uncertainty quantification in cardiovascular simulations. For computational convenience, we break this question into two successive steps. The first step, formulated as an inverse problem, uses Bayesian inference and surrogate modeling to determine the sets of outlet boundary conditions matching the preoperative PDFs. The second step propagates these uncertainties through to simulation predictions using sparse grid stochastic collocation to determine the first two statistical moments for simulation outputs of interest.

While the individual steps of the proposed approach (e.g., Kriging, MCMC estimation, forward uncertainty propagation with stochastic collocation on sparse grids, etc.) are not new, a novel combination of inverse estimation and forward uncertainty propagation is proposed in this study, and then applied in detail to determine confidence of results in virtual surgery. This goes beyond prior applications of uncertainty in cardiovascular simulation, where an assumption of i.i.d. random inputs is often made without any direct physiological or clinical motivation, leading to an increase of the total uncertainty in the system. In the proposed approach, the distributions of uncertain inputs (boundary resistances in the selected application) are instead estimated directly from the available clinical data and their random uncertainties. The other methodological interest of this paper is the selection of tools in the context of three-dimensional Computational Fluid Dynamics able to maintain the computational cost acceptable.

The main context of this work is verification & validation of translational models with indirect clinical utility given by the ability to improve the statistical significance of conclusions drawn from numerical simulations. Direct clinical use of patient-specific simulations (e.g., as complementary clinical examinations) is also possible to compare, for example, the hemodynamic effects of alternative surgical choices. In this context, quantification of confidence in the predictions will constitute an indispensable aid to draw significant conclusions.

The paper is organized as follows. The clinical history of patient MUSC2 is summarized in Section 2. This section also describes the four CFD models used to investigate pre-and postoperative surgical conditions for this patient. Section 3 is devoted to quantifying the target PDFs of right flow split ratio ρ(f s ) and average pulmonary pressure ρ(p av ) from patient-specific clinical data. Section 4 discusses estimations of the distributions of pre-operative boundary resistances from ρ(f s ) and ρ(p av ), followed by forward propagation to virtual surgery results using SC. Finally, confidence intervals on global and local hemodynamic quantities are presented in Section 5 and discussed in Section 6.

CLINICAL DATA AND MODEL CONSTRUCTION

Patient data

One patient, MUSC2, was recruited at the Medical University of South Carolina. This patient was diagnosed with Hypoplastic Left Heart Syndrome (HLHS) with mitral and aortic valve stenosis. A 3.5 mm right-modified Blalock-Taussig shunt (rmBTS) was inserted with the Norwood procedure. Development of central PA stenosis was observed following Stage I surgery. At 5 months of age, with a BSA of 0.3 m 2 , the patient underwent Stage II surgery with augmentation of the central PA using pulmonary homograft. Pressure catheter and MRI flow measurements were acquired 1.5 months prior to surgery under the same general anesthetic. All imaging and clinical data was acquired in accordance with institutional review board guidelines.

Model creation

Three patient-specific models were created to simulate pre-operative stage I hemodynamics, and predict the outcome of Stage II surgery. A baseline geometry for model A was obtained using semi-automatic segmentation from the pre-operative (stage I) MRI scans. Stage II hemi-Fontan geometries with and without LPA stenosis were successively created through virtual surgery, resulting in models B and D. Note that models B and C share the same geometry, but different boundary conditions were applied to investigate how open-loop or closed-loop formulations affect the simulation results. The three-dimensional geometries for the four models are illustrated in Figure 1.

Boundary conditions

Steady flow, open loop boundary conditions were applied to models A and B, with input flow rate of 9.3 ml/s prescribed at the shunt/SVC inlet and boundary resistances applied to each outlet.

Pulsatile closed loop boundary conditions were applied to models C and D. A lumped parameter network (LPN) model of the Stage II circulation was formulated as a system of algebraic and ordinary differential equations (ODEs) characterizing an assembly of resistance, compliance and inductance elements, consistent with an electrical circuit analogy. Five main circuit blocks are included characterizing the heart, systemic (i.e., upper and lower body) and pulmonary (i.e., left and right) circulations. Time-varying activation functions model atrial and ventricular contraction while non-linear diodes mimic atrio-ventricular and aortic valves. The lower body circulation model includes the kidneys, liver and intestine together with a venous valve in the leg block. Coronary circulation, gravity and respiratory effects were neglected [START_REF] Migliavacca | Modeling of the norwood circulation: effects of shunt size, vascular resistances, and heart rate[END_REF]. A semi-automatic registration was performed to identify the parameters of the LPN model, as discussed in our previous work [START_REF] Corsini | An integrated approach to patient-specific predictive modeling for single ventricle heart palliation[END_REF][START_REF] Kung | Predictive modeling of the virtual hemi-Fontan operation for second stage single ventricle palliation: two patient-specific cases[END_REF] to match pre-operative clinical data. Equivalent pulmonary RCR blocks for models C and D were determined by combining sets of outlet boundary resistances with pulmonary morphometry data [START_REF] Spilker | Morphometry-based impedance boundary conditions for patient-specific modeling of blood flow in pulmonary arteries[END_REF][START_REF] Troianowski | Three-dimensional simulations in Glenn patients: clinically based boundary conditions, hemodynamic results and sensitivity to input data[END_REF]. Table I summarizes the geometry and boundary conditions used for the four models employed in the present study.

Solution of the 3D Model

Simulations were performed using a custom three-dimensional Streamline Upwind Petrov-Galerkin (SUPG) finite element Navier Stokes solver from the SimVascular project [START_REF] Taylor | Finite element modeling of blood flow in arteries[END_REF], implicitly coupled to the LPN network using Neumann boundary conditions at all inlets and outlets [START_REF] Moghadam | A modular numerical method for implicit 0D/3D coupling in cardiovascular finite element simulations[END_REF]. Implementation of resistance boundary conditions for open-loop models follows the coupled multi-domain method discussed in [START_REF] Vignon-Clementel | Outflow boundary conditions for three-dimensional finite element modeling of blood flow and pressure in arteries[END_REF]. Blood was assumed to be Newtonian with a density of 1060 kg/m 3 and a dynamic viscosity of 0.004 Pa•s, and rigid vessel material properties were assumed. Stabilization was applied at coupled velocity-pressure surfaces to prevent numerical divergence due to backflow [START_REF] Moghadam | A comparison of outlet boundary treatments for prevention of backflow divergence with relevance to blood flow simulations[END_REF]. The discrete algebraic system of equations was solved iteratively using a bipartitioned algorithm with a dedicated preconditioner [START_REF] Esmaily-Moghadam | A new preconditioning technique for implicitly coupled multidomain simulations with applications to hemodynamics[END_REF].

Output quantities of interest

We consider a bounded three-dimensional region Γ ⊂ R 3 containing a subset of the human vasculature, with piecewise smooth boundary ∂Γ and parameterized in space and time by the coordinates x ∈ R 3 and scalar t ∈ R + . The vectors

Q(t) ∈ R n = (Q 1 (t), Q 2 (t), . . . , Q n (t)), P(t) ∈
R n = (P 1 (t), P 2 (t), . . . , P n (t)) represent flow rates and pressure at each outlet ∂Γ i , i ∈ J = {1, . . . , n}. If we denote as p(x, t) and v(x, t) the pressure and velocity field in ∂Γ i , respectively, the above quantities are obtained as:

Q i (t) = ∂Γi [v(x, t) • n i ] d ∂Γ i , and P i (t) = 1 |∂Γ i | ∂Γi p(x, t) d ∂Γ i , (1)
where n i is the outward normal to ∂Γ i (assumed planar) and |∂Γ i | its area. Notation is simplified in the following sections not explicitly reporting time dependence for these quantities. Left and right outlets are identified using the two index sets J l = {l 1 , l 2 , . . . , l n l } ⊂ J and J r = {r 1 , r 2 , . . . , r nr } ⊂ J, with J l ∩ J r = ∅. The total flow rate through left and right pulmonary branches is given by q l and q r , while p l and p r denote the associated average pressures. These quantities are calculated as follows:

q l = j∈J l Q j , q r = j∈Jr Q j , p l = 1 n l j∈J l P j , p r = 1 n r j∈Jr P j . (2)
RPA/PA flow split ratio (f s ), average pulmonary pressure (p av ) and pressure drop (δ p ) are then obtained as:

f s = q r q l + q r , p av = 1 2 (p l + p r ) , δp = |p l -p r |. (3)
Outlet resistances are defined using the vector R ∈ R n = (R 1 , R 2 , . . . , R n ) = ((P 1 -P a )/Q 1 , (P 2 -P a )/Q 2 , . . . , (P n -P a )/Q n ) with P a being the atrial pressure.

Assuming all other boundary and initial conditions to be fixed, a CFD analysis defines a relationship between values of boundary resistances and outlet flow and pressure. In the following sections, we will define a forward problem that determines f s , p av and δp from R, and an inverse problem that determines R given f s , p av and δp. Formally, we can introduce a non-linear operator G : R n → R 3 , and express the foward (F) and inverse (I) problems as:

F :    f s p av δp    = G(R), I : R = G -1    f s p av δp    (4)
3. ASSESSMENT OF CLINICAL DATA UNCERTAINTY Systematic errors in repeated clinical measurements affect numerical predictions and should be clearly associated with simulation results. It is therefore important to realistically estimate appropriate error distributions for the clinically measured quantities of interest, i.e., mean pulmonary pressure (p av ), RPA/PA pulmonary flow split ratio (f s ), and model geometry.

Pulmonary pressure measurements

Mean PA pressure for patient MUSC2 was determined indirectly, via measurement of the pulmonary venous wedge pressure (PVWP). Indirect assessment is a common practice in cases where it is difficult or risky to access the PAs with a catheter. For patients post Stage I single ventricle palliation surgery with aortopulmonary bypass, for example, the pulmonary arteries can only be reached by inserting a catheter through the BT shunt. This is a small cross-sectional area conduit carrying the entire pulmonary blood supply, and partial occlusion with a catheter can cause hypoxia and other severe complications.

There are various possible sources of errors in the acquisition of mean PVWP. Catheters equipped with an inflatable balloon are routinely used to provided adequate wedging. Partial vein occlusion with end-hole catheters may result in errors related to the difference between static and dynamic pressures. For single ventricle patients, respiration affects venous return [START_REF] Hsia | Effects of respiration and gravity on infradiaphragmatic venous flow in normal and fontan patients[END_REF]. As end-expiration correlates well with atmospheric pleural pressure, PVWP measurements should use end-expiratory breath hold. For patients with aortopulmonary shunts, however, errors in PVWP related to the respiration are mitigated by the arterial nature of the pulmonary veins (PV) traces and associated lower compliance.

While adequate wedging and proper accounting of the effects of respiration should produce accurate PVWP measurements, correlation between p av and PVWP may still be affected by the patient's physiology. For example, heart rate may affect the difference between PA and PV pressure due to a selective damping action of the capillary vascular bed [START_REF] Mori | Relation of pulmonary venous wedge pressures to pulmonary artery pressures in patients with single ventricle physiology[END_REF]. The vascular bed also reacts to changes in the flow conditions due to wedging, with consequences on the agreement between sequential PA-PV pressure measurements. Finally, partial occlusion of the shunt may induce falsely low p av and affect the correlation with PVWP.

The relation between PVWP and p av is investigated in [START_REF] Mori | Relation of pulmonary venous wedge pressures to pulmonary artery pressures in patients with single ventricle physiology[END_REF][START_REF] Thompson | Pulmonary venous wedge pressures accurately predict pulmonary arterial pressures in children with single ventricle physiology[END_REF] for single-ventricle patients and in [START_REF] Hill | Pulmonary venous wedge pressure provides a safe and accurate estimate of pulmonary arterial pressure in children with shunt-dependent pulmonary blood flow[END_REF] for children with shunt-dependent pulmonary blood flow. Good agreement between PVWP and p av is reported in [START_REF] Thompson | Pulmonary venous wedge pressures accurately predict pulmonary arterial pressures in children with single ventricle physiology[END_REF] where differences of -0.15 ± 1.31 mmHg are observed for patients with BT shunt not affected by pulmonary hypertension. In [START_REF] Mori | Relation of pulmonary venous wedge pressures to pulmonary artery pressures in patients with single ventricle physiology[END_REF], the correlation between PVWP and pulmonary arterial pressure (PAP) was found to be affeted by the measured pressure range, with PVWP significantly understimating PAP above 18 mmHg.

For patient MUSC2, measured left PVWP was 13.0 mmHg, right upper PVWP was 11.0 mmHg and left lower PVWP was 12.0 mmHg. Based on the above data, a value of 12.0 mmHg was clinically reported for both left and pulmonary arteries, with a maximum associated error of 1.0 mmHg, which is compatible with literature values. Because pulmonary stenosis can significantly affect the distribution of pressure in the PAs, we accounted for effects of proximal and distal stenosis in the pulmonary arteries by rounding the mean pulmonary pressure variance to 3.0 mmHg 2 (standard deviation of 1.73 mmHg).

Flow split measurement

Various sources of errors in estimation of blood velocity through phase-constrast MRI (PC-MRI) are discussed in the literature [START_REF] Fratz | Guidelines and protocols for cardiovascular magnetic resonance in children and adults with congenital heart disease: SCMR expert consensus group on congenital heart disease[END_REF]. First, the cross section of the vessel of interest should be compatible with the characteristics of the PC-MRI scanner, where an in-plane resolution of 1.0-1.3 mm and a slice thickness of 5.0 mm are typically used for children. Changes in size of the pulmonary arteries may not be negligible with the heart cycle and should also be considered. In this regard, prior studies on small animals reported a change in the PA cross sectional area of about 12% [START_REF] Loeber | How much does aortic and pulmonary artery area vary during the cardiac cycle?[END_REF]. For patient MUSC2, both the diameter of the BT shunt (3.5 mm) and the size of the LPA subject to stenosis are close to the limit of 3 cells across in diameter and 8 cells in cross section as discussed in the literature [START_REF] Fratz | Guidelines and protocols for cardiovascular magnetic resonance in children and adults with congenital heart disease: SCMR expert consensus group on congenital heart disease[END_REF]. Partial volume effects due to cells containing both blood and vessel wall may also be responsible for decreased accuracy in the flow measurements.

High velocity/turbulent flow produces signal loss which negatively affects the flow measurement through the PAs. Peak velocities were observed for patient MUSC2 equal to 43 cm/s for the RPA and 65 cm/s for the LPA. Higher velocities and more flow in the LPA may therefore lead to signal loss with possible underestimation of the clinically acquired flow rate. Also, prior studies have reported inter-observer errors in pulmonary artery MRI flow measurements in the range of 5-10%, with variation in cardiac output with repeated measurements reported in the range of ±6%.

For patient MUSC2 the measured flow was 4.5 ml/s for the RPA, and 5.23 ml/s for the LPA, resulting in a RPA/PA flow split ratio of 0.46 and total flow through the BT shunt equal to 9.73 ml/s. Similarly, measured flow rates in the pulmonary veins were 4.51 ml/s for the RPV, 5.88 ml/s for the LPV, with flow split ratio equal to 0.43 and total flow 10.39 ml/s. Flow due to development of collateral circulation was not observed. The final distribution was selected by accounting for 10% error in the measured arterial flow split due to the size of the PAs near the anastomisis with the shunt, comparable with the resolution of the PC-MRI scan. Moreover, a right skew was assigned to the RPA/PA flow split distribution as arterial flow in the LPA was likely underestimated due to signal loss. The mode of the distribution (i.e, the most likely flow split associated with the maximum of the PDF) was set to 0.43 due to the higher reliability of venous flow measurements for this patient.

Geometry

The geometry of the shunt/SVC-PAs junction was not parameterized in the current study. However, it is useful to discuss possible geometrical errors to better understand differences between clinically measured and simulated hemodynamics. Analyzing three model configurations (pre-operative, postoperative with LPA stenosis and post-operative without LPA stenosis) will provide further insight on the effects of geometry perturbations on our numerical results.

The pre-operative shunt-PAs junction model was reconstructed from three-dimensional steadystate free precession images and MR angiograms with resolution of 0.391 mm × 0.391 mm × 1 mm and 1.0 mm × 1.0 mm × 2 mm, respectively, acquired at the diastasis point (mid-diastole). MR anatomical measurements are routinely confirmed with two-dimensional planar black blood MR images (with resolution of 1.042 × 1.042 mm 2 and 4 mm slice thickness) used to enhance vessel lumen contrast. For the MUSC2 patient, both shunt and central PA stenoses were further observed using fluoroscopy images from the catheterization lab, providing measures of shunt and PAs to the nearest 0.2 to 0.3 mm.

Turbulence or high velocities in the flow may affect the measurements making the BT shunt or PA stenoses appear smaller than they are. This seems to affect the present case where high velocities and mixing are observed in the anastomosis region.

Final distribution of clinical data

From the above discussion, the error distributions associated to f s and p av were determined as follows. The distribution ρ(f s ) was selected as a skewed Gaussian distribution with range [0.41, 0.55] and mode equal to 0.43. This is characterized by three parameters ξ (location), ω (scale) and α (shape) in the expression [START_REF] Azzalini | The multivariate skew-normal distribution[END_REF]:

ρ(f s ) = 1 ωπ e -(fs-ξ) 2 2ω 2 α( fs-ξ ω ) -∞ e -t 2 2 dt (5)
The parameters ξ, ω and α were identified through Simulated Annealing with constraints µ(f s ) -2.0 σ(f s ) = 0.41, µ(f s ) + 2.0 σ(f s ) = 0.55, and arg max fs ρ(f s ) = 0.43, leading to ξ = 0.423, ω = 0.053 and α = 4.0. The symbols µ(•) and σ(•) denote mean and standard deviation, respectively. A Gaussian distribution was selected for ρ(p av ) with mean value of 12.0 mmHg and standard deviation equal to √ 3 = 1.73 mmHg. These distributions are illustrated in Figure 2. Finally, we note that the above distributions were assumed to be Gaussian, and that moments or other characteristics (e.g., mean value, mode, variance and skewness) were chosen to be consistent with clinical data availability for the specific patient pathology, clinical experience and findings in the literature.

UNCERTAINTY ANALYSIS

Applicability of the proposed approach to other fields

The general applicability of the proposed approach to problems other than single-ventricle virtual surgery is discussed in this section. The exercise of determining the probability distributions of specific measurements, can certainty be performed in other application areas, and such quantities may actually be more readily available in other settings. These tools could be applied to other problems in cardiovascular disease, for example when the uncertainty relates to a measured flow rate at the outlet of a vascular model or pressures are acquired under uncertainty. In all these case an important common feature is that clinical data are quantities that relate to model outputs (i.e., flow split or average pressure over a region, etc.), not inputs.

Once the target uncertainty is assessed, Bayesian parameter estimation (discussed in a clinical context in Section 4.5) can be performed with arbitrary models. However, difficulties may arise related to possible non identifiable parameters or excessive computational cost. Non identifiability (see, e.g., [START_REF] Miao | On identifiability of nonlinear ODE models and applications in viral dynamics[END_REF]) is typically associated with over-parameterization, i.e., inability of learning the distributions of a large number of parameters under data scarcity. In other words, equally likely solutions to a given problem may arise, without the ability to assess their difference from the available data. Identifiability depends on the parameterization of the model and, in a Bayesian context, the amount of prior information at hand. For problems that are not identifiable, an analysis of the Fisher information matrix may be useful in detecting unimportant parameters or the presence of non identifiable combinations in the neighborhood of specific parameter combinations, but other approaches to detect local or global identifiability may be used. The proposed application of inverse estimation for virtual surgery discussed in Section 4.5 is not affected by this issue as illustrated in the MCMC results of Figure 7.

Conversely, we propose a solution to the problem of an expensive model solution in Section 4.4, by replacement with an inexpensive Gaussian process interpolant. This lightweight surrogate does not need to be accurate for all the possible values of the inputs, but only in the region of interest, i.e., input values generating model solutions sufficiently close to the collected data. Similar Kriging models have been used to accelerate optimization in a wide range of engineering applications and uncertainty analysis [START_REF] Sankaran | The impact of uncertainty on shape optimization of idealized bypass graft models in unsteady flow[END_REF][START_REF] Sankaran | An efficient framework for optimization and parameter sensitivity analysis in arterial growth and remodeling computations[END_REF].

At this point, forward uncertainty propagation can be performed in full generality (with nonintrusive approaches directly applicable to existing solvers), with the only additional complication of properly transforming the MCMC input parameter samples, e.g., as realizations from an i.i.d. random vector. This is discussed in Section 4.6 for the proposed application to virtual surgery.

Section outline

After selecting clinically appropriate distributions for f s and p av in Section 3.4, two main steps are needed to propagate this uncertainty to virtual surgery simulation outputs. First, boundary conditions need to be tuned to identify parameters that produce the distributions of clinical data ρ(f s ) and ρ(p av ) identified for the pre-operative model. In other words, we aim to determine the joint distribution ρ(R) and covariance C R such that the flow split ratio and average pulmonary pressures resulting from G(R) match the clinical data distributions ρ(f s ) and ρ(p av ). We formulate this as an inverse problem using MCMC to compute ρ(R) and C R , as discussed in section 4.5.

Second, we need to propagate the uncertainty in the random vector R to statistically characterize virtual surgery outputs. This step employs SC as discussed in Section 4.6. Finally, we stress that every solution of the forward problem (i.e. 3D multiscale simulation) is expensive, therefore the number of evaluations of G(R) needs to be kept to a minimum throughout both inverse and forward propagations. One approach is to construct a surrogate model that approximates G sufficiently well for the range of boundary resistances associated with non zero ρ(R), as discussed in the next section.

Construction of a surrogate model for f s and p av : semi-empirical sampling

Initial samples in the region of non-negligible ρ(R) are determined using a semi-empirical approximation based on a 0D fluid dynamics model that mimics G -1 , starting from sets of f s and p av compatible with the clinical data. Depending on how detailed this model is and how many samples are available, our approximation of G -1 can be progressively refined. A first basic approach to generate samples, is to assume a constant pressure p av throughout the SVC-PA junction, and determine the components of Q from the following expression:

Q i = Q in S i (δ i,r f s + (1 -δ i,r )(1 -f s )) j S j δ i,j (6)
where δ i,j = 1 if outlets i, j ∈ J l or i, j ∈ J r , otherwise δ i,j = 0 and δ i,r = 1 if outlet i ∈ J r , zero otherwise. The vector S = (S 1 , S 2 , . . . , S n ) contains the surface area of each outlet. The underlying assumption in the derivation of Equation ( 6) is that the volumetric flow rate is proportional to the surface area at each outlet. Once the flow rates are determined, the vector R can be defined using the expression:

R i = (p av -P a )/Q i ∀ i ∈ {1, 2, . . . , n} (7)
The approximation in ( 6) is evident as no major pressure losses due to friction in the 3D model nor minor losses due to RPA and LPA stenoses are considered. A straightforward improvement accounts for minor pressure losses due to stenosis in the LPA and RPA, respectively, distal and proximal to the anastomosis between the rmBTS and PAs. To quantify the relationship between pressure losses and blood velocity in the region affected by hypoplasia, we adopt the formulation discussed in [START_REF] Seeley | Effect of geometry on pressure losses across models of arterial stenoses[END_REF]:

∆p s = µL s 32 S 2 v D 2 v S 2 s V + 0.76 ρ S v S s -1 2 V 2 , (8)
which quantifies the minor pressure loss through a cylindrical vessel with diameter D v and area S v affected by a discrete stenosis with length L s and minimal area S s . The quantities µ and ρ denote blood viscosity and density, respectively. If we assume ∆p l s and ∆p r s to represent the pressure losses through LPA and RPA stenosis, respectively, the resistances on the left and right PAs are equal to:

R l i = p av -0.5∆p l s -P a Q i , R r i = p av + 0.5∆p l s -P a Q i ( 9
)
when only the stenosis in the LPA is taken into account and

R l i = p av -0.5(∆p l s -∆p r s ) -P a Q i , R r i = p av + 0.5(∆p l s -∆p r s ) -P a Q i ( 10
)
accounting for both LPA and RPA stenosis. The initial samples were obtained by selecting a set of pairs (f s , p av ) belonging to the region of interest in the flow split -average pressure plane. Rather than adopting uniform or stratified sampling (see, e.g., [START_REF] Mckay | Comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF]), we preferred to determine these locations on a two-dimensional sparse grid assembled from nested univariate quadrature rules. An order 7 grid was selected consisting of 65 points. Examples of two-dimensional quadrature grids suited for uniform [START_REF] Petras | Smolyak cubature of given polynomial degree with few nodes for increasing dimension[END_REF] and Gaussian [START_REF] Genz | Fully symmetric interpolatory rules for multiple integrals over infinite regions with Gaussian weight[END_REF] probability measures, respectively, are illustrated in Figure 3. We determined 65 sets of n = 20 resistances for the pre-operative model using expression [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF] and another 65 sets using expression [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF]. After performing m = 130 model evaluations, the values of f s and p av were computed and compared to the initial sparse grid locations. Figure 4 shows the agreement between desired and computed values for these two quantities.

Results are reported for both semi-empirical formulations with stenosis in the LPA only and both LPA and RPA. Inclusion of RPA stenosis improves the accuracy in the proposed 0D semiempirical formulation for increasing RPA/PA flow split ratios and for average pulmonary pressures (Figures 4a and4b). Figure 4c quantifies the unsteadiness in computing the pre-operative RPA/PA flow split ratio by plotting the 90% confidence interval for this quantity, i.e., 2 • 1.645 • σ fs , where

σ fs = tmax t0 [f s (t) -µ fs ] 2 dt 1/2
and σ fs , µ fs represent the standard deviation and time-averaged values for the RPA/PA flow split ratio, respectively. The quantity t max denotes the total simulation time and t 0 was selected to exclude the initial transient.

The amount of flow through the LPA affects σ fs , and lower values of RPA/PA ratios produce oscillations with amplitude of 0.03 -0.04. This phenomenon, which is mainly produced by the geometrical complexity of the flow domain, will be further discussed when assessing the accuracy for distributions of pre-operative simulation outputs.

Denoting R = {R i , i = 1, . . . , m} as the sets of resistances used as pre-operative simulation inputs, we can form their empirical covariance matrix C R using the following expression: As expected, the resistances in R belong to a two-dimensional manifold generated by r 1 , r 2 ∈ R 20 , i.e., the eigenvectors of C R associated with non-zero eigenvalues. All other resistances can be obtained from linear combinations of r 1 and r 2 using coefficients β 1 and β 2 . This reduces the dimensionality of the space of outlet resistances from n = 20 to only two variables, namely β 1 and β 2 , a direct consequence of assuming only right flow split ratio and average pressures as independent variables in our 0D semi-empirical formulation. This is physically consistent with the fact that there are two main mechanisms to change the resistances and produce an effect on the average pressure and flow split ratio in the shunt-PA junction model. The first mechanism is to increase or decrease left and right pulmonary resistances simultaneously, thus mainly impacting the average pulmonary pressure. The second mechanism is to change the relative difference between left and right pulmonary resistance, keeping the same total value, which directly affects flow split ratio. Other variations outside the above patterns produce minor effects on f s and p av , primarily affecting the relative distribution to individual branches. To better illustrate this point, Figure 5a plots the average left vs. right pulmonary pressures for all 130 simulations. The sparse grid in Figure 3 is reproduced in Figure 5a, rotated 45 degrees. This confirms the correlation between the independent variables of our semi-empitical approximation f s and p av and variations in the left and right preoperative pumonary artery pressures, respectively. Finally, Figure 5b shows that pressure values are nearly uniform across outlets in all simulations, confirming that average pressure at the outlets is sufficiently well approximated by the sample mean.

Construction of a Kriging surrogate

From the results of the pre-operative simulations, and after a principal component analysis of R, we have determined m evaluations of a relationship Gi : R 2 → R, i = 1, 2, 3 mapping realizations of the coefficients β j = (β j 1 , β j 2 ), j = 1, . . . , m into values for the right pulmonary flow splits (i = 1), average pressures (i = 2) and mean pressure drop across the LPA (i = 3). The next problem we need to address is how to determine the joint distribution ρ(β) that generates the distributions of interest for f s , p av and δp. Monte Carlo strategies have proven successful to solve inverse Bayesian estimation problems (see,e.g., [START_REF] Gilks | Markov chain monte carlo[END_REF]) but they generally require a significant number of evaluations for Gi . To overcome this problem, thus significantly reducing the computational burden of estimating ρ(β), we replaced Gi with a Kriging surrogate. Consider the generic quantity y j ∈ R with y 1 = f s , y 2 = p av , or y 3 = δp and y = (y 1 , y 2 , y 3 ). Kriging selects a surrogate model for each variable y j (β) from a stochastic process that satisfies the following linear regression model:

y j (β) = p i=1 c i,j φ i (β) + e j (β), j = 1, 2, 3. ( 12
)
where the regression error is assumed to be a Gaussian process with zero mean and covariance C ej (e j (β i ), e j (β k )) = C ej (θ, β i -β k ). Among the covariance models proposed in the literature (see, e.g. [START_REF] Acar | Effects of the correlation model, the trend model, and the number of training points on the accuracy of Kriging metamodels[END_REF]), we adopt the typical choice of an exponential covariance formulation with

C ej (θ, β i -β k ) = m l=1 exp(-θ l β i -β k ).
Our generic simulation result y j is affected by intrinsic unsteadiness and therefore our covariance model was selected to guarantee the fastest possible decay to a smooth mean polynomial trend. Accordingly, a family of second order polynomial functions was selected for φ i (β), i = 1, . . . , p. The vector of parameters θ is determined using maximum likelihood estimation and the coefficients c i,j , i = 1, . . . , p by solving a generalized least squares problem with a possibly dense correlation matrix. We refer the interested reader to [START_REF] Stein | Interpolation of spatial data: some theory for Kriging[END_REF] for an in-depth discussion on the theoretical background and computational implementation of Kriging. Figure 6 shows a graphical representation of the Kriging surrogates for the three quantities of interest. Finally, we note that a combination of Kriging surrogates and derivative-free optimization has proven successful in robust optimization of virtual surgery scenarios [START_REF] Marsden | A computational framework for derivative-free optimization of cardiovascular geometries[END_REF][START_REF] Yang | Constrained optimization of an idealized y-shaped baffle for the Fontan surgery at rest and exercise[END_REF][START_REF] Yang | Hepatic blood flow distribution and performance in conventional and novel y-graft Fontan geometries: a case series computational fluid dynamics study[END_REF].

Inverse Bayesian analysis and MCMC

We formulate our estimation problem for a generic scalar model result y j using the following nonlinear statistical model:

y j = G * j (β) + e j (13)
and assume each error component e j to be Gaussianly distibuted with zero mean and variance σ 2 j . In [START_REF] Tierney | Markov chains for exploring posterior distributions[END_REF] the symbol G * j referes to the Kriging surrogate for Gj . We use MCMC (see, e.g., [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF][START_REF] Gilks | Markov chain monte carlo[END_REF]) to sample from the posterior distribution ρ(β|y) defined as:

ρ(β|y) = ρ(y|β)ρ(β) ρ(y) (14)
Due to the non-linearity in G * j and the consequent non-Gaussianity in the likelihood ρ(y|β) we adopted a Metropolis-Hastings sampler [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF] with a symmetric multivariate Gaussian trial distribution. In practice, samples belonging to the posterior ρ(β|y) are progressively generated according to a Markov Chain β k , k = 1, . . . , k max , where the probability of transition to a new point β k+1 from β k is:

ρ M H = min 1, ρ(β k+1 |y) ρ(β k |y) (15)
Figure 7 shows the MCMC tracings, ergodic statistics and marginal PDFs for all accepted β k . Both distributions are unimodal and good mixing can be observed from MCMC traces. Satisfactory convergence can also be observed from plots of ergodic mean and standard deviation. It should also be noted that, in general, β 1 and β 2 resulting from MCMC are correlated.

To verify that we had obtained the correct estimation of ρ(β 1 , β 2 ), we determined the histograms of f s and p av resulting from applying the Kriging model to the accepted β k MCMC samples. As expected, we retrieved back ρ(f s ) and ρ(p av ) as defined in Section 3, and shown in Figure 2. We thus confirmed that we had obtained a distribution of boundary resistances that correctly reproduced the known distribution of clinical data.

To summarize, in this section we have discussed an MCMC approach to sample from ρ(β 1 , β 2 ) such that forward propagation of associated boundary resistances provides pre-operative model results compatible with ρ(f s ) and ρ(p av ).

Forward propagation to multiscale model results

After obtaining the distributions of outlet resistances, parametrized through β 1 and β 2 , we proceeded to propagate them forward to the results of our four CFD simulations. Consider a complete probability space (Ω, F, P) in which Ω is the sample space, F is a σ-algebra of possible events, and P denotes a probability measure on F. The input uncertainty is characterized by the vector of independent random variables β = ( β1 , β2 ), where each random variable βi : +∞] has distribution ρ( βi ). In the following we assume a Gaussian distribution for the generic βi . Let Γ ⊂ R D , D ∈ N, be the spatial domain with boundary ∂Γ. The quantities x ∈ Γ and t ∈ [0, T ] represent the spatial and temporal variables, respectively. We consider the problem of computing expectations (in particular first and second order moments) of the solution u(x, t, β) : Γ × [0, T ] × Σ → R h , h ∈ N, to our stochastic CFD problem. These expectations, for a particular location in space x * and time t * , can be computed through numerical integration as follows:

Ω i → Σ i ≡ [-∞,
E u(x * , t * , β) = Σ u(x * , t * , β)ρ( β) dΣ ≈ ng i=1 u(x * , t * , βi ) w i V u(x * , t * , β) = Σ u(x * , t * , β) -E u(x * , t * , β) 2 ρ( β) dΣ ≈ ng i=1 u(x * , t * , βi ) -E u(x * , t * , βi ) 2 w i , (16)
where V [•] denotes the variance operator and w i , i = 1, . . . , n g are the integration weights. The numerical integrals in [START_REF] Ghanem | Stochastic finite elements: a spectral approach[END_REF] result in accurate estimates for a sufficiently smooth u(x * , t * , β), i.e., this solution can be approximated sufficiently well by only the first few terms of the following expansion:

u(x * , t * , β) = p i=1 c i H i ( β) (17)
for an arbitrary spatial location x * and time t * . In [START_REF] Nobile | A sparse grid stochastic collocation method for partial differential equations with random input data[END_REF] the family of tensor product polynomials H i ( β) is selected to be orthogonal with respect to the distribution ρ( β), i.e., Σ H i ( β) H j ( β)ρ( β) dΣ = 0, if i = j. The input random vector β needs to be introduced, as β (resulting from MCMC estimation) contains, in general, non-normalized (i.e. with zero average and unit variance) and correlated random variables. The burn-in (selected as 10% of the samples) is first removed and the empirical covariance matrix C β computed and factorized in the form C β = S T S, where S is the upper triangular Cholesky factor. The intermediate decorrelated random vector β is then computed as β = S -T β and independence is achieved by truncation to the first two moments, i.e., µ( β1 ), µ( β2 ), σ 2 ( β1 ), σ 2 ( β2 ).

4.6.1. Integration on Smolyak sparse grids Selection of the numerical integration scheme, i.e., the set of points { βk , k = 1, . . . , n g } to be used in ( 16) constitutes the next step of our analysis. Following the Smolyak construction [START_REF] Smolyak | Quadrature and interpolation formulas for tensor products of certain classes of functions[END_REF], the integration grid consists of a partial tensorization of univariate nested quadrature rules. Consider an n-dimensional function f ( β1 , . . . , βn ) : Σ ⊂ R n → R where the stochastic domain Σ = Σ 1 × • • • × Σ n is generated by Cartesian products of one dimensional intervals. Define Σ j as the interval resulting from keeping all dimensions fixed except β j . A univariate nested quadrature rule of order q is defined as follows: Σj f ( β1 , . . . , βj , . . . , βd )ρ( βj ) dΣ j ≈ nq i=1 f ( β1 , . . . , βi j , . . . , βd )

w i = Q q j [f ] N k j [f ] = Q k+1 j -Q k j [f ] where Q 0 j [f ] = 0. ( 18
)
By tensorization of nested univariate quadrature formulae, we assemble the Smolyak quadrature rule by writing:

S k [f ] = N k1 1 ⊗ • • • ⊗ N k d d [f ] and Σ f ( β) ρ( β) dΣ ≈ |k|<q+n-1 S k [f ], (19)
where |k| = n i=1 k i . The truncation of ρ( β) to a multivariate Guassian distribution allows the quadrature locations { βi , i = 1, . . . , n g } to be selected according to well known formulae, but unfortunately the support of the distribution becomes infinite. Thus, values of resistances determined by first de-normalization βi = µ( βi ) + βi σ( βi ), and then correlation with S -1 and linear combination R k = d j=1 β j r j may become negative and therefore not acceptable as boundary conditions. To overcome this problem, the standard Gaussian distribution ρ( β) is truncated to the interval [-4, 4] × [-4, 4] and a Gauss-Legendre sparse grid is used. To make sure that a unit constant function is correctly integrated over this interval for any order of the sparse grid, we also modify the integration weights such that:

ng i=1 w i d j=1 1 √ 2 π e -( βi j ) 2 2 = 1. (20)
4.6.2. Computing marginal distributions for models results In this section we discuss possible strategies to determine marginal distributions for the quantities of interest from pre-and postoperative model solutions at the n g input quadrature locations. A combined Monte Carlo and interpolation strategy is used. Interpolation can be perfomed by Kriging as discussed in Section 4.4 or alternatively using numerical integration to compute the coefficients of a multivariate Legendre polynomial expansion:

f ( β) ρ( β) ≈ p i=1 α i L i ( β), α i = f ( β) ρ( β), L i ( β) = Σ f ( β) L i ( β) ρ( β) dΣ, (21)
where the orthogonality of the polynomial family L i ( β) with respect to the uniform probability measure is used. Both approaches lead to very similar results for all models considered in this study.

Computational tasks and associated cost

A graphical layout for the tasks discussed in the previous sections is illustrated in Figure 8. The initial clinical data manipulation task (i.e., T-C) is not computationally expensive; only optimization (performed using Simulated Annealing) was required to identify the parameters of the flow split skewed Gaussian distribution, as discussed in Section 3.4. Interrogation of the 0D surrogate model to determine the set R and its principal component analysis (i.e., T-E1 and T-E2) can also be performed at a reduced cost with many, freely available, software libraries. The task of constructing the Kriging pre-operative model approximant is computationally expensive as it requires multiple model evaluations. Note that the number of model evaluations is the main contributor to computational cost as expressed in Table II, where this quantity is reported to illustrate the overall cost of the proposed approach. For task T-E3 we performed 65+65 model evaluations using outlet resistances determined through two 0D models accounting for stenosis in the LPA only and both LPA and RPA. Note that the number of required model solutions in this tasks is problem dependent and mainly affected by the smoothness of G i (β). An additional 65 model evaluations were used for forward uncertainty propagation using the pre-operative model to check the agreement between computed and prescribed clinical output distributions, as shown in Figure 9. This cost was added to Table II, even if not strictly necessary.

Finally, forward propagation was performed using 65 post-operative model solutions. We also note that the total execution time for simultaneous solution of multiple computer models results in only a modest increase in wall clock time in modern computer clusters with respect to a single model evaluation, subject to resource availability. In our application, we ran 65 simulations simultaneously using a large parallel cluster through the XSEDE program. T-E4: Run MCMC using the Kriging surrogate. See Section 4.5.

Parameter Estimation

T-P1: Transform inputs (i.e., BC parameter set) to an i.i.d. random vector. See Section 4.6.

T-P2: Use stochastic collocation on sparse grids to compute the output distributions. See Section 4.6.1. Table II. Cost of computationally intensive tasks.

Forward Propagation

RESULTS

Uncertainties in global circulation hemodynamic indicators

The agreement between ρ(f s ) and ρ(p av ) resulting from the preoperative model and those selected in Section 3.4 should be discussed first, as this step determines the sets of resistances used for all other models. Figure 9 shows good agreement between the desired and computed distributions, in particular for the average pulmonary pressure, comfirming the trend discussed in Figure 4. While the accuracy of the Kriging surrogate could be further improved by iteratively training with additional realizations {β j , j = 1, . . . , n g } and associated pre-operative simulations results, Figure 4c shows that the maximum accuracy in the RPA/PA flow split ratio is directly related to the instrinsic unsteadiness in the simulations. In the lower range of admissible flow split ratios the variability of this result (quantified in time using a 90% confidence interval) is equal to 0.03 -0.04. Therefore, in this case, iterative refinements of the Kriging approximant will likely not lead to better agreement between target and simulated pre-operative distributions of f s and p av . Convergence of ρ(f s ) and ρ(p av ) to the target distributions for increasing sparse grid orders is also shown in Figure 9. Changes in f s , p av and δp for all models A, B, C, D are illustrated in Figure 10. The RPA flow split ratio decreases as a consequence of virtual hemi-Fontan surgery. This is mainly explained by the change in the model geometry and marginally related to the choice of boundary conditions downstream of the pulmonary artery branches. A change from imposed inlet flow and resistance boundary conditions to a coupled closed-loop LPN model with pulmonary morphometry-based RCR blocks does not produce a noticeable difference in this quantity. Augmentation of the LPA stenosis leads to an increase in the mean LPA flow split ratio, as expected, also producing an increase in the variance.

The average pulmonary pressure is not strongly affected by the local geometry of the shunt/SVC-PA junction. The role played by the boundary LPN model is more relevant in this case, as this prescribes pressure-flow relationships at all inlets and outlets. The smaller variance in average pressure for models C and D confirms that changes in pulmonary RCR blocks produce a smaller effect if embedded in a global circulation model.

The mean pressure drop shows a significant reduction with Stage II surgery, after removal of the high dissipation in the region downstream of the BT shunt. A further reduction in mean pressure drop with model C is observed, which can be explained by the smaller inlet flow associated with a coupled LPN circulation. Finally, by applying arterioplasty to the LPA, δp is greatly reduced. The small variance associated with this situation also highlights the robustness of this prediction. III.

An alternative visualization of the results associated to f s , p av and δp is shown in Figure 11, where changes in the 3D model power losses are also included.

Uncertainties in local hemodynamics

Before presenting the uncertainty in local hemodynamic results, including pressures, velocities and wall shear stresses, we briefly address the accuracy associated with computing their statistical moments. Figure 12 shows average values and standard deviations of the above quantities computed using plain Monte Carlo from all results or resulting from ( 16) for sparse grid orders 3 and 5. The three approaches generate similar results, suggesting that the variation of local quantities can be well captured with a relatively small integration order. Large average pressures at the anastomosis between the shunt and PAs are observed in Figure 12 and result from the high velocity jet exiting from the shunt and impinging on the PA vessel walls. Also, changes in boundary resistance induced by clinical variability cause significant pressure changes in this region, producing standard deviations higher than 2.0 mmHg. Post-operatively, average pressures are mainly affected by the transition from open-loop to closed-loop boundary conditions, leading to both a lower average pressure and variance. The mean pressure drop across the LPA for the various models is consistent with the global results illustrated in Figure 10.

The high dissipation region downstream of the shunt-PA anastomosis is clearly observed by examining the velocity field in Figure 14. The high δp observed in the pre-operative model is therefore only partially caused by the distal LPA stenosis. The accuracy in capturing the geometry of this region may also help to explain the differences in pressure drop between clinical observations and pre-operative model results. Wall shear stresses are extremely high in the region downstream of the shunt and greatly reduced after hemi-Fontan surgery, as shown in Figure 15. Differences in SVC inflow are responsible for the observed decrease in WSS comparing model B to C. Finally, wall shear stress is significantly reduced with arterioplasty, as expected.

DISCUSSION

We have presented a general framework (see discussion in Section 4.1) for propagation of uncertainties in clinical data to uncertainties in model parameters, and ultimately simulation output predictions. The proposed strategy allows for reporting of confidence intervals on the results, which can aid in clinical decision making and quantify trust in simulation outputs. Additionally, the quantification of uncertainty in the results induced by variability in the boundary conditions enables direct comparison of modeling practices.

While some have criticized the use of LPN models due to the large number of parameters required, this work shows that the number of parameters itself is not responsible for increasing the overall uncertainty in the results. In practice, LPN parameters are chosen to match clinical data following a set of rules, and therefore should not be assumed as independent variables. On the other hand, the inclusion of the most important parameters (those significantly affecting the hemodynamic response) is a prerequisite to properly estimating their marginal contributions. This improved perspective on cardio-vascular (CV) simulation results is achieved at the expense of an increase in the overall computational cost, as shown in Table II. In our case, determination of the initial 0D semi-empirical surrogate for MCMC required 130 solutions of the pre-operative steady state model. In the uncertainty propagation stage, the pre-operative model was solved 65 times subject to various sets of boundary resistances to check the expected resulting distributions. Finally, 65 solutions of the other three models were required for the selected maximum sparse grid order. While this is an increase over the deterministic case, further improvements in efficiency may be possible in the future and it is clear that statistical results are a necessity for future clinical adoption of these tools.

Changes in the confidence associated with systemic hemodynamic indicators are illustrated in Table III. Stage II surgery with arterioplasty produces a noticeable average increase in LPA flow. However, the removal of the stenosis increases the sensitivity of the flow split ratio due to a variation in outlet resistances. Thus uncertainty of our simulation predictions is increased following LPA plasty, making a clinical determination from this result less certain. The opposite trend is observed for the pressure difference between RPA and LPA whose estimate becomes extremely precise after post-operative augmentation of the LPA.

In addition to reporting confidence intervals, UQ assessment of confidence in simulation results may also shed light on why certain stages of single-ventricle surgery seem to be associated with large physiological variability. In our study, a large variance in pressures and velocities at the shunt-PA junction for Stage I is observed; this observation is in agreement with the well known clinical observation that Stage I physiology is particularly delicate. However, a deterministic inlet flow rate was assumed at the shunt, leaving a comprehensive exploration of the variability of systemic hemodynamic quantities to future study.

The clinical estimation of δp, distributed around an average value of 1.0 mmHg, does not seem to agree well with pre-operative simulation results. This can be seen both in Figure 10 where average RPA/LPA pressure drop was observed around 4.93 ± 0.66 mmHg and in Figure 16, where values of δp are reported from the 130 pre-operative simulations with boundary resistances determined through semi-empirical sampling. Only simulations with right pulmonary flow splits and average pressures compatible with clinically observed ρ(f s ) and ρ(p av ) are considered in Figure 16. The resulting δp are contained in the interval [3.67, 6.10], with lower values obtained for higher flow splits. One possible explanation is that values of clinical flow splits determined through PC-MRI may have been underestimated. This is likely due to both the high velocities and associated variances (see Section 5.2) at the shunt-PA junction that may produce signal loss during clinical acquisition. A further explanation relates to the uncertainty in model geometry due to the reasons expressed in Section 3.3. Peak velocities of 43 and 65 cm/s were observed at Stage I for the RPA and LPA, respectively. These velocities are smaller than those obtained through CFD simulation, consistent with the above observations. Note that the proposed framework allows for stronger statements on the agreement between simulated and measured data than those possible using deterministic results.

Table III summarizes our findings. The average values for the three quantities of interest decrease from model A through D, and the change in boundary conditions seems to significantly affect mPAP and δp across the LPA, but has limited consequence on the RPA/PA flow split ratio. Uncertainties in RPA/PA flow split increase with LPA plasty, while smaller standard deviations result for mean and pressure drop estimates.

Change in the presented results due to a different solver formulation (e.g., finite volume or finite differences instead of finite elements) is, in our opinion, mitigated by two factors. First, it is important to perform preliminary mesh adaption/refinement based on error indicators to minimize discretization errors independently of the underlying numerics. This, together with the convergence properties of the chosen numerical approach should guarantee satisfactory agreement between different solvers (see, e.g., [START_REF] Arbia | Numerical blood flow simulation in surgical corrections: what do we need for an accurate analysis[END_REF]). Second, distributions of resistance boundary conditions are determined (using Kriging plus MCMC) in order to match the expected clinical targets, independently of the model formulation. This should shift the differences between alternative modeling approaches to differences in the adopted boundary conditions, leaving the final predictions marginally affected.

We recognize several limitations to our study, mostly related to various sources of uncertainty that should be added to the proposed framework. First, the effects of geometrical model uncertainty need to be investigated. This is also particularly relevant in virtual surgery applications, where the final post-operative anatomy may be differently implemented by different surgeons and therefore subject to additional uncertainty. We made a choice to focus this paper on uncertainty related to the clinical acquisition of hemodynamics data and how this propagates to the model predictions. Deformable wall simulations (including uncertainty on material properties) will also provide improved accuracy. While we recognize that considering all possible sources of uncertainty is a very challenging task, inclusion of contributions from segmented geometry and material properties will be an important step forward to understand how the uncertainty in these factors affects the outputs.

Second, computations have been presented for coupled 0D-3D models where only variations in equivalent pulmonary RCR blocks were considered. Random inputs for other LPN parameters (resistances, capacitances and inductances) should also be inferred from the available clinical data, as well as all assumptions used to assign relative values of LPN parameters. Third, the proposed framework should be applied to multiple patients to better characterize inter-patient variability of global and local hemodynamic statistics. Finally, physiologic variation is a formidable source of

CONCLUSION

A combined inverse and forward uncertainty propagation approach is presented in this study to quantify the effects of uncertainty in clinically estimated RPA/PA flow split ratio and mean pulmonary pressure on the results of stage I and stage II single-ventricle simulations. The proposed approach combines an efficient inverse modeling approach to determine the distribution of model parameters arising from distributions in clinical data, and then propagates these uncertainties to determine confidence intervals and statistics on simulation predictions. To the best of the author's knowledge this study is the first to quantify the uncertainty in virtual surgery predictions directly from estimates of uncertainty in patient-specific clinical data and proposes a new paradigm of presenting CV simulation results that we hope will challenge the common practice of providing only deterministic results.

On a concluding note about the applicability of the proposed UQ approach to aid decision making in the clinic, while we feel that a full UQ approach to cardiovascular simulation is a possibility today, due to the continuous increase in available computational resources, studies on larger populations are also possible to establish guidelines on which quantities are consistently associated with accurate predictions in virtual surgery and those that are less to be trusted on average.

Figure 1 .

 1 Figure 1. (a) Pre-and post-operative model configurations. Model A is a stage I pre-operative model with resistance boundary conditions (BC) and prescribed shunt inflow. Model B is a stage II post-operative model with resistance BC and prescribed SVC inflow. Model C is a stage II post-operative model (same geometry as model B) with closed loop BCs. Model D is a stage II post-operative model with closed loop BCs and augmented LPA stenosis. (b) Volumetric flow rate through SVC/shunt for the four models. An error bar of one standard deviation is shown for pulsatile simulations (models C and D) indicating the variation in the mean flow rate due to uncertainty in the pulmonary RCR blocks.

8 D

 8 . E. SCHIAVAZZI ET AL. DOI:10.1002/CNM.2737

10 DFigure 2 .

 102 Figure 2. Selected distributions of (a) RPA/PA flow split ratio ρ(fs) and (b) mean pulmonary arterial pressure (mPAP) ρ(pav).

1 Figure 3 .

 13 Figure 3. Two-dimensional quadrature grids suited for (left) Gaussian[START_REF] Genz | Fully symmetric interpolatory rules for multiple integrals over infinite regions with Gaussian weight[END_REF] and (right) uniform[START_REF] Petras | Smolyak cubature of given polynomial degree with few nodes for increasing dimension[END_REF] probability measures.

Figure 4 .

 4 Figure 4. Discrepancy between values of (a) flow split ratios and (b) average pressures used as inputs of the semi-empirical approximation and obtained as a result of the pre-operative model. (c) The degree of unsteadiness is represented by plotting the 90% confidence interval vs. the average RPA/PA flow split ratio.

  a) Left vs. right average outlet pressures for 130 pre-operative simulations with boundary resistances determined through the proposed semi-empirical 0D approach. (b) Pre-operative flow rates and pressures at LPA branch outlets.

Figure 6 .

 6 Figure 6. Graphical representation of the two-dimensional Kriging surrogates for right flow split ratio (left), mean pulmonary pressure (centre) and RPA/LPA pressure discrepancy (right).

Figure 7 .

 7 Figure 7. Results of the MCMC simulations. (a) Traces of β k . The picture shows a good mixing of samples for the two distributions. (b) Plots of ergodic average and standard deviation for β 1 and β 2 . (c) Resulting marginal PDFs of β 1 and β 2 .

Figure 8 .

 8 Figure 8. General outline of the various tasks proposed and associated section number in the manuscript. Specific to the proposed application to single-ventricle pathologies are the relevant clinical data and the type of boundary conditions (BC) for the 3D CFD simulations.

Figure 9 .

 9 Figure 9. Distributions of fs (top left) and pav (top right) resulting from the preoperative model using resistance distributions identified through inverse Bayesian estimation. The convergence in the first two moments for increasing sparse grid orders is also shown (bottom).

Figure 10 .

 10 Figure 10. Distributions of fs, pav and δp for all pre-and post-operative models. Target distributions for preoperative fs and pav are also indicated using dashed lines. For clarity, only the first two statistical moments are represented. Values of 90% confidence intervals are also reported in TableIII.

Figure 11 .

 11 Figure 11. Alternative visualization of changes in fs, pav, δp and power loss for models A,B,C,D. Error bars are representative of the variability produced by one standard deviation.

Figure 12 .

 12 Figure 12. Statistics of local hemodynamic results in terms of pressures, velocities and wall shear stresses resulting from straightforward Monte Carlo average or numerical integration using increasing sparse grid orders.

Figure 13 .

 13 Figure 13. (a) Average values and (b) standard deviations of local pressure results (in mmHg). Contours are limited to 16.0 mmHg and 2.0 mmHg for maximum average pressure and standard deviation, respectively.

Figure 14 .

 14 Figure 14. (a) Average values and (b) standard deviations of blood velocity results (in mm/s).

Figure 15 .

 15 Figure 15. (a) Average values and (b) standard deviations of wall shear stress results (in Pa). The maximum value for the average and standard deviation of the wall shear stress is set to 20 mmHg and 15 mmHg, repsectively.

Figure 16 .

 16 Figure 16. Distributions of pre-operative model results in terms of (left) pav (centre) fs and (right) δp as obtained from 130 simulations with boundary resistances determined through semi-empirical sampling. The range in RPA to LPA pressure differences appears in the range [3.67, 6.10].

Table III .

 III 90% Confidence intervals, averages, standard deviations and pre-op to post-op percentage change for pre-and post-operative model results.

Copyright c 2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Biomed. Engng. (2010) Prepared using cnmauth.cls DOI: 10.1002/cnm

8. ACKNOWLEDGEMENTS Comments and feedback from the anonymous Reviewers are gratefully acknowledged. This work was supported by the Leducq Foundation as part of a Transatlantic Network of Excellence for Cardiovascular Research, a Burroughs Wellcome Fund Career award at the Scientific Interface, NSF CAREER OCI-1150184, a British Heart Foundation Clinical Research Fellowship FS/12/35/29566 and used computational resources from the Extreme Science and Engineering Discovery Environment (XSEDE), supported by National Science Foundation grant number ACI-1053575. We also acknowledge the open source SimVascular project at simtk.org.

uncertainty for cardiovascular simulation predictions and further studies are warranted to better quantify its contribution.