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-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-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-of-the-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-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 vector-valued function [35], heteroscedastic GPs for input-dependent noise scenarios [36-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.

Figure 1

A schematic diagram of the ISR3D model. After the initial deployment of the stent with the IC model, the SMC and BF models run iteratively until the end of the simulation. At every time step, the SMC model passes the current lumen geometry to the BF model, which then updates the blood flow and sends the wall shear stress back to the SMC model. The SMC model then computes neointimalumen growth for the next time step, based on the wall shear stress.

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 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 post-deployment 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.

Figure 2

The simulation outcomes of ISR3D model (a) and corresponding QoIs measured over time (b). In (a), the blue part denotes the vessel wall, the beige part is the stent and the red part denotes the neointima. In (b), the average cross-sectional area and maximum relative area loss of the vessel at each day are measured and computed.

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 judge whether the restenosis happens or not. Note that both QoIs are evaluated as a function of time. The values at consecutive time steps are highly correlated with each other.

Surrogate modelling

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 , where N 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 consisting of N number of the model responses ,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. denotes a diagonal matrix with singular values σ, where i = 1, …, N and .

The objective of POD is to find out a set of orthogonal bases from the space containing all possible orthogonal bases, such that the error introduced by the projection to low dimensional space could be minimized,where denotes the L2 norm. By the Eckart–Young theorem [51], the orthogonal basis with the basis vectors 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 Ren is higher than 99.9%, the approximation reconstructed by the first k bases performs well enough. Since the values of σ 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: , where α are the projection coefficients.

Gaussian process regression

Assume that a model response is generated by the function y = f() + ε with a corresponding input , and ε denotes the noise of the measurement or stochasticity of the model and assumes that a normal distribution is followed: . A GP can be defined as a collection of random variables and any finite number of the random variables follows joint Gaussian distribution [29],where denotes a GP prior over the space of functions specified by its mean function m() and covariance functions k(, ′). 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],where is the signal variance and ℓ denotes the length scale for each input dimension. For a regression problem, an independent Gaussian kernel with variance 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 ,where = {σ, σ, ℓ1, …, ℓ} and = k(X, X). To predict the model response at an unevaluated location *, the GP prior can be rewritten intoConditioning on the observed data, the predictive distribution of the new point also follows a normal distribution,whereandThe 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 number of samples using ISR3D and collect the training data .

2. Construct the snapshot matrix and perform SVD to obtain orthogonal bases based on the relative energy threshold.

3. Project the output of the training data to each basis and compute the projection coefficients: , where .

4. Train th single-output GP with uncertain inputs and projection coefficients .

Prediction

1. For an unevaluated point , use GPs to predict its projection coefficients , where .

2. Reconstruct the corresponding model response .

Uncertainty quantification

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 in our UQ experiment. Note that the measurements in [55,56] started from an unstrained sample, while in our model the vessel is pre-strained 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 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 .

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.

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.3) − 1).

uncertain parameters	ranges (min)	ranges (max)	unit	CV	reference
endothelium regeneration	10	20	day	0.19	[53]
blood flow velocity	0.133	0.399	m s−1	0.29	[55,56]
relative threshold strain	0.446	0.785	–	0.16	[59]
percentage of fenestration	2	10	%	0.38	[60]

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() 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, V, 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 ∼ denotes a vector of all uncertain parameters in except x. The total sensitivity indices take all the relevant contributions of a uncertain input into account,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 L2 norm,where N 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 L2 errors of the POD approximation with a different number of snapshots of both QoIs are shown in figure 3. The average relative L2 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.

Figure 3

Relative L2 error of the POD with different numbers of snapshots for two QoIs, average cross-sectional area (a) and maximum relative area loss (b) using cross-validation. The error bar stands for the standard deviation computed from 100 replications.

To test the performance of the GP regression, another fourfold cross-validation was performed with 100 repetitions. 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 L2 error is 0.52% for the average cross-sectional area and 2.52% for the relative maximum area loss.

Figure 4

A comparison between expected QoIs and predicted QoIs over all time steps from a cross-validation, where the diagonal lines denote the precise predictions of expected values. ACSA, average cross-sectional area; MRAL, maximum relative area loss.

After the validations of surrogate models, the UQ experiments for both QoIs were performed. We applied the qMC method to draw 105 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 cross-sectional area after stenting was around 3.17 mm2. 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 re-endothelialization time, while a few other cases did not stop before 22 days.

Figure 5

Mean, 50%, 75% and 95% percentiles of the average cross-sectional area of the lumen over time with qMC sampling and corresponding histograms and PDFs at days 5, 10, 15, 20 and 30.

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 mm2 and a small number of the samples have a lower average cross-sectional area up to 2.63 mm2. 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 mm2), 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.

Figure 6

Mean, 50%, 75% and 95% percentiles of the maximum relative area loss over time with qMC sampling and corresponding histograms and PDFs at days 5, 10, 15, 20 and 30.

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.

Table 2

Mean, standard deviation (SD), CV (in percentage) and percentage of restenosis at days 5, 10 ,15, 20 and 30 for both QoIs computed from 100 repetitions of the UQ estimation. Owing to the large number of samples used, the confidence interval of the evaluations for each estimate is extremely small and is not shown here.

	average cross-sectional area (mm2)					maximum relative area loss (%)
estimates	day 5	day 10	day 15	day 20	day 30	day 5	day 10	day 15	day 20	day 30
mean	2.774	2.303	2.062	2.002	1.991	19.848	33.674	40.962	42.721	42.958
SD	0.046	0.098	0.179	0.224	0.226	1.023	2.745	5.672	7.037	7.124
CV	1.658%	4.255%	8.681%	11.189%	11.351%	5.154%	8.152%	13.847%	16.472%	16.591%
restenosis	–	–	–	–	–	0%	0%	5.123%	16.047%	17.873%

Apart from the uncertainty estimations, sensitivity analysis has also been performed. The sensitivity analysis was performed with 5 × 105 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 time for both QoIs are shown in figure 7. The confidence interval can hardly be seen in the figure, indicating extremely small uncertainty in our sensitivity estimations. The first-order 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.

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.

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.

Figure 8

Restenosis sample scatter distributions and corresponding histograms at days 15, 20 and 30.

In the left column of figure 8, scatter distributions of samples in terms of fenestration percentage and re-endothelialization time are shown. The range of re-endothelialization time falls between days 14 and 20, meanwhile a clear degression tendency can be observed from the corresponding histogram. The range of fenestration percentage shows that the restenosis can happen even with the lowest fenestration percentage, but the probability decreases slightly as the percentage drops. The middle column demonstrates the scatter distributions of re-endothelialization time and blood flow velocity. At day 15, only the cases with rather low blood flow velocity (under 0.27 m s−1) reached the restenosis threshold. However at the end of the simulations, the upper bound rose to 0.38 m s−1. Unlike the left and right columns, a clear separation can be found between restenotic samples and the rest of the domain. The right column is based on fenestration percentage and blood flow velocity. Same patterns could be observed. The influence of fenestration percentage is rather minor while the value of blood flow velocity significantly affects the possibility of restenosis.

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 , training data generation and surrogate training . 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 is negligible compared with , we find that the speed-up equals NUQ/512 ≈ 976.6.

Table 3

The computation cost of the ISR3D model and the surrogate model. denotes the core hour to finish one single run of the simulation. Each simulation was performed exclusively on a node with 2 × 12-core 2.6 GHz Intel Xeon E5-2690 v3 (Haswell) CPUs of Dutch supercomputer Cartesius. and stand for the sample generation time and training time for the surrogate model, respectively. NUQ denotes the number of samples used in a UQ experiment.

UQ method	T_ISR (core hour)	T_train (core hour)	T_sample (core hour)	NUQ	speed-up of UQ
qMC	585.1	–	–	–	1
qMC (surrogate)	6.1 × 10−8	3.1 × 10−2	585.1 × 512	5 × 105	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 L2 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. 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 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, finite-element 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.