Machine learning and reduced order modelling for the simulation of braided stent deployment

Endoluminal reconstruction using flow diverters represents a novel paradigm for the minimally invasive treatment of intracranial aneurysms. The configuration assumed by these very dense braided stents once deployed within the parent vessel is not easily predictable and medical volumetric images alone may be insufficient to plan the treatment satisfactorily. Therefore, here we propose a fast and accurate machine learning and reduced order modelling framework, based on finite element simulations, to assist practitioners in the planning and interventional stages. It consists of a first classification step to determine a priori whether a simulation will be successful (good conformity between stent and vessel) or not from a clinical perspective, followed by a regression step that provides an approximated solution of the deployed stent configuration. The latter is achieved using a non-intrusive reduced order modelling scheme that combines the proper orthogonal decomposition algorithm and Gaussian process regression. The workflow was validated on an idealized intracranial artery with a saccular aneurysm and the effect of six geometrical and surgical parameters on the outcome of stent deployment was studied. We trained six machine learning models on a dataset of varying size and obtained classifiers with up to 95% accuracy in predicting the deployment outcome. The support vector machine model outperformed the others when considering a small dataset of 50 training cases, with an accuracy of 93% and a specificity of 97%. On the other hand, real-time predictions of the stent deployed configuration were achieved with an average validation error between predicted and high-fidelity results never greater than the spatial resolution of 3D rotational angiography, the imaging technique with the best spatial resolution (0.15 mm). Such accurate predictions can be reached even with a small database of 47 simulations: by increasing the training simulations to 147, the average prediction error is reduced to 0.07 mm. These results are promising as they demonstrate the ability of these techniques to achieve simulations within a few milliseconds while retaining the mechanical realism and predictability of the stent deployed configuration.


Introduction
Intracranial aneurysms (IAs) are local dilations of the arteries in the brain caused by a degenerative weakening of the arterial wall. Saccular, or blister-like, are the most common IAs. Their prevalence among the general population is estimated to be around 2%-3% (Rinkel et al., 1998). With an incidence of 10/100,000 person-years, IAs rupture leads to subarachnoid haemorrhage, a life-threatening type of stroke with high morbidity and mortality (Wiebers, 2003;Vlak et al., 2011). Therefore, when IAs with a diameter larger than 5/7 mm are detected, they are often recommended for early treatment to exclude the aneurysm sac from the cerebral circulation (Rahman et al., 2010). Nowadays, endovascular options have become the preferred intervention thanks to the lower rate of complications with respect to invasive techniques (e.g., clipping) (Liu and Huang, 2015).
Endoluminal reconstruction through flow diverters represents a novel paradigm for the mini-invasive treatment of IAs, as an alternative to endosaccular occlusion through coiling (Durso et al., 2011;Pierot, 2011). Flow diverters are self-expanding devices consisting of a low-porosity braided stent. Thanks to this structure, they are highly flexible and resistant to kinking. Accordingly, they are well suited for tortuous vessels and wide-neck IAs. These stents come in different sizes, in terms of radius and length. Nowadays, surgeons choose the size based only on clinical experience and measurements taken on medical volumetric images (e.g., 3D rotational angiography or computed tomography angiography), acquired shortly before surgery.
However, the configuration assumed by these devices once deployed within the parent vessel is not easily predictable and routine 3D medical images alone may be insufficient to plan the treatment satisfactorily. The related difficulties can lengthen the interventional times. Moreover, higher doses of angiographic contrast agents and anaesthetic drugs are needed. All these factors increase the risk of postoperative complications for the patient. Therefore, there is a compelling need to develop computational models capable of simulating, in real time, the deployment of flow diverters within patient-specific vessels to assist practitioners in the planning and interventional stages (Karmonik et al., 2005).
The mechanical behaviour of braided stents is typically simulated using a finite element (FE) model where beam elements are used to discretise the wires (Auricchio et al., 2011;Shiozaki et al., 2020;Zaccaria et al., 2020;McKenna and Vaughan, 2021); however, only a few studies modelled numerically flow diverters (Ma et al., 2012;Fu and Xia, 2017). Due to the large amount of degrees of freedom (DOFs) and the necessity to solve the contact between the device and the vessel wall, the computational time required by these traditional techniques is very high. To overcome this limitation and make computational models suitable for clinical use in real time, fast virtual stenting (FVS) methods have been reported in the literature (Larrabide et al., 2012;Spranger et al., 2015;Zhong et al., 2016). They predict the stent deployed configuration by simplifying its mechanical behaviour and/or the contact against the vessel wall (e.g., by using simplex deformable meshes, mass-spring models or active contour models).
Reduced order modelling is also gaining interest in computeraided surgery thanks to its capability of reducing the computational complexity and cost of numerical problems while preserving their inner physics (Niroomandi et al., 2012a;Niroomandi et al., 2012b;Mena et al., 2015;Santo et al., 2020). One of the most powerful and widespread techniques to build reduced-order models (ROMs) is the reduced-basis (RB) method. The RB method is adapted to non-linear problems which need to be solved a large number of times for different parameter values (Hesthaven et al., 2016;Quarteroni et al., 2016). In biomechanics, the parametrization can concern boundary and initial conditions (Chang et al., 2017;Bridio et al., 2022;Girfoglio et al., 2022), loads (Biancolini and Valentini, 2018) and the geometrical domain under investigation (Ballarin et al., 2016;Biancolini et al., 2020;Kardampiki et al., 2022). The latter is sometimes handled with a statistical shape model, which enables the comparison of different anatomies in terms of the same characteristics (Lauzeral et al., 2019;Cosentino et al., 2020;Buoso et al., 2021).
The real-time simulation of such parametrized problems is achieved by exploiting the intrinsic similarities between their solutions. The most important features of the original, full-order model (FOM), i.e., the RBs, can be extracted through proper orthogonal decomposition (POD) from a set of high-fidelity (HF) solutions. The construction of such a dataset and the extraction of the RBs is referred to as offline stage. Thereafter, approximated solutions for unseen parameter values are determined as linear combination of the RBs (online stage). The powerful advantage of these methods is that, if the online stage is completely decoupled from the offline one, the computations performed in the online stage are independent of the dimension of the FOM. RBs-based methods differ in the implementation of the online stage: Intrusive methods rely on a projection onto the RBs space to generate the ROM (Ballarin et al., 2015;Ballarin et al., 2016); non-intrusive methods employ a regression model trained to learn the mapping from parameters to the solution expressed in the RBs space on the HF dataset (Guo and Hesthaven, 2018;Hesthaven and Ubbiali, 2018;Guo and Hesthaven, 2019;Fresca and Manzoni, 2022). Non-intrusive methods outperform intrusive methods in terms of efficiency, as they do not require solving a system of non-linear equations in the online phase, but only evaluating the regression model. However, they require a very large training dataset to ensure accurate solutions. Recently, physics-informed neural networks (PINNs) emerged as a promising alternative (Liu et al., 2020;Buoso et al., 2021;Chen et al., 2021).
In this study, we present the first implementation of a machine learning (ML)-based ROM scheme for the prediction of the stent deployed configuration. To assess its feasibility, we created a parametric synthetic geometry that allows control of the vessel radius, curvature and aneurysm size. Surgical decisions on the deployment site are also considered when creating the HF dataset. If given the deployment conditions the stent does not land inside the aneurysm and is well positioned against the vessel wall, the deployment is considered successful from a clinical perspective. Since there is no clinical need to predict the stent configuration after an unsuccessful deployment, the virtual framework here proposed consists of two steps (Figure 1): A first classification step that allows a priori determination of whether a simulation will be successful or not, followed by a regression step that provides an approximated solution of the deployed stent configuration. Moreover, in continuation with our previous work (Bisighini et al., 2022), we propose a fast strategy to perform FE simulations of braided stent deployment to reduce the computational time required to build the HF dataset for training.
The remainder of the paper is organized as follows. Section 2.1 summarizes the methods used to generate the HF dataset with particular emphasis on the description of the flow diverter model (Section 2.1.1), the synthetic aneurysm model (Section 2.1.2) and the scheme followed to perform FE simulations of stent deployment (Section 2.1.3). Section 2.2 describes the classification model, Section 2.3 the ML-based ROM. In Sections 3, 4, the results of testing these models against unseen scenarios are shown and discussed. Finally, Section 5 presents some concluding remarks.

Braided stent modelling
The braided stent is modelled as a tubular net of thin wires with circular cross-section (Figure 2A). For a stent with radius R s , length L s , composed by N w wires with radius R w and presenting N cells repetitive units, the nodal positions are defined by the following set of equations: with n ∈ [0, N w /2] and i ∈ [0, N cells ] and where orient is either 1 or -1 respectively for clockwise and counter-clockwise wires, dθ = 2π/(N w /2), Θ n = n ⋅ 2π/(N w /2) and ϕ = L s /N cells is the pitch angle.
In this work, the values assigned to these geometrical parameters refer to a generic flow diverter: N w = 48, R s = 2.6 mm, R w = 0.014 mm, L = 15 mm, N cells = 70. The stent is made of Phynox, a cobalt-chromium alloy, which is modelled as a linear elastic material with Young modulus E = 225 GPa, Poisson coefficient ν = 0.33 and density ρ = 9.13 ⋅10 3 kg/m 3 .

Artery and aneurysm modelling
The stent is released within a parametric idealised model of an intracranial artery presenting a saccular aneurysm. The vessel centerline is defined using a planar quadratic Bézier curve: where P 0 , P 2 are fixed and P 1 is the control point which will be included in the parametrisation. These points lie in a 2D plane, so the only DOFs of this spline are the y, z coordinates of P 1 . The curve B(t) starts from P 0 in the direction of P 1 , then bends to reach P 2 . Bézier curves allow defining smooth, continuous curves that resemble the curvature of intracranial arteries (Zyłkowski et al., 2018;Danu et al., 2019). The Visualization Toolkit (VTK) software is employed to build the model (Schroeder et al., 2006). The artery is created using the vtkTubeFilter() function that generates a tube around a line. Its diameter (D v ) is considered constant along the centerline. Following, a spherical idealised aneurysm with diameter D a is created using vtkSphereSource(), the sphere centre C a is positioned in the middle of the vessel centerline and the relative y-distance is parametrised. Through a Boolean union, the sphere is added to the artery model ( Figure 3A). In summary, the vessel geometry is fully parametrised by 5 parameters: y P 1 , z P 1 , D v , D a , y C a .

FE simulation of braided stent deployment
The simulations are performed using EndoBeams.jl, an inhouse and open-source FE modelling framework for the numerical simulation of frictional contact interactions between beams and rigid surfaces (Bisighini et al., 2022). This software has been validated against some public benchmark tests and the commercial software Abaqus (Simulia, Dassault Systems, Providence, RI, United States).
The stent structure is discretised using beam elements and the resulting mesh is composed of 3,408 nodes and 3,360 beam elements. To avoid rendering beam-to-beam contacts, a penaltybased constraint is imposed at each interconnection between the interlaced wires so that the position of nodes pairs is forced to be the same (Zaccaria et al., 2020). The contact algorithm implemented in EndoBeams.jl relies on an implicit representation of the contact surface, the signed distance field (SDF) (Aguirre and Avril, 2020). The SDF of the artery model is computed using the Julia package SignedDistanceField.jl (SignedDistanceField.jl, 2022), which only requires as input the triangular mesh in the form of a. stl file ( Figure 3B). The vessel wall is assumed rigid, thus the SDF is constant along the simulation. The normal contact pressure is treated using a viscoelastic model evaluated by means of a regularised penalty method while the tangential force is determined with a regularised Coulomb friction law considering slip-stick behaviour (Wriggers, 2006;Aguirre and Avril, 2020;Shiozaki et al., 2020).
The strategy followed to perform simulations of stent deployment was inspired by the work proposed in (Perrin et al., 2015;Spranger et al., 2015;Hemmler et al., 2018). This approach is compatible with the use of the SDF to manage the vessel-stent contacts and represents an alternative to the use of a virtual catheter, as typically done before in the literature (Bock et al., 2012). The simulations are carried out in three steps: 1) crimping; 2) positioning and 3) deployment. Since our goal is to obtain the stent configuration at the end of the deployment simulation, we focus on quasi-static simulations. Moreover, as commonly done in literature, the stent is assumed to be positioned along the vessel centerline at the beginning of the deployment (Auricchio et al., 2011;Bock et al., 2012;Hemmler et al., 2018;Leng et al., 2018).
First, the braided stent is crimped by blocking its circumferential DOFs and imposing a radial displacement equal to (R stent − R crimped ) to all its nodes, where R crimped is the stent target radius after the crimping (Figure 2A).
The braided stent nodes lying on the same z-plane can be grouped in N r rings and a central point can be introduced for each of these rings, R i with i = 1 … N r ( Figure 2B). These points define the initially straight centerline of the crimped stent (C 0 ). The stent is displaced so that R 0 coincides with the chosen deployment site along the vessel centerline. The final centerline of the stent (C T ) is computed as the projection of C 0 along the vessel one ( Figure 2C) so as to maintain constant the total length of the stent centerline. C 0 , and so C T , can be subdivided into segments, i.e., vectors connecting subsequent points along the centerline (seg i = ‖R i − R i−1 ‖). The first point R 0 is considered aligned with the z-axis. A rotation is associated with each of these segments to realign it to the preceding one. The angle (θ i ) and axis (ax i ) of rotation for each segment can be computed as follows: Thus, a rotation matrix M i can be defined for each segment. By cumulatively applying these rotations, we can align C T with the z-axis, obtaining C 0 : The advantage of this technique is the possibility of obtaining intermediate positions (C t ) between C T and C 0 by dividing θ i by the number of desired configurations and applying only this partial rotation to C T ( Figure 2D). These intermediate positions are used to drive a physical simulation where kinematics constraints are implemented between all the stent nodes lying on the same ring and their corresponding points R i on C 0 ; the difference of the coordinates between C t and C 0 is computed and the resulting displacements are applied to each point of C 0 in a successive way leading the stent to bend ( Figure 2E). Finally, the stent is allowed to freely deform within the vessel and activate the contact against the wall ( Figure 2F). The simulations are stopped when the kinetic energy falls below a certain threshold (10 -12 mJ), which represents static mechanical equilibrium.

Creation of the high-fidelity dataset
The impact of anatomical characteristics and surgical decisions on the final stent configuration is studied. Therefore, we consider a set of geometric features describing both the artery and aneurysm geometry and the stent deployment site along the vessel centerline; we refer to this generic set of parameters as simulation parameters, collected in the vector μ.
For the creation of the HF dataset, we considered the following simulation parameters: where y P 1 and z P 1 are the y and z coordinates of the middle Bézier curve point, y C a the y coordinate of the aneurysm centre point and η the stent position along the vessel centerline. Since the impact of stent misplacement is one of the objectives of this study, we consider deployment sites in which one of the ends of the stent falls in the aneurysm neck area. A Latin hypercube sampling (LHS) method is used to generate N s different values for μ B ; the corresponding artery models are created and a stent deployment simulation is performed within each model as explained in Section 2.1.3. The LHS plan is created using the Julia package LatinHypercubeSampling.jl (LatinHypercubeSampling.jl, 2020). The simulation parameters are evaluated within a range that resembles that observed in the literature (Krejza et al., 2006;Li et al., 2013): for D v , we consider a range of [2,4] mm; for D a , a range of [5,10] mm.
Alternatively, the simulation parameters μ B are substituted as input of the ML models with μ cl where, instead of the middle point of the Bézier curve (P 1 ) and the stent deployment site (η), the y and z coordinates of N cl points (Q i ) on the positioned stent centerline (C T ) are used: The μ cl vector components are calculated geometrically from the deployment conditions and no FE simulation is required.

Binary classification
Within the HF dataset, deployment solutions are considered "successful" from a clinical perspective if the stent extremities are in contact with the artery wall and do not land within the aneurysm sac; otherwise, the simulation outcome is labelled as "failure" (Figure 4). This classification is done automatically by checking if, at the end of the simulation, one or more nodes from the stent extremities are inside the aneurysm. Being it modelled as a sphere, the SDF of the aneurysm can be computed analytically and, thus, a node x p is located inside the aneurysm if its distance to the surface is positive, i.e.,: Once the results are labelled, a supervised machine-learning algorithm is trained to learn the relationship between the simulation parameters μ and the corresponding solution output: when the model is built, it allows the assignment of new, unseen scenarios to one of the two categories. The classification is done in Matlab (MathWorks, Natick, MA, United States). For this purpose, the performance of the following classifiers is compared (Singh et al., 2016): • A logistic regression (LR) model is created and fitted using the fitctree() function; • A k-Nearest Neighbour (k-NN) model is created and fitted using the fitcknn() function; • A naive Bayes (NB) model is created and fitted using the fitcnb() function; • A decision tree (DT) model is created and fitted using the fitctree() function; • An artificial neural network (ANN) model with three layers of size [10, 10, 10] and the hyperbolic function tanh as activation function is created and fitted using the fitnn() function; • A support vector machine (SVM) model with a polynomial kernel of order 2 is created and fitted using the fitsvm() function.
The architecture and hyperparameters values are chosen for each ML model using the ClassifierLearnerApp.

Metrics
Five metrics are considered to evaluate the performance of the trained classifiers. They are computed by counting true positives (TPs), true negatives (TNs), false positives (FPs) and false negatives (FNs), which are collected in the confusion matrix. Their expressions are.
• Accuracy = TP+TN TP+FN+TN+FP : it tells how good is the classifier, regardless of the label meaning; • Sensitivity (or recall or TP rate) = TP TP+FN : it tells how good the classifier is at predicting successful deployment conditions; • Specificity (or TN rate) = TN TN+FP : it tells how good the classifier is at predicting failed deployment conditions; • Precision = TP TP+FP : it tells how close predicted values are to each other; • F1-score = 2 precision⋅recall precision+recall : it is a more general measure of accuracy that combines precision and recall in a single metric.
Frontiers in Physiology 05 frontiersin.org  Another useful tool to asses the classification performance is the receiver operating characteristic (ROC), a plot showing the performance of the classifier in terms of TP rate (sensitivity) against FP rate (1-specificity) as a function of the cut-off threshold. The metric related to the ROC curve is the area under the curve (AUC). The closer the ROC curve is to the upper left corner of the graph (and thus the higher the AUC value), the more accurate the classifier is.

Reduced order modelling
In case of successful deployment, a ML-based reduced order modelling method is employed to compute an approximated solution for the stent deployed configuration within the considered vessel. At first, one might consider training a ML algorithm to learn the relationship between the simulation parameters and the vector of nodal displacements at the end of the simulation. However, given the large number of DOFs of the stent model (in our case, ∼10,000 DOFs), the size of such an output vector is very large and would result in a very long training time. Reduced order modelling allows for the reduction of the original problem dimension by extracting the most important features, the RBs, from the training dataset through POD (Guo and Hesthaven, 2018;Han et al., 2020;Bridio et al., 2022;Fresca and Manzoni, 2022). A supervised ML algorithm is then used to establish the relationship between simulation parameters and the solution expressed in the RBs space. Finally, an approximated solution of the stent deployed configuration can be recovered in real time for any combination of simulation parameters.

Proper orthogonal decomposition
The theory behind the POD algorithm is here only introduced: if interested, the reader is suggested to refer to (Chatterjee, 2000). Once the HF dataset is computed, the snapshots matrix S is built by arranging the HF solutions u h (μ) as columns of a matrix: Each of these vectors represents one snapshot and contains the nodal displacements at the end of the quasi-static deployment simulation: where N n is the number of nodes in the stent mesh, thus the dimension of the HF problem is N h = 3 ⋅ N n . The POD algorithm relies on performing the singular value decomposition (SVD) of S: where U = [ u 1 |u 2 |⋯|u N h ] is the left singular vectors matrix, Z = [ z 1 |z 2 |⋯|z N s ] is the right singular vectors matrix and Σ = diag(σ 1 , σ 2 , …, σ N s ) contains the singular values of S, sorted from the largest to the smallest (σ 1 ≥ σ 2 ≥ ⋯ ≥ σ N s ≥ 0). The Schmidt-Eckhart-Young theorem states that the columns of S, Col(S), can be well approximated by the first L left singular vectors of S, i.e., Col(U), if the singular values decay rapidly. Thus, given a tolerance ɛ POD , L can be found as the minimum integer such that: The column vectors [ u 1 |u 2 |⋯|u L ] represent the RBs of the model and are assembled in the matrix V.
The HF solution u h (μ) can be now projected onto the reduced space defined by V: where UU T = I since U is orthogonal, u L are the L projection coefficients associated with the column bases of V and u rb (μ) is the solution projected onto the reduced space.

Gaussian process regression
As proposed by Guo in (Guo and Hesthaven, 2018), Gaussian process regression (GPR) is adopted to approximate the HF solutions for any simulation parameters combination. A theoretical introduction to GPR is provided in Appendix 1, here we simply introduce the implementation aspects and variables necessary for a complete interpretation of the results.
In this work, a GPR modelf is constructed using the Matlab function fitrgp() with the Matérn 5/2 kernel from a set of training data where the predictors are the simulation parameters μ and the outputs are the projection coefficients u L (μ) computed from the HF solutions as: Once the model is trained, the function predictExactWithCov() is used to predict the projection coefficients for any desired unseen value of the simulation parameters μ*, i.e.,û L (μ * ). This allows recovering the full-order solution as follows:

Metrics
For the test cases, the ML-based ROM results are evaluated against the HF solutions by computing three absolute errors at each node of the stent mesh: • The order reduction error, E rb = ‖u rb (μ) − u h (μ)‖; • The prediction error, E p = ‖u p (μ) − u h (μ)‖; • The GPR error, E gpr = E p − E rb .
At the nodal level, these errors correspond to the distance between pairs of nodes. The accuracy of a single solution can be evaluated as the average error (AE) or maximum error (ME) on the mesh nodes: e.g., for the order reduction error, we get AE rb and ME rb . Global variables can be computed by averaging these values on the test dataset: e.g., for the order reduction error, we refer to AE rb as the average of AE rb over the test solutions and to ME rb as the average of ME rb over the test solutions.

Binary classification: Prediction of deployment success
In this section, the results of the classification models are presented. To study the influence of the dataset size on the prediction capability of the ML models, we first created a dataset of N s = 900 simulations and then, by using the subLHCoptim() function from the LatinHypercubeSampling.jl package, we defined three optimal subspaces with N s = {150, 300, 600}. Following an approach similar to (Hesthaven and Ubbiali, 2018), we considered a fixed number of samples (N test = 100) for testing the four different datasets. In Table 1, the number of samples considered for ML training and testing is reported. The input data were standardised before training. The values of the metrics for the four different dataset sizes are given in Table 2. For N train = 200, the confusion matrices of the different classification models are reported in Figure 5. The NB, ANN and SVM models all show an AUC larger than 96% when using N train = 200. The performance of all classification models improves when the training dataset is expanded. In general, the LR model shows the worst performance: regardless of the dataset size, its metrics stay below 90%. Precision is especially low (highest value = 32.3%): as visible in Figure 5A, the model fails mostly in predicting false cases (failures). For the smallest dataset, the SVM model shows the best performance with accuracy, specificity and precision larger than 90%. When considering N train = 500, all the validation metrics are larger than 90% for the SVM and ANN models. k-NN, DT and NB models also show good performance in terms of accuracy and specificity (between 83% and 96%); however, when compared with SVM and ANN, they have poorer precision, which is reflected in a lower F1-score (maximum 83.3%). The ANN and SVM models achieve the highest value of specificity (97%), precision (93.9%) and F1-score (92%); instead, the highest sensitivity is reached by the NB model (95.8%). When considering the largest dataset N train = 800, only a few improvements are observed, e.g., the specificity of the NB model increases from 95.8% to 100%. However, the majority of the metrics stabilise or even decrease in value: most likely, the models are undergoing overfitting. The results of the classification are related to the use of μ B as predictors; no improvement is observed using μ cl for classification.

Reduced order modelling: Prediction of the stent deployed configuration
In this section, the results of the ML-based ROM are presented. As proposed for the classification models, we started with a dataset with N s = 900 simulations and then defined three optimal subspaces with N s = {150, 300, 600}. Then, we built four datasets with only the successful deployment simulations, N success = {97, 196, 396, 555}, and considered a fixed number of samples (N test = 50) for testing. In Table 1, the number of samples considered for ML-based ROM training and testing is reported. The input and output data were standardised before training. The spatial resolution of imaging techniques currently used for IAs detection and treatment was considered for evaluating the prediction errors (Hacein-Bey and Provenzale, 2011;Maupu et al., 2022): magnetic resonance angiography (MRA) = 0.6 mm-1 mm, computed tomography angiography (CTA) = 0.4 mm-0.7 mm, digital subtraction angiography (DSA) = 0.2 mm, 3D rotational angiography (3DRA) = 0.15 mm.
As already mentioned in Section 2.3, the choice of the RBs number L is carried out considering the singular values of the snapshots matrix. In Figure 6A, the cumulative sum of the singular Frontiers in Physiology 09 frontiersin.org values normalised with respect to their total sum is reported: this plot shows the fraction of total variance retained by the first Lsingular values for N train = {47, 146, 346, 555}. As reported in Table 3, regardless of the dataset size, the first 4 RBs cover more than 90% (ɛ POD = 0.1) and the first 10 RBs more than 99% (ɛ POD = 0.01) of the total variance in the input data. As shown in Figure 6B, the prediction does not benefit considerably by considering a number of RBs greater than 15: in fact, analogously for any value of N train , the average order reduction error decreases towards 0 as more RBs are considered, while the average prediction and GPR errors reach a stable plateau around L = 15. Therefore, henceforth L = 15 is considered in this analysis. The results presented so far relate to the use of μ B as predictors. As reported in Table 4 and shown in Figure 7, a strong reduction of the prediction error is evident when using μ cl instead: with N cl = 3, the average prediction error is 5.65× lower than when using μ B while the maximum prediction error is 2.84× lower. Increasing the number of considered points N cl , a further but slight decrease of the prediction error is observed: the average and maximum values are respectively 1.11× and 1.03× lower with N cl = 5 and 1.08× and 1.04× lower with N cl = 8. Using μ cl , the percentage of test samples where the average prediction error is greater than the spatial resolution of 3DRA is zero; the same is true for the number of test samples where the maximum prediction error is greater than the spatial resolution of CTA. The number of test samples where the maximum prediction error is greater than the spatial resolution of 3DRA is reduced from 64% with μ B to 10% with μ cl . The improvement gained by using μ cl is not only reflected in a lower variance of the prediction error within the test dataset but also within a single test solution: this is noticeable when sampling the multivariate normal distributions corresponding to the GPR outputs, as explained in Appendix 1. In Figure 8, we reported the position of 100 points sampled within the distribution predicted by the GPR for the displacement of the first node of the stent mesh. The points are more widely distributed with μ B while they are highly concentrated around the mean node value with μ cl , which means that the uncertainty on the predicted displacement value when using μ cl is lower than when using μ B .
Finally, we analysed the influence of the number of train samples on the ML-based ROM performance. Considering N train = 47 instead of N train = 146, the average prediction error is 1.92× lower while the maximum prediction error is 1.84× lower. The average and maximum prediction error slightly decrease when considering N train = 346 instead of N train = 146, respectively 1.14× and 1.04×. It can be observed in Figure 9 that further expanding the training dataset, the prediction error converges around the value reached when N train = 346 is used.

Discussion
In this work, a fast and accurate method for braided stent deployment analysis has been proposed. It consists of a twostep workflow where FE simulations are used, first, to train a ML model that classifies successful and unsuccessful simulations and, next, to train a ML-based ROM that approximates the stent deployed configuration within the considered vessel. This approach was validated by studying the effect of a combination of geometrical and surgical parameters on the outcome of the stent deployment simulation: to this end, we employed a parametric idealised model of an intracranial artery characterised by a saccular aneurysm. The presented model relies on previous work to develop an optimised framework for numerical simulations of frictional Bézier 0.164 ± 0.120 26 (52%) 0.233 ± 0.147 32 (64%) 4 (8%) N cl = 3 0.029 ± 0.020 0 (0%) 0.082 ± 0.058 5 (10%) 0 (0%) N cl = 5 0.026 ± 0.017 0 (0%) 0.080 ± 0.059 5 (10%) 0 (0%) N cl = 8 0.024 ± 0.017 0 (0%) 0.077 ± 0.058 6 (12%) 0 (0%) contact interactions between wire-like structures discretised using beam elements and rigid surfaces [Aguirre and Avril. (2020); Bisighini et al. (2022)]. Based on this, here we proposed an efficient approach to simulate braided stent deployment, which allowed us to reduce the computational time required to collect the simulations needed for ML training. As explained in Section 2.1, stent positioning is performed by connecting a subset of the stent nodes to a virtual centerline built along its main axis and imposing a pre-computed displacement field on this centerline. These kinematic constraints lead the stent to follow the centerline and bend. Thanks to this technique, the computational time to perform a FE simulation of stent deployment takes, on average, 15 min. This makes it possible to build even large datasets in an acceptable amount of time: e.g., the largest dataset N s = 900 used in this analysis could be created in 12 h on 20 nodes of a cluster with 2.6 Ghz Intel Xeon Gold 6,132 CPUs and 6 Gb of RAM for each node. We trained six ML models on a dataset of varying size and obtained classifiers that were 80%-91% accurate in predicting the deployment outcome even with a relatively small dataset (N train = 50). Only increasing its size to N train = 500 and using SVM and NN models, we were able to perform binary classification with all the validation metrics between 89% and 97%. The surgical needs addressed by this model require us to favour the presence of FNs over that of FPs, i.e., we prefer to mislabel a successful simulation as a failure rather than the opposite. To minimize the presence of FPs then, high-specificity models are preferred. Therefore, for N train = 50, we would select the SVM model and for N train = 200, the ANN model. As shown in Figure 10, FNs misclassification may be explained by the fact that, in these cases, the stent is in a "boundary" situation in which one of its extremities lands immediately after the aneurysm neck. Introducing deformable walls, these cases would most probably fall into the "failure" condition. A similar situation is observed also for some FPs cases, i.e., one of the extremities of the stent lands within the aneurysm immediately before the aneurysm neck, but no such clear explanation was found for all the FPs.
For the regression, the POD algorithm was employed to extract the RBs while a GPR model to establish a mapping between the simulation parameter values and the projection coefficients. The results showed that GPR is strongly influenced by the input parameters: maintaining the same output data and changing the predictors from μ B to μ cl , the average prediction error could be decreased by more than 5×. This reduction can be explained by the less non-linear relationship present between these alternative predictors and the output variables. We carried out a sensitivity analysis on the number of RBs to be considered for the dimensionality reduction of the problem and on the training

FIGURE 8
Two examples from the test dataset: focus on 100 points sampled within the multivariate normal distribution predicted by the GPR for the displacement of the first node of the stent mesh. (A) HF (black) and predicted (green) solution using μ B as predictors. (B) HF (black) and predicted (red) solution using μ cl as predictors.
dataset size. The results showed that the prediction error stabilises with a relatively low number of RBs (15). On the other hand, with 47 training samples, we were able to achieve a maximum prediction error slightly lower than the spatial resolution of 3DRA; with 147 training samples, we were able to reduce it to 0.07 mm, half the spatial resolution of 3DRA. In our analysis, the average prediction error is never greater than 0.15 mm (3DRA spatial resolution) while the maximum prediction error is never greater than 0.4 mm (CTA spatial resolution) and is lower than 0.15 mm (3DRA spatial resolution) in 90% of the total test cases. By looking at the test examples reported in Figure 11, it can be observed that the prediction error is evenly distributed and the stent configuration accommodates very well the vessel curvature. In light of these results, we can conclude that a prediction error in the range of 0.01 mm-0.05 mm leads to an approximated stent that is very close to the HF results and such an error is achievable even with a relatively small dataset (47 simulations). In the cases characterised by larger prediction errors ( Figures 11C, E), the greatest differences can be observed at the stent extremities. These cases correspond to the same "boundary" situation that also affects the classification models, where one of the stent extremities lands immediately after the aneurysm neck. As mentioned above, by introducing deformable walls, the stent would most likely slip within the aneurysmal space in this situation, making the deployment unsuccessful. Thanks to the complete decoupling between the offline stage (extraction of RBs and training of the regression model) and the online one (prediction of new solutions), the computational cost is decreased from ∼15 min, using the FE simulation scheme proposed, to a few milliseconds.
The results here presented are promising as they demonstrate the ability of ML and reduced order modelling techniques to account for the non-linearities of the stent deployment problem and accurately model its outcome. Compared to other fast virtual stenting methods proposed in the literature, reduced order modelling has the  potential of reducing the computational cost of stent deployment simulation without simplifying the underlying biomechanical model [Larrabide et al. (2012); Spranger et al. (2015); Zhong et al. (2016)]. However, the current approach presents some limitations that need to be addressed. The first limitation is the ability to consider only one stent geometry at a time: due to the analytical formulation we used to build the braided stent mesh (Eq. 1), variations in the stent radius result in a different number of nodes and, therefore, number of DOFs. This in turn would render displacement vectors of different lengths that could not be assembled in the same snapshot matrix S. Considering stents of various sizes would require the construction of separate ROMs. Nevertheless, flow diverters are available in a finite and limited number of sizes, so creating different models could be feasible. An alternative solution could be the introduction of a sampling step to transform the stent mesh into a point cloud with a fixed number of points. Besides, the independent single-output GPR method used in this project does not model the correlation between outputs (Appendix 1): a multi-output GPR approach should be considered in the future (Liu et al., 2018). Secondly, both the stent and vessel FE models here considered are based on some simplifications: In particular, the penalty-based constraints at the wires interconnections for the stent and the rigid-wall hypothesis for the vessels. Given the high braiding angle and number of wires of flow diverters (Kim et al., 2008), the simplified contact technique employed for modelling wire-to-wire interactions reduces its radial and axial flexibility limiting the application within tortuous vessels (Zaccaria et al., 2020) and prevents modelling clinical practice as the "push-pull" technique (Ma et al., 2014). Therefore, with the aim of applying the same workflow to patient-specific geometries, the current approach Frontiers in Physiology 13 frontiersin.org

FIGURE 11
Six examples from the test dataset. The nodal absolute error E p between HF and predicted solution is shown as a colourmap. As comparative scale, for each solution we reported the diameter of the vessel and the aneurysm. The wire thickness is magnified (2×) to better visualize the stent. should be replaced by a more general contact formulation (Ma et al., 2012;Shiozaki et al., 2020). On the other hand, it is still unclear if the introduction of wall deformability would have a strong impact on the stenting simulation. Since intracranial arteries are more rigid than extracranial ones [Hayashi et al. (1980)], most studies assume rigid walls finding good agreement between the results of numerical simulations and experimental tests [Ma et al. (2012); Ma et al. (2013); Ma et al. (2014)]. However, some preliminary comparisons on idealised geometries highlighted the presence of differences in the stent deployed configuration, in particular in the radial direction (Fu and Xia, 2017;Cai et al., 2020). Moreover, deformable walls would enable modelling vessel straightening following the insertion of endovascular guides and stent sheaths (King et al., 2012;Gindre et al., 2017). Concerning the ROM construction, both of these improvements to the biomechanical model would increase the computational cost of the HF simulations, introduce a higher level of complexity in the modelled problem and, thus, require a larger training dataset. Moreover, deformable walls would most likely change the outcome of the deployment simulation, in particular in the "boundary" situations discussed above. Finally, the vessel geometry considered for this work is rather far from patient-specific geometries. Introducing more parameters will result in the need for a larger number of bases (hence, a larger training dataset) to predict the solution well. Therefore, efficient algorithms to construct the training dataset should be explored to reduce its size as much as possible as well as alternative ML algorithms for regression, e.g., neural networks.

Conclusion
This work represents the first attempt to combine finite element simulations with machine learning and reduced order modelling for the analysis of braided stent deployment. Its feasibility was demonstrated using an idealised vessel model, where a set of geometrical features can be controlled. Surgical decisions were also taken into account in the creation of the high-fidelity dataset. The two-step workflow allows the classification of deployment conditions with up to 95% accuracy and real-time prediction of the stent deployed configuration with a maximum prediction error always lower than the spatial resolution of computed tomography angiography (0.4 mm) and lower than that of 3D rotational angiography (0.15 mm) in 90% of test cases. Despite the simplified vessel shape and the assumption of rigid walls, this study is representative of the clinical scenario and can be extended to more realistic applications without modification. Current efforts are focused on understanding how many parameters are needed to fully describe patient-specific models. To represent such geometries in the reduced-order model, a statistical shape model will be used instead of the parametrization presented here. Future developments of the presented framework will also include the improvement of the clinical criterion used to assess the effectiveness of treatment: for example, the porosity of the stent at the neck of the aneurysm should be taken into account as it affects the outcome of flow diversion. In the future, a similar computational tool could be used by practitioners before and during intracranial aneurysm surgery to rule out conditions that would lead to unsuccessful deployment, visualize the stent configuration depending on the deployment site and check whether the chosen device covers one or more side branches.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.