Non-intrusive and semi-intrusive uncertainty quantification of a multiscale in-stent restenosis model
Dongwei Ye, Anna Nikishova, Lourens Veen, Pavel Zun, Alfons G. Hoekstra
NNon-intrusive and semi-intrusive uncertaintyquantification of a multiscale in-stent restenosis model
Dongwei Ye a, ∗ , Anna Nikishova a, ∗ , Lourens Veen b , Pavel Zun a,c,d ,Alfons G. Hoekstra a a Computational Science Lab, Informatics Institute, Faculty of Science, University ofAmsterdam, The Netherlands b Netherlands eScience Center, Amsterdam, The Netherlands c ITMO University, Saint Petersburg, Russia d Erasmus University Medical Center, Rotterdam, The Netherlands
Abstract
Uncertainty estimations are presented of the response of a multiscale in-stentrestenosis model, as obtained by both non-intrusive and semi-intrusive uncer-tainty quantification. The in-stent restenosis model is a fully coupled multiscalesimulation of post-stenting tissue growth, in which the most costly submodelis the blood flow simulation. Surrogate modelling for non-intrusive uncertaintyquantification takes the whole model as a black-box and maps directly fromthe three uncertain inputs to the quantity of interest, the neointimal area. Thecorresponding uncertain estimates matched the results from quasi-Monte Carlosimulations well. In the semi-intrusive uncertainty quantification, the most ex-pensive submodel is replaced with a surrogate model. We developed a surro-gate model for the blood flow simulation by using a convolutional neural net-work. The semi-intrusive method with the new surrogate model offered efficientestimates of uncertainty and sensitivity while keeping a relatively high accu-racy. It outperformed the result obtained with earlier surrogate models. It alsoachieved the estimates comparable to the non-intrusive method with a similarefficiency. Presented results on uncertainty propagation with non-intrusive and ∗ Corresponding authors
Email addresses: [email protected] (Dongwei Ye),
[email protected] (Anna Nikishova),
[email protected] (Lourens Veen),
[email protected] (Pavel Zun),
[email protected] (Alfons G. Hoekstra)
Preprint submitted to ArXiv September 2, 2020 a r X i v : . [ s t a t . C O ] S e p emi-intrusive metamodeling methods allow us to draw some conclusions on theadvantages and limitations of these methods. Keywords:
Uncertainty Quantification, Sensitivity Analysis, Surrogatemodelling, Semi-intrusive method, Gaussian process regression, Convolutionalneural network, Multiscale simulation
1. Introduction
Numerical simulations of real-world phenomena contribute to a better un-derstanding of these phenomena and to predicting the dynamics of the under-lying systems. Many natural phenomena occur across scales in space and time[1, 2, 3, 4, 5, 6]. As a result, multiscale models and simulations are widely used[7, 8, 1, 9, 10]. These multiscale models couple mathematical models of relevantprocesses on different spatial or temporal scales together relying on suitablescale bridging methods [11]. However, multiscale simulations can suffer fromsubstantial computational cost because of the high computational demands of,usually, the microscale simulations [10]. Uncertainty quantification (UQ) anal-ysis [12] applied to multiscale simulations adds additional substantial compu-tational burden since thousands of runs are required for good estimates of theuncertainties. Therefore simulations can become extremely time-consuming orimpractical, even on current state-of-the-art supercomputing infrastructure.To reduce the computational cost of multiscale UQ, we have recently pro-posed a set of semi-intrusive algorithms for multiscale UQ [13] and demonstratedtheir effectiveness for several multiscale UQ scenarios [13, 14]. Usually the out-put of a multiscale model is derived from a macroscale submodel, which in turnis implicitly determined by microscale dynamics to which it is coupled. Oneapproach from [13] relies on performing a Monte Carlo UQ on the macroscalesubmodel while replacing the most costly microscale submodel by a surrogatemodel. The surrogate model is trained on existing data to learn the mappingbetween input and output and then makes a prediction for a new experimentbased on the learned pattern. Replacing the expensive part of a model with2 relatively cheap surrogate can often significantly improve computational ef-ficiency, but comes at the cost of reduced accuracy. In [14], a physics-basedand a interpolation-based surrogate model were constructed to implement thesemi-intrusive UQ for the in-stent restenosis multiscale model [15, 16, 17]. Theresults were compared to black-box Monte Carlo results to demonstrate theefficiency improvement. However, a comparison between the semi-intrusive al-gorithm and non-intrusive UQ with a surrogate model replacing the completemultiscale model were not explored in that study. In this work, we first improvethe surrogate model of micro-scale model (the blood flow simulation) from [14]by using a convolutional neural network (CNN), due to its capability of patternrecognition and feature extraction [18, 19]. Additionally, a surrogate model fornon-intrusive UQ is designed to directly map the input parameters to the quan-tity of interest. The UQ estimations with both these methods are then carriedout and compared in terms of the UQ estimation and computational efficiency.The paper is arranged as follows. The two-dimensional multiscale model ofin-stent restenosis is shortly introduced in Section 2. The new surrogate modelfor the blood flow simulation and the surrogate model for in-stent restenosismodel are described in Section 2.2. The approach to estimate and analysethe uncertainty of the response of the in-stent restenosis model is explained inSection 3. The results of surrogate modelling, uncertainty quantification andsensitivity analysis are presented in Section 4. Sections 5 and 6 compare anddiscuss the UQ performance and summarise the obtained results.
2. Model
An arterial stenosis is the abnormal narrowing of an artery, usually due toaccumulation of fatty material in the walls and intimal thickening (atheroscle-rosis). In ischemic heart disease, a stenosis in a coronary artery limits bloodflow to the heart muscle, which can result in reduced heart function, shortnessof breath, chest pains or a heart attack. Coronary stenosis can be treated using3alloon angioplasty, in which a balloon is inserted into the artery via a catheterand inflated, which compresses the fatty plaque against the arterial wall. Duringthis procedure, a wire-mesh stent is deployed to keep the artery from recoilingback to the narrowed state. This procedure damages the vessel wall, and inparticular the endothelium, the innermost lining of the artery. This triggers ahealing response involving (amongst other processes) growth and proliferationof smooth muscle cells (SMCs) on the inside of the artery. In some cases, anexcessive growth response occurs, leading to a significant renewed narrowing ofthe artery inside of the stent. This is known as an in-stent restenosis (ISR), andis considered an adverse treatment outcome [20, 21, 22].The ISR2D model is a two-dimensional simulation of the post-stenting heal-ing response of an artery [16, 17], which is used here to test the proposedsemi-intrusive multiscale UQ algorithm. Note that a more realistic, but alsocomputationally much more expensive three dimensional version of the modelis available [23, 24]. The ISR2D model used in this paper consists of three sub-models: the IC submodel, which simulates initial conditions in the form of thestate of the artery immediately after stenting, the SMC submodel, which is anagent-based simulation of smooth muscle cell growth and endothelium recovery,and the blood flow submodel, which uses the Lattice Boltzmann method (LBM)to simulate blood flow through the artery. The structure of ISR2D is shown inFigure 1.Sufficiently high wall shear stress (WSS) at the arterial wall triggers anypresent endothelium to produce nitric oxide, which in turn inhibits the growthof the SMCs if it crosses a threshold value. Blood flow thus affects SMC growth,but in turn is also affected by it, as the proliferating SMCs change the geometryof the artery. Note that due to random placement of daughter cells when thegrowing SMCs divide, variability in the length of the cell cycle, and a randomspatial pattern of endothelium recovery, the SMC model is stochastic. Themain output of the model is the cross-sectional area of the neointima (the newtissue formed due to SMC proliferation) as a function of time after stenting. Aclinically recognised in-stent restenosis occurs if more than 50% of the original4 igure 1: Diagram of the ISR2D model: the initial conditions (IC) submodel provides theinitial geometry to the smooth muscle cells (SMC) single scale model where the geometry isfurther updated. SMC model calls the blood flow (BF) simulation, which provides the wallshear stress (WSS) for the updated geometry. The cycle continues until the final time step isreached, when the SMC model yields the final output of the cross-sectional area. cross-sectional area of the artery is covered by the neointima [25].The SMC growth occurs over a period of weeks, while blood flow adaptsmuch more rapidly to the changing geometry. ISR2D is therefore a multiscalemodel exhibiting time-scale separation. For every time step of the SMC model(one time step simulates one hour resulting in 1440 steps per 60 days of thetotal simulation time), the BF simulation is run to convergence for the currentgeometry, and the resulting WSS values are sent back to the SMC model, whichuses those in the model of nitric oxide production by endothelial cells.Figure 2 shows an input (domain map) and corresponding output (shearstress field) of the BF submodel. The domain map represents the geometry ofthe artery in the form of a binary grid, with 1 (shown in black) representing asolid grid cell, and 0 (shown in white) representing a fluid grid cell. The output5 .0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 mm −0.6−0.4−0.20.00.20.40.6 mm Input mm Output
Figure 2: Input and output of blood flow simulation. Left plot: 150 ×
150 binary geometrymap as the input for blood flow simulation. The black part is the vessel wall and the whitepart is the lumen (the fluid domain). Right plot: Simulated shear stress in the domain (in
P a ). of the BF simulation is a corresponding grid of shear stress values, which areset to 0 for solid domain cells, and set to the shear stress values computed byLBM for fluid domain cells. The shear stresses at the fluid-solid boundary layerare taken as the wall shear stresses and passed to the SMC model.The blood flow simulation is the computationally most expensive componentof the ISR2D model. It takes around 80% of the computational time, and thepotential gain in performance obtained by replacing it with a surrogate model inthe semi-intrusive UQ scenario is therefore highest. Nikishova et al. performeda semi-intrusive UQ analysis for the same model, comparing surrogates basedon nearest-neighbour interpolation and on simplified physics to a non-intrusiveblack box Monte Carlo approach [14]. The UQ estimates with physics surrogatehas improved the computational efficiency but the means of the cross-sectionalarea of the neointima resulting from the surrogate models were substantiallylower than the black-box Monte Carlo result. On the other hand, the UQestimates with nearest-neighbour interpolation is better but the correspondingspeedup was relatively poor. In this paper, an accurate surrogate model usingconvolutional neural network is proposed and demonstrated to result in accurateestimates of the uncertainties, while at the same time leading to the desired6eduction in computational cost for the multiscale UQ. As mentioned in the previous section, the blood flow simulation in the ISR2Dmodel is computationally expensive. To reduce the computational cost of themodel, a surrogate model to compute the required wall shear stresses is usedto replace the original blood flow model. Convolutional neural networks havebeen applied to fetch the features from irregular geometries in fluid dynamicsprediction [18, 26, 19]. In ISR2D application, the aim is to learn the latentnonlinear function between vessel wall geometry, blood flow velocity and wallshear stress.The mapping between input and output of the BF simulation can be con-sidered as a function f , which takes the geometry matrix ζ and the inlet bloodvelocity v in as input and produces a 2 × k -dimensional vector of WSS magni-tudes, τ wss as the output: τ wss = f ( ζ , v in ) , (2.1)where k = 150 is the grid size along x axis, ζ = ( ζ ij ) ∈ R k × k . The geometrymatrix ζ was used for CFD simulation. The surrogate model ˆ f replaces theoriginal blood flow model f ( ζ , v in ) and offers an approximate prediction of wallshear stress in a reduced amount of time.The CNN model follows the network structure proposed by [18] and wasoptimized to fit our application. The model consists of three parts: shape en-coding, nonlinear mapping and stress decoding, as shown in Figure 3. Theshape encoding layers extract the features of the geometry to the shape code.A fully connected (FC) layer then maps the shape code together with the bloodflow velocity to the stress code. The stress decoding part is responsible for amapping from the stress code to wall shear stress. In this surrogate model, thegeometry input was transformed from a binary map to a 2 × k array which indi-cates the locations of upper and lower fluid-solid boundaries. The convolutionlayers then take the information from both boundaries into account and predictthe shear stress on these boundaries. There are three convolution layers, a fully7onnected layer and three deconvolution layers deployed between the input layerand output layer. Each of them is followed by a rectifier linear unit (ReLU) asthe activation function. Besides, the output of each convolution layer is con-catenated to the corresponding deconvolution layer to help with the decodingprocess.The training data for the surrogate model comes from the runs of the ISR2Dmodel [27]. One run of ISR2D calls the LBM solver 1440 times (once per hourof simulated time). This means that with only a few runs of the simulation, aconsiderable amount of flow data for training is already available. We trainedthe surrogate model with the data from four runs of ISR2D simulation, hence5760 blood vessel geometries and wall shear stress distributions were used fortraining. The training optimization was based on the mean squared error lossfunction: L ( τ wss ) = 1 n n (cid:88) i =1 (cid:107) τ ( i ) wss − ˆ τ ( i ) wss (cid:107) (2.2)where n denotes the number of samples of the training set. The Adam optimizer
Figure 3: Diagram of the CNN model
This surrogate model is constructed to replace the whole ISR2D multiscalemodel for uncertainty quantification and sensitivity analysis. Let the multi-scale model function be defined by g ( ξ ) with ξ denoting a n -dimensional vectorconsisting of the stochastic variables and uncertain inputs of the model. Theresponse of the model is: y = g ( ξ ) (2.3)As mentioned before, the ISR2D model is a stochastic model and thus includesboth aleatory uncertainty and epistemic uncertainty. In the surrogate modelfor non-intrusive UQ, we assume that the aleatory uncertainty can be separatedfrom the function g , hence the expression for the surrogate model can be writtenas: ˆ y = h ( ξ ∼ ξ ∗ ) + (cid:15) ( ξ ∗ ) (2.4)where ˆ y denotes the response of the surrogate model, ξ ∗ is the aleatory uncer-tainty and ξ ∼ ξ ∗ is the parameter vector without stochastic variable. Assumingthat the aleatory uncertainty (cid:15) follows a normal distribution N (0 , σ ∗ ), such astochastic model can be easily quantified by Gaussian process regression (GPR).Note that ISR2D is a time-evolving model, and that we are not only interestedin the response at the final timestep, but also in the process. To avoid an extradimension of input, more precisely, a time t that will significantly increase thecomputational cost of training and prediction [29], a local surrogate model foreach time step is constructed. Therefore the expression can be rewritten as:ˆ y t = h t ( ξ ∼ ξ ∗ ) + (cid:15) t ( ξ ∗ ) , t = 1 , , ..., T, (2.5)9PR is based on the assumption of the joint Gaussian distribution betweentraining data and prediction mean: y train t h pred t ∼ N , K t + σ ∗ t I K ∗ t ( K ∗ t ) T K ∗∗ t + σ ∗ t I , t = 1 , , ..., T (2.6)where covariance matrices K t , K ∗∗ t and K ∗ t denote the correlation within train-ing data, within test data, and between these two respectively at each time step.Applying Bayesian inference, the posterior probability distribution of h pred t giventraining set( ξ train ∼ ξ ∗ , y train t ) follows a Gaussian process: P ( h pred t | ( ξ train ∼ ξ ∗ , y train t ) , ξ pred ∼ ξ ∗ ) = GP (¯ h pred t , Var( h pred t )) , (2.7)where ¯ h pred t = K ∗ t ( K t + σ ∗ t I ) − y train t and Var( h pred t ) = K ∗∗ t − K ∗ t ( K t + σ ∗ t I ) − ( K ∗ t ) T . The mean value offers the prediction and the variance represents the uncertaintyof this prediction. The radial basis function kernel and white noise kernel werechosen for the computation of covariance matrices. In each local GPR model,there are three hyperparameters associated to the kernels, length scale l t , signalvariance σ ft and noise variance (stochasticity) σ ∗ t . These hyperparameters aretrained by optimizing log marginal likelihood function:log p ( y train t | ξ ∼ ξ ∗ ) = −
12 ( y train t ) (cid:62) ( K t + σ ∗ t I ) − y train t −
12 log | K t + σ ∗ t I | − n π. (2.8)The training data comes from the qMC result from [27] and the size of thetraining data was chosen to match the speedup of semi-intrusive UQ. As thesemi-intrusive method gains a speedup of around 7, the speedup of the non-intrusive method is also designed to be around 7 for comparison. Therefore,150 ISR2D simulations are used for training. The Gaussian process surrogatemodel used in non-intrusive UQ was built using GPy [30].10 . Uncertainty quantification and sensitivity analysis Uncertainty Quantification is the analysis of the precision of a computationalmodel [31], including uncertainty in the model response ( forward problem ), andits input ( inverse problem ). Here, the forward propagation of uncertainty andthe precision of the ISR2D model output are studied.The semi-intrusive metamodelling method involves replacing a computation-ally expensive single-scale submodel with a surrogate, which produces an ap-proximation to the original single scale model result in a reduced time. Theinput uncertainty is propagated via the surrogate model in the same way asfor the non-intrusive method, by running an ensemble with parameter valuessampled from a distribution. Using the obtained samples of the model response,uncertainty in the model is estimated, analysing the probability density func-tion, mean and variance, as well as estimating the Sobol sensitivity indices[32, 33]. The mean and variance of the model responses at time step t are thenestimated by: E (ˆ y t ( ξ )) ≈ N N (cid:88) j =1 ˆ y t ( ξ j ) , Var (ˆ y t ( ξ )) ≈ N N (cid:88) j =1 ˆ y t ( ξ j ) − N N (cid:88) j =1 ˆ y t ( ξ j ) , (3.1)where N is the number of samples. The Sobol sensitivity index for the i -thparameter is defined by S totalξ i = Var totalξ i (ˆ y t ( ξ ))Var (ˆ y t ( ξ )) , (3.2)where the partial variance in the numerator is approximated by [34, 35]Var totalξ (ˆ y t ( ξ i )) ≈ N N (cid:88) j =1 (cid:16) ˆ y t ( ξ j ) − ˆ y t ( ξ j ∼ i , ξ j + Ni ) (cid:17) , (3.3)where ξ ∼ i is a vector of all elements of ξ except the i -th element, and g ( ξ j ∼ i , ξ j + Ni )denotes the model response with the same values of all inputs as for g ( ξ j ) exceptthe i -th input ξ i . 11t is important to note that since ISR2D is a stochastic model, the estimatesof the partial variance for the uncertain inputs include the aleatory uncertaintyas well. The aleatory uncertainty was estimated on its own using Eq. 3.3, takingadvantage of the fact that two model runs with the same input values give g ( ξ j )and g ( ξ j ∼ ξ ∗ , ξ j + Nξ ∗ ).
4. Results
First, the quality of the blood flow surrogate model is evaluated. Then, theresults of UQ and sensitivity analysis based on non-intrusive and semi-intrusivemethods are compared to a previously obtained reference solution reported in[27].
Before applying the blood flow surrogate model to semi-intrusive UQ anal-ysis, the quality of the surrogate model was evaluated. We used normalizedmean absolute error (NMAE) on the test dataset to evaluate the quality of thesurrogate model: NMAE = k (cid:107) τ wss − ˆ τ wss (cid:107) max {| τ wss |} × , (4.1)where (cid:107) · (cid:107) denotes L norm and max {| τ wss |} is the peak stress in the LBM.Figure 4 visualizes the wall shear stress prediction of the CNN surrogate modeland the LBM solver on a test dataset (CFD simulation results in one ISR2Dsimulation). The averaged NMAE of the surrogate model on the whole testdataset is around 3 . igure 4: (a) - (f) Predictions of wall shear stress distribution along the upper boundary byLBM and CNN surrogate models at 0, 5, 10, 15, 20, 30 days. The corresponding NMAEs are4.14%, 4.15%, 5.2%, 2.67%, 2.49% and 2.88% ncertain Parameter Range (min) Range (max) Unitinlet blood flow velocity 0.432 0.528 m/smaximum deployment depth 0.09 0.13 mmendothelium regeneration time 15 23 days Table 1: List of UQ parameters of ISR2D and their min and max values.
Uncertainty in the ISR2D model response is due to the model stochasticityand uncertainty in three model parameters. The ranges of the uncertain param-eters are shown in Table 4.2. These three uncertain inputs are assumed to beuniformly distributed within the given ranges. The model output of interest wasthe neointimal area as a function of time after stenting and its uncertainty wasestimated using non-intrusive method (NI) and semi-intrusive (SI) method withsurrogates of the blood flow micro model. We compared the uncertainty quan-tification result and sensitivity analysis result for similar values of UQ speedup.We also show the quasi-Monte Carlo (qMC) result from [27] as the referencesolution. The total number of model runs in both qMC and SI experiments was1024.Figure 5 shows the estimates of the mean and standard deviation with theqMC, the SI and NI methods. The mean is approximated especially well for thefirst 10 simulated days. After this point, slightly less average growth is observedin both SI and NI estimates than in the qMC results. The NI estimates slightlyoutperformed the SI estimates. The shape of the mean value of the neointimalarea is well approximated by both SI method and NI method. The results ofthe standard deviation from the SI and NI methods also show approximatelysimilar value to the qMC estimator.The comparison of the probability density functions (pdf) obtained withthree UQ methods at 5, 10, 15, 20 and 60 days after stenting are shown inFigure 6. A good fit of the pdfs is obtained at the early time points. For day5 and day 10, the two sample Kolmogorov-Smirnov (K-S) test for qMC and SI14dfs produces statistics of 0.04 and 0.07 respectively, while the statistics betweenqMC and NI stays at 0.05. At later time points, the K-S statistic is alwayssmaller for the results of NI estimates than with SI estimates, which means abetter fit. These plots also indicate the ratio of restenosis cases defined by 50%occlusion of the original lumen area [25] (shown as the vertical line). Non-zerovalues of this ratio are only observed at day 20 and 60. As expected, prediction ofa smaller growth in SI resulted in a relatively smaller predicted binary restenosisrate. The NI’s prediction is closer to the qMC result. At the final simulationtime point, NI predicted 10.3% restenosis occurrence while SI predicted 7.7%occurrence. A summary of uncertainty estimation at the final time step bydifferent methods is presented in Table 2. The NI estimations have the smallesterror in the estimation of the mean and the restenosis ratio. The SI with CNNresults have a smaller error than some other methods for each estimator. All theSI and NI results show a statistically significant underestimation of the meanvalue (two-value t-test, p < . Time ( days ) N e o i n t i m a l a r e a ( mm ) Mean and Standard Deviation qMC MeanqMC Mean ± SDSI MeanSI Mean ± SDNI MeanNI Mean ± SD
Figure 5: Mean and standard deviation of the ISR2D model output on the neointimal areawith quasi-Monte Carlo (qMC) and with the semi-intrusive (SI) method and non-intrusive(NI) method. Q Method Micro Model Mean Estimationx 10 − (mm ) Standard Deviationx 10 − (mm ) Restenosis Ratio(%)Value Error Value Error Value ErrorqMC LBM Table 2: Comparison of the estimates of means and standard deviation of neointimal growthand restenosis ratio with qMC, SI and NI methods. The indicated error is the absolutedifference from the reference qMC value. The four surrogate models for SIUQ are data-driven model I (DD I), data-driven model II (DD II), physics surrogate model (Phys) andconvolutional neural network model (CNN). See [14] for details on the Phys, and DD I andDD II surrogates. preserved. Since both the overall and the partial variances are overestimated,the error is significantly smaller in the estimation of the sensitivity indices,and the indices obtained by NI method with GP are close enough to the oneestimated by qMC. The variances estimated by the SI method are also closeto the qMC result and all within its confidence interval. As a result, the totalsensitivity indices are also well approximated with this type of method.
The main purpose of applying SI and NI methods is to speed up the simula-tion and reduce the computational cost while getting good enough estimates ofthe uncertainties. The speedup of the UQ analysis of these advanced methodswas computed as follows: S = N T ISR N T ISR ∗ + T train + T sample (4.2)where N is the number of runs of ISR2D simulation in the UQ analysis. T ISR is the execution time of an ISR2D simulation with LBM solver. T ISR ∗ is theexecution time of either an ISR2D simulation with a blood flow surrogate model16 .02.55.07.510.012.515.0 P r o b a b ili t y d e n s i t y K-S: 0.04 Restenosisthreshold (Rt)qMC:0.0%>Rt SI:0.0%>Rt0246810 P r o b a b ili t y d e n s i t y
10 days
K-S: 0.07 qMC:0.0%>Rt SI:0.0%>Rt02468 P r o b a b ili t y d e n s i t y
15 days
K-S: 0.11 qMC:0.0%>Rt SI:0.0%>Rt0246 P r o b a b ili t y d e n s i t y
20 days
K-S: 0.11 qMC:4.3%>Rt SI:1.1%>Rt0.0 0.2 0.4 0.6 0.8
Neointimal area (mm ) P r o b a b ili t y d e n s i t y
60 days
K-S: 0.13 qMC:12.2%>Rt SI:7.7%>Rt K-S: 0.04 Restenosisthreshold (Rt)qMC:0.0%>Rt NI:0.0%>RtK-S: 0.04 qMC:0.0%>Rt NI:0.0%>RtK-S: 0.04 qMC:0.0%>Rt NI:0.0%>RtK-S: 0.05 qMC:4.3%>Rt NI:4.1%>Rt0.0 0.2 0.4 0.6 0.8
Neointimal area (mm ) K-S: 0.05 qMC:12.2%>Rt NI:10.2%>Rt
Figure 6: Probability density function of the neointimal area at different simulation timesobtained with the quasi-Monte Carlo (QMC), the semi-intrusive method (SI) and non-intrusivemethod (NI). igure 7: Partial and total variances (left column) and total Sobol sensitivity indices (rightcolumn) with qMC in solid lines and the SI (top row) and NI (bottom row) in dashed lines.Each of the quantities for the uncertain inputs includes the aleatory uncertainty. The areaaround the qMC results is the 95% confidence interval obtained by a bootstrap test [36]. or the execution time of the surrogate model of non-intrusive method. T train denotes the training time for the surrogate model. T sample denotes the time togenerate the training data for the surrogate model. In Table 3, the executiontimes and resulting speedups of the SI and NI methods relative to the qMCmethod are evaluated, including previously reported SI results from [14] for18 QMethod MicroModel T ISR (min) T ISR ∗ (min) T micro (min) T train (min) T sample (min) N Speedupof UQqMC LBM / 50.9 37.8 / 894 1024 1.72SI DD II / 14.6 2.05 / 447 1024 5.94SI Phys / 11.9 0.08 / / 1024 7.51SI CNN / 12.8 0.26 9.9 357.6 1024 6.79NI / / 0.17 / 4.3 1 . × Table 3: Comparison of the computational time and corresponding speedup of different ap-proaches. The time value indicates the mean computational time obtained over N = 1024samples. T micro is the execution time of micro model (LBM/surrogate models) in one ISR2Dsimulation. The all computations were performed on the Distributed ASCI SupercomputerDAS5 [37] with Intel Haswell E5-2630-v3 CPU. comparison. Because of the light surrogate model, the SI approach with CNNwas approximately seven times faster than black-box qMC, an improvement ofmore than a factor three over the nearest-neighbour interpolation based surro-gate model. The simplified physics model was even faster, but was also the leastaccurate one, while the SI with CNN based surrogate model provided the bestuncertainty quantification and sensitivity estimates among the four surrogates(see [14] for details on the Phys, and DD I and DD II surrogates).
5. Discussion
The CNN surrogate model performed well regarding the wall shear stressprediction for the micro model. It takes advantage of convolution layers tofetch latent features in the geometry input and then uses the FC layer anddeconvolution layers to map the features to the wall shear stress prediction.Although the CNN surrogate model was able to predict the wall shear stressquite accurately, a small error still exists. This error introduced by the surrogatemodel then propagated through the iteration and led to the error in the UQestimation and restenosis prediction as shown in Figure 5 and table 2. Theaccuracy of the SIUQ estimation depends not only on the quality of the surrogate19odel but also on the structure of the multiscale simulation. However, the UQresult from the SI method suggests that the error is small enough to produceuncertainty and sensitivity estimates close to the ones obtained by qMC. Ofcourse, the UQ result can be further improved by a better CNN surrogate model,e.g by training with a larger dataset or constructing a deeper CNN structure.But such improvement in the surrogate model does not necessarily guaranteea significant improvement in the UQ estimation. A trial experiment has beenperformed on ISR2D with a better CNN surrogate model (NMAE ≈ T ISR / T macro which is around six for ISR2D. However the speedup of the SIUQ with CNNsurrogate model (Table 3) is even higher than the maximum expected value.This is because the computational cost of the SMC (macro) model varies. Itmainly depends on the number of agents in the model. As the ISR2D with a sur-rogate model underestimates the neointimal growth which means fewer agentsin the SMC model, the corresponding computational time is reduced from 14minutes to 12 minutes and even less. Thus, the computational burden of theSIUQ method comes mostly from the cost of the simulation multiplied by thenumber of samples and it is not affected visibly by the cost of obtaining thesurrogate model. Although the deployment of the surrogate model reduces thecomputational cost, the maximum speedup is still limited by the computationalcost of the other parts of the multiscale simulation. The main computationalcost of NIUQ stems from the generation of training data. Consequently, there is20o upper limit for the speedup of NIUQ as long as the surrogate model can beefficiently and effectively trained with a small amount of data, e.g by applyingactive learning [38, 39]. Of course, the exact amount of data required to train asurrogate model for NIUQ depends not only on the dimensionality of the inputbut also on the desired accuracy. Therefore, selecting which method to apply forachieving an efficient UQ analysis has to take into consideration the multiscalemodel itself and the requirement of the UQ analysis.By comparing the SI and NI results at the same computational efficiency,one can see that uncertainty and sensitivity analysis of SI were as good as theNI. SI has the additional advantage of granularity, since only part of the modelis replaced by the surrogate. This means that the parameters of the submodelsnot replaced by the surrogate can be varied and studied without changing thesurrogate, as long as the replaced micro model is not affected. For example, incase of the ISR2D model, different parameters and rule sets for cell behaviourcan be used with the existing surrogate model for flow. On the other hand, usinga NI model for a different biological ruleset would require essentially building anew NI surrogate, which would incur a significant computational cost.In general, both SI and NI approaches performed well. The SI approach ismore suitable for cyclic multiscale simulations as it retains the framework ofthe simulation and can obtain the training data for the surrogate model at arelatively low cost. Another semi-intrusive UQ strategy is to run the simulationsfor UQ and build the surrogate model on the fly. Such a dynamic system shouldbe capable of constructing and updating the surrogate during the simulationprocess, for instance, as done by Leiter et al. [40]. In the UQ scenario, thesurrogate model can also be validated dynamically [13]. We aim to apply thesetechniques to the three-dimensional versions of the ISR model in future work.
6. Conclusion
In order to implement uncertainty quantification and sensitivity analysisefficiently for the ISR2D model, a new surrogate model based on a convolutional21eural network was developed and applied in semi-intrusive UQ estimation. TheUQ estimate with the new surrogate model was compared with the result ofprevious work and a non-intrusive UQ estimates based on Gaussian processregression. The result shows that SI with the convolutional neural networksurrogate model outperformed the previous result. The result is also comparableto non-intrusive estimates. Both SI and NI are valid methods to perform UQin an efficient way for the ISR2D model.
7. Funding
This work was supported by the Netherlands eScience Center under thee-MUSC (Enhancing Multiscale Computing with Sensitivity Analysis and Un-certainty Quantification) project. This project has received funding from theEuropean Union Horizon 2020 research and innovation programme under grantagreements
ReferencesReferences [1] A. Hoekstra, B. Chopard, P. Coveney, Multiscale modelling and simulation:a position paper, Phil. Trans. R. Soc. A 372 (2021) (2014) 20130377. doi:10.1098/rsta.2013.0377 .[2] D. Groen, S. J. Zasada, P. V. Coveney, Survey of multiscale and multi-physics applications and communities, Computing in Science & Engineering16 (2) (2014) 34–43. doi:10.1109/MCSE.2013.47 .223] A. Mizeranschi, D. Groen, J. Borgdorff, A. G. Hoekstra, B. Chopard,W. Dubitzky, Anatomy and Physiology of Multiscale Modeling and Simu-lation in Systems Medicine, Springer New York, New York, NY, 2016, pp.375–404. doi:10.1007/978-1-4939-3283-2_17 .[4] S. Alowayyed, D. Groen, P. V. Coveney, A. G. Hoekstra, Multiscalecomputing in the exascale era, Journal of Computational Science 22 (2017)15 – 25. doi:10.1016/j.jocs.2017.07.004 .URL [5] B. Chopard, J.-L. Falcone, P. Kunzli, L. Veen, A. Hoekstra, Multiscalemodeling: recent progress and open questions, Multiscale and Multidis-ciplinary Modeling, Experiments and Design 1 (1) (2018) 57–68. doi:10.1007/s41939-017-0006-4 .URL https://doi.org/10.1007/s41939-017-0006-4 [6] A. G. Hoekstra, B. Chopard, D. Coster, S. P. Zwart, P. V. Coveney,Multiscale computing for science and engineering in the era of exascaleperformance, Philosophical Transactions of the Royal Society A: Mathe-matical, Physical and Engineering Sciences 377 (2142) (2019) 20180144. doi:10.1098/rsta.2018.0144 .URL https://royalsocietypublishing.org/doi/abs/10.1098/rsta.2018.0144 [7] M. Praprotnik, L. D. Site, K. Kremer, Multiscale simulation of soft mat-ter: From scale bridging to adaptive resolution, Annual Review of Phys-ical Chemistry 59 (1) (2008) 545–571, pMID: 18062769. doi:10.1146/annurev.physchem.59.032607.093707 .[8] P. M. Sloot, A. G. Hoekstra, Multi-scale modelling in computationalbiomedicine, Briefings in bioinformatics 11 (1) (2010) 142–152. doi:10.1093/bib/bbp038 . 239] S. Karabasov, D. Nerukh, A. Hoekstra, B. Chopard, P. V. Coveney,Multiscale modelling: approaches and challenges, Phil. Trans R. Soc. A372 (2021) (2014) 20130390. doi:10.1098/rsta.2013.0390 .URL https://royalsocietypublishing.org/doi/abs/10.1098/rsta.2013.0390 [10] S. Alowayyed, D. Groen, P. V. Coveney, A. G. Hoekstra, Multiscalecomputing in the exascale era, Journal of Computational Science doi:10.1016/j.jocs.2017.07.004 .[11] B. Chopard, J. Borgdorff, A. G. Hoekstra, A framework for multi-scalemodelling, Phil. Trans R. Soc. A 372 (2021) (2014) 20130378. doi:10.1098/rsta.2013.0378 .[12] R. C. Smith, Uncertainty quantification: theory, implementation, and ap-plications, Vol. 12, Siam, 2013.[13] A. Nikishova, A. G. Hoekstra, Semi-intrusive uncertainty propagation formultiscale models, Journal of Computational Science 35 (2019) 80 – 90. doi:10.1016/j.jocs.2019.06.007 .URL [14] A. Nikishova, L. Veen, P. Zun, A. G. Hoekstra, Semi-intrusive multiscalemetamodelling uncertainty quantification with application to a model ofin-stent restenosis, Phil. Trans R. Soc. A 377 (2142) (2019) 20180154. doi:10.1098/rsta.2018.0154 .URL https://royalsocietypublishing.org/doi/abs/10.1098/rsta.2018.0154 [15] D. J. Evans, P. V. Lawford, J. Gunn, D. Walker, D. Hose, R. Smallwood,B. Chopard, M. Krafczyk, J. Bernsdorf, A. Hoekstra, The application ofmultiscale modelling to the process of development and prevention of steno-sis in a stented coronary artery, Phil. Trans R. Soc. A 366 (1879) (2008)3343–3360. 2416] H. Tahir, A. G. Hoekstra, E. Lorenz, P. V. Lawford, D. R. Hose, J. Gunn,D. J. W. Evans, Multi-scale simulations of the dynamics of in-stent resteno-sis: impact of stent deployment and design, Interface Focus 1 (3) (2011)365–373. doi:10.1098/rsfs.2010.0024 .[17] H. Tahir, C. Bona-Casas, A. G. Hoekstra, Modelling the effect of a func-tional endothelium on the development of in-stent restenosis, PLoS One8 (6) (2013) e66138. doi:10.1371/journal.pone.0066138 .[18] X. Guo, W. Li, F. Iorio, Convolutional neural networks for steady flowapproximation, in: Proceedings of the 22Nd ACM SIGKDD InternationalConference on Knowledge Discovery and Data Mining, KDD ’16, ACM,New York, NY, USA, 2016, pp. 481–490. doi:10.1145/2939672.2939738 .URL http://doi.acm.org/10.1145/2939672.2939738 [19] L. Liang, M. Liu, C. Martin, W. Sun, A deep learning approach to es-timate stress distribution: a fast and accurate surrogate of finite-elementanalysis, Journal of The Royal Society Interface 15 (138) (2018) 20170844. doi:10.1098/rsif.2017.0844 .URL https://royalsocietypublishing.org/doi/abs/10.1098/rsif.2017.0844 [20] J. W. Jukema, J. J. W. Verschuren, T. A. N. Ahmed, P. H. A. Quax,Restenosis after PCI. Part 1: pathophysiology and risk factors., Naturereviews. Cardiology 9 (1) (2012) 53–62. doi:10.1038/nrcardio.2011.132 .URL [21] J. W. Jukema, T. A. N. Ahmed, J. J. W. Verschuren, P. H. A. Quax,Restenosis after PCI. Part 2: prevention and therapy., Nature reviews.Cardiology 9 (2) (2012) 79–90. doi:10.1038/nrcardio.2011.148 .URL [22] J. Iqbal, P. W. Serruys, D. P. Taggart, Optimal revascularization for com-plex coronary artery disease., Nature reviews. Cardiology 10 (11) (2013)2535–47. doi:10.1038/nrcardio.2013.138 .URL [23] P. S. Zun, T. Anikina, A. Svitenkov, A. G. Hoekstra, A comparison offully-coupled 3d in-stent restenosis simulations to in-vivo data, Frontiersin Physiology 8 (2017) 284. doi:10.3389/fphys.2017.00284 .URL [24] P. S. Zun, A. J. Narracott, C. Chiastra, J. Gunn, A. G. Hoekstra,Location-Specific Comparison Between a 3D In-Stent Restenosis Modeland Micro-CT and Histology Data from Porcine In Vivo Experiments,Cardiovascular Engineering and Technology 10 (4) (2019) 568–582. doi:10.1007/s13239-019-00431-4 .URL http://link.springer.com/10.1007/s13239-019-00431-4 [25] P. W. Serruys, J. A. Ormiston, Y. Onuma, E. Regar, N. Gonzalo, H. M.Garcia-Garcia, K. Nieman, N. Bruining, C. Dorange, K. Miquel-Hbert,et al., A bioabsorbable everolimus-eluting coronary stent system (absorb):2-year outcomes and results from multiple imaging methods, The Lancet373 (9667) (2009) 897910. doi:10.1016/s0140-6736(09)60325-1 .[26] J. Tompson, K. Schlachter, P. Sprechmann, K. Perlin, Accelerating eule-rian fluid simulation with convolutional networks, in: Proceedings of the34th International Conference on Machine Learning - Volume 70, ICML’17,JMLR.org, 2017, pp. 3424–3433.URL http://dl.acm.org/citation.cfm?id=3305890.3306035 [27] A. Nikishova, L. Veen, P. Zun, A. G. Hoekstra, Uncertainty quantificationof a multiscale model for in-stent restenosis, Cardiovascular Engineeringand Technology 9 (4) (2018) 761–774. doi:10.1007/s13239-018-00372-4 .[28] F. Chollet, et al., Keras, https://keras.io (2015).2629] C. E. Rasmussen, C. K. I. Williams, Gaussian Processes for Machine Learn-ing (Adaptive Computation and Machine Learning), The MIT Press, 2005.[30] GPy, GPy: A gaussian process framework in python, http://github.com/SheffieldML/GPy (since 2012).[31] W. L. Oberkampf, S. M. DeLand, B. M. Rutherford, K. V. Diegert,K. F. Alvin, Error and uncertainty in modeling and simulation,Reliability Engineering & System Safety 75 (3) (2002) 333 – 357. doi:10.1016/S0951-8320(01)00120-X .URL [32] I. M. Sobol, On sensitivity estimation for nonlinear mathematical models,Matematicheskoe Modelirovanie 2 (1) (1990) 112–118.[33] S. Kucherenko, M. Rodriguez-Fernandez, C. Pantelides, N. Shah, MonteCarlo evaluation of derivative-based global sensitivity measures, ReliabilityEngineering & System Safety 94 (7) (2009) 1135 – 1148, Special Issue onSensitivity Analysis. doi:10.1016/j.ress.2008.05.006 .[34] T. Homma, A. Saltelli, Importance measures in global sensitivity analysisof nonlinear models, Reliability Engineering & System Safety 52 (1) (1996)1 – 17. doi:10.1016/0951-8320(96)00002-6 .[35] A. Saltelli, P. Annoni, I. Azzini, F. Campolongo, M. Ratto, S. Tarantola,Variance based sensitivity analysis of model output. Design and estimatorfor the total sensitivity index, Computer Physics Communications 181 (2)(2010) 259–270. doi:10.1016/j.cpc.2009.09.018 .[36] G. E. B. Archer, A. Saltelli, I. M. Sobol, Sensitivity measures, ANOVA-likeTechniques and the use of bootstrap, Journal of Statistical Computationand Simulation 58 (2) (1997) 99–120. doi:10.1080/00949659708811825 .[37] H. Bal, D. Epema, C. de Laat, R. van Nieuwpoort, J. Romein, F. Seinstra,C. Snoek, H. Wijshoff, A medium-scale distributed system for computer27cience research: Infrastructure for the long term, Computer 49 (5) (2016)54–63. doi:10.1109/MC.2016.127 .[38] T. Collet, O. Pietquin, Optimism in active learning with gaussian processes,in: S. Arik, T. Huang, W. K. Lai, Q. Liu (Eds.), Neural Information Pro-cessing, Springer International Publishing, Cham, 2015, pp. 152–160.[39] D. R. Jones, M. Schonlau, W. J. Welch, Efficient global optimization ofexpensive black-box functions, Journal of Global optimization 13 (4) (1998)455–492.[40] K. W. Leiter, B. C. Barnes, R. Becker, J. Knap, Accelerated scale-bridgingthrough adaptive surrogate model evaluation, Journal of ComputationalScience 27 (2018) 91 – 106. doi:10.1016/j.jocs.2018.04.010 .URL