Full text
A low dimensional surrogate model for a fast estimation of strain in the thrombus during a thrombectomy procedure Sara Bridio a , * , Giulia Luraghi a , Francesco Migliavacca a , Sanjay Pant b , Alberto GarcíaGonz ´alez c , Jose F. Rodriguez Matas a a Laboratory of Biological Structure Mechanics (LaBS), Department of Chemistry, Materials and Chemical Engineering “Giulio Natta”, Politecnico di Milano, Milan, Italy b Faculty of Science and Engineering, Swansea University, Swansea, Wales, UK c Laboratori de C`alcul Num `eric (LaC `aN), E.T.S. de Ingeniería de Caminos, Universitat Polit `ecnica de Catalunya - BarcelonaTech, Barcelona, Spain Keywords: Acute ischemic stroke Thrombectomy, Surrogate modeling Principal components analysis Kriging, Finite element method ABSTRACT Background: Intra-arterial thrombectomy is the main treatment for acute ischemic stroke due to large vessel occlusions and can consist in mechanically removing the thrombus with a stent-retriever. A cause of failure of the procedure is the fragmentation of the thrombus and formation of micro-emboli, difficult to remove. This work proposes a methodology for the creation of a low-dimensional surrogate model of the mechanical thrombectomy procedure, trained on realizations from high-fidelity simulations, able to estimate the evolution of the maximum first principal strain in the thrombus. Method: A parametric finite-element model was created, composed of a tapered vessel, a thrombus, a stentretriever and a catheter. A design of experiments was conducted to sample 100 combinations of the model pa-rameters and the corresponding thrombectomy simulations were run and post-processed to extract the maximum first principal strain in the thrombus during the procedure. Then, a surrogate model was built with a combination of principal component analysis and Kriging. Results: The surrogate model was chosen after a sensitivity analysis on the number of principal components and was tested with 10 additional cases. The model provided predictions of the strain curves with correlation above 0.9 and a maximum error of 28%, with an error below 20% in 60% of the test cases. Conclusions: The surrogate model provides nearly instantaneous estimates and constitutes a valuable tool for evaluating the risk of thrombus rupture during preoperative planning for the treatment of acute ischemic stroke. when a blood clot obstructs a cerebral artery and prevents the perfusion of downstream cerebral tissues. This causes the development of an infarct zone in cerebral tissues, which is irreversible if not re-perfused in a short time (Saver, 2006). The main treatment for AIS due to large vessel occlusions is intra-arterial thrombectomy (IAT) (Goyal et al., 2016), an endovascular technique which aims at mechanically removing the thrombus from the patient’s artery. The IAT procedure can be performed with stent-retrievers, or aspiration catheters, or in combined techniques involving both devices (Kang and Park, 2017). In the case of stent-retriever IAT, the crimped stent is minimally invasively navigated to the occlusion site through a micro-catheter, aided by angiographies to correctly position the device. Once the occlusion is reached, the catheter is unsheathed, and the stent is deployed. Stent-retrievers are made of super-elastic Nickel–Titanium (NiTi) alloys, hence in this phase the stent tends to recover its original non-crimped configuration and can entrap the blood clot. The stent-thrombus complex is then retrieved out of the patient to restore cerebral artery perfusion. A factor determining the interaction with the stent-retriever is thrombus composition. Thrombi with higher red blood cell content are softer and usually well integrated in the stent struts, while thrombi with higher fibrin and platelet content are stiffer and tend to remain between the stent and the vessel walls, preventing the recanalization (Ospel et al., 2021). Despite being the standard of care for large-vessel occlusion AIS, the IAT procedure still requires optimization to improve the clinical *Corresponding author. Computational Biomechanics Laboratory – LaBS, Department of Chemistry, Materials and Chemical Engineering “Giulio Natta”, Politecnico di Milano, Piazza L. Da Vinci 32, 20133, Milano, Italy. E-mail address: [email protected] (S. Bridio). 1. Introduction Acute Ischemic Stroke (AIS) is a neurovascular pathology occurring
2 established methodology used to lower problem complexity. Lataniotis et al. (2020) combined Gaussian Process modeling (Kriging) and polynomial chaos expansion surrogates and kernel Principal Component Analysis (kPCA) to perform uncertainty quantification in two engineering problems. This methodology has been also applied for the uncertainty quantification and Bayesian calibration of a hydrological model (Nagel et al., 2020). Rocas et al. (2021) used a combined kPCA with metamodeling to carry out uncertainty quantification analysis in crashworthiness. Li et al. (2020) integrate dimension reduction (Principal Component Analysis, PCA) together with Kriging surrogate models to perform sensitivity analysis in models with high-dimensional output. In this work, a method combining dimensionality reduction and surrogate modeling techniques is implemented and applied for a fast and accurate estimation of the outcomes in thrombectomy procedures, avoiding costly FEM simulations. As a first feasibility study, a low dimensional surrogate model of high-fidelity virtual simulations of the IAT procedures is developed, able to predict the evolution of the maximum first principal strain in the thrombus due to the interaction with the stent-retriever, following the scheme in Fig. 1. 2. Methods 2.1. Parametric FEM model of the thrombectomy procedure Fig. 2 shows the FEM model created for the high-fidelity IAT procedure simulations. To keep a limited number of parameters, a simplified tapered vessel geometry was chosen. The part of the vessel with large diameter corresponds to the internal carotid artery (ICA) segment, while the part with small diameter corresponds to the distal tract of the middle cerebral artery (MCA). The vessel, of total length of 200 mm, was discretized with rigid quadrilateral shell elements (average element size 0.35 mm). A blood clot was positioned in the smaller straight part of the vessel, as thrombus occlusions mostly occur in the distal segments of the MCA (Dutra et al., 2019), with a diameter equal to 90% of the vessel diameter. The clot was discretized with an average of 12400 linear tetrahedral elements (average element size 0.2 mm) and modeled with a quasi-hyperelastic foam material proposed by (Kolling et al., 2007). The model approximates the principal components of the Kirchhoff stress, τ i, i=I,II,III, as: τ i=f(λi) − f(J− ν 1−2 ν )(1) where ν is the initial Poisson’s ratio (0.3 from (Luraghi et al., 2021a)), and f(•) is a function identified from a uniaxial stress-strain curve, τ = g(λ): f(λ) = λg(λ) + λ− ν g(λ− ν ) + ⋯+λ(− ν )ng(λ(− ν )n).(2) The values of f(λ)are calculated at the beginning of the analysis and do not require an analytical expression for function f(λ)since it is calculated from the uniaxial stress-strain data. This material formulation is available within the finite-element solver LS-DYNA (ANSYS, Canonsburg, PA, USA) used to perform the simulations. Unconfined compression tests conducted on ex vivo thrombi retrieved from stroke patients (treated at the Erasmus University Medical Center in Rotterdam, The Netherlands) were used to calibrate the model (Luraghi et al., 2021b). The samples were subjected to 80% compression of their initial height, at a strain rate of 10% per second (full details about the testing protocol can be found in (Luraghi et al., 2021b)). Clots of different compositions were tested, ranging from fibrin-rich to red blood cell-rich clots (Fig. 3). Linear interpolation was used to approximate the stress-strain behavior of clots of compositions between the tested ones. The chosen stent-retriever was a Trevo ProVue (Stryker, Kalamazoo, MI, USA) with 4 mm diameter and 40 mm length. The geometry of the stent was reconstructed by identifying the repeating cell and creating a outcome of the patients (Ospel et al., 2021), in particular to reduce the recanalization time and the vascular damage (Kühn et al., 2020). A frequent complication of the IAT procedure is the thrombus embolization, caused by the fragmentation of the blood clot and the occlusion of smaller distal vessels, more difficult to recanalize (Kaesmacher et al., 2017; Georgakopoulou et al., 2021). The improvement of clinical outcome could be achieved by providing tools for performing an accurate and fast biomechanical analysis as pre-operative planning (which is now based only on clinical imaging), to identify the most suitable procedure and device, considering the vascular geometry, the occlusion site and the blood clot characteristics of the specific patient. Besides, the procedure itself and the design of the devices are subjected to improvement in maximizing the grip of the stent to the clot and minimizing the clot fragmentation and vascular trauma. Recent studies in the literature proposed models for high-fidelity simulations of the IAT procedure with stent-retriever. In (Luraghi et al., 2021a), a finite-element model (FEM) of the thrombectomy procedure was created, composed of rigid-walls vessels, a quasi-hyperelastic foam model for the thrombus and a stent-retriever model discretized with beam-elements. The model was successfully validated with in vitro experiments replicating the clinical procedure. The same authors used the model developed in (Luraghi et al., 2021a) to reproduce a patient-specific case (Luraghi et al., 2021b), where the thrombus model, including fracture properties, developed by Fereidoonnezhad et al. (2021) was integrated, and to perform a study on the impact of patient-specific vascular geometry on the outcome of the IAT procedure (Bridio et al., 2021). Nevertheless, in clinical applications the use of high-fidelity numerical simulations of IAT is unaffordable, or even impossible, due to both the required time and computational cost. In (Luraghi et al., 2021b), the IAT simulation in a patient-specific vessel geometry is reported to run in approximately 24 h (on 40 CPUs of an Intel Xeon64 with 256 GB of RAM). The required time is incompatible with clinical pre-operative planning, where the treatment to AIS should be provided during the initial hours after symptoms onset (Saver, 2006). From another perspective, FEM simulations are also impractical when a big volume of virtual procedures needs to be analyzed, as it happens in the case of an in silico trial, whose possibility of substituting or integrating clinical trials has become in recent years a promising field of research (Viceconti et al., 2016). A possible solution to this limitation is the creation of surrogate models of high-fidelity simulations of the IAT procedure. Surrogate modeling consists of finding a mapping between a set of input parameters describing a model and an output function representing the response. Once the mapping function is found, the surrogate model can nearly instantaneously estimate the output, given a new set of input parameters. In the literature, few studies with biomedical application used the surrogate modeling technique to speed up the obtainment of results with respect to the use of high-fidelity simulations, with neural networks (NN) being the main approach. In (Liang et al., 2020), deep NN were used to estimate pressure and flow velocity distributions inside the thoracic aorta, trained on computational fluid dynamics simulations. The same authors used FEM simulations to train deep NN to predict stress distributions in the aortic walls (Liang et al., 2018). In (Madani et al., 2019), NN were trained with FEM simulations to predict the maximum Von Mises stress in the walls of atherosclerotic arteries, which is used to estimate the risk of plaque rupture. Finally, in (Pellicer-Valero et al., 2020) a combination of PCA and Feedforward NN, based on FEM simulations, was used to provide real-time inference on the outcome of liver surgery, with the final aim of having a computer aided surgery. These models demonstrate the great potential of the use of surrogate models trained on FEM simulations and able to substitute them when a real-time response is needed. All the above NN-based models required several hundreds or even thousands of simulations for the training phase. In different disciplines and context of applications, the combination of dimensionality reduction with surrogate modeling is a well-
3 parametric model in Matlab (Math-Works, Natick, MA, USA) (Fig. 4A). The cross-section of the stent struts was identified with the use of a confocal laser microscope. The reconstructed stent geometry was then discretized with 1032 linear beam elements (average element length 0.2 mm). The Hughes-Liu formulation with cross-section integration was chosen. The NiTi material of the stent-retriever was modeled in LSDYNA with the available shape-memory material formulation (Livermore Software Technology Corp, 2019). The material parameters were calibrated by numerically reproducing an experimental uniaxial tensile test performed on a Trevo ProVue stent. The material parameters used for the simulations are listed in Fig. 4. The developed FEM model was parametrized with 5 relevant features (green text in Fig. 2A): 1) D_ICA, the maximum diameter of the tapered vessel, 2) D_MCA, the minimum diameter of the tapered vessel, 3) L_clot, the length of the blood clot, 4) %FIB, the percentage of fibrin content of the blood clot (the composition of the thrombus determines different mechanical properties (Duffy et al., 2019)), and 5) X_stent, the distance between the proximal end of the thrombus and the head of the stent, giving an indication of the relative position of the stent-retriever with respect to the clot. The last parameter is a derived parameter: it is obtained with the formula R⋅(L_stent – L_clot), where L_stent is the length of the stent-retriever and R is a parameter which can take a value between 0 and 1, ensuring that the stent always covers the entire length of the clot. The 5 parameters vary in ranges as listed in Table 1. The parameters referring to the vessel and the thrombus were determined on the basis of literature values (Boodt et al., 2020; Rai et al., 2013). Fig. 2B shows examples of models with different parameters. Fig. 1. Schematic representation of the steps for the construction and testing of the surrogate model for the prediction of the maximum first principal strain (MaxFP strain) in the thrombus during a thrombectomy procedure. Fig. 2. A) FEM model for the IAT procedure and the 5 model parameters: small and large diameters of the tapered vessel, D_MCA and D_ICA, respectively (MCA: middle cerebral artery; ICA: internal carotid artery); length (L_clot) and composition in terms of fibrin content (%FIB) of the thrombus; positioning of the stent-retriever with respect to the thrombus (X_stent). B) Examples of models with different parameters. Fig. 3. Stress-strain curves obtained from unconfined compression tests on ex vivo thrombi with different fibrin (FIB) and red blood cells (RBC) content.
4 2.2. Design of experiment (DoE) and creation of the training dataset A DoE was conducted to create the training dataset for the surrogate model, which is constituted by a set of FEM simulations of IAT procedure. The parameter space was sampled by 100 points from the Sobol sequence, which produces a quasi-random space-filling DoE while adhering to uniformity of samples’ distribution (Garud et al., 2017). The DoE resulted in the creation of 100 different FEM models of IAT. The geometrical features (D_ICA, D_MCA, L_clot, X_stent) were varied thanks to the parametrization of the geometry of the FEM model. The different clot compositions (%FIB) were implemented by adjusting the material parameters for the blood clot model, based on experimental tests, as detailed in Subsection 2.1. Once prepared the FEM models, 100 high-fidelity simulations of the IAT procedure were run with the FE-solver LS-DYNA. Each simulation is composed of 3 steps (Fig. 5): 1) catheter and stent tracking: the guide catheter (0.5 mm diameter, modeled as rigid) is positioned inside the vessel, pushing the blood clot against the vessel walls. A friction-less soft penalty contact is defined between the clot and the catheter. Then, the stent-retriever is crimped inside the catheter, by imposing the movement of the tip of the stent along the centerline of the catheter. At the end of the crimping the stent is positioned at the selected location with respect to the occlusion. A hard penalty contact is defined between stent and catheter; 2) stent deployment: the deployment of the stent-retriever is performed simulating the unsheathing the catheter with the progressive removal of the contacts between the stent and progressive portions of the catheter; 3) stent and thrombus retrieval: the retrieval of the stent-thrombus complex is simulated by imposing the movement of the stent tip along the vessel centerline. A rough soft penalty contact is defined between thrombus and vessel wall, with friction coefficient of 0.1. Between the stent and the thrombus, a soft penalty contact is defined, with friction coefficient of 0.2, while a hard penalty contact is defined between stent and vessel wall. A mass proportional damping of 10 s −1 was applied to the thrombus to achieve numerical stability without an excessively reduced time-step (Luraghi et al., 2018). The setting of the simulations was derived from the work in (Luraghi et al., 2021a), where the thrombectomy simulation was also validated with in vitro experiments. The results of the FEM simulations were post-processed to extract the Fig. 4. A) Model of the Trevo ProVue stent-retriever. B) Stress-strain curve for the NiTi material. C) Parameters used for modelling the stent material in LS-DYNA (Livermore Software Technology Corp, 2019). Table 1 Ranges of the 5 parameters of the FEM model of the IAT procedure. Parameter Range D_ICA [4–7.5] mm D_MCA [2.5–4] mm L_clot [5–35] mm %FIB [5–95] % X_stent [0–35] mm Fig. 5. Steps of the IAT procedure (the catheter is removed from visualization to show the stent-retriever): 1) the crimped stent-retriever is positioned through the catheter at the occlusion location; 2) the catheter is unsheathed and the stent is deployed to entrap the thrombus; 3) the stent is retrieved to remove the thrombus from the vessel.
5 X=UΣVT(3) where the matrix U∈Rd×d contains the orthonormal eigenvectors of the matrix XXT, the matrix V∈Rns×ns contains the orthonormal eigenvectors of the matrix XTX, and the matrix Σ∈Rd×ns has non-zero values in the diagonal only. The values on the diagonal of Σ represent the singular values of the matrix X (Golub and Reinsch, 1971). The singular values in Σ appear in descending order: the first singular value is the one explaining the highest amount of variance. The dimensionality of the system was then reduced to the first k singular values, depending on the percentage of the system information that needed to be explained. The results obtained with 4 different percentages of explained variances were compared, namely 60%, 70%, 80% and 90%, which in the following will be referred to as Model 1, Model 2, Model 3 and Model 4, respectively. The choice of the percentage of explained variance determines the number k of singular values for the dimensionality reduction: the initial data contained in X was projected onto the reduced space, by performing the operation: A= UTX(4) where the matrix U∈Rd×k contains the first k orthonormal eigenvectors of U, and the matrix A∈Rk×ns contains the components of the projections of each MaxFP strain curve onto the new reduced space. Kriging models. After performing the reduction of the dimensionality of the system, k independent Gaussian Process surrogate models (also known as Kriging models) were trained to identify the mapping function between the set of 5 parameters of the FEM model and each of the components contained in A∈Rk×ns for the respective strain curve. For a detailed explanation of Kriging models and background the reader is referred to (Keane et al., 2008), whereas the main concepts are here reported. A Kriging model identifies a response surface by interpolating between points constituting realizations of a system. It assumes that the observed output is the result of a stochastic process with a Gaussian distribution. The mapping between the output Y and a vector of input data z can be expressed as: Y(z) = β+G(z)(5) where β is an unknown hyperparameter and G(z)is a Gaussian process with zero mean and covariance given by: Cov(G(z,z’)) = σ 2 zk(z,z′)(6) where σ 2 z is the variance of the process and k(z,z′)is a correlation (kernel) function. In this case, a Gaussian kernel function was chosen: k(z,z’) = exp[−1 2∑ p m=1(zm−z′ m)2 θ2 m](7) where p is the number of parameters of each input vector and θ= {θ1, Fig. 6. Example of MaxFP strain curve and contour plot of the maximum first principal strain at the end of stent deployment: the red areas, in contact with the stent struts, are the ones at risk of fracture initiation. output of interest. Due to the simplified geometry of the vessel, in only 2 simulations the thrombus was not removed from the vessel (due to a loss of engagement between stent and thrombus, while in the other 98 simulations the procedure was successful. The chosen output of interest was the evolution of maximum First Principal (MaxFP strain in the thrombus during the procedure (Fig. 6. More precisely, at each time step of the simulation, the 10 clot elements with the highest MaxFP strain were identified and averaged: the output extracted from each simulation is a curve reporting the evolution of the value during the procedure. The choice of the MaxFP strain was based on the finding in (Tutwiler et al., 2020, where the fracture of thrombi is demonstrated as strain-driven. Each simulation is made of 217 time steps, leading to an output dataset made of 217 × 100 MaxFP strain values. 2.3. Creation of the surrogate model A surrogate model was built for the prediction of the MaxFP strain curves. The creation of the surrogate model was entirely performed in Matlab and was made of two main steps. Principal Component Analysis (PCA). The first step consisted in the reduction of the dimensionality of the output dataset obtain from the FEM simulations. PCA (García-Gonz´alez et al., 2020 is a technique that reduces the dimensionality of a system belonging to a space Rd×ns (where ns is the number of samples of dimension d, by projecting it into a reduced space of dimension Rk×ns , with k≪d. The dimension k must be chosen such that the information contained in the original data is properly kept. The principal components are orthonormal vectors defining a new basis for describing the data. An eigenvalue or singular value is associated to each principal component, whose magnitude, with respect to the others, indicates the amount of variance (i.e. information of the system explained by that principal component. Therefore, the analysis of the eigenvalues or singular values allows to determine the optimal number of dimensions k of the reduced space. In this case, the system to reduce is the matrix X ∈ Rd×ns containing the MaxFP strain values, where ns = 100 is the number of thrombectomy simulations, and d = 217 is the number of time steps in each simulation. After centering the data, a PCA was carried out by performing a Singular Value Decomposition (SVD of the matrix X. SVD is a technique which allows to avoid matrix diagonalizations required by standard PCA, by providing a factorization of the X matrix of the form:
6 Error =‖x−x‖ ‖x‖×100 (8) where the vector x∈Rd contains the MaxFP strain calculated by the FEM simulation during the IAT procedure, and the vector x∈Rd contains the MaxFP strain predicted by the surrogate model. The goodness of the predictions of the MaxFP strain curves were also assessed by calculating the Pearson correlation coefficient (Cohen et al., 2009) between the predicted curve and the one obtained from the FEM simulation. In addition, to better characterize the ability of the surrogate model in describing the problem, an analysis of variance (ANOVA) (Gelman et al., 2005) was performed to assess the influence of the input parameters and their interactions on the outcome of both the high-fidelity and the surrogate model. For each parameter or pair of parameters, the ANOVA analysis provides a p-value indicating its significance. A p-value<0.05 indicates that the parameter has significative influence on the determination of the output. The ANOVA analysis was performed on the 100 samples used for building the surrogate, with a leave-one-out approach. The interest in this case is not on the accuracy of the prediction, but to see the ability of the surrogate model to describe the influence of the input variables on the output. For this analysis, the 5 input parameters constitute the groups, while the chosen output is the mean MaxFP strain value for each realization. Thus, two ANOVA analyses were performed to determine the influence of each of the 5 input parameters and of their pairwise interactions on the mean MaxFP strain values, obtained with the simulations and with the surrogate model. The results of the analyses are compared to assess if the relevance of the input variables in the FEM simulations is respected in the surrogate model. 3. Results 3.1. Dimensionality reduction The first step of the creation of the surrogate model consisted in the reduction of the dimensionality of the problem by means of a PCA. In Fig. 7 a graph of the individual and cumulative explained variance of the eigenvalues is shown. The individual explained variance is the amount of the system information explained by each eigenvalue, while the cumulative explained variance indicates the amount of the system information explained by the sum of a progressive number of eigenvalues (which sums to 1 when all the eigenvalues are considered). Four percentages of system information were considered: to keep 60% of the information, 3 principal components need to be used (Model 1); to keep 70%, 80% and 90% of the information, 6, 12 and 25 principal components are respectively needed (Models 2, 3 and 4, respectively). 3.2. Surrogate model predictions Four surrogate models were then built to predict the MaxFP strain evolution during the thrombectomy procedure, one for each selected percentage of system information. To test the predictive ability of each model, 10 additional simulations of the IAT procedure were run, with random combinations of the 5 input parameters of the FEM model. Fig. 8 shows the corresponding parameters together with the result of the FEM simulation (dashed blue line). The prediction error in each test case was calculated with Eq. (8), obtaining the results shown in Fig. 9, top left. The increase in the number of principal components does not produce significative changes in the prediction error. In the bottom part of Fig. 9, the MaxFP strain curves predicted by the surrogate models are compared to the true curve from the FEM simulation of test cases 5 and 8 (the curves obtained with Model 4 are omitted as they were superimposed to the ones obtained with Model 3). The two examples show that the predicted curves are very similar with all the tested number of principal components, with differences only in the initial tract. The choice of Model 2 (surrogate model built with 6 principal components and 70% of the system information) appears a reasonable trade-off between the Fig. 7. Individual and cumulative explained variance of the eigenvalues. The first 3 principal components (PCs) explain 60% of the system information (Model 1), the first 6 PCs the 70% (Model 2), the first 12 PCs the 80% (Model 3) and the first 25 PCs the 90% (Model 4). θ2, …, θp} is a vector of p unknown hyperparameters. In the current application, the vector of input data z contains, for each sample, the 5 parameters of the FEM model (D_ICA, D_MCA, L_clot, %FIB, X_stent, and the observed outputs Y(z) are the corresponding components of the projections of the MaxFP strain curve in the reduced space, contained in the columns of matrix A ∈ Rk×ns . A Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS algorithm (Liu and Nocedal, 1989, a quasi-Newtonian optimization algorithm, was selected to find the optimal values of the hyperparameters (β, σ z, θ) of the Kriging model. This procedure was repeated k times to obtain k independent Kriging models: each of them is able to estimate one of the singular components ai where i = {1, 2, …, k}, given a new set of input parameters (D_ICA, D_MCA, L_clot, %FIB, X_stent. The vector a ∈ Rk of estimated singular components, multiplied by the matrix U ∈ Rd×k, provides the vector x ∈ Rd, which is the estimation of the MaxFP strain curve for the new set of input parameters. 2.4. Surrogate model evaluation The strain curves predicted by the surrogate model were compared with the curves extracted from the output of the simulations by means of the normalized Euclidean distance between corresponding curves with the formula:
7 reduction to a low number of principal components and the amount of preserved information of the system. Fig. 9, top right, shows the values of the Pearson correlation coefficient between the predicted strain curve and the one obtained from the FEM simulations, using Model 2: for all the test cases the coefficients are very close to 1, demonstrating a strong correlation between the true and predicted curves. Fig. 8 shows the comparison of the MaxFP strain curves calculated with the FEM simulations (dashed blue line) and predicted by Model 2 (red line) for all the 10 test cases, showing the ability of the surrogate model to predict the shape of the strain curves, but also demonstrating the smoothing effect of the model. In terms of computational time, the surrogate model provides nearly instantaneous predictions of the strain curves, while each high-fidelity FEM simulation of the IAT procedure required on average 40 h on 20 CPUs of an Intel Xeon64 with 120 GB of RAM memory. 3.3. Analysis of the influence of the input parameters An ANOVA analysis was performed to assess if the chosen surrogate model (Model 2) can replicate the behavior of the high-fidelity model in terms of influence of the input parameters (as described in Section 2.4). The results of the ANOVA analyses for the FEM model and the surrogate model are summarized in Table 2, which reports the input parameters, or pairs of parameters, found as significative for the determination of the mean value of the MaxFP strain in the two models. The significant parameters for the FEM output were found to be significant also for the surrogate model, for which two additional parameters, namely X_stent and D_ICA, were also found as relevant. Fig. 8. Surrogate model predictions (Model 2, using 6 principal components) of the MaxFP strain curves (red lines) for the 10 test cases and comparison with the MaxFP strain curves obtained with the FEM simulations (dashed blue line). Each panel shows the values of the 5 input variables used to parametrize the FEM model described in Section 2.1.
8 4. Discussion In the present work, a method for creating a low dimensional surrogate model of high-fidelity simulations of the IAT procedure was proposed. The surrogate model can provide an estimation of the evolution of the maximum first principal strain in the thrombus during a thrombectomy procedure performed on a simplified tapered vessel geometry, given as input parameters the minimum and maximum diameters of the vessel, the length of the thrombus, its composition and the relative position between the thrombus and the stent-retriever. The presented model constitutes the first attempt of using the surrogate modeling techniques for the estimation of the outcome of an IAT procedure. The creation of the model was based on a combination of PCA and Kriging interpolation rather than on a NN approach to limit the number of samples required for the training of the model. The complexity of the interaction between the stent-retriever and the thrombus requires a fine discretization, both spatial and temporal, that implies a high computational cost for each simulation (approximately 40 h on a system with 20 CPUs and 120 GB of RAM). A sensitivity analysis was conducted on the number of principal components to be considered for the reduction of the dimensionality of the problem. The results (Fig. 9) showed that good predictions of the MaxFP strain curves can be obtained with a very limited number of principal components (up to 3, corresponding to Model 1, with 60% of the information of the original system). However, Model 2 (using 6 principal components and 70% of information) was chosen to present the final results, as a reasonable trade-off between the low dimension of the reduced space and the amount of preserved information of the system. With the chosen surrogate model, the shapes of the MaxFP strain curves obtained with the FEM simulations were well approximated for each test case (Fig. 8). In particular, the model was able to accurately capture the initial tract of each curve, up to the end of stent deployment. This is the phase with a sudden increase in the MaxFP strain in the thrombus, that may be decisive for the initiation of the fracture of the thrombus. When considering the point-to-point error (Fig. 9, top left), 60% of the cases had error below 20%, the minimum error was 11% (test 8), while the maximum error was 28% (test 5). Despite these errors suggest there remains scope for improvement, the model proved able to capture well the shapes of the strain curves. Indeed, the analysis of the Pearson correlation coefficient between the predicted curves and the ones obtained from the FEM simulations (Fig. 9, top right) demonstrated a strong correlation between the two curves, with values of the coefficient close to 1 for each test case. Results from the ANOVA analysis of Model 2 indicated that the surrogate model is able to reproduce the behavior of the high-fidelity model in terms of influence of the input parameters (Table 2). The analysis on the FEM outputs indicated four parameters or pair of parameters as significant: the interactions between L_clot and X_stent, between L_clot and D_MCA, between %FIB and X_stent, and L_clot. The results of this analysis seem reasonable, as the interaction of L_clot and X_stent determines the amount of contacts between the thrombus and the stent struts, which are areas of high strains (Fig. 6); D_MCA determines how much the stent is free to expand, and therefore the force exerted on the thrombus (in combination with L_clot that again is an Fig. 9. Top left: prediction errors of the MaxFP strain curve in the 10 test cases obtained with the surrogate models using 60% (Model 1), 70% (Model 2), 80% (Model 3) and 90% (Model 4) of the system information (using 3, 6, 12 and 25 principal components (PCs) respectively). Top right: Pearson correlation coefficients between the predicted curves and the curves from FEM, shown for the case of Model 2. Bottom: tests 5 and 8 are shown as examples of the obtained strain curves (the ones obtained with Model 4 are omitted as they were superimposed to the ones obtained with Model 3). Table 2 Results of the ANOVA analysis: input parameters, or pairs of parameters, found as relevant for the determination of the mean value of the MaxFP strain obtained with both the FEM simulation and the surrogate model, in decreasing order of influence (p =p-value). FEM simulation Surrogate model 1. L_clot – X_stent (p =0.0015) 1. L_clot – X_stent (p <0.001) 2. L_clot – D_MCA (p =0.0085) 2. %FIB – X_stent (p =0.0039) 3. %FIB – X_stent (p =0.0209) 3. L_clot (p =0.0098) 4. L_clot (p =0.026) 4. L_clot – D_MCA (p =0.0221) 5. X_stent (p =0.0312) 6. D_ICA (p =0.0412)
9 5. Limitations The obtained results are encouraging for the possibility of applying this surrogate modeling method for the estimation of strain in the thrombus during a thrombectomy procedure, but they indicate that there are some limitations that need to be addressed. A first limitation lies in the choice of the MaxFP strain as output of interest. The MaxFP strain is calculated, at each time step, as the average value of the 10 elements with highest strain values. This means that the 10 elements are not necessarily the same in each time step, determining a very irregular shape of the strain curves (see Fig. 6 as an example). The choice of a smoother variable would provide more regular curves, which would be more easily approximated by the surrogate model. Secondly, the model itself may be improved by increasing the number of samples in the training dataset or by building an ad hoc kernel in the Kriging model (Eq. (7)) which better adapts to the prediction of the principal components’ coefficients. Another limitation is the simplified vessel geometry. For exploring the feasibility of building this kind of surrogate model a straight tapered vessel was chosen. To extend the application of the modeling technique to patient-like cases, a more realistic vessel geometry can be considered. The parametric FEM model will include curved vessels to represent the carotid siphon and a bifurcation to represent the so-called T-junction (bifurcation of the ICA into MCA and anterior cerebral artery). Once the new training dataset is obtained, the technique for the creation of the surrogate model will not be different from the one implemented in this study. Additionally, the inclusion of the fibrin content of the thrombus as an input parameter for the prediction of the strain level is a source of uncertainty. At the moment, the resolution of the images acquired before the clinical intervention does not usually allow to determine the thrombus composition before it is retrieved. Therefore, large uncertainties would be associated to this parameter. The developed model can provide a probability of failure by interrogating the model with a range of different possible thrombus compositions. Finally, some limitations are associated to the FEM model used for the training of the surrogate model, also discussed in previously published papers (Luraghi et al., 2021a, 2021b, 2021c). Specifically, the inclusion of a compliant model of the vascular wall and of a cohesive model between the blood clot and the vessel walls (as proposed for example by Oyekole et al. (2021)) may be accounted for to obtain more accurate strain (or stress) curves for the thrombus. However, the methodology presented in this paper will apply the same, regardless of the fidelity of the FEM model. 6. Conclusion The developed methodology for the prediction of the maximum strain in the thrombus during a thrombectomy procedure (which can be adapted to the prediction of other variables of interest, e.g. the maximum stress) appears valuable for a better understanding of the thrombus mechanics during the treatment and to predict the thrombus evolution and potential fragmentation. The surrogate model provides nearly instantaneous predictions of the strain level, making it suitable for performing a pre-operative planning. By building a model for each of the most used commercial devices, when a new patient needs to be treated it would be possible to make few measurements on the patient’s vasculature and thrombus to interrogate the surrogate models and obtain an estimation of the level of strain produced in the clot by each device, choosing the one with the lowest risk of fragmenting the thrombus. A similar application could be used for future in silico trials. The creation of a surrogate model would allow to estimate the level of strain produced in the thrombus by a new device in a big number of virtual patients, without the need of running a computationally demanding and time-consuming high-fidelity simulation for each of them. Funding This project has received funding from the European Union’s Horizon 2020 research and innovation program under grant agreement No 777072 and from the MIUR FISR-FISR2019_03221 CECOMES. CRediT authorship contribution statement Sara Bridio: Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. Giulia Luraghi: Writing – review & editing, Validation, Software, Methodology, Formal analysis, Conceptualization. Francesco Migliavacca: Writing – review & editing, Supervision, Resources, Methodology, Funding acquisition, Conceptualization. Sanjay Pant: Writing – review & editing, Software, Methodology, Conceptualization. Alberto García-Gonz´ alez: Writing – review & editing, Writing – original draft, Validation, Supervision, Software, Project administration, Methodology, Formal analysis, Conceptualization. Jose F. Rodriguez Matas: Writing – review & editing, Writing – original draft, Validation, Supervision, Software, Resources, Project administration, Methodology, Formal analysis, Conceptualization. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability No data was used for the research described in the article. Acknowledgments The authors would like to thank Mattia Belloni for his help in the initial stage of this study. References Boodt, N., Compagne, K.C.J., Dutra, B.G., Samuels, N., Tolhuisen, M.L., Alves, H.C.B.R., Kappelhof, M., Lycklama ` A Nijeholt, G.J., Marquering, H.A., Majoie, C.B.L.M., Lingsma, H.F., Dippel, D.W.J., Van Der Lugt, A., 2020. Stroke etiology and thrombus computed tomography characteristics in patients with acute ischemic stroke: a mr clean registry substudy. Stroke 1727–1735. https://doi.org/10.1161/ STROKEAHA.119.027749. Bridio, S., Luraghi, G., Rodriguez Matas, J.F., Dubini, G., Giassi, G.G., Maggio, G., Kawamoto, J.N., Moerman, K.M., McGarry, P., Konduri, P.R., Arrarte Terreros, N., Marquering, H.A., van Bavel, E., Majoie, C.B.L.M., Migliavacca, F., 2021. Impact of the internal carotid artery morphology on in silico stent-retriever thrombectomy outcome. Front. Med. Technol. 3, 1–13. https://doi.org/10.3389/ fmedt.2021.719909. Cohen, I., Huang, Y., Chen, J., Benesty, J., 2009. Pearson correlation coefficient. In: S, S., B.V, B.M. (Eds.), Noise Reduct. Speech Process, pp. 1–4. https://doi.org/10.1007/ 978-3-642-00296-0_5. Duffy, S., Mccarthy, R., Farrell, M., Thomas, S., Brennan, P., Power, S., O’Hare, A., Morris, L., Rainsford, E., Maccarthy, E., Thornton, J., Gilvarry, M., 2019. Per-pass analysis of thrombus composition in patients with acute ischemic stroke undergoing indication of the number of contacts with the stent, while %FIB is related to the mechanical properties of the blood clot, determining the level of strain. The same four parameters were found as significant for the determination of the outputs of the surrogate model. In this case, two other parameters (X_stent and D_ICA) were found as significant, although with higher p-values with respect to the others. Therefore, besides the accuracy loss associated with the reduced nature of the model, the statistical analysis shows that the description of the behavior of the output provided by the surrogate model is in good agreement with high-fidelity FEM simulations.