Geometric Uncertainty in Patient-Specific Cardiovascular Modeling with Convolutional Dropout Networks
Gabriel Maher, Casey Fleeter, Daniele Schiavazzi, Alison Marsden
GGeometric Uncertainty in Patient-Specific Cardiovascular Modeling withConvolutional Dropout Networks
Gabriel D. Maher , Casey M. Fleeter , Daniele E. Schiavazzi , Alison L. Marsden ∗ Institute for Computational and Mathematical Engineering, Stanford University, Stanford, CA, USA Department of Applied and Computational Mathematics and Statistics, University of Notre Dame, Notre Dame, IN, USA Departments of Pediatrics and Bioengineering, Stanford University, Stanford, CA, USA
Abstract
We propose a novel approach to generate samples from the conditional distribution of patient-specific cardiovas-cular models given a clinically aquired image volume. A convolutional neural network architecture with dropoutlayers is first trained for vessel lumen segmentation using a regression approach, to enable Bayesian estimation ofvessel lumen surfaces. This network is then integrated into a path-planning patient-specific modeling pipeline togenerate families of cardiovascular models. We demonstrate our approach by quantifying the effect of geometricuncertainty on the hemodynamics for three patient-specific anatomies, an aorto-iliac bifurcation, an abdominalaortic aneurysm and a sub-model of the left coronary arteries. A key innovation introduced in the proposed ap-proach is the ability to learn geometric uncertainty directly from training data. The results show how geometricuncertainty produces coefficients of variation comparable to or larger than other sources of uncertainty for wallshear stress and velocity magnitude, but has limited impact on pressure. Specifically, this is true for anatomiescharacterized by small vessel sizes, and for local vessel lesions seen infrequently during network training.
Results from cardiovascular models are affected by a number of uncertainty sources, including material properties,image-data resolution, and boundary condition selection to match clinical target data. A rigorous determination ofsimulation uncertainty and the development of numerical approaches to efficiently quantify its effects on patient-specific models is necessary to increase clinical adoption of simulation tools and improve their effectiveness for earlytreatment planning and non-invasive diagnostics.Prior studies have investigated a range of methods for characterizing the influence of specific sources of uncertaintyin cardiovascular modeling [5, 37]. To reduce the computational burden with respect to Monte Carlo sampling, astochastic collocation approach was proposed in [29], while a multi-fidelity approach was proposed by [1] in thecontext of mechanical stress analysis of abdominal aortic aneurysms. Additionally, a generalized polynomial chaosexpansion is presented in [24] and applied to two pathological anatomies, i.e., an abdominal aortic aneurysm andan arteriovenous fistula, respectively. A generalized multiresolution expansion for uncertainty quantification wasdeveloped in [30] to better handle uncertainty in the presence of non-smooth stochastic responses, while mitigatingthe exponential complexity of multi-dimensional multi-wavelet refinement. Combined uncertainty in vessel wallmaterial properties and hemodynamics are investigated in [40] for several patient-specific models of coronary arterybypass grafting, leveraging a novel submodeling approach to focus the analysis only on venous and arterial bypassgrafts. This and other studies focusing on the coronary circulation, (see, e.g., [33]) contributed to show a loosecoupling between hemodynamics and wall mechanics for such anatomies. One dimensional models have been usedto better understand main pulmonary artery pressure uncertainty in mice due to material property and imagesegmentation uncertainty [3, 4].Finally, multifidelity simulations based on approximate control variate variance reduction in Monte Carlo sam-pling, were thorougly analyzed in the context of deformable cardiovascular models in [6].While the above contributions focus on the propagation of uncertainty from model inputs to outputs, an end-to-end (or clinical data to simulation results) uncertainty analysis pipeline is proposed in [31] in the context of virtualstage II single ventricle palliation surgery. Additionally, the solution of inverse problems is discussed in [39], whereautomated Bayesian estimation is applied to tune close-loop boundary condition parameters for patient-specificmulti-scale models of the coronary circulation, in order to match a number of non-invasive clinical measurements. ∗ corresponding author: [email protected] a r X i v : . [ q - b i o . Q M ] S e p he vast majority of studies in the literature focus on uncertainty in the boundary conditions and mechanicalproperties of the vascular walls. A third major source is geometric uncertainty which results from errors and operatorsubjectivity in vessel segmentation from image data, which constitutes a fundamental step in the generation of cardio-vascular models. Acquisition of medical image volumes is inherently noisy, has limitations related to the achievableresolution as well as artifact, motion, and aliasing errors. Construction of patient-specific model geometries fromimage volumes is therefore affected by image uncertainty. In the literature, analysis of the effects of geometricaluncertainty on the results of high-fidelity cardiovascular models has remained elusive due to the complexity of as-sembling end-to-end pipelines for automatic model generation and analysis. As discretization approaches invariablyrequire the geometry to be represented through a discrete surface mesh, a popular technique within the biomechanicscommunity is mesh morphing [36]. Other methods have focused on modeling the variation of geometry via segmen-tation approaches. For example, for a given input image the STAPLE algorithm generates a distribution of possiblesegmentations, but requires a set of ground-truth segmentations as input. Gaussian processes have also been used tomodel preexisting segmentation variation in [16]. Segmentation priors and multivariate sensitivity analysis proveduseful for segmentation variability estimation in [13], however it is unclear how to extend the method to multiplesimultaneous images.Only a few studies consider the effect of geometrical uncertainty for cardiovascular models. Sensitivity of hemo-dynamics to geometry variation in patient-specific cerebral aneurysms was investigated in [8], which considered twomodel samples generated using heuristic smoothing techniques. Manual segmentation uncertainty was shown to havevarying influence on FFR-CT calculations in [43], where uncertainty depended on the mean FFR-CT value. In anaortic flow simulation, geometric uncertainty was shown to be a dominant factor when compared to computationalfluid dynamics (CFD) model parameter uncertainty and boundary condition uncertainty [46]. Geometric uncertaintywas also investigated in [26] for coronary artery simulations, and obtained through local perturbations of an idealizedstenosis model. Effects on entire cardiovascular models were investigated in [25, 27, 28] by perturbing the area andsurface points of selected vessel segments using a spatial Gaussian function with uniform parameterization. Thevariation in geometry was found to produce sensitivities of up to 10% in simulated FFR-CT measurements. Wewould like to point out that geometric uncertainty is assumed a-priori in the above studies, instead of being directly learned from the image data.More recently, Bayesian Neural Networks, neural networks that are able to learn uncertainty from data, have beenincreasingly adopted in applications where it is crucial to quantify confidence in predictions [7, 12]. In particular, [7]showed that augmenting neural networks with dropout layers enables them to learn uncertainty from the trainingdataset. In the medical imaging field dropout networks have been used to model segmentation uncertainty for MRIvolumes [23]. In particular, the network’s prediction uncertainty was found to be a useful marker for detecting humanexpert prediction error.In this work, we use Bayesian deep learning to develop a cardiovascular model generation technique that learnsthe geometry distribution from a dataset of existing geometries and images. We then use this network along withMonte Carlo sampling and numerical blood flow simulation to characterize the change in model outputs due togeometric uncertainty.In section 2 we discuss our dropout network architecture and path-planning cardiovascular model generationprocess. Sections 3 and 4 provide an overview of the anatomical benchmarks we selected and the results we obtained.Finally, Section 5 contains a discussion and 6 presents our conclusions. Given a medical image volume X and a set of vessel pathlines V , our method produces samples from the distributionof patient-specific cardiovascular models, Y ∼ P ( Y | X , V , z ), compatible with both the image data, pathline anda collection of latent random variables z . To better explain how this is accomplished, we summarize in Figure 1 atypical two-dimensional segmentation or path-planning approach, i.e., a widely used method to generate anatomicalsurfaces developed in a prior work [19] and based on the cardiovascular model format developed for SimVascular[41] and the Vascular Model Repository (VMR) [45]. This requires to first define a vessel pathline (we will usethe term centerline interchangeably) by connecting user-specified point locations inside the lumen of the vessel ofinterest (e.g., the aortic arch in Figure 1). The tangent vector to this centerline is then used to generate a continuouscollection of local 2D images slices. Two-dimensional vessel lumen segmentation on this slice is accomplished througha parametric estimator, trained using a large collection of 2D cross-sectional images and corresponding ground truthlumen boundary. The resulting two-dimensional segments, representing the intersection between the lumen wallsurface and the cross-section plane, are then lofted into a three-dimensional lumen surface and the final modelgenerated by boolean union of multiple vessels. 2igure 1: Proposed model building pipeline. (a) Image data and vessel pathline are supplied by the user. (b)
Path information is used to extract local 2D cross-sectional images in the plane orthogonal to the vessel path. (c)
Two-dimensional images are extracted along vessel pathlines and fed to the CNN as inputs. (d) and (e)
Theproposed CNN processes the cross-sectional images and directly outputs an array of point coordinates, characterizinga two-dimensional lumen segmentation. (f )
The collection of two-dimensional points is transformed back to three-dimensional coordinates on the image volume. (g)
The cross-sectional segmentations are lofted along the pathlineto form the final lumen surface.Note how our vessel lumen estimator also depends on a collection of latent variables z , where two different real-izations z (1) , z (2) ∼ p ( z ) will produce two distinct anatomical surfaces. This way, we can naturally generate familiesof cardiovascular models Y := { Y (1) , . . . , Y ( n ) } , where Y i ∼ P ( Y | X , V , z ) , i = 1 , . . . , n . Given this ability to sample from distributions of cardiovascular anatomies, we adopt Monte Carlo sampling and a computational fluid dynam-ics (CFD) solver to estimate the changes in hemodynamics induced by geometrical uncertainty in the segmentedanatomy. Specifically, each sample { Y (1) , . . . , Y ( n ) } provides a computational domain where we numerically solvethe Navier-Stokes equations using the finite element method, as further discussed in Section 2.3. This Monte Carloprocess is described in Algorithm 1 and outlined in Figure 2. Algorithm 1
Monte Carlo Sampling of Geometrically Uncertain Hemodynamic Solutionsmedical image volume X vessel pathlines V element size h parametric estimator m θ x := { x , ..., x n } ← extract ( X , V ) U := {} for i = 1 , ..., K doz ∼ P ( z ) Y := { ˆ y , ..., ˆ y n } ← { m θ ( x ; z ) , ..., m θ ( x n ; z ) } Y i ← model ( Y , V ) Y hi ← mesh ( Y i , h ) U hi ← simulate ( Y hi ) U ← U ∪ U hi end for return U extract ( X , V ) creates a collection of cross-section images by slicing the image X orthogonal to the vessel pathlines V . model () generates a cardiovascular model by lofting the two-dimensional lumen segments along the vessel pathline, and merging multiplevessels together by boolean union. mesh ( Y i , h ) generates a tetrahedral mesh of the domain Y i , using an element size h (see details in [42, 41]). simulate ( Y h ) computes the solution of the Navier-Stokes equations on Y h , with appropriately chosen initial and boundary conditions. In this section we provide a more formal description of the path-planning process (Figure 3). The first inputconsists of a gray-scale medical image volume with
H, W and D voxels in the axial, sagittal and coronal direction,respectively, i.e. X ∈ R H × W × D . The second input is a single pathline which consists of a collection of spline3igure 2: Generation of geometrically uncertain cardiovascular model solutions following Algorithm 1segments V := { v ( s ) , ..., v N v ( s ) } . Each segment, v i ( s ) : [0 , → R , i = 1 , . . . , N v is a function that maps therelative length along the segment to three-dimensional locations within the image volume, and is obtained throughspline interpolation from a collection of user-specified points locations [41]. A collection of N s local two-dimensionalcross-sectional images of the vessel lumen are then extracted for each pathline at the discrete image space locations S := { v i ( s ) , ..., v i Ns ( s N s ) } and the local tangent and normal vectors to the pathline at point v i j ( s k ) are used toconstruct a planar grid where the voxel intensities from X are interpolated, to create a gray-scale cross-sectionalimage of the vessel lumen x i , i = 1 , . . . , N s . Repeating this process for all selected N s locations along the pathlineproduces a set of two-dimensional images x := { x , ..., x N s } .Realizations from a Bernoulli random vector z and the images x i , i = 1 , . . . , N s constitute the inputs to anartificial neural network designed to produce two-dimensional segmentations of the form y := { y , ..., y N s } = { m θ ( x ; z ) , ..., m θ ( x N s ; z ) } , (1)at all the N s locations along the pathline. It is important to note how the same realization from the Bernoulli vector z is used to set the dropout layer in the network across all the N s segmentation instances for each anatomical surfacerealization. In contrast to generating the segmentations with independent dropout vectors, our process ensures thatthe same network weights are used to segment all the x i , i = 1 , . . . , N s images, leading to a consistent bias across thewhole cardiovascular model geometry. Finally, the normal and tangent vectors to the pathline are used to re-orientthe 2D lumen segments back to image space (Figure 3c), where they are interpolated and joined together [42] to forma triangular surface mesh of the full cardiovascular model [41]. A volumetric finite element Delaunay triangulationis then generated using TetGen [35] for finite element analysis. For the parametric vessel lumen estimator, m θ ( x ; z ), we use a convolutional neural network. In particular, m θ mapsthe input 2D gray-scale image slice, x ∈ R H × H , to a vector of K normalized radii y ∈ [0 , K ⊂ R K . The radiicorrespond to the distance of the vessel lumen from the center of the image along rays oriented according to angularintervals φ := { φ , ..., φ K } . This allows the radius y ji ∈ y i , i = 1 , . . . , N s , j = 1 , . . . , K to be converted to a singlelocation in the cross-sectional slice x i using the expression p ji = (cid:16) y ji H cos φ j , ˆ y ji H sin φ j (cid:17) , i = 1 , . . . , N s , j = 1 , . . . , K, (2)and the full lumen segmentation from image slice x i in the set of points p i := { p i , ..., p Ki } . Even though the literaturehas witnessed an explosion in new layouts and arrangements in recent years [17, 32], a CNN generally consist of a4 a) (b) (c) (d) (e) Figure 3: Cardiovascular model construction workflow used in SimVascular [41]. Starting from (a)
Image data, (b) pathlines are manually generated by the users, (c) two-dimensional lumen segmentations are generated at each crosssection x i , i = 1 , . . . , N s , (d) the entire vessel lumen surface is reconstructed by lofting, (e) a Boolean union ofmultiple vessels is meshed to generate a 3D cardiovascular model.collection of layers, each applying a mathematical operation, such as a linear transformation or convolution, followedby an elementwise nonlinear activation. In our case, the transformation in layer l is expressed as o ( l ) = m ( l ) ( a ( l − ; Θ ( l ) ) , a ( l ) = g ( l ) ( o ( l ) ) , (3)where m ( l ) ( · ) is the specific mathematical transformation occurring through layer l , a represents the generic inputvector from the previous layer and output to the next one, and g ( l ) ( · ) the selected non linear activation. Thelearnable parameters for the l -th layer are denoted by Θ ( l ) ⊂ θ . In this study, we employ a CNN combining denseand convolutional layers. Dense layers operate on vector inputs and outputs through the linear transformation o ( l ) = Θ ( l ) a ( l − + b ( l ) , (4)where Θ ( l ) and b ( l ) are a weight matrix and bias term, respectively. The convolutional layers instead transform athird order tensor input using o ( l ) ijk = (cid:88) o (cid:88) p (cid:88) q Θ ( l ) opqk a ( l − i + o,j + p,q + b ( l ) , (5)where Θ ( l ) is a fourth order tensor of trainable weights. Note how the output from convolutional layers are flattenedinto one-dimensional vectors before being fed to dense layers.The activation functions g ( l ) ( · ) allow the neural network to learn nonlinear relationships in the data [10] anddetermine the types of output it can produce. Since the output in our case is a vector of radii in [0 ,
1] we use theelementwise sigmoid activation function g ( l ) ( x ) = 11 + e − x . (6)For the intermediate layers we instead use Leaky Rectified Linear Units (Leaky-RELU) because they avoid theproblem with vanishing gradients when optimizing the network weights using gradient-descent [18] g ( l ) ( x ) = (cid:40) x, x > α · x, x ≤ . (7)We augment our convolutional network to sample from the distribution of vessel lumens for a given image by addingdropout layers to the network, which sets the outputs of the previous layer to zero through Hadamard (elementwise)products by a vector of Bernoulli random variables z ( l ) . In practice, this is implemented as a ( l ) = m ( l ) θ ( o ( l − ) (8) o ( l ) = a ( l ) (cid:12) z ( l ) , z ( l ) ∼ B (1 − p ) (9) a ( l +1) = m ( l +1) θ ( o ( l ) ) , (10)5here l denotes the layer number after which a dropout layer is applied, B is a multivariate distribution withindepenent Bernoulli components, and p the selected dropout probability. The inclusion of dropout layers inducestochasticity to the vessel lumen segmentation process, resulting in random collections of points y obtained as y = m θ ( x ; z ) , z ∼ B (1 − p ) , (11)where z is a vector containing all Bernoulli dropout variables throughout the network. As discussed above, thesevariables are kept the same for every two-dimensional segmentation in a single cardiovascular model instance.For the network architecture, we build on our previous work [19] and use the GoogleNet architecture [38] appro-priately modified for vessel lumen regression, which consists of a CNN encoder followed by fully-connected layers totransform the encoded vector into the vessel lumen space. A dropout layer is applied to the output of the penulti-mate layer in the network in order to inject stochasticity (see Figure 4). A GoogleNet architecture was selected forcomputational efficiency, as the proposed algorithm is distributed as a SimVascular plug-in, targeting users withoutspecialized hardware. The GoogleNet architecture is computationally efficient due to the use of convolutional andpooling layers with different dimensions to compress the input image while still retaining necessary input information.For more details the reader is referred to [38]. While more recent networks have been developed, earlier studies weconducted [19] showed that the GoogleNet network achieved accuracy comparable to human experts on a 2D vessellumen segmentation task and so is sufficient for the purposes of this work.Figure 4: CNN based vessel lumen regression with and without dropout sampling. During training, the network’s weights are initialized according to the variance-scaling approach discussed in [9] andoptimized using stochastic gradient-descent ADAM algorithm [15]. We apply the angular distance transform to eachground-truth lumen surface p i , . . . , p Ki to transform it to a ground-truth vector y i and create a training dataset ofsize N d consisting of the following collection of image-radii pairs D = { ( x , y ) , ..., ( x N d , y N d ) } . (12)Additionally, we employ a l loss of the form l ( y i , (cid:98) y i ) = (cid:118)(cid:117)(cid:117)(cid:116) K (cid:88) j =1 ( y ji − (cid:98) y ji ) , L ( x , y ; θ , z , N b ) = 1 N b N b (cid:88) i =1 l ( y i , (cid:98) y i ) = 1 N b N b (cid:88) i =1 l [ y i , m θ ( x i ; z )] , (13)where y i and (cid:98) y i = m θ ( x i ; z ) represent the i -th ground-truth collection of normalized lumen radii and neural networkprediction, respectively, while L ( x , y ; θ , z , N b ) represents the expected loss over a given batch of training examples.We also pre-process each input gray-scale image x i ∈ R H × W , by computing a normalized image (cid:101) x i having zero meanand unit variance pixel intensities expressed as (cid:101) x jki = x jki − µ x σ x , j = 1 , . . . , H, k = 1 , . . . , W, µ x = 1 H H (cid:88) j =1 W (cid:88) k =1 x jki , σ x = H H (cid:88) j =1 W (cid:88) k =1 ( x jki − µ x ) / , (14)6here µ x and σ x are the mean and standard deviation pixel intensities, respectively. Finally, the training dataset D is augmented by randomly rotating and cropping each pair of image slice and 2D vessel lumen segmentation. Our dataset consists of 50 CT and 54 MR contrast-enhanced 3D medical image volumes, all publicly available from theVascular Model Repository (VMR) [45]. For each image volume, the VMR contains vessel pathlines, segmentations,3D patient-specific models and hemodynamic simulation results (see Figure 3) created in SimVascular by expert users,in many cases with supervision from a radiologist. To avoid anisotropic voxel spacing, all image volumes were re-sampled keeping an isotropic voxel spacing of 0.029 cm, which ensures the largest vessel diameter to be around 100pixels, a relatively small window size which reduces the network computation and memory requirements. Specifically,we used a window size H × W of 160 ×
160 pixels to allow the full range of vessel sizes to be represented with sufficientresolution by each two-dimensional slice. Finally, we split the data into training, validation and testing sets, of 86, 4and 14 volumes, respectively. This resulted in 16004, 239 and 6317 cross-sectional images and vessel lumen surfacepoint labels for the training, validation and testing sets, respectively.
The cardiovascular model generation process discussed above results in three-dimensional tetrahedral meshes whichprovide a domain Ω ⊂ R where the incompressible Navier-Stokes equations are solved. These are the equationsdescribing the evolution in time of a Newtonian fluid of constant density ρ , in a domain Ω with boundary ∂ Ω =Γ D ∪ Γ N , partitioned according to the application of Dirichlet and Neumann boundary conditions, respectively ρ ∂ u ∂t + ρ ( u · ∇ ) u − ∇ · τ = f in Ω × [0 , T ] ∇ · u = 0 in Ω × [0 , T ] u = g on Γ D × [0 , T ] τ · (cid:98) n = h on Γ N × [0 , T ] u (0) = u in Ω × { } , (15)in which u is fluid velocity, p is the fluid pressure, f is a given forcing term, (cid:98) n is the outward directed unit normalvector to Γ N , and τ the viscous stress tensor defined as σ = − p I + 2 µ (cid:15) ( u ). Let µ be the dynamic viscosity ofthe fluid, I the second order identity tensor and (cid:15) ( u ) the strain-rate tensor defined as (cid:15) ( u ) = ( ∇ u + ∇ u T ). Thefunctions g and h are given Dirichlet and Neumann boundary data, while u is the initial condition.We numerically solve the system (15) using a Streamline Upwind Petrov Galerkin (SUPG) finite element methodimplemented in the SimVascular flow solver (svSolver) [2], which contains specialized routines for cardiovascularCFD such as backflow stabilization [20], algebraic system solvers and preconditioners [34] and a large collection ofphysiologic boundary conditions (see, e.g., [20, 22, 21, 34]). The numerical solution is integrated in time using asecond-order generalized- α method [11]. We also apply RCR boundary conditions for generic outlets and a coronarylumped parameter boundary condition for coronary artery outlets, respectively (see, e.g., [22, 44, 14]). Finally, werestrict our attention to simulations with rigid walls. Our model generation procedure generates a set of discrete meshes Y h := { Y h , ..., Y hN y } (here the superscript h is used to indicate the size of the discrete mesh). Numerical solution of the Navier-Stokes equations on each meshsubsequently produces a set of velocity fields U h := { U h , ..., U hN y } and pressure fields P h := { P h , ..., P hN y } with U hi : Y i × [0 , T ] → R and P hi : Y i × [0 , T ] → R .Our objective is to calculate relevant Monte Carlo (i.e., sample) statistics using the ensembles U h and P h .However, the mesh geometry and hence the solution domain Ω is not constant for different realizations of the flowand pressure fields. This precludes us from considering quantities of interest defined at specific point locations inΩ, and we focus instead on output quantities that do hold meaning in the context of varying geometry. Consider ageneric model result r i ( ω , t ) for the i -th geometry realization, and the cross-sectional area A ji corresponding to the j -th slice location. We define the quantity q ji ( t ) = 1 | A ji | (cid:90) A ji r i ( ω , t ) d Γ , or the quantity q ji ( t ) = 1 | ∂A ji | (cid:90) ∂A i ( s ) r i ( ω , t ) d Γ , (16) ∂A ji (e.g., wall shear stress). The secondtype of quantities are time-average versions of the q ji ( t ), such as, q ji = 1( T − T ) (cid:90) T T q ji ( t ) dt. (17)In this study, we focus on quantities r i ( ω , t ) such as the pressure as well as the wall shear stress and velocitymagnitudes. Finally, for our Monte Carlo trials we choose to report a relative measure of variability, i.e., thecoefficient of variation, defined as CoV = σµ , (18)where σ and µ are the sample mean and the standard deviation of the quantities of interest q ji ( t ) or q ji , respectively.We further report confidence intervals for our Monte Carlo estimates. By the Central Limit Theorem the samplemean tends to a Gaussian random variable with standard deviation˜ σ = σ √ N , (19)where N is the sample size. In many cases the absolute value of the sample mean is small, therefore, for the sake ofclarity, we report the confidence interval as a percentage of the sample mean, that is CI = ± α ˜ σµ , (20)where α is a constant depending on the confidence level (e.g. α = 2 for 95% confidence). The CoV is thus an estimateof the variability of a particular QOI and CI is a measure of the accuracy of our sample mean estimates. To better understand the shape variation in our cardiovascular model samples we use the Principal ComponentAnalysis (PCA) algorithm to compute a low rank factorization of the matrix X ∈ R M × N constructed from thesamples. In particular we compute X − ¯ X = U K Σ K V TK , (21)where ¯ X is the mean of the model sample, Σ K is a diagonal matrix containing the singular values on its diagonaland the columns of U ∈ R M × K represent a reduced-order basis for the deviation of the cardiovascular models fromthe mean model. We use U to study the modes of variation in our cardiovascular model samples.By the properties of the PCA factorization, the columns of U K are ordered starting with the modes that cap-ture the most variance of X . The first column thus represents the most significant mode, the second column thesecond most significant mode etc. revealing the dominant modes in which segmentation uncertainty is causing thecardiovascular model geometry to vary. The first anatomy we consider consists of the bifurcation of the abdominal aorta in the two iliac arteries, presentedthrough a model having one inlet and two outlets (see Figure 5b). The inlet boundary condition is chosen to bea typical physiological waveform, corresponding to an average inflow of 6 L/min (see Figure 5c), while outflowRCR boundary conditions are applied, where the resistance and compliance parameters were preliminarily tuned toproduce a realistic outlet pressure range of 80-120 mmHg (see Table 5a). Initially, we conducted a mesh convergencestudy with meshes comprising 100 , ,
000 and 1 , ,
000 tetrahedral elements and boundary layer mesh with5 layers, and compared these to a reference mesh with 3 , ,
000 elements. Mesh convergence was assessed by firsttime-averaging the QOI and then using the mean absolute error (cid:15) h = 1 | Ω | (cid:90) Ω | ¯ f h ( ξ ) − ¯ f ∗ ( ξ ) | dΩ , (22)where ¯ f h and ¯ f ∗ are the time-averaged QOI on the investigated and reference mesh respectively. The 1 , , . .
7% error for the pressure, WSS magnitude and velocity magnituderespectively (see Figures 5d, 5e and 5f) and is employed in all additional experiments. The size of the Monte Carloensemble is finally selected equal to 150 for this anatomy.8essel R C R aorta 237.0 0.00115 2370left iliac 376.0 0.00085 3760 (a) RCR boundary condition parameters. Resis-tance in dyne · s/cm , capacitance in cm /dyne .(b) Patient-specific model surface. . . . . . . time ( s )050100150200250300 i n fl o w ( c m / s ) (c) Aortic inlet waveform . . . . . . . . | p − p ∗ | (d) Pressure convergence study. . . . . . . . . | T A W SS − T A W SS ∗ | (e) WSS convergence study. . . . . . . . . . || V | − | V ∗ || (f) Velocity convergence study. Figure 5: Aorto-iliac bifurcation model with boundary conditions (a,c), lumen surface (b) and mesh convergenceanalysis (d,e,f).
The second anatomy considered in this study includes the aorta and its main branches from an abdominal CT imageof a patient with an abdominal aortic aneurysm (AAA, see Figure 6a), subject to the same aortic inflow used for theprevious anatomy (see Figure 5c). RCR boundary condition parameters are reported in Table 6b. For this anatomy,a family of 110 geometries was generated through Monte Carlo sampling. The sample size was reduced as 110 modelswas sufficient for statistical convergence. We conducted a mesh convergence study, using the mean absolute error(22) and with meshes having roughly 500 , ,
000 and 3 , ,
000 elements and boundary layer mesh with 5layers and compared these to a reference mesh with 7 , ,
000 elements. The 3 , ,
000 element mesh showed aless than 0 . The third model we consider includes the left anterior descending (LAD) and left circumflex (LCx) coronary arteriesextracted from a CT image volume, also studied in [34]. Coronary lumped parameter boundary condition valueswere selected to produce physiological pressure ranges (Fig. 7a). The coronary simulations use a pulsatile coronaryinflow waveform (Fig. 7c). A sample size of 110 models was used for the Monte Carlo trials. We conducted amesh convergence study, using the mean absolute error (22) and with meshes with roughly 500 , , ,
000 and1 , ,
000 elements and boundary layers with 5 layers and compared these to a reference mesh with 3 , , , ,
000 element mesh showed a less than 0 . a) Patient-specific AAA model. Vessel R C R Aorta 660 0.0004 6600Celiac Hepatic 1600 0.00017 16000Celiac Splenic 910 0.0003 9100Ext Iliac Left 1155 0.00024 11550Renal Left 691 0.00039 6910Renal Right 1220 0.00022 12200SMA 1040 0.00026 10400 (b) RCR boundary conditions for AAA model,resistance in dyne · s/cm , capacitance in cm /dyne . .
14 0 .
16 0 .
18 0 . . . . . . . | p − p ∗ | (c) Pressure convergence study. .
14 0 .
16 0 .
18 0 . . . . . . | T A W SS − T A W SS ∗ | (d) WSS convergence study. .
14 0 .
16 0 .
18 0 . . . . . . . . || V | − | V ∗ || (e) Velocity convergence study. Figure 6: AAA model anatomy (a) with boundary conditions (b) and mesh convergence analysis (c,d,e).
Before commenting on the model results, we first investigate the statistical properties of the segmentations producedby our dropout network. In addition, we compare network samples and lumen segmentations produced by a numberof expert SimVascular users. To do so, we selected four representative image volumes with vessel centerlines fromthe Vascular Model Repository [45], including a cerebrovascular anatomy imaged by MR, a coronary anatomy withaneurysms caused by Kawasaki disease imaged by CT, a coronary anatomy following bypass graft surgery imaged byCT, and an pulmonary anatomy imaged by MR. Slices were selected at discrete intervals along the vessel centerlines,and each location was segmented by three individual SimVascular experts (see Figure 3), resulting in a total of 290segmentations per expert. For the same slices, 50 neural network lumen samples were generated for various dropoutprobabilities, i.e., p = 0 . p = 0 . p = 0 .
4, and used for statistical analysis.Segmentation radius CoV observed for SimVascular expert users is separated into two classes, i.e., large vessels(r > ≤ p = 0 . p = 0 . essel R C R C R LAD . · − . · − LAD − D . · − . · − LCx . · − . · − LCx − OM . · − . · − LCx − OM . · − . · − LCx − OM . · − . · − (a) RCRCR boundary condition parameters for leftcoronary artery model, resistance in dyne · s/cm , ca-pacitance in cm /dyne .(b) Left coronary artery model anatomy. . . . . . . time ( s )12345 i n fl o w ( c m / s ) (c) Coronary inlet waveform .
035 0 .
040 0 .
045 0 .
050 0 . . . . . . . . | p − p ∗ | (d) Pressure convergence study. .
035 0 .
040 0 .
045 0 .
050 0 . . . . . . | T A W SS − T A W SS ∗ | (e) WSS convergence study. .
035 0 .
040 0 .
045 0 .
050 0 . . . . . . . . . || V | − | V ∗ || (f) Velocity convergence study. Figure 7: Left coronary artery model with boundary conditions (a,c), lumen surface (b) and mesh convergenceanalysis (d,e,f).For large vessels, the network with dropout p = 0 . p = 0 . p = 0 . p = 0 . p = 0 . p the variability in the segmentations is reduced,and they all converge closely to the mean lumen profile (see Figures 10b-10d and 10f-10h). For the experiments inthe remainder of this work a dropout probability of 0 . . . . . . . . . . . . . r a d i u s C o V SV experts (a) SimVascular experts . . . . . . . . . . . . r a d i u s C o V dropout p = 0 . (b) CNN dropout p = 0 . . . . . . . . . . . . . r a d i u s C o V dropout p = 0 . (c) CNN dropout p = 0 . . . . . . . . . . . . . r a d i u s C o V dropout p = 0 . (d) CNN dropout p = 0 . Figure 8: Radius CoV against radius for SimVascular experts and GoogleNet network with different levels of dropout.11 . . . . . r a d i u s C o V SV expertsdropout p = 0 . p = 0 . p = 0 . (a) Large vessels r > . cm . . . . . . . r a d i u s C o V SV expertsdropout p = 0 . p = 0 . p = 0 . (b) Small vessels r < . cm Figure 9: Comparison of mean radius CoV for large and small vessels between SimVascular experts and the proposednetwork, with varying dropout probability. The distributions of CoV shown above were obtained by collecting oneCoV for each cross-section x i . (a) SimVascular experts (b) CNN dropout p = 0 . p = 0 . p = 0 . p = 0 . p = 0 . p = 0 . Figure 10: Vessel lumen segmentation generated for large (top) and small (bottom) vessels by expert users and theproposed network with varying dropout probability.
The lumen generated from our dropout network show good qualitative agreement with the depicted vessel lumen forthe left iliac artery (Figure 11a) and is able to correctly identify the relevant main branch even in the presence ofsurrounding tissue noise and branching vessels.Variation in the segmentation radii appears to be limited, with standard deviation σ r between 0.005 cm and 0.01cm (see Figure 11b and Figure 11c). In addition, variability in σ r appears to increase with decreasing vessel size,likely due to the typically poorer resolution of smaller vessels. The roughly constant σ r also results in a CoV thatincreases with decreasing vessel size, as typically seen for segmentations performed by expert operators [19], albeitwith larger magnitude.Dominant modes from PCA appear to perturb whole vessels, such as the right iliac or large segments of the aorta(Figures 12a-12e) and to increase with the reduction in the vessel radius. Conversely, higher modes involve local12eometrical changes such as the bulbous region at the proximal end of the left iliac. These latter modes, however,contribute less, being associated with much smaller singular values (e.g. the 19th and 20th mode have singular valuesthat are a factor of approximately 5.7 smaller than the first mode). Despite generating segmentations independently,this confirms that the geometric variability produced by the proposed dropout network is distributed across theentire model, and increases with small vessel radii.Convergence of Monte Carlo statistics such as the mean, standard deviation and CoV for the pressure, TAWSSand velocity magnitude integrated over the aorta appears to be satisfactory (Figure 13), with 95% confidence intervalsfor the mean showing faster convergence for the pressure, then velocity magnitude and finally TAWSS.The larger sensitivity of TAWSS to geometric uncertainty can be observed from their time histories over the lasttwo cardiac cycles, as shown in Figures 14a, 14b. Pressure variability tends to decrease towards the distal end of thevessel, whereas TAWSS and velocity magnitude show an opposite trend. This relates to an increase in the wall shearstress and velocity after the bifurcation which, for this model, amplifies the effect of the geometric uncertainty forthese two QoIs. Pressure uncertainty instead depends on the variability of the vessel resistance which cumulates thecontributions of each uncertain segmentation along the vessel.CoVs were found to be approximately 0.4%, 1.5% and 3% for pressure, velocity magnitude and TAWSS, respec-tively. Thus, CoVs for TAWSS and velocity magnitude are roughly a factor of 10 and 5 larger compared to thepressure CoV, highlighting the increased sensitivity of TAWSS and velocity to geometry variation (see Table 1). (a) Aorta and Iliacs case - Left Iliac . . . . . . . . . . . . r a d i u s C o V ( % ) (b) Radius CoV .
50 0 .
75 1 .
00 1 . . . . . . . . . . σ r ( c m ) (c) Radius standard deviation Figure 11: Lumen segmentation samples and radius CoV/standard deviation for aorto-iliac bifurcation test case,computed over cross-sectional slices x i , i = 1 , . . . , n indicates the number of cross-sectional slices for the associated vessel. Path Aorta ( n = 95) left iliac ( n = 61) Path Aorta left iliacRadius mean [cm] 0.84 0.61
TAWSS mean [dyne/cm ] 40.24 46.55 Radius CoV
TAWSS CoV
Radius conf.
TAWSS conf.
Pressure mean [mmHg] 98.43 96.02
Velocity mean [cm/s] 42.84 41.78
Pressure CoV
Velocity CoV
Pressure conf.
Velocity conf. a) λ = 16 .
86 (b) λ = 13 .
80 (c) λ = 3 .
00 (d) λ = 2 . . . . . . . . . S i n g u l a r V a l u e (e) Singular values Figure 12: PCA modes overlayed on mean aorto-iliac bifurcation model geometry. . . . . . . p r e ss u r e ( mm H g ) m e a n . . . . . . T A W SS m ag n i t ud e ( d y n e / c m ) m e a n . . . . . V e l o c i t y m ag n i t ud e ( c m / s ) m e a n . . . . . p r e ss u r e ( mm H g ) C o V . . . T A W SS m ag n i t ud e ( d y n e / c m ) C o V . . . . . V e l o c i t y m ag n i t ud e ( c m / s ) C o V Figure 13: Monte Carlo moment traces for aorto-iliac bifurcation model QoIs. p r e ss u r e ( mm H g ) W SS m ag n i t ud e ( d y n e / c m ) V e l o c i t y m ag n i t ud e ( c m / s ) (a) Aorta p r e ss u r e ( mm H g ) W SS m ag n i t ud e ( d y n e / c m ) V e l o c i t y m ag n i t ud e ( c m / s ) (b) Left Iliac Figure 14: Outlet QoIs and ± σ interval for aorta-iliac bifurcation model. The lumen generated by the network agrees well with the vessel lumen images, even in the presence of noise andbifurcations (see Figure 17a). Values of σ r from the dropout network are generally in the range of 0.005-0.01 cm, withsome outliers, particularly for larger vessel lumens. Segmentation of small vessels is affected by increased uncertainty,with radius CoV of 3%, versus 1% for the largest vessels. PCA shape analysis shows the first two modes affectinglarge scale model features such as the aortic aneurysm as well as the entire aorta and iliac branches (see Figure 18a).14 . . . . . r a d i u s ( c m ) . . . p r e ss u r e ( mm H g ) . . . T A W SS m ag n i t ud e ( d y n e / c m ) . . . V e l o c i t y m ag n i t ud e ( c m / s ) (a) Aorta . . . . . r a d i u s ( c m ) . . . p r e ss u r e ( mm H g ) . . . T A W SS m ag n i t ud e ( d y n e / c m ) . . . V e l o c i t y m ag n i t ud e ( c m / s ) (b) left Iliac Figure 15: Time averaged QoIs and ± σ interval for aorto-iliac bifurcation model, plotted along the vessel centerline. (a) (b) (c) (d) Figure 16: Nearest neighbor interpolation of cross-sectional time-averaged CoVs for aorto-iliac bifurcation model.The 19 th and 20 th modes (with singular values a factor of four smaller than the first mode) are instead associatedwith local features, like celiac branches or by asymmetric aneurysm perturbations (see Figure 18c).The relative confidence intervals of the Monte Carlo estimates for the mean pressure, TAWSS and velocitymagnitude were found to be within the ranges 0.07-0.8%, 0.44-0.79% and 0.28-0.65%, indicating a satisfactoryconvergence, particularly compared to the observed CoV for the same vessels (see Table 2 and Figure 19). Outletprofiles show increased variability in branch vessels for all QoIs with respect to the Aorta (see Figures 20a-20g).Time average flow results show increasing TAWSS and velocity magnitude variability towards the distal end of thevessel (see Figures 21a-21g).As expected, radius CoV is inversely proportional to the vessel size (see Figure 22a). Relative pressure variabilityis approximately uniform along the path of each vessel, and particularly elevated in the celiac hepatic branch, due tothe fact that it branches off of the celiac splenic which itself branches off of the aorta (see Figure 22b). The TAWSSCoV appears to be significant for small branches (see Figure 22c), and is equal to 0.06 in the aneurysm region and0.02 in the proximal regions of the aorta. Similarly, the velocity magnitude CoVs is found to be 0.04 within theaneurysm and 0.02 in the upstream aorta. This illustrates how diseased regions like aneurysms can lead to increasedgeometric uncertainty, most likely due to increased ambiguity of the vessel lumen shape and increased surroundingnoise sources in the input image volume, leading to higher neural network output variability. We note again thatthis is directly learned from the image volume by the neural network. Even for the smallest coronary arteries and in the presence of significant surrounding heart tissue, the vessel lumenshape was qualitatively well captured by the proposed dropout network (see Figure 23a). A PCA quantification ofsegmentation uncertainty shows dominant modes distributed along major arteries, and higher modes inducing localchanges to smaller branches (see Figure 24a). While the relative variance is higher in the small vessel branches whencompared to the larger branches, PCA determines modes based on absolute variance. The absolute variance is largerin the large vessels and explains their presence in the dominant PCA modes.A satisfactory convergence is observed for the Monte Carlo statistical moments after 110 model evaluations, with15 a) Abdominal Aortic Aneurysm Case - CeliacSplenic . . . . . . . . . . r a d i u s C o V ( % ) (b) Radius CoV . . . . . . . . . . . . σ r ( c m ) (c) Radius standard deviation Figure 17: Lumen segmentation samples and radius CoV/standard deviation for abdominal aortic aneurysm testcase, computed over cross-sectional slices x i , i = 1 , . . . , (a) λ = 24 .
60 (b) λ = 21 .
41 (c) λ = 6 .
28 (d) λ = 5 . S i n g u l a r V a l u e (e) Singular Values Figure 18: PCA modes overlayed on mean abdominal aortic aneurysm model geometry.asymptotic traces for more than 50 samples (see Figure 25). The 95% relative confidence intervals for the MonteCarlo estimates of the mean are approximately equal to 0.4%, 3% and 1.5% for pressure, TAWSS and velocitymagnitude, respectively, and are well below the CoV computed for the same QoIs (see Tables 3).The outlet time histories reflect the diastolic nature of the coronary flow (see Figures 26a-26f), where the WSSexhibits the largest uncertainty followed by velocity magnitude and pressure. Unlike the other two anatomies con-sidered in the previous sections, geometrical uncertainty significantly affects hemodynamic model outputs due tosmaller vessel sizes, as previously suggested in the literature in the context of coronary artery disease [28]. Timeaveraged quantities over the vessel length show a similar pattern (see Figure 27a-27f). In particular, TAWSS andvelocity uncertainty appear to be very similar and correlated with the vessel radius. Specifically, smaller radii (withlarger radius variability) produce larger TAWSS and velocity uncertainty. A different behavior is instead observed forthe
LCx − OM branch, where large TAWSS and velocity uncertainty are associated with a larger radius. However,this phenomenon is localized at the proximal end of the vessel and probably triggered by the bifurcation nearby.Pressure, TAWSS and velocity magnitude CoVs were approximately equal to 2%, 10-20% and 6-15%, respectively16
50 100sample size94 . . . p r e ss u r e ( mm H g ) m e a n . . . . . T A W SS m ag n i t ud e ( d y n e / c m ) m e a n . . . . . V e l o c i t y m ag n i t ud e ( c m / s ) m e a n . . . . p r e ss u r e ( mm H g ) C o V . . . . T A W SS m ag n i t ud e ( d y n e / c m ) C o V . . . . V e l o c i t y m ag n i t ud e ( c m / s ) C o V Figure 19: Monte carlo convergence - Abdominal Aortic Aneurysm caseTable 2: Monte Carlo sample mean, coefficient of variation (CoV) and 95% relative confidence interval for all QoIsin abdominal aortic aneurysm model. n indicates the number of cross-sectional slices for the associated vessel. Path Aorta Celiac hepatic Celiac splenic Ext. iliac left Renal left Renal right SMA ( n = 160) ( n = 30) ( n = 69) ( n = 169) ( n = 51) ( n = 35) ( n = 61) Radius mean [cm] 0.63 0.26 0.36 0.38 0.32 0.29 0.38
Radius CoV
Radius conf.
Pressure mean [mmHg] 96.45 89.11 94.51 94.35 95.29 92.35 99.26
Pressure CoV
Pressure conf.
TAWSS mean [dyne/cm ] 47.66 125.70 85.11 58.47 87.28 117.50 35.35 TAWSS CoV
TAWSS conf.
Velocity mean [cm/s] 38.05 79.84 52.21 43.88 51.43 66.67 31.98
Velocity CoV
Velocity conf. (see Table 3). Notably, the
LAD − D , LCx − OM , and LCx − OM branches showed larger TAWSS and velocityCoVs, equal to 16%, 24.5% and 92.9%, still explained by the small range of vessel sizes present in the coronaryanatomy, which amplifies the segmentation uncertainty, particularly towards the distal end of each vessel (see Table 3and Figure 28a). The pressure CoV appears to be larger in the LCx and its branches
LCx − OM and particularly LCx − OM . This is explained by the relatively small radius of such vessels, which increases resistance and amplifiesgeometric uncertainty and by the variations in flow split caused by the relatively smaller LCx branching off the
LAD trunk (see Figure 28b). Finally, the largest WSS and velocity CoVs are located near bifurcations, due to higher localflow variability and more ambiguity in the definition of the vessel lumen (see Figure 28c and Figure 28d).
Our experiments in the previous sections illustrate how the proposed Bayesian dropout network generates vessellumen segmentations characterized by a σ r between 0.005 cm and 0.01 cm. This translates in radius CoVs in therange 1-5%, where smaller vessel sizes are associated with larger radius variability, as expected. Furthermore, ourcardiovascular model sampling process resulted in realistic whole model variation as observed from the PCA modesacross model realizations.Amongst the output QoIs considered in this study, wall shear stress was the most impacted by geometry uncer-tainty, followed by velocity magnitude and then pressure. Wall shear stress CoV ranged from 3% to 20%, whereas17 p r e ss u r e ( mm H g ) W SS m ag n i t ud e ( d y n e / c m ) V e l o c i t y m ag n i t ud e ( c m / s ) (a) Aorta p r e ss u r e ( mm H g ) W SS m ag n i t ud e ( d y n e / c m ) V e l o c i t y m ag n i t ud e ( c m / s ) (b) Celiac Hepatic p r e ss u r e ( mm H g ) W SS m ag n i t ud e ( d y n e / c m ) V e l o c i t y m ag n i t ud e ( c m / s ) (c) Celiac Splenic p r e ss u r e ( mm H g ) W SS m ag n i t ud e ( d y n e / c m ) V e l o c i t y m ag n i t ud e ( c m / s ) (d) Left external Iliac p r e ss u r e ( mm H g ) W SS m ag n i t ud e ( d y n e / c m ) V e l o c i t y m ag n i t ud e ( c m / s ) (e) Left Renal p r e ss u r e ( mm H g ) W SS m ag n i t ud e ( d y n e / c m ) V e l o c i t y m ag n i t ud e ( c m / s ) (f) Right Renal p r e ss u r e ( mm H g ) W SS m ag n i t ud e ( d y n e / c m ) V e l o c i t y m ag n i t ud e ( c m / s ) (g) SMA Figure 20: Outlet QoIs and ± σ interval for abdominal aortic aneurysm model.Table 3: Monte Carlo sample mean, coefficient of variation (CoV) and 95% relative confidence interval for all QoIsin left coronary artery model. n indicates the number of cross-sectional slices for the associated vessel. Path LCx LCx-OM LCx-OM LCx-OM LAD LAD-D ( n = 29 ( n = 58) ( n = 20) ( n = 12) ( n = 48) ( n = 55) Radius mean [cm] 0.14 0.12 0.09 0.10 0.13 0.11
Radius CoV
Radius conf.
Pressure mean [mmHg] 93.46 92.37 91.78 91.92 96.02 95.14
Pressure CoV
Pressure conf.
TAWSS mean [dyne/cm ] 39.27 48.43 26.37 11.91 12.78 26.17 TAWSS CoV
TAWSS conf.
Velocity mean [cm/s] 21.91 28.04 13.69 6.66 8.56 15.50
Velocity CoV
Velocity conf. velocity magnitude and pressure resulted in CoV ranges equal to 1.4-15% and 0.2-3%, respectively. Larger variabilityis observed in the left coronary artery model, due to the prevalence of vessels with small radius.18 . . . . . . r a d i u s ( c m ) . . . p r e ss u r e ( mm H g ) . . . T A W SS m ag n i t ud e ( d y n e / c m ) . . . V e l o c i t y m ag n i t ud e ( c m / s ) (a) Aorta . . . . . . r a d i u s ( c m ) . . . p r e ss u r e ( mm H g ) . . . T A W SS m ag n i t ud e ( d y n e / c m ) . . . V e l o c i t y m ag n i t ud e ( c m / s ) (b) Celiac Hepatic . . . s/L . . r a d i u s ( c m ) . . . s/L p r e ss u r e ( mm H g ) . . . s/L T A W SS m a g n i t ud e ( d y n e / c m ) . . . s/L V e l o c i t y m a g n i t ud e ( c m / s ) (c) Celiac Splenic . . . s/L . . r a d i u s ( c m ) . . . s/L p r e ss u r e ( mm H g ) . . . s/L T A W SS m a g n i t ud e ( d y n e / c m ) . . . s/L V e l o c i t y m a g n i t ud e ( c m / s ) (d) Left external Iliac . . . s/L . . r a d i u s ( c m ) . . . s/L p r e ss u r e ( mm H g ) . . . s/L T A W SS m a g n i t ud e ( d y n e / c m ) . . . s/L V e l o c i t y m a g n i t ud e ( c m / s ) (e) Left Renal . . . s/L . . r a d i u s ( c m ) . . . s/L p r e ss u r e ( mm H g ) . . . s/L T A W SS m a g n i t ud e ( d y n e / c m ) . . . s/L V e l o c i t y m a g n i t ud e ( c m / s ) (f) Right Renal . . . s/L . . r a d i u s ( c m ) . . . s/L p r e ss u r e ( mm H g ) . . . s/L T A W SS m a g n i t ud e ( d y n e / c m ) . . . s/L V e l o c i t y m a g n i t ud e ( c m / s ) (g) SMA Figure 21: Time averaged QoIs and ± σ interval for abdominal aortic aneurysm model, plotted along the vesselcenterline.To quantify the relative importance of geometric uncertainty with respect to other sources relevant in hemody-namic simulations, we compare the output variability found in our study to those found in studies investigating othersources of uncertainty. Uncertainty due to coronary pressure waveform, intramyocardial pressure, morphometry expo-nent and vascular wall Young’s modulus was recently investigated in [34], in the context of coronary artery modeling.Coronary pressure waveform uncertainty resulted in a 7% CoV for the average pressure and <
7% for TAWSS andvelocity magnitude. Intramyocardial pressure uncertainty produced CoV of roughly 25% for TAWSS and velocitymagnitude. Morphometry exponent uncertainty resulted in a 2% CoV for velocity magnitude and negligible impacton TAWSS magnitude and pressure. Finally, vascular wall Young’s modulus uncertainty had negligible impact onthe hemodynamics, as suggested in the literature [40]. The CoV of 1.4-15% and 3-20% for velocity magintude andTAWSS produced by geometric uncertainty in this work is thus comparable to the CoV due to coronary pressureand intramyocardial pressure uncertainty, but larger than the CoV due to morphometry exponent and vascular wallYoung’s modulus. The pressure CoV of 0.2-3% due to geometry uncertainty was smaller than that produced byintramyocardial pressure uncertainty.A multifidelity uncertainty quantification approach was used in [6] to investigate simulation uncertainty due tomaterial property and boundary condition uncertainty, for both healthy and diseased aortic and coronary anatomies.The uncertain parameters were uniformly distributed with ±
30% variation around their means. Resulting CoV were ∼
3% for pressure, <
1% for velocity magnitude and ∼ a) (b) (c) (d) Figure 22: Nearest neighbor interpolation of cross-sectional time-averaged CoVs for abdominal aortic aneurysmmodel. (a) Left Coronary Artery case -
LAD − D .
10 0 .
15 0 . r a d i u s C o V ( % ) (b) Radius CoV .
10 0 .
15 0 . . . . . . . . . . σ r ( c m ) (c) Radius standard deviation Figure 23: Lumen segmentation samples and radius CoV/standard deviation for left coronary artery test case,computed over cross-sectional slices x i , i = 1 , . . . , a) λ = 8 .
09 (b) λ = 6 .
72 (c) λ = 1 .
91 (d) λ = 1 . S i n g u l a r V a l u e (e) Singular Values Figure 24: PCA modes overlayed on mean left coronary artery model geometry. . . . . p r e ss u r e ( mm H g ) m e a n . . . . . T A W SS m ag n i t ud e ( d y n e / c m ) m e a n . . . . . V e l o c i t y m ag n i t ud e ( c m / s ) m e a n . . . . . p r e ss u r e ( mm H g ) C o V . . . . . T A W SS m ag n i t ud e ( d y n e / c m ) C o V . . . . . . V e l o c i t y m ag n i t ud e ( c m / s ) C o V Figure 25: Monte Carlo moment traces for left coronary artery model QoIs.to the other two anatomies. This suggests that geometric uncertainty might play a dominant role for anatomiescharacterized by small vessel sizes (the coronary circulation is a particularly relevant case), especially for stenoticlesions, associated with substantial radius uncertainty.Table 4: Comparison between the hemodynamic effect geometry uncertainty and other sources of uncertainty.
Model type CoronaryUncertainty source
Pressure waveform Intramyocardial pressure Morphometry Exponent Wall Young’s modulus
Input distribution
Uniform Uniform Uniform Uniform
Reference [34] [34] [34] [34]
Pressure CoV
7% negligible negligible negligible
TAWSS CoV <
7% 25% negligible negligible
Velocity CoV <
7% 25% 2% negligible
Model Type Aorta/Coronary Coronary bypass graft Aorta/AAA/CoronaryUncertainty source
Boundary conditions Multiple clinical targets Full model geometry
Input distribution
Uniform ±
30% Assimilated from clinical data Dropout sampling
Reference [6] [39]
This paper
Pressure CoV ∼
3% - 0.2-3%
TAWSS CoV ∼ Velocity CoV <
1% - 1.4-15%
Even though this is the first systematic study in the literature combining machine learning and high-fidelity21 . . . p r e ss u r e ( mm H g ) . . . W SS m ag n i t ud e ( d y n e / c m ) . . . V e l o c i t y m ag n i t ud e ( c m / s ) (a) LCx . . . p r e ss u r e ( mm H g ) . . . W SS m ag n i t ud e ( d y n e / c m ) . . . V e l o c i t y m ag n i t ud e ( c m / s ) (b) LCx − OM . . . p r e ss u r e ( mm H g ) . . . W SS m ag n i t ud e ( d y n e / c m ) . . . V e l o c i t y m ag n i t ud e ( c m / s ) (c) LCx − OM . . . p r e ss u r e ( mm H g ) . . . W SS m ag n i t ud e ( d y n e / c m ) . . . V e l o c i t y m ag n i t ud e ( c m / s ) (d) LCx − OM . . . p r e ss u r e ( mm H g ) . . . W SS m ag n i t ud e ( d y n e / c m ) . . . V e l o c i t y m ag n i t ud e ( c m / s ) (e) LAD . . . p r e ss u r e ( mm H g ) . . . W SS m ag n i t ud e ( d y n e / c m ) . . . V e l o c i t y m ag n i t ud e ( c m / s ) (f) LAD − D Figure 26: Outlet QoIs and ± σ interval for left coronary artery model. . . . . . r a d i u s ( c m ) . . . p r e ss u r e ( mm H g ) . . . T A W SS m ag n i t ud e ( d y n e / c m ) . . . V e l o c i t y m ag n i t ud e ( c m / s ) (a) LCx . . r a d i u s ( c m ) p r e ss u r e ( mm H g ) T A W SS m ag n i t ud e ( d y n e / c m ) V e l o c i t y m ag n i t ud e ( c m / s ) (b) LCx − OM . . . r a d i u s ( c m ) p r e ss u r e ( mm H g ) T A W SS m ag n i t ud e ( d y n e / c m ) V e l o c i t y m ag n i t ud e ( c m / s ) (c) LCx − OM . . . r a d i u s ( c m ) . . . p r e ss u r e ( mm H g ) T A W SS m ag n i t ud e ( d y n e / c m ) V e l o c i t y m ag n i t ud e ( c m / s ) (d) LCx − OM . . r a d i u s ( c m ) p r e ss u r e ( mm H g ) T A W SS m ag n i t ud e ( d y n e / c m ) V e l o c i t y m ag n i t ud e ( c m / s ) (e) LAD . . r a d i u s ( c m ) p r e ss u r e ( mm H g ) T A W SS m ag n i t ud e ( d y n e / c m ) V e l o c i t y m ag n i t ud e ( c m / s ) (f) LAD − D Figure 27: Time averaged QoIs and ± σ interval for left coronary artery model, plotted along the vessel centerline.cardioascular models to study the effect of geometric uncertainty, we recognize several limitations. First, uncertaintypropagation is performed using standard Monte Carlo sampling. Even though a number of approaches in theliterature have shown promise to accelerate convergence to the true statistical moments (such as, e.g., stochasticcollocation [29], or generalized multiresolution expansions [30] in the context of cardiovascular flow), the Bernoullirandom vectors used in the dropout layer in this study have a dimension of around 10,000, which is extremelychallenging for approaches based on stochastic spectral expansion. Second, our study includes only one diseased22 a) (b) (c) (d) Figure 28: Nearest neighbor interpolation of cross-sectional time-averaged CoVs for left coronary artery model.anatomy, and showed that geometrical uncertainty is amplified in the aneurysm region of the abdominal aorta. Wealso found such amplification in a healthy left coronary artery model due to the typically smaller vessel radii. Thissuggests how cases of stenosed or calcified coronary arteries and vascular lesions, that are not often seen duringnetwork training, may offer new insights on the role of geometric uncertainty. Third, we assume that dropoutnetworks are able to learn output uncertainty from their training data [7], while, in practice, they provide only anapproximate representation of the true distribution of vessel lumen for a given image. Additionally, our segmentationsare generated at discrete slices along the centerline path, and thus the training data might act as a filter on the wholegeometric variability, leading the proposed algorithm to underestimate the true underlying geometric uncertainty.Removing these limitations would require new three-dimensional vessel segmentation paradigms, which is an activearea of research. Fourth, the path planning approach intrinsically limits the uncertainty at the bifurcations, as itrequires users to adjust pathlines so they originate within the parent vessel. This introduces user bias into thecardiovascular model samples and may constrain the underlying geometric uncertainty. Future work will be devotedto produce improved estimates for bifurcation uncertainty.
We have developed a Bayesian dropout network to generate families of two-dimensional lumen segmentations fromslices of a clinically acquired image volume. Of particular note is the fact that our neural network learns lumensegmentation uncertainty directly from the image training data and is thus able to generate lumen samples witha realistic uncertainty distribution. This was combined with vessel centerlines and a path-planning model buildingworkflow to create realizations of high-fidelity cardiovascular models with uncertain lumen surface. Finally, wecharacterized simulation output variability due to geometric uncertainty using Monte Carlo sampling.We have also analyzed the principal components of the lumen surfaces we generated, showing how, despite seg-menting slices independently and analyzing geometries characterized by a wide range of vessel radii, dominant modesappear to be equally distributed on the entire model, without amplifying any particular local feature. Addition-ally, our network generated vessel lumens with relatively constant radius standard deviation that were found to beindependent of vessel size. This resulted in increasing relative uncertainty for smaller vessels, similar to manualsegmentations generated by expert users.Experiments on an aortic bifurcation model, an abdominal aortic aneurysm model and a left coronary arterymodel showed that geometry uncertainty primarily resulted in wall shear stress and velocity magnitude uncertainty.This was true in particular for the coronary anatomy, characterized by smaller vessel sizes. Moreover, while TAWSSand velocity magnitude were impacted by geometrical uncertainty, especially near the distal ends of small vesselsand near bifurcations, pressure was only marginally affected. Compared to other sources of uncertainty, for example,the boundary conditions or material properties, the relative importance of geometry uncertainty was found to bedetermined by the particular patient-specific geometry being investigated and the vessel radius.Our method still requires one to manually create the vessel centerlines, which may be laborious, time consumingand introduce additional uncertainty. Automated methods to predict vessel centerlines and corresponding uncertaintycould therefore be combined with the proposed dropout network to further improve performance and increase modelbuilding efficiency. Finally, we have explored geometry uncertainty independently from other sources, disregarding23heir interaction. Future work will be devoted to combine all three sources of simulation uncertainty, i.e., boundaryconditions, material properties and geometry. acknowledgements
This research was funded by the National Science Foundation (NSF) Software Infrastructure for Sustained Innova-tion (SSI) grants 1663671 and 1562450, Computational and Data-Enabled Science and Engineering (CDSE) grant1508794, NSF CAREER award 1942662, National Institute of Health grants R01HL123689 and R01EB018302 andthe American Heart Association Precision Medicine Platform.
References [1] J. Biehler, M. W. Gee, and W. A. Wall. Towards efficient uncertainty quantification in complex and large-scalebiomechanical problems based on a Bayesian multi-fidelity scheme.
Biomechanics and Modeling in Mechanobi-ology , 14, 2014.[2] A. N. Brooks and T. J. R. Hughes. Streamline upwind/Petrov-Galerkin formulations for convection dominatedflows with particular emphasis on the incompressible Navier-Stokes equations.
Computer Methods in AppliedMechanics and Engineering , 1982.[3] M. J. Colebank, L. M. Paun, M. U. Qureshi, M. Chesler, D. Husmeier, M. S. Olufsen, and L. E. Fix. Influenceof image segmentation on one-dimensional fluid dynamics predictions in the mouse pulmonary arteries.
RoyalSociety Interface , 2019.[4] M. J. Colebank, M. U. Qureshi, and M. S. Olufsen. Sensitivity analysis and uncertainty quantification of 1-dmodels of pulmonary hemodynamics in mice under control and hypertensive conditions.
International Journalfor Numerical Methods in Biomedical Engineering , 2019.[5] V. G. Eck, W. P. Donders, J. Sturdy, J. Feinberg, T. Delhaas, L. R. Hellevik, and W. Huberts. A guide touncertainty quantification and sensitivity analysis for cardiovascular applications.
International Journal forNumerical Methods in Biomedical Engineering , 32, 2015.[6] C.M. Fleeter, G. Geraci, D.E. Schiavazzi, A.M. Kahn, and A.L. Marsden. Multilevel and multifidelity uncertaintyquantification for cardiovascular hemodynamics.
Computer Methods in Applied Mechanics and Engineering ,365:113030, 2020.[7] Y. Gal and Z. Ghahramani. Dropout as a Bayesian approximation: representing model uncertainty in deeplearning.
International Conference on Machine Learning , 2016.[8] A. M. Gambaruto, J. Janela, A. Moura, and A. Sequeira. Sensitivity of hemodynamics in a patient specificcerebral aneurysm to vascular geometry and blood rehology.
Mathematical Biosciences and Engineering , 2011.[9] K. He, X. Zhang, S. Ren, and J. Sun. Delving deep into rectifiers: surpassing human-level performance onImageNet classification.
International Conference on Computer Vision , 2015.[10] K Hornik. Approximation capabilities of multilayer feedforward networks.
Neural Networks , 4, 1991.[11] K. E. Jansen, C. H. Whiting, and G. M. Hulbert. A generalized- α method for integrating the filtered Navier-Stokes equations with a stabilized finite element method. Computer Methods in Applied Mechanics and Engi-neering , 190, 2000.[12] R. Jena and S. P. Awate. A Bayesian neural net to segment images with uncertainty estimates and goodcalibration.
International Conference on Information Processing in Imaging , 2019.[13] L. Joskowicz, D. Cohen, N. Caplan, and J. Sosna. Automatic segmentation variability estimation with segmen-tation priors.
Medical Image Analysis , 50, 2018.[14] H. J. Kim, I. E. Vignon-Clementel, J. S. Coogan, C. A. Figueroa, K. E. Jansen, and C. A. Taylor. Patient-specific modeling of blood flow and pressure in human coronary arteries.
Annals of Biomedical Engineering , 38,2010. 2415] D. Kingma and J. Ba. Adam: A method for stochastic optimization.
International Conference on LearningRepresentations , 2015.[16] M. Le, J. Unkelbach, N. Ayache, and H. Delingette. Sampling image segmentations for uncertainty quantification.
Medical Image Analysis , 34, 2016.[17] Y. LeCun, Y. Bengio, and G. Hinton. Deep learning.
Nature , 521(7553):436–444, May 2015.[18] A. L. Maas, A. Y. Hannun, and A. Y. Ng. Rectifier nonlinearities improve neural network acoustic moels.
ICML ,2013.[19] G. Maher, D. Parker, N. Wilson, and A Marsden. Neural network vessel lumen regression for automatedsegmentation in cardiovascular image-based modeling.
Cardiovascular Engineering and Technology , 2020.[20] M. E. Moghadam, Y. Bazilevs, T. Y. Hsia, I. E. Vignon-Clemente, A. L. Marsden, and (MOCHA) Alliancefor the Modeling of Congenital Hearts Alliance. A comparsion of outlet boundary treatments for prevention ofbackflow divergence with relevance to blood flow simulations.
Computational Mechanics , 48, 2011.[21] M. E. Moghadam, Y. Bazilevs, and A. L. Marsden. A bi-partitioned iterative algorithm for solving linear systemsarising from incompressible flow problems.
Computer Methods in Applied Mechanics and Engineering , 286, 2015.[22] M. E. Moghadam, I. E. Vignon-Clementel, R. Figliola, A. L. Marsden, and (MOCHA) Alliance for the Model-ing of Congenital Hearts Alliance. A modular numerical method for implicit 0D/3D coupling in cardiovascularfinite element simulations.
Journal of Computational Physics , 244, 2013.[23] T. Nair, D. Precup, D. L. Arnold, and T. Arbel. Exploring uncertainty measures in deep networks for multiplesclerosis lesion detection and segmentation.
Medical Image Analysis , 59, 2020.[24] S. Quicken, W. P. Donders, M. J. van Disseldorp, K. Gashi, B. M. E. Mees, F. N. van de Vosse, R. G. P. Lopata,T. Delhaas, and W. Huberts. Application of an adaptive polynomial chaos expansion on computationally ex-pensive three-dimensional cardiovascular models for uncertainty quantification and sensitivity analysis.
Journalof Biomechanical Engineering , 138, 2016.[25] S. Sankaran, L. J. Grad, and C. A. Taylor. Real-time sensitivity analysis of blood flow simulations to lumensegmentation uncertainty.
Conference on Medical Imaging and Computer Assisted Intervention , 2014.[26] S. Sankaran, L. Grady, and C. A. Taylor. Fast computation of hemodynamic sensitivity to lumen segmentationuncertainty.
IEEE Transactions on Medical Imaging , 34, 2015.[27] S. Sankaran, L. Grady, and C. A. Taylor. Impact of geometric uncertainty on hemodynamic simulations usingmachine learning.
Computer Methods in Applied Mechanics and Engineering , 297, 2015.[28] S. Sankaran, H. J. Kim, G. Choi, and C. A. Taylor. Uncertainty quantification in coronary blood flow simulations:impact of geometry, boundary conditions and blood viscosity.
Journal of Biomechanics , 49, 2016.[29] Sethuraman Sankaran and Alison L. Marsden. A stochastic collocation method for uncertainty quantificationand propagation in cardiovascular simulations.
Journal of Biomechanical Engineering , 133(3):1–12, 2011.[30] D. E. Schiavazzi, A. Doostan, G. Iaccarino, and A. L. Marsden. A generalized multi-resolution expansion foruncertainty propagation with application to cardiovascular modeling.
Computer Methods in Applied Mechanicsand Engineering , 314:196–221, 2017.[31] D.E. Schiavazzi, G. Arbia, C. Baker, A.M. Hlavacek, T-Y. Hsia, A.L. Marsden, I.E. Vignon-Clementel, andThe Modeling of Congenital Hearts Alliance (MOCHA) Investigators. Uncertainty quantification in virtualsurgery hemodynamics predictions for single ventricle palliation.
International journal for numerical methodsin biomedical engineering , 32(3):e02737, 2016.[32] J. Schmidhuber. Deep learning in neural networks: an overview.
Neural Networks , 61:85–117, January 2015.[33] J. Seo, D.E. Schiavazzi, A.M. Kahn, and A.L. Marsden. The effects of clinically-derived parametric datauncertainty in patient-specific coronary simulations with deformable walls.
International Journal for NumericalMethods in Biomedical Engineering , n/a(n/a):e3351, 2020.2534] J. Seo, D.E. Schiavazzi, and A.L. Marsden. Performance of preconditioned iterative linear solvers for cardiovas-cular simulations in rigid and deformable vessels.
Computational Mechanics , 64(3):717–739, 2019.[35] H. Si. TetGen, a Delaunay-based quality tetrahedral mesh generator.
ACM Transactions on MathematicalSoftware , 41, 2015.[36] I. A. Sigal, M. R. Hardisty, and C. M. Whyne. Mesh-morphing algorithms for specimen-specific finite elementmodeling.
Journal of Biomechanics , 41, 2008.[37] D. A. Steinman and F. Migliavacca. Editorial: Special issue on verification, validation, and uncertainty quan-tification of cardiovascular models: Towards effective VVUQ for translating cardiovascular modelling to clinicalutilitys.
Cardiovascular Engineering and Technology , 9, 2018.[38] C. Szegedy, W. Liu, Y. Jia, P. Sermanet, S. Reed, D. Anguelov, D. Erhan, V. Vanhoucke, and A. Rabinovich. Go-ing deeper with convolutions. In
Proceedings of the IEEE conference on computer vision and pattern recognition ,pages 1–9, 2015.[39] J.S. Tran, D.E. Schiavazzi, A. Ramachandra Bangalore, A.M. Kahn, and A.L. Marsden. Automated tuningfor parameter identification and uncertainty quantification in multi-scale coronary simulations.
Computers andFluids , 142:128–138, 2017.[40] Justin S. Tran, Daniele E. Schiavazzi, Andrew M. Kahn, and Alison L. Marsden. Uncertainty quantification ofsimulated biomechanical stimuli in coronary artery bypass grafts.
Computer Methods in Applied Mechanics andEngineering , 345:402–428, 2019.[41] A. Updegrove, N. Wilson, J. Merkow, H. Lan, A.L. Marsden, and S.C. Shadden. SimVascular: an open sourcepipeline for cardiovascular simulation.
Annals of Biomedical Engineering , 61, 2013.[42] A. Updegrove, N. M. Wilson, and S. C. Shadden. Boolean and smoothing of discrete surfaces.
Advances inEngineering Software , 95, 2016.[43] P. Venugopal, X. Li, L. Cheng, R. Bhagalia, and P. M. Edic. Sensitivity of FFR-CT to manual segmentation.
Proceedings of SPIE , 2019.[44] I. E. Vignon-Clementel, C. A. Figueroa, K. E. Jansen, and C. A. Taylor. Outflow boundary conditions for3D simulations of non-periodic blood flow and pressure fields in deformable arteries.
Computer Methods inBiomechanics and Biomedical Engineering , 13, 2009.[45] N.M. Wilson, A.K. Ortiz, and A.B. Johnson. The vascular model repository: a public resource of medicalimaging data and blood flow simulation results.
Journal of medical devices , 7(4), 2013.[46] H. Xu, D. Baroli, and A. Veneziani. Global sensitivity analysis for patient-specific aortic simulations: the roleof geometry, boundary condition and les modeling parameters.