Uncertainty quantification of a three-dimensional in-stent restenosis model with surrogate modelling

In-stent restenosis is a recurrence of coronary artery narrowing due to vascular injury caused by balloon dilation and stent placement. It may lead to the relapse of angina symptoms or to an acute coronary syndrome. An uncertainty quantification of a model for in-stent restenosis with four uncertain parameters (endothelium regeneration time, the threshold strain for smooth muscle cell bond breaking, blood flow velocity and the percentage of fenestration in the internal elastic lamina) is presented. Two quantities of interest were studied, namely the average cross-sectional area and the maximum relative area loss in a vessel. Owing to the high computational cost required for uncertainty quantification, a surrogate model, based on Gaussian process regression with proper orthogonal decomposition, was developed and subsequently used for model response evaluation in the uncertainty quantification. A detailed analysis of the uncertainty propagation is presented. Around 11% and 16% uncertainty is observed on the two quantities of interest, respectively, and the uncertainty estimates show that a higher fenestration mainly determines the uncertainty in the neointimal growth at the initial stage of the process. The uncertainties in blood flow velocity and endothelium regeneration time mainly determine the uncertainty in the quantities of interest at the later, clinically relevant stages of the restenosis process.


Introduction
Coronary heart disease is mainly due to the accumulation and development of atherosclerotic plaques, which narrow the vessel lumen and reduce the flow of blood. It can cause ischaemia or further evolve into a myocardial infarction. The most common treatment is percutaneous coronary intervention with stent deployment [1,2]. However, in addition to displacing the plaque from the lumen and restoring the blood flow, the balloon dilation for stent placement also denudes the endothelium layer and damages the vessel wall. It then triggers smooth muscle cell (SMC) activation, proliferation and migration and extracellular matrix formation, as well as other processes, e.g. inflammation and platelet aggregation [3,4]. This leads to the growth of neointima, which is newly formed tissue composed mainly of smooth muscle cells and extracellular matrix, in the vessel lumen. The excessive growth of neointima can result in a renarrowing of the vessel, a condition known as in-stent restenosis (ISR).
To study the mechanism of restenosis, a multiscale model for ISR was proposed [5] and a first two-dimensional version of that model (named ISR2D) was developed and studied in detail [6][7][8]. The model consists of three submodels: an initial condition (IC) model, an agent-based SMC model and a blood flow (BF) model. The IC model simulates balloon expansion and stent deployment and provides the input configuration for the other two models. The agent-based SMC model simulates the biological and mechanical states of each cell of the vessel, while the BF model provides the haemodynamics information as a function of the current vessel lumen shape. This multiscale model has been applied to investigate the effect of functional endothelium regeneration and the impact of stent deployment and design on restenosis [6,7,9,10]. Most recently, the effects of local blood flow dynamics with scenarios of adaptive and non-adaptive coronary vasculature on restenosis were studied based on the ISR2D model [11]. The two-dimensional model is, however, a simplification of the actual pathology. Therefore, a more comprehensive three-dimensional model (named ISR3D) was developed and compared with in vivo experimental scenarios [12,13].
Uncertainty quantification (UQ) is widely applied to study the effect of uncertainties in initial or boundary conditions and of other parameters of computational models on their simulated quantities of interest. Common UQ methods, such as those based on the Monte Carlo method [14][15][16], require a large number of simulations to provide enough data for the numerical integration of the statistical estimator [17]. However, it might be prohibitive for computationally expensive models, such as ISR3D, to achieve this. One solution could be to adopt surrogate modelling, by which a surrogate model (or metamodel) is developed to approximate the response of the original model at a relatively low cost. Subsequently, this surrogate model replaces the original simulation to realize the evaluations required for the UQ.
The construction of a surrogate model can be categorized into three types: simplified models, projection-based methods and data-fit methods [18]. Simplified models refer to a rough approximation based on simplifications of the simulated system such as spatial dimensionality reduction [19,20] or coarse-grid discretizations [21,22]. The projection-based methods proceed by identifying a low-dimensional subspace that is constructed to retain the essential character of the system input-output mapping. Stochastic collocation [23,24] and polynomial chaos expansion [25,26] are two state-ofthe-art projection-based methods for UQ analysis. Finally, the data-fit methods map out latent functions between input and output. Common methods for these types of surrogates are support vector machines [27], neural networks [28] or Gaussian processes (GPs) [29].
GP regression is widely applied in uncertainty estimation and reliability analysis owing to its non-parametric and Bayesian inference nature [30][31][32][33]. It was first proposed by Krige for geostatistics [34], and later extensively studied and extended to solve the regression problem under different scenarios, such as multi-task/multi-output GPs for vectorvalued function [35], heteroscedastic GPs for input-dependent noise scenarios [36][37][38], sparse GPs with inducing inputs for efficient training of large datasets [39,40] or deep GPs with a hierarchical structure to capture more complex processes [41].
Generally, GPs are designed for a scalar output and become cumbersome when multi-output is required as a result of the large kernel used for co-regionalization. The complexity of multi-output Gaussian process (MOGP) is associated with the dimension of output and the number of training samples. The computational cost of MOGP can easily become prohibitively expensive if the desired output dimension is high. One alternative solution is to apply dimensionality reduction techniques, such as proper orthogonal decomposition (POD) [42], to the model response before regression. The regression prediction is hence no longer the model response but the projection coefficients of the response. Owing to the limited amount of projection coefficients required for the reconstruction of the output space, several single-output GPs are sufficient in this case. This method has been widely applied for time-dependent problems [43,44], computational fluid dynamics [45], etc.
Here, the uncertainty propagation due to four uncertain parameters of the ISR3D model (endothelium regeneration time, the threshold strain for SMC bond breaking, blood flow velocity and the percentage of fenestration by area in the internal elastic lamina (IEL)) is investigated. The quantities of interest (QoIs) are the average cross-sectional area of the lumen and the maximum relative area loss as a function of time. We applied POD to reduce the dimension of the output and used GP regression as the surrogate model to map the uncertain inputs to the projection coefficients of the POD. With this computationally efficient surrogate model, uncertainty estimations and sensitivity analysis of the restenosis process are conducted and analysed.
The paper is arranged as follows. The details of the ISR3D model are introduced in §2. The construction of the surrogate model with POD and GPs is described in §3. The uncertain parameters and uncertainty estimations are presented in §4. The results of uncertainty estimates and sensitivity analysis are presented in §5 followed by a discussion in §6 and the conclusion in §7.

In-stent restenosis three-dimensional model
ISR3D is a multiscale computational model simulating the post-stenting neointima growth in a coronary artery [12,13]. It mainly consists of three single-scale submodels: the IC model, the SMC model (including details of the vascular wall, such as the lamina and the endothelial cells) and the BF model. A schematic diagram of ISR3D is shown in figure 1.
The SMC model has two parts: one deals with the biomechanics of the vessel wall post-stenting, while the other deals with the SMC biology, mainly in relation to proliferation and production of extracellular matrix. The mechanical part of the SMC model simulates the mechanical response of the vessel wall, based on cell-cell pairwise repulsive and attractive forces and calculating the cell displacements. Each SMC of the vessel wall is modelled as a spherical agent, and the interactions between them are provided by potential and bond forces. The effective radii of particles represent the radii of corresponding cells and changes during the growth governed by the biological solver [13].
The biological model of SMCs describes the cell cycle dynamics. Cell life cycle is a sequence of growth, replication and division of the cell; at the end of the life cycle, the cell divides into two daughter cells. The processes that influence the cell life cycle take place in the 30 μm neighbourhood around the cell; the time scale of one cycle is around 24-48 h.
The growth of individual SMCs is modelled by a cell cycle model, similar to the one described in [46]. Each cell can be in a state of growth (G1), synthesis/secondary growth/mitosis (S/G2/M) or a quiescent state (G0). Cells evolve from one state to the next according to their internal clock, and stop or die under the influence of external factors such as contact royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20210864 inhibition (the mechanical stresses in between SMCs) or the concentration of nitric oxide. The biological model provides new radii, states of the cells as its output and also the initial coordinates for newly formed cells. Growth of the neointima takes several dozens of cell cycles and stops several weeks after the stenting procedure [13].
The BF model is a pressure-driven fluid dynamics model, which provides relevant ranges of shear stresses on the vessel walls. The solver receives the lumen geometry every time step from the SMC solver, simulates the steady-state blood flow and returns the wall shear stress information to the SMC model. The blood is assumed to be incompressible and Newtonian, and is modelled by the lattice Boltzmann method (LBM) [47] in a three-dimensional rectangular mesh (D3Q19). The inlet boundary condition for velocity is set to a parabolic profile and its maximum velocity is defined as one of the uncertain parameters. A Dirichlet pressure boundary condition is assigned at the outlet and the vessel wall is defined as a non-slip condition. The simulation is implemented with Palabos [48].
The initial stent deployment is performed by the IC model. The stent is expanded radially with a capsule-shaped balloon until it reaches a predefined deployment depth. As there is no uncertain input of the UQ experiment related to the IC model and all the simulations start from exactly the same postdeployment state, we exclude the IC model from the execution of the UQ. For further details about the ISR3D, see [12,13]. A public version of the ISR3D model, which is studied in this paper, can be found on Github. 1 In the UQ experiments described here, the scenario of stenting a small porcine coronary vessel with 2 mm diameter is simulated. The simulation domain is limited on the outside by the middle layer of the artery wall, the tunica media. The entire length of the vessel is set to be 18 mm with a tunica media thickness of 0.35 mm and 1 mm lumen radius. The entire vessel is assumed to be slightly curved to obtain a more realistic blood flow pattern in the vessel. The stent applied in the simulations is made of intersecting spiral elements (shown in figure 2a). It can be viewed as a simplified version of the NIR stent [49] and the deployment depth is set to be 0.25 mm. The model is set to simulate the restenosis process up to 30 days after stenting.
The computational cost of ISR3D with a vessel and a stent of this size is rather expensive. A single run of the ISR3D simulation takes 500-600 core hours on a supercomputer node (a node with 2 × 12-core 2.6 GHz Intel Xeon E5-2690 v. 3 CPUs), depending on the total amount of neointima growth. For non-intrusive UQ methods, a large number of evaluations of the model are required for the statistical analysis and this becomes impractical for such a computationally expensive model. Therefore, to perform the UQ efficiently, a data-driven surrogate model based on GPs and POD is developed to learn the latent function between the uncertain inputs and the QoIs, and applied to evaluate the model response in the UQ.
Two QoIs are measured in the UQ experiment: the average cross-sectional area of the vessel lumen and the maximum relative area loss. The lumen cross-sectional areas along the centreline of the vessel are obtained using an open-source toolkit VMTK. 2 The average values of this area over the considered vessel model at each time step are used to evaluate how the uncertain parameters influence the total amount of neointima growth over time (shown in figure 2b). The relative area loss of the vessel shows the relative amount of neointima growth compared with the initial post-stenting cross-sectional area. Clinically, the restenosis is defined as the renarrowing of the lumen to more than 50% occlusion [3]. The maximum value of relative area loss of a vessel offers us a criterion to

Proper orthogonal decomposition on model response
Assume the response of the model is a series of responses (here, average cross-sectional areas of the lumen) over time y [ R Nt , where N t is the dimension of the output vector. The POD method can be applied to approximate the model responses by projecting the response to a low-rank space. The POD can be realized in three schemes: Karhunen-Loeve decomposition, principal component analysis and the singular value decomposition (SVD). In this work, the SVD method is applied for the decomposition [50]. Consider a snapshot matrix S [ R NtÂNs consisting of N s number of the model responses fy 1 , y 2 , . . . , y Ns g, The snapshot matrix can be decomposed into three matrices using singular value decomposition, where U and V denote the left and right orthonormal matrices, respectively. S denotes a diagonal matrix with singular values The objective of POD is to find out a set of orthogonal basesF ¼ ff 1 , f 2 , . . . , f k g from the space L ¼ fF [ R NtÂk :F T F ¼ Ig containing all possible orthogonal bases, such that the error introduced by the projection to low dimensional space could be minimized, where Á k k L 2 denotes the L 2 norm. By the Eckart-Young theorem [51], the orthogonal basis with the basis vectors fu i g k i¼1 taken from the ith column of U is the solution to such an optimization problem. The relative energy captured by the projection to such low-dimensional space consisting of the first k columns of U can be evaluated by [52] We assume that, if the relative energy R en is higher than 99.9%, the approximation reconstructed by the first k bases performs well enough. Since the values of σ i decay rapidly, a small k would be sufficient to achieve the relative energy threshold. Once the basis vectors are obtained, any model response can be approximated by: y %ŷ ¼ P k i¼1 a i f i , where α i are the projection coefficients.

Gaussian process regression
Assume that a model response y [ R is generated by the function y = f (x) + ϵ with a corresponding input x [ R d , and ϵ denotes the noise of the measurement or stochasticity of the model and assumes that a normal distribution is followed: N ð0, s 2 n Þ. A GP can be defined as a collection of random variables and any finite number of the random variables follows joint Gaussian distribution [29], f ðxÞ GPðmðxÞ, kðx, x 0 ÞÞ, ð3:5Þ where GP denotes a GP prior over the space of functions specified by its mean function m(x) and covariance functions k(x, x 0 ). Generally, the mean function is set to be zero without loss of generality. The kernel functions specify how the random variables are correlated with each other and also imply the smoothness of the functions. One of the common choices is the radial basis function kernel with automatic relevant determination (ARD) [29], ð3:6Þ where s 2 f is the signal variance and ℓ i denotes the length scale for each input dimension. For a regression royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20210864 problem, an independent Gaussian kernel with variance s 2 n is used to specify the noise in the function. These hyperparameters in the kernel will be determined via the optimization of likelihood function with observed data collection ðX, . To predict the model response at an unevaluated location x*, the GP prior can be rewritten into Conditioning on the observed data, the predictive distribution of the new point x* also follows a normal distribution, The y Ã stands for the mean prediction of the response and Var( y*) is the predictive variance indicating the uncertainty of the prediction. Generally, the GP regression is applied as a surrogate model to infer the latent function between uncertain inputs and QoIs. However, after the decomposition of the model response by SVD, both evaluated and unevaluated model responses can be represented by the projection coefficients on the chosen orthogonal bases, therefore the GP is now used to learn the mapping between uncertain inputs and projection coefficients of POD and predicts the new coefficients for unevaluated points. The details of the procedure are shown in algorithm 1.
Algorithm 1. Constructing a surrogate model for ISR3D with GP and POD. Training 1. Evaluate N train number of samples using ISR3D and collect the training data ðx i , y i Þ È É Ntrain i¼1 . 2. Construct the snapshot matrix ½y 1 jy 2 j Á Á Á jy Ns and perform SVD to obtain k orthogonal basesF based on the relative energy threshold. 3. Project the output of the training data to each basis and compute the projection coefficients:

Uncertainty quantification 4.1. Uncertain parameters
The four epistemic uncertainties considered in the forward uncertainty propagation of ISR3D include endothelium regeneration time, blood flow velocity, the threshold strain for SMC bond breaking and the percentage of fenestration by area in the IEL. Note that all the uncertain parameters except the blood flow velocity are parameters of the SMC submodel.

Endothelium regeneration time
The endothelium regeneration starts right after the denudation caused by the balloon dilation and stent deployment.
With sufficiently high wall shear stress from blood flow, the endothelium releases nitric oxide, which behaves as the inhibitor of the proliferation of SMCs. Therefore, the rate of endothelial regrowth significantly influences the growth of neointima. In the ISR3D, the regeneration of endothelium cells is modelled to increase linearly up to a coverage of 59% after 3 days, followed by a full recovery to 100% after a certain number of days given by the uncertain input [12]. This setting is based on experimental results from Nakazawa et al. [53]. However, the exact time for re-endothelialization may vary with many factors, such as the severity of vessel injury, the types of stenting and the degrees of inflammatory response [54]. In order to study this uncertain parameter, we consider an average endothelium regeneration time of 15 days based on the experimental data from Nakazawa et al. [53], and vary it from 10 to 20 days in the UQ.

Threshold relative strain
The threshold strain is the maximum strain that can be obtained before the bonds between SMCs break. Generally, during the stenting process, the vessel wall is overstretched in the circumferential direction, and therefore the connections between the SMCs (e.g. collagen fibres) are possibly broken and cause microfractures in the tissue. These microfractures may cause inflammation and contribute to the proliferation of SMCs after stenting. Our choice on the uncertainty of the breaking strain is inferred from stretching experiments [55,56] in which the mechanical responses of the coronary arteries under the stretch condition were gauged. The result demonstrated that the first intimal rupture occurred at around 110% strain, and the strain-stress curve became non-smooth when strain reached approximately 120%. Therefore, we consider the threshold strain around the experimental rupture value 1.1 with an uncertainty of +20% in our UQ experiment. Note that the measurements in [55,56] started from an unstrained sample, while in our model the vessel is prestrained by 30% owing to it being pressurized by the flowing blood inside it.

Blood flow velocity
Blood flow, as one of the mechanical factors, also plays an important role in the growth of neointima [57,58]. High enough wall shear stress in the vessel accelerates the production of nitric oxide in endothelial cells, which acts as an inhibitor of SMC proliferation.
As mentioned before, the blood flow in the simulation is modelled as a steady flow with a constant parabolic inlet boundary condition, since we are mainly interested in the time-averaged values of the wall shear stress. The velocity data from [59] were applied to compute time-averaged velocity and converted to the parabolic profiles, the maximum velocity of which is 0.266 m s −1 . Owing to the royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20210864 measurement error and potential variety of velocity for individual vessels, we presume a large uncertainty in the data and vary it by 50% based on the average values 0:266 m s À1 .

Fenestration percentage
The IEL is modelled in ISR3D as a layer of agents on the inner surface of the vessel wall [12]. The fenestrations on IEL significantly affect the initial growth of SMCs as they allow SMCs to migrate into the blood vessel and start proliferating there. However, the SMCs in ISR3D are not able to change shape to migrate through the fenestrations, unlike real SMCs. Therefore, a certain percentage of IEL agents are switched to SMCs in ISR3D, to obtain a smaller amount of very large fenestrations, with the same total surface area as in the experiment. The fenestration percentage is calculated as the total area of all fenestrations divided by the total area of the IEL. The uncertainty ranges for this parameter are obtained from [60], where the percentage of fenestration in the hypercholesterolaemic group is approximately 7.5% and in the control group is approximately 3.5%. To include and study the scenarios for both cases, we consider the parameters to vary from 2% to 10%.
The ranges of all the uncertain parameters mentioned above are given in table 1 and the distributions of the uncertainties are all assumed to be uniform.

Uncertainty estimations and sensitivity analysis
For the UQ we applied the quasi-Monte Carlo (qMC) sampling method with Sobol sequence [61]. The method allows the samples to be more evenly distributed in the domain, which leads to a better convergence rate than the standard random sampling.
To investigate the uncertainty propagation of the uncertain inputs through the model, mean, variance, probability density function (PDF) and coefficient of variation are estimated. In addition, global sensitivity analysis has been conducted to study how much each uncertain input has contributed to the uncertainty of QoIs. The variance-based method (Sobol method) [62] is applied, which assumes that the latent functions f (x) can be decomposed into a combination of functions of individual uncertain inputs and their higher order interactions, which also leads to the following decomposition of the variance [62]: where V i , V ij , V 12...d stand for the partial variance contributed by the ith uncertain input, by the interactions between the ith and jth uncertain inputs and by the higher order interactions. The first-order Sobol indices indicate the independent contributions from the partial variance of each single uncertain input, where x ∼i denotes a vector of all uncertain parameters in x except x i . The total sensitivity indices take all the relevant contributions of a uncertain input into account, Var fðxÞ : ð4:3Þ All the sensitivity indices mentioned above are computed by Saltelli's method [62].

Results
To train the surrogate models, 512 samples were generated by the qMC method and evaluated by the ISR3D model. Before the surrogate model was deployed to the UQ experiment, the surrogates were validated with a fourfold cross-validation. We measured the approximation error of both POD and GP regression with the relative L 2 norm, where N cv denotes the number of samples used for cross-validations. In the cross-validation of POD, a certain number of snapshots were randomly taken from the training dataset and used to construct the snapshot matrix for SVD. The validation dataset was used to measure the approximation error. The relative L 2 errors of the POD approximation with a different number of snapshots of both QoIs are shown in figure 3. The average relative L 2 error gradually decreases to around 0.07% and 0.3%, respectively, with the number of snapshots reaching 100. The tendency of the curve shows that the error has almost converged to a limit; a further increase in the number of snapshots will not greatly improve performance. The low standard deviation of the error shows that there is no significant influence on the choice snapshots. Therefore, we randomly chose 100 snapshots from the output of the training data for the POD in the construction of the surrogate model. To test the performance of the GP regression, another fourfold cross-validation was performed with 100 repetitions. Table 1. Ranges, units and coefficients of variation (CVs) of uncertain parameters of the ISR3D model. Note that the relative threshold strain is calculated with 30% pre-strained. The measurements in [55,56] started from an unstrained sample while in our model the vessel is pre-strained by 30%. Therefore, strain listed by Holzapfel et al. σ abso is scaled to obtain the relative deformation of our pre-strained tissue by σ rela = (((σ abso + 1)/1.  [59] percentage of fenestration 2 10 % 0.38 [60] royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20210864 The predicted projection coefficients were first used to reconstruct their original model responses and subsequently compared with the expected output from the validation dataset. Comparisons of the predicted QoIs versus expected QoIs over all time steps are demonstrated in figure 4. The resulting points are clustered around the diagonal line, indicating that the GP has inferred the underlying functions well. The average relative L 2 error is 0.52% for the average cross-sectional area and 2.52% for the relative maximum area loss. After the validations of surrogate models, the UQ experiments for both QoIs were performed. We applied the qMC method to draw 10 5 samples from the uncertain input domain and fed them to the surrogate models. The mean and 50%, 75% and 95% percentile estimations of average cross-sectional area over time are shown in figure 5. The corresponding histograms and PDFs of days 5, 10, 15, 20 and 30 are also shown in the same figure. The initial average crosssectional area after stenting was around 3.17 mm 2 . With the evolution of time, the cross-sectional area gradually reduced owing to the neointimal growth. The mean estimation of the average cross-sectional area shows that the neointimal growth was slow at the beginning but started to accelerate after day 1. An almost linear growth between day 1 and day 10 was observed followed by a descending growth rate until all the growth stopped at around day 22. The upper boundary of the 95% percentile shows that some samples stopped growing shortly after day 10 owing to the short reendothelialization time, while a few other cases did not stop before 22 days.
The PDFs and histograms in figure 5 show the details of the distributions of days 5, 10, 15, 20 and 30. On day 5, most of the samples cluster around 2.8 mm 2 and a small number of the samples have a lower average cross-sectional area up to 2.63 mm 2 . A certain number of samples already stopped growing between days 10 and 20. The early stop usually means a small amount of neointima and contributes to the right tail of the distributions (around 2.4-2.6 mm 2 ), while the rest of the samples still shifted towards the left owing to the growth. The difference between days 20 and day 30 is minor, indicating that the growth in most of the samples had stopped before day 20.
Similar patterns can be observed for the maximum relative area loss in figure 6. The distribution at day 30 shows that most of the simulations ended up with 30-60% area loss. Assuming that the restenosis happened when the area loss reached 50%, about 5%, 16% and 18% of the simulations had reached the restenosis threshold at days 15, 20 and 30, respectively. Table 2 provides detailed information on the mean, standard deviation (SD), coefficients of variation (CVs) at days 5, 10, 15, 20 and 30 of both QoIs computed by 100 replications of the UQ experiment. Around 11.3% and 16.6% uncertainties are observed from the average cross-sectional area and maximum relative area loss, respectively.
Apart from the uncertainty estimations, sensitivity analysis has also been performed. The sensitivity analysis was performed with 5 × 10 5 samples using the Sobol sequence and was repeated 100 times to compute the confidence interval. The first-order indices of the four uncertain inputs over    figure 7. The confidence interval can hardly be seen in the figure, indicating extremely small uncertainty in our sensitivity estimations. The firstorder indices quantify the direct effect of each uncertain input on the total variance of the QoIs. For both QoIs, a dominant influence of the fenestration percentage on the variance can be observed at the initial stage of the process and keeps decreasing over time. It has almost no impact after 10 days. The blood flow velocity is a critical factor on the growth throughout the entire process and shows significant influences on variance in between day 5 and day 10, and gradually falls to around 0.2, while the endothelium regeneration times show an increasing effect and play the most important role after 13 days. The sensitivity to threshold strain is relatively small compared with the other uncertain inputs. The total order indices measure both the independent and higher order interaction effect of an uncertain input on the total variance of QoIs. The total order indices of both QoIs are very similar to their first-order result, meaning that there is little higher order interaction between the uncertain inputs.
To further investigate the relations between uncertain inputs and restenosis, scatter distributions and histograms of the samples which reached the restenosis threshold at days 15, 20 and 30 are shown in figure 8. Note that the threshold strain is not shown in the figure since the sensitivity analysis result suggested that it is relatively not important in the process. The speed-up of the entire UQ experiment using the surrogate model has also been estimated. Table 3 shows the details of the computational cost, including the average core hour for model evaluation with ISR3D and surrogate model T ISR , training data generation T sample and surrogate training T train . Both training and prediction of a surrogate model were extremely fast. The most computationally expensive part was the generation of training data with ISR3D. The average core hour for each evaluation was around 585. Since in this case N UQ T ISR þ T train is negligible compared with T train , we find that the speed-up equals N UQ /512 ≈ 976.6.

Discussion
The result of surrogate modelling shows that the combination of POD and GP regression performs well. The decomposition and reconstruction of the model response with POD save the computational effort for regression and provide a convenient and consistent way to cover the entire model response over time. In this work, the snapshot matrix was constructed by 100 randomly chosen snapshots from the training dataset generated by qMC sampling. An adaptive sampling method [63] can be used to choose more representative snapshots with error estimations; however, this was unnecessary as a relatively large training dataset was available and the approximation error could be properly controlled.
The GP was then applied to infer the latent functions between uncertain inputs and projection coefficients of POD. In the cross-validation of the surrogate model, the relative L 2 error of the maximum relative area loss is slightly larger than the other QoI. This is mainly due to its way of computing relative area loss, which required a division of the initial cross-sectional area. The initial cross-sectional areas at each slice of the lumen are different and thus introduced the noise into the data. Therefore, the regression performance of such a QoI was slightly worse than the others.
For the UQ, around 11% and 16% of uncertainty are observed from the average cross-sectional area and maximum relative area loss, respectively. The uncertainties in the output are mainly contributed by fenestration percentage, blood flow velocity and endothelium recovery time. The fenestration percentage is important at the beginning because a larger amount of fenestrations allows more SMCs to migrate to the vessel lumen and proliferate. However such an impact drops sharply to almost 0 in the first 5 days, as the SMCs form a continuous layer over the IEL. Meanwhile the blood flow velocity starts to dominate the variance between day 5 and day 10. During day 5, re-endothelialization coverage varied from 63% to 67% and increased up to 73-87% by day 10, which means that if the wall shear stress is sufficiently high, a large percentage of cells at the lumen surface could already have their growth inhibited by nitric oxide. After day 10, the influence of the blood flow velocity drops gradually and is replaced by re-endothelialization. Figure 8 shows that, at the end of the simulations, the influences of fenestration percentage are relatively minor compared with the effect of blood flow velocity and endothelium regeneration time. This suggests that the scenarios with a high fenestration percentage, such as hypercholesterolaemia, might not have a high impact on restenosis probability if other parameters such as endothelium regeneration time can be strictly controlled.
In the design of the UQ experiment, the choices of uncertain ranges and distributions are described in §4.1. For all uncertain parameters, the average values and some of the uncertainty ranges are directly obtained or calculated from experimental results. However, owing to the lack of data on real population distributions of these parameters, the uncertainty ranges and their corresponding distributions are mainly defined by our approximations. They are simplified, and are likely to be different from the actual distribution.  Figure 7. The first-order and total Sobol sensitivity indices of both QoIs (ACSA, average cross-sectional area; MRAL, maximum relative area loss) and corresponding 95% confidence interval of each estimate based on 100 replica computing. Changing the input distribution naturally will lead to a different output distribution in the UQ. Our choice of uniform distributions and approximation of the ranges can illustrate the role of the input parameters and their effect on the restenosis process to a certain extent, but further UQ experiments based on actual distribution data of these parameters are still required and the outcome of which can also be applied to verify against clinical restenosis data. If such a realistic input distribution can be obtained, e.g. from a synthetic population, the trained surrogate can be used to perform the UQ quickly, as long as the population data lie in the same range used for generating the training sample.
In this work, we studied four biological uncertain parameters. We quantified their uncertainty propagation and sensitivity for two QoIs adapted for in silico models from clinically recognized metrics. This helps us to better understand the underlying contribution of these parameters to restenosis. In addition to the investigated biological factors, other factors and scenarios can also be studied via ISR3D; for example, variability in the stenting procedure, such as deployment depth or malapposition of the stent. Through the UQ analysis, the potential effect of such factors can be quantified and studied. Additionally, different scenarios, such as small/large vessel diameters and the tortuosity of   royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 18: 20210864 the stented vessel, can also significantly influence the outcome of a simulation. We leave the study of these factors, which all affect the initial shape of the stented vessel, to our future work. The ISR3D model itself has several limitations. First, it does not account for the inflammation processes, which are important during the early stages of post-stenting. Second, the geometry used in the UQ experiment is not based on any particular vessel, and instead is a piece of a perfectly cylindrical tube, and the stent fits the curvature of the vessel perfectly and is radially expanded in a uniform way. Additionally, a steady flow approximation was used for our BF model. The inlet was set to a parabolic profile, assuming the vessel's shape upstream of the inlet is unknown. It should be noted that the curvature affects the flow pattern to a great extent, leading to helical flows and affecting the wall shear stress also inside the stented region. The effects of the curvature are studied to some extent in our earlier publications; see, for example, fig. 9 in [13]. Similarly, finiteelement computational studies confirm the same effect of the curvature in [64]. These studies demonstrate that it is possible to capture the helical flow with a parabolic inlet condition in the regions of developed flow, but the flow pattern may be unrealistic close to the inlet.
A combination of the limitations listed above may result in the underestimation of restenotic growth. For example, Morton et al. [65] reported an area loss of 62% for porcine vessels of a similar diameter and deployment depth with an NIR stent, which is very close to the upper bound of the distribution predicted by the model. Nevertheless, the experimental values lie within the distribution, further confirming that the ranges selected for UQ reasonably overlap with the physiological ranges. There are also other limitations in the model we use, which are discussed in detail in [12,13].

Conclusion
The UQ and sensitivity analysis of a multiscale model ISR3D was performed. The uncertainty propagation from four parameters-endothelium regeneration time, threshold strain, percentage of fenestration and blood flow velocity-to two QoIs-the average cross-sectional area and the maximum relative area loss-are investigated. Owing to the high computational cost of ISR3D, surrogate modelling techniques were applied. The QoIs over time were, first, decomposed by POD and the resulting projection coefficients were learned by a GP regression model. Cross-validations are applied to validate the performance of the surrogate model. The surrogate model was subsequently deployed in the UQ experiment to replace the original model. The UQ and sensitivity analysis results showed that the blood flow velocity and endothelium regeneration time have a significant influence on the neointima growth and result in restenosis, while the impact of the fenestration percentage is limited and the threshold strain barely has any influence on the process.