Biomechanical rupture risk assessment of abdominal aortic aneurysms using clinical data: A patient-specific, probabilistic framework and comparative case-control study

We present a data-informed, highly personalized, probabilistic approach for the quantification of abdominal aortic aneurysm (AAA) rupture risk. Our novel framework builds upon a comprehensive database of tensile test results that were carried out on 305 AAA tissue samples from 139 patients, as well as corresponding non-invasively and clinically accessible patient-specific data. Based on this, a multivariate regression model is created to obtain a probabilistic description of personalized vessel wall properties associated with a prospective AAA patient. We formulate a probabilistic rupture risk index that consistently incorporates the available statistical information and generalizes existing approaches. For the efficient evaluation of this index, a flexible Kriging-based surrogate model with an active training process is proposed. In a case-control study, the methodology is applied on a total of 36 retrospective, diameter matched asymptomatic (group 1, n = 18) and known symptomatic/ruptured (group 2, n = 18) cohort of AAA patients. Finally, we show its efficacy to discriminate between the two groups and demonstrate competitive performance in comparison to existing deterministic and probabilistic biomechanical indices.


Introduction
An abdominal aortic aneurysm (AAA) is a slowly progressing vascular disease, causing an enlargement of the infrarenal aorta and is considered pathological if the aortic diameter exceeds 30 mm [1]. AAA prevalence has been reported within a range of 1.2% to 3.3% in men older than 60 years based on several studies in western societies [2]. In most cases, AAAs develop asymptomatically over several years, but they can rapidly turn into a serious clinical emergency in case of rupture. More than 50% of patients with a ruptured AAA die before reaching the hospital [1] and perioperative mortality rates range from 40% to 60% [3]. To prevent such a disastrous scenario, the clinical guidelines from the US-based Society for Vascular Surgery recommend elective repair for AAA patients with an aortic diameter greater or equal to 55 mm, regular screening intervals for patients with smaller-sized AAAs and onetime screenings for AAAs in men and women above a certain age and based on established risk factors [1]. This maximum diameter recommendation is based on a risk assessment, where the risk of rupture is weighed against the mortality risk of an elective repair. While the latter risks are relatively well-known, aneurysm rupture is a complex biomechanical failure event. With the increasing use of endovascular repair (EVAR) over open surgical repair (OSR) [4], however, which can be attributed to the significant short term mortality benefit of EVAR (1.4% compared to 4.2%) [1], interventional risks have become a less important factor in the risk assessment process.
Nonetheless, a biomechanical rupture risk assessment can provide an additional important piece of information. It enables the possibility to provide patient-specific screening guidelines, avoid unnecessary interventions [5] and support the clinical decision process for cases that are not covered by the clinical guidelines. The Society for Vascular Surgery's 55 mm recommendation, e.g., only holds for patients "at low or acceptable surgical risk with a fusiform AAA" [1]. Furthermore, there are no clear or only weak recommendations for women with AAAs of size 50-54 mm, aneurysms with non-fusiform geometries, smaller AAAs [6], or patients at higher surgical risk. In addition to that, not all AAAs are suitable for EVAR, with higher complication rates for AAA cases that are not covered by the instructions for use [7]. Lastly, recent meta-studies (e.g. [8,9]) on the long term outcomes of EVAR versus OSR could not detect any differences with regards to the all-cause mortality or even concluded in favor of OSR.
In this paper, we present a highly personalized, probabilistic framework for the biomechanical quantification of AAA rupture risk. The framework builds upon a comprehensive database, consisting of tensile experiments that were carried out on 305 AAA tissue samples from 139 patients and corresponding non-invasively and clinically accessible patient data. The approach consistently incorporates the available statistical information in terms of probability distributions in order to account for patient-specific uncertainties about relevant vessel wall properties. We emphasize the importance of accounting for these uncertainties and demonstrate that this leads to a more accurate individualized rupture risk assessment as compared to deterministic approaches.
Our work builds upon previous efforts by our group and collaborators regarding the biomechanical modeling and characterization of AAA in-vivo behavior [10][11][12][13][14][15][16], as well as several previous studies indicating that biomechanical indices are more accurate predictors for AAA rupture risk than the clinically established maximum diameter criterion [17][18][19][20][21][22][23][24]. In contrast to the approaches in [17-20, 22, 24], however, we advocate a probabilistic treatment to account for uncertain vessel wall properties. Our work thus goes along the lines of [21], but with the key difference that it includes the stiffness parameters of the AAA vessel wall as statistical quantities, uses patient-specific vessel wall properties and accounts for statistical correlations among these properties.
The paper is organized as follows. Motivated by a failure-based criterion, our rupture risk index is formulated in Section 2.1 incorporating patient-specific statistical information. Section 2.2 defines the biomechanical AAA model and specifies the probabilistic regression model to obtain personalized vessel wall properties. In Section 2.3, a method for the efficient evaluation of the rupture risk index is proposed and in Section 3, the framework is applied on a total of 36 retrospective, diameter matched asymptomatic (group 1, n = 18) and known symptomatic/ruptured (group 2, n = 18) cohort of AAA patients.

Materials and methods
2.1 Failure-based probabilistic quantification of rupture risk 2.1.1 Rupture as an event of material failure. From a mechanical point of view, rupture represents an event of local material failure at a point x in the aneurysm wall, which motivates its definition via a failure function φ(x) and the failure criterion φðxÞ > 0; at any x: We limit ourselves to stress-based failure and define rupture as an event where the local wall stress measure σ(x) exceeds the local wall strength σ γ (x). This results in the failure function φ(x) = σ(x) − σ γ (x), or the criterion sðxÞ > s g ðxÞ; at any x: ð2Þ Using the equivalent von Mises stress σ vm (x) as the local stress measure σ(x) and an assumed spatially constant wall strength σ γ , this criterion can be evaluated as where s max vm is the maximum von Mises stress s max vm ¼ max x s vm ðxÞ. It is important to note that the above definition does not incorporate any aspect about failure over time. In order to be able to include time in the analysis, i.e. to make a statement about the risk of rupture in the next year, one would require knowledge about the future progression of the AAA for this patient, such as a model for the aneurysm growth and change in vessel wall properties. Since there is hardly any knowledge about these aspects, we limit the further discussion to a rupture risk assessment at the point of time of the acquired data. While there are sudden events like calcification-induced formation of saccular aneurysms, we assume that in most cases an AAA is a slowly progressing disease and thus our approach has, at least for the near future, sufficient predictive capability.
2.1.2 Existing criteria and rupture risk indices. Rupture risk estimation for AAAs has been an ongoing research topic over several decades, with many attempts to establish decision criteria for clinical practice. The maximum diameter criterion [1] still represents the most widely used criterion for decision making today. It is often justified by Laplace's law, which states that the vessel wall stress is proportional to its diameter in spherical geometries. Based on this and with data obtained from several clinical studies, a very simple criterion, has been formulated, relating the patient's AAA diameter d to a critical maximum diameter d max . While established in clinical practice and easy to apply using CT or ultrasound imaging, this criterion has often been criticized [25] and is an ongoing subject for discussion [6]. With growing computational resources and advances in the modeling of biomechanical material behavior, the simulation of patient-specific AAA models has been advanced by several research groups. Experiments on harvested AAA samples were able to reveal material parameters and failure properties. In addition with regression models [13,26,27] for the prediction of the individual wall strength, this enabled the definition of biomechanics-based indices [19,20,22,28], such as the rupture potential index (RPI) relating the von Mises stress to the wall strength. Furthermore, it could be shown [19,20,24] that these indices can be better rupture risk indicators than the maximum diameter criterion. Experimental testing [13,26,27] also revealed significant inter-and intra-patient variabilities in the mechanical properties of AAA tissue, motivating a probabilistic approach to rupture risk estimation [14,16,21] and resulting in the probabilistic rupture risk index (PRRI) [21] PRRI ¼ where the authors used distributions for the wall thickness and wall strength that were fitted on cohort data published by our group [13].

A novel probabilistic approach.
In this work, we propose a novel failure-based, probabilistic rupture risk indicator that consistently incorporates all available statistical information and accounts for correlations among vessel wall properties. Fig 1 (left) illustrates the rationale for our approach, showing how part of the available data is directly involved in the estimation of the risk of rupture, while another part affects the evaluation of the computational model. In general, this data will be correlated, resulting in correlated quantities for the evaluation of rupture risk and necessitating a reformulation.
To that end and recalling the rupture criterion from Eq (3), we can calculate the probability of rupture over the joint probability distribution pðs max vm ; s g Þ as where 1 s max vm >s g is the indicator function defined as This formulation can be easily extended to, e.g., spatially varying vessel properties using Eqs (1) or (2) as failure events. Furthermore, it includes the PRRI in Eq (6) as a special case, when choosing pðs max vm ; s g Þ ¼ pðs max vm Þpðs g Þ. Lastly, it allows for a straightforward visual interpretation as illustrated in Fig 1 (right). The plot shows the joint probability distribution pðs max vm ; s g Þ and visualizes the rupture event area in red. The blue area implies a high probability for the joint occurrence of the corresponding stress and strength values. The probability of rupture P rupt is simply the volume of this density within the triangular rupture event area. Thus, the larger the overlap between pðs max vm ; s g Þ and the red area, the higher P rupt .
2.2 Data-informed patient-specific AAA models 2.1.1 Geometry creation from CT imaging and meshing. Patient-specific 3D AAA geometries are reconstructed via a semi-automatic segmentation process from CT imaging data using the software ScanIP (Synopsys, Mountain View, California) and based on a protocol as described in [12]. The minimal requirement for the spatial resolution of CT scans was 1 mm and for the slice thickness 3 mm. The upper boundary for the segmentation was the branching of the renal arteries and the lower boundary below the bifurcation at the iliac arteries. Due to the small thickness of the AAA wall, its low contrast and the limited resolution of the CT images, it is only possible to extract the blood lumen and intraluminal thrombus (ILT) geometries. After segmentation, the ILT geometry is exported as a surface model for meshing.
In a next step, we use the software Trelis (csimsoft, American Fork, Utah) and bi-linear quadrilateral elements to mesh the abluminal ILT surface. From this surface mesh, the arterial wall layer is extruded with a specified, spatially constant thickness t, resulting in a tri-linear, single layer, hexahedral mesh for the AAA wall. Finally, linear tetrahedral elements are employed for the meshing of the complex ILT geometry and a layer of linear pyramid elements as a transition for mesh compatibility between AAA wall and thrombus. Element sizes were set to 1.6 mm, corresponding to the median of measured thicknesses of AAA wall specimens in our database and leading to hexahedral elements of shape 1.6 mm × 1.6 mm × t for the AAA wall. A mesh convergence study has been performed to assess that the chosen spatial mesh resolution is sufficient in the context of our application. The meshing procedure is also described in [29] in more detail.

Biomechanical modeling.
Previous studies have shown that in order to accurately describe the biomechanical behavior of AAAs, a sufficient model complexity is required [20,30], while results by [31,32] indicate that also simpler models might be appropriate. For our purposes, we employ the finite deformation boundary value problem of nonlinear elasticity: where O 0 is the reference configuration of the AAA, u denotes the displacement field, F = I+ ru the deformation gradient, S the second Piola-Kirchhoff stress tensor and σ the Cauchy stress tensor. On the Neumann boundary γ σ , i.e. the luminal ILT surface, an orthonormal loadt ¼ À pn is applied, with the pressure value p and the unit outward surface normal n in the current configuration. Furthermore, at the proximal and distal end surfaces of the AAA model, Γ u , we employ a Robin-type boundary condition with spring supports following [29,33]. The stiffness parameter k s is per unit reference area and set to 100 kPa/mm in this study, while N is the unit outward surface normal in the reference configuration.
To model the constitutive behavior of the ILT, we use the strain energy function proposed in [34] C ILT ð � I 1 ; � I 2 ; JÞ ¼ cð � I 2 1 À 2 � I 2 À 3Þ þ C vol ðJÞ ð12Þ and a linearly decreasing stiffness c from the luminal to the abluminal ILT surface [12]. � I 1 and � I 2 are the first and second invariants of the modified right Cauchy-Green deformation tensor [29]. The strain energy function employed for the AAA wall material is [35,36] with stiffness parameters α and β. Both strain energy functions are equipped with an additive volumetric component including the bulk modulus with parameters � k ILT ¼ 8c and � k wall ¼ 2a for the employed ILT and wall material models and a Poisson's ratio of ν = 0.48 [12].
To obtain a pressurized in vivo configuration of the AAA, the MULF prestressing method [10,11] is used, where the applied load corresponds to the mean arterial pressure (MAP = 1/3 systolic pressure + 2/3 diastolic pressure). From this prestressed configuration, the pressure is raised by 50% to simulate elevated blood pressure conditions [21]. Following [19] and for comparability reasons, the values for the systolic and diastolic pressures were set to 121 mmHg and 87 mmHg for all cases, respectively, resulting in a MAP of 98.33 mmHg.
With the finite element discretization from Section 2.2.1, a nonlinear system of equations is obtained, which is solved using an in-house finite element code. We note that in this study we neglect the effect of calcifications in the AAA for simplicity and assume constant vessel wall thickness t and stiffness parameters α and β throughout the aneurysm. Furthermore, we evaluate the maximum von Mises stress as the 99th percentile of the von Mises stress field in the aneurysm.
For the remainder of this work, we will use the parameter to quantity of interest (QoI) map with parameter vector θ ¼ ½t; a; b� T 2 R 3 þ and QoI s max vm 2 R þ to denote the forward problem. Thus, calculating s max vm ðθÞ for one realization of t, α and β will involve one evaluation of the nonlinear finite element model.

Patient database.
The modeling of patient-specific vessel wall properties here is based on data that has been collected during several research projects between 2008 and 2017 on the mechanobiological behavior of AAAs [13,15]. The study was approved by the ethics committee of the University Hospital rechts der Isar, Technical University of Munich. AAA patients undergoing elective OSR (including emergency repair due to rupture) at the University Hospital rechts der Isar in Munich, Germany, were added to the database, whenever it was possible to extract tissue samples for mechanical testing. Apart from anamnesis and CT imaging data, hemograms were evaluated and one or more AAA tissue samples harvested during OSR. These samples were mechanically and histologically investigated, resulting in an exhaustive retrospective AAA database. Further information on data collection and experimental testing can be found in [13,15]. To date, the database contains a total number of 305 entries from an equal number of tissue samples that were collected from 139 patients.
The data can be split into two groups. Invasive properties (cf. Table 1), denoted as Θ ¼ ½t; a; b; s g � T 2 R 4 þ , are properties, which have been determined retrospectively from AAA tissue samples and cannot be obtained for a prospective patient by using clinically established methods. They are, however, essential for the biomechanical modeling and simulation of AAAs and the calculation of the probability of rupture using Eq (7). Non-invasive properties (cf. Table 2), denoted by ξ, on the other hand, can be determined with standard methods in the clinic. The subrenal diameter in Table 2 is measured directly below the renal arteries. If the aneurysm reached the renal arteries, the aortic diameter between the celiac artery and the superior mesenteric artery minus 2.5 mm was used instead [12].
Based on correlations between the invasive and non-invasive properties [13], the goal is to construct a statistical model for the patient-individualized prediction of vessel wall properties Θ(ξ) for a prospective new patient with non-invasive properties ξ. While this process is described in Section 2.2.4, a preprocessing step for the dataset is essential, since values are missing both in the invasive and non-invasive properties for several cases in our database. Moreover, the relatively small number of available data, but high number of non-invasive properties, requires a feature selection process to identify the most important properties in ξ. Similar to [15], we conduct the following preprocessing steps.
Non-invasive features in ξ, where more than 30% of the data points had missing values and patients with more than 30% of missing features were excluded and all other missing non-invasive properties imputed with the corresponding median value across the population. As a consequence, the four parameters calcium, high-sensitivity C-reactive protein (hsCRP), creatine kinase and fibrinogen were disregarded. Afterwards, all non-invasive features were normalized.
Based on a correlation analysis using Spearman's rank correlation coefficient (cf. S1 Table), the total number of features was reduced to a final selection of 8 variables: maximum AAA diameter, maximum thrombus thickness, AAA length, subrenal diameter, thrombocytes, hemoglobin, mean corpuscular hemoglobin (MCH), mean corpuscular volume (MCV). The restriction was done using a sequential forward selection algorithm similar to [15]. In an attempt to keep the number of non-invasive parameters small, we iteratively added the highest correlating non-invasive parameters to the GP model (see Section 2.2.4) until no further improvement in the leave-one-out cross-validation (LOOCV) scores could be observed. We note, however, that this does not imply that other non-invasive features such as sex, medication or anamnesis parameters do not have an influence on the biomechanical properties of the AAA wall. The resulting dataset D ¼ fξ i ; Θ i g n data i¼1 , that was used for the analysis in Section 3, consisted of n data = 251 data points from 113 individual patients and is available as supplementary information to this study (cf. S2 and S3 Tables).

Prediction of invasive vessel wall properties.
Previous approaches to create models for the AAA wall thickness, stiffness parameters or strength were either deterministic [12,26],

PLOS ONE
based on cohort statistics [21], or did not account for correlations among the vessel wall quantities [15]. In the following, we make use of a multivariate Gaussian process regression model [37][38][39] to address these shortcomings and achieve the following desiderata: 1. Patient-specific modeling: obtain personalized estimates for the vessel wall quantities Θ based on correlations with the non-invasive properties ξ of a specific, prospective patient.
2. Probabilistic treatment: take into account the uncertainties in the predictions for Θ (do not ignore statistical information). These can be determined with standard methods in the clinic and will be used as feature variables to predict the invasive properties of a prospective AAA patient. https://doi.org/10.1371/journal.pone.0242097.t002 3. Dependencies: model the correlations among the invasive properties Θ in order to obtain a more accurate probabilistic description and avoid physically implausible parameter configurations.
As a result, given the non-invasive properties ξ of a prospective AAA patient, the logarithm (acting as a positivity constraint) of the corresponding prediction Θ(ξ) will follow a multivariate Gaussian distribution with predicted mean μ logΘ and covariance matrix S log Θ , i.e.
log ΘðξÞ � N ðm log Θ ; S log Θ Þ ¼ pðlog ΘÞ: As we will see in Section 3.2, our approach leads to more accurate estimates for Θ and also a lower variance in the predictions. All relevant details regarding this model are provided in Appendix A.1.

A Kriging surrogate model for the maximum stress
2.3.1 Estimating the probability of rupture. Since the calculation of the probability of rupture P rupt from Eq (7) using the high-fidelity, nonlinear finite element model from Section 2.2.2 is infeasible for a clinical application, we propose a Kriging surrogate model to speed up computations [40][41][42]. The surrogate model will effectively serve as a proxy for the maximum von Mises stress s max vm ðθÞ in the AAA vessel wall (cf. Eq (16)) and allows to make computationally cheap predictions at an arbitrary combination of θ = [t, α, β] T  To that end, we adopt and extend the Active Learning-MacKay (ALM) strategy from [43] and choose points for high-fidelity model evaluations such as to minimize a density-and stress-weighted predictive standard deviation objective function where p(log Θ) is the patient-specific probability distribution for the invasive model parameters Θ = [t, α, β, σ γ ] T from the regression model in Section 2.2.4. The reasoning behind this choice follows from the ALM approach, where only the predictive standard deviations d log s max vm ðθÞ are considered in the objective function. In our case, we are equipped with a probability distribution, p(log Θ), so we can attribute a higher weight to the more probable regions in Θ. Additionally, we pay special attention to points in the input space, where the predicted maximum von Mises stresses m log s max vm are high to ensure the surrogate model accurately replicates the full model in these regions. The problem of choosing an appropriate point θ next for evaluation results in the optimization problem which is approximated by creating a grid fΘ i g n grid i¼1 over the input space, calculating fcðΘ i Þg n grid i¼1 using the Kriging surrogate and determining The next evaluation point θ next = [t next , α next , β next ] T can then simply be extracted from Θ next . During the active learning, we monitor the averagê and stop the training process, when there are no more significant changes inĉ with an increasing number of high-fidelity model evaluations.

Framework summary
Based on our retrospective AAA database of non-invasive and invasive data pairs and a multioutput Gaussian process model fitted to this dataset (cf. Section 2.2.4), the necessary steps to estimate the probability of rupture for a prospective patient are: • Step 1: Data generation in the clinic: CT imaging, determination of the non-invasive parameters ξ from Table 2 • Step 2: Geometry creation: segmentation and meshing of the AAA geometry While CT imaging is essential for geometry creation, the rupture risk analysis can also be carried out if no non-invasive properties ξ are available for a prospective patient by using cohort statistics (cf. Model 1, Section 3.2) without personalization. The computational procedure is summarized in Algorithm 1. In practice, it has proven feasible to choose n init = 8 (where it makes sense to include the predicted mean μ log Θ in the set of initial samples), n grid = n eval = 10, 000 and tol = 1.0 × 10 −4 .

Regression model benchmark
Before demonstrating the framework in full detail, a brief comparison between the multi-output Gaussian process regression model (cf. Section 2.2.4) with existing probabilistic modeling approaches used in the context of AAA rupture risk is provided. To that end, we employ leaveone-out-cross-validation (LOOCV) on our dataset D (cf. Section 2.2.3) to test the predictive capabilities of three different models for p(logΘ): • Model 1: assuming all variables are log-normally distributed and independent, the joint distribution is obtained, where the means and variances are calculated across the whole population using the dataset D, that is with κ 2 {t, α, β, σ γ }. This corresponds to the approach chosen in [21].
• Model 2: by training single-output Gaussian processes for each output variable separately following [15], the same decomposition of Gaussian distributions as in Eq (25) is obtained, however, with means and variances predicted individually for each patient.
In addition to the mean of the patient standardized mean square error (PSMSE) [15], we also report the mean of the patient predictive entropy (PPE), where H½pðlog ΘÞ� is the entropy of the distribution p(logΘ) and a measure of uncertainty or variance for multivariate distributions. With regards to the different measures, it is desirable for both PSMSE and PPE to be small, corresponding to a model which is accurate and produces low-variance estimates. For conciseness, values for the mean of the PSMSE are averaged over the four predictive variables Θ. We refer to [15] for an exhaustive discussion of the LOOCV and calculation of the PSMSE. The obtained results for the three models are shown in Table 3. We note that our proposed model (Model 3) was able to consistently achieve the lowest scores, although the differences are rather small.

Framework demonstration for AAA Pat17
To illustrate the application of our proposed framework we demonstrate all steps in detail below, following the outline as presented in Section 3.1. We assume we are provided with CT imaging data and non-invasive properties ξ for one specific prospective AAA (Step 1), referred to as Pat17 in the following.  Table 4 shows the relevant 8 non-invasive properties ξ that are used by the regression model (cf. Section 2.2.4) to obtain the predictive distribution p(logΘ(ξ)), which is specific to averaged over the four predictive variables Θ as well as the mean of the patient predictive entropy (E½PPE�). https://doi.org/10.1371/journal.pone.0242097.t003

Fig 2. AAA Pat17 as seen via CT imaging (I), a 3D rendering of the segmentation result (II), the generated finite element mesh (III) and a visualization of the von Mises stress field corresponding to the mean μ logΘ of the predictive distribution p(logΘ) for that AAA (IV).
https://doi.org/10.1371/journal.pone.0242097.g002 Pat17. Along with that, means and standard deviations based on all 113 patients in D are provided. Based on this data, we can predict the mean μ logΘ and covariance S logΘ for this patient ( Step 3). The obtained distribution is visualized in Fig 3 and the predictive means and standard deviations are provided in Table 5 along with reference values from the cohort. The entropy of p(logΘ) is 3.3050 and thus slightly lower than the LOOCV mean (cf. Table 3). Highest correlations among the invasive properties for Pat17 can be found between t and σ γ (r t;s g ¼ À 0:3214), β and σ γ (r b;s g ¼ 0:2274), t and β (ρ t,β = −0.1966) as well as α and β (ρ α,β = 0.1413).  Given p(logΘ), the forward model (cf. Eq (16)) for Pat17 is defined. The probability of rupture for this AAA is approximated using a Kriging surrogate model (Step 4). Lastly, the probability of rupture can be estimated using the Kriging surrogate (Step 5), which amounts to 0.47% for Pat17 (cf. Fig 5 for a visualization). We stress that this value must not be compared to the operative risks associated with OSR or EVAR in order to make decisions. Rather, it needs to be put into context with results for other AAA patients that have been computed using the same methodology, which is discussed below in Section 3.4.

Comparative case-control study using diameter matched groups
To test the efficacy of the framework as a rupture risk indicator and to compare it with existing biomechanical indices, we consider diameter matched groups of asymptomatic (group 1, n = 18) and known symptomatic/ruptured (group 2, n = 18) AAA patients from our database. The groups were chosen such that their maximum diameter mean and standard deviation approximately match (group 1: 62.17±7.18 mm, group 2: 63.06±7.56 mm), rendering a differentiation between the groups based on the maximum diameter criterion ineffective.
For a detailed overview regarding the selection of the two groups, we refer to Table 6. After preprocessing of our original dataset (cf. Section 2.2.3), we restricted the cohort to AAAs with a maximum diameter between 50 and 80 mm in order to obtain an intermediatesized group of patients. As a result, 64 patients remained, of which 47 had asymptomatic and 17 had symptomatic or ruptured AAAs. The latter were put into one group, since symptomatic AAAs are known to be at an elevated risk of rupture [44]. The reason for the much lower number of symptomatic/ruptured AAAs is that these AAAs often have very large diameters (>80mm). We included AAA patients from a previous case-control study by our group [19], which examined 13 asymptomatic and 12 symptomatic AAA patients. Finally, we manually selected 18 asymptomatic and 18 symptomatic/ruptured patients based on the following criteria:

PLOS ONE
• Find two groups with the best match in diameter.
• Preferably include cases where non-invasive data is available and thus patient-specific invasive properties can be predicted.
• Disregard cases, where CT images are not available or lack a sufficient image quality to create simulation models.
Detailed information for all AAAs of both groups is provided in Tables 7 and 8 and a visualization of their rupture risk indices, P rupt , in Figs 5 and 6 (cf. Appendix A.3). No patient had known connective tissue disorders. For 10 out of 18 AAAs in group 1 and for 9 out of 18 AAAs in group 2 we had non-invasive data and were thus able to use the multi-output regression model to determine a personalized input density p(logΘ). For the remaining 8 (group 1) and 9 (group 2) AAAs, we used cohort statistics, i.e. Model 1 from Section 3.2. We apply our framework to all 36 AAAs using an individual prospective scenario, i.e. before starting the analysis for one AAA, this patient is removed from the database, while the other 35 AAAs are included. In order to provide a comparison of P rupt with other biomechanical indices, we calculate the following additional quantities: • Maximum von Mises stress at the input parameter mean (neglects any statistical information): • Rupture potential index [28] at the input parameter mean (neglects any statistical information, but takes into account the wall strength):

PLOS ONE
• Probabilistic rupture risk index [21] (takes into account cohort-based uncertainties in the wall thickness and wall strength according to Model 1, Section 3.2): Comprehensive results for all patients are listed in Tables 7 and 8 (cf. Appendix A.3). The average number of high-fidelity model evaluations to train the Kriging surrogate was 11. Based on these results and to evaluate the performance of the individual quantities, we provide: 1. Relative mean and median differences between group 1 and group 2 (cf. Table 9).
3. Receiver operating characteristic (ROC) curves and the area under the ROC curve (AUC) (cf. Fig 8) [45]. Computed true positive rates (TPR), false positive rates (FPR) and corresponding threshold values are provided for P rupt as supplementary information (cf. S4 Table). Relative differences for a quantity q between the asymptomatic group result q a and the symptomatic/ruptured group result q s/r are calculated as Δq = |q s/r − q a |/q a . https://doi.org/10.1371/journal.pone.0242097.t009

Discussion
The obtained values for the relative mean and median differences in Table 9 confirm that group 1 and group 2 are indistinguishable based on the maximum diameter criterion. While the relative differences are higher for s vm max and RPI, PRRI and in particular our proposed index P rupt feature a significantly larger mean and median difference. Recalling that the maximum diameter, d max , is one important non-invasive parameter in our framework (cf. Section 2.2.3), we emphasize that its influence has been rendered ineffective through the study design. A similar trend as in Table 9 can be observed in Fig 7, with RPI and PRRI providing a slightly better separation between the two groups than s vm max , while for P rupt the interquartile ranges of the two groups are non-overlapping. Finally, in Fig 8 we can observe that P rupt outperforms the remaining classifiers and achieves the best performance among all quantities in terms of the AUC score, followed by PRRI, RPI and s vm max . We further note that from the 18 patients in the symptomatic/ruptured group, 11 had ruptured AAAs (Pat19, Pat23, Pat24, Pat26, Pat27, Pat28, Pat29, Pat30, Pat32, Pat34, Pat35). The mean P rupt scores for the 11 ruptured AAAs is 6.57 and thus slightly higher than the mean 5.89 for the 7 symptomatic AAAs. To summarize our key observations: • The maximum diameter criterion, by design, clearly fails to separate the two groups in all our comparisons.
• The proposed index P rupt consistently achieves the best separation.
• The results indicate that the more statistical information taken into account, the better the capability to distinguish between group 1 and group 2.
Before translating these findings into any clinical application, however, there are several limitations that have to be kept in mind. First, this is a non-randomized, retrospective casecontrol study with a relatively small cohort size (group 1: n = 18, group 2: n = 18) and the database described in Section 2.2.3. Second, there was no matching based on other risk factors such as sex, age or family history, which could be a confounder. Third, since we only have access to electively repaired or symptomatic/ruptured AAAs for mechanical testing, the mean diameters of the two groups (group 1: 62.17 mm, group 2: 63.03 mm) are larger than the Society for Vascular Surgery's decision criterion for elective repair (55 mm) [1]. In the future, due to the increasing use of EVAR, it will be even harder to obtain representative tissue samples from AAAs of relevant size for a database. As a result, caution is advised when interpreting the results presented here for smaller AAAs, e.g. of size 45 − 55 mm. Furthermore, all discussed approaches are unable to make any prediction about the future development of the AAA, such that the rupture risk assessment only holds for the point in time of data generation. In addition to that, the biomechanical model does not take into account factors like growth, calcifications and surrounding organs, which might be important for the analysis.

Conclusion
We presented a novel data-informed, highly personalized, probabilistic framework for the quantification of abdominal aortic aneurysm (AAA) rupture risk and demonstrated competitive performance in comparison to existing approaches. Our framework results in the calculation of a rupture risk index, P rupt , which can be introduced as a relevant additional piece of information in the clinical decision process for AAA cases that are not or not unambiguously covered by existing guidelines and recommendations. In view of our results it is suggested to incorporate personalized, or at least cohort-based, statistical information and choose a probabilistic approach for the biomechanical rupture risk assessment. Deterministic indices were shown to be less accurate and do not account for possible sensitivities due to uncertain vessel wall quantities.
In order to advance this framework to a clinical application, several further aspects need to be examined. Challenges lie especially in the fully automatic segmentation of the CT imaging data, which at the moment requires manual steps by a trained expert and can be time consuming. In view of the limitations discussed in Section 4, a larger, randomized study with risk factor matched groups is desirable to confirm this study's findings regarding its clinical use. Future work will also address how further model parameters such as the blood pressure influence the rupture risk index and whether this quantity should be treated probabilistically as well. Lastly, to be able to make predictions over time, it is required to incorporate AAA growth [46,47] into the framework and analyze its effect on the biomechanical rupture risk assessment.

A.1 Multi-output Gaussian process regression
The data generation process for logΘ is assumed to underly a function logΘ ξ ð Þ that is contaminated by additive Gaussian noise, such that It is further postulated that the vector logΘ ¼ ½logt; logã; logb; logs g � T � N 0; O ð Þ follows a multivariate Gaussian distribution with the positive semi-definite covariance matrix O 2 R 4�4 and it is assumed that 2 � N ð0; SÞ, with the diagonal matrix S and noise levels S dd 2 R þ ðd ¼ 1 . . . 4Þ. Demanding that every entry of the vector logΘðξÞ corresponds to the same zero mean Gaussian process with covariance function k(ξ, ξ 0 ), i.e. logtðξÞ; logãðξÞ; logbðξÞ; logs g ðξÞ � GPð0; kðξ; ξ 0 ÞÞ; ð32Þ logΘðξÞ can be expressed as a multivariate Gaussian process [39] logΘðξÞ � MGPð0; kðξ; ξ 0 Þ; OÞ: ð33Þ As a result, the collection flogΘ i g n data i¼1 follows a matrix-variate Gaussian distribution ½logΘ 1 ; . . . ; logΘ n data � T � MN ð0; K; OÞ; ð34Þ with the covariance matrix K and entries K ij = k(ξ i , ξ j ), modeling the covariance between two inputs ξ i and ξ j . Expressing the matrix Gaussian distribution as a multivariate Gaussian distribution and incorporating the additive noise �, one obtains where � denotes the Kronecker product and I n data the n data × n data identity matrix. For our purposes, we choose the covariance function with hyperparameters z 1 , z 2 , z 3 and z 4 . Following [37][38][39], and the predicted covariance where k ? denotes the vector of covariance function evaluations between ξ ? and the data fξ i g n data i¼1 , i.e. k ? ¼ ½kðξ ? ; ξ 1 Þ; . . . ; kðξ ? ; ξ n data Þ� T : Finally, the log marginal likelihood is and can be optimized with respect to its hyperparameters z.

A.2 Kriging surrogate incorporating explicit basis functions
Kriging can be regarded as a special case of a Gaussian process, where data points are assumed noise-free to interpolate the high-fidelity model at the provided high-fidelity evaluations. To find an adequate function for this purpose, we use a Kriging interpolation model that incorporates explicit basis functions as described in [42]. Using such a model, it is possible to exactly represent functions that can be described by the provided basis. Ensuring positive predictions via a log transformation, we approximate the high-fidelity model as log s max vm ðθÞ � logs max vm ðθÞ þ hðθÞ T Z; ð42Þ where logs max vm ðθÞ � GP 0; kðθ; θ 0 Þ ð Þ is a zero mean Gaussian process with covariance function k(θ, θ 0 ) and h(θ) denotes the chosen basis functions with coefficients η. A simple squared exponential kernel is chosen, where the matrix L ¼ diagðz 2 ; z 3 ; z 4 Þ 2 R 3�3 is diagonal, leading to the vector of hyperparameters z = [z 1 , z 2 , z 3 , z 4 ] T , with z 2 R 4 þ . Furthermore, trilinear basis functions, i.e.  (46) and (47) consist of a contribution from the zero mean Gaussian process predicted mean and variance, m logs max vm ðθ ? Þ and d 2 logs max vm ðθ ? Þ, and additional terms involving the provided basis functions, respectively. Finally, the marginal log likelihood is L z ð Þ ¼ log pðσjθÞ ¼ À 1 2σ T K À 1σ þ 1 2σ T Cσ À 1 2 log jKj À 1 2 log jAj À n eval À m 2 log 2p; ð50Þ where A = HK −1 H T C = K −1 H T A −1 HK −1 and m is the rank of H T . Maximizing the marginal log likelihood, optimal values for the hyperparameters z of the Kriging covariance model (cf. Eq (43)) and for the providedσ andθ can be determined.
Supporting information S1