{'fraction_non_alphanumeric': 0.038205752757327786, 'fraction_numerical': 0.00897670981690796, 'mean_word_length': 4.565031982942431, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 19, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Uncertainty and Sensitivity Analyses Methods for Agent-Based Mathematical Models: An Introductory Review

Multiscale, agent-based mathematical models of biological systems are often associated with model uncertainty and sensitivity to parameter perturbations. Here, three uncertainty and sensitivity analyses methods, that are suitable to use when working with agent-based models, are discussed. These methods are namely Consistency Analysis, Robustness Analysis and Latin Hypercube Analysis. This introductory review discusses origins, conventions, implementation and result interpretation of the aforementioned methods. Information on how to implement the discussed methods in MATLAB is included.


Introduction
Mathematical models of biological systems are abstractions of highly complex reality. It follows that parameters used in such models often are associated with some degree of uncertainty, where the uncertainty can be derived from various origins. Epistemic uncertainty refers to uncertainty resulting from limited knowledge about the biological system at hand, whilst aleatory uncertainty stems from naturally occurring stochasticity, intrinsic to biological systems [1,2,3]. Model parameters may thus be naturally stochastic, theoretically unknown, and unfeasible or impossible to measure precisely (or at all). Further magnifying the contributions of uncertainty in mathematical models of biological systems, in particular, is the fact that one parameter in the mathematical model may correspond to a multitude of underlying biological mechanisms and features in the real, biological system. This is especially true for minimal parameter models, i.e. mathematical models that aspire to be as non-complex as possible whilst still capturing all biological features of interest [4]. There exist multiple method papers that describe how to perform uncertainty and sensitivity analyses methods, authors Alden et al. even provide a free R-based software package (Spartan [2]) that enables the user to perform different such methods, including the three methods discussed in this review. However, as these methods have been developed across multiple research fields, both inside and outside of the natural sciences, it is difficult to find one comprehensive review that discusses not only how to perform these methods, but also where these methods come from, and why certain conventions are proposed and/or used. To this end, we have in this paper gathered such information for three uncertainty and sensitivity analyses techniques, namely Consistency Analysis, Robustness Analysis and Latin Hypercube Analysis. Our aim is that this will allow the reader to better evaluate uncertainty and sensitivity analyses presented by other authors, and encourage the reader to consider performing these methods when suitable.
In order to understand the impact that parameter uncertainty and parameter perturbations have on results produced by a mathematical model, uncertainty and sensitivity analyses can be used. A mathematical model that comprises a set of uncertain model parameters (or inputs), is able to produce a range of possible responses (or outputs). Uncertainty analysis assesses the range of these outputs overall, and provides information regarding how certain (or uncertain) we should be with our model results, and the conclusions that we draw from them [5]. Sensitivity analysis describes the relationship between uncertainty in inputs and uncertainty in outputs. It can be used to identify which sources of input uncertainty (i.e. which model parameters) significantly influence the uncertainty in the output and, equally importantly, which do not [5]. Assessing how sensitive the output is to small input perturbations is a healthy way to scrutinise our mathematical model [6]. Moreover, for a well-formulated model, knowledge regarding how input uncertainty influences output uncertainty can yield insight into the biological system that has not yet been empirically observed [2]. Furthermore, if the uncertainty in some input parameter is shown to not affect output uncertainty, the modeller may consider fixing that parameter, and thus reducing model complexity in accordance with a minimal-parameter modelling approach. In local sensitivity analysis techniques, model parameters (inputs) are perturbed one at a time whilst other parameters remain fixed at their calibrated value. In global sensitivity analysis techniques, all model parameters are simultaneously perturbed [7].
There exist several sensitivity and uncertainty analyses techniques, but here we will focus on three such techniques that are suitable to use in conjunction with agent-based mathematical models. These techniques are namely Consistency Analysis, Robustness Analysis and Latin Hypercube Analysis, which all answer important, and complementary, questions about mathematical models and their corresponding in silico responses [1,2].
Note that Consistency Analysis is only meaningful when analysing models with stochastic variables. Note that there exist other uncertainty and sensitivity analyses techniques, suitable for agent-based models, that are outside the scope of this review [2,8].
The statistical techniques described in this work have been developed and applied across multiple academic disciplines, both inside and outside of the natural sciences. Consequently, terminology and notations vary in the literature. The aim of this work is to combine pertinent literature from various academic fields whilst keeping terminology and mathematical notations consistent, unambiguous and tailored towards a mathematical and scientific audience. Therefore, when needed, certain algorithms from the literature are here reformulated into expressions that a mathematician would consider to be conventional. This review is intended to provide gentle, yet comprehensive, instructions to the modeller wanting to perform uncertainty and sensitivity analyses on agent-based models. Thorough directions on how to perform Consistency Analysis (Section 3), Robustness Analysis (Section 4) and Latin Hypercube Sampling and Analysis (Section 5) are provided. Consistency Analysis utilises the measure of stochastic superiority, which is therefore discussed in Section 2. Throughout this work, we have included some historical information that elucidates why certain statistical conventions are used. Each section also contains pictorial, step-by-step instructions on how to perform the aforementioned techniques. As a case study example, all methods discussed in this review are implemented in one of our recent agent-based, multiscale, mathematical oncology studies [9].
2 The measure of stochastic superiority

The Common Language Statistics
In 1992, McGraw and Wong introduced the common language statistics (CL) as an intuitive way to compare two distributions of data [10]. The CL was initially introduced as a tool to compare data from normal distributions, but was later on approximated for use on any continuous distributions. The CL describes the probability that a random data sample from one of the distributions is greater than a random data sample from the other distribution. For example, if we have two continuous data distributions B and C, and we are comparing the distributions with respect to some variable X, then the CL is simply given by where standard probability notations have been used so that P (X B > X C ) denotes the probability that a random data sample X B from distribution B is greater than a random data sample X C from distribution C [10]. Thus the subscript of X here signifies the distribution from which the data sample X was taken.

The A measure of stochastic superiority
The CL was developed to compare continuous data distributions, but Vargha and Delaney [11] introduced the A measure of stochastic superiority (or A-measure for short) as a generalisation of the CL that can directly be applied to compare both continuous and discrete distributions of variables that are at least ordinally scaled. When comparing two distributions B and C, with respect to the variable X, the A-measure A BC (X) is given by where P (X B = X C ) denotes the probability that a random data sample from distribution B is equal to a random data sample from distribution C. By comparing Equations 1 and 2, it is clear that in the continuous case, where P (X B = X C ) = 0, the A-measure reduces to the CL.
If two distributions that are identical with respect to the variable X are compared, then P (X B > X C ) = P (X C > X B ) and we say that the distributions B and C are stochastically equal with respect to the variable X. On the other hand, if P (X B > X C ) > P (X C > X B ), then we say that the distribution B is stochastically greater than distribution C, and accordingly, that distribution C is stochastically smaller than distribution B [11]. If distribution B is stochastically greater than distribution C with respect to the variable X, it simply occurs more often that the sample X B is greater than the sample X C when two random samples X B and X C are compared. Likewise, if distribution C is stochastically smaller than distribution B with respect to the variable X, it occurs more often that the sample X C is smaller than the sample X B when comparing two random samples X B and X C . These definitions of stochastic relationships (stochastically equal to, stochastically greater than, stochastically smaller than), used by Vargha and Delayney [11], amongst others, are weaker than definitions used by some other authors, but sufficient and appropriate for the our current purposes: comparing distributions of discrete data samples produced by in silico simulations based on stochastic, individual-based mathematical models.
When comparing two samples X B and X C , the possible outcomes are (i) that X B is greater than X C , (ii) that X B is equal to X C and (iii) that X B is smaller than X C . These three possible outcomes must sum up to one so that, In the continuous case, P (X B = X C ) = 0 as previously stated, and thus it follows that and thus it suffices to know only one of the values P (X B > X C ) or P (X C > X B ) in order to determine the stochastic relationship between the distributions B and C with respect to X.
For example: if P (X B > X C ) = 0.4, then it is clear that P (X C > X B ) = 0.6 and thus that P (X B > X C ) < P (X C > X B ), or equivalently, that distribution B is stochastically smaller than distribution C.
However in the discrete case, P (X B = X C ) is not generally equal to zero and therefore, Consequently, one single value P (X B > X C ) or P (X C > X B ) alone can generally not be used to determine the stochastic relationship between the distributions B and C.
For example: if, again, P (X B > X C ) = 0.4, it follows that P (X C > X B ) = 0.6 − P (X B = X C ). This does not give us enough information to determine the stochastic relationship between the two distributions B and C.
In order to proceed to compare the distributions B and C in this case, the stochastic difference δ is introduced where δ is given by Via a linear transformation, the transformed stochastic difference, δ ∈ [0, 1], can be obtained using Equation 5 so that from which we can see that the A-measure, A BC (X) (Equation 2), measures the stochastic difference between P (X B > X C ) and P (X C > X B ) under a linear transformation [11].
In order to estimate the A-measure using samples from two distributions, the point estimate of the A-measure, here denoted theÂ-measure (with a hat), is used. (In the Spartan package [2], this is referred to as the A test score). For example, if we want to compare two discrete distributions B and C, where B comprises m data samples (of some variable X) so that B = {b 1 , b 2 , .., b m } and C comprises n data samples (of some variable X) so that C = {c 1 , c 2 , .., c n } then where i = 1, 2, .., m and j = 1, 2, .., n and #(event) is the 'counting function' that simply denotes the number of time that a certain event occurs when comparing all possible pairs of data samples (b i , c j ). For clarity, Figure 1 provides an example of how theÂ-measure can be computed by simply counting events. Using more conventional mathematical notation, theÂ-measure is given bŷ where H(x) is the Heaviside step function such that IfÂ BC (X) = 0.5, then the distributions B and C are stochastically equal with respect to the variable X. TheÂ-measure can thus be used to measure 'how equal' two discrete distributions B and C are, by assessing how much theÂ-measure (∈ [0, 1]) deviates from equality, i.e. the value 0.5. The closer theÂ-measure is to 0.5, the 'more equal' the two compared distributions are [11]. In many applications, we are only interested in 'how equal' two distributions B and C are, and it is not important which distribution is the stochastically greater one. In such cases we are only interested in how much theÂ-measure deviates from stochastic equality (i.e. the value 0.5) but the direction is not important. Or in mathematical terms: the magnitude of the difference between theÂ-measure and stochastic equality is important but the sign is not. The magnitudalÂ-measure, here denotedÂ with an underscore, ignores the sign of deviation from equality and is given byÂ The statistical significance is used to describe the effect of the stochastic difference between two distributions B and C. If two distributions B and C are 'fairly equal' (i.e. if they yield an A BC -measure close to 0.5) then the statistical significance is classified as small. The statistical significance is classified using the magnitudalÂ-measure and, using guidelines from Vargha and Delaney [11], the statistical significance is classified to be small, medium or large with respect to X according to the following threshold values forÂ BC (X) , These threshold values (that might appear somewhat arbitrary) were first introduced by psychologist and statistician Cohen [12,13] in the 1960s when comparing normal distributions, but then in terms of another statistical measurement: the effect size (Cohen' and σ is the standard deviation of either B or C (as B and C here are assumed to have the same standard deviation) [13,14]. Omitting details from statistics, a small d-value essentially corresponds to a big overlap between distributions B and C, whilst a large d-value corresponds to a small overlap between distributions B and C, as is illustrated in Figure 2. Cohen decided to use the threshold d-values for describing 'small', 'medium' and 'large' effect sizes to be 0.2, 0.5 and 0.8 respectively [13]. If we hold on to the assumption that B and C are two normal distributions with the same variability, and furthermore say that they contain the same number of data samples, we can use measures of overlap to get a further 'feel' for the previously discussed effect sizes, as illustrated in Figure 2. Cohen's d value can also be converted into 'the probability that a random data sample X B from (normal) distribution B is larger than a random data sample X C from (normal) distribution C [10], but that is exactly what theÂ-measureÂ BC (X) measures! So this is where the threshold values for the descriptors 'small', 'medium' and 'large' statistical differences listed in Equation 12 come from. Now, Cohen motivated his choice of the d-value thresholds using a blend of intuitive 'everyday' examples and mathematical reasoning [13], but he did issue a warning regarding the fact that the threshold values should be determined based on the research methodology at hand. Thus the (modeller) should not blindly use Cohen's suggested thresholds, but instead reason what constitutes a small enough statistical significance in the study at hand. The (modeller) must also decide how fine the data samples in the data distributions should be before performing consistency analysis. In many applications, it is likely the amount of data samples required in order to achieve a small statistical significance increases with the fineness of the data. Nonetheless, scientific conventions are useful (no need for citations) and thus in the remainder of this chapter we will use the threshold values suggested by Cohen, as is done in other mathematical biology studies [2].

Consistency Analysis
In silico simulations based on mathematical models with built-in stochasticity will not produce the same output data every simulation run. Consistency Analysis (also called aleatory analysis) is a stochastic technique that answers the question: how many data samples do we need to produce in order to mitigate uncertainty originating from intrinsic model stochasticity? In our case, one data sample is the product of one in silico simulation, so an equivalent question is: how many in silico simulations should we run before describing our results in terms of, for example, average values, standard deviations or similar?
Let us say that one in silico simulation produces one data sample of some output variable X. This data sample can for example correspond to 'the population size at time point T ', or something similar. It is up to the modeller to identify and decide what the meaningful output variable(s) should be, and consistency analyses can be performed on multiple output variables at multiple time steps, for comprehensiveness. Before we begin, note that when performing Consistency Analysis, we always use the calibrated model parameters.
The first step involved in performing Consistency Analysis it to produce multiple distributions of data of various sizes. We say that a distribution with n data samples has a distribution-size n, and the goal of Consistency Analysis is to find the smallest n-value (here denoted n * ) that yields a small stochastic significance. To do this, we create various distribution groups that all contain 20 distributions each of some distribution-size n, as is shown in Step 1 in Section 3.1. Following the methodology described in previous work by Alden et al., and the Spartan package that they developed, [2], we create one distribution group that contains 20 distributions of size n = 1, one distribution group that contains 20 distributions of size n = 5 and so on. Here, the n-values 1, 5, 50, 100 and 300 are evaluated [2] and thus we must produce a total of 20 · (1 + 5 + 50 + 100 + 300) = 9120 in silico runs. (Note that if the desired accuracy is not achieved for the highest investigated n-value, here n = 300, higher values of n can be explored).
We here let a distribution D n,k denote the k-th distribution of distribution size n so that D n,k = {d 1 n,k , d 2 n,k , .. , d n n,k } where d h n,k is the the h-th data sample in distribution D n,k and h = 1, 2, .., n. TheÂ-measure resulting from comparing two distributions D n,k and D n,k with respect to the variable X is denoted byÂ n k,k (X). Now, within every distribution-group, we compare the first distribution (k = 1) to all other distributions (k = 2, 3, .., 20) using theÂ-measure. This yields 19Â-measures per distributiongroup (as is shown in Step 2 in Section 3.1. The maximum scaledÂ-measure with respect to X, occurring in a distribution-group g that contains distributions of size n g , is denotedÂ ng max (X). The smallest value n g for whichÂ ng max (X) ≤ 0.56 is denoted n * . In other words: n * corresponds to the smallest distribution-size for which all of the 19 computedÂ-measures yield a small stochastic significance, as is shown in Step 3 in Section 3.1. This answers the question that we set out to answer via Consistency Analysis: n * data samples (or in silico runs) are needed in order to mitigate uncertainty originating from intrinsic model stochasticity. The procedure on how to perform Consistency Analysis is outlined Section 3.1.

Quick Guide: Consistency Analysis
Here follows a quick guide for how to perform Consistency Analysis.

Robustness Analysis
Robustness Analysis answers the question: how robust are model responses to local parameter perturbations? Robustness Analysis investigates if, and how, perturbing the value of one input parameter significantly changes an output X. Using theÂ-measure, data distributions containing output data produced by perturbed input parameters, are compared to a data distribution containing output data produced by the calibrated input parameters. All perturbed data distributions are here of size n * , where n * is decided in the Consistency Analysis process, previously described in Section 3, when analysing stochastic models.
Before commencing the Robustness Analysis, we must identify the uncertain model parameters that we want to investigate the robustness of. We denote these parameters p i , where i = 1, 2, .., q, and thus we have a total of q parameters whose robustness we will investigate. Now, as illustrated in Step 1 in Section 4.1, we let each such parameter p i be investigated at r(p i ) different parameter values (including the calibrated value), and thus we need to generate a total of P distributions of sample size n * where Note that the number of investigated parameter values, r(p i ), need not be the same for every input parameter p i . Investigated distributions of sample size n * are here denoted D n * ,p i j , where i = 1, 2, .., q denotes which parameter is being perturbed and j = 1, 2, .., r(p i ) denotes the specific perturbation of parameter p i . For some perturbation j = C, the parameter value p i j equals the calibrated value for input parameter p i . For each parameter that we are investigating, theÂmeasure is used to compare the calibrated distribution D n * ,p i C to all distributions D n * ,p i j . Note that when j = C, the calibrated distribution is compared to itself and thus theÂ-measure equals 0.5. TheseÂ-measures provide information regarding the statistical significance, specifically if it can be described to be small, medium or large under parameter perturbations. Plotting the correspondingÂ-measure over the parameter value p i j for each parameter p i , paints an informative picture of local parameter robustness, as shown in Step 2, in Section 4.1. Another descriptive way to demonstrate the influence that parameter values p i j have on some output variable X is to use boxplots. As is illustrated in Step 3 in Section 4.1, boxplots can be used to clearly show the median, different percentiles, and outliers of some data distribution D n * ,p i j as a function of the parameter value p i j . The methodology to perform Robustness Analysis is outlined in Section 4.1. Note that Robustness Analysis does not pick up on any non-linear effects, between an input parameter p i and an output X, that occur when more than one model parameter is simultaneously perturbed [7]. Such effects can however be identified using a global sensitivity analysis technique, such as Latin Hypercube Analysis, as described in Section 5.

Latin Hypercube Sampling and Analysis
Latin Hypercube Analysis answers the question: how robust are model responses to global parameter perturbations? Latin Hypercube Analysis is a type of global sensitivity analysis that investigates the relationship between input parameters and output responses when all input parameters are simultaneously perturbed. The parameters that we want to perturb are (as in Section 4) denoted p i , where i = 1, 2, .., q. Thus the parameters p 1 , p 2 , ..., p q together span a parameter space of dimension q. It is impossible to test every possible combination of input parameter values if they are picked from continuous ranges. In fact, even if we select a finite number of parameter values r(p i ) to test for each parameter p i , or if we pick discrete parameter values, comparing every possible combination of parameter values may require us to produce an impractically large number of simulation runs. Thus performing in silico simulations for all possible combinations of input parameters will in many cases be at worst impossible, and at best impractical. In order to circumvent this issue, Latin Hypercube Sampling can be used [2]. It is a sampling technique that ensures comprehensive testing coverage over the parameter space whilst keeping the number of tested parameter combinations low enough to be applicable in practice [15,16]. After Latin Hypercube Sampling (Section 5.1), Latin Hypercube Analysis (Section 5.2) is used in order to assess global sensitivity.

Latin Hypercube Sampling
In the two-dimensional case, a Latin Square is an × square grid containing (traditionally Latin, hence the name) different symbols such that each symbol occurs exactly once in every row and exactly once in every column [17], as illustrated in Figure 3. Analogously, in the Latin Hypercube Sampling framework, consider two parameters p 1 and p 2 , spanning a parameter space of dimension q = 2, where both p 1 and p 2 are sectioned into intervals. We then pick combinations of input parameter values (or sampling-points) (p 1 j , p 2 j ), where j = 1, 2, ..., , such that every p 1 -interval is sampled from exactly once and every p 2 -interval is sampled from exactly once. Within the parameter range an interval, the sampled parameter value p i j is randomly selected (unless of course the interval contains only one possible value p i j ). Note that the j index denotes the coordinate combination that p i j belongs to, not the interval from which the parameter value p i j was taken. Thus there is no condition demanding that the values p i j are ordered in a way such that p i 1 < p i 2 < ... < p i . The analogy between a Latin Square and Latin Hypercube Sampling from a two-dimensional parameter space is illustrated in Figure 3. The Latin Square can be extended to higher dimensions to form a Latin Cube (dimension = 3) or a Latin Hypercube (dimension > 3) and, analogously, the two-dimensional sampling space illustrated in Figure 3 can be extended to q dimensions, spanned by the input parameters p 1 , p 2 , .., p q [17]. Figure 3: Left: An × Latin Square in which each Latin symbol occurs times, exactly once in each row and exactly once in each column. Right (analogously): A two-dimensional parameter space spanned by the input parameters p 1 and p 2 that are both sectioned into intervals. Using Latin Hypercube sampling, parameter combinations (p 1 j , p 2 j ) are sampled where j = 1, 2, .., and each p 1 -interval is sampled from exactly once and each p 2 -interval is sampled from exactly once. .
For each parameter p i , the total investigated parameter range is [min(p i ), max(p i )], where min(p i ) and max(p i ) respectively denote the minimum and maximum values of p i to be investigated. Now each parameter range [min(p i ), max(p i )] is sectioned into N intervals, and we denote these intervals by u 1 p i , u 2 p i , ...., u N p i . Note that all input parameters p i must be sectioned into the same number of intervals. If the intervals are uniformly distributed, then the size of an interval, w(p i ), is and the r-th interval u r p i has a parameter range such that where r = 1, 2, ..., N .
Note that there are more than one way to populate Latin symbols in a Latin Square, this can be realised by regarding Figure 3 and noticing that the A-symbols and the B-symbols cover the Latin Square in different ways. Analogously, and by extension, there are multiples ways to populate sampling coordinates in a Latin Hypercube Sampling framework. Some of these ways provide better coverage of the parameter space than do others [17], but details regarding such sampling-optimisation are outside the scope of this study. In this study, we use the built-in MAT-LAB function lhsdesign [18] to select which parameter combinations to use according to a Latin Hypercube Sampling approach, details about the implementation are available in the Appendix. Note that, in our case, all N intervals u 1 p i , u 2 p i , ...., u N p i for a parameter p i are uniformly spaced, but the choice of spacing can be adjusted to the specific application at hand [18]. Now let us address the choice of intervals N , as this is not straightforward. Using the Latin Hypercube Sampling framework, every parameter p i , where i = 1, 2, .., q, is partitioned into N intervals and, consequently, N combinations comprising q parameter values are sampled and tested. Compared to a small N -value, a large value of N will provide more data to use, and draw conclusions from, in the Latin Hypercube Analysis stage, however, it will also increase the computational cost in the Latin Hypercube Sampling stage. There is no strict rule for how to choose N , but suggested values for N in the literature are N = 2q for large values of q (i.e. high-dimensional parameter spaces) or N = 4q/3 which has been described to be 'usually satisfactory' [19,20]. Authors of the Spartan package use a lot larger numbers in their provided examples [2]. In this example study, we decide to use N = 100 uniform intervals. At the end of the day, the choice of N is up to the modeller, who must outweigh the (computational) cost of producing a large number of data samples, with the advantage of having a vast amount of data, and thus plentiful information, in the analysis stage. Details regarding quantitative choices of N are outside the scope of this review.

Latin Hypercube Analysis
During the Latin Hypercube Sampling process, N different points in the q-dimensional parameter space spanned by the input parameters p 1 , p 2 , ..., p q are selected as sampling-points, as shown in Step 1 in Section 5.3. One such sampling-point, C j , can be described by its coordinates in the parameter space so that C j = (p 1 j , p 2 j , ..., p q j ). Each sampling-point C j is used to generate n * output values X(C j ), where n * is determined using Consistency Analysis. Subsequently, the median output value, here denotedX(C j ), is computed for every C j . Now, our overall aim is to investigate the relationship between an input parameter p i and an output response X. We investigate this input-output relationship in two steps, one of which is qualitative and one of which is quantitative. In the first, and qualitative, step we produce two-dimensional scatterplots in which median output data, X(C 1 ),X(C 2 ), ...,X(C N ) =X(p 1 1 , p 2 1 , ..., p q 1 ),X(p 1 2 , p 2 2 , ..., p q 2 ), ...,X(p 1 N , p 2 N , ..., p q N ), are plotted over parameter values p i 1 , p i 2 , ..., p i N , for one of the input parameters p i . We do this for every input parameter i = 1, 2, ..., q and thus q scatterplots are created. By simply visually analysing the data in the scatterplots, we are able to make qualitative observations regarding the relationship between the input and the parameter. Examples of such observations are provided in Step 2 in Section 5.3 As a second step, we use a quantitative measure, such as the Pearson Product Moment Correlation Coefficient (or the correlation coefficient for short), to quantitatively describe the correlation between input parameters and output responses, as done in Step 3 in Section 5.3. The correlation coefficient is denoted r, where r ∈ [−1, +1], describes the linear association between the input parameter and the output response in terms of both magnitude and direction. A positive (linear) correlation between p i andX(C j ) means that if either the input value or the output value increases, so does the other one, and thus r is positive. Conversely, a negative correlation means that if either p i orX(C j ) increases, the other one decreases, and thus r is negative. The magnitude of r describes the strength of the correlation, where a magnitude of 1 corresponds to a strong linear association, and a small magnitude corresponds to a weak correlation. An r-value of approximately zero indicates that there is no linear correlation between the two investigated variables. Note that the Pearson Product Moment Correlation Coefficient picks up linear associations only, thus there may exist other, non-linear correlations that are not captured by the correlation coefficient r.Therefore it is important to, not only quantitatively compute input-output correlations, but to also qualitatively assess the relationships between inputs and outputs, via data visualisation in scatterplots as previously described.
The correlation coefficient, r i , describing the correlation between an input parameter p i , and an output response X (in median form) is given by [21], When it comes to interpreting quantitative input-output relationships based on the correlation coefficient r, there are no all-encompassing threshold values to use for descriptors such as 'weak', 'moderate', 'strong' [21,22,23]. Relationships quantified by correlation coefficient values close to the extrema 0 or 1 may be easy to describe as 'negligible' or 'strong', respectively but correlation coefficient values in the middle of the [0,1] range are more difficult to label. Various 'rule of thumbs' have been suggested in the literature but, at the end of the day, it is up to the researcher to appropriately judge what constitutes a 'weak', 'moderate' or 'strong' input-output relationship in the specific (modelling) application at hand, taking into account the research area, the number of data samples, and the range of investigated input values [22]. However, even without rigid descriptor threshold values, we can compare the correlation coefficient values for all input-output pairs and see which input values are the most influential within the ranges of regarded input values. As a guide, suggested correlation coefficient descriptor threshold values presented in the literature are listed in Table 5

Conclusion
This review is intended as a gentle, introductory review to three sensitivity analyses methods, namely, consistency analysis, robustness analysis and latin hypercube analysis. Information on how to implement these methods in MATLAB are available in the Appendix. Alternatively, all methods discussed in this review can be implemented using the R-based software package Spartan, developed by Alden et al. [2]. In fact, many of the proceedings and conventions used in this review follow those suggested by Alden et al. in order to allow the reader to, as easily as possible, use Spartan [2], if desired. Scrutinising mathematical models using uncertainty and sensitivity analyses methods is an important part in model development and, in many applications, knowledge about a model's robustness is crucial [24]. In the context of quantitative pharmacology, for example, a mathematical model may be used to guide preclinical or, ultimately, clinical proceedings. In such cases, understanding how confident we can be with model results, and how sensitive a model is to parameter perturbations, is of the utmost importance. With this in mind, the methods discussed in this review are implemented and interpreted in one of our recent studies, in which an agent-based, cancer cell model is used to predict in vivo treatment responses (post in vitro calibration) to an anti-cancer drug that may inhibit DNA damage repair [9].

Creating boxplots
The MATLAB function boxplot can be used to create boxplots. The input data in one column is represented by one box in the boxplot. For details regarding labeling and style alternatives, please see the MATLAB documentation [18].

Choosing Latin Hypercube Sampling points
A Latin Hypercube Sampling matrix can be created using the MATLAB function lhsdesign, which returns a matrix of size n × q, where n denotes the number of samples to be tested, and q denotes the number of input parameters to investigate (and thus perturb).

LHC _ Matrix=lhsdesign(n,q)
Each row i, in the created matrix (here denoted LHC_Matrix), corresponds to the ith sampling point. Each element (i, j) corresponds to the parameter value of the jth input parameter in sampling point i, where each parameter ranges between 0 and 1. For different criteria on how to chose the specific parameter values within each sampled interval, please refer to the MATLAB documentation [18]. Sampling points can, for example, be chosen in a way that maximises the distance between sampling points in the q-dimensional sampling space.

Qualitative and Quantitative Latin Hypercube Sampling Analysis
In order to qualitatively asses the correlation between an input parameter p, and an output response X, one can use the MATLAB function scatter. In the below listings, p and X are two data vectors.

scatter(p,X)
Further, to quantify the linear correlation between p and X, the MATLAB function corrcoef can be used to compute correlation coefficients.