Multilevel and quasi-Monte Carlo methods for uncertainty quantification in particle travel times through random heterogeneous porous media

In this study, we apply four Monte Carlo simulation methods, namely, Monte Carlo, quasi-Monte Carlo, multilevel Monte Carlo and multilevel quasi-Monte Carlo to the problem of uncertainty quantification in the estimation of the average travel time during the transport of particles through random heterogeneous porous media. We apply the four methodologies to a model problem where the only input parameter, the hydraulic conductivity, is modelled as a log-Gaussian random field by using direct Karhunen–Loéve decompositions. The random terms in such expansions represent the coefficients in the equations. Numerical calculations demonstrating the effectiveness of each of the methods are presented. A comparison of the computational cost incurred by each of the methods for three different tolerances is provided. The accuracy of the approaches is quantified via the mean square error.


Introduction
The Monte Carlo (MC) method is a widely used and effective approach for uncertainty quantification (UQ) in systems of ordinary/partial differential equations (ODEs/PDEs) with random coefficients [1,2]. The implementation of this method is straightforward, it can be applied to any type of problem including nonlinear problems, it is possible to compute an estimate of the error as part of the solution process and it 2017 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.
where h (m) denotes the pressure head, K (m s −1 ) the hydraulic conductivity, q (m 2 s −1 ) the Darcy flux and g represents the source terms.
The process considered in this study is the flow of an incompressible liquid in a horizontal confined aquifer. For this problem, we consider a square flow domain R = [0, 1] × [0, 1] ⊂ R 2 , and the source terms are set to zero for simplicity, i.e. g ≡ 0. The governing equations defined in (2.1) are coupled to yield a single equation for the pressure head: ∇ · (K(x)∇h(x)) = 0, x = (x, y) ∈ R. (2. 2) The QoI to be considered in this problem is the travel time τ that a convected particle released at the centre of the domain R takes to reach the boundary of the domain, ∂R, i.e. from the point (x 0 , y 0 ) = (1/2, 1/2) to ( To compute the travel time τ , we let x = ζ (t) = (ζ 1 (t), ζ 2 (t)) be the location of a particle. After the pressure is calculated from (2.2), the trajectory ζ (t) is computed by solving the transport equation (2.4) subject to the initial condition ζ (0) = ( 1 2 , 1/2). We then determine the time τ for which ζ 1 (τ ) = 1, i.e. the convected particle lies on the right boundary, by solving [11,12] dζ (t) dt = − K(ζ ) φ ∇h(ζ ), (2.4) where φ is the rock porosity (dimensionless), i.e. the ratio of void volume in a rock to total volume. To solve equation (2.2) in R, we used a numerical code based on the standard cell-centred finitevolume method. After the pressure field h is computed, for simplicity, the spatial gradient of heads is approximated by using the central finite difference (h i+1,j − h i,j )/|x i+1,j − x i,j |, where x i,j denotes the centroid of each cell in the computational mesh (see [2] for full details). Equation (2.4) was solved by direct Euler integration. In this application, the uncertain inputs for the simulator will be the values of the hydraulic conductivity K at each of the nodes of the computational domain. The simulator output will be the travel time. It is common in groundwater flow studies [13][14][15][16] to model K as a log-Gaussian random field, i.e. to replace the conductivity by a scalar valued field, K(x), whose log is Gaussian, Z(x) := logK(x) or K(x) = exp(Z(x)). By doing this, we also guarantee that K > 0 in R. Several studies [17][18][19] have shown that although the conductivity values can exhibit large spatial variations, these are spatially correlated. A correlation function that has been extensively used [2,11,19,20] for modelling the correlation of Z(x) is the following exponential covariance function: (2.5) where λ denotes the correlation length and σ 2 is the process variance. In groundwater flow applications, the geostatistical/variogram parameters, in this case, λ and σ 2 , must be chosen according to the geostatistics of the considered porous medium. In this work, the parameters representing the conductivity fields have been selected from ranges gathered from the literature. Values around 0.3 for λ and around 1.0 for σ 2 appear to be the preferred in similar studies (e.g. [14][15][16]). Thus, in this paper, we will use λ = 0.3 and σ 2 = 1.0. Note also that an appropriate discretization scheme for this type of models must be designed according to the value of λ; in other words, the size of the computational domain has to be chosen significantly larger than the value of λ and also allow λ to be large enough to be taken into account in the numerical formulation [2], i.e. in our case, larger than the distance between centroids in adjacent cells.
To generate samples of K(x) at the nodes of the computational domain, first, we need to generate samples of Gaussian field Z(x) at such nodes. One of the most popular methods to generate different (Gaussian distributed) Z(x) is the Karhunen-Loéve (KL) expansion method [2,[13][14][15][16]21,22]. This method provides an approximation (due to the truncation of an infinite series) of the permeability fields at all the points of the continuous domain, which can be sampled afterwards on any grid. In order to avoid adding extra errors (arisen from the truncation of the KL expansion) to the model and produce more accurate representations of the hydraulic conductivity, alternative methods might be considered, for instance, the circulant embedding algorithm [23][24][25]. The circulant embedding method provides fast and exact representations of the Gaussian field but requires the use of the fast Fourier transform method, and thus, it is not straightforward to implement. Two alternatives to the circulant embedding method for producing exact decompositions of the covariance matrix associated to the correlation function given in (2.5) are the Cholesky method [11,25,26] and the KL decomposition [22,27]. These methods are not recommended for covariance functions that are differentiable at zero lag distance, e.g. the square exponential (or Gaussian) correlation function [22,28]. In those cases, the associated covariance matrix is likely to become extremely ill-conditioned [29,30]. They could be also inappropriate for problems in which the simulator necessitates an extremely fine discretization of the computational domain [30], but this does not apply to the problem considered in this paper. Conversely, the main advantages of this approach is that it only requires a single eigen-decomposition of the covariance matrix, the results of which are stored and used to generate new realizations of the permeability field very cheaply, and furthermore, its implementation is very simple and straightforward. In this paper, we will opt for the KL decomposition method, which is described briefly next (for full details, e.g. [22,25]).
Let {x j } M j=1 ⊂ R be the set of nodes for a given discretization of the problem domain R. To generate samples of Z(x), we let C be the positive semi-definite covariance matrix associated to the function c, This matrix admits an eigen-decomposition [26]: Λ is the M × M diagonal matrix of ordered decreasing eigenvalues λ 1 ≥ λ 2 ≥ · · · ≥ λ M ≥ 0, and Φ is the M × M matrix whose columns φ i , i = 1, . . . , M, are the eigenvectors of C. Let ξ i ∼ N (0, 1), i = 1, . . . , M, be independent and identically distributed (i.i.d.) random variables. We can draw samples from Z ∼ N (m, C) at the points x j using the KL decomposition of Z := (Z 1 , . . . , Z M ) using the following [25]: Without loss of generality, we will set m ≡ 0 in (2.6), and thus, the discrete random permeability field is therefore given by The terms ξ i ∼ N (0, 1) above will be called KL coefficients. Now, for each new ensemble {ξ j 1 , . . . , ξ j M }, j ∈ N, of random variables ξ j i ∼ N (0, 1), we can generate a new realization of the conductivity K j ∈ R M . Note that this method only provides values of the conductivity at the nodes and not in the whole continuum R.
In the following section, we include a review of the literature related to the implementation of the four MC methods.

Monte Carlo simulation methods
Let (Ω, F , P) be a probability space. Let X M := (ξ 1 , . . . , ξ M ) be the random vector formed with a given ensemble of M KL coefficients which yields to a discrete random permeability field K ∈ R M . Let T M = f (X M ) ∈ R, where f denotes the travel time simulator, be the approximation of the travel time obtained with the simulator based on a computational domain of M nodes {x j } M j=1 . We denote by T the true (underlying) travel time random variable T : Ω → R solution of (2.2), and assume that the expected value , as M → ∞, and that (in mean) the order of convergence is α > 0 (see [2,3] for further details), i.e.

|E[T
for some constant C α . We are interested in estimating E[T]. Thus, given M ∈ N sufficiently large, we compute approximations (or estimators)T M of E[T M ] and quantify the accuracy of our approximations via the root mean square error (RMSE): We will denote by C ε the computational ε-Cost used to achieve an RMSE e(T M ) ≤ ε. This ε-Cost is quantified by the number of floating point operations that are needed to achieve an RMSE of e(T M ) ≤ ε. As it could be well known to the reader (but it is important to remark here before discussing the MLMC method), when solving a system of PDEs, a computer model needs to retain all the important features of the physical domain (a continuum medium) of the problem and reduce them into a simplified form, called the computational domain (a discrete set of points). Throughout this paper, we will use the term grid for the structured distribution of points, called nodes, that form the computational domain used by the computer model to solve the equations, and M will denote the number of nodes which form the corresponding grid. According to this, given two grids M i and M j with i < j, i, j ∈ N, we will say that M i is a subgrid of M j , and we will write M i < M j , if all the nodes contained in M i are also contained in M j . We will then say that M i is coarser than M j and conversely that M j is finer than M i . For solving efficiently a system of PDEs, choosing M sufficiently large corresponds to choosing a fine enough grid that guarantees that the computer model is providing an accurate approximation of the true solution of the problem. Figure 1 shows an example of two grids for the same physical domain used by a computer model.
In the following sections, we describe how to implement each of the MC methods. Note that while for all of the methods the QoI is E(T M ), in each of the methods we use a different estimator to approximate E(T M ).

Classical Monte Carlo simulation method
We define the standard MC estimator for estimating E(T M ) as follows: where T We assume that the cost to compute one sample T and hence the total cost of the MC estimator satisfies [2] The MSE ofT MC M,N can be expressed as follows [2]: where V[T M ] is the variance of the MC estimator, which represents the sampling error and decays inversely with the number of samples. The second term on the right-hand side is the square of the error in mean between T M and T, also called the discretization error or the bias.  Two samples of the same random permeability field in two consecutive levels (a) and + 1 (b) to be used as input in the MLMC method. In this example, we used = 3.

Multilevel Monte Carlo simulation
The main idea behind the MLMC simulation method is to start obtaining approximations of T from several grids, starting by the coarsest and stopping when the given MSE has been numerically achieved. For a detailed description of the method, the reader is referred to [2,3]. In this section, we only give a brief summary of the practicality of the approach.
, and so the MSE is We see that the MSE for the multilevel estimator consists of two terms, the variance of the estimator and the approximation error. To bound the RMSE by ε, we can seek to bound each term above by ε 2 /2. Note that the second term is exactly the same as in equation ( Thus, to then achieve an overall RMSE of ε, the first term of e(T ML M ) 2 is also bounded by ε 2 /2. The computational cost of the MLMC estimator is [2]: where C := C(Y (i) ) represents the cost of a single sample of Y .
The variance of the MLMC estimator can be minimized [2] for a fixed computational cost by choosing with the constant of proportionality chosen so that the overall variance is ε 2 /2. So, the total cost on level is proportional to V[Y ]C and hence we can write [31]: In practice, optimal values for L and {N } L =0 can be computed from sample averages and the unbiased sample variances of Y . If we assume that sufficiently large, providing us with a computable error estimator to determine either whether M is sufficiently large or whether the number of levels L needs to be increased.
The above conditions and statements are formally presented in the following theorem [2]: Then, for any ε < e −1 , there exist a positive constant C ML , a value L (and corresponding M ≡ M L ) and a Proof. The proof is given in [2].
The MLMC algorithm can be implemented in practice as follows: (iii) Calculate the optimal N , = 0, . . . , L, using (3.12), Remember that C := C(Y (i) ) represents the cost of a single sample of Y . This step aims to make the variance of the MLMC estimator (3.9) less than 1 2 ε 2 . (iv) Evaluate extra samples at each level as needed for the new N .
If not converged, set L = L + 1 and go back to 2. The parameters, α, β and γ that can be estimated empirically as follows: For γ , we assume that the number of operations to compute one sample on level is C = cM γ for some constant c independent of . For β, we can use as an approximation the slope of the line for log For α, we can use as an approximation the slope of the line for log

Quasi-Monte Carlo simulation
Many of the existing variance reduction methods built upon pseudo-random sequences, e.g. MLMC, are focused on reducing the overall computational cost of a numerical simulation. QMC methods aim to accelerate the rate of convergence of MC by using deterministic (also called quasi-random or low-discrepancy) sequences instead of pseudo-random. The discrepancy of a sequence is a measure of its uniformity and it is computed by comparing the actual number of sample points in a given volume of multidimensional space with the number of sample points that should be there assuming a uniform distribution. These methods normally achieve a convergence rate of O((log N) M /N). Hence, the convergence rate is directly related to the dimension of the space. This dependence on the dimension of the space together with the correlation of the points generated by the QMC method yields sometimes non-accurate and biased results. That is the main reason why, during the past two decades, QMC methods have been mostly applied to models defined over low-dimensional probability spaces [8,32,33].
In recent years, there has been an increasing interest in tackling the problem of UQ in models of physical processes, for instance, transport in porous media or carbon capture and storage in deep saline aquifers. As discussed in §2, the uncertainty in those models is often represented by truncated KL expansions of log-Gaussian random fields defined in high-dimensional probability spaces. The truncation of these KL expansions adds more uncertainty to the model and this affects the accuracy of the results. Although QMC methods have already been successfully applied to problems defined in high-dimensional spaces by employing different representations of the random inputs [34][35][36][37], to the best of our knowledge they have not yet been used in models represented by direct KL decompositions. In this section, we apply the QMC method to an extremely high-dimensional problem with log-Gaussian distributed inputs and present numerical evidence of the acceleration of the MC rate of convergence. Before introducing the QMC simulation method, let us describe in more detail the MC integration procedure. The MC method uses pseudo-random number sampling algorithms, i.e. during the generation process, uniformly distributed pseudo-random numbers are generated and transformed into the KL coefficients which jointly form random input vectors in R M , and these are distributed according to a certain probability distribution, in our case, N (0, I). Let us see an illustrative example: let (Ω, F , P) be a probability space and let g : [0, 1] M → R, and Y = g(Z), where Z is an uniformly distributed random vector in [0, 1] M . Suppose that we wish to compute I = [0,1] M g(ξ ) dξ with the MC method. Let p denote the uniform probability density function and letting ξ be uniformly distributed in [0, 1] M , we can apply MC quadrature to approximate I, for a given N ∈ N, in the following way: where ω j ∈ Ω and the values ξ (ω j ) ∈ R M are i.i.d. random vectors sampled uniformly by sampling the components ξ i (ω j ) independently and uniformly on the interval [0, 1]. Some examples of quasi-random sequences are: digital nets [38], rank-1 lattice rule [39], Faure sequences [40] or Sobol sequences [41]. From a deep review of the literature, Sobol sequences seem to be the most popular for being used by the QMC method in mathematical models with random inputs [4,6,8], and thus, we will opt for Sobol sequences in this paper. The biggest difference to pseudo-random sequences is that the sample values are chosen under consideration of the previously sampled points, thus avoiding the occurrence of spatial clusters and gaps, as we can observe in figure 3. Figure 3a shows 100 pseudo-random numbers sampled from a uniform distribution in the unit square. Figure 3b shows the same number of points generated by using a Sobol sequence. It can be observed that the sampling space is filled in a more uniform manner in figure 3b. Figure 3c,d show, respectively, the spatial distribution of 2000 points with pseudo-randon numbers generation and Sobol sequences.

Multilevel Quasi-Monte Carlo simulation
Although there are currently many researchers using MLQMC, there are still very limited works (most of them still in press) in the literature [39,42,43]. Thus, to the best of our knowledge, the application of MLQMC to the case of direct KL decompositions for log-Gaussian random fields is also new. This method is a consequence of combining the MLMC algorithm with a randomized QMC estimator instead of the MC estimator. In this paper, we use the MLQMC algorithm developed by Giles & Waterhouse [5]. In order to obtain unbiased estimators for the variances, we need to induce some randomness to the QMC points, this process is known as QMC randomization. There are several ways of QMC randomization, depending on the type of low-discrepancy sequence used. In this study, we use the digital scrambling technique described in [44]. This consists in building a set of n (we will use n = 16 in this study) scrambled   Sobol' sequences to obtain averages for the quantityŶ at level , i.e.Ŷ is the average of the quantitieŝ T 0 andT −T −1 (for > 0) over the 16 sets of N QMC points. The MLQMC algorithm (described in [5]) can then be summarized as follows: In the following section, numerical results from the application of the above methods to the model problem described in §2 for several discretizations of the physical domain are discussed.

Numerical results
The procedure followed for conducting the experiments is as follows: firstly, we check (empirically) from which level (i.e. the value of M 0 ) the asymptotic hypotheses of theorem 3.1 are satisfied (this assures that the simulations at the coarsest grid are reliable approximations of the QoI); secondly, we set the tolerance (MSE) for which we wish the MC algorithms to stop; thirdly, we use the conclusions of theorem 3.1 to implement the four methods as discussed earlier in §3; and finally, the performance of each of the methods is tested by comparing their computational ε 2 -Cost, i.e. the number of floating point operations that are needed to achieve the given MSE.
The three tolerances employed for all the comparisons are: 0.01, 0.0064 and 0.0025. The average travel times estimated with MC, MLMC, QMC and MLQMC methods will be denoted, respectively, by T MC , T MLMC , T QMC and T MLQMC . The sequence of levels will start with M 0 = 81. This enables one to get a minimal level of resolution of the problem [2,3]. The maximum level considered will be The conditions of theorem 3.1 for the mean and the variance of the MLMC and MLQMC estimators will be numerically confirmed for each of the cases. The estimates of the parameters α and β will be computed 'on the fly' from sample averages. The dominant cost will rely on the PDE solution, and an algebraic multi-grid method is used as the iterative linear solver. Hence, γ = 1 in all the simulations. To quantify the cost of the algorithms, we will assume that the number of operations to compute one sample on level is C = cM for some fixed constant c, and thus, in the results presented in this paper, we will show the standardized costs, scaled by 1/c, i.e. the cost is defined as L =0 N M .

Comparison between classical Monte Carlo and multilevel Monte Carlo
In this section, we compare the performance of MC and MLMC methods for the tolerances above. As could be expected from similar works in the field and after reviewing the theory related to both methods, the MLMC method clearly outperforms the standard MC. The MLMC results are presented in tables 1-3. The MC results are given in table 4. Thus, by looking at tables 1-4, we observe that while the MLMC method reduces the computational cost of MC for the same degree of accuracy at a rate of 20-26 for tolerances of MSE = 0.01 and MSE = 0.0025, the reduction reaches its peak at the rate of 80 for a tolerance of MSE = 0.0064. Henceforth, in this application, MLMC performs best for a grid of M 4 = 16 641 elements. Tables 1-3 show the number of samples, N , used by the MLMC method in each level, , for the given MSE, ε 2 , the final computational ε 2 -Cost (cost for that given tolerance ε 2 ), the value of the average travel time, T MLMC , and the corresponding bounds for the estimation (T MLMC − ε, T MLMC + ε). Figure 4 shows the expected value of T and Y = T − T −1 and how the slope of the line for E[T − T −1 ] has a decreasing tendency. It also shows how E[T ] is approximately constant on all levels, numerically verifying condiction (i) of theorem 3.1. Figure 5 shows the behaviour of the variance of T and Y = T − T −1 for each level , and how the condition (ii) of theorem 3.1 is numerically satisfied on the levels shown.

Comparison between Monte Carlo and quasi-Monte Carlo
In this section, we compare the performance of MC and QMC methods for the same tolerances used in the previous sections. In this case, low-discrepancy sequences clearly outperform pseudo-random for all the tolerances. Similarly to what happened with MLMC, the reduction in cost provided by the QMC method when compared to MC reaches its peak at level 4. The reduction rate achieved at this level is 9.    This could indicate that after the discretization error has been adequately reduced, and consequently, a fine resolution of the QoI is being obtained in each simulation, there is not much additional gain by reducing the sample variance (or sampling error). The latter can be also deduced from figure 9, where after level 4 (or tolerance 0.0064) the slope of the cost for standard MC is nearly constant.  Table 4 provides the data comparison of the computational ε 2 -Cost for the MC and QMC. These quantities are obtained in the corresponding MLMC and MLQMC simulations. To calculate the costs for the MC and QMC methods, we use the estimator provided in [3]: where N * = 2ε −2 V[T ], so that the variance of the MC (3.3) and QMC (3.14) estimators is 1 2 ε 2 as with the corresponding MLMC and MLQMC methods.
In addition to this ε 2 -Cost comparison, we will analyse the convergence of the MC and QMC methods at each of the levels where the multilevel methods converged. Figure 6 shows the convergence analysis,   Table 5 gives the values of the MC and QMC estimators at each level based on N = 25 000.

Comparison between quasi-Monte Carlo and multilevel quasi-Monte Carlo
In this section, we compare the performance of QMC and MLQMC methods for the same MSEs as above. In this case, unlike in the comparison between MC and MLMC, MLQMC outperforms QMC in a monotonic order, i.e. the reduction in the cost follows an increasing rate along with the increase in the degree of accuracy (or reduction in tolerance).  the logic of deterministic sequences generation, and they seem to be (as one could expect) a direct consequence of the ordered (deterministic) way in which the MLQMC estimator is built.
We illustrate next the same tables and figures shown in the previous section for the MC and MLMC methods. Tables 6-8 give the number of samples, N , used by the MLQMC method in each level, , for the given MSE, ε 2 , the final computational ε 2 -Cost incurred by using the given tolerance, the value of the average travel time, T MLQMC , and the corresponding bounds for the estimation, (T MLQMC − ε, T MLQMC + ε). Figure 7 shows the expected value of T and Y = T − T −1 and how the slope of the line for E[T − T −1 ] has a decreasing tendency. It also shows how E[T ] is approximately constant on all levels. Figure 8 shows the behaviour of the variance of T and Y = T − T −1 for each level , and how the condition (ii) of theorem 3.1 is numerically satisfied on the levels shown.

Comparison of Monte Carlo, quasi-Monte Carlo, multilevel Monte Carlo and multilevel quasi-Monte Carlo
The overall picture with the performance of all the methods is shown in figure 9. We can see how the MLQMC method produces a lower computational cost for all the tolerances. MLMC is performing better than MC and QMC, and in conclusion, MC seems to be the least efficient method.

Conclusion and further work
In this paper, we analysed the efficiency of MC, MLMC, QMC and MLQMC in achieving a desired error level on the estimation of the average travel time during the transport of particles in heterogeneous porous media. The analysis was focused on employing the four methods to solve, under the same conditions, a stochastic model defined in a high-dimensional probability space, and in comparing the computational costs incurred by the four different approaches. The improvements were related to the use of low-discrepancy (Sobol) sequences for the space filling design (QMC) and variance reduction in the multi-grid schemes (MLMC). One conclusion that can be drawn from the review of the literature and the results obtained in this paper is that, on one hand, for 'smooth' uncertain model parameters defined in high dimensions, e.g. the log-Gaussian representation of the hydraulic conductivity in Darcy's Law, we can rely on QMC methods to significantly reduce the computational cost in an uncertainty analysis, while providing accurate results when compared with other methods like MC. On the other hand, in cases where the uncertain parameters are not smooth enough (e.g. with discontinuities), the QMC method reviewed in this paper may yield inaccurate and biased results. In this case, the use of unbiased randomized QMC estimators as the one used in the MLQMC method might be an alternative, although this would lead to a loss of the deterministic control offered by the standard QMC. A description of such randomized QMC methods is provided in [5].
We provided a detailed comparison of the accuracy and efficiency between the different methods. From the numerical results obtained in the model problem studied in this paper, the QMC and MLMC methods provided the same order of accuracy that the classical MC with considerably less computational runs. The combination of both methods led to the MLQMC method, which was proved to provide the optimal computational effort for the simulator while retaining the same accuracy in the calculations.
In terms of practicality, the multilevel schemes require additional work on the simulator's numerical code in order to carry out the corresponding multi-grid approach, and this could be impractical for users of Engineering commercial packages for instance. Although the multilevel approaches could also be used for non-nested grids, for non-uniform shapes of the computational domain, methods like the multi-index Monte Carlo [45] could be a better choice.
Further research may include testing the performance of the methods by considering alternative pseudo-random sequences to Sobol when building the QMC and MLQMC estimators, for instance, rank-1 lattice rule [39] or Faure sequences [40]. Refining the MLQMC method discussed in this paper, and therefore reducing its computational cost, is also possible by exploiting the deterministic way in which the estimation of the QoI is conducted, i.e. we could design an algorithm that returns the minimum number of samples needed at each level that makes the statistical error be lower than the given tolerance, instead of using just an exceeding estimation.
Data accessibility. Our data are deposited at: http://dx.doi.org/10.5061/dryad.ft1d5 [46]. Authors' contributions. Both authors conceived and designed the study. D.C.-G. carried out the implementations, analysis and drafted the manuscript. Both authors gave their final approval for publication.
Competing interests. We declare we have no competing interests. Funding. This research was funded by the EU Panacea project, FP7, grant agreement no. 282900. D.C.-G. also acknowledges financial support from EPSRC, grant no. EP/L027682/1.