Bifurcation and sensitivity analysis reveal key drivers of multistability in a model of macrophage polarization

In this paper, we present and analyze a mathematical model for polarization of a single macrophage which, despite its simplicity, exhibits complex dynamics in terms of multistability. In particular, we demonstrate that an asymmetry in the regulatory mechanisms and parameter values is important for observing multiple phenotypes. Bifurcation and sensitivity analyses show that external signaling cues are necessary for macrophage commitment and emergence to a phenotype, but that the intrinsic macrophage metabolism is equally important. Based on our numerical results, we formulate hypotheses that could be further investigated by in vitro experiments to deepen our understanding of macrophage polarization. 2010 MSC 93A30, 34C23, 49Q12, 92C42, 37N25


Introduction
Macrophages are highly versatile immune cells which, among other roles, eliminate pathogens and damaged cells through phagocytosis. They play a critical role in innate immunity and help to initiate the adaptive immune response through antigen presentation and cytokine signaling. Due to their diverse func- 5 tions and plasticity, macrophages are able to exhibit markedly different phenotypes, depending on the external signals they receive, e.g., microbial products, damaged cells, or cytokines. The continuum of macrophage activation and the diverse spectrum of pro-and anti-inflammatory phenotypes result in nuanced immune regulations [31].
A conceptual framework has been developed for the description of macrophage activation with two polar extremes being the most widely studied and best understood. On one end of the phenotype spectrum, M1-like macrophages are classically activated by the cytokine interferon γ (IFNγ) or by an endotoxin directly [30]. Once activated, M1-like macrophages release cytokines that inhibit the proliferation of nearby cells (including cancer cells) and initiate inflammation and an immune response.
Mixed phenotypes also exist, which share some (but not all) significant features with the M1-or M2-like phenotypes [1]. The existence of mixed pheno- 25 types has been particularly demonstrated in the tumor microenvironment [43].
Macrophage polarization is mediated in part, through the canonical Janusor TYK2-kinases (JAK)-Signaling signal transducers and activators of transcription (STAT) signaling pathway. Activation of STATs is primarily driven by ligand-stimulated cytokine receptors whereby STATs become phosphorylated 30 at a critical tyrosine residue leading to their release from the receptor complex where they then cross the nuclear membrane and reach chromatin. There they bind specific cognate DNA elements and participate in complex gene regulation processes.
Therefore, the phenotype expressed by a macrophage is identified through 35 the specific STAT activation. M1 polarization is associated with STAT1 activity, whereas M2 polarization is associated with STAT6 activity [29]. The M1 and M2 polarization process is dynamic and can be reversed under certain conditions. Individual macrophages can change their phenotype in response to local signaling cues [46,22,51]. This can be especially pronounced in 40 the tumor microenvironment and manifests in tumor associated macrophages, which can demonstrate both pro-tumoral and anti-tumoral activities [36].
Therefore, a better understanding of the polarization process of macrophages has the potential to guide the development of targeted cancer therapy to redirect the polarization towards a tumor suppressing microenvironment [47,51,7]. 45 Mathematical modeling is a useful tool to better understand macrophage polarization by validating or testing hypothesis, and making predictions about possible dynamics. To our knowledge, three previous studies based on ordinary differential equations (ODEs) have modeled macrophage polarization and plasticity [38,32]. While the authors in [32] showed bistable dynamics 50 of macrophage phenotypes when exposed to external signaling cues, the authors in [38] could show that after initial differentiation into M1 and M2, the M2 phenotype was ultimately dominating. Finally, the authors in [50] used a systems-level approach to present the complexity of signaling pathways and intracellular regulation which describe macrophage differentiation under IFN-55 γ, IL-4 signaling, and cell stress (hypoxia). With their model, the authors in [50] could replicate experimental results on macrophage phenotype markers and transcription factor regulations upon external perturbations, also for the tumor microenvironment.
All three models are built using generic formulations of self-stimulation and 60 mutual inhibition, which are also common building blocks in immune cell differentiation models [5,49]. Similar modelling approaches as for T-cell differentiation have been used for macrophages in e.g., [32,38], as T-helper cells differentiate in a similar manner [27,29]. Our goal is to use mathematical modeling to shed light on the polarization and regulatory signaling dynamics related to activation of macrophage phenotypes. We aim to build a simple model, which includes less parameters than the previous models by [32,38,50], but which shows similar complex dynamics. We conduct bifurcation and stability analyses to study its dynamical diversity, and relate these dynamics to biological observations. In addition, Global 70 Sensitivity Analysis (GSA) is employed to i) guide the model reduction and ii) to identify the most sensitive drivers of the system dynamics. Finally, sensitive parameters and input signals values are altered to study their effect on the dynamics.
For the rest of this paper, Section 2 describes our mathematical model in

Mathematical Model
Our mathematical model is based on the interactions specific to the macrophage lineage commitment signaling network. For this purpose, we simplify the network of macrophage functions in the liver from [37], and consider only IFNγ 85 (input signal S1) and IL-4 (input signal S2) as relevant cytokine signals. The levels of activated STAT1 (variable x 1 ) and STAT6 (variable x 2 ) are used in our model as proxies for the two macrophage activation states.
A schematic diagram of our model is given in Fig. 1. We model the dynamics of activated STATs with a pair of coupled nonlinear differential equations, described in equations in (1)-(2).
where ≡ d dt . All parameters are assumed to be constant, positive and real numbers, except n 1,2 , which are integers.

90
The description of all model parameters is provided in Table 1. Parameter Description a 1,2 Strength of self-stimulation b 1,2, Basal activation rates n 1,2 Exponents in the Hill functions for self-stimulation k 1,2 Thresholds in the Hill functions for self-stimulation l 1,2 Exponents in the Hill functions for mutual inhibition p 1,2 Thresholds in the Hill function for mutual inhibition q 1,2 Deactivation rates S 1,2 Input signal strength

Model formulation
The equations in (1)-(2) were adapted from the T-cell model in [49]. The equation for x 2 is based on the assumption that both type I and type II interferons inhibit IL-4-induced STAT6 activation in human monocytes in a SOCS-1-95 dependent manner [10], and therefore differs from the model formulation in [49]. This change results in an asymmetry in our equations in that STAT6 inhibits both the input signal and self-stimulation, but STAT1 affects only the input signal [44]. Furthermore, we reduced model complexity by fixing the Hill coefficient in equation (4) to 1. Also, the signal input function in [49] was simplified 100 to a single parameter (S 1 , S 2 , respectively) for each phenotype in our model.
In our model equations, the parameters α i represent the maximal activation rate of STAT due to self-stimulation. STATs are however also activated at low background levels (b i ) in the absence of cytokine stimulation [9]. STATs are also inactivated by dephosphorylation, and we assume this rate is linear (terms The fact that STAT1 and STAT6 are autocrine [48,15], is captured by the stimulating Hill functions in the model equations (1)- (2). Finally, we assume respective activation of STAT1 and STAT6 via IFNγ (S 1 ) and IL-4 (S 2 ) [33].
We use stimulating (equation (3)) and inhibiting (equation (4)) Hill func-110 tions to describe STAT self-stimulation and mutual inhibition [42], respectively. The rationale behind the choice of these generic functions is that self-stimulation and inhibition are complex, non-linear processes, which consist of several individual steps. For example, in the process of self-stimulation, cytokines from the macrophage are secreted to stimulate helper T-cell differentiation [23]. Dif-115 ferentiated helper T-cells then secrete cytokines which in-turn stimulate the macrophage differentiation. However, detailed knowledge about these individual steps is unknown, which makes it difficult to derive mathematical equations for each step. In addition, we assume that the response in self-stimulation is sigmoidal, depending on the "dose" of input signals. Therefore, the Hill function 120 is used and replaces the need to model the steps individually [42]. A similar argument was used for the inhibitory Hill function.
In the Hill function of equation (3), k i represents the signaling level at which STAT stimulation is half-maximal and the Hill coefficient n i governs the steepness of the Hill function in that as this value grows, the function becomes more 125 switch-like. For the inhibitory Hill function, the parameters play a similar role.

Numerical Methods
In this section we provide the detailed description of numerical methods we employed. 130 We explore parameter variations and analyze how the different parameter sets affect variability in the system states by using three parameter sets: the initial set Θ 0 , and two variation sets, Θ 1 and Θ 2 . The parameters in the initial set Θ 0 are adapted from [49], while the variation sets Θ 1 and Θ 2 are derived using nullclines. All three parameter cases are presented in Table 2.

Bifurcation and stability analysis
We expect our model, for all three case scenarios, Θ 0,1,2 , to exhibit at least bistable dynamics, similar to the original model. Thus, we first conduct bifurcation analysis to further investigate the impact of different parameter sets on model dynamics.

140
Bifurcation analysis aims to detect critical points of the bifurcation parameters, where the system dynamics change qualitatively in the long-term [17]. Given the biological importance of external signaling cues (INF-γ and IL-4) in the macrophage polarization process [46], we are primarily interested in determining how the system dynamics change based on varying input signals (i.e., S 1 145 and S 2 ). We therefore consider S 1 and S 2 as main bifurcation parameters, with the other parameters set to their values in Table 2. The bifurcation diagrams from equations (1)-(2) were obtained using the software package XPPAUT [12]. Details on numerical settings to draw bifurcation diagrams can be found in Appendix A.1.

150
We define states of STAT activation based on model-specific thresholds. An activation level is defined as low, if S 1,2 ≤ 1.0, and as high, if S 1,2 > 1.0. It is then the ratio of STAT1 to STAT6 activation, that characterizes a macrophage phenotype. The threshold levels are chosen to allow a consistent classification of phenotype cases in our model, although they only represent relative levels.

155
Stability analysis was performed by numerical simulations in Matlab.

Sensitivity analysis
We perform sensitivity analysis to identify parameter sets that have the greatest influence on the model outputs (e.g., STAT1 and STAT6 activation), and act as key drivers of macrophage polarization. Local sensitivity analysis 160 quantifies changes in the model with respect to perturbation of a single parameter at-a-time, but is not recommended for non-linear models [52]. In contrast to local sensitivity, GSA methods explore the effects of changes in parameter values on model outcome by varying all parameters simultaneously. Although several global methods exist (e.g., extended Fourier Amplitude Sensitivity Test, 165 Latin Hypercube Sampling and Partial Rank Correlation Coefficient [45,28]), we employ Sobol's method [39], a variance decomposition technique. We chose this method because it makes no assumptions about the relationship between model inputs and outputs in contrast to, for example, the Partial Rank Correlation Coefficient method, which requires monotonicity. Additionally, Sobol's 170 method considers interactions between parameters. A detailed description of Sobol's method can be found in Appendix Appendix A.2.
We implemented Sobol's sensitivity analysis using the SALib package [18]. We varied parameters 15% in each direction from their baseline values (i.e., parameter sets Θ 0,1,2 in Table 2. We consider these scenarios separately. In all cases, we generated 300, 000 parameter set samples. The selected outcome of interest for the analysis is the ratio of STAT1 to STAT6 activation, which is responsible for macrophage polarization to specific phenotypes.

Perturbation in sensitive parameters
Based on results of the GSA, we explore the effect of perturbations in sensi-180 tive parameters on macrophage polarization dynamics. For illustrative purposes, we will only consider perturbations in the most sensitive parameter (q 2 ) on case Θ 0 . Understanding the effect of dephosphorylation on system dynamics is especially important as deactivation rates change often in biological settings [19]. By changing q 2 and keeping all other parameters fixed, we hope to better un-185 derstand its individual effect on the relation between external input signals and activation of transcription factors.

Model Results
The results of the numerical simulations are presented in this section. 190 We observe bistability, tristability, and quadstability for different combinations of the S 1 and S 2 based on the three parameter cases Θ 0,1,2 , respectively.

Bistable case
With the initial parameter set Θ 0 we observe two stable fixed points, exhibiting bistable behavior. These steady states represent state variable ratios 195 (x 1 /x 2 ) with i) high/low and ii) low/low levels.
Though the detailed bifurcation diagrams are omitted here, we validate this bistable behavior by numerically solving equations (1)-(2) with the parameter set Θ 0 . The most interesting behavior observed is that x 1 and x 2 go through a switch before converging to their respective stable fixed points, as shown in 200 Fig. 2(a). The solution trajectory of this switch behavior (in solid black) in the phase plane is provided in Fig. 2(b). Note that only two fixed points are present even though there seems to be another fixed point on the upper left part in the phase plane because of the proximity of the x 1 -and x 2 -nullclines. The bistable behavior is further confirmed by the basin of attraction shown in Fig. 2(c).

Tristable case
With parameter set Θ 1 , three stable steady states of (x 1 /x 2 ) are observed with i) low/high, ii) high/low, iii) low/low levels. The third steady state represents a situation where both STAT1 and STAT6 are present at similar levels.
Numerical solutions that converge to different stable fixed points are shown 210 in Figs. 3(a)-(c). The respective solution trajectories in the phase plane are shown in Fig. 3(d).
Because of the increased values S 1 = S 2 = 4 for this case, there are two additional intersections between the x 1 and x 2 -nullclines compared to the bistable case, as can be seen in the phase plane of Fig. 3(d). This results in the addition 215 of two fixed points, one of which is stable and the other is unstable. Thus, if we start with the same initial condition used in Fig. 2(a), the trajectory converges to the new stable fixed point with high x 2 /low x 1 , which was not observed in the bistable case. As further confirmed by the basin of attraction of Fig. 3(e), the other two stable fixed points remain as before.

220
It is the ratio of STAT1 (x 1 ) to STAT6 (x 2 ) activation levels that defines the polarization of a macrophage into the M1 or M2 phenotype [46,32]. In our results, a high level of activated STAT1 in presence of low activated STAT6 levels defines the M1 phenotype [13], while low levels of activated STAT1 and high levels of activated STAT6 define the M2 phenotype. Low STAT1 and 225 STAT6 activation levels represent a "hyporesponsive" phenotype that has not been described in the current literature. This phenotype might however have biological relevance (e.g., for cancer therapy), as an intermittent phenotype between M1 and M2. For example, recent studies by [3,6,25] have shown that tumors are initially characterized by M1 or an intermittent phenotype state, 230 while advanced cancer is defined by M2 phenotype. It is therefore possible that this "hyporesponsive" phenotype describes another intermittent phenotype that appears during this transition.

Quadstable case
Using the last parameter set Θ 2 , our model demonstrates quadstable behav- The situation where both STAT1 and STAT6 have high activation status is, however, unique to the quadstable case. High activation of both STAT1 and STAT6 shows the existence of an intermittent phenotype [1], which bears characteristics of both the M1 and M2 types. Several of such intermittent states have been identified, for example, M2a, M2b, M2c and M2d [34]. The 245 intermittent phenotype can also represent a transformation state, in which M1 branches to M2, and vice versa [8].
To understand how a varying input signal changes the activation of STAT1 and STAT6, we illustrate, based on Figs. 4(b)-(c), how one should read the bifurcation diagram: Figs. 4(b)-(c) have to be read simultaneously, starting from 250 S 1 = 0 and then increasing the S 1 value while following the bifurcation trend. Note that while S 1 is varied, all other parameters values are kept unchanged. By varying S 1 from 0 to around 12, x 1 is on the lowest stable branch while x 2 is on the highest stable branch. Increasing S 1 input signal beyond 12, x 1 and x 2 will follow the bifurcation trend up and down, respectively, to the next stable branch with x 1 activation level between 1 and 2.2, and x 2 activation level between 1.8 and 1.3. To reach the third stable branch, input signal S 1 is decreased (to follow the bifurcation trend) until x 1 and x 2 jump from the second red branch to the third branch. The third branches spans values between 0.3 and 1 for x 1 , and values between 0.3 and 0.7 for x 2 . When on the third branch, 260 S 1 input signal will be increased again, at an input signal of around 7, both x 1 and x 2 will jump onto the respectively highest and lowest branch. Figs. 4(a)-(d) can be read similarly.
In Figs. 4(b)-(c), we observe furthermore that for high, S 1 > 18 levels, the state variables x 1 and x 2 are committed to highest and lowest activation levels, respectively.
It is interesting that in the case of quadstability, the system is committed to the high/low state (see Figs. 4(b)-(c)) for high S1 values, while this could not be observed for bistable or tristable situations. Biologically, an irreversible switch into the M1 phenotype means that the macrophage will no longer be 270 able to change its phenotype when exposed to changing input signals. This suggests that for high self-stimulation in the presence of high INFγ and low IL-4 signals, the system can commit to M1 phenotype and stay reversible for the M2 phenotype. In parameter set Θ 2 , STAT1 has higher self-stimulation than STAT6, i.e., a 1 > a 2 and n 1 > n 2 . This might be a crucial driver for

Identification of key drivers of macrophage dynamics through global sensitivity analysis
The most sensitive parameters for the bistable case using total sensitivity as a metric are, in descending order, q 2 , q 1 , k 2 , S 1 , k 1 , a 2 , S 2 (see Fig. 6(a)).

285
The four most sensitive parameters for bistable and tristable cases, shown in Figs. 6(a)-(b), respectively, agree and the next three most sensitive for each case are common (k 1 , a 2 , S 2 ) but reordered. Fig. 6(c) shows that the most sensitive parameters in the quadstable case are consistent with results from the previous two cases. 290 In terms of the pathways, this indicates that deactivation rates of both STAT1 and STAT6 (q 1 and q 2 , respectively) are highly sensitive, as well as the input signal for M1 polarization, INFγ (S 1 ). Parameters k 2 and k 1 are also sensitive, and both relate to the response of the Hill functions for selfstimulation. These parameters govern the concentration at which the switch takes place. In all cases, k 2 is more sensitive than k 1 . Parameters S 2 and a 2 are the signaling input for M2 polarization (IL-4) and the maximum rate at which STAT6 stimulates its own activation via a regulative feedback mechanism.    In all instances, the parameter q 2 has the highest total sensitivity index. The cases of bistability and tristability have the same most sensitive seven parameters q 2 , q 1 , k 2 , S 1 , k 1 a 2 , S 2 with only the ordering of the last three altered. For the quadstable case, q 2 is also the most sensitive, with k 2 and a 2 moving up in the ordering compared to the previous two cases.

Effect of perturbation in sensitive parameters
states, and to increase the external stimuli needed to evoke a fate change. This example indicates that deactivation rates can contribute to the robustness of the 305 dynamical system to variations in external stimuli. In particular, it illustrates that deactivation of STAT1 and STAT6 plays an essential role in macrophage polarization, as deactivation rates indirectly affect inhibition of external input signals on the opposite state variable, while self-stimulation affects its own state variable.

Discussion
In this work, we develop and explore a novel mathematical model for the dynamics of macrophage polarization and identify key parameters of the multistable dynamics. We validate that macrophage polarization is not strictly bipo- lar, but can consist of multiple phenotypes. Our macrophage polarization model 315 is the first to show that a two-dimensional and simple model can predict bistable, tristable and quadstable phenotypes. We could validate known phenotypes (i.e., M1 and M2) and have uncovered unknown, intermittent ones (i.e., low/low and high/high) with a mixed phenotype expression. From a biological perspective, the intermittent pheno-320 types might be more likely in in vivo settings than the extreme M1 and M2 cases, which are studied in cell cultures. Although, we cannot rule out that there exist more than four different phenotypes for our system, our findings are supported by those in [26], where the authors identified maximal four stable states given a similar model formulation. To our knowledge, only one previ-325 ous study by [32], which studied a more complex model and applied also twoand three-dimensional bifurcation analyses, could identify a broader spectrum of known (e.g., M2a and M2b) and unknown macrophage phenotypes. Our identified unknown phenotypes can however not be compared directly to those in [32], because the authors classified STAT activation into high, medium and low levels, while we only made a distinction between high and low. In addition, such classification states are model dependent.
Sensitivity analysis of our model revealed the high impact of the deactivation rates, q 2 and q 1 , on the ratio of STAT1 to STAT6 activation at steady state, used as a proxy for M1 and M2 phenotype, respectively. Parameters k 2 and α 2 335 were also identified as sensitive because both of these parameters are related to the self-stimulation of STAT6 activation. Our most sensitive model parameters are similar to those identified in [41,50]. These sensitive parameters agree with results of our bifurcation analysis, where parameters of self-stimulation and deactivation seemed to have a profound impact on the dynamics. For example, in the quadstable case, parameters of self-stimulation might explain the observed system commitment and the emergence of an additional phenotype, while results of varying deactivation rates changed the response to external signaling cues, as can be seen in Fig. 7. The consistency in identifying sensitive parameters from bifurcation and sensitivity analyses is however expected, because a properly 345 designed analysis should reveal bifurcation parameters to be sensitive [28].
In summary, bifurcation and sensitivity analyses showed that external signaling cues are necessary for macrophage commitment and emergence to a phenotype, but that the intrinsic macrophage metabolism (represented by selfstimulative factors and deactivation) is equally important [14,2]. It should be 350 noted that the intrinsic metabolic pathways, which enabled fate commitment in the quadstable situation, are masked by the generic nature (i.e., Hill function) of our model. Intrinsic metabolic pathways in macrophages are in general variable [14], and one can distinguish, for example, between glucose, lipid, iron or amino acid metabolic pathways [2].

355
Our results support the expectation from the model diagram ( Fig. 1 that the system's outcome also depends crucially on the self-stimulation of x 2 . Because the equations are not symmetric (i.e., in the second equation the stimulatory and inhibitory Hill functions are additive, not multiplicative as in the first equation), the parameters associated with STAT6 have a stronger impact on the 360 model outcome. This observation is also reflected in the asymmetric values of a 1 , n 1 and a 2 , n 2 in Θ 2 . The asymmetry illustrates that lower values of a 2 , n 2 have the same effect on systems dynamics as higher values of a 1 , n 1 . The parameters in Θ 0,1 however are symmetric, because they were adapted from the mathematical model in [49], which has a symmetric model structure. The need 365 for an asymmetry in self-stimulation dynamics of STAT1 and STAT6 might be explained by the experimental finding that the signaling pathway induced by IFN dominates over the signaling pathway induced by IL-4, according to the authors in [35]. This explanation is furthermore in accordance with our finding of an irreversible switch to the M1 phenotype for high concentrations of INF-γ.
Furthermore, we illustrated how STAT deactivation impacts macrophage polarization by influencing the robustness to external stimuli. The authors in [40] pointed out that the effects of deactivation are, however, not well understood for macrophages. Finally the knowledge of sensitive parameters for macrophage polarization might guide the conduction of future in vitro experiments and thus deepen our understanding of macrophage polarization. We therefore suggest that the following hypotheses, which resulted from our analyses, to be tested experimentally: H1: The response-time and sensitivity of STATs to cytokine signaling levels can be altered by changing deactivation rates.

380
H2: Once macrophages are committed to a phenotype, further stimulation via cytokines leaves them unchanged. H3: Intrinsic metabolic characteristics, which correspond to aspects of selfstimulation and deactivation, determine the range and variability of observable macrophage phenotypes.

385
H4: There exist intermittent phenotypes with equal STAT activation levels (i.e., defined by STAT phosphorylation) in in vitro settings.

Model limitations and future work
A clear advantage of our model is its simplicity and its ability to exhibits complex dynamics in terms of multistability. The price of simplicity is, however, 390 the elimination of several key biological signaling components. One example is NF-κB, a protein complex which interacts with type 1 interferons, among other signals [11]. Future work could inspect a more refined signaling network, based on our model formulation.
Another limitation of this work is that our model considers a single macrophage 395 whereas in reality there are entire populations of macrophages which influence each other. However, understanding how a single macrophage reacts to its microenvironment is a first step to understanding population level behavior. Mathematical models are needed to address macrophage polarization on population level and to consider input signals beyond IFN-γ and IL-4, while 400 incorporating knowledge of dynamics of a single macrophage.
Finally, the absence of empirical data in this manuscript could be considered a limitation, especially as there exist data on the activation status of STATs. The primary focus of this manuscript has been to understand the qualitative characteristics of the proposed model. Hence extending the analysis to include 405 empirical data is beyond this scope.
Our model represents also a solid first step towards analyzing stochastic gene expression in macrophages. In future work, we will make use of the chemical master equation and analyze how switching probabilities between different phenotypes change with variations in extrinsic and intrinsic noise levels.

Declaration of interest
The authors declare that they have no known competing financial interests or 415 personal relationships that could have appeared to influence the work reported in this paper.

Funding Sources
The research of SR was funded by Trond Mohn Foundation, Grant No. BFS2017TMT01.

Data statement
This manuscript uses no biological or experimental data. All Figs. can be reproduced using the mathematical model, the parameters and numerical specifications presented in this manuscript.   1 Details can also be found at http://www.math.pitt.edu/~bard/bardware/tut/xppauto.html Abbreviations: N tst, number of mesh intervals for discretization of periodic orbits, N max, maximum number of steps taken along any branch, N pr, give complete info every Npr steps, Ds, initial step size for bifurcation calculation, Ds min , minimum step size, Ds max , maximum step size, P ar min , left-hand limit of the diagram for principal parameter, P ar max , right-hand limit of the diagram for the principal parameter the effect of varying x i and x j simultaneously, additional to the effect of their individual variations, termed a second-order sensitivity. Higher order terms have analogous interpretations.
Assuming that f (x) is square integrable, the functional decomposition may be squared and integrated and the total variance D can be defined as The partial variances from squaring and integrating the right hand side of Appendix A.1 are of the form D i1,i2,...i k = · · · f 2 (x i1 , x 12 , . . . , x is dx i1 dx i2 . . . dx i,s (Appendix A.3) These integrals can then be approximated with Monte Carlo integration, and the Sobol sensitivity indices are calculated by the ratio of partial to total variance, representing the fraction of total variance which is attributed to individual model parameters or to combinations of parameters.
Furthermore, the total effect sensitivity index was proffered as an extension of the Sobol sensitivity index to quantify the overall effect of a parameter alone and in combination with any other parameters on model output [20]. This is defined to be S Ti = S i + S ci (Appendix A.5) where S ci is the set of sensitivity indices in which parameter x i appears.