Uncertainty Quantification of a Multiscale Model for In-Stent Restenosis

Purpose Coronary artery stenosis, or abnormal narrowing, is a widespread and potentially fatal cardiac disease. After treatment by balloon angioplasty and stenting, restenosis may occur inside the stent due to excessive neointima formation. Simulations of in-stent restenosis can provide new insight into this process. However, uncertainties due to variability in patient-specific parameters must be taken into account. Methods We performed an uncertainty quantification (UQ) study on a complex two-dimensional in-stent restenosis model. We used a quasi-Monte Carlo method for UQ of the neointimal area, and the Sobol sensitivity analysis (SA) to estimate the proportions of aleatory and epistemic uncertainties and to determine the most important input parameters. Results We observe approximately 30% uncertainty in the mean neointimal area as simulated by the model. Depending on whether a fast initial endothelium recovery occurs, the proportion of the model variance due to natural variability ranges from 15 to 35%. The endothelium regeneration time is identified as the most influential model parameter. Conclusion The model output contains a moderate quantity of uncertainty, and the model precision can be increased by obtaining a more certain value on the endothelium regeneration time. We conclude that the quasi-Monte Carlo UQ and the Sobol SA are reliable methods for estimating uncertainties in the response of complicated multiscale cardiovascular models.


INTRODUCTION
Cardiac diseases are a leading cause of mortality in developed countries. 16 Coronary artery stenosis, or abnormal narrowing, is a particularly widespread and potentially fatal cardiac disease. This narrowing is of-ten corrected by stenting the affected artery. There are multiple types of stents currently in use, ranging from simple bare metal stents, to more advanced stents, for example ones eluting growth-inhibiting drugs, ones designed to capture endothelial progenitor cells, and to bioresorbable polymer scaffolds. 10 In-stent restenosis (ISR) is the process of excessive neointima formation in an artery following balloon angioplasty and stenting, leading to a renewed narrowing of the artery. In 5% to 10% cases it requires a repeat revascularization of the target lesion. 7 Restenosis is caused by damage to the vessel wall and by disturbed flow patterns in the stented segment. 10 Restenosis is an important complication of the stenting procedure, and because of that, it has been studied extensively in clinical trials, reviewed e.g., in Ref. 8, invivo animal experiments, reviewed e.g., in Ref. 9, and also by using computational models. 3,5,13,18,34 The computational models of ISR usually represent cells by agents, either freely moving 5,18,34 or placed on a lattice. 3,13 These agents take cues from the blood flow, concentration of drugs eluted from the stent, and from the vessel damage and mechanics, which affect the growth and proliferation of the cells.
We perform an Uncertainty Quantification (UQ) of an off-lattice in-stent restenosis model previously developed by Tahir et al., 30 where we measure the precision of the model response, not the accuracy. The model is subject to both aleatory and epistemic uncertainty. It includes random variables which represent the stochastic nature of the system. It also includes uncertain parameters, which theoretically can be known precisely, although this is a rare case in practice. Here, we study our two-dimensional version of the ISR model. 5,29,30 We do this to provide a proofof-concept, using a model that is computationally relatively cheap, in preparation for a more in-depth study of our later, more physiological and more computationally expensive three-dimensional ISR model. 34 Our ultimate goal was to access the model result sensitivity to the uncertainties in the inputs and stochastic parameters, in order to be aware of the result uncertainty and the inputs which cause a greater effect on it. This will allow us to be aware of the result uncertainty quantities, when model results are analyzed, and of the inputs, which value should be as precise as possible.
A model description is given in Sect. 2. Next, in Sect. 3, we give a brief overview of the UQ method we used for measuring uncertainties, as well as of the Sobol method for the analysis of the effect of uncertain inputs on the model result. In Sect. 4, we present the results, and we finish with a discussion and conclusions in Sect. 5.

THE ISR MODEL
The ISR2D model is a two-dimensional representation of the ISR process. 6,30 Its domain is a longitudinal section of the artery, in which five subprocesses take place: stent deployment, post-deployment smooth muscle cell proliferation into the lumen, blood flow through the lumen, re-endothelialization, and diffusion of antiproliferative drugs from the stent into the tissue. As these processes take place at different temporal scales, ISR2D is a multiscale model. [4][5][6] Initially, the simulated artery consists of lower and upper arterial walls (see Fig. 1a). These comprise a tunica media, consisting of multiple layers of smooth muscle cells (SMCs), and a layer of Internal Elastic Lamina (IEL) elements, 30 all represented by freely moving agents. The adventitia is not modelled explicitly. The initial thickness of the tunica intima is assumed to be negligible, and the endothelial cells are not modelled explicitly. However, an implicit representation of the regenerated endothelium as an attri-bute of the SMCs (see below) is an important part of the model. 31 As the only function of the endothelium in this model is to inhibit SMC growth, the lack of endothelium cover on the IEL agents does not affect the results. Square stent struts, also consisting of freely moving agents, are located within the lumen, at some distance from the wall.
The stent deployment process gradually moves the struts outwards to the desired deployment depth, deforming the wall through force-based interaction between the agents. Touching agents are subject to an adhesive Hertzian contact force, while non-touching but nearby agents are attracted by a Hookean force representing the extracellular matrix. Cells at the left and right edges of the domain are constrained horizontally as a boundary condition. The strut agents are fixed in place at a stepwise increasing distance, and forces are equilibrated, forcing the wall to conform. If an agent sustains excessive mechanical stress, or also strain for IEL agents, it is removed. The endothelium is assumed to be completely removed by the angioplasty and stent deployment procedure. While we have not validated this 2D model against known mechanical properties of the arterial wall, our 3D version of the cell interaction model has been validated, see Ref. 14. Figure 1b shows the resulting initial conditions for the second phase of the simulation. Stent deployment has deformed the arterial wall, and IELs have been removed due to excessive hoop strain (near the strut) as well as longitudinal strain (around x = 0.4 mm), exposing the SMCs to the blood flow. No SMCs are covered by functional endothelium at this point.
During the second phase, the exposed SMCs have changed from contractile to synthetic, and traverse the cell cycle, proliferating into the lumen. The neointima is modelled as consisting mainly of SMCs. At each step of the model, each endothelium-free SMC becomes covered by functional endothelium with a time-varying probability calculated to achieve a certain coverage scenario. 29 Figure 1c shows the model state after some growth. In this example, to the left of the top strut, the endothelium has recovered quickly and proliferation has nearly stopped. On the right and at the bottom, proliferation is ongoing. Green cells are in the cell growth and proliferation cycle, with dark green cells having endothelium cover. Blue cells are endotheliumcovered and have been inhibited by nitric oxide concentration as a result of high wall shear stress (see below). Orange cells are quiescent.
After each step of the cell model, a fixed latticebased representation of the domain is constructed, with agent-covered nodes marked as solid. A constant parabolic velocity profile is set on the inlet, and a Lattice Boltzmann solver is used to compute the wall shear stress (WSS) at each SMC adjacent to the flow under the assumption that the total flux through the arterial segment remains constant as the restenosis develops. At locations where an intact and functional endothelium is present, nitric oxide (NO) is produced, which in turn inhibits SMC growth if its concentration is high enough. 31 The drug diffusion submodel, which operates on the same lattice, simulates antiproliferative drugs eluting from the stent (if any) and diffusing through the tissue, inhibiting proliferation. 5 SMC growth is also inhibited if the SMC is surrounded by other cells.
Following one of the approaches in experimental studies, we base our assessment of restenosis on the remaining cross-sectional lumen area. 24 We consider a restenosis to have taken place if the final cross-sectional area of the neointima is more than 50% of the original cross-sectional lumen area of the vessel. To estimate this area at a given simulation time, we measure the mean width of the lumen at that time, and estimate the cross-sectional area of the lumen as that of a circle with that diameter. We then subtract this from the initial post-stenting lumen area, to obtain the area of the neointima.

Model Set-up
Our simulation set-up is essentially the same as in the work of Tahir et al. 31 We configured the model to represent a small section of a healthy porcine coronary artery, with a length of 1.5 mm, a lumen diameter of 1 mm, and a tunica media consisting of five layers of SMCs. A healthy porcine artery was selected to better facilitate the model's validation, since the histological data points from the corresponding in vivo experiments are readily available. We used a bare metal stent with a mean deployment depth of 110 lm. The blood flow model was set to dynamic viscosity l = 4 mPa s, density q = 1000 kg m À3 and Re = 120, resulting in a mean steady flow velocity of 0.48 m s À1 .
We considered two re-endothelialization scenarios (Fig. 2). In the first scenario, the endothelium recovers to a coverage of 59% after 3 days, followed by a full recovery after 15 to 23 days (determined by an uncertain input parameter, see Table 1). In the second scenario, the initial fast recovery is absent, with endothelium cover increasing linearly from zero to 100% after 15 to 23 days. These scenarios correspond to the three endothelium recovery cases of Tahir et al. 31 The uncertainty range in Scenario 1 was chosen to cover the space between their cases 1 and 2.
In the model, reendothelialization is driven by three parameters: the initial fast recovery time, initial fast recovery degree, and the time to total recovery. Rather than having two scenarios, we could of course have taken the initial fast recovery time and degree as additional uncertain parameters, which would be appropriate as the measurement they are based on has a high relative error. 33 However, reducing the number of uncertain parameters reduces computational cost, and with the present set-up, we can demonstrate a comparison between scenarios in the presence of uncertainty.

Aleatory Uncertainty in the Model
The model is stochastic because the length of the cell cycle is chosen randomly for each SMC (normal distribution, l ¼ 32 h, r ¼ 2 h), the relative orientation of daughter cells is chosen randomly during mitosis, and the pattern of reendothelialization is random as well. For a detailed description of the model, we refer to previous publications. 5,31 Epistemic Uncertainty in the Model For our uncertainty quantification and sensitivity analysis of the ISR2D model, we consider the uncertainty in three input parameters (Table 1). Stent deployment depth and endothelium regeneration time were shown to strongly affect simulated neointimal area by Tahir et al., 29 and were therefore included. Additionally, we included blood flow velocity, because its effect on the behavior of the model has not yet been evaluated, and because of the potentially interesting interaction with the endothelium regeneration time. Blood flow in the coronaries is also variable depending on oxygen demand in the heart muscle, and by treating the blood flow as an uncertain parameter, we capture this effect into our analysis.
Flow velocity was varied by 10%. The other ranges were chosen to correspond to previously published results, 31 so that we could validate our findings.

UNCERTAINTY QUANTIFICATION
Here we provide some relevant details of the Uncertainty Quantification methods we applied. We describe our uncertainty measures of the total model response uncertainty, estimation of aleatory and epistemic uncertainty in the model result, and sampling scheme. In this section, the ISR2D model, described in ''The ISR Model'', is denoted by function fðn; xÞ, which depends on a vector of stochastic model parameters n described in ''Aleatory Uncertainty in the Model'' and uncertain model inputs x from Table 1 (Section ''Epistemic Uncertainty in the Model'').

Uncertainty Measures
In this study, we used the variance (Varðfðn; xÞÞ), standard deviation ( ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Varðfðn; xÞÞ p ), and the coefficient of variation (CV) as measures of model response uncertainty: where Eðfðn; xÞÞ is the mean value of the model response. These measures were estimated using a Monte Carlo method (Fig. 3) by running the model with different values of x and n, and collecting samples of the model results, in order to estimate the result uncertainty.

Aleatory Uncertainty Estimation
Aleatory uncertainty is the type of uncertainty arising from natural variability of the model output. 19 We collectively denote the stochastic variables that describe such variability in the model, as n. We can apply Saltelli's method 22 in order to measure the part of variance (Var T n ), which would remain even if all uncertainty in the model inputs is removed: where fðx j ; n j Þ and fðx j ; n jþM Þ are the model results with the uncertain inputs having the same values x j , but stochastic parameters having different values n j and n jþM . M is the number of evaluations of the difference in 2.
Additionally, we are interested in the partial standard deviation (r T n ), since it has the same units as the mean value. In general, we can approximate this value using a brute force Monte Carlo approach, but this is computationally expensive. However, Jensen's inequality 11 provides us a way to control from above the partial standard deviation with the square root of the partial variance Var T n : Epistemic Uncertainty Estimation Epistemic uncertainty is model imprecision, which arises due to lack of knowledge. 19 In this study, uncer- tain values of model inputs x generate such epistemic uncertainty. In order to estimate the parts of uncertainty arising due to imprecision in the value of each of the model inputs, we applied the Sobol variance-based sensitivity analysis method. 26 The Sobol Sensitivity Indices (SIs) are the ratio between the partial variance and the variance of the model function fðn; xÞ. In this work, we compute the first-order (S x i ) and total (S T n;x i ) SI for each of the uncertain model inputs, which are where Var x i and Var T n;x i are the partial variances, and n is the number of uncertain inputs. These partial variances for a parameter x i can be expressed by Refs. 22 and 23 Var x i ¼ Var x i ðE n;x Ài ðfðn; xÞjx i ÞÞ; which is the expected reduction in variance, if x i were known precisely, and Var T n;x i ¼ E x Ài ðVar n;x i ðfðn; xÞjx Ài ÞÞ; which is the expected remaining variance, if all the parameters except x i were known exactly. 1 x Ài denotes a vector of all inputs except x i . We can approximate the partial variances for the first-order and total SI for the ith parameter by Ref. 27 Var T n; where fðn jþM ; x jþM Ài ; x j i Þ is the result with input x i having the same value as for fðn j ; x j Þ, but the rest of the model parameters having different values. fðn jþM ; x j Ài ; x jþM i Þ denotes the model outputs with the same values of all inputs as for fðn j ; x j Þ, except for the values of n and x i .
The total number of the model evaluations required for computing both aleatory and epistemic uncertainties is Mð2n þ 2Þ, which is equal to N in this work.

Sampling
Instead of uniform random numbers, we used the Sobol sequence. 25 This sequence is called quasirandom, since it is not random, but preserves enough properties of random numbers to be used in Monte Carlo methods. The Sobol sequence is low-discrepancy, which means that it covers input space X more evenly.
For quasirandom sequences, convergence is of the Var(f (ξ, x)) (for 1 j n with n number of uncertain inputs) take a random value according to their distribution pðx j Þ. The model is run, and the stochastic parameters n take some random values during the simulation. Then, the model produces a value of the output f ðn ðiÞ ; x ðiÞ Þ. We run model in parallel N times. Using these N samples, we are able to estimate the uncertainty measures of the model result. 1 Note that the total sensitivity indices includes the contributions from the stochastic variable, since aleatory uncertainty is irreducible, and remains when values of other inputs are known precisely. 2 Here, the mean value, variances and standard deviation. For the sensitivity indices and partial variances, instead of sample size N, we test the number M, which is equal to the number of estimations of the expressions inside the sums in Eqs. (2), (5) and (6). evaluations of the estimator, we compute the estimators' mean ( I) and standard deviation (eð IÞ): Thus, the precision of the estimator obtained with N samples is I AE eð IÞ: Reliable values of I and eð IÞ can be obtained with the bootstrap when K is a large number. We chose K equal to 10,000.
In our experiment, we first set the sample size N to 480. The bootstrap test of the confidence in the mean value and standard deviation estimators showed that the estimators' standard deviations (eð IÞ) were up to 10%. We doubled the number of samples to N ¼ 960, resulting in standard deviations not exceeding 2.3%. This is deemed precise enough, and all results are based on using 960 samples.
We ran the ISR2D model once for each sample, for a total of 960 runs. The model is configured via a configuration file that contains a description of the model structure as well as values for its input parameters. We used a custom Python script to generate these configuration files, using the SOBOL library 21 to generate the Sobol quasirandom sequences. 28 As the runs are independent, they were run in parallel, using a single core for each run for maximum efficiency. Runs took on average 2 h and 15 min to complete, for a total of 4300 core hours. Approximately 100 more core hours were spent on postprocessing the results.

RESULTS
We observe a single output parameter of the model, the cross-sectional area of the neointimal growth, as a function of time. We compare the two endothelial recovery scenarios, taking into account the uncertainty of the model output, and estimate the contributions of the aleatory and epistemic uncertainty. We then investigate the relative importance of the input parameters as well as the inherent model stochasticity, in a sensitivity analysis. Finally, we show the spatial distribution of (explained) uncertainty. Figure 4 shows the probability of SMC presence over the domain at the end of the simulation run, giving an idea of the shape and size of the final neointimal growth produced by the model. Colors show the fraction of samples in which the given location was covered by neointima. On the whole, the results are similar for the two scenarios, an ovalshaped neointima surrounding the stent strut, shifted somewhat in the direction of the blood flow. As expected, the slower endothelium recovery in Scenario 2 leads to larger overall growth. Since this is a longitudinal view, the area difference is distorted; see below for a quantitative analysis. Note that in Scenario 2, it appears that the vessel can become fully closed. In fact, the top and bottom halves of the neointima sometimes form asymmetrically, and these different asymmetric runs are responsible for the innermost portions of neointima. In reality, the wall shear stress-induced growth inhibition will keep the lumen from closing up entirely in individual instances of the model.

Uncertainty Propagation
Figures 5(i) show the neointimal area over time as well as its estimated uncertainty. Initially, growth is slow as only a few SMCs are active, but as the size of the neointima increases, so does its growth rate (black line). Growth slows down again as the endothelium recovers and SMCs are increasingly inhibited, until a final state is reached. As expected, Scenario 2 leads to higher growth than Scenario 1 (% 0.38 mm 2 and % 0.28 mm 2 at the final time step). For the 1 mm vessel we simulated, a 50% cross-sectional area reduction corresponds to a neointimal area of approximately 0.39 mm 2 (horizontal dashed line). For both scenarios, the mean neointimal area is below this restenosis threshold. Overall uncertainty is fairly high however, with a standard deviation of 0.07 mm 2 for Scenario 1, and 0.09 mm 2 for Scenario 2 (red dashed lines). There are two sources of this uncertainty: the three uncertain input parameters give rise to epistemic uncertainty, while the internal stochastic variables result in aleatory uncertainty. Our results show that if the uncertain parameters were known exactly, the remaining aleatory standard deviation in the final state would be at most % 0.04 mm 2 in both scenarios (purple dashed lines). With respect to the occurrence of restenosis, we can conclude that in Scenario 1 a restenosis occurs in less than 10% of cases, while in Scenario 2 the probability is close to 50%.
The relative uncertainty (coefficient of variation, Eq. (1)) is depicted in Fig. 5(ii). Total relative uncertainty (red lines) is initially around 30%, then drops to % 25% as the growth slows and the slower growing runs catch up to some extent. The serrations on the left are due to all SMCs starting their cell cycle simultaneously; they later smooth out due to the variability in cycle length.
As stent deployment is only affected by deployment depth and not by any internal stochasticity, relative aleatory uncertainty (blue lines) starts out at zero. The aleatory uncertainty increases as SMCs start to proliferate and the stochastic variables come into effect, then drops slightly as the total uncertainty does to settle at % 14% and % 9% respectively. The final absolute aleatory uncertainty is actually nearly the same for both scenarios, but as the total area is larger in Scenario 2, the relative uncertainty is lower.
The probability density functions for the final neointimal area (Fig. 5(iii)) show nearly normal shape, validating our choice to use the variance as a measure of uncertainty.

Sensitivity Analysis
The sensitivity analysis results over simulation time are shown in Fig. 6. For each of the uncertain inputs, we computed the first order Sensitivity Index S x i , as well as the total sensitivity index S T n;x i , the latter including the combined stochastic model variables n. Shaded bands indicate a one standard deviation uncertainty interval of the estimate. Although at the number of samples we took there is a fairly large amount of uncertainty in the estimates, the relative importance of the input parameters is clear. Figure 6(i) show the variance in the neointimal area as a function of time. As we saw before, the total variance (black) is larger for Scenario 2 than for Scenario 1. The total (absolute) aleatory uncertainty (blue) is the same however. This measure includes the aleatory uncertainty by itself, and all combined effects of the aleatory uncertainty and the uncertain inputs. If all uncertain inputs were known exactly, this would be the remaining variance of the neointimal area.
For the final state of the simulation, the endothelium recovery time (red) is the most important parameter in both scenarios, followed by stent deployment depth (cyan) and flow velocity (green). This order is different in the beginning of the simulation, but this is difficult to see as the overall variance is close to zero. Figure 6(ii) show the first order sensitivity indices of the input parameters. As these are relative to the total variance, the relative importance over time can be studied more easily, although there is some noise for the early stages. Still, for both scenarios, it is now clear that the deployment depth explains almost all variance at the start of the growth process, while regeneration time becomes more important during the later stages. This makes sense, as the deployment depth affects the initial activation of the SMCs, while the regeneration time influences their inhibition. The importance of the flow velocity remains low throughout, which would suggest a low sensitivity to physiological variability of the blood flow. The importance of the regeneration time relative to the other parameters appears to be larger for Scenario 2; this is mostly due to the smaller relative share of the aleatory uncertainty however.
The black line in these graphs is the sum of the three first-order sensitivity indices, plus the total aleatory uncertainty. It thus represents all the sources of uncertainty except for the higher order combined effects of the input parameters. This shows that these combined effects explain at most in the order of 20% of the uncertainty.
In Fig. 6(iii), the total sensitivity indices (TSIs) for the input parameters are shown, which include all first-order and combined effects, including combined effects with the stochastic variables. Here we see the same pattern, with deployment depth dominating in the early stages, and regeneration time in the late stages of the simulation. A bump in the aleatory TSI is visible at days 10-15, and since this is included in the other plotted values, shows up in them as well. This is likely due to the stochastic variation in the SMC cell cycle lengths, which results in some simulation runs growing faster than others. Near the end of the growth phase, the slower runs catch up, and the variance decreases again. For Scenario 2, the flow velocity now shows a clear effect. This must be a combined effect, presumably with the regeneration time, as the timing is correct and the mechanism clear. For Scenario 1, flow velocity only shows a small combined effect at around day 7. Uncertainty   Figures 7 and 8 show the variance and SI in space for the final time step of the simulation for Scenarios 1 and 2, respectively. We can conclude that none of the uncertain inputs bring high uncertainty in a particular area of the observed space solely. In the results for the first scenario, we observe that the total SI for the deployment depth shows especially high value on the outflow side of the SMC growth. Moreover, we see that the SI of the deployment depth in average over space shows a higher value in comparison of the SI for other parameters.

DISCUSSION AND CONCLUSIONS
We performed uncertainty quantification and sensitivity analysis for a 2D in-stent restenosis model, demonstrating how the variation of a few selected parameters of the model affects the simulation results as well as whether knowing the exact value of a parameter is important for the precision of the simulation results.
From the sensitivity analysis we find that the largest contribution to the uncertainty comes from the endothelium regeneration speed, followed by the deployment depth. This is similar to results obtained in an earlier studies for a similar 3D model. 34 There it was also demonstrated that a difference in re-endothelialization speed changes growth to a large extent, while the effects of the deployment depth are smaller, but still quite pronounced. As the re-endothelialization has such significant impact on the neointimal area, this part of the ISR model requires further study and validation. We are currently undertaking more detailed modelling, based on controlled experiments of endothelial cell migration on substrates. The results once more confirm our hypothesis 31 that development of a restenosis is very much driven by the inability to quickly regenerate a functional endothelium.
However, these results cannot be considered final. The model parameters were based on a previous publication, 31 which considered a largely simplified model geometry, and some other parameters also were assumed without a thorough investigation on their distribution or variability range. For instance, the effect of the flow velocity requires a deeper study, looking in more detail at the actual physiological variability that may be expected (instead of the 10% variability that we now assumed). Additionally, it is hard to determine physiological values and draw conclusions about real systems for the 2D model, since the model considers a very simplified short and straight segment of the vessel. In real arteries, the curvature of the artery also plays a big role, causing uneven flow and also causing side effects such as hinging during the stent deployment. 12 Also, the model calculates a steady flow in the vessel, while in reality pulsatile flow results in a time-varying WSS and sometimes even WSS reversal. 32 Still, the results do show that once we have a reasonable understanding of the model and its parameters, we can extract important observations from the model (i.e. would a certain stent deployment lead to an actual restenosis) with a well defined uncertainty, even though this is a relatively involved multiscale model showing a complex biomechanical response. This method for uncertainty quantification and sensitivity analysis will therefore contribute strongly to the validation of these kinds of models.
On the other hand, we need to acknowledge the limited validity of the 2D model we used. We should therefore apply the methods described in this paper to the 3D version of the ISR model, 34 which is closer to the real system and also better validated. However, the 3D model contains millions of cell agents compared to a few thousand agents in the 2D model, and the 3D flow calculations are much more expensive as well. This makes the computational cost of Uncertainty Quantification and Sensitivity Analysis for the 3D model much higher or even prohibitive. We have therefore started to develop UQ algorithms for multiscale models that are capable of reliable estimation of the uncertainties while reducing the computational costs. 17 We are planning to test these algorithms with the ISR3D model.
We conclude that quasi-Monte Carlo UQ and Sobol sensitivity analysis are reliable methods for estimating uncertainties in the response of complicated multiscale cardiovascular models such as the ISR. They can be used to determine which parameters affect the model results the most, and are therefore the most important ones to obtain precise measurements on. Given measurements of the parameters, output uncertainty can then be assessed. We believe that using these methods to assess the quality and usability of a simulation model is a crucial step on the way to certification and use in for instance in-silico clinical trials for coronary stenting.

ACKNOWLEDGMENTS
The simulation was developed and tested on the Distributed ASCI Supercomputer (DAS-5). 2 The authors would like to thank the journal editor and the anonymous reviewers for their valuable suggestions and admonitions.

CONFLICT OF INTEREST
All authors declare that they have no conflict of interest.

ETHICAL APPROVAL
This article does not contain any studies with human participants or animals performed by any of the authors.

OPEN ACCESS
This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.