Fast uncertainty quantification of tracer distribution in the brain interstitial fluid with multilevel and quasi Monte Carlo

Abstract Efficient uncertainty quantification algorithms are key to understand the propagation of uncertainty—from uncertain input parameters to uncertain output quantities—in high resolution mathematical models of brain physiology. Advanced Monte Carlo methods such as quasi Monte Carlo (QMC) and multilevel Monte Carlo (MLMC) have the potential to dramatically improve upon standard Monte Carlo (MC) methods, but their applicability and performance in biomedical applications is underexplored. In this paper, we design and apply QMC and MLMC methods to quantify uncertainty in a convection‐diffusion model of tracer transport within the brain. We show that QMC outperforms standard MC simulations when the number of random inputs is small. MLMC considerably outperforms both QMC and standard MC methods and should therefore be preferred for brain transport models.


| INTRODUCTION
Mathematical models in biology involve many parameters that are uncertain or in some cases unknown. Over the last years, increased computing power has expanded the complexity and increased the number of degrees of freedom of many such models. For instance, full scale modeling of the brain requires very high spatial resolution to resolve the details on the curvature of the surface. For this reason, efficient uncertainty quantification algorithms are now needed to explore the often large parameter space of a given model. Advanced Monte Carlo methods such as quasi Monte Carlo (QMC) and multilevel Monte Carlo (MLMC) have become very popular in the mathematical, engineering, and financial literature for the quantification of uncertainty in model predictions. However, applying these methods to physiologically relevant brain simulations is a difficult task given the typical complexity of the models and geometries involved.
Efficient elimination of waste products from the brain is crucial for a functional central nervous system. In humans, the brain is responsible for around 20% of the total oxygen consumption and 15% of the cardiac output. 1 Despite this high energy demand, the brain lacks lymphatic vessels, which carry waste products in the rest of the body, thus introducing the need for an effective alternative mechanism. Proposed theories of waste clearance all include a combination of diffusion and convection to clear substances out of the brain. There is also consensus that interstitial fluid (ISF) within the brain can exchange with cerebrospinal fluid (CSF) surrounding the brain to transfer waste products. However, there is disagreement on where the CSF/ISF-exchange occurs, in which direction it occurs, the forces involved in transport and how much of the transport within the brain that can be explained by diffusion alone. [2][3][4][5][6][7] Interchange between CSF and ISF is also an important mechanism for effective drug delivery to the brain by substance administration in the SAS. The blood-brain-barrier is one of the main obstacles for this type of drug delivery to the brain. 8 The contrast agent gadobutrol has been used as a CSF tracer to assess transport into brain tissue, and possibly across the blood-brain-barrier. 9,10 However, human studies are sparse, the number of subjects are limited, and MRI cannot directly capture phenomena on the microscale (eg, transport across the blood-brain-barrier).
Mathematical modeling is therefore a compelling tool to investigate CSF/IFS-exchange related to brain clearance or drug delivery. Computational studies [11][12][13] investigating CSF flow in paravascular spaces (PVS) have not been able to reproduce the average flow velocities reported in experimental studies with microspheres. 14,15 Within the interstitium, both experimental, 7 and computational 16 studies have concluded that diffusion dominates convection. Still, tracer experiments in human beings 9,10 show transport of gadobutrol to the brain that is unlikely to be explained by diffusion alone. 17 The aforementioned mathematical models involve parameters that are usually unknown to some degree. Therefore, evaluating the models' parameter sensitivities is vital to give confidence in the model predictions and conclusions. For instance, brain tissue permeability range over several orders of magnitude in the literature. 16,18 Recently, Fultz et al 19 found coherent oscillations of electrophysiological waves, blood flow and cerebrospinal fluid flow during non-rapid eye movement sleep, suggesting that these processes are all interlinked. Full-scale modeling of CSF/ISF flow in the brain thus introduces several more interlinked parameters related to these additional processes. Previous studies in the literature that have included extensive parameter sensitivity analyses have been computationally cheap, 20,21 allowing for parameter space exploration within a reasonable amount of time. However, when the model of interest is expensive to simulate, more advanced uncertainty quantification (UQ) methods are needed.
Standard Monte Carlo (MC) methods have successfully been applied to simulations of, for example, the cardiovascular system, [22][23][24] and to brain solute transport. 17 However, when working with physiologically realistic, MRIderived geometries, standard MC methods are typically prohibitively costly and more advanced methods, for example, quasi Monte Carlo (QMC) or multilevel methods (multilevel, multilevel quasi, or multi-fidelity Monte Carlo), are to be preferred. We remark that each of these methods bring their own additional requirements: QMC requires either the output functionals of interest to have low effective dimensionality with respect to the random input and/or the input dimensions to be ordered in order of decaying importance. 25 Whether this assumption is satisfied is strongly problem-dependent. On the other hand for Multilevel Monte Carlo (MLMC) to work well, a good statistical coupling must be enforced between the levels and a mesh hierarchy on which a good rate of bias and variance decay can be appreciated is needed. 26 Even though such a hierarchy can always be obtained by refining and/or coarsening a given mesh in theory, this is far from trivial in practice. A multi-fidelity Monte Carlo (MFMC) 27 approach can be beneficial in this case, since it can incorporate in the hierarchy low-fidelity models that are still correlated with the fine-mesh simulations, but do not require a computational grid. However, suitable low-fidelity models are not always available in practice. Finally, multilevel quasi Monte Carlo (MLQMC) methods combine the advantages and the requirements of both QMC and MLMC and, as such, may be advantageous only in the cases in which both QMC and MLMC perform well.
In light of the above considerations, the efficiency, feasibility and performance of enhanced Monte Carlo methods such as QMC, MLMC, MFMC, and MLQMC is highly problem-dependent. Moreover, determining a-priori which method that will bring the largest performance benefits for a given problem is highly nontrivial. As the model complexity and dimension grows, choosing the most efficient UQ method becomes even more important. Additionally, knowing when QMC and MLMC both work well provides a model domain for which MLQMC methods can bring substantial additional computational improvements. Generally speaking, correctly setting up a QMC or MLMC method in complex geometries is a difficult task, and to our knowledge QMC and MLMC methods have not yet been employed for UQ in simulations of brain physiology.
In this study, we therefore investigate different enhanced Monte Carlo approaches for quantifying uncertainty in a model with random field coefficients deriving from MRI-studies in humans. 9,10 The model, and its connection to tracer experiments in humans, has previously been described in detail, 17 as well as the clinical relevance of the simulated tracer experiments in Ref. 10. We introduce and evaluate the relative performance of different Monte Carlo methods, namely standard, quasi and multilevel Monte Carlo methods 29,30 as implemented in the FEMLMC library. The software platform used is comprised of generally available open-source software libraries (FEniCS, 31 libsupermesh, 32,33 and FEMLMC 34 ).
Our findings show that QMC outperforms standard Monte Carlo, but only when the number of (in this case infinite-dimensional) random inputs is small, bringing a 10-fold improvement in the computational cost. On the other hand, MLMC always outperforms both standard Monte Carlo and QMC, and is two orders of magnitude faster than standard Monte Carlo, showing that MLMC should be the method of choice when solving brain transport problems. Our results further suggest that problems with a small number of random field inputs might be amenable to a MLQMC treatment.
This paper is organized as follows. In Section 2, we give a brief overview of mathematical, numerical and algorithmic background. In Section 3, we present the baseline stochastic model for brain tracer movement, as originally introduced in Ref. 20, and in Section 4 we formalize two UQ test problems building on this baseline model. We detail the numerical and computational solution algorithms in Section 5 before presenting numerical results in Section 6. We discuss our findings in Section 7 and conclude in Section 8.

| Preliminaries
In what follows we indicate by x ∈ G & R 3 a given spatial point in a domain G, with t > 0 a time point, and with ω a given probabilistic event, living in a sample space Ω. For example, we might indicate with f(t, x) a generic function of time and space and with z(ω) a generic random variable of expected value  z ½ and variance  z ½ . Additionally, we will consider random functions of space, varying both in space and across random realizations. These are called random fields and we use the notation, for example, u(x, ω). More formally, a random field u(x, ω) is a collection of random variables such that each point value u(x, Á) is a random variable for every x and u(Á, ω) is a function of space for fixed ω.
For Gaussian fields, all point values (which are random variables) are jointly Gaussian, and a Gaussian field u(x, ω) is uniquely determined by prescribing a mean function μ x ð Þ =  u x, ω ð Þ ½ and a covariance function . A ubiquitous family of Gaussian fields is the Matérn family: a Matérn field is a Gaussian field with covariance of the type where σ 2 , ν, λ > 0 is the variance, smoothness parameter and correlation length of the field, respectively. Here Γ(x) is the Euler Gamma function and K ν x ð Þ is the modified Bessel function of the second kind. The σ 2 , ν and λ parameters determine the characteristics of the random field. For instance, in our model the random field represents CSF/ISF-flow and the correlation length was chosen to match approximately the distance between larger arteries and veins, to be able to model CSF inflow along arteries and outflow along veins. 4 In this paper we consider the solution of a convection-diffusion-reaction equation with random fields as coefficients for the movement of tracer within the brain. The convection-diffusion-reaction equation reads as: find the tracer concentration c = c(t, x, ω) for x ∈ G, ω ∈ Ω and t ≥ 0 such that Here, the domain G & R 3 represents the brain, the superimposed dot represents the time derivative, D * is the effective diffusion coefficient of the tracer in the tissue, v represents a convective fluid velocity and r ≥ 0 is a drainage coefficient. Further details, including boundary and initial conditions, will be presented in Section 3.
The solution c is random and varies in time and space as well, while the coefficients v(x, ω) and D * (x, ω) are random fields. A typical assumption is that one is interested in computing one (or more) output functional of interest Q ω ð Þ= Q c ð Þ,for example, Q could be the spatial average of c over a region of interest. Typically, solving (2) then means to compute the expected value of Q(ω). Note that most other statistics (such as the variance or the cumulative density function [CDF]) can be rewritten as an expectation and thus computed analogously. A basic method for solving (2) is the standard Monte Carlo (MC) method, for which the expectation is approximated by the sample average obtained with 0 < N ∈ N samples of Q(ω). Note that each sample of Q requires computing samples of the coefficients v (x, ω), D * (x, ω) and the corresponding solution of (2). Alas, the MC method converges slowly in terms of the number of samples N, namely at a rate of O N − 1=2 À Á . This makes standard MC quite expensive: assuming an ε 2 tolerance for the mean square error (MSE)  Q −Q À Á 2 h i and a cost per sample of ε −q for some positive problem-specific number q, standard MC has a total cost complexity of O(ε −2 − q ), see Ref 26. In practice, more advanced methods, such as Quasi Monte Carlo (QMC), multilevel Monte Carlo (MLMC) and multilevel quasi Monte Carlo (MLQMC), are to be preferred. In what follows we give a brief description of the QMC and MLMC methods.

| Quasi Monte Carlo
We now give a quick overview of the quasi Monte Carlo method and we refer the reader to the book by Lemieux 35 for a more in-depth description. The advantage of QMC with respect to standard MC is that in QMC methods, the convergence rate with respect to the number of samples N is improved from O(N −1/2 ) to O(N −1 + ε ) for any ε > 0. At the heart of QMC is the approximation or reinterpretation of the expectation as a high-dimensional integral over the unit hypercube: where P is a suitable probability measure, H(x) for x ∈ R s is some suitable function and s is typically the dimensionality of the random input needed to sample Q. The QMC method can then be expressed as a quadrature rule over [0, 1] s with equal weights, approximating the right-hand side integral with with x n ∈ R s . Choosing the x n f g N n = 1 uniformly at random results in a convergence rate of O(N −1/2 ). However, there exist deterministic point sequences, so-called low-discrepancy point sequences that can yield an improved rate of O (N −1 (logN) s ) . 36 This is the key idea of QMC methods: estimating the integral I by a low-discrepancy sequence. Now, let  Q ½ be the variance of Q. In standard MC, the statistical error is O  Q ½ =N ð Þ 1=2 (owing to the central limit theorem) and  Q ½ can be estimated by taking the sample variance of Q. However, in standard QMC, the point sequence x n f g N n = 1 is deterministic and therefore there is no notion of statistical error available. For this reason, a practical estimate for the approximation error introduced by QMC is not available. A common solution is to randomize the sequence (while still preserving the hypercube-covering properties) in order to retain a measure of estimator variability (see eg, Chapter 6 of Ref. be M independent randomizations of a low-discrepancy sequence x n f g N n = 1 , then randomized QMC estimates I as where each of the I m N ω ð Þ are now random and a confidence interval can therefore be estimated provided M is large enough. In this paper we use M = 32. Assuming a fixed M, a given MSE tolerance ε 2 , a O(ε −q ) cost per sample and a QMC convergence order of O(N −1 + ε ) for any ϵ > 0, the total cost of randomized QMC is O(ε −q − 1/(1 − ε) ). If we take ϵ to be small, this is almost ε −1 times better than standard MC. Remark 1. The (logN) s can dominate the early convergence behavior of QMC methods, in which case the suboptimal O (N −1/2 ) convergence is observed initially. This behavior is typically exacerbated when the random input dimensionality s is particularly large (eg, in applications with infinite-dimensional random inputs such as random fields), to the extent that the faster rate is never observed in practice. However, this is not always the case: if the QMC integrand H(x) is inherently lower dimensional in the sense that it can be well approximated by a function only depending on s ( s input dimensions, then it is possible to replace s with s and the transition to a faster rate will occur earlier. This is the principle of low-effective dimensionality, first introduced by Caflish et al in Ref. 13, and is fundamental when using QMC in high dimensions. Thus, for good QMC convergence it is important to either reduce the integrand dimensionality or to order the integration variables in order of decaying importance such that the improved convergence rates can be attained.

| Multilevel Monte Carlo
The multilevel Monte Carlo method was first introduced by Heinrich in Ref. 77 for parametric integration and subsequently popularized by Giles for stochastic differential equations in Ref. 32. MLMC works under the assumption that one can compute a hierarchy Q ℓ ω ð Þ f g L ℓ = 1 of different approximations of the output functional Q of increasing accuracy. For instance, one could consider solving (2) on a hierarchy of computational meshes of size decreasing with ℓ. At the heart of MLMC is the expansion of  Q ½ into the telescoping sum Approximating each expectation in the sum on the right-hand side with standard MC then yields the MLMC estimator: in which a key element is that the term is sampled by the same event ω n ℓ . This aspect is referred to as the MLMC coupling and is essential for the improved performance of MLMC over MC. In fact, while the variance of each Q ℓ might be large, the variance of the difference Q ℓ − Q ℓ − 1 is typically much smaller due to the strong correlation between the two terms. For this reason, while the estimation of each  Q ℓ −Q ℓ − 1 ½ term still occurs at a O(N −1/2 ) rate, the number of samples needed to achieve a given statistical error tolerance is smaller and decreases with ℓ. More formally, the convergence of MLMC is ensured by the following Theorem 39 : Theorem 1. (Ref 26, Theorem 1). Let Q(ω) be a random variable with finite second moment and let Q ℓ be its level ℓ approximation. Let Y ℓ be the (unbiased) MC estimator of  Q ℓ −Q ℓ− 1 ½ on level ℓ, and let C ℓ and V ℓ be the expected cost and variance, respectively, of each of the N ℓ Monte Carlo samples needed to compute Y ℓ . If the estimators Y ℓ are independent and there exist positive constants α, β, γ, c 1 , c 2 , c 3 , such that α≥ 1 2 min β, γ ð Þ and then there exist a positive constant c 4 such that, for all ε < e −1 , there is a level number L and number of samples N ℓ , such that the mean square error (MSE) of the MLMC estimatorQ = P L ℓ = 1 Q ℓ is bounded: and its total computational complexity is bounded: Remark 2. The MLMC parameters α, β, and γ must be estimated if not known a priori. However, for PDE applications, as it is the case here, a priori error estimates are typically available. The optimal number of samples on each level N ℓ and the maximum level L needed can similarly be estimated.
Remark 3. Theorem 1 assumes that one can increase the maximum level L at will. In practical applications this might not always be possible, for example, when the computational resources available are limited.
In this study, we modify the original MLMC algorithm by weighting the relative importance of bias and statistical error as introduced by Haji-Ali et al. 38 The MSE of the MLMC estimator is given by whereQ is the MLMC estimator with varianceV . To ensure that MSE ≤ ε 2 , we enforce the bounds, where θ ∈ (0, 1) is a weight balancing the two terms so as to obtain comparable error reduction. Small values of θ reduce the number of samples needed and are therefore preferred when the bias is small. Conversely, large values are beneficial when the bias is large as they allow to achieve smaller tolerances without adding finer levels to the hierarchy.

| Gaussian field sampling techniques
When Q(ω) depends on a random field (eg, through the coefficients of (2)), an efficient sampling technique is also needed as part of any Monte Carlo method. Unfortunately, sampling generic random fields from a given distribution is normally a prohibitive task. In the specific case in which the field to be sampled is Gaussian, however, the sampling problem becomes tractable, albeit still computationally expensive. Different Gaussian field sampling methods are available in the literature. The most common are: (a) Factorization of the covariance matrix, including Hierarchical matrix approximation [40][41][42][43] (b) Karhunen-Loève expansion of the random field (cf. section 11.1 in Ref. 79 and eg, Ref. 78) and (c) the circulant embedding method which uses FFT. [46][47][48][49] In this paper we use an alternative method known as the stochastic PDE (SPDE) approach for the sake of efficiency. 29,30,50,51 This sampling strategy is based on approximately drawing realizations of a Matérn field by solving the following elliptic PDE 50,51 : where ν is the smoothness parameter of the Matérn covariance to be sampled (cf. Equation (1) , and the equality has to hold almost surely and be interpreted in the sense of distributions. When k is an integer, as we shall assume here, solving (14) is equivalent to solving a second-order PDE k times. The constant η > 0 is a scaling factor that depends on the Matérn covariance parameters σ, λ, and ν, cf. Ref. 30. The solution to u(x, ω) is a Matérn field only ifĜ R d . Otherwise, for generalĜ & R d , u is still a Gaussian field, but its covariance is only approximately Matérn, with the error in the covariance decaying exponentially away from the boundary ofĜ . 52 For this reason in this paper we will always assume that the field sample is needed on a domain G &Ĝ and that we solve (14) on a domainĜ large enough so that this source of error is negligible over G. We impose homogeneous Dirichlet boundary conditions on ∂Ĝ.
Note that the term _ W is random and for each sample of u needed, we must sample _ W from its distribution and solve (14). The sampling of _ W is non-trivial in itself, especially when samples are needed within a MLMC or QMC framework: in the QMC case, a good variable ordering is required to achieve good convergence with respect to the number of samples while in the MLMC case the multilevel coupling must be enforced correctly to obtain a considerable variance reduction. In this paper, we use the sampling techniques developed in Ref. 19 for the MLMC case and in Ref. 18 for the (ML)QMC case to address these requirements. Both these MLMC and (ML)QMC sampling strategies are efficient with optimal cost complexity, making the sampling of u a linear cost and memory operation. This optimal complexity is especially advantageous when dealing with physiologically relevant applications with complex geometries in 3D.

| A STOCHASTIC MODEL OF TRACER TRANSPORT IN BRAIN TISSUE
This paper aims to evaluate the numerical and computational performance of different Monte Carlo methods for quantifying uncertainty in stochastic models of tracer transport in brain tissue. In particular, we focus on two comprehensive coefficient models with random diffusion and velocity fields. In addition to reflecting existing hypotheses on ISF tracer transport, these models are designed to offer a suitable challenge to the MLMC and QMC algorithms presented in Refs. 18, 19.

| The ISF tracer transport equation
We begin by considering the baseline model derived by the authors in Ref. 20, and briefly introduced in Section 2.1, describing the evolution of tracer concentration within the brain parenchyma under uncertainty: find the tracer concentration c = c(t, x, ω) for x ∈ G, ω ∈ Ω and t ≥ 0 such that As before, G & R 3 is the brain parenchyma domain comprised of white and gray matter from the Colin27 human adult brain atlas FEM mesh 53 version 2 ( Figure 1). The superimposed dot represents the time derivative, D * is the effective diffusion coefficient of the tracer in the tissue (depending on the tracer free diffusion coefficient and the tissue tortuosity), 54 v represents a convective fluid velocity and r ≥ 0 is a drainage coefficient potentially representing, for example, capillary absorption 55 or direct outflow to lymph nodes. 10 The parenchymal domain is assumed to not contain any tracer initially: c(0, x, ω) = 0. This model reflects the MRI-study of Ringstad et al 10 in which 0.5 mL of 1.0 mmol/ mL of the radioactive tracer gadobutrol was injected in the spinal canal (intrathecally) of 15 hydrocephalus patients and eight reference subjects.

| Boundary conditions
Let the brain boundary ∂G = ∂G S [ ∂G V , with ∂G S representing the interface between the brain parenchyma and the subarachnoid space (SAS), and ∂G V representing the interface between the brain parenchyma and the cerebral ventricles, respectively. We assume the tracer concentration on the SAS interface to be known and impose no ventricular outflux. As boundary conditions for (15), we thus prescribe where n is the unit normal vector pointing outward from ∂G. The function g(t, x) models the movement of tracer starting from the lower cranial SAS and traveling upward in the CSF surrounding the brain as observed in the study by Ringstad et al, 10 and the form where x = [x 1 , x 2 , x 3 ] and c CSF (t) is the average tracer concentration in the SAS, while h(t, x) represents its spatial distribution. The variable υ z = 1.5 × 10 −5 m/sec is the speed of tracer movement upwards in the SAS, while a = 20 m −1 reflects the gradient of tracer concentration from the lower to the upper cranial SAS. The value z 0 = − 0.2 m is the initial distance from the lateral ventricles reached by the tracer at t = 0. The average SAS tracer concentration c CSF is modeled as follows. Let n 0 = 0.5 mmol be the total amount of tracer initially injected in the CSF and let V CSF = 140 mL be the total CSF volume in the human SAS and ventricles. 56 Then, the average concentration in the SAS right after injection is given by c CSF (0) = 0.5 mmol/140 mL = 3.57 mol/m 3 . Assuming conservation of tracer molecules, the total amount of tracer in the brain and in the SAS plus or minus the tracer otherwise absorbed or produced is constant in time, and is equal to the initial amount n 0 . This observation gives the approximate relation where, for simplicity, c is given by a deterministic solution of (15) with boundary conditions (18) in which all the random coefficients are replaced by their average. Solving for c CSF , we finally let The computational domains and points of interest. A, The coarsest computational mesh, embedded in the auxiliary box used to sample the Matérn fields. B, The domain representing the brain parenchyma. The lateral ventricles are shown in red, while the parenchyma is pink. Two smaller regions of interest are marked in green: The leftmost green region, S w , is within the white matter, while the rightmost green region, S g , is within the gray matter

| Quantities of interest
We consider different output quantities (functionals) describing the characteristics of tracer movement within the brain. For each time τ = 30k min for k = 1, …, 48 (from half an hour from injection to 1 day after), we consider the total amount of tracer in the gray matter Q g and in the white matter Q w : Additionally, we consider two localized concentration measures: the average tracer concentrations q g and q w in two smaller regions, one within the gray matter S g and one within the white matter S w respectively: where V g and V w is the volume of the gray and white matter subregions, respectively. The size and location of these subregions are shown in Figure 1B.

| COEFFICIENT MODELS
The effective diffusion coefficient of a solute (or tracer) and the velocity field are heterogeneous 6,57-59 within the parenchyma and also vary from individual to individual. To account for both these types of variation and for the uncertainty in the coefficient magnitude we model them as random variables or fields.

| Diffusion coefficient
Let D * Gad = 1:2 × 10 − 10 m=s 2 be the average parenchymal gadobutrol diffusivity. 9 We model the effective diffusion coefficient as where D * γ is a random field such that for each fixed x ∈ G, D * γ x, ω ð Þ is a gamma-distributed random variable with shape k = 3 and scale θ = 0:75 × D * Gad =k . This choice of parameters ensures the positivity of D * with probability 1. Furthermore, it reflects the average value and variability reported in the literature, namely we have  D * Â Ã = D * Gad , with values larger than 3 times the average being attained with very low probability. 9,54 The probability distribution of D * (x, Á) is shown in Figure 2B.
To sample D * γ from its distribution, we first sample a Matérn field X(x, ω) using the techniques presented in Refs. 18, 19, and then transform it into a gamma random field by using a copula. 60 This consists in setting where F −1 is the inverse cumulative density function (CDF) of the target (gamma) distribution, Φ is the CDF of the standard normal distribution and X(x, ω) is a standard (zero mean, unit variance) Matérn field with smoothness parameter ν = 2.5 and correlation length λ = 0.01 m, cf. (1). Note that Φ maps any standard normal random variable to a standard uniform random variable and that F −1 maps any standard uniform random variable to the target distribution, hence the function F −1 (Φ(x)) maps standard random variables to the target gamma distribution. Samples of D * γ x, ω ð Þ obtained this way will preserve the same Spearman correlation and smoothness properties of X (x, ω), but will present a different covariance structure. 60

| Velocity and drainage coefficients
We now turn to define two models (Model 1, 2) for brain tissue fluid movement and clearance. Both models are combined with the random diffusion field defined by (23). For further physiological considerations, see, for example,. 17

| Model 1-Glymphatic velocity model with directionality
For Model 1, we model fluid transport along paravascular spaces in direct communication with the SAS 6 under the following assumptions. (a) Substantial changes in the velocity field occurs at a distance proportional to a characteristic distance (correlation length λ) between an arteriole and a venule. (b) The velocity field can be represented as a sum of a glymphatic velocity field associated with arterioles and venules and a velocity field represents a directional movement due to larger blood vessel structures. (c) Almost no fluid is filtrated or absorbed by capillaries, thus v base is divergencefree, while v dir only has a small net flow of 0.007 mL/min out of the parenchyma. We then let the drainage term r = 0 in (15) and define the stochastic velocity field The stochastic glymphatic velocity field v base is given by where η λ ð Þ = λ= ffiffi ffi 8 p is a scaling constant chosen such that the magnitude of v (denoted jvj) satisfies  v j j 2 Â Ã 1=2 = v avg , U(ω) is an independent standard uniform random variable and X(x, ω), Y(x, ω) and Z(x, ω) are standard i.i.d. Matérn fields with ν = 2.5 and correlation length λ = 1020 μm. Taking the curl of the random vector field [X, Y, Z] T ensures that v base is divergence free. A sample of the streamlines of v base is shown in Figure 2A while its velocity magnitude distribution is shown in Figure 2C. The deterministic second term v dir represents a directional velocity field induced by large vascular structures 17 and is given by where v f = 2 × 10 −6 m/s.
In fact, note that the partial derivative ∂X/∂x i of a zero-mean Gaussian field X(x, ω) with a twice differentiable covariance C(x, y) is still a zero-mean Gaussian field with covariance given by ∂ 2 C(x, y)/(∂x i ∂y i ) (cf. section 2.3 in Ref. 76). The curl components of the field within the brackets in (25) are therefore just sums of independent Gaussian fields, and hence Gaussian as well. Similarly, this also shows that the covariance of v has the same correlation length as the Matérn fields X, Y and Z since the second derivative of a Matérn covariance (cf. (1)) has the same correlation length as the original covariance function, although it is not Matérn anymore.

| Model 2-Capillary filtration model with arterial inflow and sink term
For Model 2, we consider an alternative velocity representation in which CSF enters the brain parenchyma along spaces surrounding penetrating arteries. 6,14,15,61 In this case, the velocity field is taken to be representing a net flow of CSF into the brain. The flow field is radially symmetric and directed inwards from the outer SAS to a spherical region of radius R = 8 cm around a center point x c within the lateral ventricles. Here v ω ð Þ is a gamma-distributed random variable chosen such that the probability distribution of the velocity magnitude is comparable to that of Model 1. The shape parameter is k = 2 and the scale parameter is set such that again  v j j 2 Â Ã 1=2 = v avg . Note that in this model  v ½ ≠0 and the main source of uncertainty is the random variable ( v ω ð Þ) rather than the spatially dependent random field.
Finally, we set a nonzero sink coefficient of r = 1 × 10 −5 s −1 , to model the assumption that ISF is drained along some alternative outflux route within the brain parenchyma. The value of r chosen corresponds to a 40% drainage of the injected tracer over 48 hours.

| NUMERICAL SOLUTION OF THE STOCHASTIC MODELS
More advanced Monte Carlo methods, such as QMC, MLMC and MLQMC, are known to outperform standard Monte Carlo methods, given an appropriate problem setting. In particular, QMC requires either the output functional to be of low effective dimensionality with respect to the random input, or the input dimensions to be ordered in order of decaying importance. 25 Even though our QMC method 29 is designed to expose the leading-order dimensions in each random field input, hence providing a suitable variable ordering, it is not known whether this strategy would prove to be effective for the problem at hand. In fact, two main challenges arise: (a) The state Equation (15) must be solved on a complicated 3D geometry, and (b) the random input dimensionality is large. The latter point is especially relevant in connection with Model 1 in which 4 infinite-dimensional random fields appear as coefficients. Both these challenges significantly increase the input dimensionality, possibly affecting QMC performance.
On the other hand, MLMC brings a different set of performance requirements. A good MLMC coupling is ensured by our coupling technique, 30 but the technique also hinges on the construction of an appropriate mesh hierarchy on which a good rate of decay of bias and variance can be appreciated. Even though such a hierarchy can always be obtained by increasing the mesh refinement in theory, the availability of computational resources and time may limit the maximal refinement level in practice. As it is a priori unclear whether QMC and/or MLMC actually bring any significant advantages with respect to standard MC when solving (15), we here evaluate both algorithms to determine which if any method performs the better.
We refer the reader to the book by Lemieux 35 and to the review article by Giles 26 for further information about QMC and MLMC methods. In what follows, we detail the numerical approach adopted for the solution of (15).

| Weak form and discretization
We solve (15) numerically using the finite element method (FEM). Let H 1 S G ð Þ be the standard Sobolev space of weakly differentiable functions that are zero on the boundary S . 62 For c, s ∈ H 1 (G), we define a c, s ð Þ= div vc ð Þ,s ð Þ+ D * rc, rs whereby (Á, Á) we indicate the standard L 2 (G) inner product: After a (second-order) implicit mid-point time discretization with time step Δt, the continuous weak form of (15) reads: find c n ∈ H 1 (G) such that for all s∈H 1 S G ð Þ and a.s., where c n thus represents an approximation to c(t n , Á) with t n = nΔt for n = 0, …, n T − 1 and n T − 1 = T/(Δt) and c n CSF is an approximation of c CSF (t n ), defined in (20). We approximate c n CSF by approximating the time integral in (20) with the trapezoidal rule: We note that the expression on the right-hand side is known as c can be pre-computed once. This decoupling results in a second-order scheme in time for c. To compute c we adopt the same discretization, but we approximate c n CSF explicitly (by replacing n with n − 1 in the right-hand side of (32)) thus avoiding the non-local, implicit boundary condition. This approximation gives to a first-order scheme for c, for which we compensate by computing the approximation using a very small time step (Δt c = 30 × 2 − 6 min) and we solve for c on the finest mesh available (see later in this section for a description of the meshes and time step sizes used).
We discretize (30) in space using the FEM. Given a FEM approximation subspace V h & H 1 S G ð Þ, the fully discrete weak form of (15) reads: find c n h ∈V h such that, for all s ∈ V h and a.s., where c n,h CSF is given by (32) in which c n − 1 and c i are replaced by c n − 1 h and c i h respectively. We let V h be composed of piecewise linear continuous Lagrange elements defined relative to a mesh T h of G of mesh resolution h.

| Meshes and time steps
We discretize the domain G by using various refinements of the Colin27 human adult brain atlas simplicial mesh 53 (version 2). We construct a multilevel hierarchy in which the coarsest level is given by one uniform refinement of the original brain mesh and the other 2 levels are obtained through uniform refinement. On level ℓ, we fix (Δt) ℓ = 15 × 2 −ℓ min and we terminate the simulations after T = 1 day. The numbers of cells and vertices of the meshes in the hierarchy are given in Table 1. The Matérn fields are sampled on an extended domain, i.e. a mesh of a larger box domainĜ of size sufficiently large to make the domain truncation error negligible (dimensions 0.16 × 0.21 × 0.17 m 3 ), 52 and then restricted to the brain mesh. The outer box together with the embedded Colin27 mesh is shown in Figure 1A. Each outer box mesh is constructed with the meshing software Gmsh 63 (dev version 4.2.3-git-023542a) such that the corresponding brain mesh is nested within. Furthermore, the box meshes are graded such that the cell size gradually gets larger away from the brain boundary (Matérn field values are only needed in the brain domain).

| Numerical stability considerations
For the parameter regimes considered here, the problem (15) is only mildly convection-dominated, with an upper estimate of the Péclet number of whereL≈0:084 m is half the diameter of the computational domain, v avg = 0.17μm/s, and D * Gad = 1:2 × 10 − 10 m=s 2 . Given the fine computational meshes, we obtain low-probability worst-case cell Péclet numbers of ≈43 × 2 −ℓ on level ℓ of the MLMC hierarchy used. In the numerical experiments, convection-related numerical instabilities were not observed.
However, in initial investigations, we observed that the FEM solution undershoots near the boundary, attaining negative concentration values. This is a known phenomenon in the literature and it does not depend on the velocity field, but it is typical of diffusion problems with Dirichlet-type boundary conditions. 64 We address this problem by a mass-lumping technique, which is known to reduce this effect. 64 This ill-behavior disappears as the mesh is refined to the extent that no undershoots were observed on the finer levels of the MLMC hierarchy. In our numerical experiments, MLMC convergence was resilient to non-physical behavior on the coarser levels. In fact, we noted that these coarse level solutions would still act as a good control variate for the finer levels.

| Solver and software
For the computations, we combined the University of Oxford Mathematical Institute computing servers and the University of Oslo Abel supercomputing cluster. The total amount of CPU hours (CPUh) available on the Abel cluster for the project were approximately 400 000 CPUh with the rule that computations with high memory requirements (>4 GB) would cost an amount of CPU hours given by the formula The computation was memory bound, therefore causing the overall computational resources to be limited. We used the Abel cluster exclusively for the samples on the finest MLMC level.
We used the FEniCS FEM software. 65 The linear systems of equations were solved using the PETSc 66 implementation of the GMRES algorithm preconditioned with the BoomerAMG algebraic multigrid algorithm from Hypre. 67 We use the MLMC and QMC routines from the open-source software FEMLMC, 34 which contains the implementation of the algorithms presented in Refs. 18,19. For the FEMLMC Matérn field sampling solver, we declare convergence when the absolute size of the preconditioned residual norm is below a tolerance of 10 −8 . We use the same stopping criterion for the GMRES solver.
Convergence of the numerical solver was verified with a convergence test comparing different mesh refinements and time steps for a set of deterministic worst-case models (with large velocities and small diffusion coefficients), see, for example, Ref. 20 (Supplementary Material).

| QMC and MLMC algorithms
We now describe the QMC and MLMC algorithms used more in details. We adopt the following (standard) randomized QMC algorithm 35 : QMC algorithm 1. Set the required tolerance ε, θ ∈ (0, 1). Set the mesh size h and Δt to ensure that the bias estimate is lower than ffiffi ffi θ p ε; 2. Get an initial estimate ofV QMC =  I N ω ð Þ ½ (cf. (6)) with N = 1 samples and M = 32 randomizations; The total QMC cost is given by C tot = CNM, where N is the final number of samples taken and C is the (expected) computational cost of computing one QMC sample.

Remark 5.
In the random PDE case, typically the bias stems directly from the deterministic PDE solver (the FEM error in our case) 68,69 and as such it can be estimated empirically, that is, by experimenting with the mesh size and time step, and/or by using theoretical error estimates. 70 Here, we estimate the bias empirically (cf. description of the MLMC algorithm) and we check that the empirical bias convergence order with h and Δt is second-order as expected given our discretization.
For MLMC, we first need to estimate the optimal number of samples on each level for a given tolerance ε. Let C ℓ , V ℓ be the cost and variance of one sample Q ℓ (Á) − Q ℓ − 1 (Á) respectively. The total MLMC cost and variance are It is possible to minimize the estimator variance for a fixed total cost. 26 For a fixed MSE tolerance ε 2 , the optimal number of samples for each level and related total cost are, We can now describe the MLMC algorithm used in this work: MLMC algorithm (taken from Ref. 33) • Set the required tolerance ε, θ ∈ (0, 1), the maximum level L max , the initial number of levels L and the initial number of samples N ℓ to be taken on each.
• while extra samples need to be evaluated (9ℓ : N ℓ > 0): 1. for each level, evaluate all the samples scheduled to be taken; 2. compute/update estimates for the level variance V ℓ , ℓ = 1, …, L; 3. compute optimal number of samples N ℓ by using (38) and update the numbers of extra samples needed N ℓ accordingly; 4. check whether the weak error (ie, the bias) j P −P Â Ã j is below the required tolerance ffiffi ffi θ p ε; to see how the bias is estimated in practice without knowing the exact solution, see section 3 in Ref. 33. 5. if not: if L = L max report failed convergence; otherwise set L ≔ L + 1, update N ℓ and N ℓ and compute N L = N L (again using (38)).
Remark 6. The standard way 26,71 to compare the efficiency of these algorithms is to fix the same tolerance for both methods and compare their total costs, given by C tot = CNM for QMC and by C tot = P L ℓ = 1 N ℓ C ℓ for MLMC (cf. (37)). Furthermore, to ensure that the bias level is the same for both methods, the QMC routine must be run on the finest mesh and time step size of the MLMC hierarchy. This gives C = C L . In practice, for the sake of comparing methods, the actual costs C ℓ , ℓ = 1, …, L can be replaced by pseudo-costs, i.e. by setting C ℓ ≈ĉ 3 2γ ℓ , whereĉ 3 andγ are the values, estimated with CPU timings, of the constants c 3 and γ appearing in Theorem 2.3. We use this latter approach.

| NUMERICAL RESULTS
We now compare the efficiency of standard MC, QMC and MLMC when employed to solve Models 1 and 2. In what follows, we let T = 30kmin f g k = 48 k = 1 and define to be the set of all the output functionals of interest considered (cf (21) and (22)).

| Estimation of MLMC parameters
We first estimate the MLMC parameters α, β and γ of Theorem 1. Since we are considering the estimation of multiple output functionals, we estimate α and β by monitoring the bias and variance at each level ℓ: We expect α = 2, β = 2α = 4 and γ = 4 since for random PDEs we have that the bias convergence is typically the same as the deterministic FEM convergence order, β = 2α , 68,69 the numerical method is second-order in both time and space and we are using a multigrid-preconditioned Krylov method (cf. Section 5). The value for γ stems from the fact that the number of time steps on level ℓ are proportional to 2 ℓ and the linear solver used has (essentially) linear complexity in the number of degrees of freedom, which in turn scale proportionally to 2 3ℓ . We therefore get a cost per level proportional to 2 (3 + 1)ℓ and a γ = 4.
To estimate the bias and variance in (40) we take N = 4000 samples on the first two levels and N = 100 on the finest level. The choice of the number of samples on the finest level is motivated by the following considerations. The computational resources available were limited, cf. Section 5.4. The number of vertices on the finest level of the MLMC hierarchy is large (cf. Table 1), resulting in a memory burden of ≈50 GB to load the mesh, the box mesh in which the brain mesh is embedded (cf. Section 2) and the FEM subspaces. Additionally, solving one instance of (15) on this mesh takes more than 24 hours in serial. The overall sampling procedure took around 2 weeks of combined use of the Oxford Mathematical Institute servers and the University of Oslo Abel cluster. The same computation performed in serial would have taken roughly 2 years. Adding a finer level to the hierarchy would have therefore been prohibitive. Figures 3 and 4 show the bias and variance versus refinement level for Models 1 and 2, respectively. We observe that the numerical estimates closely match the theoretical expectations. The estimated variance convergence order for both models isβ≈4:2, which is just above the theoretical value of β = 4. For Model 2, the estimated bias order isα = 2:09 which again closely matches the theoretical estimate of α = 2. However for Model 1, we observe that the bias decays F I G U R E 3 Convergence behavior of the FEM approximation to the solution of Model 1. The estimated convergence order for the variance agrees with our predictions and with what expected by the theory in the diffusion-only case. 69 The bias convergence order observed is instead higher than expected. Estimated parameters via linear regression:α≈4:16,β≈4:15 F I G U R E 4 Convergence behavior of the FEM approximation to the solution of Model 2. The estimated convergence orders agree with our predictions and with what expected by the theory in the diffusion-only case. 69 Estimated parameters via linear regression:α≈2:09, β≈4:18 more rapidly than expected (α = 4:16 versus α = 2). In this case, we are likely observing a pre-asymptotic regime and the convergence rate seems to be decaying as ℓ increases; this notion would be in agreement with theoretical convergence results. 69,70 Finally, when estimating γ as taking the average wall-clock time of each sample as a proxy (data not shown), we obtainγ≈4:09, which is close to the theoretical prediction of γ = 4.

| Mean square error weighting under limited computational resources
Note that we have a finite number of meshes available and consequently the version of MLMC considered here is "weaker" than the true MLMC algorithm, described in Section 5.5, since the maximum level L is bounded, cf. Remark 3. In fact, we are unable to reduce the bias below the threshold given by the finest mesh of the hierarchy without resorting to more advanced techniques (see later Remark 8). However, we can still balance the relative weight of bias and statistical error by choosing two different values of the MLMC weight parameter θ (cf. Ref. 38 and Section 2). Based on the fact that for this problem: (a) the maximum L is restricted; (b) N L is restricted; (c) we have estimated values for the MLMC parameters α, β, and γ (cf. Theorem 1); we can estimate a priori the largest possible values of θ that we can use without making the number of samples on the finest level exceed 100. This yields the values θ = 0.041 for Model 1 and θ = 0.72 for 2, corresponding to the best achievable MSE tolerances of roughly ε = 3.3 × 10 −4 for Model 1 and of ε = 2.5 × 10 −3 for Model 2. Note that in the Model 1 case, the bias is much smaller (compare Figures 3 and 4), hence why the chosen θ is smaller as well.
Remark 7. Note that the bias (ie, the discretization error) can be estimated by the MLMC algorithm without requiring a reference or exact solution (recall the MLMC algorithm in Section 5.5) and is the same independently from which Monte Carlo method is used. Furthermore, all methods can compute their own estimates of the statistical error. For this reason, no reference or exact solution is needed to compute the MSE or to construct the figures shown in this section.
Remark 8. In practice, it is possible to reduce the MLMC estimator bias by augmenting MLMC with Richardson-Romberg extrapolation. 28,39 However, we leave this enhancement for future research.

| Tracer evolution over time
In Figures 5 and 6 we present the expected value and SD of all quantities of interest as a function of time for both models. These results have been obtained with MLMC by setting the lowest achievable MSE tolerances given the restriction on the maximum number of samples on the finest level, cf. Section 6.2. The total number of samples used for these simulations are N ℓ = [26060, 1372, 100] for Model 1 and N ℓ = [26567, 1171, 100] for Model 2 for ℓ = 1,2,3. We remark that the figure errorbars represent one SD away from the mean, and not necessarily the actual variability of these quantities. In both models, tracer first spreads to the gray matter, reaching a peak amount of tracer in 8-12 hours ( Figures 5A  and 6A). Over time, the tracer also spreads into the white matter ( Figures 5B and 6(B)). The time-to-peak amount of tracer in white matter is estimated at ≈ 24 hours in Model 2 ( Figure 6B), while for Model 1 this peak would occur considerably later ( Figure 5B). In Figure 5, we observe that for Model 1, the variances of the global amounts of tracer Q g and Q w are small. However, the variability in the local tracer concentration q g and q w is much larger and, as such, substantially contributes to the computational complexity of the Monte Carlo simulations. We observe that for Model 1, the local variability of the random coefficients "averages out" when computing global output functionals. Thus, if global quantities were the only focus, a MC simulation with small sample size could be sufficient to achieve good accuracy for this model. Conversely, we note that all quantities of interest for Model 2 (amount of tracer in gray and white matter, and tracer concentrations in local regions) have large SDs ( Figure 6).
Comparing these results with the MSE tolerances imposed we can estimate the number of correct digits of the computed expectations: the estimates obtained for Model 1 are roughly within 3 digits of accuracy for the gray matter quantities Q g and q g , 2 digits for Q w and 1 digit only for q w . Similarly, all Model 2 gray matter quantities Q g and q g are roughly correct within 2 digits, while the white matter outputs Q w and q w are just slightly less accurate.

| Comparison of MC, QMC, and MLMC performance
In Figures 7 and 8 (left) we show the optimal number of samples chosen automatically by MLMC on each level as the root mean square error tolerance ε is reduced. The maximum level chosen is increased as ε decreases in order to satisfy F I G U R E 5 Mean and SD vs time for all quantities of interest for Model 1. The continuous lines correspond to expected values and the vertical bars indicate plus or minus one SD away from the mean (negative values are excluded here). Both were computed using MLMC with a MSE tolerance of ε = 3.3 × 10 −4 , which is the best achievable tolerance given the discretization error and the restrictions on the number of samples on the finest level, cf. Section 6.2 F I G U R E 6 Mean and SD vs time for all quantities of interest for Model 2. The continuous lines correspond to expected values and the vertical bars indicate plus or minus one SD away from the mean (negative values are excluded here). Both were computed using MLMC with a MSE tolerance of ε = 2.5 × 10 −3 , which is the best achievable tolerance given the discretization error and the restrictions on the number of samples on the finest level, cf. Section 6.2 the bias tolerance. Note that the smallest values of ε considered correspond to the lowest bias tolerance that standard MLMC can achieve with an upper limit of 100 samples on the finest level (cf. Remark 8).
In Figures 7 and 8 (right), we compare the total computational cost C tot of standard MC, MLMC, and QMC for the solution of Model 1 and 2, respectively. The cost of running a full standard MC simulation for this problem is prohibitive (on the order of 100 years in serial!) so we estimate the MC cost by using the level variance and cost. This is a standard procedure and essentially derives from the application of the central limit theorem. 26 We find that QMC is significantly slower than MLMC. For this reason we only estimate the number of QMC samples needed by running QMC on the second finest mesh of the hierarchy rather than on the finest (cf. Remark 6). The difference should be minimal since the number of samples needed appears to be approximately constant on the finer levels of the hierarchy in numerical experiments. 29 Since β = γ, Theorem 1 predicts an overall MLMC complexity of ε −2 (logε) 2 for a root mean square error tolerance ε. We therefore expect a near constant cost curve for ε 2 (logε) −2 C tot versus ε in the MLMC case. The numerical results corroborate the theoretical expectations: while the MLMC cost lines oscillate some, they are F I G U R E 7 Convergence of standard MC, QMC and MLMC for the solution of Model 1 (θ = 0.041). In the plot on the left we show how the MLMC algorithm automatically selects the optimal number of samples N ℓ on each level to achieve a given tolerance ε. In the plot on the right we compare the efficiency of the methods for different tolerances. The savings of MLMC with respect to standard MC and QMC are considerable, while QMC barely improves on standard MC (see main text) F I G U R E 8 Convergence of standard MC, QMC and MLMC for the solution of Model 2 (θ = 0.72). In the plot on the left we show how the MLMC algorithm automatically selects the optimal number of samples N ℓ on each level to achieve a given tolerance ε. In the plot on the right we compare the efficiency of the methods for different tolerances. MLMC significantly outperforms QMC, which in turn considerably outperforms standard MC well-fitted by a horizontal line (estimated slope ≈0.05 for Model 1 and ≈0.1 for Model 2). Overall, for both models MLMC significantly outperforms both QMC and standard MC, with a O(100) factor of improvement with respect to standard MC.
While the qualitative behavior of standard MC and MLMC is consistent between the two models, QMC behaves differently. For Model 2 (Figure 8), the performance of QMC considerably improves on that of standard MC, especially for smaller MSE tolerances. On the other hand, for Model 1 (Figure 7), the improvement is negligible and QMC performs essentially the same as standard MC.
This behavior could be interpreted in the context of the formulations of Models 1 and 2 (cf. Section 4). While the stochastic input in Model 2 includes 1 random field and 1 random variable, Model 1 depends on 4 random fields and 1 random variable. Given the lack of performance gain for QMC applied to Model 1, we hypothesize that the higher input dimensionality affects the QMC convergence, causing the rate to decay to a standard O(N −1/2 ) MC rate. Indeed, the fact that QMC performance degrades with high input dimensions is well-known. 25 It therefore appears that the (ML)QMC method presented in Ref. 18 is not robust with respect to the number of random field inputs, at least in 3D where the dimensionality is larger. We did not observe this ill-behavior in analogous numerical tests performed on a convectiondiffusion PDE with random coefficients on a square domain.
Remark 9. Adding a coarser level to the mesh hierarchy, given by the original Colin27 human adult brain atlas mesh 53 (version 2) did not improve the performance of MLMC. This was determined by first estimating the cost and variance of the coarsest MLMC levels in the extended hierarchy using 4000 MC pilot samples and then determining which hierarchy would give the best savings in terms of cost (see eg, section 2.6 in Ref. 33).

| DISCUSSION
Most theoretical UQ research on equations with random coefficients focuses on the same model problem, 45,69,72,73 with numerical results often performed on simple boxes and in low dimensions. For this reason, and because of their different requirements and convergence properties, it is a priori unknown which method performs best for a given problem. While asymptotic complexity results might be available, 26,69,73 it might not be possible to observe the asymptotic behavior due to practical limitations that arise in large-scale problems (eg, limited computational resources, small hierarchy size, poor low-fidelity models). In fact, when considering large-scale applications with physics-inspired problems on MRI-derived geometries, currently, the only way to compare the performance of different methods is via numerical experiments. Nevertheless, the application of UQ to realistic large-scale problems is scarce in the literature, to the extent that (that we know of) this is the first study employing advanced UQ techniques to study the human brain.
The only related study that we are aware of comes from our previous work on the subject. 17 In this study we found that diffusion alone could not explain MRI data, suggesting flow velocities may impact the distribution of tracer seen in MRI data, possibly also affecting both brain clearance and drug delivery. However, in Ref. 20 standard MC was used, which forced us to resort to coarser meshes and shorter simulation time. The need for a more efficient UQ methodology motivated the current study.
In terms of limitations, in this study we considered a single brain geometry deriving from the Colin27 mesh. 74 A set of different geometries should be considered to ensure that the results obtained are patient-independent. We only modeled a single fluid flow compartment, while future research could consider more elaborate models accounting for all the different fluids flowing in the parenchyma (venous, arterial and capillary blood, interstitial fluid, cerebrospinal fluid). The computational resources available for this study were limited and with more computational power we could have obtained more accurate results by adding a finer level to the hierarchy and/or increasing the number of samples. However, we did not encounter any significant issues with setting up a hierarchy on which good MLMC convergence could be observed, a problem mentioned in recent work by Quaglino et al, 24 in which they apply MFMC for a cardiac electrophysiology application.
We also point at a limitation related to the choice of covariance structure for the diffusion and velocity fields. Here, these fields have been assigned an isotropic covariance. As such, the correlation between any two points only depends on their Euclidean distance, and not on the (domain) interior path length between the two points. Consequently, brain surface regions such as the different gyri that are closer to each other than the field correlation length, become correlated even if separated by a sulcus. Since the average sulcal width is estimated to range between 1.5 and 3 mm in humans, 44,74 this aspect does not play a role for the velocity field (λ ≈ 1 mm), but for the diffusion coefficient it may (λ = 1 cm). We expect global quantities (eg, Q g and Q w ) and local quantities depending on white matter regions which are far from the brain surface (eg, q w ) to be largely unaffected by this modeling choice. Conversely, local quantities that depend on gray matter local regions (eg, q g ) might be affected and further investigation would be needed. We remark, however, that isotropic covariances make for very efficient sampling methods, and designing and efficiently sampling a random field which is topology-aware at the high resolution scale of the sulci is highly challenging.
Overall, the random fields parameters were chosen to represent the current understanding of the relevant physiology to the extent possible, with physical considerations as the primary factor in the choice of fields and parameters. While most of the model parameters do not affect the performance of the numerical solver, QMC or MLMC in any way, the following aspects should be considered for alternative models: • Velocity magnitude and diffusion coefficient: a much larger velocity magnitude (or a much smaller diffusion coefficient) would increase the Péclet number and harm the stability of the finite element solver used (cf. Section 5.3). A more advanced numerical solver would be needed in this case. We remark that for random PDEs the convergence of MLMC directly depends on the convergence properties of the deterministic solver used. 68,69 • Random field variance: larger variability implies that more samples would be needed for an accurate result. This does not affect the efficacy of the methods, but it does increase the CPU hours needed for the computation, which might be a problem when the computational resources are limited.
• Smoothness parameter and correlation length: small values of these parameters are known to affect most random field sampling techniques and extra care must be taken. At the same time small values would also introduce sharp and abrupt spatial changes in the sampled field that are acceptable, for example, oil reservoir modeling, 75 but would be nonphysical in a healthy brain. We remark that the velocity correlation length chosen in this study is quite small, showing that the sampling techniques in Refs. 18, 19 still work well in this case.

| CONCLUSIONS
We have compared and evaluated the MC, QMC, and MLMC uncertainty quantification and sampling methods for two physiologically realistic brain transport models. Even under restriction on the maximum number of samples on the finest mesh level and on the finest mesh available, MLMC significantly outperforms all other methods, yielding an improvement factor of roughly 100 with respect to standard MC. QMC outperforms standard MC by a factor of approximately 10 for one of the models (Model 2), but yields no performance gains for the other (Model 1). Overall, for this application, MLMC achieves the best performance and should be preferred. We remark that a standard MC simulation would have taken on the order of 10 years to complete in parallel (O(100) years in serial!), while only 2 weeks are needed for MLMC.
MLQMC methods were not investigated directly in this study. However, our numerical findings indicate that no additional improvement can be achieved in the Model 1 case. On the other hand, Model 2 seems amenable to MLQMC. Thus, MLQMC could also bring additional computational gains for similar models with a small number of input random fields.
In conclusion, even large-scale problems on MRI-derived geometries are amenable to advanced Monte Carlo methods and we hope that our work can inspire further research on uncertainty quantification for biomedical engineering problems.