Numerical simulations for computational hemodynamics in clinical settings require a combination of many ingredients, mathematical models, solvers and patient-specific data. The sensitivity of the solutions to these factors may be critical, particularly when we have a partial or noisy knowledge of data. Uncertainty quantification is crucial to assess the reliability of the results. We present here an extensive sensitivity analysis in aortic flow simulations, to quantify the dependence of clinically relevant quantities to the patient-specific geometry and the inflow boundary conditions. Geometry and inflow conditions are generally believed to have a major impact on numerical simulations. We resort to a global sensitivity analysis, (i.e., not restricted to a linearization around a working point), based on polynomial chaos expansion (PCE) and the associated Sobol' indices. We regard the geometry and the inflow conditions as the realization of a parametric stochastic process. To construct a physically consistent stochastic process for the geometry, we use a set of longitudinal-in-time images of a patient with an abdominal aortic aneurysm (AAA) to parametrize geometrical variations. Aortic flow is highly disturbed during systole. This leads to high computational costs, even amplified in a sensitivity analysis -when many simulations are needed. To mitigate this, we consider here a large Eddy simulation (LES) model. Our model depends in particular on a user-defined parameter called filter radius. We borrowed the tools of the global sensitivity analysis to assess the sensitivity of the solution to this parameter too. The targeted quantities of interest (QoI) include: the total kinetic energy (TKE), the time-average wall shear stress (TAWSS), and the oscillatory shear index (OSI). The results show that these indexes are mostly sensitive to the geometry. Also, we find that the sensitivity may be different during different instants of the heartbeat and in different regions of the domain of interest. This analysis helps to assess the reliability of in silico tools for clinical applications.
As the computational hemodynamics are progressively incorporated into the clinical practice, a rigorous assessment of the reliability of the numerical predictions is crucial. The numerical solution depends on many factors, ranging from the patient-specific geometry reconstructed from medical images to the measurements used as boundary conditions, the parameters of the mathematical model, and of the associated numerical solver. We do not have a perfect knowledge of these factors and a sensitivity analysis intends to assess the dependence of the computational solution and its associated clinical conclusions on the data as well as the empirically chosen numerical parameters of the solver. Uncertainty quantification is an essential instance of sensitivity analysis, focused on the uncertainty of the physical inputs to a mathematical model, to assess the reliability and the robustness of the numerical solutions.
In this work, we present a sensitivity analysis of the computational hemodynamics in the aorta. The aorta is a peculiar site for the shape, the size, and the regime of the flow. Many pathologies are currently investigated using numerical tools (e.g., Refs. [1–5] to mention a few). An accurate sensitivity analysis in this site is of utmost interest for the reliability of numerical tools.
1.1 Methods of Sensitivity Analysis.
Sensitivity analysis can be carried out in different ways. A possible approach is an educated sampling of different values of the inputs under scrutiny to obtain a statistical assessment of the sensitivity. This approach is quite general; it requires many samples for the significance of the results. Another approach resorts to the numerical computation of the so-called sensitivity equations. This requires the calculation of the differential equations solved by the sensitivities, i.e., the (Gateaux) derivatives of the outputs as a function of the inputs. Numerical procedures can solve these equations. This analysis is “local” as the sensitivities are generally computed by Taylor expansion around a working point, i.e., a nominal value of the inputs. Here, we opted for a “global” sensitivity analysis. Global sensitivity analysis quantifies the importance of varying inputs on the outputs of interest by exploring the entire input parameters' domain . Our approach is based on the polynomial chaos expansion (PCE) and the associated Sobol' indices [7–9]. Recent applications of PCE in vascular problems can be found in Refs. [10–15]. With PCE, the inputs are regarded as stochastic processes represented by a polynomial expansion and parametrized by random input variables. The outputs are, consequently, represented by a corresponding stochastic polynomial expansion. Our analysis computes the coefficients of a truncated PCE; the number of samples of the input space required to accomplish this task is dictated by the quadrature formulas used for computing these coefficients. The need for samples is significantly less than for other strategies.
1.2 Inputs of Interest.
There is a common agreement in computational hemodynamics that the inputs having a major impact on the numerical results with clinical relevance are the geometry and the boundary conditions, more than rheological models or fluid-structure interaction: this is the experience of the authors and experts of the field (e.g., T.J.R. Hughes, personal communication, 2016). Geometry and boundary conditions are, therefore, the focus of the present global sensitivity analysis.
A sensitivity analysis for the geometry can be done in different ways, depending on the specific sources of variations under investigation. For instance, if one is interested in the variations induced by the operator's dependence on a reconstruction procedure, a simple approach consists of sampling the reconstruction with different operators and then comparing the results. The number of samples required to have a statistically significant assessment is pretty large. Here, we pursue a novel approach related to the stochastic parametrization. To obtain physically consistent variations of the geometrical domain parametrized by a random variable, we resort to a longitudinal dataset of a patient, i.e., snapshots of the patient's geometry at different times. The patient at hand is affected by an Abdominal Aortic Aneurysm (AAA), and the evolution of the disease is available over a time range of 4 years. With appropriate image-processing techniques (nonrigid registration), we manage to create a continuous representation of the geometry as a stochastic process. This enables using PCE in a physically realistic way. Also, we argue that, in this way, our analysis spans a broad range of possible variations and uncertainties of the geometry (image noise, operator dependence, reconstruction artifacts, the time evolution of the morphology).
Investigations of the impact of uncertainty of inflow and outflow boundary conditions on the simulation results of blood flows are, e.g., in Refs.  and [13–18]. For the physical nature of the problem, with a strong convective field from the proximal to the distal sites, inflow conditions are arguably more impactful than outflow ones. We postulate that an inlet flowrate is prescribed as a stochastic variable, expanded by PCE. We analyze the impact of the variations of the time-waveform on the numerical results and their hemodynamics significance, with a focus on the different instants of the heartbeat when the sensitivity is significant.
Finally, we extend our sensitivity analysis to numerical parameters required by the blood modeling. As aorta is affected by highly disturbed flows, large Eddy simulation (LES) models have been demonstrated to mitigate the high computational costs significantly [2,19]. However, LES models require numerical parameters whose tuning is largely empirical. In particular, in our model analyzed in Ref.  and tested in Refs.  and , the filter radius is expected to have a significant impact on the solution. A local preliminary analysis based on the sensitivities equations in idealized geometries confirmed this guess . Here, we borrow the same tools used for physical inputs to assess our simulations' sensitivity to the filter radius.
1.3 Outputs of Interest.
We analyze the impact of the (physical and numerical) inputs for both an idealized aortic geometry and a patient-specific setting, on the following outputs: the total kinetic energy (TKE), the time average wall shear stress (TAWSS), and the oscillatory shear index (OSI). These quantities were demonstrated to be associated with the growth of abdominal aortic aneurysms [22–26]. We use the name of quantities of interest (QoI) hereafter to indicate these model outputs generically. In this way, we assess the robustness of our numerical model for aortic hemodynamic applications. The PCE-based Sobol' indices analysis identifies the inputs that affect the most the final results, while the confidence intervals computed for the numerical simulations certify the reliability of the model predictions.
Our results rigorously confirm the quite general experience of the geometry being the most impactful input on the results. Moreover, the analysis allows a fine investigation of different phases of the heartbeats and the region of the aorta more sensitive to the input variations. In particular, we found that some outputs are more sensitive to the inputs in regions of the vessel where other indices are more robust. If confirmed, this result suggests that a combined calculation of several indices may be beneficial in terms of the overall robustness of the computational clinical indications to the input variations.
2 Materials and Methods
2.1 The Numerical Model.
Simulating the aortic blood flows is nontrivial due to the relatively large Reynolds numbers during systole and complex morphology [2,27–29]. The direct numerical simulation of these flows requires a refinement level of the space reticulation to be able to capture the flow field at the smallest significant scale. The associated computational cost might be intimidatingly high, especially for time-sensitive clinical applications. Turbulence modeling has been adopted for simulating aortic flows [28,29]. Recently, an LES deconvolution-based Leray model was shown to capture the properties of a disturbed flow using a relatively coarse discretization finite element mesh, reducing the computational cost . This model can also suppress backflow instability, which is a typical numerical drawback caused by Neumann boundary conditions at outflow branches .
Here, for the sake of the notation, the boundary of Ω is split into two nonoverlapping portions ΓD (Dirichlet conditions) and ΓN (Neumann conditions).
Here, is an indicator function that identifies the regions of the fluid domain requiring stabilization, and δ is a user-defined parameter called filter radius, which quantifies the strength of the stabilization. For (Helmholtz operator), one choice advocated in Refs.  and  is (deconvolution based indicator function). In this study we refer to the so-called Evolve-Filter-Relax implementation of this model, as detailed in Ref. . The radius δ should be selected to find the tradeoff between the stabilization (preferring a large δ) and the accuracy (suggesting a small δ).
In what follows, we refer to the solver previously developed in Ref. , where the equations are discretized with the finite element method in space (Taylor-Hood elements) and a backward difference formula of order 2 in time with a semi-implicit (Picard) linearization. The solver is implemented in the C++ Object Oriented library lifeV . We carried out simulations on Stampede2 from the extreme science and engineering discovery environment (XSEDE). For all the simulations conducted in the study, we assume blood to be a Newtonian incompressible flow with density and viscosity . The aortic wall is assumed to be rigid. Time-step, grid size (and grid independence) and the number of heartbeats run to have a reasonable periodic flow are illustrated and discussed in Refs.  and .
The study is carried out on (1) an idealized aortic arch and (2) a patient-specific aorta with a degenerated AAA, i.e., with a degradation of the aortic wall.
2.2.1 The Idealized Aortic Arch.
As a proof of concept, we first investigate the impact of the filter radius δ and inflow rate Q(t) in a simplified aortic arch (Fig. 1). The geometry is composed of a half torus and two cylinders. The detailed dimensions are in Fig. 1. In this case, we do not test the uncertainty for the geometry. The stochasticity of the inflow and the filter radius is detailed below.
2.2.2 The Patient-Specific Abdominal Aortic Aneurysm.
To parametrize realistic variations of the arterial geometry, we considered a patient affected by the AAA, available from the iCardioCloud Project supported by the Cariplo Foundation, Italy (No. 2013–1779). The patient had a significant abdominal aortic dilation from 2010 to 2016, as shown in Fig. 2. The figure reports the 3D reconstruction of the morphology from the available CT scans, performed with the vascular modeling ToolKit . These geometries provided the basis for the construction of the stochastic map needed by the PCE, as explained in Sec. 2.3.1.
2.3 Stochastic Representation of the Inputs.
To use the PCE, we associate a stochastic process to the inputs of interest. This avoids the need for many samples, requiring many data currently unavailable. To this aim, we postulate some a priori stochastic information on the inputs. The criterion for these arbitrary choices is the consistency with the physical context.
We list the specific stochastic features of the random variable for each input in Table 1.
2.3.1 Stochastic Representation of the Aneurysmal Geometry.
To the best of our knowledge, it is the first time that a sensitivity analysis on geometry is conducted by using a longitudinal dataset and a registration method. In this way we embrace a large range of possible realizations of the geometry, on a time scale of 4 years. The same procedure can be used on different time scales and with more than two scans, adjusting the scaling function (12) accordingly.
2.3.2 Stochastic Representation of the Inflow Rate.
In the absence of patient-specific inflow data, we designed a stochastic inflow by scaling the average inflow waveform q(t) reported in Fig. 3 (retrieved from Ref. ) by a normally distributed univariate random variable , where and . This results in a range () of peak inflow being from 203 to 380 ml/s, covering the flowrate from at rest to stress .
2.3.3 Stochastic Representation of the Filter Radius δ.
The filter radius δ is modeled as a normally distributed univariate random variable with mean being the minimal mesh size and standard deviation being . As we mentioned in the introduction, this parameter is empirically tuned as a tradeoff between stabilization (preferring a larger value) and accuracy (suggesting a smaller value). The standard deviation 25% of the “nominal” value is a reasonable range in our practical experience.
2.4 Model Responses (Outputs).
The model responses considered here are the TKE, TAWSS, and OSI. These are indices with potential relevance for clinical applications in aortic diseases (see, e.g., Refs.  and [36–38]).
where Ω represents the whole computational domain.
OSI of values around indicates a high occurrence of retrograde flows.
The mean and variance of these QoIs are computed as well as the corresponding Sobol' indices to quantify the relative contribution of different uncertain inputs to the uncertainty of QoIs.
2.5 Recall on Global Sensitivity Analysis.
2.5.1 The Sobol' Index.
where f0 is a constant and denotes all the possible combinations of s elements of the pool of indices .
, where δAB is the Kronecker delta ( for ; otherwise).
which are the global sensitivity indices indicating the functional structure of the model  related to the combination of the input parameters . Notice that fA = 0 if and only if VA = 0.
2.5.2 Polynomial Chaos Expansion Based Sobol' Indices.
Traditionally, the Sobol' indices are computed using Monte Carlo (MC) or other sampling methods . To be accurate, MC approaches generally require many samples, since the accuracy improves with only the square root of the number of samples. In our case, each sample requires the numerical simulation of the LES model, with a prohibitive computational cost. Therefore, we resort to the truncated PCE for a surrogate model to propagate uncertainties and compute the Sobol' indices. In this way, the computational cost reduces to that of estimating the coefficients of the PCE. Specifically, the convergence of the Sobol's indices from PCE depends on: (1) the smoothness of the output of interest, (2) the number of dimensions of the inputs under consideration [8,44]. To balance the computational cost and accuracy, we test here truncated PCE of degree 2 and 3 (see Sec. 3).
2.5.3 Truncated Polynomial Chaos Expansion.
Any second-order random process (i.e., with a finite variance) can be represented as a series of polynomials in random inputs [45,46]. Denote again by a vector of independent random variables with a joint distribution defined on . The corresponding image probability space of the random variables is , where Γ is the image of is the Borel σ-algebra on Γ and is the probability measure defined on . Consider
as the space of all polynomials in ξi with degree up to p;
as the set of all polynomials in orthogonal to with respect to the probability measure ;
as the space of polynomials spanned by .
is conventionally taken as . Due to the orthogonality, .
2.5.4 Polynomial Chaos Expansion-Based Sobol' Indices.
where is the deterministic coefficient.
where for . As shown in Eq. (36), the cost of computing the Sobol' indices corresponds to the cost of evaluating the coefficients of the truncated PCE.
2.6 Computing the Coefficients of the Polynomial Expansion.
There are different approaches to accomplishing this task. “Intrusive approaches” (i.e., Galerkin approach [40,48]) and “nonintrusive” (point collocation  and pseudo-spectral projection) methods [49–51] are used to compute the stochastic modes for the PCE. In the former, the governing equations are reformulated to target the PCE. In the latter, available deterministic solvers are combined with collocation or projection techniques. Compressed sensing is usually adopted with the point collocation method to enforce the sparsity of the coefficients .
where and wr are quadrature points and associated weights. The number of required numerical simulations (ultimately, the computational complexity of the method) and the accuracy depend on the number and the collocation of the quadrature nodes. We adopted nested sparse grid quadrature points based on Leja sequence  and generated by the open-source tool Chaospy . In this way, the number of forward simulations required for evaluating the stochastic modes for d = 2 and a PCE of degree 2 is 15 while it is 21 for PCE of degree 3, as shown in Fig. 4. Here, the number of quadrature points are subjected to an “empirical rule” as detailed in Ref. , which is implemented in Chaospy. This is much less than the number of samples required by MC approaches. Further potential improvement can be made by adopting adaptive sparse PCE [54,55], which is subject to future work. The degree of the surrogate PCE is cross-validated by comparing its results to the results of the original LES model.
3.1 Results on the Idealized Aortic Arch.
As a proof of concept, the impact of the filter radius and inflow rate are investigated in a simplified aortic arch (Fig. 1). Details of the random input parameters are listed in first two rows of Table 1.
The truncated PCE approximations of the TKE, TAWSS are computed with polynomial basis of degree p = 2 and p = 3. As anticipated, in this case, d = 2, we needed 15 and 21 simulations, respectively. The differences between the results from the PCE of different degrees are negligible.
The mean, variance, and prediction interval of the TKE for the simplified aortic are shown in Fig. 5. The Sobol' indices for the TKE is shown in Fig. 6 (left).
The mean, variance, and prediction interval for TAWSS are reported in Fig. 7, while the Sobol' indices for TAWSS is shown in Fig. 6 (right).
3.2 Results on the Patient-Specific Abdominal Aortic Aneurysm.
The truncated PCE approximations of the TKE, TAWSS, and OSI are computed with PCE of degree p = 2. We compute the stochastic modes using the nonintrusive pseudo-spectral methods with the Leja nested sparse grid quadrature points. The resulting number of simulations required is 35 (d = 3 inputs in this case). To be concrete, the realizations of the geometry in the five different values of ξG are illustrated in Fig. 8.
The mean and coefficient of variation for the TAWSS (top row) and OSI (bottom row) are reported in Fig. 9. The Sobol' indices for the TKE, TAWSS, and OSI are reported in Figs. 10 and 11, respectively.
4.1 The Idealized Case.
The variance of the TKE caused by the variation of the filter radius δ during the systole (from 0 to 0.3 s) is minimal, while the influence of the filter radius is evident during the diastole, as shown in Fig. 5 (left). However, the inflow impact is much higher than that caused by δ, as shown in Fig. 5 (right).
4.1.1 Sobol' Indices for the Total Kinetic Energy.
The relative influence of the variations of the filter radius δ as well as the inflow boundary conditions on the hemodynamic factors of this simplified aortic arch are investigated using the Sobol' indices. The Sobol' index of δ in time (Fig. 6 left) represents the contribution of δ to the TKE variation through the cardiac cycle. Consistent with the variance results (in Fig. 5), this Sobol' index is much lower than that from the inflow rate during systole, despite the higher coefficient of variation (CoV) of δ. This result demonstrates the robustness of the model during systole. However, the influence of δ exceeds that of the inflow rate during diastole.
The TAWSS Sobol' index of δ is lower than that of the inflow in most of the arch, except for the inflow and the interior bend. The higher Sobol' index of δ in these regions might be due to the influence of the δ during the diastolic phase.
4.2 The Patient-Specific Abdominal Aortic Aneurysm Case.
TAWSS and OSI are often computed together to assess the nature of the flow in the neighborhood of the wall, and the connection to possible focal diseases or degenerations of the wall. As expected, the mean TAWSS is found to be lower in the region with more growth than other regions on the aneurysm. Correspondingly, we found that OSI is relatively high in the region with more growth. These two hemodynamics indicators point out a possible disturbed oscillating nature of the flow in the region with growth.
More interestingly, the CoV for TAWSS is higher in the region with more growth (i.e., region with relatively low mean TAWSS), while the opposite result is observed from OSI. This means that the reliability of the two hemodynamic indices is somehow complementary in terms of space-dependent reliability. Their joint computation is justified by the fact that one is more reliable where the other is less.
4.2.1 Sobol' Indices for the Total Kinetic Energy.
The Sobol' indices of the filter radius δ, the inflow rate, as well as the abdominal aneurysmal geometry for the TKE are in Fig. 10. The result for the relative contribution from the filter radius and inflow rate is consistent with the finding from the simplified aortic arch, i.e., the inflow rate has relatively more contribution to the total variation during systole. In contrast, the filter radius is more influential than the inflow rate during diastole. However, the variation of the abdominal aneurysmal geometry of the aorta has a significantly dominant role for the variability of the TKE throughout the whole cardiac cycle. The partial variance caused by the variation in the geometry is larger than those caused by the other two factors combined, as shown in Fig. 10.
4.2.2 Sobol' Indices for the Time-Average Wall Shear Stress and Oscillatory Shear Index.
The Sobol' indices of the filter radius δ, inflow rate, as well as the abdominal aneurysmal geometry for TAWSS and OSI are in Fig. 11. Similarly to what we found for the TKE, the relative contribution of the three investigated factors (ranked from high to low) to both TAWSS and OSI are the aneurysm geometry, the inflow rate, and the filter radius. The TAWSS of the two extreme quadrature points for ξG is in Fig. 12 and the OSI is in Fig. 13. The filter radius and inflow rate are frozen at their mean values, respectively. We can appreciate that the TAWSS and OSI may change significantly even when the geometry is the only source of variability.
The work bears some limitations, as listed below.
Due to the constraint of computational costs, only three sources of uncertainty are considered simultaneously.
We assume simplified models for the variations in the inflow boundary condition for the lack of patient-specific data of a large population.
In the simulations, the arterial wall is rigid, and the flow is Newtonian.
Limitations (1) and (2) will be addressed in future works, while we speculate that the limitations in (3) have a minor impact on the conclusions of this study. We argue that the qualitative conclusions drawn from the results here, particularly in terms of prioritizing the impact of the different uncertainties hold despite these limitations.
We performed a global sensitivity analysis using PCE-based Sobol' indices to investigate the impact of the filter radius of the deconvolution-based Leray model, inflow rate, as well as the geometry on the hemodynamic simulations in aortas. The conclusions of the results are summarized hereafter.
Arterial morphology has always been speculated to be a critical factor affecting the computational hemodynamic predictions. The results in this study confirm this intuition: the geometry is the most influential uncertain input, comparing to the other inputs considered. Other studies report similar results on idealized geometries. In Ref. , the authors considered idealized AAA and carotid sinus geometries  with the radius of the aneurysm as an uncertain input. Thanks to our novel use of the image-registration procedure, coupled with the truncated PCE analysis, we can extend the results to patient-specific geometries with a realistic geometrical variability. We observed significant variations of wall shear stress and axial velocity fields with the variation of the geometry. For completeness, we mention that in coronary arteries, the relative importance of the minimal luminal diameter was shown to exceed that of the viscosity and boundary resistance in computing FFRCT . Our results lead to a similar conclusion in the aorta. In this respect, it is crucial to use patient-specific geometries and to quantify the associated sensitivity to predict clinically relevant results reliably.
Minimal influence of δ on clinically relevant quantities is shown by the Sobol' index of δ for TAWSS and OSI in the patient-specific AAA case, as shown in Fig. 11. These numerical results demonstrated that the clinical QoIs are quite insensitive to the filter radius and therefore we can guarantee the robustness of applying the Leray model in aortic simulations.
The third practical conclusion we want to highlight from our results is that different hemodynamic indices may show a different level of reliability in space. The CoV results of TAWSS and OSI in the patient-specific AAA case of this study seem to “complement each other,” in the sense that one is more trustworthy in a specific region where the other is less reliable and vice versa. This result suggests that it is worth computing both in practice and complementing mutually the information they provide to have a complete, reliable picture of the clinical condition of the patient, so to use one or the other QoIs in different regions of interest.
The educated (possibly automatic) combination of the different QoIs in clinical scenarios and the extension of the present global analysis to more uncertain inputs will be the subject of future works. Specifically, we target the sensitivity analysis on the aortic syndrome of type B aortic dissections in our following work and its impact on the clinical routine of diagnosis of this specific disease.
We gratefully acknowledge the help of Romarowski Rodrigo (San Donato Hospital, Milan, IT, and University of Pavia, IT) for the patient-specific CT images.
We carried out simulations on this study on the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation (Grant No. TG-ASC160069; Funder ID: 10.13039/100000001).
H. Xu and A. Veneziani acknowledge the support of the U.S. National Science project NSF-DMS 1620406/1620384 “Collaborative Research: Efficient Modeling of Incompressible Fluid Dynamics at Moderate Reynolds Numbers” (Funder ID: 10.13039/100000001).