Global sensitivity analysis based on Gaussian-process metamodelling for complex biomechanical problems

Biomechanical models often need to describe very complex systems, organs or diseases, and hence also include a large number of parameters. One of the attractive features of physics-based models is that in those models (most) parameters have a clear physical meaning. Nevertheless, the determination of these parameters is often very elaborate and costly and shows a large scatter within the population. Hence, it is essential to identify the most important parameter for a particular problem at hand. In order to distinguish parameters which have a significant influence on a specific model output from non-influential parameters, we use sensitivity analysis, in particular the Sobol method as a global variance-based method. However, the Sobol method requires a large number of model evaluations, which is prohibitive for computationally expensive models. We therefore employ Gaussian processes as a metamodel for the underlying full model. Metamodelling introduces further uncertainty, which we also quantify. We demonstrate the approach by applying it to two different problems: nanoparticle-mediated drug delivery in a multiphase tumour-growth model, and arterial growth and remodelling. Even relatively small numbers of evaluations of the full model suffice to identify the influential parameters in both cases and to separate them from non-influential parameters. The approach also allows the quantification of higher-order interaction effects. We thus show that a variance-based global sensitivity analysis is feasible for computationally expensive biomechanical models. Different aspects of sensitivity analysis are covered including a transparent declaration of the uncertainties involved in the estimation process. Such a global sensitivity analysis not only helps to massively reduce costs for experimental determination of parameters but is also highly beneficial for inverse analysis of such complex models.


Introduction
Over the past few decades, computational biomechanical models have become an essential tool in research.The goal of these models is to allow predictions so as to better understand the underlying biological system or to support decision-making in a medical context, e.g., to choose the most efficient therapy for a specific patient.Nevertheless, the output of such models is inherently subject to uncertainty for various reasons: first, the underlying biological process is stochastic-which is particularly true for oncophysics and cancer treatment models (for example, branching process models for cancer [1,2], or stochastic models for immunotherapy of cancer [3]).Second, the experimental data used to calibrate models is uncertain [4].Third, the computational model itself includes sources of uncertainty, including the assumptions made to set up the model, other simplifications, or the input parameters [5].
When analysing the uncertainty of the model output, we distinguish uncertainty analysis from sensitivity analysis [5]: uncertainty analysis quantifies the uncertainty in the model output by propagating input uncertainties, via the model, onto the output [6, p. 262].Sensitivity analysis, on the other hand, apportions the uncertainty in the model output to different sources of uncertainty in the model input [7, p. 45].Inputs of interest can generally include not only model parameters but also boundary and initial conditions, assumptions, and constraints [8].In the context of sensitivity analysis, those inputs of interest are commonly referred to as factors.Here, we only consider model parameters as sources of uncertainty and refer to those as input parameters.
In this study, the goal is threefold: • Identify the most influential parameters on which further experimental estimation should focus (called factor prioritisation or ranking).• Identify parameters with little or no effect, which can thus be set to fixed values within their range (called factor fixing or screening).• Identify and quantify the interaction between parameters.This knowledge expedites the efficient design of future computational and experimental studies, while avoiding wasting resources on determining non-influential parameters.
We propose to apply a special type of sensitivity analysis to achieve the goals just described.One way of quantifying the sensitivity of the model output Y on the input parameter X i is to calculate the partial derivative ∂Y /∂X i .In practice, this involves choosing a base point X * and then perturbing one factor at a time while keeping all remaining factors fixed.This results in a local sensitivity measure at the base point X * which only explores one point of the input space and thus results in a deficient sensitivity analysis [5].In contrast, the Elementary Effects method (also called Morris method [9]) is not limited to one single point but explores the whole input space.It thereby overcomes the major limitation of local methods, while only requiring a relatively small number of model evaluations.While the Elementary Effects method is a global sensitivity analysis method, it only provides semi-quantitative information and is typically used for factor fixing [7,6].However, it cannot detect and quantify interactions between parameters and nonlinearities [10].In this work, we focus on complex biomechanical problems in which interactions between the different parameters can be expected.We therefore need a global method that can provide more detailed information.
Our method of choice is the Sobol 1 method [11,12], which is a variance-based global sensitivity analysis method that decomposes the output variance into portions attributed to the input parameters (see Fig. 1).The downside is that it requires many model evaluations, which quickly becomes computationally prohibitive in the case of complex models.We propose to introduce Gaussian processes [13] as a metamodel for the full model to mitigate the problem of computationally expensive model evaluations.Since the use of a metamodel introduces a further source of uncertainty in the sensitivity analysis, we estimate the uncertainty following the approach presented by Le Gratiet et

3/44
Uncertain input parameters X p(X 1 ) Training data

Figure 1. Overview of sensitivity analysis with a Gaussian process metamodel.
The Gaussian-process metamodel f GP,N (X) is trained on a data set generated by the full model f (X).Uncertain input parameters X i are propagated via the metamodel and result in an uncertain output Y .The Sobol method, as an example of a global sensitivity analysis tool, decomposes the output variance into portions attributed to the input parameters.
al. [14].After the full biomechanical model is substituted by the metamodel, we can calculate the Sobol indices based on Monte-Carlo integration.The uncertainty related to metamodelling and the uncertainty related to Monte-Carlo integration are analysed both separately and in total.
So far, the approach suggested by [14] has been applied to different computational models: an individual-based model simulation of microbial communities [15], a mathematical model of renal fibrosis [16], climate change simulations [17], and numerical wind-turbine models [18].Moreover, [19] applied the same idea to the calculation of Elementary Effects of a heart model.In [15,16,19,17], the authors state that they used the method, but no analysis of the associated uncertainties was presented.Only [18] quantified the uncertainty related to the metamodel and the uncertainty related to Monte-Carlo integration separately.
Our goal is to present the complete workflow of estimating Sobol indices based on Gaussian processes as a metamodel including the uncertainties: we demonstrate how to apply the approach suggested by [14] to two different biomechanical models, and we also assess the performance of the method when applied to such complex examples.
The article is structured as follows: we first introduce the reader to the Sobol method and to Gaussian processes in general.The Gaussian-process metamodel is then used to estimate the Sobol indices, including separate estimates for uncertainty related to In its most general form, a model f is a functional representation of the relevant physical process.The model calculates an output y = f (x) ∈ R for any given realisation x ∈ R D of the uncertain input parameters X, with D being the number of parameters (see upper left side of Fig. 1).We assume that the random vector X summarises the input parameters X 1 , . . ., X D , which are independent random variables.The probability distribution of X is described by the probability density function p(x).In the following, we assume that the input parameters are uniformly distributed, with no loss of generality.Common alternatives are the normal distribution or the log-normal distribution, among many others.
Our goal is to investigate the sensitivity of the model output to the uncertain input parameters X.Because of the randomness in the input parameters, the model output Y is also a random variable defined as (see upper right side of Fig. 1).The output distribution can (partially) be described by its first two moments: the expected value E [Y ], and the variance σ 2 Y = V [Y ] (with E [•] denoting the expectation operator, V [•] the variance operator, and σ Y the standard deviation as given in Supplement A.1).
One way of characterising sensitivity is to decompose the variance of the output V [Y ] into portions ascribed to the individual input parameters.A common sensitivity analysis method based on the decomposition of variance is the Sobol method [26].The core idea is to decompose the output variance V [Y ] as with the conditional variances given by where X ∼i denotes the vector of all input parameters except X i .Thus, the output variance V [Y ] is the sum of variances contributed by input parameter X i , including interactions with other parameters.The idea now is to attribute the total variance to the individual input parameters according to their variance contribution.Note that V ij is the variance contributed by the input parameters X i and X j but not expressed 5/44 in V i nor V j .This is called the interaction of parameters X i and X j .Note that this decomposition assumes statistical independence of the input parameters is the portion of the output variance ascribed to input parameter X i , we define the first-order Sobol index S i as The numerator describes the extent to which the output variance V [Y ] would be reduced if the parameter X i was fixed.A parameter X i with a high first-order index S i should, hence, have priority when determining parameters based on experiments so as to efficiently reduce the overall uncertainty of the model.The first-order Sobol index is typically used to identify the most influential parameters, which is our first goal [10].Moreover, a parameter with a high first-order index S i is more likely to be identifiable from experiments, but can still be non-identifiable [8]: to decide whether a parameter is identifiable or not, an identifiability analysis is required, which complements the sensitivity analysis (see [27] for an overview of identifiability analysis).
The question now arises as to whether S i = 0 is also sufficient to conclude that a parameter has no influence.In fact, this is not the case because the parameter might be involved in interactions with other parameters.A parameter may have no effect if it is varied alone; however, this may be different when it is varied in combination with another parameter, or even with several other parameters.An additional sensitivity measure that includes higher-order interaction effects is needed to identify non-influential parameters.The total-order Sobol index is, therefore, defined as In this case, the numerator ] describes the expected output variance that would be left if all parameters but X i were to be determined [10].If-and only if-this expected output variance is close to zero, is the parameter X i non-influential.The total-order index describes the total contribution of the parameter X i to the output Y : this includes the first-order effect plus any higher-order effects that arise from interactions.The difference S T i − S i , then, indicates interaction effects between factor X i and any other factor [6, p. 167].As mentioned above, the total-order index is particularly helpful in the context of factor fixing: if S T i = 0 (or is in practice sufficiently small), the parameter X i is non-influential and can be fixed anywhere in its input range without affecting the output variance.So, the first-order and the total-order Sobol indices serve our first two goals: identify the most influential and the non-influential parameters.To additionally identify interactions between two specific parameters X i and X j -which is our third goal-we define the second-order Sobol index as Finally, dividing Eq. ( 2) by V [Y ] and inserting Eqs.(3) and (5) leads to All sensitivity indices, thus, sum up to 1; furthermore, they are non-negative.This leads to an interesting implication which is worth noting: even when we have a large 6/44 number of parameters, we cannot have a large number of influential parameters.If all D parameters are equally influential, each can only contribute 1/D of the variance.If, however, a few parameters have a strong influence on the output Y , the remaining parameters can contribute even less (see [6,Sec. 2.4.3]).As [28] stated: only a small subset of parameters significantly influences one specific system output ("sparsity of factors" principle).It should also be noted that the total-order indices S T i do not, in general, sum up to 1.

Numerical approximation of Sobol indices
To estimate the Sobol indices according to Eqs. ( 3), ( 4) and ( 5), we need to compute conditional variances, e.g., , which involves evaluating multidimensional integrals in the space of the input parameters R D .Numerical integration based on quadrature rules becomes prohibitively expensive as the number of input space dimensions increases.This is why Monte-Carlo integration is employed, the accuracy of which is independent of the number of input space dimensions [29].For each single integral, Monte-Carlo integration involves evaluating M Monte-Carlo samples: to compute, for example, To make the estimation of Sobol indices more efficient, [30] rewrote the multidimensional integral so that it can be computed using a single Monte-Carlo loop (summarised in Supplement A.2).
To make the best use of the model evaluations, we employ the efficient algorithms suggested by [32]: to estimate the first and the total-order indices, we generate M samples row-wise concatenated as a matrix A and further M samples concatenated as a matrix B. This results in two independent M × D matrices.We introduce a third matrix A (i) B for each input space dimension i = 1, . . ., D, where all columns are taken from A except the i-th column, which is taken from B. One sample (B) m and the corresponding sample (A (i) B ) m have X i in common but differ in all other parameters X ∼i .
To calculate the first-order index, we then use the estimator proposed by Saltelli et al. [33] V and for the total-order index, we use the estimator proposed by Jansen [34]: Alternative forms were presented in [11,34,35,36], among others.The denominator is estimated as where we estimate the variance of the output as the sample variance V of evaluations of all samples A and B. This yields better results, i.e., an estimator with lower variance, compared to [32]. 2 The error of the Monte-Carlo estimate for the expectation is proportional to , with g denoting the integrand.If we assume V [g] to be fixed, we have to increase the number of Monte-Carlo samples M , and the error of the estimate thus decreases by 1 √ M [31].

7/44
In addition to the first and total-order indices, we estimate the second-order indices as proposed by Saltelli [32]: where the matrix B (i) More details on different sensitivity-index estimators can be found in [32,33], among others.
We hence have to evaluate our model f at all samples of the triplet A, B and A A ), resulting in M (2D + 2) simulations in total (for more details see Supplement A.3 and the original publication by [32]).In practice, Quasi-Monte-Carlo (QMC) integration is often used to generate the samples because of its superior rate of convergence compared to Monte-Carlo integration [37].

Gaussian process metamodels
As just described, Monte-Carlo integration to estimate the Sobol indices requires a large number of sample evaluations and is thus computationally prohibitive if the evaluation of the underlying model is expensive.We therefore use a metamodel (also known as surrogate model or emulator) as an approximation of the full model.Classically used metamodels include polynomials, splines, neural networks, polynomial chaos expansion, support vector regression, and Gaussian processes (GPs), among others [38,39].Before a metamodel can be used for a sensitivity analysis, for example, it has to be trained to later ensure that it is a good approximation of the full model.This process consists of three steps, which we first summarise (see Fig. 2) and, then, explain in more detail below: 1. Generate N training samples summarised in X (resulting in an N × D matrix).

Evaluate the full model at the training samples to obtain the corresponding
response: Y = f (X ) (resulting in an N × 1 vector).
3. Form and train the metamodel.
First, we generate N training samples that are summarised in the matrix X (N × D matrix).The choice of training samples has to provide a good coverage of the input space to later ensure a good predictive quality of the metamodel.To this end, we use a QMC approach based on Sobol sequences [40] to generate the training samples.

Remark (Sequential design).
A commonly used alternative to a QMC approach is Latin Hypercube Sampling (LHS) [41].[38] state that optimised LHS is particularly well-suited for metamodel fitting.However, sequential design is also important in the context of metamodelling: if the original number of training samples is not sufficient to achieve a good predictive quality of the metamodel, additional samples can be added while still making use of the original training samples.Since only advanced LHS methods [42,43] enable sequentially adding new points, while this is straightforward with QMC schemes [44], we use a QMC approach. 8/44 where each row corresponds to one training point.Finally, we have to train the metamodel.We use a GP f GP as metamodel similar to [45,46,14] and summarise what training means for GPs below.Note that we only include a compact overview of the most relevant concepts used in this paper.For more details, the reader is referred e.g. to [13].
A GP defines a distribution over functions such that any finite set of function values f GP (x (1) ), f GP (x (2) ), . . ., f GP (x (n) ) has a joint Gaussian distribution [13,47].From a Bayesian point of view, we distinguish between prior and posterior: the prior GP reflects our beliefs about the metamodel before seeing any (training) data, and the posterior GP is then conditioned on the (training) data, i.e. includes the knowledge from the data (see Fig. 2).This conditioning on the data is what we refer to as training.
The prior GP f GP (X) is given by and is completely specified by its mean function m GP (X) and its covariance function k(X, X ) between all possible pairs (X, X ).The covariance function is a positive definite kernel, e.g., the squared exponential covariance function (also called radial basis function) with the characteristic length scale , variance parameter σ f and || • || denoting the Euclidean L2-norm.We assume that the prior mean function is zero: m GP (X) = 0, which is common practice and does not limit the GP model, as any uncertainty about the mean function can be included in the choice of a covariance function [47].Different covariance functions exist and can be combined, for example, through multiplication or addition.[47, Chap.2] presents a concise overview of different covariance functions for GPs and how to use them to express the structure of the data.The choice of a suitable covariance function is essential since the more a-priori knowledge goes into choosing the covariance function, the fewer data we need to train the metamodel [48].
As described above, we observe the output at N training points.The goal, then, is to predict the output at N * new points summarised in X * ; in our case, those new points (where we predict the output) will be the Monte-Carlo samples for the estimation of the Sobol indices.
Remember that Eq. ( 12) is only the prior distribution and does not yet incorporate our knowledge from the training data.To obtain the posterior, we now condition the prior GP on our set of N training data points.This conditioning results in the key predictive equations In Eqs. ( 15) and ( 16), K ε = K + σ 2 y I, where K = k(X , X ) denotes the N × N matrix we obtain when evaluating the covariance function (given in Eq. ( 13)) for all pairs of training points, similarly for K * = k(X , X * ) and K * * = k(X * , X * ).We use σ 2 y as an artificially introduced variable nugget term to alleviate numerical problems [49,50].The 10/44 hyperparameters θ = (σ f , ) are optimised by maximising the log marginal likelihood using a gradient-based optimiser.The log marginal likelihood is given by log with | • | denoting the determinant.Maximising the log marginal likelihood given by Eq. ( 17) with respect to the hyperparameters θ automatically incorporates a trade-off between model fit and model complexity: the first term in Eq. ( 17) penalises the model's failure to describe the data while the second term penalises high model complexity.Thus, this favours the least complex model that is able to explain the data [13].
One advantage of employing GPs as a metamodel is that predictions can be computed exactly in a closed form [47] and that GPs inherently provide uncertainty measures over the predictions.Moreover, one can incorporate a wide range of modelling assumptions into the choice of the covariance function.However, note that computing the inverse in the first term in Eq. ( 17) (and the determinant in the second term) is computationally expensive, i.e., on the order O(N 3 ).This cubic complexity results in slow inference as the number of training samples increases.One further challenge of using GPs as a metamodel is that they are susceptible to the curse of dimensionality: as the dimensionality of the input space increases, the number of training samples required to train the metamodel grows exponentially [51,52] and the optimisation of hyperparameters θ becomes impractical.[53] review approaches to improve the scalability of GPs to large data sets, e.g., by using stochastic variational inference [54]; [52] present an approach with built-in dimensionality reduction.

Estimation of Sobol indices and their uncertainty
To estimate the Sobol indices, we now use the estimators given by Eqs. ( 7), ( 8) and (10) and substitute the realisations of the full model f with those of the trained GP metamodel f GP,N as suggested by [14]: with M again being the number of Monte-Carlo samples.We summarise the estimates as Ŝ with ∈ {i, T i, ij} for the first, total, or second-order index estimates, respectively.Remember that we now evaluate the Monte-Carlo samples with the metamodel instead of the full model.We can, therefore, afford considerably larger numbers of Monte-Carlo samples.Since we sample realisations of the GP metamodel f GP,N , the resulting estimates Ŝ are again random variables.These include two sources of uncertainty: one related to the metamodel approximation and one related to the Monte-Carlo integration.To estimate those uncertainties, and additionally the total uncertainty, we employ the algorithm suggested by [14].The steps described in the following can equally be applied to all indices of different order.

11/44
Algorithm 1. Calculation of (N GP ×B) estimates Ŝ k,b of the Sobol index.We sample N GP realisations of the GP metamodel and subsequently resample each realisation B times using the bootstrap technique as suggested by Le Gratiet et al. [14].

Draw samples with replacement:
We visually summarise the approach in Algorithm 1; a more detailed version is included in the Supplement A.4.The core idea is to sample N GP realisations of the GP metamodel and, subsequently, resample each realisation B times using the bootstrap technique [55].This results in N GP × B estimates Ŝ k,b of the respective Sobol index.We then calculate the mean as (21) and the total variance as Since this estimator includes two sources of uncertainty (one related to the metamodel approximation and one related to the Monte-Carlo integration), we decompose the variance of S as Sampling realisations of the metamodel f GP,N (X) as opposed to using only the predictive mean m GP,N (X) allows us to take into account the covariance structure of the metamodel.The part of the variance related to the metamodel approximation can 12/44 be estimated as Alternatively, [56] present an approach to estimate an upper bound for the metamodel error based directly on the covariance function of the GP, but their approach only provides a rough upper bound [14].In addition, [57] present an approach to investigate the accuracy of Sobol indices based on a general relation between the accuracy of an arbitrary metamodel and the error of the estimated indices.
To estimate the uncertainty related to Monte-Carlo integration, we use the bootstrap technique [55].The Monte-Carlo samples A, B, A

B and B (i)
A are resampled (i.e.sampled with replacement) B times as depicted in Algorithm 1.We then calculate Ŝ according to Eqs. ( 18), ( 19) and (20) for each bootstrap sample, resulting in B estimates for the Sobol index for each realisation k of the metamodel.The part of the variance related to the Monte-Carlo integration is given by The bootstrap technique is based on the fact that sampling with replacement from a set of independent, identically distributed data equals sampling from the empirical distribution function of the data [58].It is important to note that bootstrapping does not require further model evaluations.For a general introduction to the bootstrap technique, the reader is referred to [59] or [60, Sec.5.2].
3 Application to nanoparticle-mediated drug delivery in a multiphase tumour-growth model

Model definition
An excellent example for the proposed overall approach is nanoparticle-mediated drug delivery in a multiphase tumour-growth model as all challenging motivating arguments for our approach are present in this problem class, like complex costly models and a large number of parameters.The tumour-growth model in its original form is based on the works [61,62,63].[20,21] extended the model to a five-phase model including the vasculature.Finally, [22] included and studied nanoparticle delivery.We use the tumour-growth model as a precursor to generate physically plausible results of the tumour and its microenvironment.Those results then serve as initial condition for the sensitivity analysis of nanoparticle-mediated drug delivery.

Underlying tumour-growth model
The model considers the tumour as a porous structure: the extracellular matrix (ECM) is the solid phase (denoted by superscript s) with several fluid phases filling its pore space.We include three fluid phases: tumour cells, host cells, and interstitial fluid (IF), denoted by superscripts t, h and l, respectively.In addition, the vasculature is modelled as an independent porous network and denoted by superscript v.The governing equations of the model are formulated on the macroscale by employing the Thermodynamically Constrained Averaging Theory (TCAT) [64,65].Each phase is modelled in an averaged sense based on volume fractions ε α with α denoting an arbitrary phase.phases must satisfy the equation The different phases additionally transport species.The vasculature and the IF transport oxygen with mass fractions denoted by ω nv and ω n l, respectively.Similarly, ω NPv and ω NP l denote the mass fractions of nanoparticles in the vasculature and the IF, respectively.Finally, tumour cells and host cells are divided up into living and necrotic cells.The mass fraction of necrotic tumour cells is denoted by ω N t and the mass fraction of necrotic host cells by ω N h.Fig. 3 schematically summarises all of the components of our multiphase tumour-growth model that are considered here.Note that the lymph system is not explicitly modelled.Many commonly used tumour-growth models are data driven and thus adhere to observed data.Our model, by contrast, is based on physical laws.To give just one example: we use Fick's laws to describe the motion of oxygen in IF and Darcy's law to describe the flow of IF through the pores of the extracellular matrix.Such a physicsbased approach allows us to describe the system, even under unobserved circumstances.Nevertheless, [6, Sec.1.1.4]state that physics-based models are customarily overparametrised: they include more laws and parameters than available data would support.This becomes particularly critical when model parameters are to be determined, e.g. by inverse analysis, and thus sensitivity analysis becomes a crucial part of model development [66]. 14/44

Nanoparticle-mediated drug delivery
The transport of nanoparticles in our multiphase tumour-growth model is included as described in [22].We only present a short summary in the following.Further details can be found in the original publication [22].
Nanoparticles are intravenously injected and transported in the vasculature.They then extravasate into the IF and are transported towards tumour cells and host cells.The focus here lies on the extravasation into the IF and the transport therein.Therefore, we do not explicitly include transport in the vasculature, but rather assume a constant mass fraction of nanoparticles in the vasculature ω NPv where the physical interpretation of the different transport mechanisms is included in Fig. 3. Nanoparticles extravasate from the vasculature into the IF through two different pathways: the interendothelial and the transendothelial pathway [67,68].
Those pathways are described by with the oncotic pressure difference between blood vessels and IF σ π v − π l , the surface-to-volume ration S/V , and Macaulay brackets • + .This equation is based on the Staverman-Kedem-Katchalsky equation similar to [69].
The first term describes the interendothelial pathway, which is a convective process: nanoparticles are dragged by the transvascular fluid flow through gaps in the blood-vessel wall [69].This process is governed by the hydraulic conductivity L v p of the blood-vessel wall which is defined as with the pore radius r 0 , the vessel-wall thickness t, and the fraction of pores γ p [70].
The second term in Eq. ( 28) describes the transendothelial pathway, which is a diffusive process: nanoparticles diffuse through the vessel-wall.This diffusive flux is 15/44 governed by the vascular permeability P v .
Finally, the lymph system absorbs nanoparticles from the IF NPl→NP ly governed by the lymphatic filtration coefficient L p S V ly .Lymphatic drainage is impaired above the collapsing pressure p ly coll , and thus no particles are removed.We further assume that the nanoparticles transport and release anti-cancer drugs.Those kill tumour cells and thus increase the mass fraction of necrotic tumour cells ω N t.At the same time, those drugs have adverse side effects and kill host cells.This additionally increases the mass fraction of necrotic host cells ω N h.For the sake of simplicity, we assume that the mass fraction of killed cells is directly proportional to the mass fraction of nanoparticles present in the IF at a certain position.We introduce intraphase reaction terms that increase the mass fraction of necrotic tumour cells and host cells according to where γ t kill and γ h kill characterise the strength of the drug.As quantity of interest for the sensitivity analysis, we consider the mean of the necrotic fraction of tumour cells given by where the tumour size is defined as A t = ´H(S t − 0.1) dΩ with the Heaviside function H(•), and S t denotes the saturation of tumour cells [20].We define the tumour as the part of the domain where S t > 0.1.In the context of sensitivity analysis, it is important to choose the quantity of interest carefully and to bear in mind that a parameter that is non-influential under one particular investigated condition, e.g., one particular quantity of interest, might be highly influential under a new condition [8].

Set-up of the numerical example
The set-up presented in the following is similar to the example in [22].We therefore only present a summary, and again refer the interested reader to the original publication for further details.The major addition to the original example is the nanoparticle-mediated killing of tumour and host cells as described by Eqs. ( 31) and (32).We investigate nanoparticle transport and subsequent killing of cells by nanoparticlemediated drugs.The transport of nanoparticles depends on the hydraulic conductivity of blood-vessel walls L v p and the blood-vessel wall permeability P v (both influence the transport of nanoparticles into the IF), the diffusivity D NPl of nanoparticles in the IF, and the lymphatic filtration coefficient L p S V ly (influencing the transport of nanoparticles out of the IF).Subsequently, drugs mediated by the nanoparticles kill tumour and host cells depending on the strength of the drug, characterised by γ t kill and γ h kill .Note that the amount of killed cells largely depends on the amount of nanoparticles reaching a particular region of the domain, and hence depends on the transport parameters.
The transport of nanoparticles, and thus the question of which regions nanoparticles reach and where drugs can kill cells, essentially depends on the microenvironment of the tumour.Solid tumours exhibit typical features relevant in this context: the majority of living tumour cells is located in the tumour periphery, whereas the tumour core mainly consists of necrotic cells (see Fig. 4A).In addition, Fig. 4B shows that the volume fraction of the vasculature is considerably lower in the tumour area because the growing tumour collapses blood-vessels.The inner core of the tumour even contains no vessels at all.Finally, the interstitial pressure in the tumour is increased and can reach 6 mmHg (see Fig. 4C).To sum up, the tumour has a necrotic core with collapsed blood vessels as well as an increased interstitial pressure, which is a structure typical for solid tumours and which has also been observed in experiments [71,72,73,74].We analyse a domain of 1 mm × 1 mm, but due to the symmetry of the problem, we only simulate one quarter of the domain (500 µm × 500 µm).The grown tumour has a radius of 440 µm as presented in Fig. 4. We analyse a time interval of 20 min of nanoparticle transport and killing of cells based on examples by [75,76].Assuming the intravenous infusion of nanoparticles, we prescribe a constant value of ω NPv = 2.0 × 10 −3 for the mass fraction of nanoparticles in the vasculature.
Table 1 summarises the six uncertain input parameters included in the sensitivity analysis.We assume that all input parameters are distributed uniformly within the 17/44 The given values for the hydraulic conductivity of the blood-vessel wall correspond to a pore radius of r 0 = [50; 200] nm as used in [22].
given ranges, which are based on experimental data (see references in Table 1).The uniform distribution is chosen because we lack more specific information about the input parameters: given only the range of the input parameters (and no further information such as mean or variance), uniform distributions maximise the information entropy and hence minimise the introduced bias [77,78].Note that the killing coefficient of host cells γ h kill has no influence on our quantity of interest, the mean of the necrotic fraction of tumour cells-neither directly nor indirectly through coupling terms.We nevertheless include the killing coefficient of host cells γ h kill in the sensitivity analysis to investigate how reliably we can identify a non-influential input parameter as such.
Fig. 4D and E present a result for the distribution of nanoparticles in the IF and for the mass fraction of necrotic tumour cells: for this example, we used the mean values of the six uncertain input parameters given in Table 1.Most nanoparticles accumulate at the edge of the tumour, while lymphatic drainage removes most particles outside the tumour area, and roughly 50% of the tumour cells are necrotic.
The nanoparticle-mediated transport included in the multiphase tumour-growth model is implemented in our in-house research code BACI [82].The sensitivity analysis methods, as presented above, are implemented in QUEENS [83].QUEENS is a general purpose framework for uncertainty quantification, physics-informed machine learning, Bayesian optimisation, inverse problems and simulation analytics on distributed computer systems.We use GPy [84] as a GP framework and PyTorch [85] to generate Sobol sequences.

Predictive quality of the metamodel
Since we use a GP metamodel to estimate the Sobol indices, we first assess its predictive quality.To this end, we investigate the quality of the metamodel predictions for two different covariance functions k(X, X ) used for the GP: we compare a tensorised, squared, exponential covariance function to a tensorised 5/2-Matérn covariance function (with ν = 5/2) [86,87,13].A tensorised covariance function has the form k(X, ) and as such includes a set of hyperparameters θ = (σ f i , i ) i=1,...,D for all input space dimensions, which we optimise by maximising the log marginal likelihood.
For this comparison, we consider different sizes of training sample sets N = [10,15,20,25,30,40,60,80,100,150,200], which we generate based on Sobol sequences.Additionally, we generate a set T of N T = 1000 testing samples disjoint of the training samples.Based on the training samples and the testing samples, we calculate the Nash-Sutcliffe efficiency 18/44 Q 2 [88] given by similar to [14].This is based on the predictive posterior mean m GP,N (X) of the GP with optimised hyperparameters, and thus compares the mean prediction of the posterior GP to the actual output of the full model f .A Nash-Sutcliffe efficiency close to 1 indicates good agreement and, hence, reliable predictions.Fig. 5 shows good convergence of the Nash-Sutcliffe efficiency for both covariance functions with values close to 1, even for smaller training sample set sizes.If the number of training samples is restricted due to the computational cost, [89] suggest an algorithm to improve the metamodelling accuracy and efficiency based on sequential sampling.Note that we use a set of testing samples here that is disjoint of our set of training samples; this means that we also evaluate our full model N T times, which might be infeasible if the model is computationally more expensive.In those cases, one can use cross-validation methods, such as those explained in [13,Sec. 5.3], where the training set itself is split into two disjoint sets: one is actually used for training and the other for validation.
Both covariance functions yield a very similar predictive quality.In the following, we only use the tensorised, squared, exponential covariance function because it is the default choice in most applications of GPs [47] and, moreover, because it is a universal covariance function [90].For the detail plot, we repeated the process five times with different training sample sets of sizes N = [10,15,20,25,30].For reference, the dashed lines in the detail plot are identical to the dashed lines in the main plot.

19/44
Remark (Randomness of metamodel training).In Fig. 5 we notice a small kink in the Nash-Sutcliffe efficiency for N = 25 in the case of the Matérn covariance function.Therefore, we repeat the training of the GP with other randomly chosen training sample sets for N = [10,15,20,25,30] to check whether the original training sets happen to perform exceptionally well.The grey detail plot in Fig. 5 presents the results: while the Nash-Sutcliffe efficiency is still above 0.90 in all cases, we notice that some training sample sets result in slightly worse efficiencies than others for those smaller sample sizes.This is precisely the case for the Matérn covariance function with N = 25 in the main plot.Moreover, the two different covariance functions lead to slightly different efficiencies.
One possible reason for such behaviour may be the convergence of the optimiser used to optimise the hyperparameters according to Eq. (17).As with all gradient-based optimisers, the optimisation may get stuck in a local minimum.To avoid ending the optimisation in a local minimum, one can repeat the optimisation multiple times from random initial points [52], as provided by the package GPy [84], for example, or use stochastic optimisation, such as Adam optimisation [91].
Further, we take a closer look at the underlying GP, which is depicted in Fig. 6 for 20 training samples.The projection of the GP into the input-space dimensions reveals a linear relation in most dimensions.We only see considerable nonlinearity for the blood-vessel wall permeability P v .Those characteristics make it much easier to train the GP based on a small number of training samples.

Remark (Projection of the D-dimensional Gaussian process).
The projection m GPi (X i ) of the D-dimensional posterior mean m GP,N (X) into the input space dimension X i (as presented in Fig. 6) is calculated as follows.First, we uniformly sample discrete values of the posterior mean m GP,N (X) in the D-dimensional input space.Second, we project those values over the input space dimension X i .Third, the results are binned in the X i -direction, and we calculate the mean and confidence interval for each bin.Note that we only project the posterior mean m GP,N (X) and neglect the covariance function k N here.
Plotting the model output over a specific input in the form of scatterplots-as done with the training samples in Fig. 6-helps us gain a general understanding of the magnitude of the underlying sensitivity [10].[6, Sec.1.2.7] offers a compelling interpretation of scatterplots in relation to the first-order Sobol index: if the conditional expectation E X∼i [Y |X i ]-here represented by the projection of the mean m GPi of the GP-has a large variation across X i , the corresponding input parameter has a high firstorder Sobol index.Fig. 6 reveals that the projection of the mean m GPi is almost constant for the diffusivity of nanoparticles D NPl , the lymphatic filtration coefficient L p S V ly , and the killing coefficient of host cells γ h kill .In contrast, the variation of the projection of the mean m GPi is larger for the blood-vessel wall L v p , the vascular permeability P v , and the killing coefficient γ t kill , and we therefore expect the output to be highly sensitive to those parameters.

First-order Sobol index estimates
We now assess the convergence of the first-order Sobol index estimates for increasing numbers of training samples.To calculate the mean Si based on Eq. ( 21), we use N GP = 500 metamodel realisations, and the number of Monte-Carlo samples is set to M = 10 000.We do not include bootstrapping here.The uncertainties in the estimates will be studied in the next section.

21/44
parameters-namely the vascular permeability P v , the killing coefficient γ t kill , and the hydraulic conductivity of the blood-vessel wall L v p -have considerably higher first-order Sobol indices than the remaining three parameters.Fig. 7 also allows to assess the convergence of the Sobol indices for an increasing number of training samples: even small sizes of training sample sets yield values close to the value based on N = 200.This is due to the high values of the respective Nash-Sutcliffe efficiency, as discussed in the previous section.Those results are promising, in particular for models that are computationally very expensive, and thus do not allow a large number of evaluations of the full model: the computational cost of the demonstrated approach is considerably reduced compared to an analysis based directly on evaluations of the full model, as for example presented in [92,93] for tumour-growth models.
We assessed convergence visually based on Fig. 7.In addition, [94] present a thorough definition of convergence criteria for global sensitivity analysis results.Nevertheless, the computationally limiting factor is usually the number of Monte-Carlo samples.Since we evaluate the Monte-Carlo samples on the GP metamodel, this limitation is less critical in the presented workflow.If sampling the realisations of the metamodel for very large numbers of Monte-Carlo samples becomes an issue, [14] include an efficient approach based on conditional GPs.

Uncertainties of Sobol index estimation
We now go on to not only estimate the mean Si but also include the uncertainty related to the metamodel and to the Monte-Carlo integration given by Eqs. ( 24) and (25).In addition to the first-order index S i , we also include the total-order Sobol index S T i .Again, we use different training sample set sizes N = [10, . . ., 200] and a tensorised, squared, exponential covariance function with hyperparameters optimised based on maximising the log marginal likelihood of the GP.We draw N GP = 500 realisations of the GP, B = 300 bootstrap samples, and M = 10 000 Monte-Carlo samples.We calculate 95% confidence intervals on the basis of the variance related to the metamodel σ2 GP and the variance related to Monte-Carlo integration σ2 MC .Fig. 8A presents the results for all six input parameters.Note the different scaling on the vertical axes.We first take a look at the results for the indices themselves.Fig. 8A again confirms that even for small numbers of training samples, the estimates for first and the total-order Sobol indices rapidly converge.As mentioned above, we do not expect any influence of the parameter γ h kill on the quantity of interest.Fig. 8A shows that we can identify this non-influential parameter as such even for small numbers of training samples.The parameter D NPl also leads to Sobol indices close to zero for N < 60.For larger training sample set sizes however we get a slightly higher total-order Sobol index, which is nevertheless small.Hence, we can clearly separate the three most influential parameters from the three non-influential parameters.
Moreover, the total-order index is higher than the first-order index, in particular for the hydraulic conductivity of the blood-vessel wall L v p and the vessel wall permeability P v .This leads to the conclusion that higher-order effects are indeed present.We will therefore analyse the second-order Sobol indices in the next section.In addition, Table 2 summarises the values for the first and the total-order indices for 200 training samples: the sum of all first-order Sobol indices is 0.967.Since this is close to one, we conclude that higher-order effects are present, but only play a minor role.The largest part of the output variance is covered by the first-order indices.
We now focus on the uncertainties: we assess the uncertainty related to the GP metamodel and the total uncertainty, where the latter includes both sources of uncertainty (related to Monte-Carlo integration and related to the metamodel).For small training sample set sizes, we see considerable uncertainty related to the metamodel (depicted in light blue/orange in Fig. 8A).However, the uncertainty related to Monte-Carlo 22/44

24/44
integration dominates for N > 40.We therefore present the uncertainty related to the metamodel in detail in Fig. 8B: the uncertainty rapidly decreases as the number of training samples increases for both the first-order and the total-order index and becomes one order of magnitude smaller than the uncertainty related to Monte-Carlo integration.[18] also found the uncertainty related to the GP metamodel to be much smaller than the uncertainty related to Monte-Carlo integration in their example.The total uncertainty (depicted in grey) could be reduced even further by increasing the number of Monte-Carlo samples.
Based on these results, we conclude that including the uncertainty related to the metamodel is not absolutely necessary in our example.However, the example presented in the outlook and the example presented by [14] illustrate that this is not always the case: only taking into account the Monte-Carlo uncertainty might then underestimate the confidence interval.In such cases, it is essential to consider the uncertainty related to the metamodel.Hence, this largely depends on the model, the input parameters, and the quantity of interest, and no one-size-fits-all rule can be given.
Nevertheless, even taking into account the metamodel and the Monte-Carlo uncertainty may incorrectly estimate the confidence intervals: poor optimisation of the hyperparameters may result in underestimated or overestimated confidence intervals.In such cases, one could additionally consider the uncertainty related to the estimation of the hyperparameters of the GP covariance function by using a full-Bayesian approach with hyperpriors [14].
To sum up, the demonstrated workflow not only identifies parameters with a high first-order Sobol index (necessary for factor prioritisation) but also parameters with a small total-order Sobol index (necessary for factor fixing).In both cases, small numbers of training samples suffice in our example.

Second-order Sobol index estimation
Since we concluded from the results in the previous sections that higher-order effects are indeed present in our example, the goal now is to estimate the second-order Sobol indices, and thereby identify interaction effects between the input parameters.The results in the previous section show that the uncertainty related to Monte-Carlo integration is dominant, and the uncertainty related to the GP metamodel is much smaller.We therefore estimate the second-order indices based on the predictive mean m GP,N (X) of the GP and do not take into account the uncertainty related to the metamodel.Thus, we estimate the second-order Sobol indices as which includes estimating the first-order effects Ŝi and Ŝj also based on the mean of the GP only.
Using the same number of Monte-Carlo samples as before (M = 10 000), however, results in a 95% confidence interval with the same order of magnitude as the indices themselves, and even leads to negative values for the Sobol indices (as given in Table 3).Therefore, we increase the number of Monte-Carlo samples to M = 1 000 000 to obtain reasonably small confidence intervals.
Table 3 and Fig. 9 summarise the results for the second-order Sobol indices: the highest interaction is present between P v and L v p , as we already expected based on the results presented in Fig. 8. Summing up all first and second-order Sobol indices results in 0.999.We thus (almost) completely apportioned the variance in the output to the input parameters, including interaction effects.

26/44
The large number of Monte-Carlo samples necessary to estimate the second-order Sobol indices highlights the relevance of metamodel-based estimation approaches.Evaluating the full model f several million times is computationally prohibitive for most models.Without using a metamodel, estimating higher-order Sobol indices is thus impossible in most cases.

Application to a model of arterial growth and remodelling
In the previous sections, we assessed the performance of the current workflow as applied to a model of nanoparticle-mediated drug delivery in more detail: even small numbers of training samples result in reliable estimates of the Sobol indices and a small uncertainty related to the GP metamodel.To give an outlook, we now apply the workflow to another complex biomechanical example, namely a homogenised, constrained mixture model of arterial growth and remodelling [23,24,25].
[95] performed an exhaustive global sensitivity analysis, where they estimated the first and total-order Sobol indices by evaluating the full model for all Monte-Carlo samples.This however entails a large computational burden (> 70 000 model evaluations).Therefore, the question arises as to whether we can reduce this computational cost by using a GP metamodel and still get reliable Sobol index estimates, including reliable uncertainty estimates.
For this comparison, we investigate Case 2 of the original publication [95], where the maximum diameter of an idealised cylindrical abdominal aorta was studied 15 years after spontaneous damage to elastin.In this case, the majority of samples lead to minor dilatation of the vessel d max < 3 cm.In contrast, a considerable number of samples do not stabilise and keep enlarging, leading to aneurysms with a much larger diameter d max >> 3 cm (see Fig 4b in [95]).We use the original results from [95] as a reference and compare them to our results based on the metamodel approach.
First, we take a look at the predictive quality of the GP metamodel.Once again, we use a tensorised, squared, exponential covariance function and compare the results for different numbers of training samples, N = [40,60,80,100,150,200,300,500] in this case.As an example, Fig. 10A presents the training samples for N = 300 for two parameters.The results for the remaining parameters are included in the Supplement (see Fig A .1 in Supplement).We see that the majority of samples result in a small dilatation in contrast to the fewer aneurysmatic samples with a very large diameter of up to 8 cm.This bimodal structure of the data makes training the GP metamodel more difficult compared to our previous example.Accordingly, the Nash-Sutcliffe presented in Fig. 10B is lower, particularly for small numbers of training samples, i.e.N < 80.
Second, we calculate the first and total-order Sobol indices and respective uncertainties for different numbers of training samples.To this end, we use M = 10 000 Monte-Carlo samples, N GP = 500 realisations of the GP metamodel, and B = 300 bootstrap samples.By way of example, we present the results for two parameters, the gain parameter kσ and the initial volume fraction of elastin φ el t0 , in Fig. 10C.Similar plots for the remaining eight parameters are included in the Supplement (see Fig. A.2 in Supplement).For the gain parameter kσ , the estimates based on the metamodel converge to the reference values from [95] for both the first and the total-order Sobol index.For the initial volume fraction of elastin φ el t0 , the reference values for the Sobol indices are very small (0.01 or smaller).In this case, exact estimates based on the metamodel approach are much harder to achieve: for N ≤ 300 training samples the estimates only stabilise.Nevertheless, we can still reliably separate the three most influential parameters from the non-influential parameters, even for small numbers of training samples (see Fig. 28/44 in Supplement).One further detail should be mentioned as an example: the plot for φ el t0 reveals problems in estimating the first and total-order Sobol indices for N = 80 or 100.The Sobol indices and the uncertainties are all close to zero.Similar behaviour can be observed for other parameters (see Fig. A.2 in Supplement).For N < 300, the GP has not yet converged, and hence does not capture all features of the quantity of interest.Furthermore, we note that the uncertainty related to the metamodel is much higher and in some cases even dominates the total uncertainty in this example, while the uncertainty related to Monte-Carlo integration dominated the previous example.It is therefore important to include the uncertainty related to the metamodel because considering only the uncertainty related to Monte-Carlo integration would underestimate the total uncertainty in the Sobol index estimate.
Finally, the computation of higher-order indices was not feasible with the approach chosen in the original contribution [95].In contrast, the following will show that the metamodel-based approach enables their computation.As an example, we again consider the gain parameter kσ : a closer look at Fig. 10C reveals that the total-index is considerably higher than the first-order index: S T kσ − S kσ = 0.34.This delta indicates interactions with other parameters, and thus estimating the second-order indices is of particular interest for this example.Since we see in Fig. 10C that the metamodel contributes significantly to the total uncertainty, we include uncertainty estimates for the metamodel (as opposed to relying solely on the predictive mean, as in the previous second-order estimates).The estimates indeed show interaction with two parameters: the turnover time τ 3 and the constitutive parameter k 2 (S kστ = 0.19 and S kσk2 = 0.06).However, the sum of all second-order indices (S kσj = 0.25) still does not cover the delta between the first and total-order index.Hence, we specifically estimate the thirdorder Sobol index for the three most influential parameters, kσ , τ , and k 2 , resulting in considerable third-order interaction: S kστ k2 = 0.06.Detailed results for the second and third-order indices, including confidence intervals, are included in the Supplement (see Thus, we are able to identify influential parameters for factor prioritisation based on the metamodel approach with small numbers of training samples, and we can also separate the influential parameters from the non-influential ones.Hence, the metamodelbased approach provides the same results as the approach based directly on the full model in the original publication [95].The computational cost, i.e., the number of evaluations of the full model, however, is much lower when using a metamodel.Additionally, we can quantify higher-order indices which is infeasible based on evaluations of the full model.

Conclusion
Since a global sensitivity analysis is computationally expensive, modellers often rely on local methods alone, which may be inadequate [5].The use of a metamodel-based approach, however, allows a global variance-based sensitivity analysis to be performed, even for computationally expensive biomechanical models with a moderate number of input space dimensions at a manageable computational cost.The number of training samples required to obtain reliable estimates for the Sobol indices depends largely on the problem set-up itself: our results demonstrate that we can identify the most influential input parameters and separate them from non-influential parameters with small numbers of training samples.However, quantifying the exact value of the Sobol indices requires more training samples.Moreover, the approach is able to quantify the uncertainty related to the metamodel: including this uncertainty is important, because considering only the uncertainty related to Monte-Carlo integration could underestimate the total uncertainty 3 The turnover time was denoted by T in the original publication [95].

29/44
in the Sobol index estimates.The metamodel-based approach also allows an estimation of higher-order Sobol indices, and thus a quantification of interaction effects, which is not feasible without a metamodel due to the computational costs involved.While there is no one-size-fits-all rule, the approach is general and efficient enough to allow a study of different aspects of sensitivity analysis, including a transparent declaration of the uncertainties involved in the estimation process.
We demonstrated how a rigorous global sensitivity analysis can be applied to complex, computationally expensive problems.A carefully performed sensitivity analysis is generally an integral part to ensure the high quality of any model development [5].By demonstrating the workflow and its application for biomechanical problems, we contribute to closing the gap between proposals of new sensitivity analysis methods and application papers [8].We hereby encourage sensitivity analysis in general and the metamodel-based approach in particular in the biomechanics community.In the big picture of model development, the presented workflow can be a building block towards inverse analysis, or it can be a valuable tool to better understand the model itself.us, summing up the rst-order index S i plus all second-order indices S i j and comparing this sum to the total-order index S Ti reveals whether we expect third-order or even higher order interaction e ects.For all three parameters kσ , τ and k , the total-order index is higher than the sum of rst-and second-order indices S Ti > S i +

A
if second-order indices are included).This means 2M simulations are needed for computing f (A) and f (B) plus D • M simulations needed for computing f (A (i) B ) for i = 1, . . ., D. The cost of first and total-order indices is, hence, M (D + 2) simulations.If second-order indices are included, we need an additional D • M simulations for f (B (i)

Figure 2 .
Figure 2. Schematic overview of the steps involved in training the Gaussian process metamodel.Conditioning the prior Gaussian process on the training samples results in the posterior Gaussian process.
Sample one realisation of the GP surrogate model at M Monte-Carlo samples 1 2 3 4 5 6 7 8 M ... Repeat for all bootstrap samples: b = 2,..., B Repeat for all realisations: k = 1,..., N GP return: (N GP x B) estimates for the Sobol index 5 . The governing equation for the transport of nanoparticles in IF is the mass balance equation of nanoparticles with mass fraction ω NPl ε l ∂ω NP l ∂t − k l µ l ∇p l • ∇ω NP l convective transport − ∇ • ε l D NPl ∇ω NP lwhere the effective diffusivity of nanoparticles in IF is given by D NPl .The superscript l denotes the IF as one of the fluid phases filling the pore space of the extracellular matrix.The IF is characterised by its viscosity µ l , density ρ l , permeability tensor k l , and finally the IF pressure p l resulting from the mass balance of the fluid equation as part of the tumour-growth model.The last term results from employing the product rule, see[61,20].The mass transfer of nanoparticles to and from the IF includes three terms

Figure 4 .
Figure 4. Grown tumour in its microenvironment (A, B, C) with distribution of nanoparticles (D) and necrotic tumour cells (E). A. Volume fraction of living tumour cells ε LTC ; B. Volume fraction of vasculature ε v ; C. Pressure p l in the IF; D. Mass fraction of nanoparticles in the IF ω NP l; E. Mass fraction of necrotic tumour cells ω N t; Subfigures D and E as an example present the result for the mean values of the uncertain input parameters.

Figure 5 .
Figure 5. Nash-Sutcliffe efficiency Q 2 for different training sample set sizes with a tensorised, squared, exponential covariance function and a tensorised Matérn covariance function.N training samples were randomly generated for the main plot.For the detail plot, we repeated the process five times with different training sample sets of sizes N =[10,15,20,25,30].For reference, the dashed lines in the detail plot are identical to the dashed lines in the main plot.

Fig. 7 Figure 6 .
Figure 6.Gaussian process for N = 20 with a tensorised, squared, exponential covariance function.Projected mean m GPi (X i ), projected 95% confidence interval (CI) and training samples D for the mean of the necrotic fraction of tumour cells ωN t

Figure 7 .
Figure 7. Convergence of first-order Sobol index estimate for increasing training sample set size.We use M = 10 000 Monte-Carlo samples and N GP = 500 realisations of the Gaussian process.

Table 2 .Figure 8 .
Figure 8. First-order and total-order Sobol indices and 95% confidence intervals (CI) for an increasing number of training samples.We use M = 10 000 Monte-Carlo samples, N GP = 500 realisations of the Gaussian process, and B = 300 bootstrap samples.Subfigure A shows Sobol indices with metamodel CI and the sum of metamodel and Monte-Carlo CI for the six input parameters separately.Monte-Carlo abbreviates as MC.Subfigure B details the metamodel CI.

Figure 9 .
Figure 9. Second-order Sobol indices Ŝij for M = 1 000 000 Monte-Carlo samples and B = 300 bootstrap samples.The blue circles represent the first-order Sobol indices.The orange surrounds represent the total-order Sobol indices.The grey areas connecting the nodes represent the second-order Sobol indices.All other second-order indices are Ŝij < 0.001.

Q 2 Figure 10 .
Figure 10.Summary of results for arterial growth and remodelling.A. Projected mean m GPi (X i ), projected 95% confidence interval (CI) and training samples D for the diameter d max (the y-axis labels apply to both figures).B. Nash-Sutcliffe efficiency Q 2 for different training sample set sizes.C. First-order and total-order Sobol indices and 95% confidence intervals (CI) for an increasing number of training samples.Monte-Carlo abbreviates as MC.

Fig. A.Figure
Fig. A. depicts the projection of the (trained) Gaussian process into the input-space dimensions including the projection of the mean m GPi (X i ), projected  con dence interval (CI), and training samples D for the diameter d max .

Fig
Fig. A. presents the estimates for the rst and the total-order Sobol indices for increasing numbers of training samples.e results are based on M =Monte-Carlo samples, N GP = realisations of the Gaussian process, and B = bootstrap samples.We additionally include the  con dence intervals (due to the metamodel and the sum of CI due to the metamodel and Monte-Carlo integration).
Figure A. : First-order and total-order Sobol indices and con dence intervals (CI) for an increasing number of training samples.We use M = Monte-Carlo samples, N GP = realisations of the Gaussian process, and B = bootstrap samples.

Finally, we estimate
the third-order Sobol indices for the parameter combination kσ , τ and k .e total-order Sobol index S Ti includes all contributions of parameter X i S Ti = S i + D j≠i S i j + ⋅ ⋅ ⋅ + S ...D .(A. ) one would need M samples to calculate the inner expectation and then repeat this M times to calculate the outer variance, resulting in a computational cost of O(M 2 ) [6, p. 164].Since M usually has to be large 2 , this is impractical, especially considering that we would need to evaluate the full model f for each Monte-Carlo sample.

Schematic summary of the components of the multiphase tumour- growth model.
The model comprises the ECM, as the solid phase, with three fluid phases filling its pore space: host cells, tumour cells and IF.The vasculature is included as an independent porous network.Species transported by the different phases are: oxygen, nanoparticles, necrotic host cells and necrotic tumour cells (marked in italics).
Figure 3.We consider three mass transfer terms that transport nanoparticles into and out of the IF: A. the interendothelial pathway; B. the transendothelial pathway; C. lymphatic drainage (all transfer terms are marked in bold).Drugs mediated by the nanoparticles kill tumour and host cells (marked by red lightning).The lymph system is not explicitly included in the model.

Table 1 . Probability distributions of uncertain input parameters.
We assume that all parameters are distributed uniformly within the given range.

Table 3 . Second-order Sobol indices for
M = 1 000 000 Monte-Carlo samples and B = 300 bootstrap samples.All other second-order indices are S ij < 0.001.