Design and execution of a verification, validation, and uncertainty quantification plan for a numerical model of left ventricular flow after LVAD implantation

Background Left ventricular assist devices (LVADs) are implantable pumps that act as a life support therapy for patients with severe heart failure. Despite improving the survival rate, LVAD therapy can carry major complications. Particularly, the flow distortion introduced by the LVAD in the left ventricle (LV) may induce thrombus formation. While previous works have used numerical models to study the impact of multiple variables in the intra-LV stagnation regions, a comprehensive validation analysis has never been executed. The main goal of this work is to present a model of the LV-LVAD system and to design and follow a verification, validation and uncertainty quantification (VVUQ) plan based on the ASME V&V40 and V&V20 standards to ensure credible predictions. Methods The experiment used to validate the simulation is the SDSU cardiac simulator, a bench mock-up of the cardiovascular system that allows mimicking multiple operation conditions for the heart-LVAD system. The numerical model is based on Alya, the BSC’s in-house platform for numerical modelling. Alya solves the Navier-Stokes equation with an Arbitrary Lagrangian-Eulerian (ALE) formulation in a deformable ventricle and includes pressure-driven valves, a 0D Windkessel model for the arterial output and a LVAD boundary condition modeled through a dynamic pressure-flow performance curve. The designed VVUQ plan involves: (a) a risk analysis and the associated credibility goals; (b) a verification stage to ensure correctness in the numerical solution procedure; (c) a sensitivity analysis to quantify the impact of the inputs on the four quantities of interest (QoIs) (average aortic root flow QAoavg, maximum aortic root flow QAomax, average LVAD flow QVADavg, and maximum LVAD flow QVADmax); (d) an uncertainty quantification using six validation experiments that include extreme operating conditions. Results Numerical code verification tests ensured correctness of the solution procedure and numerical calculation verification showed a grid convergence index (GCI)95% <3.3%. The total Sobol indices obtained during the sensitivity analysis demonstrated that the ejection fraction, the heart rate, and the pump performance curve coefficients are the most impactful inputs for the analysed QoIs. The Minkowski norm is used as validation metric for the uncertainty quantification. It shows that the midpoint cases have more accurate results when compared to the extreme cases. The total computational cost of the simulations was above 100 [core-years] executed in around three weeks time span in Marenostrum IV supercomputer. Conclusions This work details a novel numerical model for the LV-LVAD system, that is supported by the design and execution of a VVUQ plan created following recognised international standards. We present a methodology demonstrating that stringent VVUQ according to ASME standards is feasible but computationally expensive.


S1 Verification S1.1 Mesh convergence metrics
The root mean square error (RMSE) between the solution i and j is defined thought the L2 norm ||.|| 2 as: Once all the cases are computed and the errors ϵ 1,2 and ϵ 2,3 computed the observed order of convergence is calculated as [1]: p i,k = [1/ ln(r j,k )] [ln |ϵ i,j /ϵ j,k | + q i,k (p i,k )] Note that in our case r i,j = r j,k = r = 2.0, so q i,k (p i,k ) = 0.0 and therefore the previous system of equations is reduced to: For three velocity fields computed u 1 , u 2 , u 3 , u 1 being the coarsest, for three subdivision levels, the order of the convergence of the numerical scheme is given by: With the observed p value, the grid convergence index (GCI) can be computed as [1]: This uncertainty estimate provides an interval f ± U 95% within the true mathematical value f T falls with a probability of 95%.

S1.2 Numerical code verification
Numerical code verification is executed as per Section 2 of [2], for a 2D Poiseuille and a 3D Womersley flow problem in a cylindrical tube. These problems have non-trivial analytical solutions that are used as true value. For both cases the discretisation error is monitored as the grid is systematically refined by halving as in [3]. If the ratio between mesh subdivisions is defined as r i,j = r i /r j then, for this case r 1,2 = r 2,3 = r = 2.0, a figure considerably larger than 1.3, the minimum value recommended [2]. The velocity field is the quantity of interest (QoI) to be verified as it is also the raw variable used calculated in the numerical model. The mesh convergence metrics are described in Section S1.1.

S1.2.1 Software Quality Assurance
The finite elements method (FEM) software is developed with a continuous integration and continuous deployment (CD/CI) strategy based on Git 1 , which combines feature-driven development and feature branches with issue tracking. Git pipelines ensure continuous integration, running a series of software checks, builds, and regression tests when the developers modify the source code. The build evaluation includes 27 combinations of architectures (Intel, IBM), compilers (gnu, intel, pg, xl), and optimization options, running more than 200 regression tests executed with various MPI and OpenMP configurations, more than 4000 different executions. This method helps to detect bugs early in the development cycle and guarantee the correctness of the simulation results and the software's stability. On this manuscript Alya version 2.1 was used, available in the Unified European Applications Benchmark Suite (UEABS) [4]. Further description of the code and training material can be found in the following links 2 3 . Following the "Credible practice of modeling and simulation in healthcare" [5], a data availability statement is included. The referenced section details how to gather the data to reproduce this work.  The outlet is free with weakly imposed pressure equal to zero. The other two walls Γ w have no slip condition v i | Γw = 0.
Results: The calculation of the convergence was calculated at the nodes of the coarsest mesh, to avoid any interpolation of the solution. Results for each case are sumarised in Table 1. The observed order of convergence p was equal to 1.907, compatible with the theoretical order of convergence of 2 of the 2nd order backward differentiation formula (BDF) time scheme used.

S1.2.3 3D Womersley flow
where: For a representation of the analytical solution refer to the blue line in Fig. 1.
Results : Figure 1 shows the analytical (Eq. (9)) and the simulated solution. The maximum velocity reaches expected periodic behaviour after 3[s]. Table 2

S1.3 Numerical calculation verification
The original model was subdivided 2 times, splitting each element into 8, to obtain meshes with 6.6M , 53.1M , and 425M elements respectively. The time step used is ∆t = 0.00428 [s]. If the critical time step is ∆t c then, ∆t/∆t c = 10 −3 [−].

S1.3.1 Ventricular geometry with stationary boundary conditions
The goal is to estimate the stationary error and observed convergence order for the computational fluid dynamics (CFD) solver on the mesh created for the ventricular geometry. Geometry deformation, valve model, and the pump boundary condition are turned off as they depend on measurements obtained from the CFD solver. As such geometry does not have an analytical solution, the RMSE are calculated against the finest mesh computed. Results are summarised in Table 3. The observed order of convergence is p obs = 0.971 compatible with the theoretical order of convergence of 1 provided by the first order trapezoidal time integration.

S1.3.2 Ventricular geometry with transient boundary conditions
Here we show the RMSE for the complete model described in the main document.. As the complete model have transient boundary conditions, the RMSE is provided in time-averaged quantities. RMSEs are calculated against the finest mesh Results are summarised on Table 4. The time-averaged observed order of convergence is p obs = 0.85 compatible with the theoretical order of convergence of 1 provided by the first order trapezoidal time integration.