Multilevel and multifidelity uncertainty quantification for cardiovascular hemodynamics
Casey M. Fleeter, Gianluca Geraci, Daniele E. Schiavazzi, Andrew M. Kahn, Alison L. Marsden
MMultilevel and multifidelity uncertainty quantification for cardiovascularhemodynamics (cid:73) , (cid:73)(cid:73) Casey M. Fleeter a , Gianluca Geraci b , Daniele E. Schiavazzi c , Andrew M. Kahn d , Alison L. Marsden a,e, ∗ a Institute for Computational and Mathematical Engineering, Stanford University, Stanford, CA, USA b Center for Computing Research, Sandia National Laboratories, Albuquerque, NM, USA c Department of Applied and Computational Mathematics and Statistics, University of Notre Dame, Notre Dame, IN, USA d Department of Medicine, University of California San Diego, La Jolla, CA, USA e Departments of Pediatrics and Bioengineering, Stanford University, Stanford, CA, USA
Abstract
Standard approaches for uncertainty quantification in cardiovascular modeling pose challenges due to thelarge number of uncertain inputs and the significant computational cost of realistic three-dimensional sim-ulations. We propose an efficient uncertainty quantification framework utilizing a multilevel multifidelityMonte Carlo (MLMF) estimator to improve the accuracy of hemodynamic quantities of interest while main-taining reasonable computational cost. This is achieved by leveraging three cardiovascular model fidelities,each with varying spatial resolution to rigorously quantify the variability in hemodynamic outputs. Weemploy two low-fidelity models (zero- and one-dimensional) to construct several different estimators. Ourgoal is to investigate and compare the efficiency of estimators built from combinations of these two low-fidelity model alternatives and our high-fidelity three-dimensional models. We demonstrate this frameworkon healthy and diseased models of aortic and coronary anatomy, including uncertainties in material propertyand boundary condition parameters. Our goal is to demonstrate that for this application it is possible toaccelerate the convergence of the estimators by utilizing a MLMF paradigm. Therefore, we compare ourapproach to single fidelity Monte Carlo estimators and to a multilevel Monte Carlo approach based onlyon three-dimensional simulations, but leveraging multiple spatial resolutions. We demonstrate significant,on the order of 10 to 100 times, reduction in total computational cost with the MLMF estimators. Wealso examine the differing properties of the MLMF estimators in healthy versus diseased models, as well asglobal versus local quantities of interest. As expected, global quantities such as outlet pressure and flowshow larger reductions than local quantities, such as those relating to wall shear stress, as the latter relymore heavily on the highest fidelity model evaluations. Similarly, healthy models show larger reductionsthan diseased models. In all cases, our workflow coupling Dakota’s MLMF estimators with the SimVascularcardiovascular modeling framework makes uncertainty quantification feasible for constrained computationalbudgets.
Keywords:
Cardiovascular modeling, Uncertainty quantification, Multilevel Monte Carlo, MultifidelityMonte Carlo, Multilevel multifidelity Monte Carlo
1. Introduction
Cardiovascular disease is the leading cause of mortality worldwide, as well as the leading cause of deathfor both men and women in the United States [1, 2]. Computational models of the cardiovascular system are (cid:73) https://doi.org/10.1016/j.cma.2020.113030 (cid:73)(cid:73) c (cid:13) ∗ Corresponding author
Email address: [email protected] (Alison L. Marsden)
Accepted Manuscript submitted to Computer Methods in Applied Mechanics and Engineering April 20, 2020 a r X i v : . [ q - b i o . Q M ] A p r ncreasingly adopted due to the hemodynamic insights they provide, which can be beneficial in diagnosis,assessment of disease progression risk, and treatment planning for cardiovascular disease. These models canbe generated in a non-invasive fashion using routinely collected clinical data [3], and their applications covera wide range of diseases and anatomies. Coronary artery disease is the most prevalent cause of death in theUnited States [4], and has thus been the subject of many modeling studies. These include computing frac-tional flow reserve for patients with coronary stenoses [5] and assessing differences in the hemodynamic andmechanical response of arterial and venous grafts towards determining cause of vein graft failure [6]. Com-putational models have also been applied to congenital heart disease for prognosis and treatment planning,from testing designs for surgical interventions of single ventricle congenital heart disease [7–9] to thromboticrisk assessment for Kawasaki disease [10–12], and assessment of disease progression and mechanical stimuliin pulmonary hypertension [13, 14]. Hemodynamic studies are also used for additional disease cases, such asanalyzing blood flow patterns in abdominal aortic aneurysms [15, 16] and predicting the development andrupture of cerebral aneurysms [17–19].Realistic three-dimensional hemodynamic simulations require segmentation of vascular anatomies frommedical image data, followed by numerical solution of the equations governing blood flow in elasticallydeformable vessels. This is typically achieved through a complex and expensive workflow. However, varioussimplifying assumptions can be made to generate models of intermediate complexity. One-dimensional hemodynamic models are formulated by integrating the Navier-Stokes equations across the vessel crosssections, with additional assumptions on the properties of the fluid and the vascular walls [20–22]. Alinearization of the incompressible Navier-Stokes equations around rest conditions leads to an even simpler zero-dimensional (or lumped parameter) formulation where vascular networks are analyzed using analogouselectrical circuits [23–26]. These simplified representations are computationally inexpensive, especially whencompared to the cost of a full three-dimensional model.Unfortunately, all models of the cardiovascular system, as well as biological and biomedical systemsmore generally, are intrinsically affected by model and data uncertainty [27]. These uncertainties arise frommodeling choices, measurement errors, and intra- and inter-patient variability. Therefore, relying on a solelydeterministic framework provides limited information, as results cannot include statistical distributions orconfidence intervals, and hinders their application in the clinical setting. Here, we aim to account forthese uncertainties and provide confidence intervals through uncertainty quantification (UQ). By adoptinga stochastic framework, model parameters are sampled from appropriate probability distributions which areeither assumed, relying on existing literature and clinical data, or assimilated from available patient-specificdata.Recently, UQ has gained traction in the field of cardiovascular modeling, primarily utilizing UQ tech-niques centered around a single model complexity, primarily three-dimensional computational fluid dynamics(CFD) simulations or lower complexity one-dimensional models. Recent studies have investigated the im-pact of geometry and boundary condition parameters on coronary fractional flow reserve [28], demonstrateda multi-resolution approach to quantify boundary condition uncertainties [29] and in conjugation with ran-dom fields [30] for coronary artery bypass grafts, and performed stochastic collocation in a one-dimensionalarterial network [31], along with others [32–39].Several challenges arise when quantifying uncertainty in cardiovascular simulations. First, there aretypically multiple sources of uncertainty to account for and propagate through the model. Such uncertaintiesinclude, but are not limited to, boundary condition parameters, spatial variability of material properties invascular tissue, and operator dependent vessel lumen segmentation leading to an uncertain vascular modelanatomy. With myriad sources of uncertainty, it follows that any UQ scheme must account for a largenumber of possibly heterogeneous random inputs; this leads to the so-called curse of dimensionality , whichposes significant challenges for generalized polynomial chaos expansions [40] and stochastic collocation [41]due to the fast increase of the computational complexity in high-dimensions. Second, for any realizationof uncertain input parameters, a three-dimensional model solution is computationally expensive. Three-dimensional simulations involve discrete representations with millions of degrees of freedom that are solved inparallel over multiple computational nodes. Maintaining a reasonable computational cost therefore becomeschallenging when relying solely on three-dimensional simulations.Monte Carlo estimators [42] have, in this context, many desirable properties. They are unbiased, offer2exibility with respect to heterogeneous input sources, and the associated variance depends only on the truevariance and number of model evaluations, not on the problem dimensionality or regularity. However, onecan typically afford only a small number of these evaluations when working with expensive deterministicsolvers. Multilevel and multifidelity estimators lead to a reduced variance compared to their “vanilla”Monte Carlo analogues for the same computational cost. The interested reader is referred to [43, 44]for an in depth discussion of these methods. Another alternative to multilevel estimators, which rely onspatial resolution, and multifidelity estimators, which are based on a control variate approach, are so-calledmultilevel multifidelity (MLMF) methods [45–48] which combine the efficiency of multilevel Monte Carloacross discretization levels with embedded control variates based on a low-fidelity model at each level.In this way, MLMF estimators leverage a cascade of varying-fidelity models ranging from computationallyinexpensive, patient-agnostic lumped parameter models to computationally intensive, patient-specific modelsof blood flow. We provide an overview of these approaches in section 3.The main contribution of this paper is to introduce the integration of multilevel multifidelity estimatorsfor uncertainty quantification with cardiovascular models. This framework integrates and manages multiplehemodynamic solvers from the SimVascular open source package [49] with Sandia National Laboratories’Dakota toolkit [50] for uncertainty analysis in a manner designed to automatically compute MLMF estima-tors for a variety of quantities of interest (QoIs). Due to the much lower computational cost of one- andzero-dimensional models, the MLMF estimators allow stochastic analysis to be performed at a drasticallyreduced computational cost. The specific contributions of this paper are1. Definition of model fidelities (3D, 1D and 0D) with corresponding discretization levels obtained aftera thorough convergence study.2. Creation of a semi-automated framework enabling Dakota to directly manage and invoke the SimVas-cular modeling simulations.3. Numerical experiments utilizing the above framework to demonstrate the performance of the estimatorsfor cardiovascular hemodynamics on pulsatile simulations with realistic boundary conditions, healthyand diseased models, and a large number of quantities of interest not previously addressed in theliterature.This paper is organized as follows. Sections 2 and 3 provide background on the cardiovascular modelingand uncertainty quantification techniques used in the study, respectively. Section 4 describes the workflow,including further details on model construction as well as the selection of uncertain parameters. Section 5provides UQ results for four patient anatomies, providing comparisons of the performance of uncertaintyquantification estimators for healthy versus diseased models and local versus global QoIs, as well as com-parisons of different uncertainty quantification estimators. Sections 6 and 7 provide discussions and drawconclusions.
2. Cardiovascular Modeling
We first illustrate the three model formulations in decreasing order of complexity (3D, 1D and 0D) anddiscuss their trade-offs (Figure 1). While zero-dimensional formulations are computationally inexpensive,these models provide only a coarse spatial representation of the quantities of interest. One-dimensionalmodels, while somewhat more expensive, provide finer grained resolution for QoIs along the vessel length,but not across vessel cross-sections. Finally, while the most computationally expensive by far, the three-dimensional CFD models are well-resolved even for local flow features, allowing detailed examination ofquantities of interest such as wall shear stress (WSS) in specific regions of interest.
Three-dimensional anatomical models are typically constructed from medical image data of specificpatients. In our workflow, these images are imported into the SimVascular open source platform [49], whichprovides a pipeline for model creation, generation of tetrahedral meshes, application of physiologic boundaryconditions, and finite element solution. More specifically, the modeling pipeline is as follows. First, the user3 oarse Resolution Fine Resolution C o m pu t a t i o n a l C o s t
1D 3D0D Figure 1: Improved fine-scale resolution for quantities of interest comes at the expense of higher computational cost whensolving the model hemodynamics with different fidelity methods. obtains patient-specific medical image data. Next, centerline paths are generated for the vessels of interest.The lumen is then segmented along, but perpendicular to, the centerline path in a semi-automated fashion.Next, the model is lofted using boolean operations to merge the vessels together. Meshing then occurs,with many specific options for local mesh refinement and model smoothing. At this point, the model isready for the application of boundary conditions, which range from prescribed flow or pressure values toa geometrically multi-scale approach [22, 51] where zero-dimensional lumped parameter networks are usedto model the upstream and downstream circulation [23, 52]. Finally, the model is solved with a stabilizedincompressible Navier-Stokes finite element solver [53, 54] with linear tetrahedral elements and second orderimplicit generalized- α time integration [55]. A custom linear solver with specialized preconditioning increasessolver efficiency [56, 57]. SimVascular post-processing tools are then used to extract hemodynamic indicatorsof interest, such as pressures and wall shear stress, known to be correlated in the literature with endothelialdamage and thrombus formation [58]. Hemodynamic simulations can be run with rigid or deformable walls,i.e., with or without accounting for the mutual interaction between fluid and structure. In this study thecoupled momentum method (CMM) [59], which assumes a linear elastic constitutive model characterized bymembrane and shear stiffness, is selected to model fluid-structure interaction (FSI). By using a monolithiccoupling approach, CMM FSI results in a monolithic algorithm in which the solid momentum contributionsare embedded, using the same degrees of freedom, in the fluid equations. The CMM formulation is availablein the open source solver package svSolver, distributed as part of SimVascular.Following the derivation in [59], blood flow is approximated as flow of incompressible Newtonian fluidin a domain D ⊂ R . The boundary Γ of domain D is considered to be divided into three partitionsΓ = ∂D = ∂D g ∪ ∂D h ∪ ∂D s , where Dirichlet, Neumann, and wall boundary conditions are applied,respectively. The mesh on D consists of n el linear tetrahedral elements D e , e = 1 , . . . , n el . To preventrigid body motion, the perimeters of all inlet and outlet vessels (denoted as Γ s ) are mechanically fixed byapplying zero relative displacement and zero velocity boundary conditions. Therefore, the combined weak4orm for the coupling of the stabilized incompressible Navier-Stokes equations and the CMM FSI methodcan be written as (cid:90) D (cid:104) (cid:126)w · ( ρ ∂(cid:126)v∂t + ρ(cid:126)v · ∇ (cid:126)v − (cid:126)f ) + ∇ (cid:126)w : ( − pI + τ ) − ∇ q · (cid:126)v (cid:105) d(cid:126)x − (cid:90) ∂D h (cid:126)w · (cid:126)hds + (cid:90) ∂D h qv n ds + (cid:90) ∂D g qv n ds + n el (cid:88) e =1 (cid:90) D e ∇ q · τ M ρ (cid:126) L ( (cid:126)v, p ) d(cid:126)x + n el (cid:88) e =1 (cid:90) D e (cid:104) ( (cid:126)v · ∇ ) (cid:126)w · τ M (cid:126) L ( (cid:126)v, p ) + ∇ · (cid:126)wτ C ∇ · (cid:126)v (cid:105) d(cid:126)x + n el (cid:88) e =1 (cid:90) D e (cid:104) (cid:126)w · ( ρ ˆ (cid:126)v · ∇ (cid:126)v ) + ( (cid:126) L ( (cid:126)v, p ) · ∇ (cid:126)w )( τ (cid:126) L ( (cid:126)v, p ) · (cid:126)v ) (cid:105) d(cid:126)x + ζ (cid:90) ∂D s (cid:104) (cid:126)w · ρ s ∂(cid:126)v∂t + ∇ (cid:126)w : σ s ( (cid:126)u ) (cid:105) ds − ζ (cid:90) Γ s (cid:126)w · (cid:126)h s dl + (cid:90) ∂D s qv n ds = 0 , (1)where (cid:126)v is velocity, p is pressure, ρ is blood density, τ = µ ( ∇ (cid:126)v + ( ∇ (cid:126)v ) T ) is the viscous stress tensor where µ is the dynamic viscosity, and ζ and σ s are the prescribed vessel wall thickness and stress tensor, respectively.Vectors and matrices are denoted by (cid:126) · and · , respectively. These equations lend themselves to the choice ofuncertain parameters we will discuss in subsection 4.5. One dimensional formulations for cardiovascular hemodynamics result from a long history of researchand experiments, whose origins can be traced back to Euler’s work on blood flow in elastic arteries [60]. Forour purposes, an in-house stabilized finite element solver is used for the one-dimensional hemodynamics. Itsspecific formulation is adapted from Hughes and Lubliner [20], with implementation details discussed in [61–63]. Blood is simulated as a Newtonian fluid, assuming velocity in the axial direction ( z ) of each cylindricalbranch, constant pressure over the vessel cross section, and a non-slip boundary condition applied at the ves-sel lumen. The governing equations are obtained by integrating the incompressible Navier-Stokes equationsover the cross section of a deformable cylindrical domain. A constitutive model relates pressure to changein cross-sectional area. The resulting conservation equations for mass and momentum are, respectively [20], ∂A∂t + ∂Q∂z = − ψ∂Q∂t + ∂∂z (cid:20) (1 + δ ) Q A (cid:21) + Aρ ∂p∂z = Af + N QA + v ∂ Q∂z . (2)The solution variables are the cross-sectional area A and flow rate Q , and model parameters are blooddensity ρ , external force f , kinematic viscosity ν , and an outflow function ψ for permeable vessels. Toprescribe a velocity profile over the cross-section, parameters δ and N are defined as δ = 1 A (cid:90) A ( φ − ds and N = ν (cid:90) ∂A ∂φ∂m dl. (3)Pressure continuity and conservation of mass are enforced at the bifurcations via Lagrange multipliers. Asingle branch point is presented without loss of generality, as this case easily generalizes to multiple branchpoints. For the case of one branch point l with m inlets and n outlets, the following m + n constraints areenforced: m (cid:88) C =1 Q in C − n (cid:88) C =1 Q out C = 0 p in C − p in l = 0 C = 2 . . . mp out C − p in l = 0 C = 1 . . . n (4)5he first equation enforces conservation of mass. The second equation enforces pressure continuity amongthe m inlet vessels, with one vessel chosen without loss of generality as a reference vessel. The pressurecontinuity is enforced between this reference branch and each of the m − n outlet vessels, with one constraintfor each outlet vessel. With Lagrange multipliers, these constraints are formulated as a potential function Z as Z = λ Q (cid:34) m (cid:88) C =1 Q in C − n (cid:88) C =1 Q out C (cid:35) + m (cid:88) C =2 (cid:2) λ p ( C − ( p in C − p in l ) (cid:3) + n (cid:88) C =1 (cid:2) λ p ( C − m ) ( p out C − p in l ) (cid:3) . (5)The details of the numerical implementation of the Lagrange multipliers and potential function are availablein [63]. Branch losses are not included here, though empirical relations for branch losses have been consideredin prior work and could be easily included in future studies. Such minor loss values at branch junctions areimplemented and explained in [64].Finally, a constitutive relationship is needed to complete the system of equations, relating the pressure p to the cross-sectional area: p ( z, t ) = ¯ p ( A ( z, t ) , z, t ) . (6)Specifically, we use a linear elastic model, consistent with the CMM formulation discussed in the next section,given by ¯ p ( A ( z, t ) , z, t ) = p ( z ) + (cid:18) E hr ( z ) (cid:19) (cid:32)(cid:115) A ( z, t ) A ( z ) − (cid:33) , (7)where p ( z ) = p ( z,
0) is the initial undeformed pressure, E is the elastic modulus of the vascular tissue, h the wall thickness, and r ( z ) and A ( z ) are the undeformed inner radius and area, respectively.A weak (integral) formulation of these equations is solved on each segment using a stabilized finiteelement method in space and a discontinuous Galerkin method in time. The non-linear algebraic system ofequations is finally solved using modified Newton iterations. Pressure and cross-sectional area continuity areenforced at the joints between two segments through Lagrange multipliers. Further details on the numericalimplementation are provided in [63]. A lumped parameter network representation of blood flow is used as our lowest-fidelity model. This isa circuit layout formulated by hydrodynamic analogy in terms of flow rate (electrical current) and pressuredifferences (voltage), where each circuit element is associated with an algebraic or differential equation, i.e.resistor ∆ P = R Q, capacitor Q = C dPdt , inductor ∆ P = L dQdt . (8)These equations are assembled in a system of ODEs which can be efficiently solved, for example, with a 4 th –order Runge-Kutta scheme. Resistors are used to represent viscous dissipation on the vessel walls, capacitorsare used to represent vascular tissue compliance, and inductors are used to represent the inertia of blood.A Poiseuille flow assumption is used to determine the model parameters for these circuit elements [65], i.e. R = 8 µ lπ r , C = 3 l π r E h , L = l ρπ r , (9)where µ is the dynamic viscosity of blood, l the vessel length, r the vessel radius, E the elastic modulus, h the wall thickness and ρ the density of blood. These models have been used extensively to model theheart and circulatory system, achieving remarkably realistic flow and pressure waveforms and allowing forthe assessment of physiologic changes to treatments [7, 35, 66, 67].6 .4. Model fidelity validation and generation To take full advantage of multiple fidelities in uncertainty propagation, generation of lower fidelity modelsfrom existing three-dimensional models should be accomplished with little effort. We therefore implementeda conversion tool which automatically extracts the one-dimensional model geometry from the centerlinesor two-dimensional segmentations constructed as part of the three-dimensional modeling process, describedin subsection 2.1. This allows one-dimensional models to be obtained at a minimal cost once a three-dimensional model has been built. The same inlet and outlet boundary conditions can be applied to theone- and three-dimensional models.To verify the agreement between low fidelity and three-dimensional model results, we performed a pre-liminary validation study on the aorto-femoral model. Results from all three model formulations werecompared under steady and pulsatile flow conditions and with resistance and RCR boundary conditionsat the model outlets. One- and zero-dimensional models produce very similar results as they both can beformulated to include vessel compliance. Good agreement is observed between the one-dimensional resultsand three-dimensional simulations with deformable (compliant) vessel walls. Agreement in selected flow andpressure quantities of interest for pulsatile inlet flow and resistance boundary conditions is acceptable forall three model fidelities (Figure 2). The agreement between models leads to high correlations between themodel fidelities, which is desirable for a multifidelity approach. .
00 0 .
25 0 .
50 0 . P r e ss u r e ( mm H g ) .
00 0 .
25 0 .
50 0 . F l o w ( c m / s ) .
00 0 .
25 0 .
50 0 . F l o w ( c m / s ) .
00 0 .
25 0 .
50 0 . P r e ss u r e ( mm H g ) (a) (c)(b) (d)
0D 1D 3D
Pressures Flows
Figure 2: Agreement for outlet flow and pressure quantities of interest is seen between 3D, 1D, and 0D models with resistanceboundary conditions and vessel compliance built into the model. (a) Pressure at superior mesenteric artery. (b) Pressure atright iliac artery. (c) Flow at right renal artery. (d) Flow at left internal iliac artery. . Uncertainty Quantification In this section we briefly explain multilevel, multifidelity control variate, and multilevel multifidelityMonte Carlo methods for uncertainty quantification (schematics in Figure 3). The main goal of thesemethods is to leverage computationally inexpensive, but potentially less accurate, models to decrease thevariance of our estimators for a variety of quantities of interest (QoIs) while using a limited number ofhighest-fidelity model evaluations. (b) (c)(a) LF Fine - LF Medium DiscrepancyLF Medium - LF Coarse DiscrepancyLF CoarseCVCVCVHF Fine - HF Medium DiscrepancyHF Medium - HF Coarse DiscrepancyHF CoarseHF M - HF M-1
DiscrepancyHF
M-1 - HF
M-2
DiscrepancyHF HF LFCV
Fidelity L e v e l Figure 3: Multilevel multifidelity simulations are comprised of (a) a multifidelity control variate approach coupling high andlow fidelity models combined with (b) a multilevel approach incorporating varying mesh levels of the high fidelity model in atelescopic sum of discrepancies between each level. (c) In this study, the specific models we employ are three mesh levels of the3D model coupled by control variates with either three mesh levels of the 1D model or a fine mesh resolution 1D model andtwo different zero-dimensional models, a “full” model with multiple circuit elements and a “simple” model with only resistorelements.
Consider a complete probability space (Ω , F , P ) where Ω is a set of elementary events, F is a Borel σ -algebra of subsets in Ω, and P a probability measure with values in [0 ,
1] over events in F . A vectorof random variables ξ = ( ξ , ξ , . . . , ξ d ) ∈ R d has components ξ i : Ω → Σ ξ i , i = 1 , . . . , d with marginals ξ i ∼ ρ i ( ξ i ) and joint probability density ρ ( ξ ). These variables represent the sources of uncertainty in themodel, while x ∈ D ⊂ R and t ∈ R + are the spatial and temporal variables, respectively.We are interested in the statistical characterization of the quantity of interest Q = Q ( ξ ), which isobtained as a functional of the numerical solution in time and space. Specifically we look at the first twomoments, the expected value E [ Q ] or the variance V [ Q ], where for a given realization ξ ( i ) of the randominputs, Q is determined through the numerical solution of the incompressible Navier-Stokes equations on adiscrete domain D ⊂ R . While in general we assume our random inputs to be arbitrarily distributed andcorrelated, for this study we only consider independent inputs.8n the remainder of this work, we will always use an estimator for the mean value of our quantity ofinterest Q , with various methods employed to reduce its variance, but the proposed methodology can beextended to the estimation of higher-order statistical moments. In general, variance reduction techniquescan be applied to estimators for other statistics pertaining to the quantity of interest Q .Due to the discrete nature of the solution scheme adopted, the accuracy of the generic output Q isaffected by the discretization level M , i.e., the number of degrees of freedom in the spatial or temporalmeshes. Therefore, a choice of discretization levels for which Q M → Q as M → ∞ guarantees that E [ Q M ] → E [ Q ] as M → ∞ . The Monte Carlo estimator based on N realizations for the expected value of Q M is defined as E [ Q M ] = (cid:90) Ω Q M ( ξ ) ρ ( ξ ) d ξ (cid:39) ˆ Q MC M,N = 1 N N (cid:88) i =1 Q ( i ) M , (10)where the i -th realization is denoted as Q ( i ) M = Q M ( ξ ( i ) ). The variance of this estimator can be computedas V [ ˆ Q MC M,N ] = 1 N V [ Q M ] , (11)whereas the mean squared error (MSE) is E [( ˆ Q MC M,N − E [ Q ]) ] = V [ ˆ Q MC M,N ] + ( E [ Q M − Q ]) . (12)We see from (12) that the MSE contains two contributions, from the estimator variance V [ ˆ Q MC M,N ] anddeterministic bias ( E [ Q M − Q ]) . The deterministic bias can be reduced with a finer discretization for Q M ,i.e., a larger value of M , while the estimator variance is reduced by increasing the number of realizations N . As a consequence, obtaining a highly accurate MC estimator (i.e. an estimator with a small MSE)would potentially require a large number N of highly accurate (large M ) numerical simulations, makingthis approach impractical for high-fidelity models. In this work, we make the following assumptions. Weperformed convergence studies (not reported here for brevity) to select the three-dimensional resolutionlevel guaranteed to obtain an acceptable discretization error for the QoIs we considered. Following thisconvergence study, a fixed set of resolution levels was constructed via a comprehensive mesh convergencestudy (details in subsection 4.3). Our goal in this work is to build an unbiased estimator with respect tothe high-fidelity model. Hence, we only target the variance contribution of the MSE. We assume the biasof the high-fidelity model to already be satisfactory; that is, we do not extend the hierarchy of models tothose with finer discretizations. However, it is also possible to extend this approach in order to target thefull MSE as suggested in the MLMC literature (e.g. see [43]).Due to the large computational cost of evaluating the output Q ( i ) M of a single hemodynamic simulation,increasing the number of samples is impractical, and multilevel, multifidelity, and multilevel multifidelityapproaches (subsections 3.2, 3.3, and 3.4, respectively) offer a better alternative. These approaches allow usto combine models of multiple fidelities and computational costs into the same stochastic framework, thusproviding a means to efficiently manage the available computational resources. The main idea of multilevel MC (MLMC) approaches is to replace the quantity of interest Q with atelescoping sum of the differences between the next coarsest resolution level (Figure 3b). For an extensivereview on MLMC estimators we refer the reader to [43].Consider L discretization (or resolution) levels, i.e., { M (cid:96) : (cid:96) = 0 , . . . , L } , where M < M < · · · < M L := M . Using linearity of expectation and a telescopic sum, we re-write the expected value of our QoI afterintroducing the discrepancy between QoIs at successive resolutions Y (cid:96) , Y (cid:96) = (cid:40) Q M if (cid:96) = 0 Q M (cid:96) − Q M (cid:96) − if 0 < (cid:96) ≤ L ⇒ E [ Q M ] = L (cid:88) (cid:96) =0 E [ Y (cid:96) ] . (13)9he multilevel Monte Carlo (MLMC) estimator for E [ Q M ] is assembled from the Monte Carlo independentestimators of E [ Y (cid:96) ]. At level (cid:96) , N (cid:96) realizations are used to estimate Q ML M , and so the expected value of ourQoI is ˆ Q ML M,N = L (cid:88) (cid:96) =0 ˆ Y (cid:96),N (cid:96) = L (cid:88) (cid:96) =0 N (cid:96) N (cid:96) (cid:88) i =1 Y ( i ) (cid:96) . (14)As before, Y ( i ) (cid:96) is evaluated at the i -th, i = 1 , . . . , N (cid:96) , realization of the stochastic vector ξ drawn fromthe distribution ρ ( ξ ). The MLMC estimator has variance equal to V [ ˆ Q ML M,N ] = L (cid:88) (cid:96) =0 N (cid:96) V [ Y (cid:96) ] . (15)The advantages of the MLMC method come from the hierarchical nature of Y (cid:96) , the difference in thequantity of interest between successive resolutions. As Q M → Q as M → ∞ , Y (cid:96) → (cid:96) increases andtherefore the contribution to the overall variance from different resolution levels decreases with (cid:96) . Thisallows us to shift the computational burden to the coarser levels, which are computationally cheaper. Extrapolation can be performed to approximate the samples needed to improve estimator variance by afactor (cid:15) . The optimal sample allocation at each level can be determined by minimizing the computationalcost of MLMC subject to a fixed target (cid:15) improvement factor [43]. The minimum total computational costis C [ ˆ Q ML M,N ] = (cid:80)
L(cid:96) =0 N (cid:96) C (cid:96) , where C (cid:96) is the computational cost of one evaluation of Y (cid:96) (note that this is thecombined cost of two evaluations on successive resolution levels for (cid:96) ≥ (cid:96) givenby N (cid:96) = 1 (cid:15) (cid:32) L (cid:88) k =0 (cid:112) V [ Y k ] C k (cid:33) (cid:115) V [ Y (cid:96) ] C (cid:96) . (16) Multifidelity (MFMC) approaches represent a flavor of the better known control variate (CV) variancereduction technique in Monte Carlo estimation (Figure 3a). Two models are now defined, i.e., a low-fidelity(LF) and a high-fidelity (HF) model, at discretization levels M LF and M , respectively. Without loss ofgenerality, in the following we only use the level M , where it is implicitly assumed that the discretizationlevels of the HF and LF model are in general independent. In this approach, the generic quantity of interest Q HF M is replaced by Q CV , HF M which embeds a correction term based on the LF model. For full details ofthe approximate CV estimators used here, we refer the reader to [68, 69]. For more details of the generalmultivariate CV approach, we refer the reader to [70].The final form of the selected approximate CV MFMC estimator is obtained by combining the MCestimators as ˆ Q CV , HF M,N HF = ˆ Q HF M,N HF + α (cid:16) ˆ Q LF M,N HF − ˆ Q LF M,N LF (cid:17) , (17)where N LF = N HF + ∆ LF = N HF (1 + r ) and the additional LF realizations ∆ LF = rN HF are drawnindependently of the initial set of N HF .The regression coefficient α and the value of the parameter r > α value α = − ρ (cid:118)(cid:117)(cid:117)(cid:116) V [ ˆ Q HF M,N HF ] V [ ˆ Q LF M,N HF ] (18)yields a minimal variance V [ ˆ Q CV , HF M,N HF ] = V [ ˆ Q HF M,N HF ] (cid:18) − r r ρ (cid:19) , (19)10here ρ is Pearson’s correlation coefficient between the LF and HF estimators. Note that this multifidelityestimator always guarantees variance reduction, since ρ ∈ (0 , r (cid:63) r (cid:63) = − (cid:115) w ρ − ρ , (20)where w = C HF / C LF is the cost ratio between the two fidelities. With this, we can compute the optimalnumber of high-fidelity samples by extrapolation targeting an estimator variance improvement factor (cid:15) : N HF = V [ Q HF M ] (cid:15) (cid:18) − rr + 1 ρ (cid:19) . (21)This follows the extrapolation for MLMC estimators in subsection 3.2. For an extension of the CV MFMCestimators to multiple fidelity models we refer the reader to [44, 71]. We combine ideas from subsections 3.2 and 3.3 to further reduce the variance of our estimators ina MLMF approach (Figure 3c). We apply one control variate for each (high-fidelity) resolution level (cid:96) ,or equivalently apply the multifidelity control variate approach to the Monte Carlo estimators associatedwith the difference Y (cid:96) . MLMF schemes can be implemented for either the same or different numbers ofhigh-fidelity and low-fidelity model levels. For further details, the reader is referred to [45, 47, 48].The MLMF estimator which assumes the same number of high-fidelity and low-fidelity model levels isgiven by E [ Q HF M ] ≈ ˆ Q MLMF M = L (cid:88) (cid:96) =0 (cid:16) ˆ Y HF M (cid:96) ,N HF (cid:96) + α (cid:96) (cid:16) ˆ Y LF M (cid:96) ,N HF (cid:96) − E [ Y LF M (cid:96) ] (cid:17)(cid:17) . (22)As before, we determine the optimal sampling distribution (by level) N HF (cid:96) to compute ˆ Y HF M (cid:96) ,N HF (cid:96) andˆ Y LF M (cid:96) ,N HF (cid:96) by extrapolation. We then calculate the redistribution of computational burden to the moreinexpensive LF model as the set of (independent) samples ∆ LF (cid:96) = r (cid:96) N HF(cid:96) needed to compute ˆ E [ Y LF M (cid:96) ,N LF(cid:96) ]as in subsection 3.2. The solution for the optimal number of samples per level (cid:96) is given by N HF (cid:96) = 2 (cid:15) L HF (cid:88) k =0 (cid:18) V [ Y HF k ] C HF k − ρ (cid:96) (cid:19) Λ k ( r k ) (cid:115) (1 − ρ (cid:96) ) V [ Y HF (cid:96) ] C HF (cid:96) , (23)with the optimal redistribution of samples given by ∆ LF (cid:96) = r (cid:63)(cid:96) N HF (cid:96) where r (cid:63)(cid:96) = − (cid:115) w (cid:96) ρ (cid:96) − ρ (cid:96) (24)and the quantity Λ (cid:96) ( r (cid:96) ) = 1 − r (cid:96) r (cid:96) ρ (cid:96) (25)measures the variance reduction on each level (cid:96) by accounting for the CV effect.At each level, the allocation of samples between LF and HF is controlled by r (cid:96) , which is a functionof the correlation ρ (cid:96) and the computational cost ratio w (cid:96) = C HF (cid:96) / C LF (cid:96) . The optimal variance reductiontherefore corresponds to 1 − Λ( r (cid:63)(cid:96) ). Intuitively this is correct since for an increase in correlation and/or inthe cost ratio, more low-fidelity simulations (i.e., larger r (cid:96) ) are required. To ensure our model correlationsare sufficiently high for all fidelity and discretization levels, we employ the version of the MLMF algorithmdetailed in [47]. 11 . Methods In this study, we developed an interface coupling the three-, one-, and zero-dimensional hemodynamicsolvers with the functionality of Sandia National Laboratories’ Dakota toolkit. Dakota is an extensive open-source framework providing a library of computational tools for optimization, uncertainty quantification,parameter estimation and sensitivity analysis [50]. It allowed us to assign distributions to uncertain pa-rameters, to automatically generate parameter samples, to manage the simultaneous execution of multiplesimulations, and to characterize statistically the resulting QoIs. The MLMF estimators within Dakota havebeen successfully applied to aerospace applications [47], wind power plants [72], and other applications. Ini-tial work by these authors for the application of cardiovascular modeling [73, 74] is extended and improvedin this paper.The Dakota cardiovascular UQ interface was created to be easily adaptable to different UQ scenarios.Most files that connected Dakota to the hemodynamic solvers remained unchanged as the cardiovascularmodels were exchanged. Similarly, changing uncertain parameter distributions and QoIs was achieved withan update to a small number of interface files. The power of Dakota’s behind-the-scenes simulation manage-ment came at runtime. We developed the interface such that after providing the model-specific files and allcardiovascular solver executables, only a single executable was called to launch an uncertainty quantificationstudy on a high performance computing system. Dakota then managed the allocation of each model fidelityand resolution level to achieve a specified convergence tolerance. Upon completion, all statistics for theQoIs were automatically generated. As such, the developed framework was user-friendly and did not requireconstant management of parameter realizations or solution files. We also note that the framework was notlimited to MLMF estimators, but could interface with any of Dakota’s capabilities simply by providing adifferent input file.Specific details of the models used in this paper are discussed in subsections 4.2, 4.3, and 4.5, while theQoIs are discussed in subsection 4.6. A repository detailing the necessary user-provided files, implementation,and usage of the Dakota cardiovascular UQ interface is publicly available on GitHub [75]. The repositoryincludes more detailed information on the simulation process. A schematic representation of the UQ workflowis provided as a general overview (Figure 4).
Four models, healthy and diseased geometry of aorto-femoral and coronary anatomy, each with threefidelities (3D, 1D, and 0D) were used in this study, with original image data and models provided bythe Vascular Model Repository ( ). The realistic diseased models wereconstructed by adapting the geometry of the healthy models to introduce an abdominal aortic aneurysm(AAA) and coronary vessel stenosis (Figure 5).The healthy aorto-femoral model is a patient-specific model of a healthy abdominal aorta with iliacand femoral arteries, with nine outlet branches. The three model fidelities were generated as discussedin section 2, illustrated in Figures 5a, 5b, 5c, and 5d. The diseased aorto-femoral model was adapted inSimVascular to include an AAA at a location typically observed in the clinic. The aneurysm was inducedby increasing the aorta radius by 50%, which is the accepted threshold for classification as an AAA [76].The healthy coronary model includes the aorta, the left and right coronary arteries and ten outletcoronary branches. The fluid dynamics in the coronaries are more complex due to large area differencesbetween coronary and aortic inlets and outlets and the offset in time between pressure and flow waveforms.The three model fidelities were generated as discussed in section 2 above, illustrated in Figures 5e, 5f, 5g,and 5h. The diseased coronary model was adapted in SimVascular to include a stenosis in the LAD, acommon location of coronary stenoses [77]. Significant vessel stenosis is often defined as a reduction invessel diameter of at least 50%; we used a reduction in diameter of 75%.Two zero-dimensional models were constructed for each anatomy. The first, more complex zero-dimensionalmodel uses both RCL and RC blocks for each vessel segment between bifurcations. Including the induc-tance and capacitance elements allows these models to incorporate inertia and wall elasticity. The secondmodel uses only resistance elements for each vessel segment, creating a model simplification with the same12
AKOTA (a)Mesh Convergence (e)Stochastic Parameters(c) Deterministic Parameters(c) DAKOTA Driver(d) 0D Levels 0D Solve(f) QoI Vector (Full)1D Levels 1D Solve(f)3D Levels 3D Solve(f) Outputs(h)All inputs(b)
MediumMedium F i n e F i n e C o a r s e C o a r s e F u l l Simple
QoI selection (g)
Figure 4: The developed workflow utilizes Sandia National Laboratories’ Dakota toolkit to automate the uncertainty quantifi-cation process. (a) Dakota (official logo from [50]) automatically manages all steps within the outer box after the user providesin a single input file (b) with all problem-specific inputs (such as stochastic parameters and their distributions, deterministicparameters, the desired number of simulations, and convergence criterion). The user can easily change (c) the parameters fromstochastic to deterministic and vice versa. Dakota automates (d) the selection of models of different fidelities and mesh levels,with the relevant simulation files provided by the user after (e) a mesh convergence study is conducted to select the appropriatemodels. Dakota also (f) automatically invokes the 0D, 1D or 3D solvers, including all pre- and post-processing, and stores thecomplete set of model outputs. Finally, Dakota extracts (g) the QoIs specified by the user and delivers (h) the moments (mean,standard deviation) of the estimators for the QoIs. uncertain parameters. While this is not the only possible simplified zero-dimensional model, the resistor-only model remains highly correlated with the higher fidelity models. It also performs well as an even lesscomputationally expensive zero-dimensional model, allowing for the two reduced models to be consideredas coarse and medium “resolution levels” in the framework. In the future, additional reduced-order modelscould be added to the hierarchy. Many such examples for one- and zero-dimensional models are enumeratedin [78].
After creating three-, one-, and zero-dimensional models of all four types (aorto-femoral healthy and dis-eased, coronary healthy and diseased), a mesh convergence study was performed to select coarse, medium,and fine resolution levels for each model which demonstrated increasing accuracy compared to a referencemesh for our quantities of interest. These QoIs included time-averaged outlet flow, time-averaged outlet pres-sure, time-averaged pressure spatially-averaged throughout various regions of the model, and time-averagedwall shear stress (TAWSS) values also spatially-averaged on various regions of the model wall, as describedin subsection 4.6. For three- and one-dimensional models, we performed a spatial mesh convergence. Asmentioned in subsection 4.2, in lieu of a mesh convergence for the zero-dimensional models lacking in tradi-tional spatial resolution, a simplified LPN was constructed from only resistor elements to serve the purposeof a coarse model representation.For each of the four model geometries, our three-dimensional mesh convergence study compared theperformance of seven mesh resolutions, ranging from about 50,000 to 2 million elements, against a referencemesh of approximately 4 million tetrahedral elements. At this reference mesh resolution, the QoI values areno longer changing by more than 1% as the mesh is further refined.13 a) (b) (c) (d)(e) (f) (g) (h)
Figure 5: Healthy and diseased aorto-femoral and coronary models of three varying fidelities are used in this study. Specifically,we utilize (a) 0D, (b) 1D, (c) 3D healthy, and (d) 3D diseased (AAA) aorto-femoral models. We also utilize (e) 0D, (f) 1D, (g)3D healthy and (h) 3D diseased (75% stenosis) coronary models. The diseased model regions are circled in (d) and (h).
Deformable wall simulations were used, requiring first a rigid simulation for each mesh level followed bya deformable wall simulation. Steady inlet flow and resistance boundary conditions were used. The totalresistance was equivalent to that of the mean value of the resistance for the boundary conditions used inthe UQ study (see subsection 4.5). Non-linear finite element iterations were regarded as converged whenthe residual norm was below a 0.0001 threshold. Each simulation was run for four cardiac cycles, and theconvergence of each was confirmed to ensure a fair comparison of the mesh levels.Convergence plots for the quantities of interest were used to select the mesh levels used for the multi-levelframework for the UQ study. We confirmed that the finest mesh produced QoI values close to those of thereference mesh in the linear scaling range (Figure 6). After selecting the fine mesh, medium and coarsemeshes were chosen from among the remaining meshes in the linear scaling range.Mesh convergence was also performed for the one-dimensional models. One-dimensional models areconstructed automatically in SimVascular after the lumen segmentation phase of three-dimensional modelbuilding through the “1D-Plugin” feature, which is based on the method developed and discussed in sub-section 2.4. When the one-dimensional models were extracted from the three-dimensional geometry, linearsegments with different diameters were automatically generated to approximate changes in vessel diame-ter throughout the model. Different mesh resolutions for these one-dimensional models were generated byvarying the number of finite elements used to discretize each of these model segments along the z -axis (thevessel centerlines). The coarsest level used 1 element per segment, the finest level used 250 elements, and areference mesh used 500 elements. Details of the selected meshes are shown in Table 1. Mesh quality values14 Number of Mesh Elements05101520 A v e r ag e N o r m a li ze d P r e ss u r e E rr o r Unselected MeshesSelected Meshes
Figure 6: Sample convergence plotfor mesh selection. Selected meshesfor aorto-femoral diseased model demar-cated with red xs. (For interpretation ofthe references to color in this figure leg-end, the reader is referred to the web ver-sion of this article.)
Aorto-Femoral Aorto-Femoral Coronary CoronaryHealthy Diseased Healthy Diseased3D Fine Mesh 1 219 672 1,036 483 1 026,675 1 069 3283D Medium Mesh 446 806 425 416 411 958 401 1703D Coarse Mesh 223 927 200 610 106 681 204 7621D Fine Mesh 50 50 25 251D Medium Mesh 10 10 10 101D Coarse Mesh 5 5 5 5
Table 1: Total number of elements in each of the meshes selected for the UQstudy for the 3D and 1D models. for the aspect ratio and Jacobian were checked for all 3D meshes to ensure they were of high quality for ourUQ study.
After selecting the mesh resolution levels for the three-, one-, and zero-dimensional models of all fourtypes (aorto-femoral healthy and diseased, coronary healthy and diseased), an initial computational costof each type was determined by running each model with the mean value realizations of all stochasticparameters (see subsection 4.5), along with the deterministic parameters, including transient inlet flow.The significant cost difference between the model fidelities and resolution levels lends itself to the MLMFframework. These costs were updated to reflect the averaged cost from the 25 pilot runs for each modelafter the pilot phase (details in subsubsection 5.1.2) of our UQ study (Table 2).
Aorto-Femoral Healthy Aorto-Femoral Diseased Coronary Healthy Coronary DiseasedFidelity & Level Cost Effective Cost Cost Effective Cost Cost Effective Cost Cost Effective Cost3D Fine Mesh 870.80 h 1 667.23 h 1 2 164.61 h 1 1 198.48 h 13D Medium Mesh 228.44 h 2 . × − . × − . × − . × −
3D Coarse Mesh 98.02 h 1 . × − . × − . × − . × −
1D Fine Mesh 11.60 m 2 . × − . × − . × − . × −
1D Medium Mesh 2.95 m 5 . × − . × − . × − . × −
1D Coarse Mesh 1.90 m 3 . × − . × − . × − . × −
0D Full Model 0.49 m 3 . × − . × − . × − . × −
0D Simple Model 0.03 m 6 . × − . × − . × − . × − Table 2: Costs needed to generate simulation results using various hemodynamic solvers and resolution levels. Cost refers tothe cost of one simulation while the effective cost is scaled relative to the cost of the corresponding 3D fine mesh simulation.
Both stochastic and deterministic parameters were used for all hemodynamic models in this study. Whenconsidering healthy and diseased models of the same anatomy, identical parameters were used to isolate theeffect of adapting the model geometry. All models included eight independent uncertain inputs with valuesdrawn from a uniform distribution with limits equal to ±
30% of the literature values for each parameter(Table 3). 15 orto-Femoral Ranges Coronary RangesUncertain Parameter Min Max Min MaxBC: Total R . × . × . × . × BC: Total C . × − . × − . × − . × − BC: Ratio of R p /R total . × − . × − . × − . × − BC: Ratio of R p /R total (renal arteries) 1 . × − . × − — —Young’s Modulus 4 . × . × . × . × Young’s Modulus (coronary arteries) — — 8 . × . × Inlet waveform total flow 5 . × . × . × . × Blood Density 7 . × − . . × − . . × − . × − . × − . × − Table 3: Uncertain parameters and ranges for uniform distributions of these parameters for aorto-femoral and coronary models.Resistances have units of dyn s / cm . Capacitances have units of cm / dyn. Young’s moduli have units of Pa. Flow has unitsof cm s − . Density has units of g / cm . Viscosity has units of Pa × s. All parameters were assigned reference (mean) values chosen to match targets from clinical literature.RCR boundary conditions were applied to all outlets of both models to represent the downstream vasculature.The total resistance and capacitance in the model were tuned using a three-element Windkessel model toproduce a physiologic pressure waveform, then distributed proportional to the vessel outlet areas [79]. Themean values for the proximal ( R p ) and distal ( R d ) resistance split was consistent with those used in [80]. Asin that study, a different ratio was used for the renal arteries to account for their unique flow features in theaorto-femoral models. Further details on boundary condition assignment can be found in [80]. Mean valueschosen for the Young’s moduli were corroborated by other studies [6, 81, 82]. A physiologic inlet waveformwas applied to the aortic inlet of both models. As the aorto-femoral model inlet is in the abdominal aorta andthe coronary model inlet is in the thoracic aorta, these waveforms differed slightly, but both were consistentwith typical clinical assessments for their respective anatomy. All analysis presented in this paper was alsocarried out with steady inlet waveforms in a preliminary study [73].Though all stochastic parameters used in this study were uniformly distributed and independent, anadvantage of the MLMF method is the flexibility it affords. Generally, sampling based methods such as thisallow for heterogeneous sources of uncertainty. Parameters can be a combination of independent parameterswith distributions assumed from literature data and parameters assimilated or sampled from clinical data,for example using a Markov Chain Monte Carlo approach as in [34]. Allowing for a spectrum of methodsfor defining uncertainties is an asset of sampling-based approaches such as the MLMF method. Both local and global quantities of interest were included in this study. These were divided into fourmain categories: time-averaged flow and pressure values at model outlets, and time-averaged blood pressureand wall shear stress values spatially-averaged over individual branches, strips of interest, and the entiremodel. The strips of interest are the highlighted regions seen in Figures 12 and 13. A total of 132 and148 QoIs were computed for the aorta and coronary models, respectively. The same quantities of interestmeasured at the same regions of interest were used for the healthy and diseased geometries to facilitatecomparisons.While all QoIs used in this study are single values computed by temporally averaging the results overone cardiac cycle, the MLMF method is capable in principle of handling QoI estimators at multiple timepoints, for example to generate flow and pressure waveforms with confidence intervals over a cardiac cycle.Estimators for a single time step were used here to demonstrate the performance of the MLMF method.As outlined in section 3, high correlations between the model fidelities at each resolution level arenecessary for the success of our MLMF methods. We find high ( > .
6) averaged values for Pearson’scorrelation coefficient across all global and local QoI categories ( Table 4). As expected, the global QoIs are16ore highly correlated than the local QoIs. Within the local QoIs, higher 3D-3D mesh level correlationsare observed for the coarser resolution levels than at the fine resolution level, due to the fact that the mostspatially resolved models are less well approximated by the low-fidelity models. We look at the correlationbetween both low- and high-fidelity models as well as within the mesh resolution levels of the high-fidelitymodel.
5. Results
To justify the additional overhead of introducing lower-fidelity models for the MLMF estimators, we com-pare the performance of the MLMF estimators to Monte Carlo (MC) and multilevel Monte Carlo (MLMC)estimators, which rely only on three-dimensional models. We see that across a variety of metrics, the MLMFestimators outperform the MC and MLMC estimators, demonstrating the power of the method when subjectto constrained computational budgets. MLMF schemes using only one-dimensional low fidelity models andincluding both one- and zero-dimensional models were evaluated for all four healthy and diseased modelgeometries, as discussed in subsection 3.4. We will refer to the setup including zero-dimensional modelsas the scheme and the setup excluding zero-dimensional models as the scheme. For the3D-1D scheme, there are three spatial mesh resolutions of both the HF and LF (1D) models, leading tothree discrepancy levels Y (cid:96) . For the 3D-1D-0D scheme, the same three HF mesh resolutions are employed.Now, the two lower resolutions of the 1D LF model are replaced with 0D models. The lowest level is the 0Dmodel comprised of only resistor elements, while the middle level is a 0D model incorporating resistance,capacitance, and inductance elements (see subsection 2.3 for additional details). The highest LF resolutionlevel remains unchanged from the 3D-1D scheme, and is still the finest 1D resolution. As we observed in Table 4, the correlations between the low- and high-fidelity models are high acrossour various QoIs. For the global QoIs, we see very high ( > .
9) correlations between both the 0D and3D and the 1D and 3D models. These observed correlations, combined with the inexpensive nature of thelow-fidelity models, led us to investigate the bias in the low-fidelity models, in order to justify the use of ourmultifidelity estimators instead of simply relying on MC estimators from the low-fidelity models.We utilized a Monte Carlo method with only the 0D RCL model simulation results or only the 1D finemodel results for various quantities of interest, and compared these to the 3D fine model results. We foundvarying degrees of bias for different QoIs between the LF model MC estimators and the HF MC estimator.Box plots for representative global and local QoIs from the aorto-femoral model show that even for theglobal QoI (flow), the 0D RCL MC estimator has a significantly different mean value than the 3D or 1DMC estimators do (Figure 7). The 3D and 1D model MC estimators do agree in this case. However, for thelocal QoI (TAWSS), all three MC estimators exhibit different mean values. In this case, the 0D model failsto generate any simulations with QoI values near the 3D MC estimator mean value.The bias exhibited by the 0D and 1D MC estimators is not a simple correction value, as the bias valuediffers across QoIs. This prevents solely the 0D or 1D model from being used without any 3D simulations.For this, and other, reasons, the general applicability of the MLMF method appeals to us. Even with a smallnumber of high-fidelity evaluations we can rely heavily on a biased 0D model without having to explicitlycalculate or consider this bias. Due to the way the MLMF estimators are constructed, any bias in the LFmodel is automatically compensated for.
In our uncertainty quantification workflow, we first complete a pilot run , an a priori prescribed numberof simulations of each model level and fidelity, chosen here to be 25 samples of each fidelity discrepancy level Y HF (cid:96) and Y LF (cid:96) for (cid:96) = 0 , ,
2, for a total of 250 simulations. Recall that for all (cid:96) > Y (cid:96) is the discrepancybetween a lower-resolution and higher-resolution simulation, amounting to a total of 50 simulations foreach Y (cid:96) , (cid:96) >
0. This educated sample size was informed as a reasonable portion of the overall available17 utlet Flow Outlet Pressure Model Pressure Model TAWSSModel and Mesh 3D-1D 3D-0D 3D-1D 3D-0D 3D-1D 3D-0D 3D-1D 3D-0D A F H Coarse 0.997 0.995 0.995 0.995 0.995 0.992 0.825 0.647Medium 0.997 0.996 0.995 0.989 0.995 0.990 0.843 0.782Fine 0.997 — 0.994 — 0.995 — 0.792 — A F D Coarse 0.998 0.995 0.996 0.996 0.996 0.992 0.834 0.612Medium 0.993 0.992 0.995 0.989 0.994 0.990 0.880 0.812Fine 0.998 — 0.994 — 0.995 — 0.818 — C H Coarse 0.999 0.999 0.998 0.999 0.998 0.999 0.854 0.775Medium 0.999 0.999 0.998 0.995 0.998 0.995 0.851 0.850Fine 0.999 — 0.998 — 0.998 — 0.811 — C D Coarse 0.999 0.999 0.997 0.999 0.997 0.999 0.847 0.744Medium 0.999 0.999 0.997 0.993 0.997 0.994 0.853 0.798Fine 0.999 — 0.996 — 0.996 — 0.799 —(a) HF-LF Model Correlations
Outlet Flow Outlet Pressure Model Pressure Model TAWSSModel and Mesh 3D-3D 3D-3D 3D-3D 3D-3D A F H Coarse–Medium 0.999 0.999 0.999 0.999Medium–Fine 0.999 0.999 0.999 0.998 A F D Coarse–Medium 0.996 0.999 0.999 0.939Medium–Fine 0.995 0.998 0.998 0.924 C H Coarse–Medium 0.999 0.999 0.999 0.999Medium–Fine 0.999 0.999 0.997 0.998 C D Coarse–Medium 0.999 0.999 0.999 0.999Medium–Fine 0.999 0.999 0.998 0.998(b) HF Mesh Resolution Level Correlations
Table 4: Averaged Pearson correlation coefficients for each QoI category for all models between (a) the 3D and 1D models(3D-1D) and the 3D and 0D models (3D-0D) of the same resolution level and between (b) the 3D models of different meshresolutions levels. The former correlations are needed for the control variate portion and the latter are needed for the MLMCportion of our approach. The correlations are very high for the global (flow and pressure) QoIs, and slightly lower for the local(TAWSS) QoIs. However, even for the local QoIs, the correlations are large enough to produce good performance from theMLMF UQ methods. No 0D model fine resolution level was used in our study, so no correlation is reported. Abbreviations:aorto-femoral healthy (AFH), aorto-femoral diseased (AFD), coronary healthy (CH), and coronary diseased (CD). a) (b) Figure 7: Examining the bias between MC estimators for different model fidelities for (a) a global flow QoI and (b) a localTAWSS QoI in the aorto-femoral model. MC estimators were built individually from simulations of the 3D Fine model, the1D Fine model, and 0D RCL model. computational budget. The purpose of the pilot run is threefold: to estimate (i) the variances of the QoIs oneach level, (ii) the correlations between high and low fidelity models, and, possibly, (iii) the average modelevaluation cost per level. From the estimator variances, we compute the normalized confidence intervalof each QoI after the pilot run. After assessing the initial normalized confidence intervals, we estimatehow many additional simulations are needed to converge the solution to the desired normalized confidenceinterval value for the QoIs.
The normalized confidence interval (nCI) of a QoI estimator Q is defined here in accordance with pre-vious literature [46, 47] as six times the coefficient of variation, or the size of three standard deviationscorresponding to a 99.7% confidence interval normalized by the expected value, nCI [ Q ] = 6 (cid:112) V [ Q ] E [ Q ] = 6 σµ , (26)for each QoI. As such, smaller values of nCI are preferred as they indicate a higher level of confidence, i.e.tighter confidence intervals, in the estimator Q . For this study, a value of nCI [ Q ] = 0 .
01 was arbitrarilychosen as the target, though this may be more accurate than needed for different applications. A moreextensive study of the proper nCI target value for clinical applications will be the subject of future work.As explained in subsection 3.4, extrapolation can be performed to estimate the additional number ofsimulations of each model and fidelity needed to improve the estimator variance by a factor (cid:15) , with (cid:15) ∈ (0 , (cid:15) can be changed to target different variances and, by extension, different nCIs for the QoIs, asnCI is a function of variance. The formulas for the optimal number of samples needed for extrapolation aregiven by (16), (21), and (23) for the MLMC, MFMC, and MLMF estimators, respectively.MLMF estimators outperformed MC and MLMC methods when considering extrapolation with differenttarget variances (i.e. different (cid:15) factors) for a specific QoI. In this context, outperforming means thatthe computational cost to obtain a given variance is lower for the MLMF schemes. The superior MLMFperformance holds for both healthy and diseased models and both global and local QoIs (Figure 8).When extrapolating to a nCI of 0.01, the exact breakdown of projected costs for the same QoIs as in Fig-ure 8 shows that the MLMF methods offer one to two orders of magnitude improvement in computationalcost over the MLMC method and one to three orders of magnitude improvement over the MC method(Tables 5a and 5b for the aorto-femoral and coronary models, respectively). Lower improvement factors are19orto-Femoral Healthy (Estimator Standard Deviation) Equivalent HF Fine Runs10 − − S t a nd a r d D e v i a t i o n (a) Equivalent HF Fine Runs10 − S t a nd a r d D e v i a t i o n (b) Equivalent HF Fine Runs10 − S t a nd a r d D e v i a t i o n (c) Equivalent HF Fine Runs10 − S t a nd a r d D e v i a t i o n (d) Coronary Healthy (Estimator Standard Deviation) Equivalent HF Fine Runs10 − − S t a nd a r d D e v i a t i o n (e) Equivalent HF Fine Runs10 − S t a nd a r d D e v i a t i o n (f) Equivalent HF Fine Runs10 − S t a nd a r d D e v i a t i o n (g) Equivalent HF Fine Runs10 − S t a nd a r d D e v i a t i o n (h) Monte Carlo Multilevel MLMF 3D-1D MLMF 3D-1D-0D
Figure 8: Comparing the performance of estimators from MC, MLMC, and MLMF UQ schemes for (top) aorto-femoral and(bottom) coronary healthy models. The values shown are extrapolated expected costs (in comparable units of the equivalentnumber of 3D fine simulations) to obtain various estimator standard deviation values for representative QoIs from each category.Standard deviation cost curves are shown for representative outlet flow QoIs in figures (a) and (e), outlet pressure QoIs infigures (b) and (f), model time-averaged pressure QoIs in figures (c) and (g), and model TAWSS QoIs in figures (d) and (h) forthe aorto-femoral and coronary models, respectively. These quantities are spatially and temporally averaged over the outletface for the flow and pressure QoIs, throughout regions of the model volume for the model pressure QoIs, or over regions ofthe model surface for model TAWSS QoIs respectively. seen for local QoIs, as the local QoIs rely heavily on the detailed spatial information from the HF modelevaluations. This means that the extrapolated costs will be more similar across UQ methods for these QoIs.The extrapolated cost reported includes the initial cost of the pilot run. Extrapolated costs are reported incomparable units of the equivalent number 3D fine simulations . This is calculated as N DEq,F ine = (cid:0) N DF ine C DF ine + N DMed C DMed + N DCoa C DCoa + N DF ine C DF ine + N DMed C DMed + N DCoa C DCoa (cid:1) /C DF ine , (27)where C F idelityLevel is the cost of one simulation of that level and fidelity and N F idelityLevel includes the total numberof all pilot run simulations and all additional simulations needed for the extrapolation to nCI [ Q ] = 0 . ealthy Model Diseased ModelQoI and Method Effective Cost Effective Cost O u t l e t F l o w MC 9 758.0 9 471.0MLMC 1 293.6 2 593.8MLMF (3D-1D) 77.9 468.3MLMF (3D-1D-0D) 66.5 440.2 O u t l e t P r e ss u r e MC 19 340.0 18 374.0MLMC 2 915.4 4 588.1MLMF (3D-1D) 82.2 282.2MLMF (3D-1D-0D) 107.0 282.2 M o d e l P r e ss u r e MC 19 598.0 19 050.0MLMC 3 016.3 4 430.4MLMF (3D-1D) 119.4 321.9MLMF (3D-1D-0D) 210.5 369.2 M o d e l T A W SS MC 23 94.0 22 581.0MLMC 5 578.7 4 695.4MLMF (3D-1D) 1 151.4 1 922.4MLMF (3D-1D-0D) 581.2 847.0(a) Aorto-Femoral Model Healthy Model Diseased ModelQoI and Method Effective Cost Effective Cost O u t l e t F l o w MC 10 656.0 10 619.0MLMC 435.1 1 095.5MLMF (3D-1D) 41.0 49.4MLMF (3D-1D-0D) 39.5 43.7 O u t l e t P r e ss u r e MC 20 467.0 20 587.0MLMC 948.1 2 480.6MLMF (3D-1D) 46.6 65.1MLMF (3D-1D-0D) 40.4 44.7 M o d e l P r e ss u r e MC 20 213.0 20 330.0MLMC 934.1 2 445.8MLMF (3D-1D) 45.1 65.0MLMF (3D-1D-0D) 40.4 46.4 M o d e l T A W SS MC 21 748.0 16 391.0MLMC 2 066.7 2 391.4MLMF (3D-1D) 901.9 1 303.8MLMF (3D-1D-0D) 290.6 269.7(b) Coronary Model
Table 5: Comparison of the performance of different UQ methods. Each of the four UQ methods compared (Monte Carlo,multilevel, MLMF with 3D, 1D, and 0D models, and MLMF with only 3D and 1D models) requires differing costs to achievea nCI value of 0.01. For four representative outlet flow, outlet pressure, time-averaged model pressure, and TAWSS QoIs,the cost of each method (in units of the equivalent number of 3D fine simulations) is shown. The same representative QoI isused for both healthy and diseased models. These QoIs are spatially and temporally averaged over the outlet face, throughoutregions of the model volume, or over regions of the model surface for the flow and pressure, model pressure, and model TAWSSQoIs respectively.
To demonstrate how the cost burden is shifted to the least expensive models with our estimators, alevel-by-level breakdown of the extrapolated simulations is shown for one global and one local QoI. Fromthe aorto-femoral model we have chosen the flow and TAWSS QoIs (Table 6a and Table 7a). From thecoronary model we have chosen the pressure and TAWSS QoIs (Table 6b and Table 7b). We look at boththe breakdown in terms of pure number of simulations of each model level (Table 6), where a single LF andHF model are considered equivalent, and the breakdown in terms of contribution to the total cost of theextrapolation (Table 7).As the LF models are much cheaper to evaluate than the HF models (Table 2), this allocation keeps thecost of the MLMF methods below the MC or MLMC methods, even though more than twice the additionalsimulations are needed for the MLMF methods than for the MC or MLMC methods. Both the MLMC andMLMF methods allocate the majority of their simulations to the least expensive model. For the MLMCmethod, the majority of additional simulations are HF coarse. For the MLMF method, the majority are LFcoarse. However, in the MLMF cases the LF coarse simulations do not comprise to the majority of the totalcost of the extrapolation.
When relying on estimator values and their variances from only the pilot run, MLMF estimators out-perform the MC and MLMC methods. In this context, outperforming means that the estimator variance islower from the MLMF schemes than from the other methods. In fact, with the size of this pilot, we wouldconsider the MLMF estimators to be converged after the pilot run. The size of the pilot run can be shapedsuch that the variance is larger and a specific variance or nCI can be targeted with extrapolation followingthe pilot phase. A computational budget constraint on, for instance, the number of most expensive modelevaluations can also be enforced.We show the convergence of four quantities of interest to the final estimated mean and variance obtainedfrom the pilot run (Figure 9). As we build up the full pilot sample, we can see the differences in variance21 utlet Flow Model TAWSSMC MLMC MLMF MLMF MC MLMC MLMF MLMF3D-1D 3D-1D-0D 3D-1D 3D-1D-0D H e a l t h y LF Coarse — — 96.6% 98.9% — — 77.9% 94.7%LF Medium — — 1.7% 0.9% — — 10.8% 4.1%LF Fine — — 1.4% 0.2% — — 10.1% 1.1%HF Coarse — 98.1% 0.2% 0.0% — 84.4% 1.1% 0.0%HF Medium — 1.1% 0.0% 0.0% — 8.2% 0.0% 0.0%HF Fine 100% 0.8% 0.0% 0.0% 100% 7.4% 0.1% 0.0%
Total No. Sims. 9 758 10 386 133 843 883 690 23 940 25 469 543 309 3 060 780 D i s e a s e d LF Coarse — — 92.2% 97.1% — — 78.6% 94.0%LF Medium — — 5.3% 2.4% — — 14.1% 5.3%LF Fine — — 2.3% 0.4% — — 5.6% 0.7%HF Coarse — 89.9% 0.1% 0.0% — 84.3% 1.4% 0.0%HF Medium — 6.7% 0.1% 0.0% — 10.3% 0.1% 0.0%HF Fine 100% 3.4% 0.1% 0.0% 100% 5.4% 0.1% 0.0%
Total No. Sims. 9 471 18 626 363 442 2 240 202 22 581 27 537 673 023 3 407 848 (a) Aorto-Femoral Model
Outlet Pressure Model TAWSSMC MLMC MLMF MLMF MC MLMC MLMF MLMF3D-1D 3D-1D-0D 3D-1D 3D-1D-0D H e a l t h y LF Coarse — — 96.4% 98.4% — — 82.4% 98.6%LF Medium — — 2.0% 0.9% — — 11.6% 0.9%LF Fine — — 1.4% 0.7% — — 4.6% 0.4%HF Coarse — 99.6% 0.1% 0.0% — 91.6% 1.3% 0.0%HF Medium — 0.3% 0.0% 0.0% — 6.3% 0.1% 0.0%HF Fine 100% 0.1% 0.0% 0.0% 100% 2.1% 0.0% 0.0%
Total No. Sims. 20 467 24 799 197 712 653 682 21 748 27 298 932 759 3 717 084 D i s e a s e d LF Coarse — — 96.7% 98.5% — — 91.2% 99.0%LF Medium — — 1.3% 0.6% — — 3.5% 0.4%LF Fine — — 1.9% 0.9% — — 4.0% 0.6%HF Coarse — 99.6% 0.1% 0.0% — 94.0% 1.2% 1.2%HF Medium — 0.3% 0.0% 0.0% — 2.8% 0.0% 0.0%HF Fine 100% 0.1% 0.0% 0.0% 100% 3.2% 0.0% 0.0%
Total No. Sims. 20 587 24 188 214 491 563 280 16 391 16 615 699 517 2 894 913 (b) Coronary Model
Table 6: Comparison of the extrapolated number of simulations to obtain nCI [ Q ] = 0 .
01 for different UQ methods. Thedistribution of number of simulations by model level and fidelity is shown as a percentage of the total number of simulationsneeded for extrapolation. The total number of simulations is shown on the final row for the healthy and diseased models. Thesame representative QoI is used for both healthy and diseased models. Values marked as 0.0% contribute negligibly to theoverall percentage of simulations. between the methods arise. The MC, MLMC, and both MLMF estimators are shown at increasing sub-samples of the full set of pilot run simulations. These estimators are plotted against the equivalent cost ofthe samples making up the estimators. The final, rightmost points in each plot are the estimators obtainedfrom the pilot study.As the sample comprising the estimator grows, each method becomes more accurate. We can see thatthe MLMF methods obtain very accurate, i.e. small variance and nCI, estimators at a much lower cost thanthe MLMC or MC estimators.Only the healthy model results are presented here, as the diseased models exhibit the same trends.Similarly, specific quantities of interest from each category were used. The conclusions drawn from thesequantities are consistent with other quantities of interest of the same type (outlet flows and pressures,time-averaged model pressure, or TAWSS values). 22 utlet Flow Model TAWSSMC MLMC MLMF MLMF MC MLMC MLMF MLMF3D-1D 3D-1D-0D 3D-1D 3D-1D-0D H e a l t h y LF Coarse — — 6.0% 0.8% — — 1.3% 0.3%LF Medium — — 0.3% 0.1% — — 0.5% 0.2%LF Fine — — 0.7% 0.6% — — 1.3% 1.4%HF Coarse — 88.7% 40.5% 36.9% — 43.4% 59.8% 5.9%HF Medium — 3.2% 12.0% 14.1% — 14.0% 7.0% 8.4%HF Fine 100% 8.1% 40.5% 47.4% 100% 42.6% 30.0% 83.8%
Total Cost (Hours) 8 497 300 1 126 500 67 833 57 927 20 847 000 4 857 900 1 002 700 506 120 D i s e a s e d LF Coarse — — 2.7% 0.4% — — 1.0% 0.3%LF Medium — — 0.4% 0.2% — — 0.5% 0.3%LF Fine — — 0.6% 0.6% — — 0.7% 0.9%HF Coarse — 54.4% 7.9% 10.3% — 41.7% 42.7% 7.7%HF Medium — 15.5% 27.9% 27.6% — 19.3% 10.4% 16.3%HF Fine 100% 30.1% 60.4% 60.9% 100% 39.0% 44.7% 74.5%
Total Cost (Hours) 6 319 300 1 730 600 312 490 293 710 15 067 00 3 132 900 1 282 700 565 160 (a) Aorto-Femoral Model
Outlet Pressure Model TAWSSMC MLMC MLMF MLMF MC MLMC MLMF MLMF3D-1D 3D-1D-0D 3D-1D 3D-1D-0D H e a l t h y LF Coarse — — 3.4% 0.4% — — 0.7% 0.3%LF Medium — — 0.2% 1.1% — — 0.3% 0.%LF Fine — — 0.3% 1.2% — — 0.2% 0.6%HF Coarse — 94.7% 15.9% 4.9% — 44.0% 47.1% 4.1%HF Medium — 2.1% 14.3% 16.4% — 22.1% 15.4% 27.6%HF Fine 100% 3.2% 65.9% 76.0% 100% 33.9% 36.3% 66.4%
Total Cost (Hours) 44 303 000 3 052 200 100 830 87 551 47 076 000 4 473 700 1 952 300 628 990 D i s e a s e d LF Coarse — — 5.0% 0.6% — — 0.8% 0.5%LF Medium — — 0.2% 0.3% — — 0.1% 0.2%LF Fine — — 0.6% 1.2% — — 0.2% 0.7%HF Coarse — 97.8% 33.6% 9.5% — 65.7% 64.2% 16.8%HF Medium — 0.9% 13.1% 19.0% — 6.5% 5.0% 14.2%HF Fine 100% 1.3% 47.6% 69.4% 100% 27.7% 29.8% 67.6%
Total Cost (Hours) 24 673 000 2 972 900 78 018 53 530 19 644 000 2 866 000 1 562 600 323 170 (b) Coronary Model
Table 7: Comparison of the cost of the extrapolated simulations to obtain nCI [ Q ] = 0 .
01 for different UQ methods. The costdistribution by model level and fidelity is shown as a percentage of the total cost of all simulations needed for extrapolation.The total cost of simulations, in hours, is shown on the final row for the healthy and diseased models. The same representativeQoI is used for both healthy and diseased models.
We introduced the process of extrapolation to estimate the additional number of simulations of eachmodel and fidelity needed to decrease the estimator variance by a factor (cid:15) . The formulas for the optimalnumber of samples needed for extrapolation are given by (16), (21), and (23) for the MLMC, MFMC, andMLMF estimators, respectively. These formulas have been used successfully for applications such as enginenozzles, scramjets, and wind turbines [47, 72]. To validate these for cardiovascular modeling, we neededto ensure that the cost of converging to various target estimator variances is reflective of the projectedextrapolated costs.We utilized the following validation procedure. We first performed a pilot run, with 15 of the lowestdiscrepancy level Y (cid:96) = Y , 10 of the middle discrepancy level Y , and 5 of the highest discrepancy level Y .This is a smaller pilot than used elsewhere, producing larger estimator variances. With the QoI estimatorsand their variances computed from this data, we chose target estimator variances equal to an improvementfactor of (cid:15) = [0 . , . , . V target [ ˆ Q MLMF ] = (cid:15) V pilot [ ˆ Q ML ] , (28)where ˆ Q ML is the MLMC estimator defined in (14) and ˆ Q MLMF is the MLMF estimator defined in (22). Thisequation is consistent with the convergence criterion within Dakota. After selecting the target variances, we23orto-Femoral Healthy F l o w ( c m / s ) (a) P r e ss u r e ( mm H g ) (b) P r e ss u r e ( mm H g ) (c) T A W SS ( g / c m · s ) (d) Coronary Healthy
10 20 30 40Equivalent HF Fine Runs0.250.300.350.400.45 F l o w ( c m / s ) (e)
10 20 30 40Equivalent HF Fine Runs80.00100.00120.00140.00 P r e ss u r e ( mm H g ) (f)
10 20 30 40Equivalent HF Fine Runs80.00100.00120.00140.00 P r e ss u r e ( mm H g ) (g)
10 20 30 40Equivalent HF Fine Runs8.0010.0012.0014.0016.00 T A W SS ( g / c m · s ) (h) Monte Carlo Multilevel MLMF 3D-1D MLMF 3D-1D-0D
Figure 9: Comparing the performance of UQ estimators from four methods for global and local QoIs. The four methods areMonte Carlo using 3D models, multilevel using 3D models, MLMF using 3D and 1D models, and MLMF using 3D, 1D, and0D models. Estimators with their confidence intervals are shown converging to the values in the pilot run for representativeoutlet flow QoIs in figures (a) and (e), outlet pressure QoIs in figures (b) and (f), model time-averaged pressure QoIs in figures(c) and (g), and model TAWSS QoIs in figures (d) and (h) for the aorto-femoral and coronary models, respectively. Thesequantities are spatially and temporally averaged over the outlet face for the flow and pressure QoIs, throughout regions of themodel volume for the model pressure QoIs, or over regions of the model surface for model TAWSS QoIs respectively. then extrapolated to determine the optimal number of additional samples of each model level and fidelity.After then running the simulations to convergence in Dakota, the final costs were approximately thosepredicted by our extrapolation for the MLMF method (Figure 10).
We now examine the differences between estimators from the MLMF 3D-1D and 3D-1D-0D schemes forthe healthy aorto-femoral and coronary models (Figure 11). These results are presented for the healthymodels only for the sake of clarity, and similar trends hold for the diseased cases.
We assessed the initial nCIs of our estimators from the pilot run samples and estimated how manyadditional simulations are needed to converge the solution to a desired estimator variance or nCI valuefor the QoIs. The Dakota workflow continues to automatically manage the UQ study until convergence isachieved.Calculating the initial nCI from the pilot run, we found that while the global QoIs have similar nCIs forthe two schemes, much better nCIs are achieved for local quantities of interest with the 3D-1D-0D scheme24 Equivalent HF Fine Runs10 × S t a nd a r d D e v i a t i o n (a) Equivalent HF Fine Runs10 × × × S t a nd a r d D e v i a t i o n (b) Monte Carlo Extrapolation Multilevel Extrapolation MLMF 3D-1D-0D Extrapolation MLMF 3D-1D-0D Actual
Figure 10: Validating the extrapolation of estimators for the MLMF UQ schemes for (a) aorto-femoral and (b) coronary healthymodels. The values shown in green pentagons and triangles, respectively, are the projected and actual extrapolated costs (inunits of the equivalent number of 3D fine simulations) to obtain various estimator standard deviations for representativeQoIs. The extrapolated costs for the MC and MLMC estimators are also shown for comparison. Representative time- andspatially-averaged (a) outlet and (b) model pressure QoIs are shown for the aorto-femoral and coronary model, respectively.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) than with the 3D-1D scheme for both aorto-femoral and coronary models (Figures 11a and 11d, respectively).We demonstrate this by computing the average nCI across several global and local QoI categories for bothMLMF schemes, including the variability within each group. Local quantities of interest are more difficultto resolve than global quantities, and we find the largest variability occurs within the TAWSS categories.Between the aorto-femoral and coronary model geometries, we see similar nCI values in addition to similartrends between categories.
Extrapolating from the pilot run nCIs to obtain nCI [ Q ] = 0 .
01 for all QoIs, we see that the extrapolatedcost is lower for the 3D-1D-0D scheme than for the 3D-1D scheme across most QoI categories for both theaorto-femoral and coronary models (Figures 11b and 11e, respectively). For the other categories, the costsare equivalent within the margin of variability in each category. As before, we see the most significantdifference between the two MLMF methods for the local QoIs. This trend towards lower cost of the 3D-1D-0D scheme follows directly from the fact that this scheme has better nCIs after the pilot run (Figures 11aand 11d). The extrapolation calculation for the estimated number of additional samples required relies onthe nCIs obtained from the pilot run. Therefore, since the 3D-1D-0D scheme had smaller nCI values tobegin with, it follows that the 3D-1D-0D scheme also requires a smaller additional cost to decrease thatnCI to below 0.01. The extrapolation is also reliant on correlations and the cost ratio between high andlow fidelity models. Both the one- and zero-dimensional models have similarly high correlations with thethree-dimensional model as well as favorable cost ratios, though the 3D-1D-0D scheme has better cost ratiosfor the coarse and medium low fidelity models.Looking at the global and local QoI categories, we see that it is more expensive, on the order of 10 to 100times higher, to obtain the target nCI for the local QoIs than the global QoIs. This is also due to the fact25orto-Femoral Healthy O u t l e t F l o w s O u t l e t P r e ss u r e s B r a n c h P r e ss u r e s S tr i p P r e ss u r e s B r a n c h T A W SS s S tr i p T A W SS s QoIs0 . . . . n C I ( p il o tr un ) (a) O u t l e t F l o w s O u t l e t P r e ss u r e s B r a n c h P r e ss u r e s S tr i p P r e ss u r e s B r a n c h T A W SS s S tr i p T A W SS s QoIs01000200030004000 E x tr a p o l a t e d c o s t( fi x e dn C I ) (b) O u t l e t F l o w s O u t l e t P r e ss u r e s B r a n c h P r e ss u r e s S tr i p P r e ss u r e s B r a n c h T A W SS s S tr i p T A W SS s QoIs0 . . . . . n C I ( fi x e d H F c o s t) (c) Coronary Healthy O u t l e t F l o w s O u t l e t P r e ss u r e s B r a n c h P r e ss u r e s S tr i p P r e ss u r e s B r a n c h T A W SS s S tr i p T A W SS s QoIs0 . . . . . . . n C I ( p il o tr un ) (d) O u t l e t F l o w s O u t l e t P r e ss u r e s B r a n c h P r e ss u r e s S tr i p P r e ss u r e s B r a n c h T A W SS s S tr i p T A W SS s QoIs0500100015002000 E x tr a p o l a t e d c o s t( fi x e dn C I ) (e) O u t l e t F l o w s O u t l e t P r e ss u r e s B r a n c h P r e ss u r e s S tr i p P r e ss u r e s B r a n c h T A W SS s S tr i p T A W SS s QoIs0 . . . . . n C I ( fi x e d H F c o s t) (f) Figure 11: Comparing the performance of MLMF estimators from 3D and 1D models to the MLMF estimators from 3D, 1D,and 0D models for the (top) aorto-femoral and (bottom) coronary healthy models. Averaged nCIs for various groups of QoIscalculated after a fixed pilot run of 250 simulations are shown in plots (a) and (d). Averaged projected extrapolated cost(in units of the equivalent number of 3D fine simulations) to obtain a nCI ≤ .
01 are shown in plots (b) and (e). Averagedprojected nCI values with the computational budget for all HF simulations equivalent to the cost of 50 3D fine simulations areshown for various QoI categories in plots (c) and (f). Representative results for the healthy anatomies are shown here. that the local quantities generally had the largest (worst) nCIs after the pilot run. We also see the largestintra-category variation in the cost for local TAWSS quantities. As local quantities likely rely more heavilyon high fidelity simulations to resolve their estimators, this cost variation is expected, as the high fidelitysimulations have larger cost range amongst the levels than the low fidelity simulations.
The results presented in this section thus far may seem to suggest that the 3D-1D-0D scheme is almostalways preferable to the 3D-1D scheme. However, it is likely that many UQ studies will be constrained bysome computational budget. We needed to investigate if the 3D-1D-0D scheme continued to outperform the3D-1D scheme in this realistic setting. Extrapolation to a specific nCI may prove computationally prohibitivefor some model applications. Instead, with a fixed budget, it is possible to compute the projected nCI forestimators subject to a computational budget constraint. Here, we report the projected nCI values for a high26delity computational budget equivalent to the cost of 50 3D fine mesh simulations, which we considered tobe a realistic UQ scenario. This means that the combined cost of all simulations of all high fidelity modelinglevels must have an equivalent to 50. Note that only the HF models are constrained here, as the simulationsat all LF levels are obtained after optimization by the algorithm and there is no control over their numbernor their total cost. With that said, many thousands of LF simulations can be run for the cost of one HFsimulation. With these constrained budgets, for the aorto-femoral model we see the 3D-1D-0D stronglyoutperforming the 3D-1D scheme for local QoIs, and performing similarly for global QoIs (Figure 11c). Forthe coronary model, the 3D-1D-0D scheme again outperforms the 3D-1D scheme, with smaller nCI valuesachieved across all QoI categories (Figure 11f). Recall for nCIs, smaller values are preferable as this amountsto higher estimator certainty. When including the overall computational cost, adding in the LF cost, weexpect that indeed the 3D-1D-0D should perform better across all QoIs and models. All things equal (i.e.for the same correlation) the optimal r factor controlling the number of LF simulations for each level isequal between 1D and 0D, but the total cost will be far less for 0D simulations. We now examine the difference in estimator performance for healthy and diseased models. We considerlocal QoIs, such as the spatially-averaged TAWSS value in specific regions on the model, and global QoIs(time-averaged flows and pressures at outlets, spatial- and time-averaged pressure in various regions of themodel) separately. Here we use the estimator values and variances from the 3D-1D-0D MLMF scheme, asthese are more accurate after the pilot run (see subsection 5.3). The conclusions drawn from these specificquantities are reflective of the behavior of other QoIs from the same category.In the aorto-femoral model, which has a larger diseased region, we find different mean values, with non-overlapping 99.7% confidence regions, for the local QoIs near the diseased region in the healthy and diseasedmodels (Figure 12c and 12d). We also observe a large difference in the percentage of vessel wall surface areawith low WSS, below a 5 g/cm · s threshold (Figure 12a). Though it is a function of WSS, this QoI is moreof a global quantity as it incorporates the WSS at all vessels of the model. However, in the coronary model,which has a more localized diseased region, the differences in QoI values are only noticeable near the diseasedregion (Figure 13), and most distinct for the local QoIs. As global quantities should not be significantlyaffected by a locally diseased region, this follows our expectations. Within the model geometries, we do seeresults consistent with clinical expectations for each diseased case. For example, in the diseased AAA modelwe see lower wall shear stresses due to the aneurysm and in stenosed coronary model, near the stenosis, wesee increased wall shear stress, both consistent with clinical expectations.Using the data from the pilot run, each QoI value is shown along with its estimated 99.7% (unnormalized)confidence interval (i.e. µ ± σ , where σ is the standard deviation of the estimator value). Recall we use nCIsto facilitate comparisons between model geometries. In general, the nCIs of the diseased model QoIs areslightly larger (less accurate) than the healthy models. As diseased models require better local resolution,and therefore more three-dimensional simulations, to resolve quantities in the diseased regions, we expectedlarger nCI values to arise from the pilot run for the diseased model than for the same QoIs and pilot runsize of a healthy model.The specific QoIs shown were chosen to be representative of each of our four main QoI categories. Thediseased AAA model is shown overlaid on the healthy aorto-femoral model (Figure 12), while the healthycoronary model is overlaid on the diseased (stenosed) coronary model, with the stenosed region enlargedin the inset (Figure 13). The highlighted strips of interest are colored by the local TAWSS of the diseasedmodels. The QoI examining the regions of model surface area with low WSS (Figure 12a), is generally ofinterest in cardiovascular hemodynamics as such regions can be predictive of thrombosis [83].We can also compare the healthy and diseased geometries in the same way we compared the MLMFestimators with and without zero-dimensional models (subsection 5.3). We again averaged the nCIs fromeach QoI category for the healthy and diseased models, computed the extrapolated cost of converging toa nCI of 0.01, and predicted the nCI obtained subject to a computational budget constraint (Figure 14).For both the aorto-femoral and coronary models, we see that the diseased QoIs have similar nCIs to thehealthy QoIs across all categories (Figures 14a and 14d). In the aorto-femoral model, which we discussed asa diseased case with slightly more global effects, we do see that most categories have slightly worse averaged27 odelTAWSS ModelPressureQoIs5.010.0 T A W SS ( g / c m · s ) P r e ss u r e ( mm H g ) ModelTAWSS ModelPressureQoIs5.010.015.0 T A W SS ( g / c m · s ) P r e ss u r e ( mm H g ) ModelTAWSS ModelPressureQoIs5.010.015.0 T A W SS ( g / c m · s ) P r e ss u r e ( mm H g ) OutletFlow OutletPressureQoIs0.05.010.0 F l o w ( c m / s ) P r e ss u r e ( mm H g ) OutletFlow OutletPressureQoIs0.05.010.0 F l o w ( c m / s ) P r e ss u r e ( mm H g ) th = 5 g/cm · s QoI20.030.0 % A r e a w i t h W SS < t h (a)(c)(e) (b)(d)(f) Healthy Model Diseased Model
Figure 12: Comparison of estimators for healthy and diseased model QoIs. Model is shown with strips of interest colored byTAWSS, with the diseased model including the AAA superimposed over the healthy model. Plot (a) shows the percentage ofthe model surface area with WSS below a threshold value of 5 g/cm · s . Plots (b)-(d) show the spatially-averaged TAWSSand time-averaged pressure in (b) the strip in the descending thoracic aorta, (c) the strip superior to (above) the AAA in theabdominal aorta, and (d) the strip inferior to (below) the AAA in the abdominal aorta. Plots (e) and (f) show the spatial-and time-averaged flow and pressure QoIs at the (e) right internal iliac and (f) left internal iliac outlets. (For interpretation ofthe references to color in this figure legend, the reader is referred to the web version of this article.) nCI for the diseased model. When comparing the expected cost of extrapolation to a nCI of 0.01, we see thatthe diseased model QoIs are almost always more expensive to converge than the healthy model QoIs for theaorto-femoral model (Figure 14b). For the coronary model, many of the QoI nCIs were already below 0.01after the pilot run, which is why the extrapolated cost is identical here; it is merely the cost of the pilot run(Figure 14e). Finally, the projected nCIs for a fixed computational budget are better (lower) or very similarfor the healthy model than for the diseased models across all QoI categories of the aorto-femoral model(Figure 14c). Again, the coronary model shows similar results between the two cases for the global QoIs,but with a smaller nCI for the healthy model when considering local QoIs (Figure 14f). For the local QoIs,this is likely explained by the complex flow features seen in the diseased models which rely more heavily on28 ealthy Model Diseased Model OutletFlow OutletPressureQoIs0.480.490.50 F l o w ( c m / s ) P r e ss u r e ( mm H g ) OutletFlow OutletPressureQoIs0.460.470.48 F l o w ( c m / s ) P r e ss u r e ( mm H g ) ModelTAWSS ModelPressureQoIs5.010.0 T A W SS ( g / c m · s ) P r e ss u r e ( mm H g ) ModelTAWSS ModelPressureQoIs5.010.0 T A W SS ( g / c m · s ) P r e ss u r e ( mm H g ) ModelTAWSS ModelPressureQoIs10.020.030.0 T A W SS ( g / c m · s ) P r e ss u r e ( mm H g ) ModelTAWSS ModelPressureQoIs10.020.0 T A W SS ( g / c m · s ) P r e ss u r e ( mm H g ) (a)(c)(e) (b)(d)(f) Figure 13: Comparison of estimators for healthy and diseased model QoIs. Model is shown with strips of interest colored byTAWSS, with the healthy model superimposed over the diseased (stenosed) model. The stenosed region is shown enlargedin the inset image. Plots (a), (b), (d), and (f) show the spatially-averaged TAWSS and time-averaged pressure in (a) thestrip superior to (above) the left and right coronary artery branches, (b) the strip preceding the stenosis of the left anteriordescending (LAD) coronary artery (d) the strip within the peak stenosed region, and (f) the strip following the stenosis. Plots(c) and (e) show the spatial- and time-averaged flow and pressure QoIs at the (c) right coronary artery (RCA) and (e) stenosedLAD outlets. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of thisarticle.)
HF simulations to resolve.
In subsection 5.4, we hinted at the difference between estimators for local and global QoIs. We examinethese differences more fully here, using the healthy aorto-femoral and coronary models. In general, we seemuch better (smaller) nCI values for global QoIs than for local QoI estimators (Figure 15).We first observe a general nCI trend between local and global QoIs by examining average nCI for eachQoI category (Figures 15a and 15d for the aorto-femoral and coronary models, respectively). As before,these values are obtained by computing the sample average and variance of the estimator nCIs from the pilotrun for all QoIs of a certain type. The local QoIs have nCIs more than four times larger than the global QoInCIs calculated from the same pilot run of simulations. This relates to the ability of lower fidelity modelsto capture the change in global QoIs resulting from changes in the uncertain model parameters. Lowercorrelations and less accurate estimators are expected for local QoIs such as local flow indicators (as foundin our coronary benchmark) or for pathological anatomies. We also see that flows are, in general, the mostaccurate, followed by pressures.Since constraining our estimators to be calculated from the fixed pilot run is similar to a budget con-straint, we also expect to see better global than local estimator nCIs when calculating estimators with a fixedcomputational budget. This supports our findings from previous sections. We expect less computationalexpense to achieve a target nCI level for global QoIs than for local QoIs (Figures 14b and 14e).29orto-Femoral Models O u t l e t F l o w s O u t l e t P r e ss u r e s B r a n c h P r e ss u r e s S tr i p P r e ss u r e s B r a n c h T A W SS s S tr i p T A W SS s QoIs0 . . . . n C I ( p il o tr un ) (a) O u t l e t F l o w s O u t l e t P r e ss u r e s B r a n c h P r e ss u r e s S tr i p P r e ss u r e s B r a n c h T A W SS s S tr i p T A W SS s QoIs0500100015002000 E x tr a p o l a t e d c o s t( fi x e dn C I ) (b) O u t l e t F l o w s O u t l e t P r e ss u r e s B r a n c h P r e ss u r e s S tr i p P r e ss u r e s B r a n c h T A W SS s S tr i p T A W SS s QoIs0 . . . . . n C I ( fi x e d H F c o s t) (c) Coronary Models O u t l e t F l o w s O u t l e t P r e ss u r e s B r a n c h P r e ss u r e s S tr i p P r e ss u r e s B r a n c h T A W SS s S tr i p T A W SS s QoIs0 . . . . n C I ( p il o tr un ) (d) O u t l e t F l o w s O u t l e t P r e ss u r e s B r a n c h P r e ss u r e s S tr i p P r e ss u r e s B r a n c h T A W SS s S tr i p T A W SS s QoIs0200400600800 E x tr a p o l a t e d c o s t( fi x e dn C I ) (e) O u t l e t F l o w s O u t l e t P r e ss u r e s B r a n c h P r e ss u r e s S tr i p P r e ss u r e s B r a n c h T A W SS s S tr i p T A W SS s QoIs0 . . . . . n C I ( fi x e d H F c o s t) (f) Healthy Model Diseased Model
Figure 14: Comparing the performance of 3D-1D-0D MLMF estimators from healthy and diseased geometries for the (top)aorto-femoral and (bottom) coronary models. Averaged nCIs for various categories of QoIs calculated after a fixed pilot run of250 simulations are shown in plots (a) and (d). Averaged projected extrapolated cost (in units of the equivalent number of 3Dfine simulations) to obtain a nCI ≤ .
01 are shown in plots (b) and (e). Averaged projected nCI values with the computationalbudget for all HF simulations equivalent to the cost of 50 3D fine simulations are shown for various QoI categories in plots (c)and (f).
Examining specific QoIs reveals the range of nCI values within each QoI category. Identifying the QoIswith the best (smallest) and worst (largest) nCI in each category, we affirm the finding that global QoIshave better best-case and worst-case nCIs than the local QoIs (Figures 15b and 15e versus 15c, and 15f,respectively). More interestingly, we examine the spread in nCI values. The WSS values are not only leastaccurate, but also differ the most between the best and worst case nCI (Figures 15c and 15f). The flowsare most accurate, as seen in the aggregated trends, and there is also little difference between the best andworst case nCI (Figures 15b and 15e).We use the estimator values and variances from the 3D-1D MLMF scheme here as we want to observedifferences not in converged values of the estimators as in subsection 5.4, but in the nCIs obtained for localand global QoIs. Since the nCIs were larger (worse) after the pilot run of the MLMF 3D-1D scheme, weconsider these here. While these results are presented for the healthy models only, similar trends hold forthe diseased cases. 30orto-Femoral Healthy O u t l e t F l o w s O u t l e t P r e ss u r e s B r a n c h P r e ss u r e s S tr i p P r e ss u r e s B r a n c h T A W SS s S tr i p T A W SS s . . . . . . . .
07 QoIs n C I ( p il o tr un ) (a) O u t l e t F l o w O u t l e t P r e ss u r e B r a n c h P r e ss u r e S tr i p P r e ss u r e F l o w ( c m / s ) P r e ss u r e ( mm H g ) (b) B r a n c h T A W SS S tr i p T A W SS T A W SS ( g / c m · s ) (c) Coronary Healthy O u t l e t F l o w s O u t l e t P r e ss u r e s B r a n c h P r e ss u r e s S tr i p P r e ss u r e s B r a n c h T A W SS s S tr i p T A W SS s . . . . . .
05 QoIs n C I ( p il o tr un ) (d) O u t l e t F l o w O u t l e t P r e ss u r e B r a n c h P r e ss u r e S tr i p P r e ss u r e F l o w ( c m / s ) P r e ss u r e ( mm H g ) (e) B r a n c h T A W SS S tr i p T A W SS T A W SS ( g / c m · s ) (f) Global QoIs Local QoIs
Best nCI in QoI Category Worst nCI in QoI Category
Figure 15: Comparison of estimators for global and local model QoIs for the (top) aorto-femoral and (bottom) coronary models.All values are shown for the healthy models and come from the 3D-1D MLMF scheme. Plots (a) and (d) show the averagednCIs in both local and global QoI categories. Plots (b) and (e) show the QoI estimators with the best and worst nCI in sixglobal QoI categories. Plots (c) and (f) show the QoI estimators with the best and worst nCI in three local QoI categories.
6. Discussion
This study demonstrated the significant computational benefits which arise from utilizing multifidelityestimators for uncertainty quantification of cardiovascular hemodynamics. By incorporating less accurateinexpensive low fidelity models into the uncertainty quantification workflow, significantly improved estimatorconfidence can be obtained for the same cost as traditional UQ methods relying only on high fidelity models.31e further demonstrated that for this application, the MLMF estimators performed best, achieving the bestconfidence intervals at the lowest computational cost. Healthy and diseased models of both aorto-femoraland coronary anatomies all demonstrated the computational advantages of the MLMF methods, whichsupports continued use of these methods in our cardiovascular pipeline. MLMF estimators were shown tobe effective for distinct model geometries, suggesting that the method would also be effective for modelsfrom other anatomic regions (e.g. pulmonary or cerebral). As the MLMF estimators have proven effectivefor large problem dimensionality, future work may include many additional uncertain parameters to handlemodels with more complex anatomy, closed loop LPN boundary conditions, and other FSI models.By quantifying the effect of uncertain parameters from both boundary conditions and model materialand fluid properties on global and local QoIs, the proposed framework demonstrated robustness for differentapplications of interest. As this study was conducted for the purpose of demonstrating the power of theMLMF method with an eye toward future cardiovascular modeling applications and clinical questions, thereare limitations of this study which can be addressed in future work. Uncertain parameters were assumed tobe distributed uniformly according to literature targets. In future work, these uncertain mean values couldbe tuned to patient-specific target values. In this case, a more sophisticated uncertain distribution, such as aGaussian, may more accurately reflect prior knowledge. It would also be possible in the future to assimilateuncertain distributions directly from patient data if given repeated readings of patient-specific targets suchas heart rate and blood pressure.To possibly further reduce the number of high-fidelity simulations, future work can investigate a moreideal distribution of simulations in the pilot study. When using 25 simulations at each level as we did inthis study, we obtain nCIs on the order of 0.01 for flow and pressure QoIs, suggesting we are likely utilizingtoo many 3D simulations. We could reshape the pilot size using tools already included in Dakota. This willresult in higher initial nCIs after the pilot study, but we will eventually obtain nCIs on the order of 0.01with possibly fewer overall high-dimensional simulations (though more low-fidelity simulations) following orextrapolation and convergence simulations.Additional uncertain parameters could also be explored in future studies. The three-dimensional Navier-Stokes equations (1) include at least one parameter we did not explore in this study: the vessel wall parameterfor thickness ζ . From previous work [30], we know that vessel wall thickness uncertainties do propagateto have a significant impact on QoIs such as wall strain, though their impact on other quantities is lesssignificant. These uncertainties should be considered in the future. With that said, model geometry is likelythe source of most of our unaccounted for uncertainties. The models used in this study were constructed usingthe SimVascular workflow. This workflow requires semi-automated lumen segmentation of vessels, whichis an operator-dependent method. This leads to inherent uncertainties in the model geometry which werenot addressed in the current study. Future work should explore the effect of these modeling uncertainties.One possible approach could rely on machine learning techniques for model building to assimilate uncertaindistributions for the vessel lumina.In this study, many quantities of interest were tested, but a more targeted approach could be developedwith a specific clinical question in mind. For instance, QoIs could be concentrated in aneurysms or otherdiseased regions, rather than distributed throughout the model as chosen for the purpose of exploration inthis study. Additionally, quantities of interest were chosen to be single values obtained from temporallyaveraging over an entire cardiac cycle. In clinical applications, we are often interested in time-dependentQoIs. The current methods can be easily extended to provide estimators at multiple time points. This wouldgive full waveforms and provide enhanced clinical information. Similarly, time-averaged model pressure andTAWSS values could be computed along vessel centerlines, with uncertainties computed at each point alongthis centerline, instead of one estimator per model branch or strip of interest.We assumed a consistent parameterization among models for this work. However, in general this maynot be the case. To handle that scenario of inconsistent parameterization, we will explore the possibility ofadopting dimension reduction strategies, for instance an active subspace mapping as described in [84]. Thiswould allow for additional low-fidelity modeling approaches. Recent developments in reduced order modeling(ROM) are one way to expand our low-fidelity modeling approach, as new model generation approaches couldeasily be fit into the MLMF framework as additional fidelity levels, or used to replace one or both of ourlow fidelity models. Integration of ROMs within the multifidelity control variate framework has recently32een demonstrated [85, 86], and could be extended to the MLMF approach. One such ROM method isthe so-called “hierarchical model reduction” approach [87–89], which allows for more flexibility in fidelityand discretization levels through the use of Fourier-like expansions and has been effectively demonstratedfor general cardiovascular problems [90–92] as well as specifically for cardiovascular UQ applications [93].Another method approximated the Navier-Stokes partial differential equations with a reduced basis andneural networks [94]. Such ROMs may help address some of the limitations we see with our current one-dimensional modeling, such as extension to multi-inlet and highly tapered anatomies, and could, for example,enable reduced modeling of more complex models that include ventricular flow. Certain quantities of interestfall outside the bounds of any of our low-fidelity modeling capabilities. New methods such as [94] couldextend our UQ methods to include thrombus modeling or QoIs such as residence time and flow recirculation.In addition to development aimed at addressing such concerns, there is active development on automatedgeneration of the one- and zero-dimensional low fidelity models used in this study, including zero-dimensionalmodels constructed directly from three-dimensional lofted models as the one-dimensional models are. Anautomated process to generate low fidelity models would allow the proposed uncertainty quantificationframework to be implemented more easily across a range of patient data, anatomies, and diseases.
7. Conclusion
As clinicians continue to incorporate computational simulations into patient treatment, it is imperativethat robust methods of uncertainty quantification are readily available. To assist in widespread adoption,these UQ methods need to be adaptable to practical computational budgets and be relatively easy touse. By using MLMF estimators, we have demonstrated enormous cost savings compared to traditionalUQ approaches such as Monte Carlo, which rely only on the highest fidelity three-dimensional models. Byincorporating the existing Dakota framework, our workflow is streamlined and partially automated, avoidingthe need for continual monitoring by the user. As such, estimators with high levels of confidence can beobtained with a reasonable computational budget in a semi-automated workflow. This study demonstratedthe effectiveness of the method on multiple patient-specific model geometries for both healthy and diseasedanatomies. By substituting LF one-dimensional models with zero-dimensional LPNs which only looselyapproximate the HF model geometry, we can obtain even higher estimator confidence improvements for thesame computational cost. If we are interested in global quantities of interest, we require fewer simulationsthan if we are interested in local quantities of interest. This is due to the degraded spatial resolution of theLF models, resulting in local QoIs with reduced correlation between model fidelities. This study examinedthe effect of uncertain boundary conditions and model parameters on a wide range of quantities of interest,demonstrating the robustness of MLMF estimators. Uncertainty quantification is essential for providingclinicians with not only simulated predictions, but the confidence of those predictions. In this study weinvestigated the possibility of employing multifidelity approaches in the future to improve the confidence oftargeted clinical QoIs while using uncertainties assimilated from patient data. Several additional multifidelityapproaches more general than MLMF estimators are now available that will be the subject of future studies.These include generalized approximate control variate approaches [71, 95] and latent variables [96].
Acknowledgments
This work was supported by NIH grants R01-EB018302 and R01-HL123689 (P.I. Alison Marsden), theNational Science Foundation under grant NSF CDSE CBET 1508794, and a Center for Computing Re-search Summer Fellowship at Sandia National Laboratories (Casey Fleeter). This work used computationalresources from the Extreme Science and Engineering Discovery Environment (XSEDE) [97], which is sup-ported by National Science Foundation grant number ACI-1548562, and the Stanford Research ComputingCenter (SRCC). We thank Mahidhar Tatineni for his assistance building Dakota on the Comet cluster andadvice on efficiency in the workflow, which was made possible through the XSEDE Extended Collabora-tive Support Service (ECSS) program [98]. We thank Michael S. Eldred of Sandia National Laborato-ries for his ongoing support with Dakota. We also acknowledge the open source SimVascular project at and the open source Dakota toolkit at https://dakota.sandia.gov .33andia National Laboratories is a multimission laboratory managed and operated by National Technology& Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for theU.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525.This paper describes objective technical results and analysis. Any subjective views or opinions that mightbe expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or theUnited States Government.
Disclosures
The authors have no conflicts of interest to disclose.
References [1] GBD 2017 Causes of Death Collaborators, Global, regional, and national age-sex-specific mortality for 282 causes of deathin 195 countries and territories, 1980-2017: a systematic analysis for the global burden of disease study 2017, Lancet 392(2018) 1736–1788.[2] M. Heron, Deaths: Leading causes for 2016, Natl Vital Stat Rep 67 (2018).[3] C. A. Taylor, T. J. Hughes, C. K. Zarins, Finite element modeling of blood flow in arteries, Comput Methods Appl MechEng 158 (1998) 155–196.[4] E. J. Benjamin, S. S. Virani, C. W. Callaway, A. M. Chamberlain, A. R. Chang, et al., Heart disease and strokestatistics—2018 update: A report from the american heart association, Circulation 137 (2018) e67–e492.[5] C. A. Taylor, T. A. Fonte, J. K. Min, Computational fluid dynamics applied to cardiac computed tomography fornoninvasive quantification of fractional flow reserve: Scientific basis, J Am Coll Cardiol 61 (2013) 2233–2241.[6] A. B. Ramachandra, A. M. Kahn, A. L. Marsden, Patient-specific simulations reveal significant differences in mechanicalstimuli in venous and arterial coronary grafts, J Cardiovasc Transl Res 9 (2016) 279–290.[7] E. Kung, A. Baretta, C. Baker, G. Arbia, G. Biglino, et al., Predictive modeling of the virtual hemi-fontan operation forsecond stage single ventricle palliation: Two patient-specific cases, J Biomech 46 (2013) 423–429.[8] D. Schiavazzi, E. Kung, A. Marsden, C. Baker, G. Pennati, et al., Hemodynamic effects of left pulmonary artery stenosisafter superior cavopulmonary connection: A patient-specific multiscale modeling study, J Thorac Cardiovasc Surg 149(2015) 689–696.e3.[9] A. Verma, M. Esmaily, J. Shang, R. Figliola, J. A. Feinstein, T.-Y. Hsia, A. L. Marsden, Optimization of the assistedbidirectional glenn procedure for first stage single ventricle repair, World J Pediatr Congenit Heart Surg 9 (2018) 157–170.[10] D. Sengupta, A. M. Kahn, J. C. Burns, S. Sankaran, S. C. Shadden, A. L. Marsden, Image-based modeling of hemodynamicsin coronary artery aneurysms caused by kawasaki disease, Biomech Model Mechanobiol 11 (2012) 915–932.[11] N. G. Gutierrez, O. Shirinsky, N. Gagarina, G. Lyskina, R. Fukazawa, et al., Assessment of coronary artery aneurysmscaused by kawasaki disease using transluminal attenuation gradient analysis of computerized tomography angiograms, AmJ Cardiol 120 (2017) 556–562.[12] N. Grande Gutierrez, M. Mathew, B. W. McCrindle, J. Tran, A. M. Kahn, et al., Hemodynamic variables in aneurysmsare associated with thrombotic risk in children with kawasaki disease, Int J Cardiol 281 (2019).[13] W. Yang, A. L. Marsden, M. T. Ogawa, C. Sakarovitch, K. K. Hall, et al., Right ventricular stroke work correlates withoutcomes in pediatric pulmonary arterial hypertension, Pulmonary Circulation 8 (2018).[14] W. Yang, M. Dong, M. Rabinovitch, F. P. Chan, A. L. Marsden, J. A. Feinstein, Evolution of hemodynamic forces inthe pulmonary tree with progressively worsening pulmonary arterial hypertension in pediatric patients, Biomech ModelMechanobiol 18 (2019) 779–796.[15] G.-Y. Suh, A. S. Les, A. S. Tenforde, S. C. Shadden, R. L. Spilker, et al., Quantification of particle residence time inabdominal aortic aneurysms using magnetic resonance imaging and computational fluid dynamics, Ann Biomed Eng 39(2011) 864–883.[16] A. Arzani, S. Shadden, Characterization of the transport topology in patient-specific abdominal aortic aneurysm models,Phys Fluids 24 (2012) 81901.[17] M. Piccinelli, A. Veneziani, D. A. Steinman, A. Remuzzi, L. Antiga, A framework for geometric analysis of vascularstructures: Application to cerebral aneurysms, IEEE Trans Med Imaging 28 (2009) 1141–1155.[18] J. Cebral, F. Mut, J. Weir, C. Putman, Association of hemodynamic characteristics and cerebral aneurysm rupture,AJNR Am J Neuroradiol 32 (2011) 264–270.[19] J. R. Cebral, M. A. Castro, J. E. Burgess, R. S. Pergolizzi, M. J. Sheridan, C. M. Putman, Characterization of cerebralaneurysms for assessing risk of rupture by using patient-specific computational hemodynamics models, AJNR Am JNeuroradiol 26 (2005) 2550–2559.[20] T. J. Hughes, J. Lubliner, On the one-dimensional theory of blood flow in the larger vessels, Math Biosci 18 (1973)161–170.[21] M. Mirramezani, S. Diamond, H. Litt, S. Shadden, Reduced order models for transstenotic pressure drop in the coronaryarteries, J Biomech Eng 141 (2019) 031005–031011.[22] A. Quarteroni, A. Veneziani, C. Vergara, Geometric multiscale modeling of the cardiovascular system, between theoryand practice, Comput Methods Appl Mech Eng 302 (2016) 193–252.
23] M. E. Moghadam, I. E. Vignon-Clementel, R. Figliola, A. L. Marsden, A modular numerical method for implicit 0d/3dcoupling in cardiovascular finite element simulations, J Comput Phys 244 (2013) 63–79.[24] A. Quarteroni, S. Ragni, A. Veneziani, Coupling between lumped and distributed models for blood flow problems, ComputVis Sci 4 (2001) 111–124.[25] F. Migliavacca, G. Pennati, G. Dubini, R. Fumero, R. Pietrabissa, et al., Modeling of the norwood circulation: effects ofshunt size, vascular resistances, and heart rate, Am J Physiol Heart Circ Physiol 280 (2001) H2076–H2086.[26] M. M. Esmaily, F. Migliavacca, I. Vignon-Clementel, T. Hsia, A. Marsden, Optimization of shunt placement for thenorwood surgery using multi-domain modeling, J Biomech Eng 134 (2012) 051002–0510013.[27] D. Schiavazzi, G. Arbia, C. Baker, A. Hlavacek, T. Y Hsia, A. Marsden, I. Vignon-Clementel, T. Modeling Of CongenitalHearts Alliance Mocha Investigators, Uncertainty quantification in virtual surgery hemodynamics predictions for singleventricle palliation, Int J Numer Method Biomed Eng 32 (2015) e02737.[28] S. Sankaran, H. J. Kim, G. Choi, C. A. Taylor, Uncertainty quantification in coronary blood flow simulations: Impact ofgeometry, boundary conditions and blood viscosity, J Biomech 49 (2016) 2540–2547.[29] D. Schiavazzi, A. Doostan, G. Iaccarino, A. Marsden, A generalized multi-resolution expansion for uncertainty propagationwith application to cardiovascular modeling, Comput Methods Appl Mech Eng 314 (2017) 196–221.[30] J. S. Tran, D. E. Schiavazzi, A. M. Kahn, A. L. Marsden, Uncertainty quantification of simulated biomechanical stimuliin coronary artery bypass grafts, Comput Methods Appl Mech Eng 345 (2019) 402–428.[31] P. Chen, A. Quarteroni, G. Rozza, Simulation-based uncertainty quantification of human arterial network hemodynamics,Int J Numer Method Biomed Eng 29 (2013) 698–721.[32] S. Sankaran, A. L. Marsden, The impact of uncertainty on shape optimization of idealized bypass graft models in unsteadyflow, Phys Fluids 22 (2010) 121902.[33] V. Eck, J. Sturdy, L. Hellevik, Effects of arterial wall models and measurement uncertainties on cardiovascular modelpredictions, J Biomech 50 (2017) 188–194.[34] J. S. Tran, D. E. Schiavazzi, A. B. Ramachandra, A. M. Kahn, A. L. Marsden, Automated tuning for parameter identifi-cation and uncertainty quantification in multi-scale coronary simulations, Comput Fluids 142 (2017) 128–138.[35] D. E. Schiavazzi, A. Baretta, G. Pennati, T.-Y. Hsia, A. L. Marsden, Patient-specific parameter estimation in single-ventricle lumped circulation models under uncertainty, Int J Numer Method Biomed Eng 33 (2017) e02799.[36] J. Biehler, W. A. Wall, The impact of personalized probabilistic wall thickness models on peak wall stress in abdominalaortic aneurysms, Int J Numer Method Biomed Eng 34 (2017) e2922.[37] A. D. Marquis, A. Arnold, C. Dean-Bernhoft, B. E. Carlson, M. S. Olufsen, Practical identifiability and uncertaintyquantification of a pulsatile cardiovascular model, Math Biosci 304 (2018) 9–24.[38] A. Brault, L. Dumas, D. Lucor, Uncertainty quantification of inflow boundary condition and proximal arterial stiffness-coupled effect on pulse wave propagation in a vascular network, Int J Numer Method Biomed Eng 33 (2017) e2859.[39] A. Boccadifuoco, A. Mariotti, S. Celi, N. Martini, M. Salvetti, Impact of uncertainties in outflow boundary conditions onthe predictions of hemodynamic simulations of ascending thoracic aortic aneurysms, Comput Fluids 165 (2018) 96–115.[40] D. Xiu, G. Karniadakis, The wiener–askey polynomial chaos for stochastic differential equations, SIAM J Sci Comput 24(2002) 619–644.[41] I. Babuˇska, F. Nobile, R. Tempone, A stochastic collocation method for elliptic partial differential equations with randominput data, SIAM J Numer Anal 45 (2007) 1005–1034.[42] N. Metropolis, S. Ulam, The Monte Carlo method, J Am Stat Assoc 44 (1949) 335–341.[43] M. Giles, Multilevel Monte Carlo methods., Acta Numer 24 (2015) 259–328.[44] B. Peherstorfer, K. Willcox, M. Gunzburger, Optimal model management for multifidelity Monte Carlo estimation., SIAMJ Sci Comput 38 (2016) A3163–A3194.[45] F. Nobile, F. Tesei, A multi level monte carlo method with control variate for elliptic pdes with log-normal coefficients,Stochastic Partial Differential Equations: Analysis and Computations 3 (2015) 398–444.[46] G. Geraci, M. Eldred, G. Iaccarino, A multifidelity control variate approach for the multilevel Monte Carlo technique,Center for Turbulence Research, Stanford University, 2015, pp. 169–181.[47] G. Geraci, M. S. Eldred, G. Iaccarino, A multifidelity multilevel Monte Carlo method for uncertainty propagation inaerospace applications, in: 19th AIAA Non-Deterministic Approaches Conference, American Institute of Aeronautics andAstronautics, Grapvine, Texas, 2017.[48] H. R. Fairbanks, A. Doostan, C. Ketelsen, G. Iaccarino, A low-rank control variate for multilevel monte carlo simulationof high-dimensional uncertain systems, Journal of Comput Phys 341 (2016).[49] A. Updegrove, N. Wilson, J. Merkow, H. Lan, A. Marsden, S. Shadden, SimVascular: An open source pipeline forcardiovascular simulation, Ann Biomed Eng 45 (2016) 525–541.[50] B. Adams, M. Ebeida, M. Eldred, G. Geraci, J. Jakeman, et al., Dakota, a multilevel parallel object-oriented framework fordesign optimization, parameter estimation, uncertainty quantification, and sensitivity analysis: Version 6.6 user’s manual,Sandia Technical Report SAND2014-4633 (2014). Updated May 2017.[51] L. Formaggia, F. Nobile, A. Quarteroni, A. Veneziani, Multiscale modelling of the circulatory system: a preliminaryanalysis, Comput Vis Sci 2 (1999) 75–83.[52] I. E. Vignon-Clementel, C. A. Figueroa, K. E. Jansen, C. A. Taylor, Outflow boundary conditions for three-dimensionalfinite element modeling of blood flow and pressure in arteries, Comput Methods Appl Mech Eng 195 (2006) 3776–3796.[53] M. Esmaily Moghadam, Y. Bazilevs, T.-Y. Hsia, I. E. Vignon-Clementel, A. L. Marsden, Modeling of Congenital HeartsAlliance (MOCHA), A comparison of outlet boundary treatments for prevention of backflow divergence with relevance toblood flow simulations, Comput Mech 48 (2011) 277–291.[54] M. Esmaily-Moghadam, Y. Bazilevs, A. L. Marsden, A bi-partitioned iterative algorithm for solving linear systems arising rom incompressible flow problems, Comput Methods Appl Mech Eng 286 (2015) 40–62.[55] K. E. Jansen, C. H. Whiting, G. M. Hulbert, A generalized- α method for integrating the filtered navierstokes equationswith a stabilized finite element method, Comput Methods Appl Mech Eng 190 (2000) 305–319.[56] M. Esmaily-Moghadam, Y. Bazilevs, A. L. Marsden, A new preconditioning technique for implicitly coupled multidomainsimulations with applications to hemodynamics, Comput Mech 52 (2013) 1141–1152.[57] J. Seo, D. E. Schiavazzi, A. L. Marsden, Performance of preconditioned iterative linear solvers for cardiovascular simulationsin rigid and deformable vessels, Comput Mech (2019).[58] Y. Chatzizisis, A. Coskun, M. Jonas, E. Edelman, C. Feldman, P. Stone, Role of endothelial shear stress in the naturalhistory of coronary atherosclerosis and vascular remodeling: molecular, cellular, and vascular behavior, J Am Coll Cardiol49 (2007) 2379–2393.[59] A. Figueroa, I. Vignon-Clementel, K. Jansen, T. Hughes, C. Taylor, A coupled momentum method for modeling bloodflow in three-dimensional deformable arteries, Comput Methods Appl Mech Eng 195 (2006) 5685–5706.[60] L. Euler, Principia pro motu sanguinis per arterias determinando, Opera posthuma mathematica et physica anno 1844detecta 2 (1775) 814–823. Ediderunt PH Fuss et N Fuss Petropoli; Apund Eggers et Socios.[61] I. Vignon, C. Taylor, Outflow boundary conditions for one-dimensional finite element modeling of blood flow and pressurewaves in arteries, Wave Motion 39 (2004) 361–374.[62] I. Vignon, A Coupled Multidomain Method for Computational Modeling of Blood Flow, Ph.D. thesis, Stanford University,2006.[63] J. Wan, B. Steele, S. A. Spicer, S. Strohband, G. R. Feijo, T. J. Hughes, C. A. Taylor, A one-dimensional finite elementmethod for simulation-based medical planning for cardiovascular disease, Comput Methods Biomech Biomed Engin 5(2002) 195–206.[64] B. Steele, C. Taylor, J. Wan, J. Ku, T. Hughes, In vivo validation of a one-dimensional finite element method forsimulation-based medical planning for cardiovascular bypass surgery, volume 1, 2001, pp. 120–123.[65] V. Miliˇsi´c, A. Quarteroni, Analysis of lumped parameter models for blood flow simulations and their relation with 1Dmodels, ESAIM Math Model Numer Anal 38 (2004) 613–632.[66] F. Migliavacca, R. Balossino, G. Pennati, G. Dubini, T.-Y. Hsia, et al., Multiscale modelling in biofluidynamics: Appli-cation to reconstructive paediatric cardiac surgery, J Biomech 39 (2006) 1010–1020.[67] C. Corsini, D. Cosentino, G. Pennati, G. Dubini, T.-Y. Hsia, F. Migliavacca, Multiscale models of the hybrid palliationfor hypoplastic left heart syndrome, J Biomech 44 (2011) 767–770.[68] R. Pasupathy, B. W. Schmeiser, M. R. Taaffe, J. Wang, Control-variate estimation using estimated control means., IIETrans 44 (2012) 381–385.[69] L. Ng, K. Willcox, Multifidelity approaches for optimization under uncertainty., Int J Numer Methods Eng 100 (2014)746–772.[70] R. Y. Rubinstein, R. Marcus, Efficiency of multivariate control variates in monte carlo simulation, Operations Research33 (1985) 661–677.[71] A. A. Gorodetsky, G. Geraci, M. Eldred, J. D. Jakeman, A Generalized Approximate Control Variate Framework forMultifidelity Uncertainty Quantification, J Comput Phys (2020) 109257.[72] D. C. Maniaci, A. L. Frankel, G. Geraci, M. L. Blaylock, M. S. Eldred, Multilevel uncertainty quantification of awind turbine large eddy simulation model, in: 6th European Conference on Computational Mechanics—7th EuropeanConference on Computational Fluid Dynamics, International Centre for Numerical Methods in Engineering, Glasgow, UK,2018, pp. 2747–2758.[73] C. M. Fleeter, G. Geraci, D. E. Schiavazzi, A. M. Kahn, M. S. Eldred, A. L. Marsden, Multilevel multifidelity approachesfor cardiovascular flow under uncertainty, in: Sandia Center for Computing Research Summer Proceedings 2017, A.D.Baczewski and M.L. Parks, eds., volume Technical Report SAND2018-2780O, Sandia National Laboratories, 2018, pp.27–50.[74] D. Schiavazzi, C. M. Fleeter, G. Geraci, A. L. Marsden, Multifidelity approaches for cardiovascular hemodynamics, in:6th European Conference on Computational Mechanics—7th European Conference on Computational Fluid Dynamics,International Centre for Numerical Methods in Engineering, Glasgow, UK, 2018, pp. 2759–2770.[75] C. Fleeter, Dakota-Simvascular Cardiovascular UQ Interface, 2019. doi: .[76] K. C. Kent, Abdominal aortic aneurysms, N Engl J Med 371 (2014) 2101–2108.[77] G. D. Giannoglou, A. P. Antoniadis, Y. S. Chatzizisis, G. E. Louridas, Difference in the topography of atherosclerosis inthe left versus right coronary artery in patients referred for coronary angiography, BMC Cardiovasc Disord 10 (2010) 26.[78] J. Peir´o, A. Veneziani, Reduced models of the cardiovascular system, Springer Milan, Milano, 2009, pp. 347–394.[79] Y. Zhou, G. Kassab, S. Molloi, On the design of the coronary arterial tree: a generalization of Murray’s law, Phys MedBiol 44 (1999) 2929–2945.[80] A. S. Les, S. C. Shadden, C. A. Figueroa, J. M. Park, M. M. Tedesco, et al., Quantification of hemodynamics in abdominalaortic aneurysms during rest and exercise using magnetic resonance imaging and computational fluid dynamics, AnnBiomed Eng 38 (2010) 1288–1313.[81] J. Coogan, J. Humphrey, C. Figueroa, Computational simulations of hemodynamic changes within thoracic, coronary,and cerebral arteries following early wall remodeling in response to distal aortic coarctation, Biomech Model Mechanobiol12 (2013) 79–93.[82] S. Roccabianca, C. Figueroa, G. Tellides, J. Humphrey, Quantification of regional differences in aortic stiffness in theaging human, J Mech Behav Biomed Mater 29 (2014) 618–634.[83] B. A. Zambrano, H. Gharahi, C. Lim, F. A. Jaberi, J. Choi, W. Lee, S. Baek, Association of intraluminal thrombus,hemodynamic forces, and abdominal aortic aneurysm expansion using longitudinal ct images., Ann Biomed Eng 44 (2015) Communications in Computer and Information Science , SpringerVerlag, Germany, 2016, pp. 3–13., SpringerVerlag, Germany, 2016, pp. 3–13.