Strongly Coupled Morphological Features of Aortic Aneurysms Drive Intraluminal Thrombus

Over 75% of abdominal aortic aneurysms harbor an intraluminal thrombus, and increasing evidence suggests that biologically active thrombus contributes to the natural history of these potentially lethal lesions. Thrombus formation depends on the local hemodynamics, which in turn depends on morphological features of the aneurysm and near vasculature. We previously presented a hemodynamically motivated “thrombus formation potential” that predicts where and when thrombus might form. Herein, we combine detailed studies of the three-dimensional hemodynamics with methods of sparse grid collocation and interpolation via kriging to examine roles of five key morphological features of aneurysms on thrombus formation: lesion diameter, axial position, length, curvature, and renal artery position. Computational simulations suggest that maximum diameter is a key determinant of thrombogenicity, but other morphological features modulate this dependence. More distally located lesions tend to have a higher thrombus formation potential and shorter lesions tend to have a higher potential than longer lesions, given the same aneurysmal dilatation. Finally, movement of vortical structures through the infrarenal aorta and lesion can significantly affect thrombogenicity. Formation of intraluminal thrombus within an evolving abdominal aortic aneurysm thus depends on coupled morphological features, not all intuitive, and computational simulations can be useful for predicting thrombogenesis.

of activated platelets alone appears to initiate thrombus formation on an undamaged endothelium 5 . Herein, we use the TFP to evaluate how different combinations of five key morphological features -maximum lesion diameter, axial position of the lesion within the infrarenal aorta, projected lesion length, relative renal artery positions, and aortic curvature -contribute to the potential formation of a thrombus within an AAA. Toward this end, we introduce an analytical descriptor of idealized lesion geometries to facilitate systematic variations of the five morphological features and then perform 179 full unsteady 3-D hemodynamic simulations for those lesions that are identified as potentially revealing by an adaptive sparse grid collocation method. Then, using kriging, we interpolate results across this five-dimensional morphometric space to estimate thrombus formation more generally, thereby effectively extending significantly the number of lesions studied without necessitating the many additional computational simulations that would have been needed. In this way, one can more easily investigate parameter sensitivity and correlations among morphometric parameters, and thereby begin to build intuition otherwise not possible.

Results
Hemodynamic Predictions. Figure 1 shows our simple yet broad parameterization scheme using five key geometric variables g i (i = 1 … 5) to characterize the morphology of each AAA: (1) the maximum diameter of the aneurysm (g 1 ) is that value within the largest cross-section of the lesion orthogonal to the centerline, (2) the axial position of the aneurysm (g 2 ) is the vertical distance from the iliac bifurcation to the geometric center of the largest cross-section, (3) the length of the aneurysm (g 3 ) is the Euclidean distance between the proximal and distal entrances into the lesion, which are identified by a radius of the infrarenal abdominal aorta greater than 110% of the baseline value, (4) the spinal bending magnitude (SBM or g 4 ) of the aneurysm is the maximum distance between the centerline of the lesion and the vertical axis of the aorta at the level of the iliac bifurcation, and (5) the renal artery offset (RAO or g 5 ) is the signed distance between the ostia of the right and left renal arteries. These five parameters can be assessed from standard clinical imaging (e.g., computed tomography or magnetic resonance angiography) and are often considered by vascular surgeons when evaluating treatment options. An associated analytical representation (see Methods) enabled a straightforward generation of myriad idealized aneurysms suitable for systematic evaluation of the consequences of individual or combined morphological features. Figure 2, panel A, shows a prototypical lesion geometry as well as morphologies that resulted from the extreme ranges of the parameters chosen (Table 1), noting that our focus was on thrombus formation and thus lesions early in their evolution. Shown, too, are representative "points" from our five-dimensional parameter space collapsed onto a three-dimensional space (panel B) to visualize multiple levels (0-3) of the collocation scheme that was used to identify those lesion geometries of most interest (see Methods). Finally, for further ease of visualization, we show select idealized geometries by superimposing a particular lesion geometry on the associated collocation point in two of the many possible two-dimensional projections (panel C). Overall, this adaptive collocation scheme suggested that we examine 179 individual lesions to evaluate how diverse hemodynamic and thrombus-related metrics varied as a function of changes in the five key morphological features. Figure 3, panel A, illustrates key hemodynamic characteristics (e.g., distributions of wall shear stress (WSS), particle shear history, and vortical structures) via corresponding scalar metrics for a representative axisymmetric AAA (g 1 = 4.0 cm, g 2 = 3.0 cm, g 3 = 7.5 cm, g 4 = 0.0 cm, g 5 = 0.75 cm) that has a high TFP (>3 within the lesion). Some features of the distributions of WSS-based indices seen in this AAA are common to most aneurysms. As expected, luminal enlargement reduces the time-averaged wall shear stress (TAWSS), particularly in central and distal regions of the lesion (e.g., 0.09 Pa when averaged over the central and distal lesion vs. 0.45 Pa when averaged over the healthy portion of the infrarenal aorta); see the first column in panel A. The oscillatory shear index (OSI) was high within the lesion in regions of low TAWSS (OSI = 0.33 averaged over the central and distal aspects) as well as in the healthy segment of the infrarenal aorta (OSI = 0.27) where the renal split induced oscillatory (i.e., backflow) flow during diastole (second column). Retrograde flow occurs in the infrarenal aorta in both healthy subjects and AAA patients, and is thought to favor atherogenesis and significantly influence AAA hemodynamics. Not surprisingly then, an Endothelial Cell Activation Potential (ECAP) index, which combines OSI and TAWSS after opportune scaling (see Methods), highlighted central and distal portions of the lesion as more likely to include susceptible endothelial cells (third column). Particle tracking revealed lower values of a platelet activation potential (PLAP) index, which accounts for the flow-induced shear history accumulated by blood-borne particles (see Methods), in the healthy portion of the infrarenal aorta (average value of 0.21), but higher values of PLAP throughout the AAA (average of 0.64), particularly on the sides of the proximal and distal necks of the lesion (average of 0.71 in those regions), as seen in the fourth column. These results suggest that lateral portions of the necks might collect flowing particles that have accumulated a significant shear history during their recent hemodynamic path, and may thereby convey activated platelets to the aneurysmal wall. The interplay between near-wall hemodynamics and particle advection is recapitulated by the distribution of TFP (rightmost column). The upper 99 th percentile of this index (denoted TFP-99p) for this lesion was 4.3, a value well above the thrombogenic threshold of 2.5 to 3.0 suggested by our previous patient-specific modeling that compared thrombogenic AAAs to non-thrombogenic carotid arteries and healthy infrarenal aortas 5 . TFP was high throughout the central part of the aneurysm, where ECAP was highest, and also laterally at the AAA necks, where PLAP was highest and ECAP indicated a potentially susceptible endothelium. In contrast, TFP was low everywhere else, suggesting that the combined occurrence of low and oscillatory WSS plus delivery of platelets having experienced a high shear history is peculiar to AAA hemodynamics, which renders the TFP index distinct among the metrics considered.
A more detailed interpretation of the distribution of multiple hemodynamic indices can be achieved by analyzing vortical structures through plots of negative λ 2 iso-contours 7 . Figure 3, panel B, reveals a distinct hairpin vortical structure (VS) that originated during diastole at the start of retrograde flow within the infrarenal aorta. This structure was fully formed at about 80% of the time into diastole (i.e., 0.75 s into the 0.9 s cycle), but broke-up just before the end of the cardiac cycle. Small remnants of the larger VS remained throughout systole of the following beat, during which they were advected to the iliac bifurcation. While snapshots of systole are not shown directly, these events are discernable through iso-contour plots of the minimum λ 2 field, which collects traces left by VSs throughout a cardiac cycle (see the rightmost column). Our simulations suggested that the position and geometry of VSs as they break-up differ significantly for different AAA geometries (not shown), which could help explain hemodynamic mechanisms underlying different values of TFP for different lesions. For the lesion highlighted in Fig. 3, break-up of the largest VS occurred proximally at about one-third of the aneurysm length (hence VS Depth = 30%), and its area at break-up was 0.69 cm 2 , which corresponded to ~6% of the largest AAA cross-sectional area. As expected for an axisymmetric lesion with a zero SBM (i.e., g 4 = 0.0 cm), the center of the VS aligned perfectly with the centerline of the lesion, thus VS Bending was 0.
In contrast, Fig. 4 shows five lesions (selected from the 179 full simulations) that exhibited some unexpected values of the thrombogenic indices. Together, these AAAs illustrate, in part, why hemodynamic simulations can be helpful in predicting AAA thrombogenicity, particularly in cases that cannot otherwise be intuited or explained by simple correlations to individual morphometric features of an aneurysm. The first two rows of Fig. 4 show 4-cm diameter lesions (i.e., g 1 = 4 cm) that are positioned centrally within the infrarenal aorta (g 2 = 6 cm) and have the maximum allowed length (g 3 = 12 cm). The first lesion, AAA 33 (1 st row), had a mid-range value of Figure 1. Parallel display of a reconstructed patient-specific AAA (left) and an idealized representation (right) that is defined geometrically by 5 constants (c i , i = 1, 2, …, 5); see text for the associated equations. Importantly, both lesions are characterized by 5 morphometric features: maximum diameter (g 1 ), projected distance from the iliac bifurcation to the maximum diameter (g 2 ), overall projected axial length (g 3 ), a spinal bending magnitude (SBM or g 4 ), and a renal artery offset (RAO or g 5 ). A cross-section is shown for both geometries at the maximum diameter (top right) to highlight that parameter g 1 is a mean value. The curvature and length of an idealized axisymmetric AAA (middle right) can be defined as a function of three constants (c 1 , c 2 , c 3 ); for example, c 2 defines the length of the lesion (g 3 ). As seen in the lower right, the spinal bending magnitude (g 4 ) is generated by two additional constants (c 4 and c 5 ). See text. spinal bending (g 4 = 2.5 cm) but a maximum value of renal artery offset (g 5 = 2.0 cm); we thus denote this lesion by = .
. g [4, 6, 12, 2 5, 2 0] 33 . This lesion was selected because it showed the highest upper percentile values of PLAP and TFP (4.79 and 8.55, respectively) among all idealized aneurysms considered. For reference, mean  values of ECAP, PLAP, and TFP were 3.98, 1.01, and 3.46, respectively, for all lesions having a 4.0 cm diameter; the proposed TFP threshold of 2.5 to 3.0 would thus suggest that 4.0-cm diameter lesions would be weakly thrombogenic on average. Moreover, given the low diameter-to-length ratio (4:12), one might not have expected high thrombogenicity in AAA 33 despite its large lesion volume, and given its slight tortuosity one might have expected flow stagnation (and therefore higher TFP) antero-laterally. Yet, the peculiar coupling of morphologic features in this AAA generated pronounced flow mixing and a heterogenous ensemble of VSs of varying sizes, most of which broke-up just below the axial center (VS Depth = 51%). As for all lesions with a forward SBM (g 4 = 2.5 cm), the VSs shifted towards the back of the aneurysm close to the posterior wall (VS Bending = 0.77 cm). ECAP was thus higher towards the front of the aneurysm, and high PLAP manifested across the lesion wall, resulting in high TFP regions on the anterior proximal neck, anterior center of the AAA, and lateral and posterior portions of the distal neck. AAA 107 (2 nd row) was fully axisymmetric (g 4 = 0 cm) and had a small RAO (g 5 = −0.5 cm), but otherwise the same maximum diameter, axial position, and overall length as AAA 33; namely, = − . g [4, 6, 12, 0, 0 5] 107 . The lack of tortuosity suggests that the only asymmetrical bias would arise due to the RAO, which was minimal. This different combination of features again resulted in pronounced mixing, with high ECAP occurring only at the anterior center of the lesion while PLAP was high almost everywhere within the central region of the lesion. The largest VSs broke-up proximally (VS Depth = 29.5% of aneurysm length), releasing high shear history particles that were advected to central and distal portions of the lesion (i.e. where PLAP was highest) during the following systole. TFP was again high on the anterior surface and, with a TFP-99p of 6.95, this lesion had high thrombogenic potential. . It is thus a much shorter lesion (g 3 = 3 cm) though with the same diameter (g 1 = 4 cm) and a more proximal location (g 2 = 7.1 cm) with a moderate SBM (g 4 = 2.5 cm). This lesion has a small volume, so one might anticipate a low TFP. Due to the tortuosity, one might expect flow stagnation (and therefore higher TFP) to occur antero-laterally. Due to the high axial position and immediately proximal branching vessels, one might expect disturbed flows and therefore high TAWSS and low TFP. AAA 143 showed high ECAP values, however, driving high thrombogenicity (upper percentile values of 5.74 and 5.07 for ECAP and TFP, respectively). Minimum λ 2 iso-contours (rightmost column of 3 rd row) suggest that a high aspect ratio can induce VSs that do not deviate significantly from their pre-lesion path, resulting in lower TAWSS. In this case, spinal bending (SBM) exacerbated this effect on the anterior portion of the lesion, which showed the highest ECAP and PLAP, with a TFP-99p well above the expected threshold.
Finally, the last two rows of Fig. 4 show counter-intuitive examples for two larger AAAs (g 1 = 4.0 cm and 5.0 cm, respectively, 4 th and 5 th rows) with thrombogenic indices that are either lower than or on the lower end of our proposed threshold for TFP. The fourth row depicts AAA 144, = .
. . g [4, 7 06, 7 5, 0, 0 75] 144 , a more proximally-located lesion having a standard diameter, length, and RAO, but no tortuosity. This lesion thus has a typical volume and aspect ratio, which might naively suggest a more "typical" TFP presentation; recall that the average TFP-99p for lesions with a 4.0-cm diameter and 7.5-cm length was 3.10, so one might expect this lesion  Fig. 3). These cases were selected given the a priori non-intuitive findings. See text for details but note that [g 1 , g 2 , g 3 , g 4  to be thrombogenic. The lack of tortuosity implies that the only asymmetric bias would be due to RAO, which was mid-range, so one might expect biasing of high TFP to the left side. Yet, this lesion highlights how a proximal axial positioning can lower TFP, a recurring finding of our parametric study. In particular, AAAs located closer to the renal bifurcation tend to be exposed to higher shear stress (as expected near a flow discontinuity), and can be traversed by heterogeneous vortical structures that often maintain their integrity to the distal neck of the lesion. Indeed, for AAA 143 (3 rd row), VS break-up was close to the distal neck (VS Depth = 80%), well beyond the high ECAP islet on the front side of the proximal neck. Not surprisingly, then, with little overlap of ECAP and PLAP, AAA 144 (4 th row) had an upper percentile TFP of only 2.0, below our proposed threshold for thrombogenicity.
The lesion in the 5 th row, AAA 175, = .
. . . g [5, 7 0607, 7 5, 2 5, 0 75] 175 , had the same axial position as AAA 144, but maximum diameter and spinal bending within the range considered. From diameter alone, one might expect this lesion to be highly thrombogenic (high TFP-99p); due to the high aspect ratio, one might also expect greater flow stagnation and thus higher TFP. Yet, the more proximal axial position lowered ECAP and thus TFP, and because of spinal bending, the VSs tended to remain in the posterior half of the lesion (rightmost column of the 5 th row) before breaking-up at about 50% of the aneurysm length. Plots of TFP suggest that particles released by the breaking-up of a VS may find a susceptible endothelium on the anterior distal neck, where ECAP was highest (with an upper percentile value of 3.25). The relatively far distance between the high ECAP region and VS break-up could explain a TFP-99p of only 2.53, which was lowest among all lesions with a maximum diameter of 5.0 cm. Taken together, then, Fig. 4 suggests again that it need not be a single dominant morphological feature that dictates where and when a thrombus will form. Statistical Inference. Although it is useful to consider individual cases, as in large patient-specific data sets, it can be more revealing to systematically assess the influence of each geometric parameter on thrombogenicity while determining possible coupling amongst the parameters. Figure 5, panel A, depicts combinatorial maps of TFP-99p for metrics of geometry, hemodynamics, and vortical structures. Kriging was used to interpolate TFP-99p results from the 179 full simulations over the five-dimensional space Ξ formed by g i over the specified ranges (see Methods and Supplement). Panel A shows 2D slices of Ξ for TFP-99p resulting from multiple combinations of two of the g i , with all other parameters fixed at their median values within each slice. Results from the 179 full simulations that resided in that plane are depicted by open circles. Panel B shows similar maps, except among diameter (g 1 ), axial position (g 2 ), VS Bending , VS Area , and VS Depth . These parameters are most responsible, in descending order (40%, 30%, 10%, 10%, 10%, respectively), for variability in TFP-99p, as determined by ANOVA of the 179 simulations. Panel 5B does not show slices of Ξ, but rather similar slices of another parametric space, Ξ', composed of this particular set of parameters. All 179 points are projected onto each slice in this case.
The first column of Fig. 5A relates lesion diameter (g 1 -abscissa) to the other four morphological features (g 2−5 -ordinates). Effects of maximum diameter tend to combine with the axial position (g 2 -1 st row) and length (g 3 -2 nd row) of a lesion to increase thrombogenic potential (i.e., TFP-99p > 2.5), but less so with SBM (g 4 -3 rd row) and RAO (g 5 -5 th row). Consider first the first row/first column (maximum diameter and axial position). While TFP-99p tends to increase monotonically with diameter, the gradient increases as the lesion location becomes more distal. Importantly, this coupling again suggests that more proximal lesions may be protected from thrombogenesis if the maximum diameter is 5 cm or less, with smaller lesions having a more distal location also appearing to be less susceptible to thrombogenesis. Results in the second row/first column show combined effects of maximum diameter and axial length, which are highly nonlinear with several regions of interest. TFP-99p again tends to increase monotonically with diameter, but the effect of length on this gradient is not monotonic. One region of interest is the lower right of the map. These lesions have a large diameter and short length (high aspect ratio) and are extremely thrombogenic, tending to have TFP-99p > 5. Indeed, AAA 143, depicted in the 3 rd row of Fig. 4, falls within this region. A second region of interest is the upper right of the map. These lesions have a large diameter and length (large volume) and tend to be thrombogenic, though not as thrombogenic as those with a high aspect ratio. Interestingly, there is a valley between these two regions where the Gaussian curvature is negative, approaching zero (recall that these maps are projections onto a plane for easier visualization). This implies that for large diameters, particular aspect ratios may be more protective against thrombus, especially at certain axial positions. This behavior is not seen at small diameters. The third row/first column shows coupling between maximum diameter and SBM, which yields minimal to modest effects on TFP-99p, hence this map primarily shows the effect of diameter. The fourth row/first column shows coupling between maximum diameter and RAO, which again yields small TFP-99p in many regions though slightly higher thrombogenicity at large diameters and large renal artery offsets. Interestingly, the gradient with diameter appears most right-shifted when RAO is zero and the geometry is bilaterally symmetric, suggesting that this state may be more protective against thrombus.
The second column of Fig. 5A relates effects of axial position (g 2 -abscissa) with three other g i (ordinate) metrics (see the first row/first column for its relation to maximum diameter). As noted before, more proximal lesions (i.e., those a greater distance from the iliac bifurcation) tend to have lower thrombogenicity, and below a certain axial position, the behavior and coupling of other parameters appear to dominate. The second row/second column shows coupling between axial position and lesion length, which is highly nonlinear. Short lesions (g 3 < 5) have high values of TFP-99p that appear to be relatively insensitive to axial position, whereas a saddle point appears to develop at moderate axial position and length, again implying that a certain aspect ratio (diameter is fixed at 4.0 cm in this slice) may be protective against thrombus. Longer lesions tend to be thrombogenic only if located more distally. The third row/second column shows coupling between axial position and SBM, which appears to be minimal, with a slight rightward shift of the axial position TFP-99p gradient for lesions with no tortuosity. The fourth row/second column shows coupling between axial position and RAO, which is not large. Here, TFP-99p appears to decrease monotonically with axial position, though with a slight rightward shift in the gradient at low and high RAO.
The third and fourth columns show remaining couplings, namely between lesion length and either SMB or RAO and between SBM and RAO (all at a maximum diameter of 4.0 cm). Generally, lesion length appears to couple nonlinearly with the other g i (all fourth-row panels), though this coupling is weaker with SBM and RAO. Interestingly, longer or shorter lesions can be more thrombogenic while those of intermediate length can be less so. The saddle points within the length maps are non-intuitive, but aspect ratio and volume, which impact ECAP, play roles. Certain aspect ratios may be protective at higher diameters. Finally, the fourth row/fourth column shows the effect of Spinal Bending Magnitude (g 4 ) on RAO, which is nonlinear but generally weak. The lower left region of the map contains lesions that have no tortuosity and a more distal left renal artery origin (and are close to being bilaterally symmetric). These lesions appear to be more thrombogenic, which may seem counter to the idea that both symmetry and distance from proximal branching vessels protect against thrombus formation. Figure 5B shows similar maps, though focusing on vortical structures. The first column shows the coupling of diameter (abscissa) with axial position as well as three metrics of vortical structure (VS): VS Bending , VS Area , and VS Depth . Again, maximum diameter is a dominant morphological feature, though modulated by vortical structures. The first row/first column shows the coupling between diameter and axial position, which is similar to that in Fig. 5A except that the interpolations and mappings are within a different space, Ξ'. The first column/second row shows the coupling between maximum diameter and VS Bending . Recall that VS Bending represents the deviation of the center of the vortical structure from the centerline of the vessel, which, due to the anterior deflection of SBM, is always posterior. Similar results are in the third and fourth rows/first column, except for VS Area and VS Depth . Low values of these metrics of the vortical structures combined with large lesion diameters result in high values of TFP-99p. It seems, for example, that high VS Bending , that is, deflection of the VS from the centerline towards the walls, is mildly protective. Recalling that VS Area represents the cross-sectional area of the VS orthogonal to the centerline at the point of break-up, VS Area appears to be insignificant at low values of maximum diameter but important at high values. When VS Area is large, the VS is close to the wall and decreases ECAP. Beyond a VS Area of 1.0 cm 2 , the diameter of the lesion no longer appears to be relevant and the lesion is protected from thrombus. Finally, recall that VS Depth refers to the distance the VS traverses within the lesion prior to break-up. Results (fourth row/first column) suggest that the VS passes through the lesion at high VS Depth , not breaking up and not releasing shear-activated platelets to regions of high ECAP, thus resulting in low TFP. Indeed, diameter becomes irrelevant at high VS Depth and the gradient is intensified at low VS Depth .
The second column shows couplings between axial position and vortical structures. As it can be seen, lower values of the VS metrics tend to associate with a higher potential for thrombosis. When assessing VS Area and VS Depth in terms of axial position of the lesion, we also see local TFP maxima close to (but not coinciding with) the bottom left corner of the maps, thus more distal lesions are slightly more vulnerable to thrombosis than proximal lesions given similar values of the VS metrics. This finding may be because distal lesions tend to exhibit high ECAP regions, especially close to the iliac bifurcation. Activated platelets that are released proximally (at minimal VS Depth ) or farthest from the aneurysmal wall (at minimal VS Area ) may need to travel a relatively far distance to reach susceptible regions, thus resulting in lower TFP. Finally, the third and fourth columns show coupled effects among different metrics of the vortical structures. These also exhibit local maxima and saddle points, but changes are minimal despite the complex topography. Again, lower values of VS Bending , VS Area , and VS Depth are more thrombogenic, consistent with other maps. Figure 6 provides another means of visualizing possible correlations amongst the different geometric and hemodynamic metrics, including the five key geometric parameters g i , lesion aspect ratio (a derived metric), the three VS metrics, and finally ECAP, PLAP, and TFP. Panel A shows results for all simulations. As it can be seen, these five geometric parameters were mutually independent (the largest correlation was between SBM and diameter, r = 0.055), as desired of potentially useful clinical metrics. Of these morphological features, consistent with findings in Figs 3-5, the maximum diameter of a lesion has the strongest direct correlation with ECAP (r = 0.66) and TFP (r = 0.59), followed by the aspect ratio (r = 0.50 for ECAP, r = 0.45 for TFP), which was determined largely by lesion length. Axial position (from the aorto-iliac bifurcation) correlated negatively with ECAP (r = −0.26) and TFP (r = −0.26), consistent with more proximal lesions tending to be less thrombogenic (Fig. 5), which may be due to the disturbed flow that occurs near branching vessels. For example, the entrance length was insufficient to reestablish laminar flow in the more proximal lesions and the coherent vortical structures seen in more distal lesions were disrupted, leading to greater mixing, higher ECAP, and lower TFP. Indeed, axial position correlated weakly with VS Area (r = 0.32) and inversely with VS Depth (r = −0.15), implying that vortical structures in proximal lesions tend to deflect from the centerline, widen towards the walls, and break-up early, all indicative of highly disturbed flows. VS Area correlates weakly with VS Bending (r = 0.16) and inversely so with VS Depth (r = −0.33), consistent with the above. VS Depth correlates with VS Bending (r = 0.42), suggesting that VSs that penetrate deeper into a lesion tend to move off-center before breaking apart (often a direct consequence of SBM). Of note, ECAP correlated strongly with TFP (r = 0.86), while PLAP correlated less well (r = 0.33), suggesting that ECAP, which is less difficult to compute, may be an acceptable proxy for TFP if computational resources are limited.
Panel B in Fig. 6 shows similar results except for the larger lesions (g 1 > 4.  function of the three most influential geometric parameters (i.e., diameter, axial location, and lesion length) as predicted by the kriging surrogate at 12,000 locations of the hyperparameter space. The plots highlight the aforementioned trends of TFP increasing with lesion diameter and decreasing for more proximally located aneurysms; the more complex relationship between TFP and lesion length is also seen in the rightmost plot. The scatter in  Fig. 6, generating HDMR contour plots did not require any explicit assumptions on the values of the remaining 3 input parameters that are not shown explicitly. Somewhat surprisingly, the plots resemble strongly the corresponding contour plots in Fig. 6, panel A. Finally, to evaluate whether the number of full simulations was sufficient, we ran 5-fold cross-validation tests on training sets of varying sizes built with randomly selected entries from our 179 full simulations. The top panel in Fig. 8 shows average relative errors when approximating TFP via kriging rather than from full hemodynamic simulations. As expected, larger training pools led to better accuracy, with initially fast improvements that tended toward a plateau. Expanding the training set from 20 to 60 simulations yielded a ~5% decrease in average error (from ~30% to ~25%), while an improvement of only ~1.5% (from 22% to 20.5%) was observed when 39 simulations were added to a training pool with 140 entries. We also assessed performance of the kriging model as a classifier of high-versus low-TFP lesions separated by a TFP threshold value of 2.5. Such categorization subdivided the 179 simulations into a low-TFP group with 44 members and a high-TFP one with 135 members. The bottom panel shows the percentage of aneurysms misclassified by the kriging approximation for both of these categories, and for training pools of varying sizes. As expected, the largest training set yielded the fewest misclassifications, amounting to ~10.5% of the considered lesions. Most of those (~9% of the total) were falsely classified as high-TFP aneurysms, a finding that persisted even when increasing the TFP threshold from 2.5 to 3.0 (whereby 14% of the lesions were erroneously labeled as high-TFP out of 18% misclassified aneurysms). Predictions from the kriging model could therefore be considered conservative, with only a few lesions (i.e., only 1.5% or 4% of all aneurysms for the 2.5 or 3.0 thresholds, respectively) erroneously labeled at a low-risk of thrombogenesis.

Discussion
There is now tremendous information regarding the histopathology, biomechanics, and clinical history of abdominal aortic aneurysms 3,8-10 , yet clinical decisions regarding intervention continue to rely primarily on the maximum dimension of the lesion or its rate of enlargement. Increasing evidence suggests, however, that ILT plays important roles in lesion enlargement and possible rupture 1-3 and there is motivation to include predictions of thrombus formation within an overall clinical decision rubric. Notwithstanding many prior studies on the biosolid or biofluid mechanics of ILT in AAAs, both idealized and patient-specific 5-7,11-14 , the pronounced patient-to-patient variability in lesion geometry has delayed general understanding. In contrast, studying idealized lesions is advantageous for one can prescribe and assess carefully individual and coupled morphometric features and hemodynamic characteristics that might otherwise be masked by confounding variations within patient-specific geometries.
In particular, we reasoned that intra-aneurysmal hemodynamics, and thus thrombus formation, depends largely on coupled effects of multiple morphological features of the diseased abdominal aorta. As a first study of such effects, we selected five parameters that collectively account for considerable patient variability (cf. Fig. 1). Maximum lesion diameter is an obvious metric to consider, as is the overall length of the lesion (thus addressing the key non-dimensional parameter of length/diameter). Vessel curvature, or tortuosity, can have significant effects on the hemodynamics, including generation of secondary flows, and was deemed important to include. Similarly, proximal branch sites and angles can influence flows distally, with effects expected to be more dramatic in proximally versus distally located AAAs. Finally, of many possible considerations of proximal branch geometry, we focused on renal artery positioning since asymmetric renal arteries can break the symmetry of the flow into the lesion. Even with only five parameters, the curse of dimensionality leads to considerable computational challenges. Like patient-specific simulations 5 , calculating hemodynamics (involving unsteady 3-D Eulerian and particle-tracking Lagrangian analyses) within idealized lesions placed within a representative human abdominal aorta is computationally expensive, with each simulation requiring approximately 1000 CPU-hours on a high-performance cluster (128 CPU cores) and generating nearly 300 GB of raw data each. A straightforward parametric study was therefore impractical and so too a Monte Carlo approach was deemed not viable.
We thus sought to reduce the computational expense by minimizing the number of simulations that would be needed while enabling reliable interpolations across the 5-D parameter space. Adaptive sparse grid collocation was selected based on prior successes in studying geometric effects on hemodynamics 15 . As outlined in the Supplemental Information, multiple levels of selection led to full simulations at 179 collocation points, which in turn enabled reliable interpolation via kriging with a predicted low classification error (less than 4%) for false positives of low TFP. Not surprisingly, maximum diameter correlated well with high TFP, consistent with the clinical observation that larger lesions tend to harbor an ILT that often forms within the region of greatest aneurysmal bulging, namely ILT exists in up to 95% of lesions 3.5 cm or more in diameter 16 . Nevertheless, careful evaluations of select groups of simulations revealed couplings amongst multiple parameters (Figs 4 and 5), some not necessarily intuitive. For example, lesions with equal maximum diameters had very different thrombus formation potential depending on axial position, which was likely influenced dramatically by the integrity or break-up of the vortical structures 7 .
Importantly, interpolation-based maps from kriging suggested that, given the same diameter, increased thrombus formation potential may arise in lesions that are located more distally and/or are shorter; spinal curvature and renal artery offset affected the results much less. Because multiple morphometric features contribute to the hemodynamics proximal to (which influences platelet activation) and within (which influences endothelial vulnerability) the aneurysmal region, single metrics such as maximum diameter are unlikely to be predictive. Indeed, it has similarly been suggested that, despite its continued clinical usage, a simple metric such as maximum diameter is not as good of a predictor of aneurysmal rupture potential as local wall stress 17 , which can incorporate information on wall geometry, properties, and applied loads. Given that AAAs rupture when wall stress exceeds wall strength, multi-factorial metrics have been introduced to predict rupture potential based on wall strength 18,19 . Using multiple linear regression, such models can include geometric metrics such as normalized maximum diameter but also effects such as medical history and sex. Until we understand precise mechanisms, and perhaps even thereafter, there is similarly a need to investigate potential effects on ILT due to age, sex, past medical history, blood composition, prescription drugs, and potential co-morbidities such as hypertension or diabetes. Of course, consideration of such metrics will demand large patient-specific data sets, which was beyond the present scope. This study was clearly but a first assessment, with consideration of idealized lesions simplifying the parametric evaluation. We computed and compared peak ECAP and TFP by taking the 99 th percentile of luminal surface fields. Both metrics tended to be highest within the aneurysmal domain, rendering both ECAP-99p and TFP-99p potentially useful indices. Other metrics or reductions of field quantities could be envisioned. Moreover, other geometric parameters of potential importance could be considered. Among these, cross-sectional eccentricity (cf. Fig. 1), the degree and direction of bulging of the aneurysm from the natural curvature of the infrarenal aorta, the angle of the iliac bifurcation, and the presence of an extra renal artery (which manifests in approximately one third of all cases) 20 could be considered. There is similarly a need to consider possible effects of peripheral vessel disease, such as stenosis or occlusion, and effects of spinal cord injury or lower limb amputation, all of which can adversely affect AAA hemodynamics [21][22][23] . Such simulations were simply beyond the present scope, particularly noting that one would need to appropriately adjust outlet conditions to account for remodeling of the contralateral iliac/femoral in cases of unilateral disease or amputation, hence increasing the parametric space considerably.
Increasingly sophisticated multiscale and/or systems-based models are becoming available for studying thrombus formation and deposition in complex hemodynamic fields [24][25][26][27] . Such models offer the promise of understanding mechanisms, which in turn could point to targeted therapeutics, and thus should be pursued in the future. Nevertheless, such studies require identification of many model parameters, and particle-based models remain computational impractical for large computational domains such as the human infrarenal aorta and distal common iliac arteries. Hence, we employed a simpler phenomenological model of thrombogenicity that has been validated against results from 9 patients 5,6 and allows reasonable computational solutions (still ~1000 CPU hours each) to explore a 5-D geometric parameter space. We submit that, given that thin ILTs appear to incite a dramatic proteolytic insult to the aneurysmal wall [28][29][30] , having a computational model that depends primarily on easily obtained imaging based information and predicts well when (and where) an ILT will form 5 could be sufficient to time clinical intervention during potential periods of rapid aneurysmal remodeling and expansion 4 . Mechanistic models, in turn, could help inform that clinical intervention and may not need to be patient-specific. Clearly, much remains to be accomplished before such models become part of the standard of care.
In conclusion, this paper presents the first melding of methods of computational hemodynamics, stochastic collocation, and Kriging to examine potentially coupled effects of common morphological features on ILT formation in AAAs. We have also proposed novel methods to extract and interpret quantitative hemodynamic features that drive ILT formation, revealing for example the importance of vortical structure dynamics in differently sized and shaped lesions. Simple maps (Fig. 6) suggest that ECAP correlates well with TFP, which is important from a practical perspective. ECAP can be determined via a standard (Eulerian) solution of the Navier-Stokes equations, not requiring the less straightforward (Lagrangian) particle tracking that is essential for calculating TFP. Hence, consistent with findings of thrombus formation in iliac artery aneurysms 31 , ECAP may be a computationally less expensive alternative for predicting possible thrombus formation. Regardless, computational studies are clearly critical for predicting thrombogenicity given the identified nonlinear couplings amongst parameters and they suggest that simple multilinear regressions will not be sufficient to determine a single scalar metric having clinical utility.

Methods
Morphological Quantification. Information from previous studies reported in the literature 32,33 as well as from direct measurements on nine patient-specific reconstructions that we analyzed previously 5,6 suggested reasonable ranges for each of the five parameters ( Table 1) that are revealed by Fig. 1. Idealized geometric model. Unsteady, 3-D simulations of patient-specific hemodynamics within the abdominal aorta are computationally expensive, particularly when computing TFP, which requires standard Eulerian simulations via Navier-Stokes solutions plus Lagrangian particle tracking in post-processing to infer platelet strain histories. Spatial discretization of the former often requires 5+ million elements and temporal resolution often requires 3600+ time steps per cardiac cycle, with four or more cycles common to ensure convergence of the unsteady solution 6 . Total CPU time can often exceed 1000 hours per patient-specific model, which can be reduced in real time using parallel computing but is nonetheless expensive. Notwithstanding the difficulty of obtaining and reconstructing large numbers of patient-specific medical images to generate models of the hemodynamics, the number of geometries that can be investigated practically within a standard parametric study is limited by the computational expense. To facilitate a larger parametric study, we combined an analytical descriptor to define classes of AAA geometries and a statistical approach (sparse grid collocation to identify optimal parameter combinations and kriging to interpolate results) to reduce the total number of full simulations needed to characterize the five-dimensional parametric space (see below).
Inspired by previous work linking AAA morphology with the wall biomechanics 34 , each idealized aneurysm was modeled while assuming circular cross-sections orthogonal to a specified centerline, with luminal radius r given as a function of axial location z along that centerline: where r IAA is the radius of the healthy segment of the infrarenal abdominal aorta (IAA), r AAA is the maximum radius of the AAA (that is, g 1 /2), and z 0 is the axial position of the maximal radius (cf., g 2 ), all in mm. Together, the model parameters c 1 (unit of length), c 2 (unit-less when axial distance z is normalized per unit length, l o = 1 mm), and c 3 (unit-less) dictate the length (g 3 ) and curvature (g 4 ) of the aneurysmal dilatation. To approximate aneurysmal curvatures for multiple lesion lengths and diameters, we let c 1 = 0.2 mm and c 3 = 6.25⋅10 −4 , while c 2 was varied according to a best-fit regression to match the desired length of the lesion (see below). The centerline of the lesion was parameterized by the axial coordinate z, while the x-axis points to the left of the patient and the y-axis points anterior. The centerline is given by: where parameters c 4 (unit of length) and c 5 (unit of inverse length) specify, respectively, the desired forward bulging, or deflection, and angle of deflection, θ. The parameter z 1 defines the bending offset where maximal bulging angle θ is achieved. We let z 1 = 60 mm from the iliac bifurcation. The parameters θ, c 4 , and c 5 are related via: This relation provides two degrees of freedom in ascribing curvature with aortic dilatation. To limit this to one variable, c 5 was set to π 1 6 3 mm −1 , which yields a forward deflection (i.e., c 4 , g 4 ) of 30 mm when θ = 30°. This selection allows the entire span of the anterior bulging to remain within the infrarenal aorta (i.e., prior to the iliac bifurcation) without any need to vary the curvature of the suprarenal aorta. Distal to the iliac bifurcation, centerlines for the iliac arteries follow a reflection of the above equation. The final centerline approximates the geometry observed in vivo due to the curvature of both the spine and the AAA itself.
It is important to distinguish the implicit parameters c i , which are used in the radial, centerline, and bending equations above, and the explicit parameters g i , which define the desired morphological features of each aneurysm. We used Mathematica 10.0 (Wolfram Research, Champaign IL) to perform a mapping between these two sets of parameters. With c 1 , c 3 , and c 5 fixed, c 4 = g 4 , and the desired spinal bending parameter θ determined via Equation (3), c 2 was computed via a numerical interpolation for each specified radius and length. Specifically, roots of r(z)−1.1r IAA were computed for different values of c 2 and r AAA and the Euclidean distance between these roots defined the length of the aneurysm (g 3 ), recalling that the ends of the aneurysm are defined by a 10% increase in radius from r IAA . The resulting set of 37 tuples spanned maximum diameters g 1 ∈ [30,50]  In addition to the idealized lesion, the overall computational domain included the iliac arteries, renal arteries, superior mesenteric artery (SMA), and celiac artery (Fig. 1) since each of these branches and associated flows can affect flow within an AAA. In particular, the axial offset between the ostia of the two renal arteries (g 5 ) was included for it can break the symmetry of the flow into the lesion and thus affect vortical structures. To minimize entrance effects, the computational domain was extended from the suprarenal abdominal aorta (SAA) into the distal descending thoracic aorta using segments of appropriate length and diameter (Fig. 1). The inferior mesenteric artery, which typically stems from the anterior portion of the IAA, was omitted because it is rarely identifiable in AAA scans. The lengths, diameters, branching origins, and branching angles of all included vessels are listed in Table 2. Lengths are computed along the centerline, with axial position representing either the ostia of a branching vessel or the axial domain of an aortic segment. All axial positions are given via the z coordinate, which is zero at the iliac bifurcation. The branch angle represents the absolute value of angular deflection from the abdominal aorta. Clock angles represent the direction of branching of a given vessel. Statistical learning of the parameter space. Our primary goal was to investigate parametrically the influence of vascular and aneurysmal morphologies on infrarenal hemodynamics and to predict potential ILT formation. Because of expected complex nonlinear coupling of such morphological features (e.g., lesion position and curvature should modulate the effects of maximal diameter), we sought to examine myriad combinations of five key parameters on the intra-aneurysmal TFP. Yet, appropriate spanning of this five-dimensional space over pathophysiological ranges (Fig. 2A 5 5 = 3125 full simulations. One approach to reduce the burden of spanning high-dimensional parameter spaces is the Chebyshev sparse grid stochastic collocation method [39][40][41] . Let Ξ represent the stochastic geometric parameter space, with dimensions defined by g i (i = 1, 2, …, 5). Rules provided by the Chebyshev sparse grid method yielded points ξ l,j (l = 0 ... 3, j = 1... n l ) that would efficiently probe Ξ. As is common in Smolyak-type constructions 15 , our sparse grid was built hierarchically, meaning that the set of collocation nodes was progressively expanded to reach the desired accuracy. We constructed 4 levels (l = 0 ... 3) of progressive refinement (see Supplemental Information), which resulted in collocation points at the roots of 5-D Chebyshev polynomials of degree 1 to 4.
After running simulations at collocation points for each refinement level, multiple hemodynamic metrics were extracted in post-processing, including wall shear stress indices, thrombus formation potential, and characteristics of vortical structures (see Sections 2.5-2.7). A 5-D Gaussian process regression (kriging) model with a radial basis function kernel was then trained to interpolate each of the hemodynamic features over Ξ. Our implementation used the open-source package GPy 42 . From this, the value and uncertainty of each metric was estimated at any point in Ξ, and thus for any complex geometry within the parameter space.
To minimize further the number of simulations needed, we employed the adaptive refinement strategy introduced by Sankaran and Marsden 15 . This allowed us to discard 62 out of 180 collocation points initially prescribed for refinement level 3, since results from levels 0, 1, and 2 showed that certain areas of the parameter spaces were already characterized sufficiently. Thus, we considered a total of 179 lesion morphologies. Recall that Fig. 2B shows a 3-D projection of Ξ to visualize collocation points while compressing the g 4 and g 5 dimensions. The nodal distribution reveals that Chebyshev grids tend to be denser along the mid-range and around boundary vertices of the parameter space. Also, asymmetries in the collocation points resulted from the adaptive refinement scheme, which revealed, for example, that limited new information is obtained when varying aneurysm length when at the smallest lesion diameter. Figure 2C details two slices of the 3-D projection along g 1 − g 2 and g 1 − g 3 mid-planes.
Computational hemodynamics simulations. Full 3-D finite element solutions of the Navier-Stokes equations were computed for each of the 179 model geometries. Inlet and outlet boundary conditions were prescribed for rest conditions, which yields worst case scenarios for ILT since exercise may limit thrombus formation 5 . The volumetric flow rate over a cardiac cycle at the inlet of the descending thoracic aorta was obtained from Les et al. 43 and the associated inlet velocity profile was computed as a Womersley pulsatile laminar flow 44 , with 15 Fourier series terms. The cardiac cycle was assumed to have a period of 0.9 seconds, corresponding to a heart rate of ~67 beats per minute. Outlet boundary conditions were specified using a lumped-parameter circuit model (Windkessel resistance-capacitance-resistance boundary conditions, see Fig. S1 in supplemental material), as is now common 45,46 . Values of normal outlet conditions were prescribed from Xiao et al. 47 . The blood was modeled as a Newtonian fluid, with a constant viscosity and density of 0.004 g/mm•s and 0.0016 g/mm 3 , respectively. All vessel walls were assumed rigid similar to our previous studies wherein potential effects of this assumption were discussed and accepted 5 .
The incompressible Navier-Stokes equations were solved in SimVascular, using a stabilized finite element flow solver, to generate unsteady pressure and velocity fields within the 179 model aneurysms as well as in the normal aortic segments just proximal and distal to each lesion. The cardiac cycle was subdivided into steps of 0.001 seconds, with computational results saved every 5 steps (0.005 seconds). Simulations were performed for four cardiac cycles having a uniform period T = 0.9 seconds. To allow transient flow patterns to dissipate and to establish periodicity, we used results from the last cardiac cycle for post-processing and Lagrangian particle tracking simulations as described previously 5 .
Post-processing allowed computation of several biomechanical metrics from the velocity fields. Of primary relevance were two luminal surface metrics: a Time-Averaged Wall Shear Stress magnitude TAWSS ( ) that is normalized by the mean value of TAWSS in the descending thoracic aorta, and a Oscillatory Shear Index (OSI) that measures changes in the WSS field due to complex oscillatory flows. The quotient of these metrics defines the Endothelial Cell Activation Potential (ECAP), which delineates regions of high OSI and low TAWSS where the vessel is thought to be particularly susceptible to thrombus formation 5,48 . TAWSS, OSI, and ECAP are given by the following point-wise metrics defined on the luminal surface: Vortical structure analysis. Shear stress indices such as TAWSS, OSI, and ECAP are highly affected by the dynamics of vortical structures near the wall. To investigate how morphological features of a lesion might influence the path taken by vortical structures throughout the cardiac cycle, we examined λ 2 iso-surfaces 7 . λ 2 is the second largest eigenvalue of the tensor D 2 + Ω 2 , where D = 0.5(l + l T ) is the symmetric part of the spatial velocity gradient tensor l and Ω = 0.5(l − l T ) is its skew-symmetric counterpart. A surface of a negative λ 2 encloses a pressure maximum, which typically indicates the presence of a vortex 49 . To analyze paths and other key features of vortical structures during a cardiac cycle, we built iso-surfaces directly on minimum values of λ 2 at each location. As such, our plots merged the dynamic evolution of the vortical structure into a single snapshot showing envelopes of the traces of the vortices. While this representation tends to hide temporal features of vortex dynamics (e.g., it does not discriminate whether a vortical structure moves on quickly or stations longer within a given region), it simplifies the analysis of vortical paths and their locations. Specifically, we measured for each case: (1) VS Depth , namely how deep the main vortical structures penetrated into the aneurysm before breaking into smaller structures; (2) VS Area (in mm 2 ), the cross-sectional area enclosing vortical structures when they break-up; and (3) SCiENTifiC REPORTS | (2018) 8:13273 | DOI:10.1038/s41598-018-31637-6 VS Bending (in mm), the distance between the centerline of the aneurysm and the center of the cross-sectional area where the vortex breaks-up.
Thrombus formation potential. Standard hemodynamic simulations allow one to compute ECAP, a measure of wall susceptibility to thrombus deposition. A second significant contributor to thrombogenesis is the presentation, then aggregation, of activated platelets at the susceptible site. This activation can depend highly on the mechanical shear history 7,50-53 . To determine which regions of the luminal surface were presented with shear-activated platelets, we used Lagrangian particle tracking simulations, as described previously 5 . Briefly, a secondary mesh with nominal resolution of 0.5 mm and increased curvature refinement was generated to fill the computational domain and a particle was placed at each node of the mesh. We then advected the platelets backward through time, looping over 8 cardiac cycles using flow results to track the trajectories and shear history of model platelets that ended up adjacent to the vascular wall. Ten simulations were run in parallel, each with a different evenly-distributed injection time in the cardiac cycle, and results were averaged. The particle tracking simulations were used to calculate the PLatelet Activation Potential (PLAP), proposed by Shadden and Hendabadi 53 and employed by Di Achille et al. 5 . PLAP is a non-dimensional volumetric index of the cumulative shear rate history experienced by particles that arrive at a particular point in the domain (i.e., the starting point prior to advection backward through time). It is defined as where |D(x(τ),τ)| is the Frobenius norm of the symmetric part of the spatial gradient of the velocity tensor, t is the time of injection of the particle, and T is the period of the cardiac cycle. A surface PLAP is generated by averaging the volumetric PLAP as described in Di Achille et al. 5 and normalizing by its average value in the distal descending thoracic aorta. PLAP provides significant insight into stagnant flows, particle attractors, and the typical delivery distribution of activated platelets to the endothelium, all of which are relevant to thrombogenesis. The ECAP and PLAP were combined to generate a Thrombus Formation Potential (TFP):

TFP ECAP PLAP
OSI PLAP TAWSS (6) Co-localization of wall priming (high ECAP) and advection of shear-activated platelets (PLAP) thus creates a surface index that predicts thrombus initiation at a given location and time 5,6 . This metric can be used as a proxy for thrombogenicity in a given lesion and thus has clinical applicability.
Global sensitivity analysis. Sobol's variance-based sensitivity analysis provides an effective strategy to decompose the variance of output metrics in terms of relative direct and indirect contributions of the input parameters 54 . This analysis nevertheless requires extensive probing of the hyperparameter space to rule out potentially spurious correlations and to quantify objective sensitivity measures. To render the global analysis feasible in our case, we applied Sobol's analysis on predictions from optimized kriging interpolations (see Section 4.4). Samples were collected following Saltelli's scheme 55 , which in its most robust form requires r(2k + 2) kriging evaluations, where k is the number of input parameters and r indicates the resolution chosen to characterize each input dimension. Following common practice, we chose r = 1000 for a total of 12000 kriging evaluations. Sample selection and sensitivity indices were calculated using the open-source python library SAlib 56 . In particular, we computed and compared first-order and total effect sensitivity indices. We also used analysis of variance theory to generate alternative visualizations of the combined effects of morphological parameters on TFP. Specifically, we built high-dimensional model representations of TFP as functions of selected couples of inputs according to [57][58][59][60] = | = … ( ) ( ) TFP g g E TFP g g i j , , , i j i j where g i , g j represent values of a selected input features, E(.) is the expected value operator and | indicates a conditional probability. The conditional expected values of TFP were computed via inference on 2-D Gaussian process regressions optimized on the 12000 extracted samples.

Data Availability
All input data needed to reproduce the results are provided in the text or tables.