Variance decomposition in stochastic simulators

It has long been recognized that stochastic simulators^1–7 offer an attractive computational means for simulating the chemical master equation (CME). These algorithms are based on describing the evolution of the number density of individual chemical species on the basis of propensity functions, which provide a probabilistic characterization of the rate of progress of elementary reactions. In its elementary formulation, the reaction is advanced by sampling the time for next reaction to occur, whereas the system remains frozen in its current state in the mean time. After a reaction has occurred, new samples are drawn, and the procedure is iterated in order to determine successive elements of a Markov process associated with the evolution of the state vector of species number densities. By computing multiple realizations of such chains, the algorithm affords in a straightforward fashion the possibility to determine various statistical moments of the solution.

One of the difficulties associated with stochastic simulators is that it is not generally straightforward to assess the impact of individual reaction channels on the evolution of the solution, nor to quantify the effect of uncertainties on the variability of the solution. This is especially the case for stochastic systems that do not involve parametric uncertainty, which are the focus of the present study. The central question we aim to address is whether it is possible, without perturbing the reaction rates (propensity functions), to exploit the stochasticity inherent in stochastic simulations for the purpose of quantifying the impact of individual channels on the evolution of particular species of interest, and also potentially characterizing the role of channel interactions. As outlined in Sec. II, our approach to this question is based on reformulating classical Stochastic Simulator Algorithms (SSAs, see Refs. 1 and 2) using independent standardized Poisson processes similarly to the next reaction method proposed in Ref. 8. This reformulation affords a straightforward identification of independent random realizations for individual channel dynamics, and consequently enables us to exploit sampling strategies in order to quantify the impact of the stochasticity associated with elementary channels, or groups of channels. As discussed in Sec. III, the Sobol-Hoeffding (SH) decomposition is applied for this purpose. In particular, this enables us to adopt well-established Monte-Carlo sampling strategies to obtain estimates of the contributions of individual channels, and of channel interactions, to the variance of specific quantities of interest. Implementation of the resulting methodology is illustrated in Sec. V through applications to well-known systems, including the birth-death (BD), Schlögl, and Michaelis-Menten models. Major conclusions are summarized in Sec. VI.

The present work focuses on the analysis of the respective importance of inherent sources of variance in stochastic simulations of chemical reaction networks. Stochastic effects arise predominantly when the reacting system involves a small number of reactant molecules, due to discrete evolutions. Specifically, we consider a well-stirred chemical system involving M_s molecular species {S₁, …, S_Ms} and K_r reaction channels {R₁, …, R_Kr}. We denote $X ( t ) = ( X 1 ( t ) , … , X M s ( t ) ) T ∈ N 0 M s$ the state of the system at time t, where X_i(t) is the instantaneous number of molecules of species S_i. The dynamics of the system are prescribed by the definition of the propensity functions a_j and state-change vectors ν_j ∈ ℤ^M_s associated to the reaction channel R_j. The propensity functions,

are defined such that, for infinitesimal time dt, a_j(x) dt is the probability that one jth reaction occurs in the time interval [t, t + dt) given that X(t) = x, whereas ν_j = (ν_1,j, …, ν_Ms,j) describes the change in the molecular population due to a unique reaction R_j.

The system then obeys the CME,⁹ 

where P(x, t|x₀, t₀) is the probability that X(t) = x given that X(t₀) = x₀.

Solving (1) is impractical even for systems with small M_s and K_r. The Gillespie SSA exactly emulates the chemical master equation by generating samples of X(t) with probability law P(x, t|x₀, t₀). In Sec. II A we briefly recall the original SSA. In Sec. II B we introduce a slightly modified version of the SSA, which is suitable for the analysis of the variance.

Gillespie’s SSA^1,2 consists in the construction of a Markov-chain X(t) having for distribution P(x, t|x₀, t₀). The first idea of the SSA is that, given X(t) = x, the probability of the next reaction (irrespective of the channel) to occur in the infinitesimal time interval [t, t + dt) is

In other words, the time elapsed to the next reaction, τ, follows an exponential distribution with mean 1/a₀(x). Second, the relative probability of reaction channel R_j to fire first is p_j(x) = a_j(x)/a₀(x). Gillespie’s Algorithm I proceeds along these two stages to advance the state of the system: first, drawing at random a time increment to the next reaction event, and second, selecting at random the firing reaction channel to update the state. In more detail, given the current state X of the system, the SSA first evaluates the K_r propensity functions a_j and their sum a₀. Drawing r₁ uniformly in (0, 1) the algorithm sets τ = − log(r₁)/a₀, and selects the index of the firing channel using the largest l such that $∑ j < l a j < a 0 r 2$⁠, where r₂ is another random number, independent of r₁, drawn uniformly from (0, 1). Time and system state are then advanced, t←t + τ and X←X + ν_l, and the procedure is repeated until the final time is reached.

It is thus seen that the SSA advances in time using a consolidated next reaction time distribution and relative reactions probabilities, so the method inherently depends on the whole state of the system. In other words, the algorithm fully couples the reaction dynamics. It requires drawing 2 random numbers from the unit interval for every elementary event (occurrence of a reaction). In practice, a pseudo-random number generator RG is used to simulate the sequence of random numbers r₁ and r₂. Classically, the pseudo-random number generator RG provides the exact same sequence of pseudo-random couples (r₁, r₂) when properly seeded (initialized). As further discussed below, this property of pseudo-random generators is attractive for our purpose, because it enables us to generate the same sequence of pseudo-random events multiple times.

However, the SSA has several drawbacks that make it unsuitable for the separation of the variances caused by the different channels. The central drawback is that SSA combines all channels to determine the next reaction time, while the analysis of the variability inherent to specific channels, as proposed in Sec. III, requires us to distinguish and control the dynamics of individual channels. An alternative way to understand the issue is that for a fixed sequence of pseudo-random numbers, the resulting trajectory of X(t) is not stable to permutation of the reaction indexation. This drawback motivates the use of the Next Reaction Algorithm (NRA) outlined below.

Following the next reaction method proposed in Ref. 8, the dynamics of the stochastic state vector X(t) can be expressed through the following random time-change:¹⁰ 

where N_j(τ) are independent Poisson processes with unit rates, and the scaled time t_j are given by

Thus, the system state can be expressed as a linear combination of K_r standard Poisson processes, with individual time scaling given by the time-integrals of the corresponding propensity functions from t₀ to t. Expression (2), together with (3), results from the fact that Poisson processes have independent increments distributed according to

Consequently, the probability of R_j to fire once in the infinitesimal time interval [t, t + dt) is

which is precisely the definition of the propensity function a_j.

The representation of X(t) through (2) and (3) suggests an alternative version of the original SSA, hereafter referred as the NRA, which is a slightly modified version of the next reaction method proposed in Ref. 8. Algorithm II introduces two additional vectors, (τ₁, …, τ_Kr) and $( τ 1 + , … , τ K r + )$⁠, respectively, to store the unscaled time of the Poisson processes and their next reaction times. These two vectors are used to keep track of the trajectories of the K_r independent Poisson processes. Given the current state X of the system, normalized times (τ₁, …, τ_Kr) and next reaction times $( τ 1 + , … , τ K r + )$⁠, one can deduce which channel will fire first from scaled times advancement rules (3). Specifically, if no other channel fires before it does, the jth channel will react after a time lapse $d t j = ( τ j + − τ j ) / a j$⁠. Letting $l=arg min j d t j$ be the first channel to fire, the procedure selects at once both the firing channel and the time lapse (dt_l) to the next reaction, unlike the SSA. Having found l, time is advanced (t←t + dt_l), system state X is updated (X←X + ν_l), the unscaled times vector is updated (τ_j←τ_j + a_jdt_l), and the unscaled next reaction time $τ l +$ is drawn at random (⁠$τ l + ← τ l + −log ( r l )$⁠, where r_l is drawn at random from (0, 1)). The procedure is then repeated until t reaches the final time. The NRA (Algorithm II) involves drawing a unique random variable per time step (whereas SSA draws two) but requires the initialization, computation, and storage of the scaled-times and next jump-times associated to the K_r Poisson processes associated to the reaction channels. However, both algorithms result in states X that are equal in distribution, that is,

Our interest in considering the NRA lies in the possibility of maintaining constant realizations of selected standardized Poisson processes N_j while varying others. It amounts to keeping invariant the sequences ${ τ j + }$ of unscaled jump-times. This possibility is enabled by the introduction of the K_r individual pseudo-random number generators, one for each reaction channel, when in fact a unique random number generator would have been sufficient in NRA (Algorithm II). The use of K_r independent pseudo-random number generators distinguishes NRA (Algorithm II) from the next reaction method in 8. As for SSA, the storage of these K_r independent sequences is not required as they are generated on the fly using the pseudo-random number generators RG_j. Doing so, the realizations of the Poisson processes are fully determined by imposing the seeds of the K_r pseudo-random number generators. As a result, we now have a clear understanding of the inherent source of stochasticity which is kept invariant when using the same set of seeds, namely, the underlying Poisson processes, and a fine control over this source of stochasticity. This control over the inherent stochasticity enables a decomposition of the variance of the state vector, and the analysis of the variability induced by distinct channels, as discussed in Sec. IV.

In parametric sensitivity analysis, one is interested in char acterization of the variability of the deterministic output (or quantity of interest) of a model when some parameters are changed. Local sensitivity analyses characterize the dependence of the deterministic output about a particular value of the parameters, usually by computing the first derivatives of the output with respect to the parameters. Contrary to local methods, global sensitivity analyses aim to elucidate the global contribution of a particular input on the total output variability, assuming a certain distribution of the parameters. In the ANOVA (analysis of the variance) or Sobol’s variance decomposition, the sensitivity indices are associated to the input parameters via a decomposition of the output variance.^11,12 When the input parameters vary independently, Sobol’s de composition is orthogonal and the definition of the sensi tivity indices is immediate from the so-called partial variances (see below). The case of dependent inputs requires appropriate definition and interpretation of the sensitivity indices (see, for instance, Refs. 13 and 14). Sobol’s decomposition is also related to High Dimensional Model Representation (HDMR), which corresponds to a truncated SH decomposition; such decompositions have in particular been used for the construction of surrogate models (see, for instance, Refs. 13 and 15–17). In the context of reaction networks presenting an inherent stochastic dynamics, analyses have been restricted to averaged functionals of the stochastic model solution: the analyses characterize the sensitivity of such averages with respect to parameters, for instance coefficients in the propensity functions defining the stochastic network.^18–22 

Parametric sensitivity analyses have to be contrasted with the use of the Sobol decomposition we are proposing in the present work. Here, we assume no variability in any model parameter. Instead, we seek to quantify the respective contributions of different reaction channels to the variance of a given functional of the stochastic model solution.

In this section, we summarize relevant aspects of the Sobol-Hoeffding decomposition, provide a brief outline of how these concepts can be applied to define useful sensitivity indices, and detail a Monte-Carlo sampling procedure for their estimation.

Consider a vector N = (N₁, …, N_D) of D independent random quantities defined on an abstract probability space $P= ( Ω , Σ , μ )$⁠. Let F : N ↦ F(N) ∈ ℝ be a second-order random functional in N, that is,

Let $D$ be the power set of {1, …, D}. For $u∈D$⁠, denote |𝔲| = Card(𝔲) and $u ∼ =D∖u$⁠, such that $u∪ u ∼ =D$ and 𝔲∩𝔲_∼ = 0̸. Given $u∈D$ we denote N_𝔲 the sub-vector of N with components (N_𝔲1, …, N_𝔲|𝔲|), so N = (N_𝔲, N_𝔲∼).

For the assumptions considered, the function F(N) has a unique orthogonal decomposition of the form²³ 

Equation (6) is called the SH decomposition of F. For instance, in the case D = 3 the decomposition is expressed as

The orthogonality condition for the functions implies that ∀𝔲≠𝔰,

and the SH functions F_𝔲 are recursively defined according to¹¹ 

where $E F n u$ is the conditional expectation of F(N) given N_𝔲 = n_𝔲, namely,

For instance, in the previous example with D = 3, we have

Clearly, $F 0̸ =E F$⁠, and the orthogonality of the decomposition implies $E F u =0$ for 𝔲≠0̸. Decomposition (6) being orthogonal, the variance $V F$ is decomposed into the sum of partial variances corresponding to the variances of the SH functions,

Recalling that $E F u ≠ 0̸ =0$⁠, and using definition (7), the partial variances $V F u$ can be expressed as

where we have used $E E F N u =E F = F 0̸$ and the orthogonality of the SH functions.

The partial variance $V F u$ measures the contribution to $V F$ of the interactions between the variables N_𝔲. Since there are 2^D such partial variances, the sensitivity analysis is usually reduced to a simpler characterization, based on first and total-order sensitivity indices associated to individual variables N_i = N_{i} or group of variables N_𝔲.

The first-order sensitivity index S_𝔲, defined as¹² 

represents the fraction of the variance $V F$ due to the variable N_𝔲 only, i.e., without any interaction with variables in N_𝔲∼. For instance,

The total-order sensitivity index T_𝔲, defined as¹² 

accounts in addition for the interactions between variables N_𝔲 and N_𝔲∼. For our example with D = 3, we have

In fact, it is immediate to show from the definition of the total-order sensitivity index that S_𝔲 ≤ T_𝔲 ≤ 1. In addition, the following inequalities hold for the case of sensitivity indices associated to a single variable:

while

In the latter case, F is said to be additive since there is no interaction between the N_i’s.

From (10) and (11) it is seen that the computation of the first and total order sensitivity indices essentially amounts to the determination of the variance of conditional expectations. For this purpose, we shall rely on the Monte Carlo method proposed in Ref. 24.

Consider two independent random sample sets $N I$ and $N I I$ of M realizations of N. The conditional variance $V E F N u$ can be estimated using the following average:²⁴ 

Here, we have denoted N^I,(i) and N^II,(i) the ith elements of $N I$ and $N I I$⁠, respectively. Introducing the classical (⁠$N I$-)sample estimator for the mean,

and variance of F,

the Monte Carlo estimator of the first-order sensitivity index is

A similar expression can be derived for the estimation of the total-order sensitivity index T_𝔲 in (11),

It is seen that the MC evaluation of the sensitivity indices requires M evaluations of F (for the elements of $N I$⁠), and M new function evaluations for each S_𝔲 or T_𝔲. In particular, the MC estimation of the first and total-order sensitivity indices S_{i} and T_{i} for the D input variables requires a total of (2D + 1) M model evaluations.

Our objective is now to propose a decomposition of the variance for a functional g(X(t)), where X(t) is the stochastic output of a stochastic simulator. In all the examples of this section, the functional g is simply a component of X(t) at some time t = T, but the proposed methodology can be readily extended to more complex functionals including path integrals and exit times, provided that g(X(t)) is a second-order random variable.

To decompose $V g ( X ( t ) )$ into contributions arising from the different K_r reaction channels and their interaction, and accordingly define the first and total-order sensitivity indices, associated with the variability of an individual channel (or group of reaction channels), an appropriate definition of the conditional expectations is required. Specifically, we need to provide a clear meaning to the expectation of g(X(t)) conditioned on a realization of a particular (or group of) reaction channel(s).

The central idea proposed in the present work is to identify a realization of the reaction channels with the realization of the underlying vector of independent standard Poisson processes,

To this end, let us denote X(t, N(ω)) the solution of the Poisson process formulation of the stochastic dynamics given by (2), that is,

where

Then, defining F(N(ω))≐g(X(t, N(ω))), with X(t, N(ω)) given by (2), and assuming F has finite second-order moment, the variance decomposition with respect to the Poisson processes in N can proceed as discussed previously in Sec. III, because the K_r channels are associated with independent standard Poisson processes N_k(ω). To this end, let $K$ be the power set of {1, …, K_r}; for $u∈K$⁠, the conditional expectation of g(X(t, N)) given N_𝔲 = n_𝔲 becomes

Observe that given n_u, the physical times and sequence of firing reaction channels remain in general stochastic, that is, dependent on N_𝔲∼, even for the channels R_j with index j ∈ 𝔲, owing to the definition of the scaled time t_j: the conditioning fixes the sequence of firing channels only in the unscaled time. However, $V E g ( X ( t , N ) ) N u$ does measure the variability in g(X(t)) induced by the underlying stochastic processes N_j(ω) with j ∈ 𝔲, and so characterizes the variance caused by the inherent stochasticity in channels R_j, j ∈ 𝔲. As a result, the sensitivity indices defined above can be used to quantify and characterize the impact of individual or group of reaction channels on the variability of g(X(t)). Even though the indices are not characterizing a sensitivity with respect to some changes in the definition of the propensity functions associated to the channels (case of parametric studies), but rather characterize a contribution to the inherent stochasticity of the system, we shall continue to refer to them as sensitivity indices.

We complete the section by discussing the practical implementation of the Monte Carlo procedure for the estimation of the sensitivity indices of g(X(t, N)), and drop the time dependence of X for simplicity.

In addition to the classical mean and variance estimations for g(X(N)), the Monte Carlo estimation of sensitivity indices involves correlations of the form (see (13) and (14))

where N^I,(i) and N^II,(i) are elements of two independent sample sets of N(ω).

Direct sampling of the vector of K_r independent standard Poisson processes is generally not practical, as it would require the storage of a prohibitively large number of firing times for each channel. Instead, it is more convenient to construct the Poisson processes “on the fly,” as shown in Algorithm II, where the next unscaled firing time $τ j +$ for channel j is determined only when it actually fires (that is, when the unscaled time τ_j reaches $τ j +$⁠).

However, for the correct MC estimation of the correlation in (16), it is crucial to ensure that the sequence of (unscaled) firing times $τ j +$ remains the same for the channels j ∈ 𝔲 when computing $X ( N u I , ( i ) , N u ∼ I , ( i ) )$ and $X ( N u I , ( i ) , N u ∼ I I , ( i ) )$⁠. In other words, the same sequence of time increments between the successive firing events of channel j ∈ 𝔲 must be repeated; because the time increment between two successive firing events for channel j is simulated by means of the pseudo-random number generators RG_j, as shown in line 20 of Algorithm II, this implies that the generator RG_j must reproduce the same sequence of pseudo-random numbers r_j when computing $g ( X ( N u I , ( i ) , N u ∼ I , ( i ) ) )$ and $g ( X ( N u I , ( i ) , N u ∼ I I , ( i ) ) )$⁠, whenever j ∈ 𝔲. On the contrary, for a channel index j ∉ 𝔲, two independent sequences of random numbers r_j have to be used. The control of sequence of pseudo-random numbers delivered by RG_j can be classically enforced through the “seed” of the generator. In this context, the Poisson process N_j(τ) is totally determined by RG_j and its seed, s_j. For such generators, consistent solutions,

are obtained by calling Algorithm II twice, with the same initialization (seeding) RG_j(s_j) if j ∈ 𝔲, and two random and independent realizations $RG j ( s j I )$ and $RG j ( s j I I )$ for j ∈ 𝔲_∼.

The whole procedure for the Monte Carlo estimation of the whole set of K_r first and total-order sensitivity indices S_{i} and T_{i} is schematically illustrated in Algorithm III, in the case of g(X) = g(X(t = T)). Note that Algorithm III is not optimized; in particular the computationally intensive part of the algorithm, namely, the (K_r + 1) calls to NRA can be carried out in parallel. Nonetheless, the present formulation provides a general framework that can be easily adapted to the estimation of generic sensitivity indices associated to $u∈K$⁠, and to other functionals g(X).

In this section, we apply and illustrate the proposed sensitivity analysis and variance decomposition methods to well-known stochastic models.

The BD process²⁵ involves a single species S (M_s = 1) and K_r = 2 reaction channels,

with propensity functions

We set b = 200, d = 1, and use M = 1 000 000 Monte Carlo samples to compute the estimates presented in this section.

Figure 1 reports typical trajectories of the system for the initial condition X(t = 0) = 0. The plot shows the transient growth of X toward its asymptotic distribution, which is the Poisson distribution with mean b/d = 200. As depicted in Figure 2, the histogram of X(t = 8) is very close to the asymptotic distribution.

For this simple model with only two reaction channels, R_b and R_d, the analysis of the variance for g(X(t)) = X(t) yields only three terms corresponding to the two first-order sensitivity indices S_{b} and S_{d}, and the mixed contribution $S { b , d } ′ ≐V X { b , d } /V X$⁠; we have

and $S { b } + S { d } + S { b , d } ′ =1$⁠.

Figure 3 shows the evolution of the first-order and total sensitivity indices of X for t ∈ [0, 8]. To appreciate the transient dynamic, the sensitivity indices have been scaled by the variance $V X$ whose evolution is also provided. For the each reaction channel, the scaled first-order and total sensitivity indices are reported as a colored strip (red for R_b and blue for R_d) spanning vertically the interval $[ V X S , V X T ]$⁠, so the extent of the strips is equal to $V X { b , d }$⁠. The results show that at early times during the transient (t < 1) the variability in X is predominantly caused by the birth channel stochasticity, whereas the death channel induces a variability that grows in time at a much lower rate. This can be explained by the fact that R_b is a zero-order reaction with constant rate, whereas R_d is a first-order reaction and thus requires a sufficiently large population X to fire at a significant rate. Further, during this transient stage, the dynamics are essentially additive since S_{b} + S_{d} ≈ 1. For 1 ≤ t ≤ 4, the variability induced by R_d only, measured by S_{d}, continues to grow with the population size (recall that $E X =V X$ for this model) and becomes eventually close to, but lower than S_{b}. In the mean time, the mixed effects also increase, as reflected by the growing extent of the strips. Finally, for t > 4, interval for which the variance has essentially achieved its asymptotic value (⁠$V X ∼200$⁠), the first order sensitivity index S_{b} becomes constant while the mixed effect S_{d} starts to decrease with a very low rate. The slow decay of S_{d} underlines the slow emergence of an increasingly important mixed contribution of the two channels on the variability of X.

To illustrate the channels interaction, we provide in Figure 4 the time evolution up to t = 500 of the first-order sensitivity indices S_{b} and S_{d}, together with the interaction term $S { b , d } ′ =1− S { b } − S { d }$⁠. The figure shows that while the variability induced by R_b remains constant for large t and amounts to 50% of the overall variance, the sole impact of channel R_d monotonically decays to the benefit of the mixed term. Again, this behavior is explained by the different natures of the two channels. First, because channel R_b is a zero-order reaction, the stochasticity that it induces is independent of the state X and therefore becomes asymptotically constant, accounting for half the variance. Second, Eqs. (2) and (3) indicate that the contribution of channel R_d (a first-order reaction) to the dynamics of X, involves not only the intrinsic stochasticity of the standardized Poisson process N_d, but also the variability of the scaled-time t_d which involves the path-integral of X(t). As a result, the variance due to N_d(t_d) involves also the variability brought by channel R_b, through its effect on X(t), and our computation quantifies how the sole effect of the variability of N_d monotonically decays while the interaction of the two channels grows.

The Schlögl model²⁶ has K_r = 4 reaction channels, expressed according to

Species B₁ and B₂ are assumed in large excess with constant population over time: X_B1 = X_B2/2 = 10⁵. The stochastic state then reduces to the single evolving species S with M_s = 1. The propensity functions are given by

and

We set c₁ = 3 × 10⁻⁷, c₂ = 10⁻⁴, c₃ = 10⁻³, and c₄ = 3.5, so that the dynamics involve a zero-order reaction (channel R₃), a first-order reaction (channel R₄), a second-order reaction (channel R₁), and a third-order reaction (channel R₂). In addition, we use the deterministic initial condition X(t = 0) = 250.

For the present settings, the system exhibits a bifurcation with two attracting branches. This is illustrated in Figures 5 and 6, which depict selected trajectories of X(t) and the empirical histogram of X(t = 8), respectively. The bi-modal character of the solution at t = 8 is clearly evidenced, with a first peak at x ∼ 100 that is well separated from a second peak occurring at x ∼ 600. Also note the different spreads of the distributions, as the second peak is significantly broader than the first.

We first examine sensitivity indices of the K_r reaction channels on g(X(t)) = X(t). For each channel, Figure 7 shows the ranges [S_{j}(t), T_{j}(t)] scaled using the state variance $V X ( t )$ as scaling factor. Also plotted is the sum of the first-order contributions. Focusing first on the first-order indices, we observe that for t ≤ 8, we have

Since most of the variance arises from the bi-modality of the solution, this suggests that channels R₁ and R₄ are the primary contributors to the solution branch selection. Figure 7 also indicates that up to t ≈ 2 the system dynamics are essentially additive, while at later times mixed interactions between channels become increasingly significant. In fact, it is seen that the variance induced by the first-order effects, $V X × ∑ j = 1 K r S { j }$⁠, i.e., by the individual channels without interaction, is essentially constant for t > 4 while the total variance continues to grow as interaction contributions develop. At t = 8, X is close to its asymptotic distribution, and the importance of the mixed effects in all the channels, T_{j} − S_{j}, can be appreciated from the extents of the colored strips in the plots. Again, mixed effects for channels R₁ and R₄ appear to be comparable, but significantly larger than they are in channels R₂ and R₃.

To gain further insight into the structure of the variability, the characterization of high-order partial variances is needed. For the present model with K = 4, we can afford to estimate all first-order sensitivity indices S_𝔲, $u∈K$⁠, from which the partial variances $V X u$ can be retrieved. We provide in Figure 8 the second-order partial variances $V X u$ for |𝔲| = 2. It is seen that they all assume relatively low values (≲500) except for $V X { 1 , 4 }$ which asymptotically peaks at approximately 3500, confirming the preponderance of these two channels. Note also that $V X { 2 , 3 }$ is estimated to be zero, so not only these two channels have lower first-order variance, but also have no second-order interaction.

Regarding the higher-order partial variances for |𝔲| > 2, reported in Figure 9, we observe again that $V X u$ have significant values only for {1, 4} ∈ 𝔲. The partial variance $V X { 1 , 2 , 3 , 4 }$⁠, which quantifies the variability due to the interaction between all the channels, becomes asymptotically dominant over all other partial variances except $V X { 1 , 2 , 4 }$ and $V X { 1 , 3 , 4 }$⁠, underscoring the importance of complex high-order interactions between the channels. Further, unlike the third-order partial variances, $V X { 1 , 2 , 3 , 4 }$ has not yet converged at t = 8, highlighting the slower time scales at which full interactions develop.