Stochastic models of chemical systems are often subjected to uncertainties in kinetic parameters in addition to the inherent random nature of their dynamics. Uncertainty quantification in such systems is generally achieved by means of sensitivity analyses in which one characterizes the variability with the uncertain kinetic parameters of the first statistical moments of model predictions. In this work, we propose an original global sensitivity analysis method where the parametric and inherent variability sources are both treated through Sobol’s decomposition of the variance into contributions from arbitrary subset of uncertain parameters *and* stochastic reaction channels. The conceptual development only assumes that the inherent and parametric sources are independent, and considers the Poisson processes in the random-time-change representation of the state dynamics as the fundamental objects governing the inherent stochasticity. A sampling algorithm is proposed to perform the global sensitivity analysis, and to estimate the partial variances and sensitivity indices characterizing the importance of the various sources of variability and their interactions. The birth-death and Schlögl models are used to illustrate both the implementation of the algorithm and the richness of the proposed analysis method. The output of the proposed sensitivity analysis is also contrasted with a local derivative-based sensitivity analysis method classically used for this type of systems.

## I. INTRODUCTION

Stochastic models are widely used in many scientific fields, including population dynamics, queues and social-network studies, biosciences, and chemical systems. Stochastic modeling is necessary when deterministic models are not able to capture important features of the dynamics, for instance, when large-scale effects of unresolved small-scale fluctuations occur, or when systems are subjected to important inherent noise. Regarding chemical modeling, which is the focus of the present work, deterministic models are usually unsuitable for describing the evolution of small systems where reactive species are present in small amounts, such that stochastic fluctuations have great effects on the system response and carry important meaning.

Chemical system modeling requires the selection of reaction mechanisms (or pathways), that is, the set of reactions between the species constituting the system, and the definition of propensity (or intensity) functions that relate the current state of the system to the probability of occurrence of a reaction in the next infinitesimal time interval. Algorithms and stochastic simulators can then be designed for the generation of paths, or trajectories, of the system’s stochastic dynamics. Sample sets of paths can be used to compute, *à la* Monte Carlo, the statistical moments of some model output or, more generally, the expected value of some functional of the stochastic state.

Typically, the propensity functions involve kinetic parameters, such as thermodynamic variables and reaction parameters, e.g., the mixture temperature and pressure, the rate constants, the pre-exponential factors, and the activation energies of the Arrhenius law. Often, these kinetic parameters are not well known and must be considered as uncertain. In this case, it is critical to assess the impact of the uncertain kinetic parameters on the model predictions.^{1} Classically, this is achieved by performing a sensitivity analysis (SA), which, roughly speaking, measures the variability of a model prediction when the uncertain parameters are varied.

A number of SA methods have been proposed for uncertain chemical systems.^{2,3} Restricting the discussion to the case of stochastic models, SA methods have been mostly limited to the assessment of parametric variability in the first statistical moments (mean and variance) of the model output due to the uncertain kinetic parameters. Along these lines, one can distinguish between the local and global SA methods. The former essentially consists in the characterization of the variability in the statistical moments, based on their derivatives with respect to the kinetic parameters at some nominal values (see, for instance, Refs. 4–7).

Global SA (GSA) methods, on the other hand, aim at accounting for the uncertainty effects when the parameters vary over their whole range. The GSA methods are often based on second-order characterizations of the variability (with respect to the parameters) using the so-called Sobol or ANOVA decompositions.^{8,9} The Sobol decomposition can be determined efficiently relying on a projection step, where the dependences of the moments are first represented as Polynomial Chaos (PC) expansions,^{10,11} before computing the partial variances and so-called sensitivity indices.^{12,13}

More recently, we have proposed in Ref. 14 a different approach for GSA, for parametric uncertainties in stochastic systems. Compared to the classical approaches discussed above, the idea is to treat the inherent stochasticity of the system as an additional and independent source of uncertainty in the Sobol decomposition. Within this framework, one proceeds directly with the decomposition of the stochastic model output, and not with the decomposition of its statistical moments. This paradigm shift enables a much more detailed and finer GSA, such as the assessment of the partial variances associated with the interactions between parameters and inherent stochasticity, a task that is not possible for a GSA of the moments.

In the GSA method proposed in Refs. 14 and 15, the decomposition of the variance was performed by relying on PC expansions of the individual random trajectories. This was possible owing to the smoothness (differentiability) of the stochastic trajectories with respect to the parameters. Such path-wise smoothness is not to be expected in general for the systems considered in the present work, either because of strong non-linearities and bifurcations or because of the discrete nature of the stochastic state spaces (typically, the system state is integer valued). Another important challenge in extending the GSA method of Refs. 14 and 15 to stochastic simulators, as explained later, comes from the need to perform conditional re-sampling of the trajectories. Resampling trajectories conditioned on fixed kinetic parameter values is straightforward in the context of stochastic simulators. In contrast, the meaning of resampling trajectories conditioned on fixed “values” of the inherent stochasticity must be clarified and properly defined, prior to be implemented in a stochastic simulator.

This issue was recently addressed in Ref. 16, where it has been proposed to identify a particular realization of the inherent stochasticity with the corresponding realization of the standard Poisson processes appearing the so-called random-time-change representation of the dynamics.^{17} Adopting this definition, we were able, in Ref. 16, to perform an analysis of the relative influences of different reaction channels on the total variability of some model outputs of a stochastic simulator.

In the present paper, we generalize the methodology in Ref. 16 to propose a GSA method for stochastic simulators having additionally uncertain parameters. In doing so, we aim to preserve the ability to perform a global variance decomposition, namely, one that is not simply based on the SA of statistical moments. Achieving this new capability for stochastic simulators requires us to devise effective means to distinguish between the effects of parametric uncertainty and state-dependent stochastic variability and to quantify their joint contributions. The presently proposed methodology enables us to overcome associated challenges and in addition offers several key features. First, the method is global and can be applied to models having arbitrary levels of parametric uncertainty as well as inherent variability. Second, the method enables a detailed GSA, through the quantification of the partial variances due to individual parameters or groups of parameters, randomness in individual or groups of the channel dynamics (inherent stochasticity), and finally the interaction between arbitrary groupings of parameters and random channels. Finally, the method is purely based on the conditional sampling of the model and is then easily amenable to stochastic simulators relying on a random-time-change formulation as well as being trivially implemented in parallel.

The outline of the paper is as follows. In Sec. II, we briefly discuss the stochastic simulators and introduce the modified next reaction method, which is subsequently used to perform the conditional sampling of the stochastic model. Sec. III details the proposed GSA method, starting with a discussion of the parametric and inherent variabilities in Sec. III A. The Sobol-Hoeffding variance decomposition and the associated sensitivity indices are introduced in Sec. III B, starting from the simple separation of parametric uncertainty and random channels variability and leading to the computation of general sensitivity indices. A numerical algorithm performing the GSA and its practical implementation are described in Sec. III C. The proposed method is illustrated based on two simple stochastic models, namely, the birth-death process in Sec. IV, and the Schlögl model in Sec. V. The results of the GSA are contrasted with a local-derivative-based SA in order to highlight the robustness and richness of the proposed method. Concluding remarks are provided in Sec. VI.

## II. STOCHASTIC CHEMICAL SYSTEMS AND SIMULATOR

In this section, we review background materials on the stochastic modeling of chemical systems and relevant simulation methods. In Sec. II B, we introduce the Modified Next Reaction Method (MNRM) which is the basis of subsequent developments.

### A. Stochastic chemical systems

Let $(\Omega ,\Sigma ,\mu )$ be a probability space, where Ω is the event space, Σ denotes a *σ*-algebra on Ω, and *μ* is a probability measure. We consider a chemical reaction system involving $Ms$ chemical species, labeled $S1,\u2026,SMs$, whose evolutions in continuous time are discrete and stochastic. The random state of the system is $\bm{X}:[0,T]\xd7\Omega \u2192\mathbb{N}0Ms$, where the $i$th component $Xi(t)$ is the number of molecules of species $Si$ at time $t$. The state vector $\bm{X}$ evolves through $Mr$ reactions channels, $R1,\u2026,RMr$, each characterized by a propensity or intensity function $aj$ and a stoichiometric vector $\bm{\nu}j$. The stoichiometric vector $\bm{\nu}j\u2208\mathbb{Z}Ms$ prescribes the changes in the population of the reactant species due to the occurrence of a reaction in channel $Rj$ and are therefore independent of the state and time. The propensity functions, on the contrary, are dependent on the system state and thermodynamical conditions and, therefore, they change in time. For simplicity, we shall restrict ourselves to well-stirred systems at constant temperature and pressure, such that the propensity functions can be expressed as $aj(\bm{X}(t);\bm{c})$, where $\bm{c}=(c1,\u2026,cNc)\u2208\mathbb{R}Nc$ is a vector of real valued kinetic parameters which are not time dependent. For $\tau \u226a1$, $aj(\bm{X}(t);\bm{c})\tau $ is the probability that, given the current system state $\bm{X}(t)$, a reaction along channel $Rj$ occurs during the small time interval $(t,t+\tau ]$. In other words, the random process $\bm{X}(t)$ is a continuous-time Markov chain defined by the set of transition probabilities for $j=1,\u2026,Mr$,

Note that the propensity functions and stoichiometric vectors of the reaction channels must satisfy a non-negativity assumption, meaning that the dynamic cannot produce a population with a negative size.

The Chemical Master Equation^{18,19} (CME) describes the evolution of the probability distribution of the stochastic state vector $\bm{X}(t)$; it is expressed as

In the CME (2), we denoted $Pt,t0[\bm{x}|\bm{x}0]$ as the probability of the system to reach a state $\bm{x}$ at time $t\u2265t0$, knowing that it is in a state $\bm{x}0$ at time $t0$. Equation (2) has no general explicit solution, and its numerical solution is too costly to compute in most cases, especially when $Ms$ exceeds a few units. One has then resort to a simulation method to generate sample trajectories from the distribution that is the solution of Eq. (2) without having to explicitly solve the CME. One classically uses so-called exact simulation algorithms to compute such trajectories and estimate the expected value of some functional, $G$, of the random system state, through

where $\omega \u2208\Omega $ denotes a particular random event and $\bm{X}(i)$ are sampled trajectories of *𝑿*. Consequently, relying on a simulator or solving the CME is equivalent, although using stochastic simulators has an inherent sampling error as a finite number *M* of paths can be generated. Note that there also exist inexact stochastic simulation algorithms that produce trajectories according to an approximation of the CME solution. Such approximate simulators are not considered in this work.

### B. Modified next reaction method (MNRM)

One of the best-known exact simulation methods is the Stochastic Simulation Algorithm (SSA) which was first introduced by Gillespie.^{20,21} Several modifications and variants of the original SSA have been proposed to improve its computational efficiency and adapt the algorithm to specific properties and characteristics of the simulated system (see, for instance, Refs. 22 and 23). For the purpose of decomposing the variance, as pursued below, a fine distinction and control of the individual channels’ dynamics is required.^{16} The standard SSA is not suitable for the purpose of performing the variance decomposition because the algorithm combines all channels’ dynamics to determine the next reaction time and the first firing channel. For this reason, we consider the Modified Next Reaction Method^{24} (MNRM) (see Algorithm 1) which is a slightly modified version of the original Next Reaction Method (NRM).^{25} The MNRM is an exact simulation algorithm making explicit the respective contributions of the $Mr$ channels in the evolution of the system state. This characteristic is subsequently exploited to perform the variance decomposition, as discussed in Sec. III.

The algorithm uses Kurtz’ random time-change (RTC) representation^{17} to simulate exact paths, expressing the random state $\bm{X}(t)$ through

where $Yj$ are independent unit-rate Poisson processes with respective scaled times $tj$ given by

Equation (3) shows that the change in the state, from the initial condition $\bm{X}(t0)$, is proportional to the number of times ($Yj$) that each channel has fired, times their respective stoichiometric vectors $\bm{\nu}j$. The random number of times a channel has fired is identified as the value of a standard (unit rate) Poisson process, evaluated at a scaled time given by the integral of the propensity (or intensity) function along the trajectory. Hence, in the MNRM, the randomness is encapsulated in the Poisson processes describing the channels’ dynamics. This crucial point was exploited in our previous paper,^{16} where control of the realizations of the Poisson processes was exploited to perform a variance decomposition and analyze the variability induced by different reaction channels.

An algorithm for the simulation of the stochastic evolution governed by Eq. (3) can be easily set up, remarking that both the system state $\bm{X}(s)$ and the propensity functions $aj(\bm{X}(s);\bm{c})$ remain constant between successive reactions. Denote $tj+>tj$ as the scaled time at which $Yj$ jumps from $Yj(tj)$ to $Yj(tj)+1$; the main idea of Algorithm 1 is that, given the current state $\bm{X}(t)$, scaled times $tj$ and paths of Poisson processes $Yj$, we can determine from (4) the (physical) time increments $dtj=(tj+\u2212tj)/aj$ till $Yj(tj)$ jumps to $Yj+1$. The minimal $dtj$ hence determines both the time and reaction channel that will first fire. Therefore, it suffices to keep track of the set of $tj+$ and $tj$ to simulate trajectories of the system, as shown in Algorithm 1.

. | . |
---|---|

Function MNRM$(\bm{X}0,T,{\bm{\nu}j},{aj},\bm{c},\mathbf{R}\mathbf{G})$ | |

Input: Initial condition $\bm{X}0$, final time $T$, stoichiometric vector ${\bm{\nu}j}$, | |

propensity functions ${aj}$, vector of kinetic parameters $\bm{c}$, | |

and initialize pseudo-random number | |

generators $\mathbf{R}\mathbf{G}=(RG1,\u2026,RGMr)$ | |

Output: State $\bm{X}(T)$ | |

1: for $j=1$ to $Mr$ do | |

2: Draw $rj$ from $RGj$ | |

3: $tj\u21900$, $tj+\u2190\u2212log(rj)$ {set next reaction times} | |

4: end for | |

5: $t\u21900$, $\bm{X}\u2190\bm{X}0$ | |

6: loop | |

7:for $j=1$ to $Mr$ do | |

8: Evaluate $aj(\bm{X};\bm{c})$ and $dtj=tj+\u2212tjaj$ | |

9:end for | |

10: Set $l=argminjdtj$ {pick next reaction} | |

11:if $t+dtl>T$ then | |

12: Break {final time reached} | |

13:else | |

14: $\bm{X}\u2190\bm{X}+\bm{\nu}l$ {update the state vector} | |

15: $t\u2190t+dtl$ {advance time} | |

16: for $j=1$ to $Mr$ do | |

17: $tj\u2190tj+ajdtl$ {update scaled times} | |

18: end for | |

19: Get $rl$ from $RGl$ and set $tl+\u2190tl+\u2212log(rl)$ {next | |

reaction time} | |

20: end if | |

21: end loop |

. | . |
---|---|

Function MNRM$(\bm{X}0,T,{\bm{\nu}j},{aj},\bm{c},\mathbf{R}\mathbf{G})$ | |

Input: Initial condition $\bm{X}0$, final time $T$, stoichiometric vector ${\bm{\nu}j}$, | |

propensity functions ${aj}$, vector of kinetic parameters $\bm{c}$, | |

and initialize pseudo-random number | |

generators $\mathbf{R}\mathbf{G}=(RG1,\u2026,RGMr)$ | |

Output: State $\bm{X}(T)$ | |

1: for $j=1$ to $Mr$ do | |

2: Draw $rj$ from $RGj$ | |

3: $tj\u21900$, $tj+\u2190\u2212log(rj)$ {set next reaction times} | |

4: end for | |

5: $t\u21900$, $\bm{X}\u2190\bm{X}0$ | |

6: loop | |

7:for $j=1$ to $Mr$ do | |

8: Evaluate $aj(\bm{X};\bm{c})$ and $dtj=tj+\u2212tjaj$ | |

9:end for | |

10: Set $l=argminjdtj$ {pick next reaction} | |

11:if $t+dtl>T$ then | |

12: Break {final time reached} | |

13:else | |

14: $\bm{X}\u2190\bm{X}+\bm{\nu}l$ {update the state vector} | |

15: $t\u2190t+dtl$ {advance time} | |

16: for $j=1$ to $Mr$ do | |

17: $tj\u2190tj+ajdtl$ {update scaled times} | |

18: end for | |

19: Get $rl$ from $RGl$ and set $tl+\u2190tl+\u2212log(rl)$ {next | |

reaction time} | |

20: end if | |

21: end loop |

We remark that the time-marching method depicted in Algorithm 1 does not require the storage of all the sequences of jumping times $tj+$. Instead, the sequences are determined on the fly using the random number generators (see line 19). In practice, the sequences of jumping times are entirely determined by the seeds initializing the pseudo-random number generators $RGj$, and the same sequences of random (unscaled) jumping times can be reproduced by reseeding the pseudo-random number generators accordingly. Another important remark concerning Algorithm 1 is that it considers $Mr$ distinct random number generators $RGj$, one for each reaction channel. This may appear unnecessary, in view of the advancement in time performed in Algorithm 1 as a single random number generator would be enough. In fact, the use of different random number generators allows for the control of the individual Poisson processes and consequently enables the variance decomposition which requires reproducing exact sequences of jumping times, as explained below. To this end, the possibility of identifying the sequence of jumping times for a channel to the random seed of its generator is critical for the computational efficiency because it alleviates the burden of storing sequences of sampled $tj$.

## III. GLOBAL SENSITIVITY ANALYSIS

In this section, we introduce the variance decomposition method and the related global sensitivity analysis. As discussed in the introduction, our method differs from the classical sensitivity analysis Refs. 4–7 of stochastic systems that are usually based on the estimation of the local sensitivity of selected statistical moments with respect to uncertain system parameters. Our method, on the contrary, aims at providing a direct second-order (variance-based) characterization of the respective variabilities induced by the uncertain parameters, globally, and quantifying separately the inherent stochasticity due to the channels. To this end, we build the approach originally proposed in Ref. 16, which focused on the variance decomposition over the reaction channels without considering parametric uncertainty.

We start by discussing parametric uncertainties in the context of stochastic reaction networks in Sec. III A. The decomposition of the variance and the uncertainty characterization using sensitivity indices are detailed in Sec. III B. Finally, in Sec. III C, we describe a Monte Carlo sampling approach for the estimation of the sensitivity indices and we provide a brief outline of the algorithm and implementation used in the present work.

### A. Parametric and stochastic uncertainties

As mentioned before, stochastic systems can involve two sources of randomness. On the one hand, the random dynamics associated with the reaction channels lead to inherent variability that cannot be reduced. On the other hand, the definition of the reaction networks involves several kinetic parameters that are usually estimated within a finite confidence level only. To properly apprehend stochastic simulations, it is therefore important to assess the respective impacts of these inherent and parametric uncertainties on the variability of a quantity of interest, $G(X)$, related to the system output.

A typical situation of parametric uncertainty, considered in this work, is the case of propensity functions $aj$ not perfectly known, but involving some uncertain parameters. Practical examples of uncertainties in the propensity functions will be given in the following sections; we can already mention the case of uncertainties in the reacting mixture temperature and pressure, in the Arrhenius activation energies, rate constants, and pre-exponential factors of the different reaction channels. In these situations, it is relevant to treat the uncertain parameters as random quantities. To this end, we shall denote by $\bm{C}(\omega )\u2208\mathbb{R}Nc$ the uncertain vector of random parameters in the propensity functions $aj$ (recall that $\omega \u2208\Omega $ is a random event). We shall also denote $p\bm{C}:\mathbb{R}Nc\u21a6\mathbb{R}+$ as the probability density function of *𝑪*. As a result of the uncertainty in the parameters of $aj$, the dynamics given by Eq. (3) are random, because of the stochastic nature of the Poisson processes $Yj$ and also because of the scaled times in Eq. (4) which involve propensity functions with now random coefficients *𝑪*. Specifically, in the case of random parameters *𝑪*, we have

In the following, we shall refer collectively to the set of $Nc$ random parameters *𝑪*, as the uncertain kinetic parameters. We shall also assume that the components of *𝑪* are independent random variables. Consequently, the scaled times $tj$ are given by

For convenience, we shall write the stochastic state as $\bm{X}(t,\bm{Y},\bm{C})$ to underline its dependence on the random Poisson processes *𝒀* and random parameters *𝑪*. The main difference between the two sources of variability (inherent and parametric) is that while the channel variability cannot be reduced, the parametric uncertainty can eventually be reduced by a better estimation of the propensity function parameters. A fundamental assumption of the proposed method is that *𝒀* and *𝑪* are independent. This assumption will allow us to decompose the variance as outlined below; it also implies that the inherent stochastic nature of the system is not dependent on the kinetic parameters. It does not mean, however, that the stochastic dynamics are not dependent on the parameters. In the remainder of the section, we discuss how to separate and quantify the variance associated with the two sources of randomness.

### B. Variance decomposition and sensitivity indices

We begin by introducing the Sobol-Hoeffding (SH) decomposition^{8} of a functional of several independent random variables. The decomposition is further exploited to derive an expression for the sensitivity indices, which will provide us with a measure of the impact of the distinct sources of uncertainty on the process.

Recall that we consider *𝒀* and *𝑪* as independent random quantities defined on the probability space $(\Omega ,\Sigma ,\mu )$. Let $G:\Omega \u2192G(\bm{X}(\bm{Y};\bm{C}))\u2208\mathbb{R}$ be a second-order random functional of the stochastic process, i.e., $\mathbb{E}[G(\bm{X}(\bm{Y};\bm{C}))2]<\u221e$. With a slight abuse of notation, we shall denote $G(\bm{X}(\bm{Y};\bm{C}))$ simply as $G(\bm{Y},\bm{C})$. The SH decomposition separates the function $G(\bm{Y},\bm{C})$ into independent and orthogonal elements as follows:

where $G\xaf$ is a constant, $Gpar(\bm{C})$ is a random functional that depends only on the parameters, $Gch(\bm{Y})$ depends only on the Poisson processes, and $Gmix(\bm{Y},\bm{C})$ depends on both random inputs. The pairwise orthogonality between the functions on the right-hand-side of Eq. (5) ensures the uniqueness of the decomposition. These functionals are given by $G\xaf=\mathbb{E}[G]$,

In the previous expressions, the conditional expectations are defined by

where *𝒄* and *𝒚* are deterministic realizations of the parameters and the Poisson processes, respectively. Because the decomposition in Eq. (5) is orthogonal, the variance of $G(\bm{Y},\bm{C})$ can be split according to

where $Vpar=\mathbb{V}[Gpar]$, $Vch=\mathbb{V}[Gch]$, and $Vmix=\mathbb{V}[Gmix]$ are called the partial variances. These partial variances account for the variability due to parameters, the reaction channels, and the interactions between them, respectively. Finally, the Sobol sensitivity indices are obtained by normalizing partial variances,

The indices $Spar$ and $Sch$ are the fractions of the variance in $G$ explained by the parametric uncertainty and stochastic channel dynamics, respectively. The index $Smix$ accounts for the fraction of the variance due to the interactions between the parametric and channel variabilities. To measure the total impacts of parametric uncertainty and channel stochasticity, we introduce the total indices

The decomposition of the variance of $G(\bm{Y},\bm{C})$ above can be generalized to perform finer sensitivity analyses and measure, for instance, the impact of a subset of reaction channels, of a single parameter, or some joint effects of channel(s) and parameter(s), provided that the components of *𝑪* are mutually independent (recall that the Poisson processes $Yj$ of different channels are independent). In this case, the random vector $\bm{U}=(Y1,\u2026,YMr,C1,\u2026,CNc)$ has components $Ui(\omega )$ which are all mutually independent. We shall set$NU=Mr+Nc$ as the total number of variability sources, defined as the sum of the number of reaction channels and uncertain kinetic parameters. With a slight abuse of notation, we shall now denote $G(\bm{X}(\bm{Y};\bm{C}))$ as $G(\bm{U})$. The general SH decomposition can immediately be generalized by considering the variability sources in the random vector *𝑼*. Let us denote the index set $U={1,\u2026,NU}$ and $PU$ its power set. For $u\u2286U$, let $u\u223c=U\u$ be its complement in U. With these notations, the generalized SH decomposition of *G* is

where $Gu(\bm{U}u)$ are mutually orthogonal random functionals, each depending only on the subvector of *𝑼* having components $\bm{U}i$ with $i\u2208u$. These functionals are recursively defined by

with the conditional expectations given by

Due to the orthogonality of the SH decomposition in Eq. (7), and observing that $G\u2205=\mathbb{E}[G]$, the variance of $G(U)$ can be expressed as

Similarly, “first-order” and total sensitivity indices associated with variability sources in u can be defined by extending the previous expressions according to

### C. Monte Carlo algorithm

We now show how to estimate the sensitivity indices by means of the Monte Carlo sampling and we briefly discuss the practical implementation of the sampling method.

#### 1. Monte Carlo estimation of the sensitivity indices

In this work, we rely on the method of Saltelli^{26} (see also Ref. 9) to estimate the sensitivity indices. For simplicity, we restrict the exposition to the separation of the effects associated with the channels, on the one hand, and kinetic parameters, on the other hand. The task is then to estimate $Spar$, $Sch$, and $Smix$.

The Monte Carlo estimation makes use of independent sample sets for both the Poisson processes and the kinetic parameters. The estimation of the partial variances is derived from evaluations of the functional $G(\bm{Y},\bm{C})$ at selected sample points resulting from the appropriate combination of the sample sets.

Let $Y$ and $C$ denote two sample sets of $M$ independent realizations $\bm{Y}(i)$ and $\bm{C}(i)$ of the Poisson processes and kinetic parameters, respectively, and where the superscript $1\u2264i\u2264M$ refers to the sample index. As discussed at the end of Section II B, the sampling of the Poisson processes in $\bm{Y}$ simply amounts to sampling (independently) the seeds of the random number generators that yields the respective sequences of jumping times in $\bm{Y}(i)$. On the other hand, the sampling of $\bm{C}$ amounts to drawing random vectors $\bm{C}(i)$ from the probability density function $p\bm{C}$. Standard methods, such as probabilistic transformations, can be used to generate $\bm{C}(i)$. From these sample sets, the Monte Carlo estimators of the expectation and variance of *G* are

Consider now independent replicas of the sample sets $Y$ and $C$, denoted, respectively, $Y\u223c$ and $C\u223c$, with elements $\bm{Y}\u223c(i)$ and $\bm{C}\u223c(i)$, respectively, defined as discussed previously. The original sample sets and their replicas can be appropriately combined to obtain Monte Carlo estimates $V^par$ and $V^ch$ of the partial variances of $G$; we use

The sensitivity indices are finally estimated as follows:

The MC estimators above can be extended to compute sensitivity indices $Su$ associated with different groupings u the input variables, for instance, the individual effects of a variable, $\bm{U}u=(Yj)$ or $\bm{U}u=(Cj)$. The estimate $Su^$ can be immediately achieved using the appropriate combination of the original sample set and replicas. For instance, the estimation of the partial variance $Vu$ associated with $Cj$ involves the estimator

where $\bm{C}k\u223c\u223c(i)=\bm{C}k\u223c(i)$ for $k\u2260j$ and $\bm{C}j\u223c\u223c(i)=\bm{C}j(i)$. We estimate the total sensitivity indices $Tu$ using the relation $Tu=1\u2212Su$ and the Jansen formula^{27} which, while also based on combining two replicas of sample sets, was shown to be more accurate.^{28}

#### 2. Implementation details

We have just seen that the Monte Carlo estimation of the sensitivity indices involves multiple evaluations of the quantity of interest for suitable combinations of sample sets of the Poisson and kinetic parameters. Algorithm 2 provides a sketch of a practical implementation in the case of the simple channel and parametric separation. The algorithm can be easily adapted to handle the estimation of more general sensitivity indices $Su$ and $Tu$. It is adapted from the original algorithm proposed in Ref. 16 for the separation of the variance due to different channels (without parametric uncertainty). Specifically, instead of combining different sample sets for the individual Poisson processes $Yj$, Algorithm 2 combines different sample sets of the vector of Poisson processes *𝒀* and parameter *𝑪*.

. | . |
---|---|

Function MNRM-GSA$(M,\bm{X}0,T,{aj},{\bm{\nu}j},\bm{C})$ | |

Input: Number $M$ of samples, initial condition $\bm{X}0$, final time | |

$T$, propensity functions ${aj}$ and stoichiometric vector ${\bm{\nu}j}$, | |

random kinetic parameters $\bm{C}$. | |

Output: $S^par$, $S^ch$ and $S^mix$ | |

1: $EG^\u21900$, $VG^\u21900$, {init mean and variance} | |

2: $V^par\u21900$, $V^ch\u21900$ and $V^mix\u21900$ {init partial var} | |

3: for $i=1$ to $M$ do | |

4: Draw seeds $\bm{S}$ and replicas $\bm{S}\u223c$ for the random number | |

generators | |

5: Draw $\bm{c}$ and replicas $\bm{c}\u223c$ from density $p\bm{C}$ | |

6: $\bm{X}\u2190\U0001d67c\U0001d67d\U0001d681\U0001d67c(\bm{X}0,T,{\bm{\nu}j},{aj},\bm{c},\mathbf{R}\mathbf{G}(\bm{S}))$ | |

7: $\bm{X}ch\u2190\U0001d67c\U0001d67d\U0001d681\U0001d67c(\bm{X}0,T,{\bm{\nu}j},{aj},\bm{c}\u223c,\mathbf{R}\mathbf{G}(\bm{S}))$ | |

8: $\bm{X}par\u2190\U0001d67c\U0001d67d\U0001d681\U0001d67c(\bm{X}0,T,{\bm{\nu}j},{aj},\bm{c},\mathbf{R}\mathbf{G}(\bm{S}\u223c))$ | |

9: $EG^\u2190EG^+G(\bm{X})$ | |

10: $VG^\u2190VG^+G(\bm{X})2$ | |

11: $Vch^\u2190Vch^+G(\bm{X})G(\bm{X}ch)$ | |

12: $Vpar^\u2190Vpar^+G(\bm{X})G(\bm{X}par)$ | |

13: end for | |

14: $EG^\u2190EG^/M$ | |

15: $VG^\u2190VG^/(M\u22121)\u2212EG^2$ | |

16: $Vpar^\u2190Vpar^/(M\u22121)\u2212EG^2$, $Vch^\u2190Vch^/(M\u22121)\u2212EG^2$ | |

17: $Spar^=Vpar^/VG^$, $Sch^=Vch^/VG^$ and $Smix^=1\u2212Spar^\u2212Sch^$ |

. | . |
---|---|

Function MNRM-GSA$(M,\bm{X}0,T,{aj},{\bm{\nu}j},\bm{C})$ | |

Input: Number $M$ of samples, initial condition $\bm{X}0$, final time | |

$T$, propensity functions ${aj}$ and stoichiometric vector ${\bm{\nu}j}$, | |

random kinetic parameters $\bm{C}$. | |

Output: $S^par$, $S^ch$ and $S^mix$ | |

1: $EG^\u21900$, $VG^\u21900$, {init mean and variance} | |

2: $V^par\u21900$, $V^ch\u21900$ and $V^mix\u21900$ {init partial var} | |

3: for $i=1$ to $M$ do | |

4: Draw seeds $\bm{S}$ and replicas $\bm{S}\u223c$ for the random number | |

generators | |

5: Draw $\bm{c}$ and replicas $\bm{c}\u223c$ from density $p\bm{C}$ | |

6: $\bm{X}\u2190\U0001d67c\U0001d67d\U0001d681\U0001d67c(\bm{X}0,T,{\bm{\nu}j},{aj},\bm{c},\mathbf{R}\mathbf{G}(\bm{S}))$ | |

7: $\bm{X}ch\u2190\U0001d67c\U0001d67d\U0001d681\U0001d67c(\bm{X}0,T,{\bm{\nu}j},{aj},\bm{c}\u223c,\mathbf{R}\mathbf{G}(\bm{S}))$ | |

8: $\bm{X}par\u2190\U0001d67c\U0001d67d\U0001d681\U0001d67c(\bm{X}0,T,{\bm{\nu}j},{aj},\bm{c},\mathbf{R}\mathbf{G}(\bm{S}\u223c))$ | |

9: $EG^\u2190EG^+G(\bm{X})$ | |

10: $VG^\u2190VG^+G(\bm{X})2$ | |

11: $Vch^\u2190Vch^+G(\bm{X})G(\bm{X}ch)$ | |

12: $Vpar^\u2190Vpar^+G(\bm{X})G(\bm{X}par)$ | |

13: end for | |

14: $EG^\u2190EG^/M$ | |

15: $VG^\u2190VG^/(M\u22121)\u2212EG^2$ | |

16: $Vpar^\u2190Vpar^/(M\u22121)\u2212EG^2$, $Vch^\u2190Vch^/(M\u22121)\u2212EG^2$ | |

17: $Spar^=Vpar^/VG^$, $Sch^=Vch^/VG^$ and $Smix^=1\u2212Spar^\u2212Sch^$ |

To this end, within the main loop of the MC sampler (loop over $i$), two vectors of seeds are first independently drawn (*𝑺* and its replica $\bm{S}^$, see line 4). These seeds are subsequently used to initialize the random number generators of the stochastic simulator, so as to appropriately set the sequence of jumping times of the Poisson processes and individual reaction channels. Second, two independent realizations of the parameters are drawn from the density $p\bm{C}$ ($\bm{c}$ and $\bm{c}\u223c$, see line 5). In the present work, we shall consider independent $Cj$ with uniform distributions, such that their sampling is straightforward.

The trajectories $X(\bm{Y},\bm{c})$, $X(\bm{Y},\bm{c}\u223c)$, and $X(\bm{Y}\u223c,\bm{c}\u223c)$ are then determined using the MNRM Algorithm 1 (see lines 6-8); the statistic and correlations of the functional $G(X)$ between these trajectories are subsequently accumulated, before proceeding to the next sample. Finally, the sensitivity indices of *G* are estimated when *M* samples have been computed.

We remark that the procedure depicted in Algorithm 2 aims at clarifying the sampling procedure and is not designed in view of optimizing the computational efficiency. In particular, it is possible to parallelize the main loop of the algorithm.

## IV. BIRTH-DEATH EXAMPLE

### A. Birth-death model

The BD system models the dynamic of a single species *S* ($Ms=1$), evolving through $Mr=2$ reaction channels, simply indexed as $Rb$ and $Rd$, according to

The first reaction channel describes the spontaneous creation (birth) of individuals, at a time-rate *B*; the second channel describes the disappearance (death) of individual at a rate proportional to the population size times the time-rate *D*. In the following, we shall consider the birth and death rates as uncertain parameters, that is, $C1=B$ and $C2=D$, and proceed with their modeling as two independent random variables having uniform distributions. For simplicity, we shall directly denote *B* and *D* as the two uncertain parameters in the rest of this section. We denote $\mu B$ and $\sigma B$ (respectively, $\mu D$ and $\sigma D$) as the mean and standard deviation of the birth rate (respectively, death rate). The variability ranges of the two rates are characterized by their coefficients of variation, $CoVB=\sigma B/\mu B$ and $CoVD=\sigma D/\mu D$. For simplicity the two coefficients of variation will be set equal to $CoVB=CoVD=CoV$, such that the two channels have comparable amounts of parametric variability. Finally, the propensity functions are given by

By Eqs. (3) and (4), and from the stoichiometric coefficients associated with the BD model, the state can be rewritten as the difference between two Poisson processes,

where the scaled times are given by

Unless specified otherwise, we set $\mu B=\u2009200$, $\mu D=\u20091$, and $CoV=0.1$. The initial population size at $t0=0$ is set to zero.

#### 1. Channel variability

Figure 1 illustrates the effects of channel variability. Specifically, in the three plots of the figure, we condition the trajectories to a fixed value of the kinetic parameters $B(\omega )=\mu B$ and $D(\omega )=\mu D$. In Figure 1(a), we report different trajectories obtained using the MNRM algorithm and generate different realizations of the Poisson processes $Yb$ and $Yd$. The trajectories of $X(t)$ quickly grow from the initial condition $X=0$ and converge to an invariant distribution as $t$ increases. The invariant distribution corresponds to a Poisson distribution with mean population size $\mu B/\mu D$ (see below). For different realizations of the Poisson processes, all trajectories have qualitatively the same behavior and asymptotically fluctuate around this terminal expected value.

In Figure 1(b), not only we sample trajectories conditionally on the parameters’ mean values, but we also condition the dynamics on a particular realization of the Poisson process associated with channel $Rd$. Therefore, the trajectories plotted reflect the variability induced by the birth channel $Rb$ only. It can be observed that, compared to trajectories in Figure 1(a), the set of trajectories in Figure 1(b) is more correlated. A similar behavior is reported in Figure 1(c), where instead of conditioning the system on the mean parameters and a particular realization of $Yd$, it is now conditioned on a realization $Yb(\omega )$ of the birth channel. Note that the conditioning of the channels in the last two plots is based on the realizations of $Yb$ and $Yd$ corresponding to the trajectory reported with a thick magenta line in all plots.

#### 2. Parametric variability

Similarly to the previous exercise, we now condition the trajectories of $X(t)$ on the same realizations used previously to condition on $Yb$ and $Yd$; the MNRM algorithm is then applied for different realizations of the kinetic parameters *B* and *D*. We also report in Figure 2 the previous trajectory corresponding to the expected rate parameters $\mu B$ and $\mu D$, shown again with a thick magenta line.

Figure 2(a), which should be contrasted with Figure 1(a), illustrates the effect of sampling the kinetic parameters only. It is seen that this results in a variability of the trajectories that is significantly higher than the variability obtained by sampling only the Poisson processes. Moreover, about the trajectories in Figure 2(a) highlight the fact that the parametric variability induces trajectories that are significantly more correlated than for the channels variability. The parametric uncertainty is alsoseen to have a greater impact on the population size, with fluctuations of the trajectories around a time-averaged that change from a trajectory to another. The variability among trajectories of the asymptotic time-average contrasts with the previous case of the channels variability, where the trajectories fluctuate around the same value.

To gain further insight into the effects of the uncertain kinetic parameters, we subsequently also condition the trajectories on $D=\mu D$ or $B=\mu D$, see Figures 2(b) and 2(c), respectively. It is seen that the correlation structure between the sampled trajectories is significantly different depending on the sampled parameter, *B* or *D*. In particular, it is seen that considering variability in *B* only (Figure 2(b)) induces trajectories that remain very correlated in time, while for a variability in *D* only (Figure 2(c)) a time shift between the trajectories is also visible. These different effects can be understood from the structure of the solution as expressed in (12) and (13). From the plots, the variability of the trajectories reported in Figure 2(c) appears larger than that in Figure 2(b), suggesting a more important impact of the uncertainty in *D*, than in *B*, on the variability of *X*. However, the experiments reported in this section are not sufficient to confidently draw such a conclusion: the complete variance decomposition must be performed to properly assess the relative importance of each of the variability sources.

### B. Variance decomposition

The variance decomposition is performed using Algorithm 2 and its generalization. We take advantage of the BD model simplicity to accurately estimate the partial variances, using large sample sets with $M=5\u22c5107$.

We start by performing the simple parametric/channel variance decomposition, on the quantity of interest consisting of the system state, $G=X(t)$, at different times *t*. Figure 3 reports the evolution with *t* of the total variance of *X* and its decomposition in parametric ($Vpar$), channels ($Vch$), and mixed ($Vmix$) contributions. The results in Figure 3(a) focus on the early transient stage for $t<\u200910$. We see that the variability in $X(t)$ is clearly dominated by the parametric uncertainty, except at the early stages, $t<\u20090.5$, where the channels variability is larger than parametric variability. The results also indicate that for $t>\u20095$ the total and parametric variances are essentially constant, the latter accounting for roughly 80% of the former.

In contrast, the channels’ partial variance $Vch$ is seen to reach a maximum for $t\u22482.5$ and then progressively decreases while the variance associated with mixed effects monotonically increases. Inspecting Figure 3(b), which reports the evolution of the partial variances over an extended time range, we observe that $Vch$ continues to decrease and eventually becomes negligible at $t\u2248100$, while the variability of the mixed effects levels to $20%$ of the total variance. These evolutions can be understood from Eq. (13): as time advances, the impact of the uncertain parameters on the scaled times becomes more and more pronounced such that $Vmix$ increases and $Vch$ decreases.

To gain additional insight into individual effects of the parameters and the channels, we provide in Figure 4 the complete separation of the partial variances $Vpar$, $Vch$, and $Vmix$. Figure 4(a) shows the decomposition of $Vpar$ into pure contributions from $B$ and $D$, and mixed effects of $B$ and $D$. It is seen that the partial variances $VB$ and $VD$ become very close for $t>\u20096$, while the interaction between the two parameters remains negligible.

Regarding the decomposition of $Vch$ into contributions of the two channels and their interaction, we observe in Figure 4(b) that the two variances $VYb$ and $VYd$ decay at the same rate and remain sensibly equal, while the mixed contribution $VYb,Yd\u22480$ denotes a negligible interaction between the two channels. It should be remarked that this finding is seemingly in contradiction with the finding reported in our previous work,^{16} where the interaction between the two channels was increasing over time with asymptotically $VYb\u2248VYb,Yd$ and $VYd\u22480$. However, no parametric uncertainty was considered in Ref. 16.

Figure 4(c) indicates that mixed effects are still existing in the present system, but they essentially consist in interactions between parameter $B$ and the two channels: $Vmix\u2248VB,Yb+VB,Yd$. Note also that third order variances, e.g., $VB,Yb,Yd$, are negligibly small and are thus not reported.

Figure 5 reports the evolution of the total sensitivity indices $TB,D,Yd,Yd$ scaled by the total variance. The results in Figure 5(a) correspond to $CoV=0$, with all the variability due to the channels and a larger effect of $Yb$ compared to $Yd$. In Figure 5(b), corresponding to $CoV=0.015$, a low contribution of the parameters is visible with a more noticeable impact of *B* compared to *D*, but the variability due to the two channels remains predominant. The total variance is also seen to be slightly higher than for $CoV=0.0$.

Increasing CoV to 0.05, in Figure 5(c), we see that the birth rate parameter *B* becomes the dominant source of variability for $t>2.5$, while parameter *D* and the channels have comparable impact. Note also that the variance is larger than previously observed. Finally, Figure 5(d) corresponding to $CoV=0.15$ shows that in this case most of the variability is due to the kinetic parameters, which contribute roughly equally at $t\u224810$. The total contribution of the channels is small and dominated by interactions with the parameters. We also observe that the total variance is up to 10 times larger for $CoV=0.15$ compared to the case of $CoV=0$.

We have just seen that the different variabilities in $X(t)$ increase with CoV, making difficult the quantitative comparison of their relative impacts. We provide in Figure 6 the dependence on CoV of the first order (Figure 6(a)) and total (Figure 6(b)) sensitivity indices of $X$ at a fixed time $t=10$. The figures show that as CoV increases, the first order sensitivity indices associated with the channels, $SYb$ and $SYd$ decrease in favor of $SB$ and $SD$ which become dominant for $CoV\u22730.05$. A similar trend is reported for the total sensitivity indices, shown in Figure 6(b), though $TYb$ and $TYd$ decrease with a slower rate than $SYb$ and $SYd$, denoting the emergence of stronger interactions between the parametric and channel variabilities and the emergence of non-additive effects.

### C. Comparison with local sensitivity estimates

The proposed method is global, in the sense that it measures the variabilities induced by the considered sources, either parametric or inherent, over the whole events set. It also provides a natural and objective way to compare the respective influences of different variability sources. These characteristics are contrasted with the classical approaches of local sensitivity measures, involving derivatives at nominal points, and which can be limited to low variability levels. In addition, comparing local sensitivity associated with different sources may be difficult, as it usually involves derivatives with respect to different quantities, sometimes hardly comparable (temperature, pre-exponential factor, …).

Our objective in this section is to demonstrate that our approach is consistent and agrees with the local derivative-based parametric sensitivity estimates when the coefficient of variation of the kinetic parameters goes to zero. Specifically, considering again for quantity of interest $G(X)=X$, we introduce the local sensitivity index,^{4} associated with the birth rate parameter as follows:

A similar expression can be written for $LD(b,d)$. In fact, for the birth-death process, the exact expression for the conditional expectation in Eq. (14) is available

Accounting for the initial condition, $x0=0$, we see that the local sensitivity indices at the nominal parameter values are given by

and

Thus, for large enough times, one gets in our case $LB(\mu B,\mu D)\u2248\mu D\u22121=1$ and $LD(\mu B,\mu D)\u2248\mu B/\mu D2=200$.

In order to show how global and local sensitivity indices relate to each other, we consider the case of the birth rate and assume a sufficient regularity of $\mathbb{E}[X|b,d]$ with respect to $b$; we then have the following Taylor series:

From this expansion, one can easily establish the expression of the variance of the conditional expectation,

where $hot$ stands for higher order terms involving higher order moments in $(B\u2212\mu B)$, and that are assumed to vanish faster than $\sigma B2$ when $CoVB\u21920$. Therefore, it follows that

with a similar expression in the case of the uncertain death rate *D*. Note that these expressions hold in the limit of a vanishing variability of all kinetic parameters. A bias term, e.g., $\mathbb{E}[X|\mu B]\u2212\mathbb{E}[X|\mu B,\mu D]$, would be otherwise introduced.

In Figure 7 we check the validity of relation in Eq. (15) in our computations, reporting as a function of the number of MC samples *M* the ratios of the partial variances, $VB$ and $VD$ for *X* at time $t=10$, with the respective variances of the parameters, $\sigma B2$ and $\sigma D2$. Different values of CoV are reported as indicated. The plots show that the MC estimates of $VB/\sigma B2$ and $VD/\sigma D2$ are converging with increasing number *M* of samples. However, the magnitude of the fluctuations are increasing as CoV decreases, denoting an increasing sampling error. This trend can be explained by the large number of samples needed to accurately estimate the partial variances $VB$ and $VD$ when the global variance becomes more and more dominated by the channels’ partial variances.

Conversely, the sampling error is seen to converge quickly with *M* when CoV is larger, but we should not expect the ratios to necessarily converge to the exact squared local sensitivity index, $LB2$ and $LD2$ which are shown with dashed lines in the plots because of the contributions of higher order moments. This is a classical trade-off problem, between bias and sampling errors, that is also present in the numerical estimation of the local sensitivity indices. In this context, one has to select an appropriate value for *h* that leads to a small enough bias in the approximation of the derivative by the difference formula, and a small enough sampling error in the MC estimate of the difference which demands more samples as *h* decreases.

Different methods (see, for instance, Refs. 4, 6, and 7) have been proposed recently to reduce and balance the bias and sampling errors in the estimation of the local sensitivity indices. In particular, correlations between the estimates of the two expectations in the difference formula (14) and multilevel Monte Carlo constructions have been proposed to obtain accurate estimators of the local sensitivity indices at a reduced cost; we plan to adapt these techniques to improve the computational efficiency of the MC estimation of the global sensitivity indices in future works.

Anyway, we emphasize that, in the case of finite CoV, a difference between the ratios and the local sensitivities is not a limitation of our method. It underlines the importance of accounting *globally* for the parametric uncertainty effects and reveals that a local characterization is not appropriate. Further, even if a local sensitivity coefficient (say $LB$) agrees with its counterpart ($VB/\sigma B2$), the presently developed global sensitivity analysis indicates that the local index may fail to account for some of the variability actually induced by the parameter (*B*). Indeed, the parameter may be responsible for an additional variability through its interactions with the channels (and the other parameters). In this case, the *total* sensitivity index is the appropriate measure of induced variability and restricting the analysis to local indices may be misleading.

## V. SCHLÖGL MODEL

In this section, we consider the Schlögl model,^{29} which has $Mr=4$ reaction channels, ${Ri,i=1,\u2026,4}$, corresponding to the reaction mechanism

Species $B1$ and $B2$ are assumed in large excess with constant population size, set, respectively, to $x1=x2/2=\u2009105$. The stochastic state then reduces to the single evolving species *S* with $Ms=1$. The uncertain propensity functions are

The parameters $C1$ to $C4$ in the definition of the propensity functions are considered as independent random variables with uniform distributions. As previously, we define the uniform distributions by their means $\mu 1\u2264j\u2264Mr$ and coefficients of variation CoV; again the latter is taken to be the same for all the parameters. Regarding the expected values of $Ci$, we set $\mu 1=0.015,\mu 2=10\u22124/6,\mu 3=200$, and $\mu 4=3.5$. Finally the initial condition is $X(t=0)=250$.

### A. Qualitative analysis

To gain an understanding of the qualitative impact of the different sources of variability we select a reference trajectory, corresponding to the mean value of the kinetic parameters and a specific realization of the Poisson processes. This particular trajectory is plotted using a thick magenta line in the following plots. We draw samples of trajectories conditioned on the same parameter values or Poisson processes.

Figure 8(a) shows several trajectories of the state conditioned on fixed values of the kinetic parameters $Cj(\omega )=\mu j$ and sampling of the Poisson processes *𝒀*. These trajectories illustrate the random dynamics of the system, that exhibit a bifurcation: starting from the initial condition, the system selects in an initial stage one or the other branch, called lower and upper branches in the following, and later fluctuates around the selected branch. It is seen that the inherent (channels) variability acts on both the branch-selection process and the subsequent fluctuations.

Figure 8(b) depicts trajectories of the system for sampled values of the parameters, but conditioned on the reference realization of the Poisson processes. The sampling of the kinetic parameters uses $CoV=0.1$. The plot shows that although the trajectories are conditioned on the same realization of the channel dynamics, the branch selected changes from one sample of the kinetic parameters to another. This behavior highlights the sensitivity of the trajectories and of the branch selection process with regard to both sources of variability.

In addition, qualitative differences are reported between the two sets of sampled trajectories in Figure 8. First, it is seen that the variability level is much larger when sampling the kinetic parameters, in Figure 8(b), than when resampling the Poisson processes in Figure 8(a). More precisely, the fluctuation amplitudes around the selected branches appear to be essentially the same for the inherent and parametric resampling, but the latter induces an additional variability by displacing the location of the attracting branches.

The analysis of the inherent and parametric variability can be further refined by distinguishing the elementary effects of each channel and parameter. This is done by resampling the Poisson processes, $Yj$, and kinetic parameters, $Cj$, one at a time. Our previous work^{16} focused on the resampling of the Poisson processes $Yj$ to investigate the relative variability induced by the different channels in the absence of parametric uncertainty. The conclusion of this analysis is briefly summarized in Figure 9, which shows that channels $R1$ and $R4$ are the principal sources of inherent variability in the state $X$, while $R2$ and $R3$ are less important. This is explained by the selected branch which depends more heavily on the processes $Y1$ and $Y4$ than on $Y2$ and $Y3$. Indeed, one can observe that for all the resampled processes $Y2$ and $Y3$ the system has selected the same (upper) branch.

Figure 10 shows the individual effect of each kinetic parameter on the reference trajectory. The resampling of $Cj$ uses increasing coefficients of variation: $CoV=\u20090.05$ (top row), 0.10 (center row), and 0.15 (bottom row). Focusing first on the case of the largest parametric variability $CoV=0.15$, we observe that the resampling of the reference trajectory leads to very different effects depending on the considered parameter. Specifically, the resampling of $C1$ and $C2$ induces a large variability of the resulting trajectories. In fact, both the location of the attracting branch and the time scale of the selection process are evidently affected when $C1$ and $C2$ are resampled. On the contrary, resampling $C3$ and $C4$ appears to essentially affect the selected branch but not its location. Also, $C4$ seems to have a greater impact on the time scale to reach the selected branch, compared to $C3$. Finally, we observe that reducing the coefficient of variation, CoV, reduces the variability in the trajectories obtained by resampling $C1$ and $C2$. This is due to the decreasing uncertainty in the location of the attracting branches. A similar effect can be seen on the time scale for reaching an attracting branch and the exchange of attracting branch when resampling $C2$ and $C3$ with a lower CoV.

### B. Global sensitivity analysis

The previous exercises with different resampling of the channel processes and kinetic parameters highlight the complex structure governing the effects and interactions between inherent and parametric variability sources. This complexity stresses the need for quantitative and generic methodologies enabling the quantification and importance ranking of individual sources, as performed below, in order to fully understand the variability in the model output. It also underlines the importance of the global sensitivity approach, as the stochastic model response can substantially depend on the input uncertainty level (CoV).

#### 1. Variance decomposition

As for the birth-death example, we start by performing the parametric/channel variance decomposition, for the system state $G(X)=X(t)$; unless otherwise noted $CoV=0.1$. The MC estimators of the partial variances reported in this section use sample sets with size $M=106$.

Figure 11 depicts the time evolution of decomposition of $\mathbb{V}[X(t)]$, into parametric, channels and mixed contributions, for $t\u2208[0,8]$. The results show that, for this value of CoV, the variability in $X(t)$ is predominantly caused by the uncertainty in the kinetic parameters, which accounts for up to $80%$ of the total variability. The remaining fraction of the variance is essentially due to the mixed effects, particularly for $t>2$ where the channels contribution $Vch$ is negligible. The dominance of $Vpar$ is explained by the parametric impact on the attracting branch locations; changing the kinetic parameters immediately affects the asymptotic time-averages and even possibly the selected branch, causing significant variability. The reason for the low contribution of $Vch$ to $\mathbb{V}[X]$ is less obvious and requires additional analysis, as further discussed below.

To have a better perception of the structure of the total variability in the model state, we present in Figure 12 the time evolution of the first-order partial variances. Starting with the first-order partial variances associated with the kinetic parameters, $VC1$ to $VC4$, shown in Figure 12(a), we observe that the parameter having the largest effect on $X(t)$ is $C1$, accounting for roughly $50%$ of the parametric variance, followed by $C4$ ($\u224825%$ of $Vpar$) and $C2$ ($\u224810%$ of $Vpar$), while the single effect of $C3$ is relatively small. The two dominant parameters are in fact associated with the propensity functions of the two channels that were shown in Ref. 16 to be responsible for most of the state variability. The figure also reports the sum of the first order partial variances $VCi$, which is seen to amount to roughly $85%$ of $Vpar$, indicating the presence of some interactions between the parameter variabilities, independently of the interactions with the channels. In words, the effects of the parameters are not additive.

Figure 12(b) depicts the evolution of the first-order partial variances, $VYi$, associated with the reaction channels. Again, channels $Y1$ and $Y4$ are seen to dominate the two other channels in terms of first-order partial variances. It is remarked that the partial variances $VYi$ are lower than their parametric counterparts, $VCi$, and that the interactions between channels appear to be less pronounced than the interactions between parameters, as $\u2211iVYi$ is close to $Vch$. Regarding the time behavior, the partial variances associated with the channels are observed to initially grow, reaching a maximum for $t\u22481$, then decrease progressively to roughly half of their maximum values and become constant. The initial growth phase corresponds to the stage where the system evolves from the initial condition to the randomly selected attracting branch, with an increasing effect of the inherent stochasticity during the selection stage. Shortly after the branch has been selected, not only the interactions between the channels develop, as in the case of certain parameters,^{16} but the interactions with the parameters become more and more pronounced, with a resulting significant reduction of the first-order partial variances, $VYi$.

Interactions between the channels and the parameters can be assessed from the plots of Figure 13. In Figure 13(a) we report the total effects associated with each parameter $Ci$ and its interactions with the channels, that is the total contribution of $Ci$ to $Vmix$. It is seen that once more the parameters $C1$ and $C4$ have the highest contributions to $Vmix$ while $C2$ and $C3$ have contributions roughly one half smaller. Further, the sum of these total contributions largely exceeds $Vmix$ indicating the importance of not only interactions between the $Ci$ and the channels stochasticity, but also of the joint interactions between the set of parameters and channels. In fact, the partial variances due to the interaction between single parameter $Ci$, only, and all Poisson processes $\bm{Y}$ are all very low compared to $Vmix$ (not shown). Similarly, one can compute the total contribution of individual channels $Yi$ to $Vmix$. Because $Vch$ is small we directly report in Figure 13(b) the time evolution of the total contributions of $Yi$ to the total variance. Compared to the corresponding first-order partial variances $VYi$, shownin Figure 12(b), we see that the total channel contributions are significantly higher and exhibit a different behavior in time, monotonically increasing to asymptotically constant values. Again, channels $R1$ and $R4$ are found to contribute roughly equally and twice as much as channels $R2$ and $R3$.

#### 2. Total sensitivity indices

Provided in Figure 14 are the total sensitivity indices of $X(t)$, associated with the kinetic parameters $Ci$ and channels’ Poisson processes $Yi$. The total sensitivity indices are scaled by the total variance $\mathbb{V}[X(t)]$, and are reported for different coefficients of variation, $CoV=0$, 0.05, 0.1, and 0.15.

Figure 14(a), where $CoV=0.0$, corresponds to the case of a model having no parametric uncertainty as studied in Ref. 16. The total sensitivity indices $TCi$ of the parameters are then 0, and the variability is principally due to channels $R1$ and $R4$ roughly equally, and to a lesser extent to channels $R2$ and $R3$. Increasing the coefficient of variation of the kinetic parameters, we observe successively in Figures 14(b)–14(d) that the total sensitivity indices of the parameters, $TCi$, become progressively dominant compared to the total sensitivity indices of the channels, $TYi$. It is remarked that the total sensitivity indices of the parameters do not increase proportionally with the squared coefficient of variation, as one would have expected for a linear parametric dependence. Also, the total variances due to the channels are seen to slightly diminish with increasing CoV. Regarding the parametric total indices $TCi$, the plots indicate that $TC1>TC4>TC2>TC3$ for the whole range of CoV values presented and that $TC1$ becomes more dominant as CoV increases.

### C. Local sensitivity analysis

We finally contrast the outputs of the global sensitivity analysis above with those of a local sensitivity characterization. It should be clear that a local characterization of the Schlögl model is prone to be inappropriate because of the complicated non-linear behavior of the system, which involves third-order reactions, and the non-trivial dependencies on the parameters. Although trajectories of the model have complex dependencies with respect to the uncertain kinetic parameters, the statistical moments of the quantity of interest may behave well and smoothly as the parameters change; conducting derivative-based local approaches is then quite legitimate.

However, we show below that using local information does not necessarily provide a useful characterization of the variability in the quantity of interest. In particular, following the local approach previously introduced in Sec. IV C for the birth-death model, we wish to estimate by $V\u223cCi$ the first-order variances $VCi$ associated with the uncertain kinetic parameters, using (see Eq. (15))

where $LCi$ is the derivative, with respect to the *i*th kinetic parameter, of the averaged quantity of interest at the parameter mean values $Ci=\mu i$. In our case, $LCi$ is then the derivative of the mean population size, since $G(X)=X$. It should be clear that relying on moments to characterize the influence of parameters can be of limited relevance because of the branching nature of the Schlögl model. Specifically, the average trajectory and its sensitivity to parameters is not necessarily representative of the system response, especially after the system has selected one or the other branch.

To support these claims, we provide in Figure 15 a comparison of the first-order and total-order variances associated with the four uncertain kinetic parameters $Ci$ and the local estimate $V\u223cCi$ as defined above. In practice, we estimate the derivatives $LCi$ of $\mathbb{E}[G(X)]$ with a large Monte Carlo sample set (typically $M=106$ samples) using different values for *h* in Eq. (14) to optimize the trade-off between bias and sampling error. The errors in $V\u223cCi$ reported below were estimated to be roughly 5%. Figure 15(a) corresponds to the case of $CoV=0.1$ at an early time $t=0.1$, where the system has not yet clearly selected the random branch. For such a short time, the distribution of $G(X)$ for the nominal kinetic parameters is still uni-modal and the local approach is meaningful. Indeed, the estimates $V\u223cCi$ are not far from the first-order partial variances $VCi$, as one would expect. However, it is already observed that the variability due to the interactions of $Ci$ with other parameters and inherent stochasticity is not accounted by $V\u223cCi$, although these interactions amount to roughly half of the induced total variability. Figure 15(b) provides the same comparison but at a later time, $t=4$. At this time, the system has, with high probability, already selected its attracting branch: the distribution of $G(X)$ for the nominal parameters is bimodal, and the derivative of its expected value is not representative of the effects of the parameters. Specifically, we see that the estimates $V\u223cCi$ completely overestimate not only the first-order partial variances $VCi$ but also the total variances induced by parameter $Ci$, namely, $\mathbb{V}[X]TCi$. It could be argued that, though the local estimates of the variances are far off the actual partial variances associated with the parameters, the local approach is still capable of ranking the relative importance of the kinetic parameters. This is, however, not true in general, as can be seen in Figure 15(c) which reports the previous variances at $t=4$, in the case of a coefficient of variation $CoV=0.5$: the global sensitivity analysis indicates that the most important source of variability is $C2$ (because of non-linear saturation effects), whereas with the local approach the ranking is obviously the same as for $CoV=0.1$ because in the local approach it is independent of CoV.

## VI. CONCLUSIONS AND FUTURE WORK

In this work, we have extended the variance decomposition approach proposed in Ref. 16 to perform sensitivity analysis in stochastic simulators having uncertain parameters. The parametric uncertainty is treated as an additional source of variability, assumed to be independent of the inherent stochasticity of the reaction channels. This assumption permits the decomposition of the variance of second-order quantities of interest derived from the system state, distinguishing between different variability sources, in particular the channel randomness and the parametric contributions. The approach differs from the classical parametric sensitivity analyses based on the parametric dependences of some statistical moments of the quantities of interest. It offers practical means to fully characterize the complex interplays between inherent sources of noise and individual uncertain parameter effects, through the evaluation of global sensitivity indices associated with an arbitrary set of parameters and random channels.

A Monte Carlo algorithm has been proposed and implemented to illustrate the computation of these sensitivity indices. It has been applied to simple stochastic systems, namely, the birth-death and Schlögl models. Through these examples, we have shown that inherent and parametric variability sources influence the model output in a different fashion. In particular, the impact of uncertain parameters becomes more pronounced when their coefficients of variation are increased. In contrast, the inherent variability due to the reaction channels is less affected by the amount of parametric uncertainty, demonstrating the irreducibility of the channels stochasticity.

Numerical examples have also highlighted the importance of a global sensitivity approach, as opposed to local sensitivity methods using derivatives with respect to the uncertain parameters of statistical moments. In particular, our computations showed that first-order derivatives of the first moment was not able to properly account for the interactions between parameters and is not suitable to assess the global variability induced by a given parameter, even for moderate uncertainty levels. To some extent, these limitations can be remedied considering high-order and mixed derivatives, a global-derivative-based approach^{30} or more advanced representations of the parametric dependencies of the moments.^{31–33} We advocate for the simplicity and generality of the global variance decomposition method, as other sensitivity approaches may miss the characterization of the variance incurred due to the stochasticity of the channels, which is essential for a complete understanding of the stochastic dynamics. The advantages afforded by global sensitivity approach are thus regarded as highly beneficial in a range of applications. Specifically, besides providing an effective means for avoiding the limitations of local sensitivity analyses and moment based approaches, which were highlighted above based on simple model simulations, the present approach also offers the capability of quantifying the impact of parametric uncertainties, the inherent stochasticity of individual channels and of groups of channels, as well as mixed effects involving reducible and irreducible sources of uncertainty. Examples where these features would be of special importance include the analysis of system robustness, model reduction, and parameter calibration.^{4,34–38}

The application of the proposed method to complex stochastic systems remains challenging because of the sampling effort. Even though the method is trivial to implement in parallel, future works should focus on the improvement of the methods to lower its computational complexity. Specifically, the results presented in this work relied on exact simulators, such that the errors in the MC estimation of the sensitivity indices are directly related to the sampling error. As for any standard MC method, a rate in $1/M$ is expected for the convergence of the sensitivity indices, and a large enough sample set must be used to properly explore both the parametric and Poisson process domains and reduce the variance of the estimators. Note that a different number of samples *M* could be used depending on the considered sensitivity indices, as their estimators may have fairly different variances. In any case, the convergence with the number of samples is slow and we plan to reduce the computational complexity of Algorithm 2 in the future, by relying on multi-level Monte Carlo strategies.^{39–41} Specifically, we envision to alleviate the computational cost of estimating the correlations, between quantities of interest at different parameter values and Poisson processes, using multi-level estimators. The main difficulty in designing such a multilevel Monte Carlo approach for the estimation of the sensitivity indices will be the enforcement of consistent realizations for the channels’ Poisson processes in tau-leaping approximations.^{23,42–44} Indeed, one would need to sample the same realization of a Poisson process for different time step *and* parameter values and not just to generate highly correlated but distinct realizations of the Poisson process. This can be achieved via a conditional sampling of the Poisson process increments. This would necessitate devising a multi-level procedure capable of balancing the sampling and the bias errors of a set of estimators for the sensitivity coefficients. Progress along this direction will be reported elsewhere.

## NOMENCLATURE

- $\Omega ,\omega \u2208\Omega $
Event space, random event

- $\mu ,\Sigma $
Probability measure, $\sigma \u2212$algebra

- $Ms$
Numbers of chemical species

- $S1,\u2026,SMs$
Chemical species

- $Xi$
Random number of $Si$ molecules

- $\bm{X}=(X1\u2026XMs)$
Random state vector

- $Mr$
Number of reaction channels

- $R1,\u2026,RMr$
Reaction channels

- $Yj$
Poisson process of channel $Rj$

- $\bm{Y}=(Y1\u2026YMr)$
Vector of Poisson processes

- $aj$
Propensity functions of channel $Rj$

- $\bm{\nu}j$
Stoichiometric vector of channel $Rj$

- $t,tj$
Time and scaled time of channel $Rj$

- $RGj$
Random number generator for $Yj$

- $Nc$
Number of kinetic parameters

- $\bm{c},\bm{C}$
Deterministic and uncertain vectors of kinetic parameters

- $p\bm{C}(\bm{c})$
Probability density function of

*𝑪* - $NU=Mr+Nc$
Total number of stochastic sources

- $\bm{U}=(\bm{C},\bm{Y})$
Vector of all stochasticity sources

- $G(\bm{X})$
Second order random functional

$\mathbb{E}[\u22c5],\mathbb{V}[\u22c5]$ Expectation and variance operators

- $Vpar$
Variance of

*G*due to parameters (*𝑪*) - $Vch$
Variance of

*G*due to the channels (*𝒀*) - $Vmix$
Variance of

*G*due to joint effects - $Spar,ch,mix$
First order sensitivity indices

- $Tpar,ch,mix$
Total order sensitivity indices

- $u\u2208{1,\u2026,NU}$
Subset of stochastic sources

- $Su,Tu$
First and total order sensitivity indices associated with u

- $VSu,VTu$
Scaled version of the sensitivity indices

- $LCi$
Local sensitivity index

- V•^,T•^,S•^,…
Monte Carlo estimates

- $V\u223cCi$
Estimate of $VTCi$ based on $LCi$

- CoV
Coefficient of variation

- $B,D$
Uncertain birth and death rates

## ACKNOWLEDGMENTS

This work was supported in part by the SRI Center for Uncertainty Quantification in Computational Science and Engineering at King Abdullah University of Science and Technology, and by the US Department of Energy (DOE), Office of Science, Office of Advanced Scientific Computing Research, under Award No. DE-SC0008789.