Global Sensitivity Analysis of Four Chamber Heart Hemodynamics Using Surrogate Models

Computational Fluid Dynamics (CFD) is used to assist in designing artificial valves and planning procedures, focusing on local flow features. However, assessing the impact on overall cardiovascular function or predicting longer-term outcomes may requires more comprehensive whole heart CFD models. Fitting such models to patient data requires numerous computationally expensive simulations, and depends on specific clinical measurements to constrain model parameters, hampering clinical adoption. Surrogate models can help to accelerate the fitting process while accounting for the added uncertainty. We create a validated patient-specific four-chamber heart CFD model based on the Navier-Stokes-Brinkman (NSB) equations and test Gaussian Process Emulators (GPEs) as a surrogate model for performing a variance-based global sensitivity analysis (GSA). GSA identified preload as the dominant driver of flow in both the right and left side of the heart, respectively. Left-right differences were seen in terms of vascular outflow resistances, with pulmonary artery resistance having a much larger impact on flow than aortic resistance. Our results suggest that GPEs can be used to identify parameters in personalized whole heart CFD models, and highlight the importance of accurate preload measurements.

u + (u − w) · ∇u − ∇ · σ(u, ) + (u − u ) = 0 in R + × Ω( ), σn − ((u − w) · n) − = h on Γ out ow ( ), with the time dependent uid domain Ω( ) de ned as using the ALE mapping d transforming an arbitrary reference con guration Ω 0 into the current uid domain Ω( ). Here , u, and w := d represent the uid pressure, the ow velocity, and the ALE velocity respectively, is the dynamic viscosity and the density. The volume penalization term ( ,x) (u( , x) − u ( , x)) is commonly known as Darcy drag which is characterized by the permeability ( , x). Another interpretation of the Darcy drag can be given as resistance term, as for example aerodynamic drag is also interpreted as a resistance. In (1) the Darcy drag is modi ed to enforce correct no-slip conditions for obstacles moving with the obstacle velocity u (x, ). The uid stress tensor σ(u, ) and strain rate tensor (u, ) are de ned as follows: σ(u, ) = − I + 2 (u, ), (u, ) = 1 2 ∇u + (∇u) .
The ALE domain Ω( ) is split up into three time dependent sub-domains by means of the permeability ( , x), namely the uid sub-domain Ω ( ), the porous sub-domain Ω ( ) and the solid sub-domain Ω ( ).
In Ω ( ) the classical ALE-Navier-Stokes equations are recovered, while in Ω the full ALE-NSB equations describe uid owing trough a moving porous medium, u and are understood in an averaged sense in this context. In Ω ( ) the velocity u is approaching u and thus asymptotically satisfying the no-slip condition on the Ω ( )/Ω ( ) interface. Note that even in the case where → 0 + the penalization term has a well de ned limit, see [S1].

S.I.A. Hemodynamic A erload Models
Modeling of afterload for hemodynamics is modeled by using a 0D Windkessel model. This means we de ne h in (4) as h := − WK ( )n with the Windkessel pressure WK is governed by the di erential algebraic system [S16] In the case of multiple Windkessel outlets we will use the same notation for variables with an added subscript indicating multiple outlets. Tools for personalization of the individual Windkessel parameters can be found in [S33].

S.I.B. Variational Formulation and Numerical Stabilization
Following [S7, S8] the discrete variational formulation of (1) including the boundary conditions (4), (5) and (3) can be stated in the following abstract form: with the bilinear form of the NSB equations the bilinear form RBVMS , which will be explained later in Equation (21), and the right hand side contribution We use standard notation to describe the nite element function space [S 1 ℎ,g,Γ (T )] 3 as a conformal vector-valued trial space of piece-wise linear, globally continuous basis functions v ℎ over a tessellation T of Ω into nite elements constrained by v ℎ = g on Dirichlet boundaries, more precisely where for this work Γ := Γ noslip ∪ Γ in ow . For further details we refer to [S10, S40]. As previously described in [S26] we utilize the residual based variational multiscale (RBVMS) formulation as proposed in [S7, S8], providing turbulence modeling in addition to numerical stabilization. In the following we give a short summary of the changes necessary to use RBVMS methods for the ALE-NSB equations. Brie y, the RBVMS formulation is based on a decomposition of the solution and weighting function spaces into coarse and ne scale subspaces and the corresponding decomposition of the velocity and the pressure and their respective test functions. Henceforth the ne scale quantities and their respective test functions shall be denoted with the superscript . We assume u = u ℎ , quasi-static ne scales ( u = 0), as well as v ℎ = 0, u = 0 on Ω( ) and incompressibility conditions for u ℎ and u . The ne scale pressure and velocity are approximated in an element-wise manner by means of the residuals r and .
The residuals of the NSB equations and the incompressibility constraint are: Taking all assumptions into consideration and employing the scale decomposition followed by partial integration yields the bilinear form of the RBVMS formulation RBVMS (v ℎ , ℎ ; u ℎ , ℎ ), The residuals (19) and (20) are evaluated for every element Ω ∈ T . Following [S37] the stabilization parameters SUPS and LSIC are de ned as: Here G is the three dimensional element metric tensor de ned per nite element as with J being the Jacobian of the transformation of the reference element to the physical nite element ∈ T , Δ denotes time step size and is a positive constant, taken as 30, derived from an element-wise inverse estimate. For further details see [S7, S8].

S.I.C. Numerical Solution Strategy
Spatio-temporal discretization of all PDEs and the solution of the arising systems of equations relied upon the Cardiac Arrhythmia Research Package (CARPentry), see [S43]. For temporal discretization of the ALE-NSB equations we used the generalized-method, see [S23] with a spectral radius ∞ = 0.2. For updating the Windkessel pressures WK we discretized (11) with an implicit Euler method. For ease of coupling with our CFD solver the calculation of ( ) in (13) is lagged by one Newton iteration. After discretization in space as described in Section S.I.B and temporal discretization using the generalizedintegrator we obtain a nonlinear algebraic system to solve for advancing time from timestep to +1 . A quasi inexact Newton-Raphson method is used to solve this system with linearization approach similar to [S8] adapted to the NSB equations. At each iteration a block system of the form is solved with K ℎ , B ℎ , C ℎ , and D ℎ denoting the Jacobian matrices, Δu, Δp representing the velocity and pressure updates and R upper , R lower indicating the residual contributions. In this regard we use the exible generalized minimal residual method (fGMRES) and e cient preconditioning based on the PCFIELDSPLIT 1 package from the library PETSc [S4-S6] and the incorporated suite HYPRE BoomerAMG [S19]. By extending our previous work [S3, S26, S27] we implemented the methods in the nite element code Cardiac Arrhythmia Research Package (CARPentry) [S42, S43].

S.II. Obstacle Representation
Here we want to give a brief description of how we represent moving obstacles for usage in the ALE-NSB equations. This task is solved by representing obstacles using triangular surface meshes followed by element-wise calculation of the partial volume covered by the obstacle. In the rst step, all nodes within the obstacle are identi ed using the ray casting algorithm [S18, S35]. Subsequently, all elements are split into three categories and receive a corresponding volume fraction value , describing the partial volume covered by the obstacle: • Elements fully covered by the obstacle lie in Ω , consequently = 1.
• Elements outside the obstacle lie in Ω and obtain = 0.
• Elements that are split by the element surface correspond to elements in Ω , hence = 1 https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCFIELDSPLIT.html where denotes the element volume covered by the obstacle and is the total element volume.
This procedure is carried out for every time step and yields a time-dependent, element-based volume fraction distribution ( , x), that serves as a basis to provide a suitable permeability distribution, see Figure S1. In this work we de ne 1 ( ,x) := ( ,x) with being a xed penalization factor, e.g.
= 10 −8 . All permeability distributions in this work have been generated using the open-source software Meshtool 2 , see [S36] and [S17] for algorithmic details. In the case of obstacles that change from open to closed state over time we use a simple scaling function. For example, assume an obstacle representing a heart valve region will be open at a time instance op and it takes dur V time to switch from open to closed we de ne and modify 1 ( ,x) to ( ,x) .
Ω Ω Ω Figure S1: Schematic representation of the distribution associated to an obstacle,which is represented by the red line, at a xed point time t.

S.III. Computation of Residence Times on Moving Domains
Here we will give a brief outline of the methods and algorithms used to compute residence time distributions. The starting point is the following PDE describing the time evolution of a residence time distribution eld originating from [S14, S31]. Given a moving uid domain Ω( ) ⊂ R 3 [ ] and a region of interest ( ) ⊂ Ω( ) the evolution of the time ( , x) spent in by an arbitrarily small uid particle caught at point x ∈ Ω( ) at time can be described as with the uid velocity u, the ALE mesh velocity w and a small arti cial di usion parameter [m 2 s −1 ] added for numerical stability of the method. For this work we have used = 1 × 10 −12 . The uid velocity u as well as the ALE mesh velocity w are assumed to be given functions, e.g. coming from a pre-computed CFD simulation. In our applications we set Γ ( ) = ∅ and Γ ( ) = Ω( ). Furthermore, the region of interest ( ) is assumed as a time-evolving tag region assigned to a particular anatomic region, e.g.: ventricular blood pools, and left atrial appendage respectively. After discretization we have with { ( )} =0 being the time-dependent test and trials functions in the ALE domain. For regular domain movement it is safe to assume that M ℎ ( ) is invertible for all and we can rewrite (26) as Next, we apply the modi ed Crank-Nicholson scheme in time as proposed in [S15] giving where we used the following shorthand notation Equation (28) is our starting point for applying the FCT scheme similar to [S25]. Following the ideas of FEM-FCT methods we de ne the matrices The construction of L ensures zero row and column sums. Instead of (28) we consider which represents a stable low-order scheme whose solution doesn't possess any over or undershoots but su ers from to smeared layers. To correct this behavior and arti cial ux correction vector f * ( +1 , ) is added to the right hand side of (29). The de nition of f * follows from an ad-hoc ansatz with the uxes de ned as and weights ∈ [0, 1]. The representation for follows from rst subtracting (28) from (29) and applying the properties of the matrices M and D. This formulation represents a nonlinear system. In [S25] a linear variant has been proposed which we adapted to the moving-domain case. For this we use the explicit solutionτ to (29), by means of an explicit Euler scheme approximating the solution τ + 1 2 at time step + Δ 2 , readingτ

Insertingτ into (30) and rearranging terms yields
) . Additionally, as suggested in [S29], we employ prelimiting in the The computation of the weights follows the same procedure as in [S25] using Zalesak's algorithm [S44]. We also refer to [S29, S30] for a more detailed overview of the presented method. Computation of the residence time distribution elds have been included as addon in CARPentry. After computation of the residence time distribution we can calculate the residence time RT spend in ( ) over a time period ( 0 , 1 ) [S31] as (a) ( ) de ned as tag region of left ventricular bloodpool.
(b) ( ) de ned as tag region of left atrial appendage. Figure S2: Time-averaged residence time distributions AVG with ( ) de ned through di erent labels in the computational mesh. Time average taken over the nal two heartbeats with beatlength equal to 0.725 s. Figure S2 and Figure S3 show illustrations of time averaged residence time distributions that were generated for this work as part of the sensitivity analysis.

S.IV. Pope's Criterion of Turbulence Resolution
In LES type formulations the resolved velocity eld is fundamentally linked to the numerical method used, hence there is no notion of convergence to the solution of a partial di erential equation [S38]. This leads to the problem that mesh convergence often cannot be established by the classical methods [S12, S21, S24, S32, S39]. To remedy this problem [S38] proposes the use of a measure of turbulence resolution , see (32), utilizing the fraction of turbulent kinetic energy resolved by the grid in question.
In order to obtain a point-wise measure, rather than the kinetic energy itself, the kinetic energy density (x, ) is considered: The resulting point-wise measure of turbulence resolution reads: Here is the turbulent kinetic energy of the residual motions, hence of the motions not resolved by the grid, and stands for the total kinetic energy. may be written as the sum of the resolved turbulent kinetic energy ℎ and the not resolved turbulent kinetic energy : The resolved turbulent kinetic energy ℎ is calculated from (31) using the uctuating part of the resolved uid velocity u , which is given by: Figure S3: Time-averaged residence time distributions RV AVG with ( ) de ned as right ventricular blood pool in the computational mesh. Time average taken over the nal two heartbeats with beatlength equal to 0.725 s.
where u ℎ is an average with respect to time. When considering a constant in ow u ℎ is given by the standard mean over all time steps ( = 1 . . . , hence u ℎ is not time-dependent) : In the case of a pulsatile behavior however a phase average is considered: where is the number of cycles or beats and is the period or beat length. By the use of (32) a criterion for su cient mesh resolution is given: In [S38] a value of = 0.2 is proposed, which corresponds to requiring a minimum of 80% of the total turbulence energy to be resolved.
We calculated this criterion based on the simulation of four beats for di erent spatial and temporal resolutions with average edge lengths of Δ = 3000, 1500, 750 m and time step sizes of Δ = 1450, 725, 362.5 ms for the mesh of the left heart. Figure S4 shows a comparison of the average Pope criterion, with average taken over all elements of the mesh. From this we concluded that the mesh with average edge length of Δ = 750 m and a time step size Δ = 362.5 ms was su cient to resolve the ow dynamics in the left heart. For reasons of practicality we used the same average edge length and time step size for the right side heart.

S.V. Input Parameter Variance E ect on Output Features
This serves as additional interpretation for the results in the main manuscript. Figure S5 shows the extracted temporal signals of all parameter sets for the pressure di erences Δ MV1,2,3,4 in the LA. While Figure S4: Average Pope criterion for the di erent spatial and temporal discretization levels. The average is taken over all elements in the mesh.
there is no strong in uence on the output in systole, one can see a clear variation in the outputs in diastole. Table S1 and Table S2 summarize the results of the cross-validation. Features with an average 2 -score of < 0.5 were excluded from global sensitivity analysis in this work. It can be observed that the 2 scores are on average lower for left heart features. We attribute that to the increased turbulence in left heart blood ow compared to right heart blood ow.

S.VII. Comparing Influence of Time Step Size on Main Vortex Center in the Le Ventricle
To study the in uence of time step size on the accuracy of vortex formation we ran a simulation of left heart blood ow with successive re ned time step sizes Δ = 362.5, 181.25, 90.625 ms. It is known, see [S28] and references therein, that vorticity is not a suitable criterion for vortex identi cation, hence we used used the scaled -criterion de ned as := S 2 , := 1 2 ∇u − ∇u , S := 1 2 ∇u + ∇u , 2 .
The di erence to the standard -criterion is the normalization w.r.t. strain, making it easier to nd adequate thresholds for visualization. We used the threshold value of 2.5 and calculated the isosurfaces Figure S5: Extracted time signals for pressure di erences in the LA for the second heart beat. Diastole is indicated by the shaded blue area in the plots.   Table S3 shows a comparison of the relative error of the identi ed vortex centers and Figure S6 shows a visual comparison of the resulting isosurfaces. From this results we concluded that our simulations accurately resolved the vortices at the chosen time step size of Δ = 362.5 ms coinciding with the results obtained from using Pope's criterion as described in section S.IV.

S.VIII. Sca erplot of the Locations of Failing Samples
Here we include a projection of the training data set used for calibrating the left heart GPEs overlayed with the samples that led to failed simulations. From the result in Figure S7 we conclude that the failed samples appear to be random samples. Figure S6: Comparison of the major vortex centers extracted from isosurfaces of the scaled -criterion at value 2.5 at peak diastole for di erent time discretizations: a) Δ = 362.5ms (orange) and reference time step Δ = 90.625ms (green), and b) Δ = 181.25ms (orange) and reference time step Δ = 90.625ms (green). Additionally the vortex centers are depicted as colored spheres in orange and green.