The impact of shape uncertainty on aortic-valve pressure-drop computations.

Patient-specific image-based computational fluid dynamics (CFD) is widely adopted in the cardiovascular research community to study hemodynamics, and will become increasingly important for personalized medicine. However, segmentation of the flow domain is not exact and geometric uncertainty can be expected which propagates through the computational model, leading to uncertainty in model output. Seventy-four aortic-valves were segmented from computed tomography images at peak systole. Statistical shape modeling was used to obtain an approximate parameterization of the original segmentations. This parameterization was used to train a meta-model that related the first five shape mode coefficients and flowrate to the CFD-computed transvalvular pressure-drop. Consequently, shape uncertainty in the order of 0.5 and 1.0 mm was emulated by introducing uncertainty in the shape mode coefficients. A global variance-based sensitivity analysis was performed to quantify output uncertainty and to determine relative importance of the shape modes. The first shape mode captured the opening/closing behavior of the valve and uncertainty in this mode coefficient accounted for more than 90% of the output variance. However, sensitivity to shape uncertainty is patient-specific, and the relative importance of the fourth shape mode coefficient tended to increase with increases in valvular area. These results show that geometric uncertainty in the order of image voxel size may lead to substantial uncertainty in CFD-computed transvalvular pressure-drops. Moreover, this illustrates that it is essential to assess the impact of geometric uncertainty on model output, and that this should be thoroughly quantified for applications that wish to use image-based CFD models.

INTRODUCTION

Imaging modalities such as computed tomography (CT), magnetic resonance imaging (MRI), and ultrasound are used on a daily basis in the clinic. These imaging modalities are traditionally used to extract diagnostic information directly, for example, jet‐velocity or aortic valve area for patients with suspected aortic valve stenosis. , , In addition to deriving diagnostic metrics directly, these images may be used to inform image‐based computational models, an approach that has become increasingly popular over recent years. With segmentation tools, the complex three‐dimensional (3D) patient‐specific geometries are readily obtained from these images. Such patient‐specific geometries have been extensively used as input to computational fluid dynamics (CFD) models, and can provide detailed pressure and velocity fields in the blood vessel. Applications of image‐based CFD models include, amongst others: the estimation of the transvalvular pressure‐drop , ; estimation of the significance of coronary lesions , ; quantification of wall shear stresses in carotid arteries ; hemodynamic evaluation of arteriovenous dialysis grafts ; or estimating energy losses in intracranial aneurysms. However, segmentation of the patient‐specific geometry is not exact, and uncertainty in the segmented shape is inevitable. Previous CFD studies have demonstrated that uncertainty in shape may substantially affect the outcome of any subsequent modeling. , , Segmentation quality is largely determined by image quality, but may also be a consequence of the segmentation method itself. The general procedure to evaluate the performance of a (semi) automatic segmentation method is to compare it to a ground‐truth segmentation; that is, a segmentation that is performed or corrected manually by a (clinical) expert. Though these uncertainties in shape will propagate to the output of the CFD model and can have an enormous impact on model‐based decisions, output uncertainties are rarely considered, let alone quantified.

Previous studies have demonstrated that differences between automatic and manual segmentation depend on image quality and the geometry of interest. For instance, Ecabert et al. demonstrated that segmentation from CT images of the heart and surrounding vessels with a shape‐constrained deformable model may introduce errors of 0.50–0.82 mm for the heart chambers, and 0.60–1.32 mm for part of the great vessels, when compared to manual segmentations. With a similar approach for aortic valve tract segmentation from MRI images, Queirós et al. observed average errors of 0.54 0.08 mm. The thin and flexible nature of heart valve leaflets arguably makes segmentation of these constructs more difficult. Nevertheless, various authors have proposed algorithms for segmenting heart valves from CT or echocardiography images. Pouch et al. used a deformable model to semi‐automatically segment the aortic valve from 3D transesophageal echocardiography data, and reported mean and 95th percentile errors of 0.5 0.1 mm and 1.0 0.2 mm compared to manual segmentations. Similarly, performed a validation study of a valve segmentation tool on a data set consisting of 1516 4‐D transesophageal and 690 cardiac 4‐D CT images. Segmentation errors of 1.54 1.17 and 1.36 0.93 mm were reported for the aortic and mitral valves when using transesophageal echocardiography and CT images respectively. Using machine learning, Liang et al. reported mean discrepancies of 0.69 0.13 mm between automatic and manual aortic valve segmentation from CT images. These contributions demonstrate that segmentation errors are typically in the order of image resolution, that is, roughly between 0.5 and 1.0 mm with the current state‐of‐the‐art hard‐ and software. Although seemingly small, this inherent geometric uncertainty may substantially affect any subsequent modeling of the detailed flow field.

A small number of studies have demonstrated that uncertainty in shape may indeed substantially affect model output. Sen and colleagues compared three segmentation methods with ground‐truth segmentations of 45 intracranial aneurysms, and found that, on average, CFD‐computed energy losses and wall shear stresses differed by 23.2% 8.7% and 24% 8.5%, respectively. Likewise, studies that aimed to assess coronary stenosis severity identified that geometric uncertainty (lumen diameter) may contribute substantially to uncertainty in hemodynamic simulations. , , , In addition, the work by Sankaran et al. recognized the value of coronary lumen uncertainty, highlighting regions in the coronary tree where simulation output was sensitive to lumen uncertainty. These studies demonstrate that shape uncertainty may be a major source of error in image‐based patient‐specific computational models, and the impact of shape uncertainty on model output should thus be thoroughly assessed in order to strengthen confidence in such computational models. Moreover, such analyses are indispensable and an integrated part of the market approval process when aiming for clinical applications that rely in image‐based CFD models. ,

In previous work we developed a computational framework to compute the transvalvular pressure‐drop from CT images by making use of CFD. In addition, we combined statistical shape modeling, meta‐modeling, and CFD in order to obtain a cheap‐to‐evaluate meta‐model, which related changes in shape to changes in the transvalvular pressure‐drop versus flow relationship. It was shown, that with a limited number of shape modes, physically relevant geometric variation in the population was captured. However, uncertainties in valvular shape, and consequently in simulation output, were not considered. In this study, we build upon this previous work, and quantified uncertainty of the transvalvular pressure‐drop versus flow relation to uncertainty in valvular shape. Uncertainty in valvular shape was imposed by assuming uncertainty in the shape mode coefficients. In addition, we performed a variance‐based sensitivity analysis to apportion the output uncertainty to uncertainties in the input.

This manuscript is outlined as follows: first a brief introduction to variance‐based sensitivity analysis and uncertainty quantification is presented. Second, the available data, statistical shape model, simulation framework, and meta‐model are discussed. Third, specific details on how sensitivity analysis and uncertainty quantification was applied is presented. The results of the sensitivity analysis and uncertainty quantification are summarized in the results section. Finally, the discussion puts the findings in an academic and clinical perspective.

METHODS

Scope of study

The aortic valve sits between the left ventricle and the ascending aorta, and ideally these surrounding structures should be included in image‐based computational models. Moreover, flow through the aortic valve is pulsatile in nature and as such also impacts the dynamic behavior of the aortic valve. This study focused on geometric uncertainty of the aortic valve, and as such made several simplifications with respect to the physiological condition. Hence, to make the study of the impact of shape uncertainty tractable, three simplifications were considered for modeling transvalvular flow: (1) the aortic valve was isolated from its surrounding structures: the shape of the aorta and ventricle were not considered; (2) steady‐state flow was assumed; and (3) the aortic valve construct was considered rigid.

Overview sensitivity and uncertainty quantification

When reducing output uncertainty of a computational model, sensitivity analysis and uncertainty quantification can be used to assess which parameters contribute most to output uncertainty, and are thus most rewarding to measure or estimate accurately. More specifically, a computational model with stochastic input factors can be represented by an arbitrary function , which yields an uncertain output , that is: Typically uncertainty in the input parameters are non‐negligible, and can considerably affect the uncertainty in model output . Here, a single output variable is considered, but in a similar way, a multivariable output can be analyzed. Sensitivity analysis and uncertainty quantification can be used to quantify the effect of input uncertainties on model output. For instance, by using the variance in as a measure for output uncertainty. Global variance‐based sensitivity analysis can separate direct contributions of each individual uncertain input from interactions between different model inputs, for example, interactions between and , where represents all parameters except . In variance based sensitivity analysis, the contribution of the uncertain parameters to the variance of can be expressed as the main () and total () Soboĺ sensitivity indices : where and are operators for the variance, and expected value, respectively. The main indices are measures for the expected reduction in variance when would be known without uncertainty. Besides the main effects, the total sensitivity index includes interaction as well, and reflects the contribution of all terms which include . A low total sensitivity index indicates that this parameter may be fixed within its uncertainty domain. Sensitivity indices may be computed by methods such as Monte Carlo or adaptive sparse generalized Polynomial Chaos Expansion (agPCE). ,

Uncertainties in boundary conditions and shape can have a substantial influence on the output uncertainty of a three‐dimensional computational model (Equation 1), for example, due to segmentation errors. In traditional “forward‐engineering” CAD models, shape variations are easily introduced by changing the shape‐defining parameters, such as angles, diameters, ratios, or thickness. However, 3D patient‐specific computational models are generally of complex shape, difficult to parameterize, and the influence of shape variation is therefore mostly neglected. Instead of using physically meaningful parameters, this work used statistical shape modeling to parameterize the shape of the valve. The statistical shape modes provided a parameterization of the geometry, and facilitated the training of a meta‐model (Figure 1(A)). The meta‐model was trained on the output of CFD simulations, and related variations in shape to variations in simulation output. Consequently, this cheap‐to‐evaluate meta‐model was used in a sensitivity analysis and uncertainty quantification framework (Figure 1(B)) for two purposes: (1) to evaluate the importance of shape variation on a population level, and (2) to evaluate the importance of shape variation on the level of an individual patient, and how shape uncertainty affects model output uncertainty for that patient. Sensitivity analysis and uncertainty quantification was applied to aortic valve pressure‐drop versus flow relations, computed with a CFD model. The following sections elaborate on the individual components presented in Figure 1.

FIGURE 1

Schematic of the procedure for meta‐model training (A). An inner meta‐model is trained with simulation data, and yields a surrogate function that relates the simulation input (without uncertainties) to simulation output . Note that the meta‐model function is an approximation of the original function . (B) The cheap‐to‐evaluate surrogate function can be used for sensitivity analysis and uncertainty quantification to determine how uncertain inputs propagate to output uncertainty

Inner meta‐model training

Data acquisition and valve segmentation

Computed tomography imaging datasets of 74 patients were available, and corresponded to the dataset that was used in Hoeijmakers et al. Images were acquired with an in‐plane spatial resolution of 0.31–0.68 mm, and slice thickness of 0.34–0.70 mm. Images were acquired with electrocardiography gated CT, and segmentations represented the peak‐systolic state of the aortic valve. Consequently, aortic valve segmentation was performed with a Shape Constrained Deformable model framework, see, for example, previous work by Ecabert et al. , or Weese et al. In the segmentation framework, a template mesh model was iteratively adapted to an image. From the resulting segmentation mesh, a submesh () was extracted that consisted of the left ventricular outflow tract, aortic valve, sinuses, and part of the ascending aorta (Figure 2). The submesh consisted of 1808 vertices and 4223 triangular faces, and had a consistent topology (). This resulted in mesh correspondence between patients. Hence, any surface mesh was completely described by the coordinate vector :

FIGURE 2

Visualization of the first five shape modes. Note that the mean mesh is the same for each realization. captures variations in the open/closed state of the aortic valve. elongates the entire valve and its surrounding structures. seems to skew the geometry orthogonal to the centerline. mainly affects the diameter of the sinotubular junction. affects sinus size, and angle between left ventricular outflow tract and the ascending aorta

Statistical shape model

Training of the inner meta‐model (Figure 1(A)) is considerably easier when a compact parameterization of the shape of the aortic valve is available. In this study we used a statistical shape model to obtain such a compact parameterization. Segmented aortic valves were aligned by a generalized Procrustes analysis, which optimally translated, rotated, and scaled each of the meshes by minimizing the sum of squared errors between the mean and the target mesh. Consequently, statistical shape modeling was used to extract the main directions of shape variance from the 74 segmented aortic valves by performing an eigen‐decomposition of the co‐variance matrix (principal component analysis), see previous work by Hoeijmakers et al. The eigenvectors of the co‐variance matrix, typically referred to as shape modes, describe the main directions of shape variance. The eigenvalues of this eigen‐decomposition describe the amount of shape variance that was captured with the corresponding eigenvector.

Using the statistical shape model, any shape in the training set was then approximated by the mean shape of the training set () plus a linear combination of a small number of shape modes (), which were weighted by . That is: with the number of shapes that were available in the training set, and the number of shape modes used for the approximation. The optimal set of weights () for each patient's aortic valve was found by minimizing the difference between the original segmentation and the approximation : The introduction of allowed to control the weight of each vertex to the minimization problem, and was chosen such that the 216 vertices that were part of, or adjacent to the free cusp edges weighted five times stronger than all other vertices. This weighting factor was empirically established, and helped to reduce the approximation error of vertices that were part of, or close to the free cusp edges. For a more detailed description on the performance (e.g., compactness and generalizability) of the statistical shape model and for a comparison between CFD results on the original segmentation and the approximated shape by the SSM we would like to refer to the work of Hoeijmakers et al.

The shape coefficients can be regarded as a patient‐specific parameterization of the aortic valve and surrounding structures, and values for the 74 original valve segmentations were found by Equation (6). Besides approximating the surface meshes in the training set, Equation (5) was used to generate a population of “virtual” aortic valve geometries that were within three times the square root of the eigenvalue () of its respective shape mode, also see Hoeijmakers et al. or Heimann and Meinzer. Figure 2 illustrates the shape variation that was contained in the first five shape modes, and showed that the first shape mode (), weighted by , controlled the opening and closing behavior of the aortic valve. The second shape mode () stretches the aortic valve construct in axial direction. The third shape mode () seemed to skew/shear the valve orthogonal to the centerline. The fourth shape mode () mainly changed the diameter of the sinotubular junction. Finally, the fifth shape mode () mainly affected the size of the sinuses, and the angle between the left ventricular outflow tract and ascending aorta. Shape variation was extracted by principal component analysis, and as such, the shape modes were ordered from high to low explained shape variance, that is, ordered with respect to statistical relevance. The first five shape modes captured 61% of the shape variance. However, statistical relevant shape variation may not correspond to physically meaningful shape variation. Hence, a computational fluid dynamics workflow was developed that allowed the computation of the transvalvular pressure‐drop, given the boundary conditions (e.g., flow‐rate) and valvular shape (any valve realization of Equation 5).

Simulation workflow

An automated CFD workflow was developed for computing the transvalvular pressure‐drop at specified flowrates. To facilitate CFD modeling, the outflow boundary was extended by five diameters. Additionally, the volume enclosed by the surface mesh (Figure 2) was discretized by approximately polyhedral elements (ANSYS Fluent R18.2, ANSYS Inc., Canonsburg, PA). Edge lengths of the polyhedrals were chosen based on a mesh‐sensitivity study, and ranged between 0.15 and 2.0 mm. Blood was modeled as an incompressible Newtonian fluid with a density of 1060 kg/m3 and dynamic viscosity of 4 mPa s. At the inflow‐boundary (left ventricular outflow tract), a plug‐velocity profile was prescribed that corresponded to steady‐state volumetric flow‐rates between 50 and 650 mL/s. Pressure at the ascending aorta was set to zero, and no‐slip boundary conditions were assumed at the walls. Reynolds numbers at the inflow boundary depended on flowrate/shape combination, but were estimated to be between 600 and 13,000, hence a shear stress transport model (5% turbulent intensity at the inflow boundary) was used to model turbulence. The governing equations were solved with ANSYS Fluent R18.2 (ANSYS Inc., Canonsburg, PA) by making use of the semi‐implicit method for pressure linked equations. Exemplary results of the velocity field and pressure evolution along the centerline are shown in the results section, that is, in Figure 6(B),(C). The transvalvular pressure‐drop was defined as the difference between inflow pressure and outflow pressure (see also Figure 6(C)). We would like to note that this approach may lead inaccuracies for healthy cases where pressure losses caused by the valve are small. For these cases friction losses in the extended flow section may become non‐negligible. Resampled surface meshes of the 74 segmented aortic valves can be found in an online data repository. Moreover, this repository also contains the centerline definitions and exemplary centerline pressures at a flow rate of 400 mL/s.

Inner meta‐model

A meta‐model was trained on the simulation input parameters (without uncertainty) and the corresponding output parameter (the transvalvular pressure‐drop). Selecting the most suitable meta‐model, and meta‐model settings is often difficult since no universal meta‐model exists that performs well for all problems. However, it has been shown that ensemble approaches, where a weighted sum of meta‐models is considered, can yield a good approximation. , That is, the goal of an ensemble‐type meta‐model is to obtain the best weighted‐average of a selection of meta‐models: where is the prediction at of the final meta‐model, that is, the weighted ensemble of various meta‐models and their settings. To find the optimal combination of meta‐models a penalized predictive score was proposed by Ben Salem and Tomaso. This score combined three components (Equation 8): (a) optimizing the internal accuracy by evaluating the mean square error on training samples/points (); (b) use a 10‐fold cross‐validation to evaluate predictive capability on unseen samples (); and (c) minimize over‐fitting of the meta‐model by a thin‐plate spline Bending Energy Functional (). , The penalized predictive score was then constructed by weighting the contribution of each of these components: In this work the Genetic‐Aggregation meta‐model of Ben Salem and Tomaso was used to relate the shape mode coefficients (, , , , ), a global scaling parameter (), and the volumetric flow rate () to the CFD‐computed transvalvular pressure‐drop.

The quality of the meta‐model (Equation 7) should converge when the number of training points increases. Hence, the seven‐dimensional (five shape parameters, scaling, and flowrate) input space was uniformly sampled with Latin Hypercube designs (maximin) of 25, 50, 100, 200, 400, 800, 1600, and 3200 samples. Samples were excluded when the transvalvular pressure‐drop exceeded 300 mmHg, when the aortic valve was completely closed, or when simulations diverged. On average 26% of the simulation samples were excluded based on these criteria. The input‐space for the shape parameters were limited to lie within . The scaling parameter was limited to values between 0.8 and 1.25 (observed in the set of 74 valve segmentations), and volumetric flow‐rate between 50 and 650 mL/s. Consequently, the Genetic‐Aggregation meta‐model was trained on the resulting simulation data. In the remainder of this manuscript, this meta‐model will be referred to as the inner meta‐model (Figure 1(A)).

Quality of the inner meta‐model

Quality of this inner meta‐model was evaluated by the root mean square error, relative root mean square error, and the mean absolute percent error. Figure 3 demonstrates that errors for the verification samples, which were excluded from training, reduced considerably with an increase in the available training samples. Additionally, Figure 3(B),(C) suggest that beyond 1000 samples the quality of the meta‐model levels off. The meta‐model trained with the most training samples yielded the best quality, and was used as the inner meta‐model which facilitated sensitivity analysis and uncertainty quantification (Figure 1(B)).

FIGURE 3

Inner meta‐model quality as function of the number of successfully simulated training samples. An increase in available training samples improves the root mean square error (RMSE), relative root mean square error (rRMSE) and the mean absolute percent error (MAPE) between the meta‐model and verification points which were not used for training the meta‐model

Sensitivity analysis and uncertainty quantification

Changes in boundary conditions and shape affect the CFD‐computed transvalvular pressure drop. In order to quantify the sensitivity of the transvalvular pressure‐drop to changes in shape, main and total sensitivity indices (see Equations 2 and 3) were obtained by the agPCE method. , The method uses polynomial chaos expansion in order to obtain an explicit formulation of (Equation 1), and expands the stochastic model output into a series of orthogonal polynomials. The generalized polynomial chaos expansion is defined as: where represents the polynomials and represents the expansion coefficients. An adaptive algorithm (agPCE) was used to include only the polynomials that significantly increased meta‐model quality. This kept meta‐model quality high, while keeping the required training set small. Sensitivity indices were analytically derived from this expansion. , Quality of this “outer” meta‐model was evaluated with a leave‐one‐out cross‐validation coefficient (). ranges between 0 (lowest quality) and 1 (highest quality), and it has been shown that with increasing sensitivity indices convergence. For this study, Legendre polynomials—which are most suitable for uniform distributions—up to order 3 were considered.

Population‐based sensitivity

A global variance‐based sensitivity analysis was performed in order to determine the sensitivity of the transvalvular pressure‐drop to changes in shape mode coefficients and flow rate. The lower and upper limits of each of the shape mode coefficients were restricted to and , respectively. This range covered a wide range of feasible shapes in the population, whilst avoiding the corners of the input domain where meta‐model results were likely poor. Uncertainties in flow‐rate were based on the work by Namasivayam et al., who measured the average systolic flow rate in 1131 patients with severe or moderate aortic stenosis, and found mean systolic flow rates of 243 42 mL/s. It was assumed that systolic flow approximately follows a sine squared profile, hence mean‐systolic values were multiplied by two to obtain peak systolic flow. This resulted in a peak‐systolic lower limit of 382 mL/s and upper limit of 590 mL/s in flow uncertainty. The distribution of the input uncertainties were unknown, hence a uniform distribution was assumed for all input parameters, which can be considered a conservative estimate of the uncertainty distribution.

Patient‐based sensitivity analysis and uncertainty quantification

In addition to a global variance‐based population sensitivity analysis, an additional sensitivity analysis was performed for each patient separately. Recall that each individual segmentation was approximated by Equation (5) (). With variations in these shape‐mode coefficients, uncertainty in shape can be emulated. That is, the original segmentations were not exact, and some uncertainty with respect to the actual in‐vivo geometry is expected.

Previous studies that used deformable model‐based segmentation frameworks have demonstrated that segmentation errors are in the order of voxel size. , Hence, it is assumed that local uncertainty in shape is between 0.5 and 1.0 mm. This shape uncertainty may be imposed by introducing variations in the shape mode coefficients. That is, any change in leads to (localized) changes in vertex positions in a particular direction (Figure 4). The maximum () and mean () rate of change with respect to was obtained from the mean mesh and shape modes numerically (using Equation 5), and depends on the specific shape mode, see Table 1 and Figure 4. Uncertainty in shape mode coefficients were calibrated to impose a maximum () vertex‐to‐vertex displacement of 0.5 and 1.0 mm, which were considered in separate analyses. Figure 4 specifies shape uncertainty per shape mode when imposing an uncertainty of 0.5 mm. This resulted in patient‐specific uncertainty ranges, and is in Table 1 expressed as a percentage of the total feasible range (). Only shape changes were of interest, and volumetric flow‐rate was fixed and uncertainty in the scaling parameter was neglected. The uncertainties in Table 1 were applied to each patient, and using the agPCE algorithm, the patient‐specific sensitivity indices were computed. Exploratory Monte‐Carlo simulations (all cases, 0.5 mm uncertainty, flowrate 400 mL/s, 500,000 samples per case) demonstrated that the distribution was different for each case, where stenotic cases generally showed more skewed output distributions. Hence, boxplots were used to visualize uncertainty in the transvalvular pressure‐drop versus flow relation. Boxplots (min, max, 25, 50, and 75th percentiles) were constructed by making use of the samples (typically around 40) that were used to build the agPCE meta‐model of highest . Due to the low‐discrepancy Soboĺ sequence sampling of the uncertain input space around each patient, it was assumed that boxplots of the agPCE samples were representative of the expected output distribution.

FIGURE 4

Resulting uncertainty in vertex‐to‐vertex position (mapped onto ) by introducing a 0.5 mm uncertainty by considering uncertainty in the shape mode coefficients (Table 1). Each shape mode introduces variation in a specific part of the surface model (first and second column) and in a particular direction (third column). The histogram illustrates how uncertainty in vertex position was distributed over all 1808 vertices for each shape mode

TABLE 1

Patient‐based parameter uncertainties for sensitivity analysis and uncertainty quantification

Parameter	Symbol	∂ε_max∂α_m mm/α †	∂ε¯∂α_m mm/α ‡	Minimum	Maximum	Assumed change in α_m
Shape mode	α_1	85	17	α_p,1−0.018⋅6λ_1	α_p,1+0.018⋅6λ_1	3.6% of 6λ_1
Shape mode	α_2	76	21	α_p,2−0.025⋅6λ_2	α_p,2+0.025⋅6λ_2	5.0% of 6λ_2
Shape mode	α_3	58	22	α_p,3−0.039⋅6λ_3	α_p,3+0.039⋅6λ_3	7.8% of 6λ_3
Shape mode	α_4	79	21	α_p,4−0.035⋅6λ_4	α_p,4+0.035⋅6λ_4	7.0% of 6λ_4
Shape mode	α_5	84	21	α_p,5−0.035⋅6λ_5	α_p,5+0.035⋅6λ_5	7.0% of 6λ_5

Note: Note that with these parameters any shape variation is defined with respect to the patient‐specific shape found by Equation (5), and these values of uncertainty were calibrated to emulate a maximum () shape variation of 0.5 mm. mm was obtained by multiplying these values by 2.

Partial derivative indicating the maximum displacement that is observed with a unit‐change in that particular shape mode coefficient. The region where vertex displacement is maximum is illustrated—per shape mode—in Figure 4.

Partial derivative indicating the mean displacement of all vertices in the mesh with a change in that particular shape mode coefficient.

RESULTS

Main and total sensitivity indices and of the population‐based sensitivity analysis are depicted in Table 2. The outer agPCE meta‐model with highest quality reached a of 0.9998, which can be considered excellent quality. Results of this high‐quality outer meta‐model are summarized in Table 2, and suggest that the transvalvular pressure‐drop is most sensitive to changes in the first shape mode (open/closure of the aortic valve). 93% of the output variance was explained by variations in this input parameter (, ). Contribution of shape modes (, ), (, ), and (, ) to the output variance was limited to around 1%–2%. Although shape mode was ranked higher statistically compared to , variation in did not contribute to the output variance (, ). Hence, shape variation in the direction of shape mode seemed unimportant from a population perspective. Flow‐rate is another contributing factor, but only accounted for 1.6% of the total output variance (, ). Total sensitivity indices suggest that limited interaction was present between shape modes and/or flow‐rate. Additionally, the total sensitivity index of remained approximately zero.

TABLE 2

Main and total sensitivity indices of the global population‐based sensitivity analysis

Parameter	Symbol	S_i †	S_i Range ‡	S_T,i †	S_T,i Range ‡
Shape mode	α_1	0.930	[0.900–0.937]	0.953	[0.900–0.955]
Shape mode	α_2	0.010	[0.009–0.012]	0.016	[0.010–0.018]
Shape mode	α_3	0.011	[0.003–0.011]	0.018	[0.003–0.020]
Shape mode	α_4	0.000	[0.000–0.004]	0.001	[0.001–0.004]
Shape mode	α_5	0.010	[0.008–0.017]	0.018	[0.012–0.020]
Flow‐rate	Q	0.016	[0.015–0.067]	0.018	[0.018–0.067]

Sensitivity index obtained from the outer meta‐model (agPCE) with highest ().

Range of sensitivity indices considering all generated meta‐models irrespective of .

The global population‐based sensitivity analysis suggests that does not contribute to the output variance, and may thus be fixed. This is partially supported by the results of the patient‐based sensitivity analysis. Patient‐based sensitivity indices (in order of increasing aortic valve area) are illustrated in Figure 5. Uncertainty in the first shape mode coefficient accounted for more than 90% of the output variance for severely and moderately stenotic heart valves (Figure 5). Additionally, Figure 5 suggests that the remaining output variance for these two subgroups was mostly explained by uncertainty in and . Similar to the observations of the population‐based sensitivity analysis, uncertainty in seemed to play a minor role for severely and moderately stenotic valves. Interestingly, main sensitivity indices seem to strongly depend on valvular area. That is, importance of decreased, whereas importance of tended to increase for more open valves, and actually even exceeded that of and (e.g., see main indices of cases 69: and 70: ).

FIGURE 5

Sensitivity indices per case, given the uncertainty in the shape mode coefficients of Table 1. Results are sorted from small to large geometric valve area (bottom graph). With an increase in valvular area, variation in other shape modes start to have a substantial effect on the transvalvular pressure‐drop. Classification according to the basic grading criteria. Severe: valve area < 100 mm2; moderate: valve area 100–150 mm2; mild/healthy: valve area >150 mm2

Besides the main and total sensitivity indices, samples that were used to construct the agPCE model were used to construct boxplots around the meta‐model results. These boxplots provide and estimate for uncertainty in the transvalvular pressure‐drop versus flow relation for each case. In Figure 6(A) the transvalvular pressure‐drop vs. flow curves of a typical stenotic, moderate, and mild/healthy case are enhanced with these boxplots, given an imposed uncertainty of 0.5 and 1.0 mm (Table 1). Figure 6(A) exposes that the uncertainty in output was considerable when uncertainties of 0.5 and 1.0 mm were imposed. For example, for Case 7 (severe aortic valve stenosis; valve area 86 mm2), the inner meta‐model gives a transvalvular pressure‐drop of 69 mmHg at 400 mL/s. However, interquartile ranges () were 74–63 (median: 68, min–max: 55–84) mmHg and 82–58 (median: 68, min–max: 44–105) mmHg for imposed uncertainties of 0.5 and 1.0 mm, respectively. Case 38 (moderately stenotic), the inner meta‐model yields a transvalvular pressure‐drop of 25 mmHg with interquartile ranges of 26–23 (median: 25, min–max: 21–29) mmHg and 28–21 (median: 24, 17–34) mmHg, for imposed uncertainties of 0.5 and 1.0 mm, respectively. Finally, for Case 73 (healthy), uncertainty seems to deteriorate in a relative sense. That is, the inner‐meta model yielded a transvalvular pressure‐drop of 0.6 mmHg with interquartile ranges of 0.6–0.5 mmHg (median: 0.6, min–max: 0.4–0.7), and 0.7–0.4 (median: 0.5, min–max 0.2–0.9) mmHg for imposed uncertainties of 0.5 and 1.0 mm respectfully. Transvalvular pressure‐drop versus flow relations and uncertainties for all cases can be found in Figures A1, A2, A3 of the Appendix. These results clearly demonstrate that with small deviations in geometry, for example, in the direction of the first shape mode, the uncertainty in the output of the model can become substantial.

FIGURE 6

(A) Typical results for a severely stenotic case (case 07—left column), a moderately stenotic case (case 38—middle column), and a mildly stenotic case (case 73—right column). (A) Transvalvular pressure‐drop versus flow curves, filled circles represent the CFD results on the original segmentation mesh . Results of the inner meta‐model (dashed line with triangles) are augmented with boxplots. Boxplots of imposed uncertainties of 0.5 mm (blue) and 1.0 mm (red) are depicted. (B) flow field at 400 mL/s. (C) Corresponding pressure evolution over the centerline, distance normalized by outflow diameter. Pressure was averaged over cross‐sections perpendicular to the centerline

FIGURE A1

Volumetric flow‐rate versus transvalvular pressure‐drop curves for subgroup A: valve opening area <100 mm2. Filled circles represent the CFD results on the original segmentation mesh . Results of the inner meta‐model (dashed line with triangles) are augmented with boxplots. Boxplots of imposed uncertainties of 0.5 mm (blue) and 1.0 mm (red) are depicted

FIGURE A2

Volumetric flow‐rate versus transvalvular pressure‐drop curves for subgroup B: valve opening area 100–150 mm2. Results of the inner meta‐model (dashed line with triangles) are augmented with boxplots. Boxplots of imposed uncertainties of 0.5 mm (blue) and 1.0 mm (red) are depicted

FIGURE A3

Volumetric flow‐rate versus transvalvular pressure‐drop curves for subgroup C: valve opening area >150 mm2. Results of the inner meta‐model (dashed line with triangles) are augmented with boxplots. Boxplots of imposed uncertainties of 0.5 mm (blue) and 1.0 mm (red) are depicted

DISCUSSION

The aim of this study was to quantify how sensitive transvalvular pressure‐drop computations were to uncertainty in valvular shape. Two sensitivity analyses were done. First, a global variance‐based population‐level sensitivity analysis was performed, and the main and total sensitivity indices were extracted. This revealed that uncertainty in the weighting of was most important, accounting for 93% of the expected variance. The main indices of , , and were low, suggesting limited importance of uncertainty in the direction of shape modes , , and . Similarly, the contribution of volumetric flow rate to output uncertainty was low and may be fixed in the range that was considered. We would like to note that may however become important again when uncertainty in valvular shape is substantially reduced. In addition, the total sensitivity indices of suggest that geometric changes in the direction of shape mode were unimportant.

Secondly, uncertainties in patient‐specific valvular shape were imposed by considering uncertainty in shape mode coefficients. In this approach, we introduced an uncertainty in shape mode coefficients that corresponds to likely segmentation errors of up to 0.5 and 1.0 mm. With this we were able to emulate geometric uncertainties, and study how this propagated through to patient‐specific pressure‐drop vs. flow relationship. The results of this sensitivity analysis partly corroborate the population‐based sensitivity analysis. However, main sensitivity indices from the patient‐based sensitivity analysis demonstrate that the importance of shape variation is in fact patient‐specific. More specifically, for severely stenotic and moderately stenotic cases, sensitivity indices roughly correspond to the indices found with the population‐based analysis, that is, uncertainty in seems most important (). However, Figure 5 also shows that geometric uncertainty in the direction of the other shape modes tended to become more important with an increase in valvular area. The results from the population‐based sensitivity may therefore not be representative for the entire spectrum of valve configurations.

Shape mode captures valve opening and closing, and naturally would have the most substantial effect on the computed transvalvular pressure‐drop. That is, a change in the weighting of leads to a change in valvular area, and as a consequence the predicted transvalvular pressure‐drop changes. Additionally, it was shown that errors in the direction of which are in the order of voxel size can lead to substantial uncertainty in the transvalvular pressure‐drop. This suggest that accurate segmentation of the free cusp edge of the aortic valve could substantially reduce uncertainty in the computed transvalvular pressure‐drop versus flow relationship.

The direction of and seem to correspond with changes in ascending aorta diameter. The observation that the (recovered) transvalvular pressure‐drop becomes more sensitive to changes in and with an increase in valvular area may be explained by pressure‐recovery. That is, when blood is accelerated into the narrow orifice, pressure decreases. Consequently, downstream from the aortic valve flow decelerates again, and pressure is (partly) recovered. , , , When the cross sectional area of the ascending aorta is large, more kinetic energy will be converted back to pressure. When the aortic valve is fully open (e.g., when healthy) the net transvalvular pressure‐drop is small. As a consequence, the relative importance of uncertainties in ascending aorta diameter ( and ) will increase compared to moderately or severely stenotic cases. Nevertheless, uncertainty in is still most important, and accounts for at least 60% of output variance (e.g., see Case 69). Additionally, it is observed that for severely and moderately stenotic valves and generally exceed . These observations indicate that geometrical variation that is statistically relevant, may not necessarily be relevant from a fluid mechanics perspective. This is supported by the work of Wu et al., who used statistical shape modeling to introduce global deformation modes to emulate airfoil geometric uncertainty. It was shown that deformation (shape) modes of lower statistical relevance can in fact be more important for transonic aerodynamic lift and drag coefficients, than modes of higher statistical relevance.

Available literature shows that statistical shape modeling is a technique that has been applied for numerous applications. Some examples in the bio‐medical field include, organ segmentation , , and extraction of morphology from medical images , , Only a limited number of studies have tried to combine statistical shape modeling with computational models, and are mainly found for orthopedic applications. For example, statistical shape models have been used for real‐time prediction of knee joint mechanics, predicting femur bone strength, or creating parametric models to model cervical spine loading. The study by Bredbenner et al. introduced uncertainty in geometry considering cervical spines that were from the mean shape, and showed that shape variation influenced the computed axial, flexion‐extension, and lateral displacements. However, sensitivity indices were not explicitly computed, and only variations from the mean shape were considered. Khalafvand and colleagues studied the influence of shape variation on intraventricular flow variables, such as wall shear stress, vortex formation time, and the time integral of energy dissipation, but did not compute sensitivity indices of the simulation output parameters. In this study we have demonstrated that besides the accepted applications, statistical shape modeling in combination with computational modeling can also be used to determine physically relevant shape variation. Moreover, we have shown that this approach may be used to estimate how shape uncertainty affects uncertainty in the output of the computational model.

Image‐based CFD models often include complex anatomical shapes. Some examples of computational models of complex shapes include: the aortic tree, coronary tree, or lungs. Similar to valve geometries, these shapes are difficult to parameterize, and as such investigating how shape uncertainty affects the results is challenging. Sankaran et al. resolved this by splitting the coronary tree into sections, and exploring the solution in a family of probable geometries with an assumed uncertainty in radius of 0.3 mm. Similar to our approach, Sankaran et al. used a surrogate (meta) model based on machine learning, to replace the compute‐intensive blood flow simulations. This demonstrates that meta‐models, which learn the relation between physics and shape, will be crucial for uncertainty quantification.

Limitations

The Shape Constrained Deformable Model framework is yet to be validated for aortic valve segmentations, and imposed uncertainties were therefore hypothetical. Such a validation study of the segmentation framework could include manual corrections, and may yield further insight where performance of automatic segmentation is poor. Nevertheless, uncertainties of 0.5 and 1.0 mm seem plausible based on existing literature. , , , , The shape modes inherently depend on the data‐set. As a consequence, precise control over local shape variation was not possible with the method that was proposed. Furthermore, our dataset was restricted to aortic valves and their direct surroundings, and ascending aorta shape—and variation therein—was not captured by the analysis. The patient‐specific shape of the ascending aorta may have affected hemodynamics and pressure losses. For instance the work by Ha et al. suggests that the angle of the aortic‐valve jet, which is a combination of valve shape (captured in this study) and ascending aorta shape (not captured in this study), may affect the rotational direction and strength of helical blood flow in the aorta. Hence, future work should consider including ascending aorta shape (variation) to assess how much this affects uncertainty in transvalvular pressure‐drop computations. Additionally, patient‐specific flow profiles were not available, and a plug‐flow assumption was assumed. Previous work suggests that this may lead to large relative errors below 10 mmHg in particular. , These limitations would however not change our conclusion that geometric uncertainty in the order of voxel size may lead to considerable uncertainties in simulation output, and as such should be thoroughly quantified for image‐based computational models.

The inner meta‐model relates variations in shape modes (weighted by the shape coefficients) and flow rate to the transvalvular pressure‐drop. A small number of shape modes was used, and yielded an approximation of the original segmentation. Likewise, the inner meta‐model, relating shape coefficients to the transvalvular pressure‐drop, is an approximation of the transvalvular pressure versus flow relationship as well. Figures A1, A2, A3 demonstrate that for most cases the meta‐model adequately approximates this relationship. However, for some cases in Subgroup C (healthy), this method seems to break down. That is, training of the inner meta‐model seemed insufficient, or the five shape modes did not seem to adequately capture the features that were relevant for the transvalvular pressure‐drop versus flow relationship.

The agPCE method is a variance‐based sensitivity method. This assumes that variance can fully capture the uncertainty in the output parameter. However, this may be inappropriate when the output distribution is skewed or multi‐modal. Hence, the sensitivity indices (Figure 5) may not be appropriate for all cases. The exploratory Monte‐Carlo simulations with the uncertainty ranges in Table 1 demonstrated that for healthy cases the distribution in transvalvular pressure‐drop followed a normal distribution, and variance is an appropriate statistical measure. However, for moderately stenotic cases, the distribution became lightly left‐skewed. Furthermore, it was observed that left‐skewness (tendency towards a lower transvalvular pressure‐drop) increased with an increase in stenosis severity. Hence, variance‐based methods for sensitivity analysis may not be ideal in all circumstances. Density based methods which characterize the output distribution by the cumulative density function, may in those cases be more appropriate, and should be considered in future work.

In line with Hoeijmakers et al., the transvalvular pressure‐drop at peak systole was approximated by CFD simulations, and expressed as an uncertain scalar parameter . In‐vivo however, the transvalvular pressure‐drop would strongly vary throughout the cardiac cycle due to flow pulsatility. Hence, the current approach could be expanded to include uncertainty in the transvalvular pressure‐drop over time. This would require computationally more demanding pulsatile simulations. Moreover, would need to be expanded to an uncertain vector that captures temporal variation.

Previous work by Marlevi et al. suggested that the recovered transvalvular pressure‐drop is driven by turbulent dissipation. In this study this turbulent dissipation was modeled through a Reynolds‐averaged Navier–Stokes‐type turbulence model which assumes fully developed isotropic turbulence, and which may not be physiologically realistic. However, we would like to note that this choice was a practical consideration due to the limited computational cost of such models compared to more accurate—but substantially more demanding—scale resolving methods such as direct numerical simulation or large eddy simulations.

Conclusion

We have developed a method for sensitivity analysis and uncertainty quantification for transvalvular pressure‐drop versus flow relationships. This method assumes that physically relevant shape variation can be adequately captured with a statistical shape model. Consequently, a meta‐model that is trained on a limited number of these shape modes can be used to quantify how (local) geometric uncertainty affects the transvalvular pressure‐drop computations. With this method we have demonstrated that geometric uncertainties in the order of voxel size may in fact strongly influence transvalvular pressure‐drop predictions by image‐based CFD. Hence, it is indispensable to quantify the impact of shape uncertainty on model results for applications that rely on image‐based CFD models.