LDDMM stochastic interpolants: an application to domain uncertainty quantification in hemodynamics

We introduce a novel conditional stochastic interpolant framework for generative modeling of three-dimensional shapes. The method builds on a recent LDDMM-based registration approach to learn the conditional drift between geometries. By leveraging th…

Authors: Sarah Katz, Francesco Romor, Jia-Jie Zhu

LDDMM stochastic interpolants: an application to domain uncertainty quantification in hemodynamics
LDDMM STOCHASTIC INTERPOLANTS: AN APPLICA TION TO DOMAIN UNCER T AINTY QUANTIFICA TION IN HEMOD YNAMICS SARAH KA TZ 1 , FRANCESCO ROMOR 1 , JIA-JIE ZHU 2 , AND ALFONSO CAIAZZO 1 Abstract. W e in troduce a no v el conditional stochastic interpolant framework for generative mo deling of three-dimensional shapes. The method builds on a recent LDDMM-based registration approach to learn the conditional drift betw een geometries. By leveraging the resulting pull-back and push-forw ard operators, w e extend this formulation beyond standard Cartesian grids to complex shapes and random v ariables defined on distinct domains. W e present an application in the context of cardio v ascular simu- lations, where aortic shapes are generated from an initial cohort of patients. The conditioning v ariable is a latent geometric representation defined by a set of cen terline p oin ts and the radii of the corresp onding inscribed spheres. This metho dology facilitates b oth data augmentation for three-dimensional biomedi- cal shap es, and the generation of random p erturbations of con trolled magnitude for a given shap e. These capabilities are essen tial for quan tifying the impact of domain uncertainties arising from medical image segmentation on the estimation of relev ant biomarkers. Contents 1. In tro duction 1 2. Preliminaries 3 2.1. Sto c hastic in terpolants 3 2.2. LDDMM for image registration 5 2.3. ResNet LDDMM for registration of surface meshes 6 3. Conditional LDDMM sto c hastic in terpolant 6 3.1. Sto c hastic in terpolant betw een n umerical solutions on different shapes 6 3.2. Generativ e mo deling with conditional LDDMM stochastic in terpolants 7 4. In ter-patient n umerical modeling of aortic blo od flow 9 4.1. A cquisition of the dataset 10 4.2. Multi-grid ResNet-LDDMM registration of aortic geometries 10 4.3. Computational hemo dynamics 11 4.4. Preserv ation of hexahedral mesh hierarc hies structure upon registration 12 5. Application to shap e uncertain t y quan tification in cardio v ascular models 13 5.1. Numerical metho ds for geometric uncertain t y quan tification 13 5.2. Domain p erturbations for aortic geometries 14 6. Numerical results 15 6.1. Sim ulation setup 15 6.2. Shap e uncertain t y quan tification on selected biomark ers 18 6.3. Sample size study 23 7. Conclusions 27 A ckno wledgemen ts 27 References 27 App endix A. T raining of the conditional LDDMM stochastic in terp olan t 31 App endix B. Smoothing algorithms for mesh asp ect ratio optimization 31 App endix C. A dditional details on biomark er statistics 33 1. Intr oduction Numerical metho ds for the simulation of ph ysical systems rely on mathematical mo dels and on the definition and generation of computational domains that can accurately represent the problem of in terest. In many applications, the definition of the domain itself might b e affected by uncertainties due, e.g., to 1 WIAS, W eierstrass Institute for Applied Analysis and Sto c hastics, Anton-Wilhelm-Amo-Str. 39, 10117 Berlin, Germany . 2 KTH Roy al Institute of T echnology , Lindstedtsv ägen 25, 114 28 Stockholm, Sweden. 1 2 KA TZ, R OMOR, ZHU, AND CAIAZZO measuremen t noise or pro cess-dep enden t errors making it necessary to quantify and manage how these affect the results and the relev ant quan tities of interest. This work is motiv ated b y the imp ortance of these domain uncertainties in blo o d flow simulations, where CFD is often used in com bination with a v ailable medical data (e.g., medical images and measure- men ts) to infer biomarkers to diagnose and monitor cardio v ascular pathologies. Quantifying the sensitivit y of n umerical results to domain v ariability is driven by t w o main asp ects. On the one hand, cardiov ascular mo dels rely on anatomical data, suc h as MRI or CT scans, whic h are (man ually , semi-automatically , or automatically) processed to create surface and v olume meshes. P ossible errors and biases in tro duced in this process result in unav oidable aleatoric and epistemic uncertainties in the computational domain used for the simulations. On the other hand, understanding the influence of in ter-patient v ariability is extremely imp ortan t to draw conclusions relating different hemo dynamics quantities b eyond a single patien t. Classical uncertaint y quantification (UQ) metho ds aim at appro ximating the statistics of quantities of in terest dep ending on the distribution of unkno wn mo del inputs. In the case of parametric uncertain ties, this task can b e addressed defining proper sampling strategies of a suitable parameter space. Ho wev er, handling the geometrical uncertainties requires suitable methods for generating instances of anatomical shap es, which cannot b e describ ed by a low num b er of geometrical parameters. F ew w orks hav e addressed this issue limited to the v ariability of lumen diameter, in the case of coronary arteries [ 45 ] using an ide- alized geometry and multiple patient data, and for the ascending aorta [ 37 ], generating domain samples from tw o-dimensional cross sections obtained during the segmen tation. A more systematic approac h w as recen tly proposed in [ 8 ], generating multiple samples of ascending aorta from an underlying description based on few geometrical parameters, and including a surface p erturbation defined via Gaussian fields. Other methods to mo del in ter-patient shap e v ariability from a more general persp ectiv e considered sam- pling based on Statistical Shape Mo dels (SSMs, see, e.g., [ 1 ]), a technique used to c haracterize a patient p opulation with a low er-dimensional parametric shap e model which can b e com bined with a probabilistic mo del to generate new domains [ 24 , 47 , 48 ]. The idea of a machine learning generativ e mo del based on sampling in a low-dimensional laten t space has b een very recently explored in [ 46 ], learning the la- ten t representation from diffeomorphic maps b et w een an initial set of shap es, and in [ 18 ] considering a v ariational auto-enco der for the node co ordinates. Generativ e mo deling aims at learning a mapping b etw een a simple prior distribution (e.g., a Gaussian) and a complex target distribution kno wn only through samples, enabling the generation of new realistic instances. This paradigm has achiev ed remark able success in image generation, where, for instance, the conditioning v ariable can b e a text prompt guiding the synthesis pro cess. Among the most recent and effectiv e approaches, conditional flow matching [ 33 ], rectified flows [ 34 ], and sto c hastic interpolants [ 2 ] ha ve b een introduced as frameworks for generative mo deling based on learning contin uous-time transp ort maps betw een distributions. These m ethods construct a time-dep enden t velocity field (or drift) that con tinuously transforms the prior in to the target distribution, offering adv an tages in terms of training stabilit y and sampling efficiency compared to classical diffusion models. Extensions and theoretical connections to broader classes of stochastic pro cesses ha v e also been explored [ 36 , 16 ]. While these tec hniques hav e b een extensively developed and applied in tw o-dimensional settings such as image [ 43 ] or video [ 35 ] synthesis, their application to three-dimensional geometric data, and in particular to the generation of anatomical shap es, remains largely unexplored. W e prop ose a new generativ e mo deling metho d for three-dimensional shap es based on conditional sto c hastic in terp olan t and diffeomorphic shape registration. The approach represen ts shap es as random v ariables in the three-dimensional Euclidean space. Using av ailable data, the metho d learns the drift of a sto c hastic differential equation (SDE) defining a sto c hastic process ( sto chastic interp olant ) b et ween differen t distributions, depending on a conditioning v ariable in a latent space. New shap es are then generated b y sampling the conditioning v ariable and evolving the corresponding SDE. The main nov elty concerns the usage of shap e registration maps within the training phase to handle random v ariables defined on differen t physical domains. T o this purp ose, we used a multilev el ResNet- LDDMM metho d [ 44 ], which can robustly and efficien tly compute registration maps b et ween anatomical shap es. Analogously to optimal transp ort-based sto c hastic interpolant methods, the paths generated b y LDDMM sto chastic in terp olan t minimize the LDDMM energy functional. This improv es the sampling efficiency and the qualit y of the generated shapes, as the training process similarly to optimal transp ort- based sto c hastic in terp olan t metho ds, is guided by the LDDMM energy functional, which is a natural metric for shap e registration. In the context of computational fluid dynamics, a further asp ect which cannot b e neglected concerns the generation of computational meshes on the deformed domains. In the numerical examples shown in this paper forward simulations are based on a high-order discontin uous Galerkin solver [ 5 ], which relies CONDITIONAL LDDMM STOCHASTIC INTERPOLANTS 3 Conditional LDDMM stochastic interpolant F orwar d domain uncertainty quantification LDDMM minimal path V elocity and pr essur e statistics Biomark ers' statistics Application 0 5 10 15 20 25 29 0.1 0.2 0.3 OSI Figure 1. Overview of the prop osed metho d and application in forw ard domain uncer- tain ty quan tification for aortic blo od flows. on a hierarc hy of hexahedral grids. This issue is addressed by defining a suitable extension of the surface deformation on to the whole three-dimensional mesh b y iteratively solving a non-linear elasticity problem. W e test the framework in the con text of uncertaint y quantification for aortic blo od flows, v alidating the capabilit y of the LDDMM sto c hastic interpolant to sample the shape space. The rest of the paper is organized as follows. Section 2 in tro duces the background of generative mo d- eling based on sto c hastic interpolants and of shap e registration using LDDMM and ResNet-LDDMM. Section 3 fo cuses on the no v el LDDMM sto c hastic interpolant metho d, whilst section 4 describ es the mathematical mo dels and the n umerical methods for blo od flo w simulations. The application to geomet- rical uncertaint y quantification is detailed in section 5 , numerical results are presented in section 6 , and section 7 dra ws the main conclusions. 2. Preliminaries This section focuses on the bac kground for introducing the LDDMM sto c hastic interpolant framew ork from a general p oin t of view, reviewing the k ey concepts of sto c hastic in terp olan ts in se ction 2.1 and the LDDMM metho d in section 2.2 . Section 2.3 discusses the application of the recen t multilev el ResNet- LDDMM in the case of registration of three-dimensional computational meshes. 2.1. Sto c hastic in terp olan ts. Sto c hastic interpolants are a mathematical framework that connects differen t approac hes for generative mo deling, such as diffusion mo dels, flow matching, and con tinuous normalizing flows. The main idea b ehind sto c hastic interpolants is to bridge , in a contin uous wa y , tw o differen t probabilit y distributions via a time-dep enden t random pro cess, whic h, once learned, can be used to generate new samples. In what follows, let us denote the probability space as ( X , A , P ) , comp osed of a sample space X , a σ -algebra of even ts A , and a probabilit y measure P . F or a random v ariable X : ( X , A , P ) → ( R n , B ( R n )) , where B stands for the Borel σ -algebra, w e denote b y µ X the probability measure induced b y X on R n and b y p ( x ) its probability densit y function, when it exists. Definition 2.1 (Sto chastic in terp olan t) . Let F = {F t } t ∈ [0 , 1] b e a filtration supported on ( X , A , P ) , and let X 0 , X 1 : ( X , A , P ) → ( R n , B ( R n )) b e tw o random v ariables. A sto c hastic interpolant b et ween X 0 and X 1 is a sto c hastic pro cess I : ([0 , 1] × R n , B ([0 , 1]) ⊗ F ) → ( R n , B ( R n )) , adapted to F , that allo ws one to bridge t wo random v ariables X 0 and X 1 , suc h that: I t =0 = X 0 , I t =1 = X 1 , P - a.s. Let there b e given tw o data distributions, which w e assume are represen ted b y tw o random v ariables X 0 and X 1 , with probability densities p ( x 0 ) and p ( x 1 ) with resp ect to the Leb esgue measure. The sto c hastic in terp olan t b etw een X 0 and X 1 is mo deled as the solution of an SDE ( dI t = b θ t ( I t ) dt + σ t dW t , t ∈ [0 , 1] , I 0 ∼ p ( x 0 ) , I 1 ∼ p ( x 1 ) , (1) where W t is a standard Brownian motion, σ t is a time-dep enden t diffusion coefficient, and the time- dep enden t drift b θ t : R n → R n needs to b e learned minimizing the follo wing sto chastic interp olant loss 4 KA TZ, R OMOR, ZHU, AND CAIAZZO function: L ( b θ ) = Z 1 0 E I t [ ∥ b θ t ( I t ) − u t ( I t ) ∥ 2 ] dt, (2) dep ending on a parameter θ , for a prescrib ed drift u t : R n → R n . T ypically , during the training, the drift u t is conditioned on x 0 and x 1 (e.g., tw o images as in Figure 2 ) sampled from X 0 and X 1 , resp ectiv ely , and is giv en in the form: u t := u cond t ( x 0 , x 1 ) = ˙ α t x 0 + ˙ β t x 1 + ˙ σ t W t , t ∈ [0 , 1] , (3) where α t , β t , σ t : R → R are deterministic functions satisfying α 0 = 1 , α 1 = 0 , β 0 = 0 , β 1 = 1 , σ 0 = 0 , σ 1 = 0 . The aim is to approximate a stochastic in terp olan t ˜ I t b et ween X 0 and X 1 as ˜ I t = α t X 0 + β t X 1 + σ t W t , t ∈ [0 , 1] . (4) A common c hoice is α t = 1 − t , β t = t , and σ t = σ p t (1 − t ) for some fixed parameter σ > 0 . Considering the conditioned drift, the loss function ( 2 ) b ecomes L ( b θ ) = Z 1 0 E X 0 ,X 1 [ ∥ b θ t ( I t ) − u cond t ( X 0 , X 1 ) ∥ 2 ] dt, (5) with the v ector field u cond t defined in equation ( 3 ). The definitions of the sto c hastic interpolant in the form ( 4 ), and the corresp onding conditional drift ( 3 ), require computing the sums and differences of discretized realizations of tw o random v ariables X 0 and X 1 in R n . This op eration can b e naturally defined in the case of images of the same size, discretized, e.g., as tw o- or three-dimensional grids with a one-to-one corresp ondence betw een pixels (see Figure 2 , left). How ever, it cannot b e directly extended to the case of random v ariables mo deling ph ysical fields on general shapes (see Figure 2 , righ t), which might b e defined on p oin t clouds with differen t cardinalities or through unstructured (surface or volume) meshes with different num bers of vertices and elemen ts. A generalization to this case, and the application in the context of n umerical metho ds for PDEs, will b e discussed in section 3 . Figure 2. Left: On a t wo-dimensional grid, the prescrib ed velocity field u cond t ( 3 ) can b e defined by computing the difference b et ween the realizations of the random v ariables X 0 and X 1 . Right: In the case of a three-dimensional shap e (described, e.g., as a surface mesh), a definition of u cond t based on the difference ( 3 ) is no longer p ossible, since there is no direct corresp ondence b et ween p oin ts. R emark 2.2 (Equiv alent ODE flow formulation) . W e can rewrite the sto c hastic in terp olan t formulation with resp ect to its probability density p I t as ( ∂ t p I t + ∇ · ( p I t v SI t ) = 0 , t ∈ [0 , 1] , p I t =0 = p X 0 , p I t =1 = p X 1 , (6) whic h is equiv alent to equation ( 5 ) via its corresp onding F okk er-Planck equations, with v SI t = b t − σ 2 t 2 ∇ log p I t . Sampling from the sto c hastic interpolant then requires solving the ODE ( ∂ t ϕ SI t ( x ) = v SI t ( ϕ SI t ( x )) , t ∈ [0 , 1] , ϕ SI 0 ( x ) = x , (7) for a map ϕ SI t : [0 , 1] × R n → R n , and setting p I t as the push-forward of p X 0 with resp ect to the map ϕ SI , i.e., p I t = p X 0 ◦  ϕ SI t  − 1 . CONDITIONAL LDDMM STOCHASTIC INTERPOLANTS 5 The existence of the solution ϕ SI of ( 7 ) is guaranteed under suitable regularity assumptions on v SI t , e.g., when v SI ∈ L 2 ([0 , 1] , L 2 ( R n )) [ 3 ] and: Z [0 , 1] sup x ∈ B ∥ v SI t ( x ) ∥ + Lip ( v SI t , B ) dt < ∞ , ∀ B ⊂ R n compact . (8) R emark 2.3 (Regularizers for sto c hastic interpolants) . The definition ( 4 ) is limited to linear interpolants b et ween X 0 and X 1 , which ma y not b e optimal in terms of minimizing the kinetic energy densit y . F or this purp ose, the minimization can b e formulated as W 2 2 ( p X 0 , p X 1 ) = inf p I t ,v SI t Z 1 0 E I t [ ∥ v SI t ∥ 2 ] dt, s.t. p I t , v SI t satisfy ( 6 ), with E I t [ ∥ v SI t ∥ 2 ] = R R n ∥ v SI t ( x ) ∥ 2 dp t ( x ) where W 2 2 stands for the 2-W asserstein distance in the dynamic formulation (Benamou–Brenier) b et ween the probability densities p X 0 and p X 1 . When appro ximating physical fields obtained as solutions of PDEs, additional constrain ts ma y b e imp osed on the drift b t to enforce ph ysical consistency , e.g., divergence-free constrain ts for incompressible fluids or a prescrib ed kinetic energy rate of change (see, e.g., [ 39 ]). 2.2. LDDMM for image registration. The problem of image registration concerns finding a smo oth transformation b et ween a sour c e/template image (defined, e.g., as a p oin t cloud or a mesh) and a tar get image. The Large Deformation Diffeomorphic Metric Mapping (LDDMM) metho d [ 49 ] addresses this task by mo deling the transformation as the flow of an ordinary differential equation (ODE), ensuring smo oth and inv ertible mappings. Let us consider the source and the target ob jects denoted b y S and T , resp ectiv ely ( S, T ⊂ R d ), and defined by their c haracteristic functions χ S : R d → R and χ T : R d → R , resp ectiv ely . Without loss of generalit y , let us assume that b oth ob jects are contained in an op en b ounded set U , i.e., supp ( χ S ) , supp ( χ T ) ⊂ U ⊂ R d . The LDDMM metho d seeks a time-dependent v elo cit y field v LDDMM t : U → R d , such that the flo w ϕ LDDMM t : U → U , solution to ( ∂ t ϕ LDDMM t ( x ) = v LDDMM t ( ϕ LDDMM t ( x )) , t ∈ [0 , 1] , ϕ LDDMM 0 ( x ) = x , (9) smo othly transforms the source ob ject S into the target T . W e assume that v LDDMM ∈ L 2 ( I , H s ( U )) , Z I ∥ v LDDMM t ∥ 2 H s ( U ) dt < ∞ , (10) so the Sob olev embedding theorem implies v ∈ L 2 ( I , C 1 ,α ( U )) with α = s − d 2 − 1 > 0 . Hence, for eac h t ∈ [0 , 1] , the flow ϕ LDDMM t is a diffeomorphism, ϕ LDDMM ∈ C 0 ( I , Diff 1 ( U )) . The transformation of the source image at time t is given by: χ S t ( x ) = χ S (  ϕ LDDMM t  − 1 ( x )) , ∀ t ∈ [0 , 1] . The LDDMM metho d seeks to minimize the following augmen ted energy functional: E ( v ) = Z 1 0 ∥ v LDDMM t ∥ 2 H s ( U ) dt + λ Z U | χ S t ( x ) − χ T ( x ) | 2 d x , (11) where ∥ v LDDMM t ∥ H s ( U ) ensures smo othness of the v elo cit y field, and λ > 0 is a regularization parameter balancing the smo othness of the transformation and the fidelity to the target image. W e refer the reader to [ 17 ] for additional details on the existence of a minimizer with the required regularity . Definition 2.4 (Push-forward and pull-bac k op erators) . Given a diffeomorphism ϕ : U → U , we can define the push-forw ard operator ϕ # and the pull-bac k op erator ϕ # that allow us to transp ort fields defined on U via: ϕ # ( f ) = f ◦ ϕ − 1 , ϕ # ( f ) = f ◦ ϕ, for any f : U → R m . These op erators will b e useful in section 3 to define a map b et ween different computational domains. R emark 2.5 (Analogies b et ween LDDMM and diffusion mo dels) . The LDDMM metho d and sto c hastic in terp olan ts (section 2.1 ) are defined in different settings (deterministic and probabilistic) and address differen t problems. Ho wev er, both approac hes are based on flows of ODEs that aim to transform a source (sub domain or probability measure, resp ectiv ely) into a target through suitable time-dep endent 6 KA TZ, R OMOR, ZHU, AND CAIAZZO v elo cit y fields – compare equations ( 7 ) and ( 9 ), as w ell as the regularity assumptions ( 8 ) for the sto c hastic in terp olan t and ( 10 ) for the velocity field in LDDMM. 2.3. ResNet LDDMM for registration of surface meshes. Let us define a generic mesh as M = (Ω , m , E ) , where Ω ⊂ R d stands for the (physical) domain or a piecewise smo oth manifold asso ciated with it, m = { m i } N i =1 ⊂ R d is the arra y of the mesh vertices’ co ordinates, and E is the set of d -dimensional cells defined as subsets of vertices in m (for instance, E is a set of subsets of 4 v ertices for tetrahedral meshes in 3D). In what follo ws, we assume that there exist pairwise bi-Lipschitz homeomorphisms b et ween shap es. This assumption do es not affect the generality of the framework and is consistent with the setting where all the considered shapes represen t the same underlying physical system, e.g., a PDE on differen t domains with a given set of b oundaries and corresp onding b oundary conditions – see section 6 . How ever, w e do not assume an y a priori relation on the sizes of the discrete spaces (n um b er of vertices, faces, or elemen ts) of the different shap es. Instead of characteristic functions (as done in section 2.2 ), we now consider source and target shap es represen ted by computational meshes M S = (Ω S , m S , E S ) and M T = (Ω T , m T , E T ) , resp ectively . In this case, the minimization problem ( 11 ) is solv ed by minimizing a shap e dissimilarity measure b et ween the meshes. A common choice is the so-called Chamfer distance b et ween the p oin t clouds given by mesh v ertices, defined by d Chamfer ( m S , m T ) = X x ∈ m S min y ∈ m T ∥ x − y ∥ 2 + X y ∈ m T min x ∈ m S ∥ x − y ∥ 2 . (12) In this setting, the LDDMM registration problem can b e formulated as: find the time-dep enden t v elo cit y field v LDDMM t : U → R d that minimizes the energy functional: E ( v ) = Z 1 0 ∥ v LDDMM t ∥ 2 H s ( U ) dt + λ d Chamfer ( ϕ LDDMM 1 ( m S ) , m T ) , (13) where ϕ LDDMM t is the flow generated by v LDDMM t (solution of equation ( 9 )), and λ > 0 is a regularization parameter. Dep ending on the c hoice of the n umerical optimization metho d and the parametrization of the velocity field v LDDMM t , several v ariants of LDDMM can be defined. In what follo ws, we fo cus on the recent and efficient multilev el ResNet LDDMM [ 44 ] (based on [ 4 ]), in which the velocity field is parametrized using a multilev el ResNet architecture, and the minimization considers a generalization of the Chamfer distance ( 12 ), accounting for dissimilarity b et w een b oundary elements and mesh centerlines, and with regularization terms to preserve mesh qualit y during the deformation (see [ 44 ] for details). In this case, the regularity of the registration map c hanges from the space of diffeomorphisms ϕ LDDMM t ∈ Diff 1 ( U ) to the space of bi-Lipschitz homeomorphisms, which is more suitable for applications in which the shap es to b e registered are not necessarily diffeomorphic (e.g., due to the presence of corners or edges). The in vertibilit y constrain t is not implicitly imp osed on the ResNet architecture, but it can b e chec k ed a p osteriori [ 44 ]. Moreov er, the size of the mesh is refined during the training stage, allo wing for efficient computation and scalability to large datasets of 3D shap es. 3. Conditional LDDMM stochastic interpolant 3.1. Sto c hastic interpolant b etw een n umerical solutions on differen t shap es. This work is mo- tiv ated b y the application of generative mo deling in the context of the numerical solution of PDEs, in whic h the computational domain is discretized b y a mesh and the solution is sought in a discrete space, e.g., of piecewise p olynomials defined on the mesh elements. Building on the framework of sto c hastic interpolants (section 2.1 ), we consider data distributions and random v ariables represen ting domain co ordinates or numerical solutions of a giv en PDE, supported on a p o ol of contin uous, three-dimensional domains, denoted as { Ω i } N shapes i =0 . In this case, it is necessary to define a map b et w een different computational domains on which differen t fields are supp orted. This is a crucial difference with resp ect to the application to the more classical case of image generation, where data structures are naturally defined on ob jects of the same size (e.g., images). Let us assume that there exists an op en, b ounded, and connected set U ⊂ R 3 suc h that Ω i ⊂ U for i = 0 , . . . , N shapes , and let us introduce a background Hilb ert space H ( U ) , e.g., a functional space on U , on which all the numerical solutions corresp onding to the different computational domains are defined. Among the a v ailable domains, let Ω 0 ⊂ R 3 denote the template (or reference) shap e, and let M 0 = (Ω 0 , m 0 , E 0 ) b e the corresp onding computational mesh. Using the notation of the probabilistic framework in tro duced in subsection 2.1 , we consider a source random v ariable X 0 : ( X , A , P ) → H ( U ) supp orted on Ω 0 , whose realizations x 0 ∈ H ( U ) b elong to the background Hilb ert space, and whose discretizations CONDITIONAL LDDMM STOCHASTIC INTERPOLANTS 7 Multilevel R esNet-LDDMM L evel 0 L evel 1 L evel 2 Sour ce mesh 0 3000 4000 5000 Number of epochs 10 5 10 4 10 3 Chamfer loss Figure 3. Example of application of the multilev el ResNet-LDDMM [ 44 ] for the reg- istration of computational meshes of aortic shap es. In this example, meshes are refined after 3,000, 4,000, and 5,000 ep o c hs. x h 0 ∈ R n 0 ha ve the dimension ( n 0 ) of the space of piecewise p olynomial functions (e.g., basis functions of a finite element metho d) defined on the mesh M 0 . The target random v ariable X 1 : ( X , A , P ) → H ( U ) will represent an y of the remaining shap es Ω i ⊂ R 3 , i ∈ { 1 , . . . , N shapes } , with realizations x i ∈ H ( U ) , and discretizations x h i ∈ R n i supp orted on the corresp onding meshes M i = (Ω i , m i , E i ) with n i degrees of freedom. As noted in section 2.3 , no relation is assumed a priori on the sizes of the discretizations corresp onding to different shap es; in particular, n i  = n j for i  = j . R emark 3.1 (General realizations of X 0 and X 1 ) . The framework has so far b een introduced from a general p ersp ectiv e, i.e., considering that the realizations x 0 , x 1 ∈ H ( U ) of X 0 and X 1 , and the corre- sp onding discretizations x h 0 ∈ R n 0 and x h i ∈ R n i , can represent any scalar or vector fields defined on the computational mesh. In particular, the distributions could represent numerical solutions of an incom- pressible flow (e.g., velocity and pressure), a deformation field defined on the computational mesh, or the co ordinates of the supp ort p oin ts of the finite elemen t spaces. Section 5 will fo cus on the latter case, as an application of shap e generative mo deling. R emark 3.2 . In this study , for the purp ose of v alidating the conditional LDDMM sto c hastic interpolant metho d, we restrict ourselves to a fixed, single template domain. In general, it is p ossible to consider v arying domains Ω j ⊂ R 3 also for the random v ariable X 0 , which will increase the ov erall computational effort required. 3.2. Generativ e mo deling with conditional LDDMM sto c hastic in terp olan ts. W e prop ose a new generativ e mo del that combines sto c hastic interpolants (section 2.1 ) with ResNet-LDDMM registration (section 2.3 ), in order to define a suitable conditioned drift for data distributions on v ariable domains. Firstly , let us observe that, to apply the framework introduced in section 2.1 , it is necessary to extend the definition of noise to the background space H ( U ) , i.e., to define a Q -Wiener pro cess W Q t on H ( U ) , with the cov ariance op erator Q : H ( U ) → H ( U ) b eing a symmetric, positive, trace-class op erator. W e consider an extension based on the standard theory of Q -Wiener pro cesses in Hilb ert spaces; see [ 13 ] for more details. W e denote b y C : ( X , A , P ) → ( R c , B ( R c )) , c > 0 , an additional c onditioning random v ariable that will trigger the generative pro cess. In the field of image [ 43 ] or video [ 35 ] generation, for example, C t ypically represents the embedding of a text prompt. Let t w o samples x 0 and x i of t w o random v ariables X 0 and X 1 resp ectiv ely , b e given, supp orted on the meshes M 0 (template shape) and M i . The generative mo del follows the approach describ ed in section 8 KA TZ, R OMOR, ZHU, AND CAIAZZO 2.1 , considering the sto chastic interpolant solution to the SDE ( dI t = b θ t ( I t , C ) dt + σ t dW Q t , I t =0 = X 0 , I t =1 = X 1 , (14) where the co efficient b θ t ( I t , C ) is trained by minimizing the loss function L ( b θ ) = Z 1 0 E X 0 ,X 1 ,C [ ∥ b θ t ( I t , C ) − u cond t ( X 0 , X 1 ) ∥ 2 ] dt. (15) for a suitable conditional drift u cond t . Unlike the case discussed in section 2.1 , prescribing u cond : [0 , 1] × Conditional LDDMM stochastic interpolant Conditional LDDMM vector fields magnitude [m] LDDMM path arbitrary path Figure 4. Graphical ov erview of the prop osed conditional LDDMM sto c hastic inter- p olan t metho d. The map defining the conditional drift b θ as a function of the condition- ing v ariable C is learned registering the av ailable shap es onto a common reference ( x 0 ), in order to find a sto chastic in terp olan t b et ween the tw o distributions. H ( U ) × H ( U ) → H ( U ) , requires defining a time-dep enden t velocity field u cond t ( x 0 , x i ) on an interme diate domain (at time t ) b et w een Ω 0 ⊂ U and Ω i ⊂ U . This task will b e addressed by computing the LDDMM registration map ϕ 0 ,i t : R 3 → R 3 b et ween the meshes M 0 and M i (i.e., b et ween the corresponding p oin t clouds m 0 and m i ) with the m ultilev el ResNet LDDMM framew ork (see section 2.3 ). Namely , if the realizations x 0 and x i represen t sets of co ordinates within the computational domains (the case considered in this work), we set u cond t ( x 0 , x i ) = v 0 ,i t ( ϕ 0 ,i t ( x 0 )) + σ t ˙ W Q t , t ∈ [0 , 1] , (LDDMM flow matching) (16) where v 0 ,i t is the LDDMM velocity field obtained by registering M 0 to M i , equation ( 9 ) (see also figure 5 , left) for the training of b θ t ( I t , C ) ( 15 ) for a given realization of the conditioning v ariable. Next, we sample different solutions I t [ c ] of ( 14 ) and define the new deformed domain as E [ I t [ c ]] . The training and generativ e phases of the prop osed conditional LDDMM sto c hastic interpolant metho d are summarized in Algorithms 1 and 2 , resp ectiv ely . More details on the training pro cedure for the conditioning drift are pro vided in app endix A . R emark 3.3 (Generalization to ph ysical fields on meshes) . If the random v ariables describe physical fields defined on the computational mesh (e.g., piecewise p olynomial numerical solutions of PDEs), one can first introduce maps betw een the discrete spaces on the meshes M 0 and M i based on the push-forw ard ( ϕ 0 ,i t ) # and pull-back ( ϕ 0 ,i t ) # op erators induced by ϕ 0 ,i t , in order to transp ort fields from M 0 to M i and vice versa. The conditional drift can b e defined via u cond t ( x 0 , x i ) = ( ψ 0 ,i t ) # ( x i ) − ( ϕ 0 ,i t ) # ( x 0 ) + σ t ˙ W Q t , t ∈ [0 , 1] , (17) CONDITIONAL LDDMM STOCHASTIC INTERPOLANTS 9 Figure 5. The conditional drift u cond t ( x 0 , x i ) is defined based on the LDDMM regis- tration. When x 0 and x i represen t domain co ordinates (left), the conditional drift is computed from the registration field ϕ 0 ,i b et ween the tw o realizations. When x 0 and x i represen t physical fields/solutions of PDEs (right), the conditional drift is computed, at eac h intermediate step t ∈ [0 , 1] , considering the push-forward of ϕ 0 ,i (registration of x 0 on to x i ) and the pull-back of ψ 0 ,i (registration of x i on to x 0 ). Algorithm 1 Conditional LDDMM Sto c hastic Interpolant – T raining Phase Require: T emplate mesh M 0 , target meshes {M i } N shapes i =1 , conditioning v ariables { c i } N shapes i =1 Ensure: T rained drift netw ork b θ t 1: for ep och = 1 to N epo chs do 2: Sample batch of indices B ⊂ { 1 , . . . , N shapes } 3: for i ∈ B do 4: Compute registration map ϕ 0 ,i t b et ween M 0 and M i and velocity v 0 ,i t (using the multilev el ResNet LDDM M – section 2.3 ) 5: Sample time: t ∼ Uniform (0 , 1) 6: Compute in termediate state: I t = ϕ 0 ,i t ( x 0 ) 7: Compute conditional drift: u cond t ( x 0 , x i ) = v 0 ,i t ( ϕ 0 ,i t ( x 0 )) + σ t ϵ , where ϵ ∼ N (0 , I ) 8: Compute loss: L i = ∥ b θ t ( I t , c i ) − u cond t ( x 0 , x i ) ∥ 2 2 9: end for 10: Up date parameters: θ ← θ − ∇ θ 1 |B| P i ∈B L i 11: end for 12: return b θ t Algorithm 2 Conditional LDDMM Sto c hastic Interpolant – Generative Phase Require: T emplate sample x 0 on mesh M 0 , conditioning v ariable c , trained drift netw ork b θ t , time discretization { t k } K k =0 with t 0 = 0 , t K = 1 Ensure: Generated sample E [ I t [ c ]] | t =1 on deformed domain 1: Initialize: I 0 ← x 0 2: Initialize: S ← 0 3: for s = 1 to N samples do 4: Initialize: I ( s ) 0 ← x 0 5: for k = 0 to K − 1 do 6: Compute drift: b ( s ) k = b θ t k ( I ( s ) t k , c ) 7: Sample noise: ξ ( s ) k ∼ N (0 , I ) 8: Up date state: I ( s ) t k +1 = I ( s ) t k + b ( s ) k · ( t k +1 − t k ) + σ t k √ t k +1 − t k · ξ ( s ) k 9: end for 10: A ccumulate sample: S ← S + I ( s ) t K 11: end for 12: return E [ I t [ c ]] | t =1 ≈ 1 N samples S where ψ 0 ,i t : Ω i ⊂ U → Ω 0 ⊂ U is defined as the backw ard transformation ψ 0 ,i t := ϕ i, 0 − t (figure 5 , right). 4. Inter-p a tient numerical modeling of a or tic blood flo w The conditional LDDMM sto c hastic interpolant framework introduced in the previous section can b e applied to an y n umerical model to handle v ariability of computational domains related via bi-Lipsc hitz 10 KA TZ, R OMOR, ZHU, AND CAIAZZO homeomorphisms. T o illustrate its capabilities, we apply the prop osed metho d to data augmentation and geometric uncertaint y quantification for the case of aortic hemo dynamics. This section introduces the exp erimen tal setup, numerical metho ds, and data sources underpinning the results in the following sections 5 and 6 . Specifically , w e introduce the dataset of aortic geometries (section 4.1 ) and the applica- tion of ResNet-LDDMM for registration (section 4.2 ). W e then detail the mathematical mo del and solv er calibration (section 4.3 ), concluding with the strategy for transp orting high-quality hexahedral meshes (section 4.4 ). 4.1. A cquisition of the dataset. F or the training of the conditioned LDDMM stochastic in terp olan t generativ e mo del, we emplo y a dataset S = { Ω i } i =0 ,n S − 1 consisting of n S = 1261 aortic synthetic shap es ( 1209 for training and 52 for testing) generated with statistical shap e mo deling (SSM) from an initial dataset of 228 patient surfaces. The patient dataset was acquired with 3D steady-state free-precession (SSFP) magnetic resonance imaging (MRI) (vo xel resolution 2 mm × 2 mm × 4 mm, resolution used for surface reconstruction 1 mm × 1 mm × 2 mm). The SSM is built on a cen terline enco ding, describing eac h geometry by a fixed num b er of centerline p oin ts ( n cntrl = 390 ) and the corresp onding radii of the largest inscrib ed spheres with cen ter on those p oin ts. Hence, the enco ding vector has size 390 × 4 . All shap es are top ologically equiv alen t and represen t a common p ortion of the aorta, containing aortic arch, thoracic aorta, and three additional branches. F or a complete description of the pro cedure, we refer the reader to [ 24 , 47 , 48 ]. A reduced version of this dataset has been recently used in [ 44 ] to v alidate the ResNet-LDDMM method for shap e registration. Figure 6 (left) sho ws the template geometry selected for the scop e of this work (surface mesh). Figure 6 (righ t) sho ws the test and training shap es in a t w o-dimensional represen tation as the result of a clustering with a t-SNE [ 50 ] em b eddings based on the similarity matrix D i =  d i ij  of the displacements from the template, i.e., d i ij := ∥ ϕ 0 ,i 1 − ϕ 0 ,j 1 ∥ L 2 . 30 20 10 0 10 20 30 t-SNE 1 40 20 0 20 40 t-SNE 2 T rain T est t-SNE with L 2 metric on the r egistration displacements Centerline encoding T emplate geometry Figure 6. Dataset of aortic geometries. Left: template geometry , and example of cen terline enco ding with inscrib ed radii. Right: t-SNE clustering of the dataset based on the similarity matrix of the LDDMM registrations from the template geometry . 4.2. Multi-grid ResNet-LDDMM registration of aortic geometries. W e apply the m ulti-grid ResNet-LDDMM registration (section 2.3 , see also [ 44 ]), to compute the registration maps ϕ i , i = 1 , . . . , n S − 1 from the c hosen template geometry Ω 0 to eac h differen t target geometry Ω i . W e refer the reader to [ 44 ] for a detailed analysis of the registration accuracy in terms of the Chamfer distance, as w ell as for a detailed study on the robustness of the approach with resp ect to the choice of the template geometry . The multi-grid ResNet-LDDMM registration is performed on tetrahedral surface meshes, considering three levels of refinemen t during the training, as describ ed in [ 44 ]. T o p erform the numerical sim ulations, the surface registration maps are then extended in the bulk domain, and used to push forw ard the hexahedral mesh defined on the template to each target domain (Figure 9 left). This is p ossible b ecause the mappings ϕ 0 ,i 1 are defined contin uously on the ambien t space R 3 and are indep enden t of the sp ecific mesh discretization of the domains, except through the optimization problem solved during the LDDMM registration. How ever, giv en the high v ariability of the shap e surfaces, a priori it is not guaran teed that the extension based on interpolation preserves the quality and the consistency of the original mesh (e.g., CONDITIONAL LDDMM STOCHASTIC INTERPOLANTS 11 leading to degenerating elemen ts). The approach emplo y ed for the extension of the registration will be describ ed in section 4.4 . T o illustrate the v ariability of the dataset, Figure 7 sho ws the distributions of relev ant quantities of interest, suc h as the minim um and maximum deformation gradient, the ratio b et ween maximum and minimum principal stretches, and the energy of the registration map. 10 4 10 2 min x det( 0 i 1 ( x )) 0.0 0.2 0.4 0.6 Density Deter minant minima 10 0 10 1 max x det( 0 i 1 ( x )) 0 1 2 Deter minant maxima 10 1 10 2 10 3 max x max ( 0 i 1 ( x )) min ( 0 i 1 ( x )) 0.0 0.5 1.0 Anisotr opy 0.0 0.2 1 0 || v 0 i t ( 0 i t )|| 2 dt 0.0 2.5 5.0 7.5 Ener gy Figure 7. Distributions of selected quantities of interest of the registered dataset o ver the 1239 training and test geometries. 4.3. Computational hemo dynamics. Let us denote by Ω ⊂ R 3 a computational domain whose bound- ary is defined b y a generic shap e of the considered dataset. The boundary of Ω is decomposed in to a v essel w all ( Γ wall ), an inlet b oundary ( Γ in ), and four outlet b oundaries (BCA, LCCA, LSA, T A, represen ted with Γ i , i ∈ { 1 , 2 , 3 , 4 } ), as shown in Figure 8 . :B CA :LCCA :LS A :T A 0.0 0.2 0.4 0.6 0.8 1.0 t [ s ] 0.0 0.2 0.4 0.6 0.8 1.0 a ( t ) Figure 8. Left: sketc h of the computational domain for the computational blo o d flow mo del ( 18 )-( 19 ). Righ t: Time–dep endent inflo w profile. In Ω , w e mo del the blo od as an incompressible, Newtonian flow, whose dynamics is described b y the v elo cit y u : Ω → R 3 and the pressure fields p : Ω → R , solution to the incompressible Navier–Stok es equations ( ρ∂ t u + ρ u · ∇ u − µ ∆ u + ∇ p = 0 , in Ω , ∇ · u = 0 , in Ω , (18) 12 KA TZ, R OMOR, ZHU, AND CAIAZZO where ρ = 1 . 06 × 10 3 kg / m 3 stands for the blo od density , and µ = 3 . 5 × 10 − 3 s · P a is the dynamic viscosit y . Equations ( 18 ) are completed by Dirichlet b oundary conditions on the inlet, no-slip b oundary conditions on the vessel wall, and by a so-called lump ed-parameter mo del on the four outlet b oundaries:      u = u in = a ( t ) β ( x ) , on Γ in , u = 0 , on Γ wall ,  − pI + µ ( ∇ u + ∇ u T )  n = − P i n , on Γ i , ∀ i. (19) The v elo cit y prescrib ed at inlet is mo deled as a time-dependent function a ( t ) m ultiplied b y a parabolic profile β ( x ) . The outlet pressures P i , i = 1 , . . . , 4 , are computed as a function of the b oundary flow rates Q i := R Γ i u · n , i ∈ { 1 , 2 , 3 , 4 } , via a three-elemen ts (RCR) Windkessel mo del [ 51 ]:    C d,i dπ i dt + π i R d,i = Q i , on Γ i ∀ i, P i = R p,i Q i + π i , on Γ i ∀ i, (20) dep ending on an auxiliary distal pressure π i , a proximal resistance R p,i , a distal resistance R d , and a capacitance C d . These parameters mo del the effect of the neglected cardiov ascular system, in terms of the resistance of the arteries in the vicinity of the outlet b oundary , the do wnstream resistance, and the compliance of the arterial system. The choice of the Windkessel parameters strongly influences the resulting hemo dynamics and their calibration typically dep ends on the flo w regime under consideration, on sub ject-sp ecific anatomical features, and on the av ailable clinical data. In this w ork, w e adopt a calibration strategy describ ed in [ 44 ], which is based on prescribing a target flow split among the outlets, and enables the representation of b oth physiological and pathological scenarios in a consistent wa y . W e refer to [ 44 ] for the details ab out the pro cedure. The system of equations ( 18 )-( 19 ) is discretized in space with discontin uous mixed-order DG2-DG1 finite elements and in time using a high-order dual splitting scheme (see [ 19 , 23 , 31 , 30 , 32 , 20 , 22 , 21 ] for details on the numerical metho d and its implementation). F or mo deling the effect of turbulence, a σ -model is used [ 28 ]. The mo del is implemented using ExaDG [ 5 ], which is based on the op en-source library deal-II [ 6 ] and relies on a matrix-free implementation on hierarchical hexahedral meshes. The hierarc hical structure of meshes is necessary for applying the geometric multigrid preconditioners whic h are critical for the efficiency of the metho ds, and must b e preserved when registering the template meshes on to different shap es. This issue will b e discussed in the following section 4.4 . 4.4. Preserv ation of hexahedral mesh hierarchies structure up on registration. The registration maps b et ween aortic shap es (section 4.2 ) are computed to match only the surface meshes, since considering a full three–dimensional registration would considerably increase the computational burden. A more efficient alternative consists then in prop erly defining extensions of the surface m aps into the in terior domain and transp orting the mesh accordingly , using, e.g., radial basis function interpolation, or solving an auxiliary problem. Additionally , smo othing steps migh t b e required in the presence of large deformations to improv e mesh quality , av oid the degeneration of mesh cells, and preserve mesh regularit y prop erties which might b e required or affect the stability of the n umerical metho d. The considered n umerical solv er, implemented in ExaDG [ 5 ], employs geometric multigrid precondition- ers on nested hexahedral meshes. Hence, the hierarchical structure needs to b e preserved when extending the surface displacement map to the mesh vertices in the volume domain. The extension algorithm considered in this w ork consists of an iterativ e algorithm solving a non- homogeneous linear elastic problem and a smo othing step. Namely , let M 0 and M i b e tw o different meshes, the template mesh and the target mesh generated with the conditional LDDMM stochastic in terp olan t algorithm 2 resp ectiv ely , and let us denote with Ψ s : = E [ I t [ c ]] | t =1 the displacement of the surface mesh no des from M 0 to M i . T o handle large deformations, the surface displacement Ψ s is split in N steps iterations. Setting Ψ 0 = I d , for n = 1 , . . . , N steps , we solve a fictitious linear elastic problem for the intermediate displacement field.        ∇ ·  E ( x ) 2(1 + ν )  ∇ Ψ n + ( ∇ Ψ n ) T  + E ( x ) ν (1 + ν )(1 − 2 ν ) ( ∇ · Ψ n ) I  = 0 , in Ψ n − 1 (Ω 0 ) , Ψ n | ∂ Ω = e − n N steps − 1 e − 1 − 1 Ψ s , on ∂ Ψ n − 1 (Ω 0 ) , (21) where the non-homogeneous Y oung mo dulus E ( x ) is defined as the square of the finite element cell aspect ratio and the Poisson’s ratio ν = 0 . 3 is fixed. A t the end of each step, inv erted cells are detected and iterativ ely corrected adjusting the v alue of the displacement field Ψ n at the no des. CONDITIONAL LDDMM STOCHASTIC INTERPOLANTS 13 The surface solution obtained after N steps is then pro cessed by a smo othing algorithm that defines in ternal displacement limiting the num b er of cells with bad asp ect ratio (details of the smo othing step are provided in app endix B ). Displacement e xtended to the interior [m] T emplate's nested meshes T ar get's generated nested meshes Figure 9. T ransport of nested hexahedral mesh hierarchies needed for geometric multi- grid preconditioners. Figure 9 sho ws an example of the nested hexahedral mesh hierarch y transp orted from M 0 to the target generated mesh M i , asso ciated to a particular realization c i ∈ R c of the conditioning random v ariable (section 5.1 ). R emark 4.1 (Generation of the template mesh) . The describ ed mesh extension pro cedure has b een used to define also the nested hexahedral mesh M 0 on the template geometry , transporting it from another mesh [ 9 ], not contained in our original dataset. 5. Applica tion to shape uncer t ainty quantifica tion in cardiov ascular models Studying the sensitivit y of quan tities of interest with resp ect to domain uncertain ties requires the generation of random domain perturbations. In the case of anatomical shap es, which are not describ ed b y a small num b er of parameters, this step needs a suitable and robust pro cedure that preserves the ph ysiological features of the anatom y and fits in to the considered numerical scheme in terms, e.g., of mesh top ology and mesh quality . 5.1. Numerical methods for geometric uncertain ty quantification. W e apply the conditional LDDMM sto c hastic in terp olan t metho d to generate domain p erturbations that will b e used to study the impact of geometrical uncertaint y . Namely , shap es are generated through deformations Ψ : ˆ Ω → Ω of a reference domain ˆ Ω ⊂ R 3 , obtained from LDDMM flows (as describ ed in section 3 ), by v arying a conditioning v ariable c : c ∈ R c 7→ { b θ t [ c ] } t ∈ [0 , 1] 7→ Ψ c : = E [ I t | t =1 [ c ]] , (22) where I t is obtained from the SDE ( dI t [ c ] = b θ t [ c ]( I t ) dt + σ t dW Q t , I t =0 = Id . (23) dep ending on the drift b θ [ c ] . F or ease of notation, w e also denote by Ψ c the LDDMM displacement field generated by the conditional sto c hastic interpolant ( 22 )-( 23 ), extended in the interior domain as describ ed in subsection 4.4 . R emark 5.1 (Regularity of the p erturbation) . The spatial regularity of the domain p erturbation x 7→ Ψ c ( x ) and the parameter-to-perturbation map c 7→ Ψ c are gov erned b y the regularity of the vector field b with resp ect to x and the parameter θ , resp ectiv ely . When b θ and the map c 7→ { b θ t [ c ] } t ∈ [0 , 1] are discretized via neural ne t works, their regularit y is inherited from the choice of activ ation functions. Sp ecifically , while ReLU activ ations provide Lipsc hitz contin uity (often asso ciated with H s for s < 3 / 2 in this context), tanh or sigmoid activ ations yield analytic (Gevrey-1) regularity . F urthermore, the use of Gevrey bump functions allows for the recov ery of Gevrey- σ regularity with σ > 1 . Several studies ha ve b een dedicated to analyzing the regularity of parameter-to-solution maps: analytic regularity (see, e.g., [ 26 , 10 ]), ( b, ϵ, p ) -holomorphic maps [ 15 ], in particular for the stationary incompressible Na vier-Stokes equations [ 12 ], and the Gevray regularity [ 11 , 25 , 14 ]. 14 KA TZ, R OMOR, ZHU, AND CAIAZZO The v anilla Monte Carlo metho d will b e used to compute estimators of the exp ectation of a selected quan tity of interest Q ( u (Ω) , p (Ω)) , defined as a function of the Navier–Stok es equations on a p erturb ed domain Ω , i.e., E [ Q ] ≈ 1 M M X i =1 Q ( u (Ω i ) , p (Ω i )) , where M is the num b er of samples, Ω i := Ψ c i ( ˆ Ω) , and c i ∼ Λ P ( C ) , the probability distribution of the conditioning v ariable C , as defined in subsection 3.2 . R emark 5.2 . Alternative n umerical strategies to compute statistics of quantities of interest include v ari- an ts of the Mon te Carlo metho ds (e.g., quasi-Mon te Carlo, multilev el Monte Carlo), p olynomial chaos expansions, realized through sto chastic collo cation [ 40 ] (tensor quadrature rules, sparse quadrature rules, etc.), and sto chastic Galerkin metho ds [ 7 ]. R emark 5.3 . F or the sak e of completeness, it is worth mentioning that in the context of uncertain t y quan- tification, the most common numerical metho d to induce random domain p erturbations is the truncated Karh unen-Lo èv e expansion of Gaussian random fields (see, e.g., [ 8 ]) defined on the reference domain ˆ Ω with a prescrib ed cov ariance kernel (e.g., RBF kernel, Matérn kernel, etc.): Ψ( ˆ x ; c ) = ˆ x + K X i =1 p λ i ψ i ( ˆ x ) c i , c = ( c 1 , . . . , c K ) ∼ N (0 , I ) , (24) where { ( λ i , ψ i ) } K i =1 are the first K eigenpairs of the cov ariance operator associated with the chosen k ernel. In this case, the regularity of the bases { ψ i } K i =1 of the expansion determines the regularity of the p erturbations. Other techniques include the randomization of parametrized geometrical descriptions suc h as free-form deformation [ 38 ], radial basis functions, F ourier and wa v elet representations, or ad-ho c parametrized deformations [ 9 ]. 5.2. Domain p erturbations for aortic geometries. In the case of aortic shap es under study , the conditioning v ariable c is defined as a low-dimensional representation of the aortic shap es, comp osed of the co ordinates of 94 centerline control p oin ts and the asso ciated inscrib ed radii, i.e., c = ( p , r ) ∈ R 94 × 4 . The time-dep enden t conditioned v ector field b θ t [ c ] : R 3 → R 3 is approximated using a graph neural net work made of consecutive MeshGraphNet lay ers [ 42 ]. The training pro cedure follows algorithm 1 and is describ ed in more detail in app endix A . Notice that, during training, the original LDDMM sto chastic in terp olan t loss (equation ( 15 )) is approximated with mini-batc hes of size n b randomly sampled from the 1209 training shap es: L batch ( θ ) = 1 n b n b X i =1 E t,ϵ  ∥ b θ t ( I t [ c i ] , c i ) − v 0 i t ( ϕ 0 i t ( m 0 )) − σ t ϵ ∥ 2 2  , where c i = ( p i , r i ) is the conditioning v ariable asso ciated with the i -th training shap e. F or eac h batc h, the functions are ev aluated on the same set of p oin ts m 0 (corresp onding to the no des of the hexahedral template mesh), and then mapp ed onto the deformed geometries via the corresp onding registration map. Since our aim is to generate p erturbations of given aortic shap es, at the end of training, w e fine-tune the mo del on the 30 test shap es for a few ep ochs. Once the mo del is trained and fine-tuned, we can generate p erturbations of a giv en test shap e asso ciated with c by simulating the SDE ( 23 ) and setting the generated p erturbation as Ψ = E [ I t | t =1 [ c ]] . In particular, Figure 10 , left, shows an example of shap es generated from a conv ex latent interpolation. In this case, w e ha v e c hosen 4 test geometries, with corresp onding latent co ordinates c 1 , . . . , c 4 , and considered the conditioning v ariable asso ciated with the conv ex combinations c interp = 4 X i =1 α i c i , with α i ∈ [0 , 1] , 4 X i =1 α i = 1 . W e can observe that the interpolated shap es represent smo oth transitions b et ween the four chosen ge- ometries, and that the intermediate conditioning v ariables yield realistic aortic shap es, which were not con tained in the original dataset. A dditionally , it is p ossible to increase the v ariability of the generated shap es by v arying uniformly the radii of inscrib ed spheres and p erturbing the control p oin ts’ co ordinates with a Gaussian random field with RBF k ernel. Figure 11 shows some examples of shap es generated with this approach. This strategy will b e used in section 6.2 for the shap e uncertaint y quantification study . CONDITIONAL LDDMM STOCHASTIC INTERPOLANTS 15 Displacement magnitude [m] Centerline inscribed radius [m] 4.8e-05 0.001 0.002 0.003 4.2e-03 1.4e-01 0.4 0.6 0.8 9.9e-01 Conve x latent interpolation Induced parametric defor mations (1, 0 , 0 , 0 ) (2 / 3 , 1 / 3 , 0 , 0 ) (1/ 3 , 2 / 3 , 0 , 0 ) (2 / 3 , 0 , 1 / 3 , 0 ) (2 / 5 , 1 / 5 , 1 / 5 , 1 / 5 ) (1/ 5 , 2 / 5 , 1 / 5 , 1 / 5 ) (0, 2 / 3 , 0 , 1 / 3 ) (1/ 3 , 0 , 2 / 3 , 0 ) (1/ 5 , 1 / 5 , 2 / 5 , 1 / 5 ) (1/ 5 , 1 / 5 , 1 / 5 , 2 / 5 ) (0, 1 / 3 , 0 , 2 / 3 ) (0, 0 , 1 , 0 ) (0, 0 , 2 / 3 , 1 / 3 ) (0, 0 , 1 / 3 , 2 / 3 ) (0, 0 , 0 , 1 ) (0, 1 , 0 , 0 ) Figure 10. Left: Example of shape generated from a con v ex in terp olation on the laten t space of four test geometries. Righ t: Limitations of the metho d with resp ect to lo cal parametric deformations. W e imp ose a v ariation of the inscrib ed radii b y a factor of α r ∈ { 0 . 2 , 0 . 7 , 1 . 3 , 1 . 8 } , limited to the cen terline points close to a selected lo cation (red circles). This results in p erturb ed shap es with global radius v ariations. The prop osed sto c hastic interpolant LDDMM allows one to generate domain p erturbations of aortic shap es based on their latent represen tations c . On the other hand, lo cal deformations, such as enlarge- men ts or contractions at sp ecific lo cations, which are relev ant for certain applications (e.g., presence of aneurysms or stenoses), might b e difficult to achiev e. T o illustrate this issue, we selected a test geometry and generated p erturbed shap es by v arying lo cally the latent co ordinates corresp onding to the inscrib ed radii close to a giv en lo cation by a factor of α r ∈ { 0 . 2 , 0 . 7 , 1 . 3 , 1 . 8 } . The results of this exp erimen t (Figure 10 , right) sho w that the radii were globally increased or decreased with resp ect to the original radii. This b eha vior can b e explained by the fact that the training dataset do es not con tain shap es with lo cal v ariations of the inscrib ed radii, and thus the mo del is not able to learn suc h lo cal deformations. P ossible strategies to ov ercome this issue include augmenting the training dataset with shap es with lo cal v ariations of the inscrib ed radii. R emark 5.4 . No metric constraints were imp osed during the training of the generative mo del. Hence, c hanging the radius of the inscrib ed spheres by a sp ecific factor do es not guarantee that the generated shap e will reflect this change accurately . This can b e seen in Figure 11 , which sho ws that the v alue of the inscrib ed spheres’ radii and the asso ciated displacemen t magnitude are of the same order of magnitude but not exactly the same. Adjustmen ts and enhancements of the metho d to impro ve the accuracy of these parametric deformations will b e the sub ject of future research. 6. Numerical resul ts 6.1. Sim ulation setup. Numerical sim ulations for uncertaint y quantification hav e b een run based on a test dataset containing 30 geometries. F or eac h of these, we generated three additional sets of n s = 10 p erturbed domains by imp osing a radius v ariation by a factor α r ∈ { 0 . 7 , 1 , 1 . 3 } . Figure 12 shows the whole test dataset, while Figure 13 sho ws a selection from the generated batches of geometries after 16 KA TZ, R OMOR, ZHU, AND CAIAZZO Inscribed spher es' radii [m] Inscribed spher es' radii [m] Inscribed spher es' radii [m] Displacement magnitude fr om one single fix ed test geometry [m] 1.1e-04 0.005 9.7e-03 0.031 0.005 0.004 0.025 0.002 0.017 Figure 11. Shap e generation with conditional LDDMM sto c hastic interpolant. The three different rows corresp ond to a different factor applied to the reference inscrib ed radii: 0 . 7 , 1 . 0 , and 1 . 3 ; in each row, w e show the four generated samples, v arying only the Gaussian random field employ ed to p erturb the con trol p oin ts co ordinates, rep orted on the left. The displacement magnitude of the asso ciated registration maps from the template to the generated test shap e is also shown. the p erturbation. The c hoice α r ∈ { 0 . 7 , 1 . 3 } is made to test the generative model and mesh extension algorithms with high deformations: they may not produce physiologically realistic geometries, but they are useful to explore the b ehavior of the mo del ov er a wide range of v ariations. In all considered cases, the generativ e mo del yielded v alid surfaces. Ho w ever, for certain geometries, particularly when reducing the radius of already relatively narrow ones, it was not p ossible to generate a v alid volume mesh or to guaran tee mesh quality sufficien t for stable simulations. These geometries (46 in total out of the full set of 900) were remov ed from the corresp onding batches and not considered further in the study . The num ber of v alid and excluded cases in eac h batch is detailed in T able 1 . T able 1. V alid geometries for each batch (num bered from 0 to 29 ). Eac h cell contains the num b er of full simulations av ailable for eac h batc h and for the corresponding v alue of α r . F or the batches not listed here, all 10 p erturbations, for all three v alues of α , were a v ailable. The last ro w provides the total num b er of cases for each α r . Incomplete batches α r 0 2 7 16 17 18 21 23 24 25 total 0 . 7 6 4 0 1 10 8 5 10 8 10 262 1 . 0 10 10 10 10 9 10 10 5 9 9 292 1 . 3 10 10 9 10 10 10 10 10 6 10 295 In each geometry , the Navier–Stok es equations are solved n umerically as describ ed in Section 4.3 . The same p erio dic inflow profile, with a p eriod of 1 second, is used for all simulations, rescaled according to the area of the inlet surface in order to provide the same volumetric p eak of 400 cm 3 s − 1 (see Figure 8 ). No-slip b oundary conditions on the vessel walls are enforced weakly due to the discontin uous Galerkin CONDITIONAL LDDMM STOCHASTIC INTERPOLANTS 17 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 Figure 12. The 30 test geometries used to generate the different batches. metho d employ ed by the ExaDG solv er. The parameters defining the b oundary conditions are calibrated only once for each batch, dep ending on the base test geometry . The calibration is based on mo deling the flo w split b et ween the four outlets as a function of the areas, as describ ed in [ 44 ]. Sim ulations are run for three cardiac cycles in order to reach a p eriodic state, and the final cycle is considered for ev aluating the quantities of interest. T o visualize the v ariability of the p erturbed domains, Figure 14 shows the sample standard deviation of the vertex co ordinates for batches 0, 10, 20, and 29. Figure 15 gives an o verview of the v ariabilit y within eac h batc h, our source of geometric uncertain t y . Figures 16 and 17 show the v ariability of velocity and pressure fields for one selected test geometry (the base geometry for batch 1) with α r = 1 . 0 at different 18 KA TZ, R OMOR, ZHU, AND CAIAZZO Figure 13. Examples from batches 0, 10, 20 and 29 (from left to right). T op: p er- turbations with radius reduced b y 30% (factor α r = 0 . 7 ). Middle: p erturbations with similar radius (factor α r = 1 ). Bottom: p erturbations with radius increased by 30% (factor α r = 1 . 3 ). time steps. Most of the uncertain t y in the v elo cit y field is concentrated during diastole, in the region of the desce nding aorta, as observed also in [ 8 ], due to the developmen t of secondary flo w patterns. Figure 14. Meshes generated from av erage co ordinates for the batches 0, 10, 20, and 29. The surface color displays the sample standard deviation. 6.2. Shap e uncertaint y quan tification on selected biomarkers. W e will compare the impact of geometric v ariation on a n umber of quantities of interest: • Ov erall maximum velocity magnitude throughout, and mean magnitude of the velocity compo- nen ts across cross–sections orthogonal to the aorta’s centerline. • Outflo w through the descending aorta ( Q desc ). • Ov erall mean pressure ¯ P and pressure difference ∆ P b et ween the descending aorta outlet and the inlet. CONDITIONAL LDDMM STOCHASTIC INTERPOLANTS 19 0 5 10 15 20 25 29 1.0 1.5 2.0 2.5 3.0 d c [ m m ] 0 5 10 15 20 25 29 1 2 3 4 5 [ m m ] Figure 15. Shap e v ariability within each batch. T op: F or each batch index (0 to 29), the bars show mean and sample standard deviation of all pairwise Chamfer distances of the shap es within the batch. Bottom: F or each batch index (0 to 29), the bars show mean and sample standard deviation of vertex co ordinate standard deviations. 4.6e-04 3.6e00 7.2e-05 5.7e-01 5.3e-06 3.6e-01 9.2e-07 1.7e-01 V elocity mean magnitude [m/s] 1.7e-03 2.5e00 2.7e-04 5.2e-01 2.6e-05 2.3e-01 6.3e-06 1.1e-01 V elocity std magnitude [m/s] t=0.2s t=0.4s t=0.6s t=0.8s t=0.2s t=0.4s t=0.6s t=0.8s Figure 16. V elo cit y mean and standard deviation ov er 10 sampled perturbations of one selected test geometry (the base geometry for batch 1) with α r = 1 . 0 , and for t ∈ { 0 . 2 s , 0 . 4 s , 0 . 6 s , 0 . 8 s } . • W all shear stress (WSS), defined as τ w ( x ) = ν ∂ ∂ n  u − ( u · n ) n  ( x ) , 20 KA TZ, R OMOR, ZHU, AND CAIAZZO 1.40e04 2.16e04 1.36e04 1.42e04 1.18e04 1.19e04 9.61e03 9.72e03 P r essur e mean [P a] 2.2e02 2.0e03 2.3e01 1.8e02 2.1e01 5.0e01 2.1e01 2.9e01 P r essur e std [P a] t=0.2s t=0.4s t=0.6s t=0.8s t=0.2s t=0.4s t=0.6s t=0.8s Figure 17. Pressure mean and standard deviation ov er 10 sampled p erturbations of one selected test geometry (the base geometry for batch 1), with α r = 1 . 0 , and for t ∈ { 0 . 2 s , 0 . 4 s , 0 . 6 s , 0 . 8 s } where n = n ( x ) is the outer normal at a b oundary p oin t x ∈ Γ , so that homogeneous shear flo w o v er a plane corresp onds to w all shear stress in the direction opposite to the flow. The WSS will b e analyzed mostly in terms of its mean magnitude ov er Γ wall or ov er a reference patch Γ ref ⊂ Γ wall on the underside of the descending aorta. • Oscillatory shear index (OSI), defined as OSI( x ) = 1 2   1 −    R 3 2 τ w ( x )    d t R 3 2 | τ w ( x ) | d t   , whic h measures how strongly τ w oscillates ov er a given length of time during the simulated cardiac cycle, at a given b oundary point x . Note that, unlike the other quan tities in tro duced here, the OSI is not time-dep endent. W e will, as with τ w itself, consider the mean OSI ov er Γ wall or Γ ref . • Secondary flo w degree (SFD), defined, for a cross section S , as SFD( S ) = R S | u − ( u · n ) n | d x R S | u · n | d x , whic h measures the ratio of tangential to normal flow ov er a cross–section S . Low SFD cor- resp onds to strongly directed flow through S , whereas high SFD indicates a large amount of in–plane flow due to eddies or (in the aortic arch) lateral outlets. • Normalized flo w displacement (NFD), defined, for a cross-section S , as NFD( S ) = | x n ( S ) − x g ( S ) | r H ( S ) . (25) The NFD is a measure of flow e c c entricity , i.e., of the normaliz ed distance b et ween the moment of the flo w through S , x n ( S ) = R S | ( u · n ) | x d x R S | u · n | d x , and its geometric centroid, x g ( S ) = R S x d x R S d x . In ( 25 ), the normalizing factor r H ( S ) is the h ydraulic radius, which is t ypically defined as the ratio of area to (wetted) p erimeter; we use the more stable estimate r H ( S ) = 3 4 R S | x − x G ( S ) | d x R S d x , which matches the conv entional definition if S is a disc. CONDITIONAL LDDMM STOCHASTIC INTERPOLANTS 21 Figure 18 shows the cross-sections S used for the computation of the quantities of in terest for one case in batch 20. Figure 18. Left: selected geometry (from batc h 20) showing the cross–sections used for computing av erages of quantities of interests. The outer rings are in-plane circles matc hing an estimate of the outer radius of each cross–section. The inner rings show the estimated hydraulic radius r H . Figures 19 show the results for the pressure drop as a function of the av eraged cross-sectional velocity and of the outflow at the descending aorta for batches 0, 10, 20, and 29. The cases with larger v essel radii (batc hes 0 and 10), where the pressure drop is low er (to achiev e the same v olumetric flo w), sho w less sensitiv e b eha vior with respect to the geometry v ariation. Similarly , within the same batch, the highest impact of the geometry is seen for the set of simulations with reduced radius ( α r = 0 . 7 ). 0.0 0.5 1.0 1.5 2.0 2.5 | u | [ m / s ] 150 100 50 0 P [ m m H g ] 0.0 0.5 1.0 1.5 2.0 2.5 | u | [ m / s ] 150 100 50 0 P [ m m H g ] 0.0 0.5 1.0 1.5 2.0 2.5 | u | [ m / s ] 150 100 50 0 P [ m m H g ] 0.0 0.5 1.0 1.5 2.0 2.5 | u | [ m / s ] 150 100 50 0 P [ m m H g ] 0 5 10 15 Q d e s c [ l / m i n ] 150 100 50 0 P [ m m H g ] 0 5 10 15 Q d e s c [ l / m i n ] 150 100 50 0 P [ m m H g ] 0 5 10 15 Q d e s c [ l / m i n ] 150 100 50 0 P [ m m H g ] 0 5 10 15 Q d e s c [ l / m i n ] 150 100 50 0 P [ m m H g ] Figure 19. T op row: Pressure difference b et w een descending aorta and inlet v ersus mean cross-sectional normal v elo cit y magnitude (a veraged o v er the length of the domain) o ver one cardiac cycle (time go es counterclockwise). Bottom row: Pressure difference b et ween descending aorta and inlet versus outflow at the descending aorta v ersus pressure difference b et w een descending aorta and inlet ov er one cardiac cycle. F rom left to right: batc hes 0, 10, 20, and 29. Black: α r = 1 . 3 (30% radius increase), blue: α r = 1 , red: α r = 0 . 7 (30% radius decrease). Figure 20 analyzes the v ariability of the results at different instants in time. In particular, the fraction of outflo w through the descending aorta at peak systole v aries notably . This b eha vior can be seen as a consequence of the differences in flow splits, which are imp osed through the Windkessel b oundary conditions (Section 4.3 ). The cases with low er radii are also those showing the highest pressure v ariabilit y . Figure 20 illustrates the velocity–pressure relationship for all simulations at different time steps. At peak 22 KA TZ, R OMOR, ZHU, AND CAIAZZO 0.8 0.6 0.4 0.2 0.0 P [ m m H g ] 0.3 0.4 0.5 0.6 Q d e s c [ l / m i n ] 150 100 50 0 P [ m m H g ] 8 10 12 14 16 Q d e s c [ l / m i n ] 60 40 20 0 P [ m m H g ] 4 6 8 10 Q d e s c [ l / m i n ] 4 3 2 1 0 1 P [ m m H g ] 0.5 1.0 1.5 Q d e s c [ l / m i n ] 0.8 0.6 0.4 0.2 0.0 P [ m m H g ] 0.04 0.06 0.08 0.10 0.12 | u | [ m / s ] 150 100 50 0 P [ m m H g ] 0.5 1.0 1.5 2.0 2.5 | u | [ m / s ] 60 40 20 0 P [ m m H g ] 0.5 1.0 1.5 | u | [ m / s ] 4 3 2 1 0 1 P [ m m H g ] 0.05 0.10 0.15 0.20 0.25 0.30 | u | [ m / s ] Figure 20. T op row: Pressure difference b et w een descending aorta and inlet v ersus outflo w at the descending aorta for all simulations. Bottom ro w: Pressure difference b et ween descending aorta and inlet versus mean cross–sectional velocity magnitude (a v- eraged ov er the length of the domain). F rom left to right: times t ∈ { 2 , 2 . 16 , 2 . 3 , 2 . 7 } , corresp onding to the end of diastole, peak flow, end of systole, and middle diastole. Blac k: α r = 1 . 3 (30% radius increase), blue: α r = 1 , red: α r = 0 . 7 (30% radius de- crease). systole ( t = 2 . 16 s ), which corresp onds to large pressure differences, larger pressures and large v ariations are alwa ys asso ciated with narrow er geometries. Figure 21 shows that the correlation b et ween pressure and v elo cit y , as w ell as b et ween pressure and w all shear stress, is clearer when examined in terms of pressure differences rather than absolute pressure, since the latter is strongly separated b etw een batches. 40 60 80 100 120 P [ m m H g ] 2.5 5.0 7.5 10.0 12.5 | w | [ P a ] 40 30 20 10 0 P [ m m H g ] 2.5 5.0 7.5 10.0 12.5 | w | [ P a ] 40 60 80 100 120 P [ m m H g ] 0.2 0.4 0.6 0.8 | u | [ m / s ] 40 30 20 10 0 P [ m m H g ] 0.2 0.4 0.6 0.8 | u | [ m / s ] Figure 21. Mean wall shear stress (top) and mean velocity magnitude (b ottom) versus mean pressure (left) and mean pressure drop (right), for all simulations, a v eraged ov er time. Note the clear batc h separation in the absolute mean pressure compared to the pressure drop. Blac k: α r = 1 . 3 (30% radius increase), blue: α r = 1 , red: α r = 0 . 7 (30% radius decrease). The sensitivity of ov erall secondary flo w (normal vs. tangential velocity) is analyzed in Figure 22 . As for the pressure difference, cases with reduced radius ( α r = 0 . 7 ) are more sensitive to geometrical uncertain ty . Ho wev er, the secondary flow degree (SFD) (Figure 22 , middle and b ottom) is considerably less sensitive to geometric v ariation. The clearest observ ation is that, as exp ected, normal (forward/bac kw ard) flo w dominates when flow rates are high, leading to low SFD, whereas the slow er flow during diastole exhibits SFD closer to 1, corresp onding to decaying turbulence throughout the aorta. Figure 23 shows sample means and standard deviations for velocity , pressure, wall shear stress, and OSI. Additional tables with detailed statistics are provided also in app endix C . In most cases, one can observe a correlation b et w een geometries with smaller v essel radius and larger and more sensitive velocities and pressures. F or one batc h (num ber 12), we obtain rather different results. This geometry is smaller than the others (see Figure 12 ), and this b eha vior could b e a consequence of ha ving scaled the b oundary conditions (inlet flow and outlet lump ed parameters) without taking into CONDITIONAL LDDMM STOCHASTIC INTERPOLANTS 23 0.0 0.5 1.0 1.5 2.0 | u n | [ m / s ] 0.0 0.2 0.4 0.6 0.8 | u ( u n ) n | [ m / s ] 0.0 0.5 1.0 1.5 2.0 | u n | [ m / s ] 0.0 0.2 0.4 0.6 0.8 | u ( u n ) n | [ m / s ] 0.0 0.5 1.0 1.5 2.0 | u n | [ m / s ] 0.0 0.2 0.4 0.6 0.8 | u ( u n ) n | [ m / s ] 0.0 0.5 1.0 1.5 2.0 | u n | [ m / s ] 0.0 0.2 0.4 0.6 0.8 | u ( u n ) n | [ m / s ] 2.0 2.2 2.4 2.6 2.8 3.0 t [ s ] 0.0 0.2 0.4 0.6 0.8 1.0 SFD 2.0 2.2 2.4 2.6 2.8 3.0 t [ s ] 0.0 0.2 0.4 0.6 0.8 1.0 SFD 2.0 2.2 2.4 2.6 2.8 3.0 t [ s ] 0.0 0.2 0.4 0.6 0.8 1.0 SFD 2.0 2.2 2.4 2.6 2.8 3.0 t [ s ] 0.0 0.2 0.4 0.6 0.8 1.0 SFD 2.0 2.2 2.4 2.6 2.8 3.0 t [ s ] 0.0 0.2 0.4 0.6 0.8 1.0 SFD 2.0 2.2 2.4 2.6 2.8 3.0 t [ s ] 0.0 0.2 0.4 0.6 0.8 1.0 SFD 2.0 2.2 2.4 2.6 2.8 3.0 t [ s ] 0.0 0.2 0.4 0.6 0.8 1.0 SFD 2.0 2.2 2.4 2.6 2.8 3.0 t [ s ] 0.0 0.2 0.4 0.6 0.8 1.0 SFD Figure 22. T op ro w: Cross-sectional tangential v elo cit y v ersus normal velocity o v er one cardiac cycle. The time runs counterclockwise around these lo ops. During systole the flo w accelerates and the normal flo w dominates, while, during diastole, the presence of dissipating eddies might result in larger tangential v elo cit y . Middle row: SFD (av eraged o ver the length of the domain) ov er one cardiac cycle. Bottom ro w: SFD (maximum o ver the length of the domain) o ver one cardiac cycle. F rom left to righ t: batc hes 0, 10, 20, and 29. Blac k: α r = 1 . 3 (30% radius increase), blue: α r = 1 , red: α r = 0 . 7 (30% radius decrease ). accoun t the ov erall v olume. A different behavior can b e seen for the OSI. In this case, sensitivities are comparable across batches and radius v ariation. Considering the geometry v ariability of the dataset (see, e.g., the geometries in Figures 13 ), the results suggest that, although the w all shear stress sensitivity v aries across the batc hes, its directional v ariabilit y and the related sensitivity are comparable. Detailed statistics on the biomarkers are provided in App endix C . 6.3. Sample size study. The n umerical study in Section 6.2 prioritized the v ariability of the geometry within the test dataset, with the limitation of considering only a limited sample (a maximum of n s = 10 sim ulations for each individual batch). In this section, we study the influence of considering, for a single batc h (batch 3, with α r = 1 ), a larger sample with an additional n ′ s = 100 geometries. W e compare the resulting statistics with those obtained with the smaller sample size in terms of distributions of quantities of interest, as well as in terms of the W asserstein–1 distance defined, for tw o one–dimensional real–v alued random v ariables X 1 and X 2 , as W 1 ( X 1 , X 2 ) := Z R | F 1 ( x ) − F 2 ( x ) | d x, (26) dep ending on the resp ectiv e cumulativ e distribution functions F 1 and F 2 . Note that, by a change of v ariables, one can write W 1 ( cX 1 , cX 2 ) = c W 1 ( X 1 , X 2 ) for any c > 0 . If X i are concentrated at x i (i.e., F i are step functions), one obtains W 1 ( X 1 , X 2 ) = | x 1 − x 2 | . Hence, if the random v ariables X i represen t tw o distributions of a physical field in a given unit, W 1 ( X 1 , X 2 ) is interpretable as a quantit y in the same unit. Figure 24 shows the distributions of mean pressure and pressure differences for the t wo samples, sho wing that the statistics obtained in b oth cases are comparable. Quantifying this b y the W asserstein–1 distance (b ottom ro w) do es not rev eal any surprises — the distributions are clearly very close compared to the ov erall scale of the pressure. Similar conclusions can b e drawn for the mean and maximum velocity magnitude (along the centerline and ov er the whole domain), Figure 25 , although here some elev ated sensitivit y at p eak flow can b e observed. Figures 26 compare pressure and velocity quan tities of in terest o ver a full cycle. Also in this case, the samples of different sizes b eha ve very similarly . A similar agreement b etw een the tw o cases w as also observ ed for the remaining quantities of in terest. 24 KA TZ, R OMOR, ZHU, AND CAIAZZO 0 5 10 15 20 25 29 2 4 6 8 10 m a x ( | u | ) [ m / s ] 0 5 10 15 20 25 29 30 20 10 0 P [ m m H g ] 0 5 10 15 20 25 29 2.5 5.0 7.5 10.0 12.5 | w | [ P a ] 0 5 10 15 20 25 29 0.1 0.2 0.3 OSI Figure 23. Sample means and standard deviations for differen t quantities o ver the whole dataset of geometries, divided according to the corresp onding batch (from 0 to 29, in the x -axis). First row: maxim um velocity magnitude at p eak systole ( t = 2 . 16 s); second ro w: largest pressure difference betw een inlet and descending aorta o ver the cardiac cycle; third row: mean wall shear stress magnitude ov er the whole blo od v essel w all and o v er the cardiac cycle; F ourth ro w: mean oscillatory shear index o v er a reference patc h on the underside of the descending aorta. Blac k: α r = 1 . 3 (30% radius increase), blue: α r = 1 , red: α r = 0 . 7 (30% radius decrease). Analyzing the results at selected time steps (Figure 27 ), one can see that, as expected, the larger sample allows us to capture a broader range of physiological conditions. CONDITIONAL LDDMM STOCHASTIC INTERPOLANTS 25 2.0 2.2 2.4 2.6 2.8 3.0 t [ s ] 40 60 80 P [ m m H g ] 2.0 2.2 2.4 2.6 2.8 3.0 t [ s ] 20 10 0 10 P [ m m H g ] 2.0 2.2 2.4 2.6 2.8 3.0 t [ s ] 40 60 80 P ± 1 [ m m H g ] 2.0 2.2 2.4 2.6 2.8 3.0 t [ s ] 20 10 0 10 P ± 1 [ m m H g ] Figure 24. T op row: Distribution (mean and standard deviation) of mean pressure in the aorta (left) and pressure difference b et ween descending aorta and inlet (righ t) o ver time. Blac k: larger batc h of 100 simulations; red : original smaller batch. Bottom row: Join t mean of the separate distributions and W asserstein–1 distance b et ween them. 2.0 2.2 2.4 2.6 2.8 3.0 t [ s ] 0.0 0.2 0.4 0.6 0.8 1.0 | u | [ m / s ] 2.0 2.2 2.4 2.6 2.8 3.0 t [ s ] 0.0 0.5 1.0 1.5 2.0 2.5 m a x ( | u | ) [ m / s ] 2.0 2.2 2.4 2.6 2.8 3.0 t [ s ] 0.0 0.2 0.4 0.6 0.8 | u | ± 1 [ m / s ] 2.0 2.2 2.4 2.6 2.8 3.0 t [ s ] 0.0 0.5 1.0 1.5 2.0 2.5 m a x ( | u | ) ± 1 [ m / s ] Figure 25. T op ro w: Distribution (mean and standard deviation) of mean cross– sectional velocity magnitude (left) and ov erall maximum v elo cit y magnitude (right) o v er time. Blac k: larger batc h of 100 simulations; red : original smaller batch. Bottom row: Join t mean of the separate distributions and W asserstein–1 distance b et ween them. 26 KA TZ, R OMOR, ZHU, AND CAIAZZO 0.0 0.2 0.4 0.6 0.8 1.0 | u | [ m / s ] 40 60 80 P [ m m H g ] 0.0 0.2 0.4 0.6 0.8 1.0 | u | [ m / s ] 20 10 0 10 P [ m m H g ] 0.0 2.5 5.0 7.5 10.0 12.5 Q d e s c [ l / m i n ] 40 60 80 P [ m m H g ] 0.0 2.5 5.0 7.5 10.0 12.5 Q d e s c [ l / m i n ] 20 10 0 10 P [ m m H g ] Figure 26. Left: Mean pressure in the aorta versus mean cross–sectional velocity mag- nitude, a veraged ov er the length of the domain (top) and outflo w through the descending aorta (b ottom), ov er one cardiac cycle. Righ t: Pressure difference b et ween descending aorta and inlet v ersus mean cross–sectional v elocity magnitude, a veraged ov er the length of the domain (top) and outflo w through the descending aorta (bottom), ov er one cardiac cycle. Black: larger batch of 100 simulations; red: original smaller batch. 0.10 0.05 0.00 0.05 0.10 P [ m m H g ] 0.0350 0.0375 0.0400 0.0425 0.0450 0.0475 | u | [ m / s ] 16 14 12 10 8 6 P [ m m H g ] 0.85 0.90 0.95 1.00 1.05 | u | [ m / s ] 4 5 6 7 P [ m m H g ] 0.50 0.55 0.60 0.65 | u | [ m / s ] 0.20 0.25 0.30 0.35 P [ m m H g ] 0.060 0.065 0.070 0.075 | u | [ m / s ] Figure 27. Pressure difference b et ween descending aorta and inlet versus mean cross– sectional velocity magnitude (av eraged ov er the length of the domain) for all simulations at times t ∈ { 2 , 2 . 16 , 2 . 3 , 2 . 7 } , corresp onding to the end of diastole (top-left), p eak flow (top-righ t), end of systole (b ottom-left), and middle diastole (b ottom-righ t). Black: larger batch of 100 simulations; red: original smaller batch. Crosses and ellipses indicate sample means and standard deviations. REFERENCES 27 7. Conclusions W e prop ose a new approach for shape generative mo deling based on the combination of sto c hastic in terp olan ts and LDDMM registration. The metho d leverages pairwise diffeomorphic maps computed via multilev el ResNet-LDDMM shap e registration to obtain the time-dependent conditioned drift on in termediate domains required for training the sto c hastic in terpolants, allo wing, in particular, the bridging of distributions defined on different three-dimensional domains. W e v alidated the framework by considering a cohort of synthetic aortic geometries, generating v aria- tions of the av ailable shap es to p erform domain uncertaint y quantification studies for blo o d flow simula- tions, monitoring different quantities of interest (suc h as v elo cit y , pressure differences, wall shear stresses, and OSI). T o this end, the generative model was used to gene rate differen t batches of p erturb ed domains starting from 30 test geometries, considering b oth increases and decreases in v essel radius. Our results sho w that the sensitivit y with resp ect to domain changes is considerably higher for geometries with re- duced lumen size. The pattern differs only in the case of OSI, suggesting that the v ariability of the stress v ector is more affected by v ariations in the surface shap e rather than by the diameter of the vessel. F or this application, we considered a conditioning v ariable related to a latent space representation of the aor- tic shap es, i.e., a given num b er of cen terline p oin ts and the corresp onding radii of the inscrib ed spheres at each p oint. This is a natural choice for this particular setting, and it is also consistent with previous studies that consider the generation of digital shap es using statistical shap e mo dels (SSMs). How ev er, the framew ork can b e adapted to handle different latent descriptions. F or the presented numerical simula- tions, the Navier–Stok es equations are solved via a high-order, discontin uous Galerkin matrix-free solver whic h relies on hexahedral volume meshes. While the choice of the numerical metho d is indep enden t from the generativ e mo del, it is extremely imp ortan t, from the p oin t of view of the application, to b e able to rely on a p erforman t solver, particularly for uncertaint y quantification studies requiring multiple sim ulations. Hence, while prescribing the surface deformation can b e achiev ed via the generative mo del, it is necessary to preserve the internal mesh structure that guarantees solv er performance. T o this end, w e prop ose an iterative metho dology for the deformation of the volume mesh based on the solution of a fictitious elasticity problem, combined with an optimization step for controlling the mesh asp ect ratio. These preliminary results hav e b een obtained b y considering only a rather limited sample of geometries and conditioning v ariables. T o address this issue, w e performed numerical simulations using a larger sample (100 geometries instead of 10) for one particular test geometry , showing go o d agreemen t b et w een the statistics in b oth cases. A limitation of the proposed approach is that it do es not allow for the study of local domain uncertain ty , as shown in Figure 10 (right). This limitation could b e addressed by enriching the training dataset with lo cal deformations. Imposing metric constraints w ould improv e the corresp ondence b et ween the latent v ariables asso ciated with centerline p oin ts and radii, and the actual cen terline p oin ts and radii of the generated shapes. W e remark that, although in the case of the aorta a natural encoding via the cen terline is a v ailable, the metho dology can nonetheless be extended to truly nonparametric shap es, such as general anatomical organs, pro vided they b elong to the same class of bi-Lipsc hitz homeomorphisms. In this regard, the current LDDMM sto c hastic in terp olant can b e extended to generate shapes in different classes of bi-Lipschitz homeomorphisms by considering m ultiple template shap es [ 29 ]. These extensions will b e explored in future work. Curren t and future research, based on the results of this work, includes the application of this approach to differen t clinical con texts (differen t organs and pathological conditions, simulation of devices and treatmen ts) for tasks related not only to uncertaint y quantification, but also to virtual patien t generation for in silico trials, and shape optimization. Current and future research, based on the results of this work, includes the application of this approac h to differen t clinical contexts (different organs and pathological conditions, simulation of devices and treatments) for tasks related not only to uncertain ty quan tification, but also to virtual patient generation for in silico trials, and shap e optimization. A cknowledgements F unded by the Deutsche F orsch ungsgemeinsc haft (DFG, German Research F oundation) under Ger- man y ´ s Excellence Strategy – The Berlin Mathematics Researc h Cen ter M A TH+ (EXC-2046/1, EXC- 2046/2, pro ject ID: 390685689). References [1] Jadie A dams, Riddhish Bhalo dia, and Shireen Elhabian. “ Uncertain-DeepSSM: F rom Images to Probabilistic Shap e Mo dels”. In: Shap e in Me dic al Imaging . Springer In ternational Publishing, 2020, pp. 57–72. doi : 10.1007/978- 3- 030- 61056- 2_5 . 28 REFERENCES [2] Mic hael S Albergo, Nic holas M Boffi, and Eric V anden-Eijnden. “ Sto chastic in terpolants: A unifying framew ork for flows and diffusions”. In: arXiv pr eprint arXiv:2303.08797 (2023). [3] Luigi Ambrosio, Nicola Gigli, and Giusepp e Sav arè. Gr adient Flows in Metric Sp ac es and in the Sp ac e of Pr ob ability Me asur es . Lectures in Mathematics ETH Zürich. Basel: Birkhäuser, 2005. [4] Boulbaba Ben Amor, Sylv ain Arguillère, and Ling Shao. “ ResNet-LDDMM: Adv ancing the LD- DMM F ramework Using Deep Residual Netw orks”. In: IEEE T r ansactions on Pattern Analysis and Machine Intel ligenc e 45.3 (2023), pp. 3707–3720. doi : 10.1109/TPAMI.2022.3174908 . [5] Daniel Arndt et al. “ ExaDG: High-Order Discontin uous Galerkin for the Exa-Scale”. In: Softwar e for Exasc ale Computing - SPPEXA 2016-2019 . Ed. by Hans-Joachim Bungartz et al. Cham: Springer In ternational Publishing, 2020, pp. 189–224. [6] Daniel Arndt et al. “ The deal. I I library , V ersion 9.6”. In: Journal of Numeric al Mathematics 32.4 (2024), pp. 369–380. [7] Iv o Babusk a, Raúl T emp one, and Georgios E. Zouraris. “ Galerkin Finite Element Approximations of Sto c hastic Elliptic Partial Differen tial Equations”. In: SIAM Journal on Numeric al Analysis 42.2 (Jan. 2004), pp. 800–825. doi : 10.1137/s0036142902418680 . [8] Domago j Bošnjak et al. “ Geometric uncertaint y of patient-spe cific blo od vessels and its impact on aortic hemo dynamics: A computational study”. In: Computers in Biolo gy and Me dicine 190 (2025), p. 110017. [9] Domago j Bošnjak et al. “ Syn thAorta: A 3D Mesh Dataset of P arametrized Physiological Healthy A ortas”. In: IEEE T r ansactions on Me dic al Imaging (2025), pp. 1–1. doi : 10 . 1109 / tmi . 2025 . 3599937 . [10] Julio E Castrillon-Candas, F abio Nobile, and Raul F T emp one. “ Analytic regularit y and collo cation appro ximation for elliptic PDEs with random domain deformations”. In: Computers & Mathematics with Applic ations 71.6 (2016), pp. 1173–1197. [11] Alexey Chernov and T ung Le. “ Analytic and Gevrey class regularit y for parametric elliptic eigen- v alue problems and applications”. In: SIAM Journal on Numeric al Analysis 62.4 (2024), pp. 1874– 1900. [12] Alb ert Cohen, Christoph Sch wab, and Jakob Zech. “ Shap e Holomorphy of the Stationary Navier– Stok es Equations”. In: SIAM Journal on Mathematic al Analysis 50.2 (Jan. 2018), pp. 1720–1752. doi : 10.1137/16m1099406 . [13] Giusepp e Da Prato and Jerzy Zab czyk. Sto chastic e quations in infinite dimensions . V ol. 152. Cam- bridge universit y press, 2014. [14] Ana Djurdjev ac et al. Unc ertainty quantific ation for stationary and time-dep endent PDEs subje ct to Gevr ey r e gular r andom domain deformations . 2025. arXiv: 2502.12345 [math.NA] . [15] Jürgen Dölz and F ernando Henríquez. F ul ly discr ete analysis of the Galerkin POD neur al network appr oximation with applic ation to 3D ac oustic wave sc attering . 2025. arXiv: 2502.01859 [math.NA] . [16] Ric hard Duong et al. T ele gr apher’s Gener ative Mo del via Kac Flows . 2025. arXiv: 2506 . 20641 [math.AP] . [17] P aul Dupuis, Ulf Grenander, and Michael I. Miller. “ V ariational problems on flows of diffeomor- phisms for image matching”. In: Quarterly of Applie d Mathematics 56.3 (Sept. 1998), pp. 587–600. doi : 10.1090/qam/1632326 . [18] F rancesco F abbri et al. “ Graph-Con v olutional-Beta-V AE for syn thetic ab dominal aortic aneurysm generation”. In: Me dic al & Biolo gic al Engine ering & Computing (2025). doi : 10 . 1007 / s11517 - 025- 03491- y . [19] Niklas F ehn. “ Robust and Efficien t Discon tinuous Galerkin Metho ds for Incompressible Flo ws”. PhD thesis. T echnisc he Universität München, 2021. [20] Niklas F ehn, W olfgang A. W all, and Martin Kronbic hler. “ Efficiency of high-p erformance discon- tin uous Galerkin spectral element metho ds for under-resolved turbulent incompressible flows”. In: Int. J. Numer. Meth. Fluids 88.1 (2018), pp. 32–54. doi : 10.1002/fld.4511 . [21] Niklas F ehn, W olfgang A. W all, and Martin Kronbic hler. “ On the stability of pro jection methods for the incompressible Navier–Stok es equations based on high-order discontin uous Galerkin discretiza- tions”. In: Journal of Computational Physics 351 (2017), pp. 392–421. doi : 10.1016/j.jcp.2017. 09.031 . [22] Niklas F ehn, W olfgang A. W all, and Martin Kronbic hler. “ Robust and efficien t discontin uous Galerkin metho ds for under-resolv ed turbulent incompressible flows”. In: Journal of Computational Physics 372 (2018), pp. 667–693. doi : 10.1016/j.jcp.2018.06.037 . [23] Niklas F ehn et al. “ Hybrid multigrid metho ds for high-order discontin uous Galerkin discretizations”. In: Journal of Computational Physics 415 (2020), p. 109538. doi : 10.1016/j.jcp.2020.109538 . REFERENCES 29 [24] Leonid Goubergrits et al. “ CT-based analysis of left ven tricular hemo dynamics using statistical shap e mo deling and computational fluid dynamics”. In: F r ontiers in Car diovascular Me dicine 9 (2022), p. 901902. [25] Helm ut Harbrech t, Marc Schmidlin, and Christoph Sch w ab. “ The Gevrey class implicit mapping theorem with application to UQ of semilinear elliptic PDEs”. In: Mathematic al Mo dels and Metho ds in Applie d Scienc es 34.05 (2024), pp. 881–917. [26] F ernando Henríquez and Christoph Sch wab. “ Shap e Holomorphy of the Calderón Pro jector for the Laplacian in R 2 ”. In: Inte gr al Equations and Op er ator The ory 93.4 (July 2021). doi : 10 . 1007 / s00020- 021- 02653- 5 . [27] T ero Karras et al. “ Elucidating the design space of diffusion-based generative models”. In: A dvanc es in neur al information pr o c essing systems 35 (2022), pp. 26565–26577. [28] Sarah Katz et al. “ Impact of turbulence mo deling on the simulation of bloo d flow in aortic coarc- tation”. In: International Journal for Numeric al Metho ds in Biome dic al Engine ering 39.5 (2023), e3695. [29] F anw ei Kong et al. “ SDF4CHD: Generativ e mo deling of cardiac anatomies with congenital heart defects”. In: Me dic al Image Analysis 97 (2024), p. 103293. doi : https : / / doi . org / 10 . 1016 / j . media.2024.103293 . [30] Benjamin Krank et al. “ A high-order semi-explicit discon tinuous Galerkin solv er for 3D incompress- ible flo w with application to DNS and LES of turbulen t c hannel flo w”. In: Journal of Computational Physics 348 (2017), pp. 634–659. doi : 10.1016/j.jcp.2017.07.039 . [31] Martin Kron bic hler and Katharina Kormann. “ F ast matrix-free ev aluation of discon tin uous Galerkin finite element op erators”. In: ACM T r ansactions on Mathematic al Softwar e (TOMS) 45.3 (2019), pp. 1–40. doi : 10.1145/3325864 . [32] Martin Kronbic hler et al. “ A Next-Generation Discontin uous Galerkin Fluid Dynamics Solver with Application to High-Resolution Lung Airflow Sim ulations”. In: Pr o c e e dings of the International Confer enc e for High Performanc e Computing, Networking, Stor age and Analysis . SC ’21. St. Louis, Missouri: Ass ociation for Computing Machinery, 2021. doi : 10.1145/3458817.3476171 . [33] Y aron Lipman et al. “ Flo w matching for generative modeling”. In: arXiv pr eprint (2022). [34] Xingc hao Liu, Chengyue Gong, and Qiang Liu. “ Flow straight and fast: Learning to generate and transfer data with rectified flow”. In: arXiv pr eprint arXiv:2209.03003 (2022). [35] Yixin Liu et al. Sor a: A R eview on Backgr ound, T e chnolo gy, Limitations, and Opp ortunities of L ar ge Vision Mo dels . 2024. arXiv: 2402.17177 [cs.CV] . [36] Ziming Liu et al. GenPhys: F r om Physic al Pr o c esses to Gener ative Mo dels . 2023. arXiv: 2304.02637 [cs.LG] . [37] Gabriel D. Maher et al. “ Geometric uncertaint y in patient-specific cardiov ascular mo deling with con volutional drop out netw orks”. In: Computer Metho ds in Applie d Me chanics and Engine ering 386 (Dec. 2021), p. 114038. doi : 10.1016/j.cma.2021.114038 . [38] Andrea Manzoni, Alfio Quarteroni, and Gianluigi Rozza. “ Shap e optimization for viscous flows by reduced basis metho ds and free-form deformation”. In: International Journal for Numeric al Metho ds in Fluids 70.5 (2012), pp. 646–670. [39] Nik ola j T. Mück e and Benjamin Sanderse. Physics-awar e gener ative mo dels for turbulent fluid flows thr ough ener gy-c onsistent sto chastic interp olants . 2025. arXiv: 2504.05852 [cs.CE] . [40] F abio Nobile, Raúl T emp one, and Clayton G W ebster. “ A sparse grid sto c hastic collo cation method for partial differential equations with random input data”. In: SIAM Journal on Numeric al Analysis 46.5 (2008), pp. 2309–2345. [41] NVIDIA PhysicsNeMo con tributors. “ NVIDIA PhysicsNeMo: An op en-source framework for physics- based deep learning in science and engineering”. In: (2025). doi : 10.5281/zenodo.15588532 . [42] T obias Pfaff et al. “ Learning mesh-based simulation with graph netw orks”. In: International c on- fer enc e on le arning r epr esentations . 2020. [43] A dity a Ramesh et al. Hier ar chic al T ext-Conditional Image Gener ation with CLIP L atents . 2022. arXiv: 2204.06125 [cs.CV] . [44] F rancesco Romor et al. Data assimilation p erforme d with r obust shap e r e gistr ation and gr aph neur al networks: applic ation to aortic c o ar ctation . 2025. arXiv: 2502.12097 [math.NA] . [45] Seth uraman Sank aran et al. “ Uncertaint y quantification in coronary blo o d flow sim ulations: Impact of geometry , b oundary conditions and blo od viscosity”. In: Journal of Biome chanics 49.12 (Aug. 2016), pp. 2540–2547. doi : 10.1016/j.jbiomech.2016.01.002 . [46] Riccardo T enderini et al. “ Deformable registration and generative mo delling of aortic anatomies by auto-deco ders and neural ODEs”. In: npj Biolo gic al Physics and Me chanics 2.1 (2025), p. 26. 30 REFERENCES [47] Ben te Thamsen et al. “ Syn thetic database of aortic morphometry and hemo dynamics: ov ercom- ing medical imaging data a v ailability”. In: IEEE T r ansactions on Me dic al Imaging 40.5 (2021), pp. 1438–1449. [48] Ben te Thamsen et al. “ Unsup ervised learning and statistical shap e mo deling of the morphometry and hemo dynamics of coarctation of the aorta”. In: International Confer enc e on Me dic al Image Computing and Computer-Assiste d Intervention . Springer. 2020, pp. 776–785. [49] Alain T rouvé. “ Diffeomorphisms Groups and Pattern Matching in Image Analysis”. In: International Journal of Computer Vision 28.3 (July 1998), pp. 213–221. doi : 10.1023/a:1008001603737 . [50] P auli Virtanen et al. “ SciPy 1.0: F undamen tal Algorithms for Scientific Computing in Python”. In: Natur e Metho ds 17 (2020), pp. 261–272. doi : 10.1038/s41592- 019- 0686- 2 . [51] Nico W esterhof, Jan-Willem Lankhaar, and Berend E W esterhof. “ The arterial windk essel”. In: Me dic al & biolo gic al engine ering & c omputing 47.2 (2009), pp. 131–141. REFERENCES 31 Appendix A. Training of the conditional LDDMM stochastic interpolant T o approximate the conditioned drift b θ t [ c ] we employ a graph neural net work (GNN), whose architec- ture is a MeshGraphNet implemented in the op en-source framew ork physicsnemo [ 41 ]. The pro cessor size is 10 , and the hidden dimension of the no de and edge enco der and deco der is 128 . The input lay er is defined by the template hexahedral mesh, whose num b er of vertices will b e denoted by N 0 . No de features are a concatenation of a feature em b edding of I t ∈ R N 0 × 3 , considering, for n f = 64 F ourier features, • a feature em b edding φ I t : R N 0 × 3 → R N 0 × 3 · n f · 2 , defined as φ I t ( I t ) = [cos( v ⊗ I t ) , sin( v ⊗ I t )] , (27) where ⊗ is the Kroneck er pro duct of I t ∈ R N 0 × 3 with v = { 1 , 2 , . . . , 2 n f − 1 } ∈ R n f ; • a time sin usoidal enco ding φ t : R → R 3 · n f · 2 defined as φ t = φ 2 t ◦ φ 1 t , where φ 1 t : R → R d , φ 1 t ( t ) =  sin( ω 0 t ) , cos( ω 0 t ) , . . . , sin( ω n f / 2 − 1 t ) , cos( ω n f / 2 − 1 t )  , (28) where the frequencies are given by ω i = exp  − ln(10000) · 2 i d  , for i = 0 , . . . , d/ 2 − 1 , with d = 94 · 4 b eing the em b edding dimension, and φ 2 t : R d → R 3 · n f · 2 is a FNN architecture with 1 hidden lay er of dimension 128 and R eLU activ ation, while the final activ ation is a tanh activ ation; • an additional feed-forw ard neural netw ork φ c : R d → R 3 · n f · 2 with the same structure as φ 2 t . F or defining the underlying graph structure, instead of the mesh edges, we consider each node to b e connected to its k = 12 nearest neigh b ors. This choice improv es the information propagation in the GNN. The training is p erformed for 2750 epo c hs on 1209 training data with a batch size n b = 4 , using the A dam optimizer with learning rate 10 − 4 . A scheduler that reduces the learning rate by a factor of 0 . 5 ev ery time there is a plateau of 100 ep ochs in the training error is employ ed. A sigma noise schedule [ 27 ] is emplo yed for the diffusion co efficien t σ t : ρ = 1 and the initial and final v alues are set to σ max = 0 . 002 and σ min = 0 . 001 , resp ectively σ t =  σ 1 /ρ max + t ( σ 1 /ρ min − σ 1 /ρ max )  ρ . After the training a fine-tuning phase is performed on the 30 test shapes for 1750 epo c hs with the same training parameters. The L 2 -relativ e errors on the drift ev aluated on the test set of 30 shap es during the fine-tuning phases are sho wn in Figure 28 . Figure 29 shows the L 2 -relativ e error during training and fine-tuning phases for different sampled time instances, as well as the Chamfer distance b et ween the target shap es and the generated shap es with the learned drift. 0 5 0 0 1 0 0 0 1 5 0 0 2 0 0 0 2 5 0 0 Epochs 1 0 1 D r i f t m e a n L 2 - r e l a t i v e e r r o r Er r or on test dataset Initial training on 1209 shapes F ine-tuning on 52 test shapes Figure 28. L 2 -relativ e errors on the drift during training (blue curv e) and fine-tuning (red curve) phases, ev aluated only on the test set of 30 shap es. Appendix B. Smoothing algorithms f or mesh aspect ra tio optimiza tion Let Ψ denote a v olume displacement field obtained b y solving the iterative linear elastic problem ( 21 ) (section 4.4 ). F or each (hexahedral) mesh cell K after mesh deformation, w e consider the trilinear reference map F K ;Ψ ( x ) : [0 , 1] 3 → K. 32 REFERENCES Drif t L2 r elative er r or T raining over 1209 training shapes Chamfer distance Drif t L2 r elative er r or Sampled time interval T est inde x F ine tuning over 52 test shapes Chamfer distance Sampled time interval T est inde x Figure 29. Left: Drift L 2 -relativ e error on the test set during training and fine-tuning for differen t sampled time instances. Right: Chamfer distance b et w een the target shapes and the generated shap es with the learned drift on the test set during training and fine- tuning. defined by mapping the vertices of a reference cub e to the vertices of K . Under the assumption that problem ( 21 ) was solv ed successfully , the Jacobian J = ∇ F K ;Ψ ( x ) ∈ R 3 × 3 is well-defined and non-singular for each p oin t x ∈ [0 , 1] 3 . In particular, since the deformation is assumed to preserve cell orientation, the Jacobian has three p ositive singular v alues 0 < σ 1 ( x ) ≤ σ 2 ( x ) ≤ σ 3 ( x ) at each p oint x ∈ [0 , 1] 3 . W e define the asp ect ratio of the cell K under the deformation Ψ as the maximum of the p oin twise ratio b etw een the largest and the smallest singular v alue of J : Asp( K ; Ψ) = max x ∈ [0 , 1] 3 σ 3 ( x ) σ 1 ( x ) . (29) R emark B.1 . The asp ect ratio ( 29 ) can b e seen as the deviation of the shap e of K from a cub e. Indeed, Asp( K ; Ψ) ≥ 1 , and the minimum v alue of 1 corresp onds to the case of all equal singular v alues, at each p oin t x ∈ [0 , 1] 3 . In practice, we consider an approximated asp ect ratio g Asp( K ; Ψ) = max x ∈ S asp σ 3 ( x ) σ 1 ( x ) , (30) using uniformly spaced p oints S asp =  0 , 1 4 , 1 2 , 3 4 , 1  3 . Let us now define, for a given displacement field, the set I (Ψ) =  K : g Asp( K ; Ψ) > 3 4 max K ′ g Asp( K ′ ; Ψ)  of b ad c el ls , i.e., of mesh cells whose approximated asp ect ratio is larger than 75% of the maximum v alue o ver the whole domain. W e consider the optimization problem Ψ ∗ = argmin Ψ | ∂ Ω =Ψ s L (Ψ) , L (Ψ) := 1 # I (Ψ) X K ∈I (Ψ) g Asp( K ; Ψ) 2 . (31) where Ψ s is the prescribed deformation on the surface mesh (see equation ( 21 )). The minimization problem ( 31 ) has b een solv ed with PyTorch using a mo dified adaptive gradient descent metho d with an initial v alue of Ψ N steps , stopping when a threshold of L < 100 w as reached, which in each case corresponded to a maximum approximate asp ect ratio b etw een 10 and 20 . REFERENCES 33 T able 2. Mesh v ariation statistics. Here d c is the mean pairwise chamfer distance ov er eac h batch, whereas σ is the ro ot mean square sample standard deviation of the vertex co ordinates. Both are in units of millimeters. α r = 0 . 7 d c σ α r = 1 d c σ α r = 1 . 3 d c σ 0 1.419 1.79 0 1.923 2.902 0 2.043 2.875 1 1.561 2.389 1 1.871 2.901 1 1.712 2.207 2 1.223 1.694 2 1.359 1.938 2 1.735 2.791 3 1.302 1.848 3 1.483 1.949 3 1.602 1.966 4 1.42 2.031 4 1.432 1.943 4 1.631 2.674 5 1.444 2.019 5 1.904 3.128 5 1.571 2.12 6 1.497 2.058 6 1.861 2.602 6 2.035 2.974 7 — — 7 1.54 2.494 7 1.56 2.131 8 2.041 3.192 8 1.634 2.208 8 1.672 2.311 9 1.562 2.38 9 1.89 3.07 9 1.457 1.815 10 1.773 2.14 10 1.894 2.367 10 1.626 1.682 11 1.764 2.615 11 1.695 2.15 11 1.989 2.51 12 1.592 2.933 12 1.568 2.468 12 1.8 2.893 13 2.111 3.639 13 1.675 2.583 13 1.351 1.641 14 1.547 2.465 14 1.597 2.431 14 1.601 2.255 15 1.899 3.077 15 1.931 2.785 15 1.803 2.299 16 — — 16 1.586 2.767 16 1.795 2.854 17 1.611 2.341 17 1.576 2.101 17 1.93 2.789 18 2.099 3.273 18 1.785 2.889 18 1.253 1.333 19 2.122 3.53 19 1.665 2.288 19 1.539 2.013 20 1.617 2.494 20 1.622 2.382 20 1.46 1.985 21 1.453 2.074 21 1.418 1.677 21 2.059 3.355 22 1.172 1.518 22 1.674 2.566 22 1.735 2.572 23 1.412 1.898 23 1.968 2.825 23 1.956 3.214 24 1.471 2.167 24 1.797 2.832 24 1.802 2.507 25 1.591 2.413 25 1.543 2.103 25 1.572 2.445 26 1.523 2.16 26 1.799 2.814 26 1.47 2.105 27 1.595 2.069 27 2.095 3.103 27 1.984 2.948 28 1.532 2.296 28 1.467 1.796 28 1.531 2.098 29 1.586 2.534 29 1.427 1.893 29 1.467 2.454 A dditional changes were made to prev ent the reemergence of inv erted cells in the course of this pro ce- dure or conv ergence to po or local minima. Note also that as we are using a multigrid n umerical solv er, this procedure was actually p erformed in sequence from the coarsest to the finest lev el, fixing inherited co ordinates at each step. Appendix C. Additional det ails on biomarker st a tistics This supplementary section details mesh v ariation statistics (table 2 ), as well as sample means and standard deviations for v elo city (table 3 ), pressure (table 4 ), wall shear stress (table 5 ), and OSI (table 6 ) for all considered batches. 34 REFERENCES T able 3. Maximum v elo cit y magnitude (in m/s) per batch at t = 2 . 16 (peak systole), sample means and standard deviations. α r = 0 . 7 µ σ α r = 1 . 0 µ σ α r = 1 . 3 µ σ 0 4 . 108 0 . 696 0 1 . 750 0 . 177 0 1 . 368 0 . 109 1 7 . 032 0 . 828 1 4 . 886 0 . 800 1 3 . 059 0 . 250 2 5 . 395 0 . 418 2 3 . 850 0 . 393 2 2 . 815 0 . 313 3 4 . 098 0 . 142 3 2 . 398 0 . 076 3 1 . 778 0 . 072 4 5 . 198 0 . 490 4 3 . 903 0 . 568 4 2 . 590 0 . 423 5 3 . 774 0 . 260 5 2 . 402 0 . 465 5 1 . 992 0 . 346 6 5 . 424 0 . 925 6 3 . 801 0 . 579 6 2 . 712 0 . 315 7 — — 7 2 . 694 0 . 228 7 1 . 906 0 . 151 8 3 . 031 0 . 277 8 2 . 412 0 . 316 8 2 . 147 0 . 255 9 4 . 464 0 . 602 9 2 . 838 0 . 209 9 1 . 975 0 . 278 10 1 . 811 0 . 166 10 1 . 081 0 . 097 10 0 . 857 0 . 048 11 3 . 298 0 . 419 11 1 . 806 0 . 189 11 1 . 416 0 . 130 12 8 . 675 1 . 354 12 4 . 576 0 . 407 12 3 . 708 0 . 273 13 4 . 019 0 . 562 13 1 . 986 0 . 148 13 1 . 510 0 . 066 14 4 . 504 0 . 317 14 2 . 483 0 . 215 14 1 . 725 0 . 170 15 3 . 545 0 . 358 15 1 . 951 0 . 221 15 1 . 324 0 . 069 16 5 . 503 — 16 3 . 097 0 . 243 16 2 . 078 0 . 183 17 4 . 538 0 . 548 17 2 . 152 0 . 205 17 1 . 387 0 . 202 18 4 . 159 0 . 446 18 2 . 663 0 . 381 18 1 . 920 0 . 170 19 3 . 618 0 . 675 19 2 . 047 0 . 122 19 1 . 919 0 . 065 20 5 . 007 0 . 470 20 3 . 814 0 . 207 20 3 . 074 0 . 299 21 6 . 123 0 . 625 21 3 . 998 0 . 257 21 2 . 831 0 . 430 22 4 . 122 0 . 268 22 4 . 312 0 . 742 22 2 . 568 0 . 348 23 4 . 306 0 . 374 23 3 . 040 0 . 562 23 2 . 301 0 . 220 24 5 . 132 0 . 423 24 2 . 810 0 . 223 24 1 . 669 0 . 154 25 5 . 801 0 . 695 25 2 . 370 0 . 260 25 1 . 819 0 . 150 26 5 . 828 0 . 626 26 4 . 139 0 . 511 26 3 . 592 0 . 295 27 2 . 835 0 . 236 27 1 . 646 0 . 129 27 0 . 978 0 . 049 28 3 . 730 0 . 518 28 2 . 405 0 . 265 28 1 . 700 0 . 141 29 6 . 422 0 . 965 29 3 . 744 0 . 326 29 2 . 880 0 . 159 REFERENCES 35 T able 4. Largest pressure difference b et ween inlet and descending aorta ov er one car- diac cycle (in mmHg), sample means and standard deviations. α r = 0 . 7 µ σ α r = 1 . 0 µ σ α r = 1 . 3 µ σ 0 28 . 644 1 . 120 0 11 . 824 0 . 427 0 8 . 618 0 . 377 1 88 . 027 12 . 313 1 43 . 453 6 . 466 1 21 . 899 0 . 787 2 63 . 888 5 . 698 2 36 . 947 4 . 547 2 20 . 059 1 . 510 3 42 . 669 1 . 190 3 20 . 990 0 . 762 3 13 . 680 0 . 311 4 40 . 048 6 . 825 4 21 . 498 2 . 175 4 14 . 061 0 . 969 5 29 . 633 1 . 412 5 16 . 304 0 . 374 5 12 . 987 0 . 387 6 47 . 970 5 . 846 6 29 . 287 3 . 956 6 19 . 765 1 . 208 7 — — 7 12 . 421 0 . 615 7 8 . 622 0 . 392 8 31 . 649 4 . 326 8 19 . 859 1 . 453 8 13 . 917 0 . 595 9 41 . 759 6 . 948 9 23 . 142 1 . 697 9 14 . 836 0 . 796 10 20 . 012 0 . 942 10 11 . 997 0 . 444 10 9 . 477 0 . 361 11 21 . 471 1 . 368 11 11 . 015 0 . 429 11 9 . 040 0 . 379 12 111 . 399 22 . 018 12 31 . 068 4 . 557 12 22 . 435 1 . 983 13 42 . 676 6 . 047 13 20 . 369 1 . 047 13 14 . 811 0 . 413 14 45 . 675 4 . 582 14 22 . 484 1 . 245 14 14 . 173 0 . 314 15 33 . 101 2 . 233 15 17 . 511 0 . 476 15 11 . 313 0 . 386 16 57 . 824 — 16 23 . 487 1 . 808 16 14 . 042 0 . 675 17 31 . 135 2 . 935 17 18 . 584 0 . 527 17 11 . 783 0 . 375 18 37 . 853 4 . 359 18 19 . 838 0 . 644 18 14 . 422 0 . 452 19 30 . 507 4 . 736 19 15 . 413 0 . 421 19 12 . 476 0 . 321 20 59 . 358 10 . 298 20 31 . 828 2 . 152 20 20 . 841 0 . 719 21 53 . 520 5 . 594 21 30 . 180 1 . 310 21 17 . 860 1 . 253 22 46 . 991 3 . 758 22 33 . 863 3 . 379 22 18 . 590 0 . 829 23 41 . 321 4 . 160 23 24 . 068 2 . 784 23 14 . 496 0 . 500 24 41 . 528 2 . 528 24 22 . 245 0 . 666 24 14 . 692 0 . 287 25 49 . 545 7 . 475 25 23 . 881 0 . 559 25 16 . 854 0 . 383 26 70 . 156 5 . 729 26 35 . 460 3 . 147 26 25 . 416 1 . 248 27 28 . 887 2 . 554 27 18 . 003 0 . 955 27 12 . 013 0 . 337 28 34 . 776 2 . 610 28 19 . 231 1 . 041 28 14 . 515 0 . 503 29 86 . 733 19 . 557 29 39 . 778 3 . 418 29 22 . 793 0 . 496 36 REFERENCES T able 5. Mean w all shear stress magnitude ov er the whole blo od v essel wall and one cardiac cycle (in Pa), sample means and standard deviations. α r = 0 . 7 µ σ α r = 1 . 0 µ σ α r = 1 . 3 µ σ 0 3 . 671 0 . 273 0 1 . 532 0 . 133 0 1 . 137 0 . 118 1 7 . 392 0 . 968 1 4 . 333 0 . 682 1 2 . 409 0 . 248 2 7 . 223 0 . 659 2 4 . 291 0 . 337 2 2 . 669 0 . 342 3 4 . 276 0 . 172 3 2 . 173 0 . 132 3 1 . 327 0 . 086 4 4 . 931 0 . 604 4 2 . 533 0 . 192 4 1 . 699 0 . 190 5 3 . 113 0 . 155 5 1 . 478 0 . 123 5 1 . 097 0 . 081 6 5 . 560 0 . 719 6 3 . 488 0 . 590 6 2 . 341 0 . 321 7 — — 7 2 . 849 0 . 360 7 1 . 719 0 . 241 8 2 . 960 0 . 513 8 1 . 800 0 . 202 8 1 . 309 0 . 092 9 5 . 386 0 . 791 9 2 . 949 0 . 312 9 1 . 663 0 . 145 10 1 . 678 0 . 137 10 1 . 018 0 . 093 10 0 . 766 0 . 043 11 2 . 766 0 . 374 11 1 . 272 0 . 107 11 1 . 003 0 . 115 12 10 . 617 1 . 704 12 4 . 292 0 . 516 12 2 . 983 0 . 365 13 4 . 313 0 . 751 13 2 . 020 0 . 208 13 1 . 331 0 . 081 14 3 . 927 0 . 410 14 2 . 129 0 . 183 14 1 . 355 0 . 114 15 3 . 781 0 . 452 15 1 . 933 0 . 197 15 1 . 064 0 . 082 16 6 . 217 — 16 3 . 235 0 . 359 16 1 . 821 0 . 255 17 4 . 229 0 . 523 17 2 . 135 0 . 152 17 1 . 369 0 . 117 18 4 . 053 0 . 539 18 2 . 027 0 . 161 18 1 . 403 0 . 095 19 3 . 061 0 . 514 19 1 . 635 0 . 126 19 1 . 258 0 . 047 20 6 . 283 1 . 070 20 3 . 493 0 . 415 20 2 . 208 0 . 111 21 4 . 877 0 . 388 21 3 . 137 0 . 141 21 1 . 820 0 . 256 22 4 . 686 0 . 303 22 3 . 414 0 . 352 22 2 . 162 0 . 214 23 4 . 784 0 . 438 23 2 . 653 0 . 471 23 1 . 795 0 . 200 24 5 . 091 0 . 576 24 2 . 714 0 . 183 24 1 . 707 0 . 111 25 4 . 832 0 . 728 25 2 . 225 0 . 122 25 1 . 579 0 . 081 26 6 . 546 0 . 732 26 3 . 611 0 . 409 26 2 . 408 0 . 146 27 2 . 839 0 . 340 27 1 . 616 0 . 167 27 0 . 931 0 . 064 28 4 . 247 0 . 374 28 2 . 242 0 . 213 28 1 . 593 0 . 130 29 6 . 724 1 . 170 29 3 . 364 0 . 328 29 2 . 105 0 . 115 REFERENCES 37 T able 6. Mean oscillatory shear index ov er a reference patc h on the underside of the descending aorta, sample means and standard deviations. α r = 0 . 7 µ σ α r = 1 . 0 µ σ α r = 1 . 3 µ σ 0 0 . 260 0 . 019 0 0 . 209 0 . 022 0 0 . 231 0 . 014 1 0 . 119 0 . 060 1 0 . 178 0 . 057 1 0 . 250 0 . 040 2 0 . 142 0 . 020 2 0 . 081 0 . 024 2 0 . 155 0 . 016 3 0 . 307 0 . 027 3 0 . 226 0 . 039 3 0 . 197 0 . 019 4 0 . 077 0 . 045 4 0 . 044 0 . 014 4 0 . 088 0 . 028 5 0 . 114 0 . 048 5 0 . 125 0 . 017 5 0 . 239 0 . 023 6 0 . 222 0 . 030 6 0 . 243 0 . 025 6 0 . 203 0 . 040 7 — — 7 0 . 131 0 . 017 7 0 . 181 0 . 022 8 0 . 149 0 . 027 8 0 . 144 0 . 024 8 0 . 163 0 . 026 9 0 . 056 0 . 029 9 0 . 050 0 . 009 9 0 . 082 0 . 013 10 0 . 246 0 . 033 10 0 . 243 0 . 024 10 0 . 191 0 . 023 11 0 . 213 0 . 035 11 0 . 201 0 . 016 11 0 . 225 0 . 013 12 0 . 225 0 . 028 12 0 . 156 0 . 024 12 0 . 152 0 . 023 13 0 . 137 0 . 052 13 0 . 117 0 . 033 13 0 . 149 0 . 018 14 0 . 220 0 . 057 14 0 . 197 0 . 063 14 0 . 159 0 . 025 15 0 . 219 0 . 049 15 0 . 227 0 . 011 15 0 . 211 0 . 026 16 0 . 045 — 16 0 . 102 0 . 038 16 0 . 155 0 . 028 17 0 . 273 0 . 026 17 0 . 248 0 . 012 17 0 . 236 0 . 015 18 0 . 120 0 . 028 18 0 . 116 0 . 040 18 0 . 115 0 . 019 19 0 . 220 0 . 043 19 0 . 254 0 . 022 19 0 . 216 0 . 018 20 0 . 149 0 . 039 20 0 . 144 0 . 062 20 0 . 170 0 . 038 21 0 . 244 0 . 020 21 0 . 230 0 . 027 21 0 . 203 0 . 042 22 0 . 185 0 . 031 22 0 . 218 0 . 027 22 0 . 204 0 . 013 23 0 . 171 0 . 037 23 0 . 226 0 . 018 23 0 . 278 0 . 027 24 0 . 237 0 . 041 24 0 . 245 0 . 019 24 0 . 270 0 . 012 25 0 . 174 0 . 042 25 0 . 118 0 . 020 25 0 . 164 0 . 019 26 0 . 119 0 . 029 26 0 . 148 0 . 021 26 0 . 237 0 . 041 27 0 . 138 0 . 012 27 0 . 184 0 . 026 27 0 . 177 0 . 031 28 0 . 097 0 . 016 28 0 . 192 0 . 035 28 0 . 208 0 . 012 29 0 . 132 0 . 024 29 0 . 147 0 . 030 29 0 . 208 0 . 024

Original Paper

Loading high-quality paper...

Comments & Academic Discussion

Loading comments...

Leave a Comment