# Exploiting deterministic algorithms to perform global sensitivity analysis of continuous-time Markov chain compartmental models with application to epidemiology

Henri Mermoz Kouye, Gildas Mazo, Clémentine Prieur, Elisabeta Vergu

► To cite this version:

Henri Mermoz Kouye, Gildas Mazo, Clémentine Prieur, Elisabeta Vergu. Exploiting deterministic algorithms to perform global sensitivity analysis of continuous-time Markov chain compartmental models with application to epidemiology. 2024. hal-03565729v5

HAL Id: hal-03565729

https://hal.science/hal-03565729v5

Preprint submitted on 18 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Performing global sensitivity analysis on simulations of a continuous-time Markov chain model motivated by epidemiology

Henri Mermoz Kouye, Gildas Mazo, Clémentine Prieur, Elisabeta Vergu

► To cite this version:

Henri Mermoz Kouye, Gildas Mazo, Clémentine Prieur, Elisabeta Vergu. Performing global sensitivity analysis on simulations of a continuous-time Markov chain model motivated by epidemiology. Computational & Applied Mathematics, 2024, 43 (409), 10.1007/s40314-024-02897-y. http://hal-03565729v6

HAL Id: hal-03565729

https://hal.science/hal-03565729v6

Submitted on 22 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Exploiting deterministic algorithms to perform global sensitivity analysis of continuous-time Markov chain compartmental models with application to epidemiology

Henri Mermoz KOUYE(1), Gildas MAZO(1), Clémentine PRIEUR(2) and Elisabeta VERGU(1,†)

(1) Univ. Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France

(2) Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP, LJK, 38000, Grenoble, France

(†) E. Vergu died on May 13, 2023, during the writing of this paper

March 18, 2024

## Abstract

In this paper, we propose a generic approach to perform global sensitivity analysis (GSA) for compartmental models based on continuous-time Markov chains (CTMC). This approach enables a complete GSA for epidemic models, in which not only the effects of uncertain parameters such as epidemic parameters (transmission rate, mean sojourn duration in compartments) are quantified, but also those of intrinsic randomness and interactions between the two. The main step in our approach is to build a deterministic representation of the underlying continuous-time Markov chain by controlling the latent variables modeling intrinsic randomness. Then, model output can be written as a deterministic function of both uncertain parameters and controlled latent variables, so that it becomes possible to compute standard variance-based sensitivity indices, e.g. the so-called Sobol' indices. However, different simulation algorithms lead to different representations. We exhibit in this work three different representations for CTMC stochastic compartmental models and discuss the results obtained by implementing and comparing GSAs based on each of these representations on a SARS-CoV-2 epidemic model.

**Keywords:** stochastic compartmental models, continuous-time Markov chains, epidemic models, global sensitivity analysis, uncertainty quantification.

## 1 Introduction

In epidemiology, stochastic compartmental models facilitate the prediction and understanding of spreads of infectious diseases in a host population, such as humans,
Performing global sensitivity analysis on
simulations of a continuous-time Markov chain
model motivated by epidemiology

Henri Mermoz KOUYE$^{(1,2)}$, Gildas MAZO$^{(1)}$, Clémentine PRIEUR$^{(2)}$
and Elisabeta VERGU$^{(1,†)}$

(1) Univ. Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France

(2) Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP, LJK, 38000, Grenoble, France

(†) E. Vergu died on May 13, 2023, during the writing of this paper

July 22, 2024

Abstract

In this paper we apply a methodology introduced in Navarro Jimenez et al.
[2016] in the framework of chemical reaction networks to perform a global
sensitivity analysis on simulations of a continuous-time Markov chain model
motivated by epidemiology. Our goal is to quantify not only the effects of
uncertain parameters such as epidemic parameters (transmission rate, mean
sojourn duration in compartments), but also those of intrinsic randomness
and interactions between epidemic parameters and intrinsic randomness. For
that purpose, following what was proposed in Navarro Jimenez et al. [2016],
we leverage three exact simulation algorithms for continuous-time Markov
chains from the state of the art which we combine with common tools from
variance-based sensitivity analysis as introduced in Sobol' [1993]. Also, we
discuss the impact of the choice of the simulation algorithm used for the sim-
ulations on the results of sensitivity analysis. Such a discussion is new, at
least to our knowledge. In a numerical section, we implement and compare
three sensitivity analyses based on simulations obtained from different exact
simulation algorithms of a SARS-CoV-2 epidemic model.

**Keywords:** stochastic compartmental models, continuous-time Markov chains, epidemic models, global sensitivity analysis, uncertainty quantification.

# 1 Introduction

In epidemiology, stochastic compartmental models help in the prediction and un-
derstanding of spreads of infectious diseases in a host population, such as hu-
mans, animals, or plants. The output of those models usually depends on nu-
merous uncertain parameters such as transmission rate, mean infectious period
animals, or plants. The outputs of those models depend on numerous uncertain
parameters such as transmission rate, mean sojourn duration in each compart-
ment or transition probabilities. The aim of sensitivity analysis is to identify,
among these parameters, the ones which have the greater impact on the infec-
tion spread [Hanthanan Arachchilage et al., 2023, Goel et al., 2023, Massard et al.,
2022]. This is useful, e.g., to elaborate efficient control strategies [Ngonghala et al.,
2015, Yang et al., 2016] or perform models comparison [Torii et al., 2023]. In the
following, we are interested in the global sensitivity analysis (GSA) of stochastic
compartmental models used in epidemiology [Courcoul et al., 2011]. We focus on
on GSA rather than local sensitivity analysis, as the former is better adapted to
nonlinear models.

Compartmental models [Brauer, 2008] consist in dividing the host popula-
tion into compartments, each containing individuals with a similar health status.
Health statuses of individuals change over time. Transitions between compart-
ments strongly depend on the individuals’ characteristics or the patterns of con-
tacts between them. While in large populations randomness due to individual-to-
individual variability averages out, it has a large impact on the transmission pro-
cess for small populations [Britton, 2009, Bittihn and Golestanian, 2020]. Because
they incorporate stochastic effects related to biological or contact events, stochastic
compartmental models are used to analyze thoroughly the outbreak of infectious
diseases. Throughout, we thus focus on stochastic models, and, more precisely, on
continuous-time Markov chains (CTMC). CTMCs assume that inter-event times
are distributed according to exponential distributions.

GSA aims at determining the extent to which the variability of an input parame-
ter or of a set of input parameters affects the variability of the model’s output (see,
e.g., Saltelli et al. [2000], Marino et al. [2008]). Performing GSA for stochastic
models is more complex as the model output is tainted with two sources of un-
certainty: the intrinsic randomness of the model and the uncertainty on the model
parameters (such as mean sojourn duration in each compartment, transmission rate
and others). So far, different paradigms have been introduced in the literature for
sensitivity analysis of stochastic models.

To do GSAs of stochastic models, a pragmatic approach consists in performing
the analysis on both the expectation and the variance of the model output condi-
tionally on the uncertain parameters. Intrinsic randomness is thereby considered as
noise and is smoothed by the average. The stochastic model is thus reduced to two
deterministic models, for which it is possible to estimate sensitivity indices of the
uncertain parameters, such as Sobol' indices [Sobol', 1993]. This approach is of-
ten used in practice in various applications, for instance in: Courcoul et al. [2011]
to identify key parameters of a model describing the spread of an animal disease
in a cattle herd; Rimbaud et al. [2018] for a model describing the spatio-temporal
spread of plant pathogens; Richard et al. [2021] for a SARS-CoV-2 spread model;
Cristancho Fajardo et al. [2021] for a theoretical metapopulation model. However,
such an approach can suffer from inconsistent conclusions. Since GSA is per-
formed separately on the conditional mean and conditional variance, a parameter
or mean sojourn duration in the different compartments. The aim of sensitivity analysis is to identify, among these parameters, the ones which have the greater impact on the infection spread [Hanthanan Arachchilage et al., 2023, Goel et al., 2023, Massard et al., 2022]. This is useful, e.g., for elaborating efficient control strategies [Ngonghala et al., 2015, Yang et al., 2016] or performing model comparisons [Torii et al., 2023]. In the following, we focus on global sensitivity analysis (GSA) of stochastic compartmental models commonly used in epidemiology [Courcoul et al., 2011]. We focus on GSA rather than local sensitivity analysis, as the former is better adapted to nonlinear models.

Compartmental models [Brauer, 2008] consist in dividing the host population into compartments, each containing individuals with a similar health status. Health status of each individual changes over time. Transitions between compartments are highly dependent on individual characteristics or the contact pattern between individuals. While in large populations randomness due to individual-to-individual variability averages out, it has a large impact on the transmission process for small populations [Britton, 2009, Bittihn and Golestanian, 2020]. As they incorporate stochastic effects related to biological or contact events, stochastic compartmental models are used to analyze thoroughly the outbreak of infectious diseases. Throughout this paper, we thus focus on stochastic models, and more precisely on continuous-time Markov chains (CTMC). A CTMC is a continuous time and memoryless discrete event stochastic process, which means that the past history impacts on the future evolution of the system only via the current state of the system. For a CTMC, inter-event times are exponentially distributed.

GSA aims at determining the extent to which the variability of an input parameter or of a set of input parameters affects the variability of model output (see, e.g., Saltelli et al. [2000], Marino et al. [2008]). In the framework of a deterministic model, a variance-based global sensitivity analysis is often used to rank the importance of input parameters (or sets of input parameters), based on their contribution to the variance of the output quantity of interest (QoI), via the computation of the so-called Sobol' indices introduced in Sobol' [1993]. Performing GSA for stochastic models is more complex as the model output is tainted with two sources of uncertainty: the intrinsic randomness of the model and the uncertainty on epidemic model parameters (such as mean sojourn duration in the different compartments, transmission rate and others). So far, different paradigms have been introduced in the literature for sensitivity analysis of stochastic models.

A pragmatic approach for GSA of stochastic models consists in performing the analysis on both expectation and variance of model output conditionally on epidemic uncertain parameters. Intrinsic randomness is thereby considered as noise and is smoothed by averaging. More generally, this approach amounts to summarizing the output of the stochastic model by several deterministic QoIs (obtained by smoothing intrinsic random noise) and then to perform a GSA for each of these QoIs, e.g., by computing Sobol' indices. This approach is often used in practice in various applications, for instance in: Courcoul et al. [2011] to identify key parameters of a model describing the spread of an animal disease in a cattle herd;
can appear to be important for a quantity and not for the other one. Moreover, it may be difficult to efficiently estimate first order Sobol' indices for the conditional mean, as this requires a fine trade-off between explorations and repetitions [Mazo, 2021]. Often, to avoid computational burdens, GSA is performed on meta-models. In Marrel et al. [2012], the conditional mean and variance are emulated by a joint Gaussian process-based model. In Étoré et al. [2020], the conditional mean or the mean exit time from any compact domain are emulated by polynomial chaos. A second approach, which avoids the loss of information induced by averaging over intrinsic randomness, has been introduced by Hart et al. [2017]. In this latter paper, the authors interpret the Sobol'-Hoeffding decomposition (see, e.g., Hoeffding [1948], Sobol' [1993]) as a random decomposition, where the randomness is inherited from the intrinsic randomness of the model. From this decomposition they obtain Sobol' indices (variance-based sensitivity indices) which are themselves random variables. More recently, a different point of view was adopted (see, e.g., Fort et al. [2021], Veiga [2021]). Stochastic models as interpreted as deterministic models with outputs taking values in some set of probability distributions. Sensitivity indices that measure the sensitivity of the output probability distribution to variations of the input parameter values are defined. Finally, another path to handle sensitivity analysis of stochastic models consists in emulating the output probability distribution (see, e.g., Zhu and Sudret [2021, 2023]).

However, none of the aforementioned approaches enables the possibility of doing a complete GSA, in which not only the effects of the uncertain parameters are quantified, but also those of intrinsic randomness and the interactions between the two. The problem of performing a complete GSA for stochastic models is tackled, e.g., in Le Maître and Knio [2015], Jimenez et al. [2017], Navarro Jimenez et al. [2016]. In these papers, model outputs are simulated from a stochastic algorithm. Then the approach consists in reinterpreting this stochastic algorithm as a deterministic one with an augmented set of inputs comprising both the uncertain parameters and the latent variables controlling intrinsic randomness, thus allowing for the application of standard GSA tools to the deterministic algorithm defined in the augmented input space.

There exist different algorithms to simulate a CTMC. To perform complete GSAs of CTMC compartmental models, therefore, we investigate three exact simulation algorithms. The first, Gillespie Direct Method algorithm [Gillespie, 1976], is the most common in epidemiology. However this algorithm does not allow us to separate and control the different types of transitions in the model. Thus we consider two other representations based respectively on Gillespie First Reaction Method [Gillespie, 1976] and Modified Next Reaction Method [Anderson, 2007]. The latter was used for doing a complete GSA of reaction networks in Navarro Jimenez et al. [2016]. The underlying idea of Modified Next Reaction Method is a random time change introduced in Ethier and Kurtz [1986].

Each of the above simulation algorithms leads to a distinct deterministic representation of the same stochastic model. But there is no reason for the results of the corresponding sensitivity analyses to coincide. In the present paper, we aim
Rimbaud et al. [2018] for a model describing the spatio-temporal spread of plant pathogens; Richard et al. [2021] for a SARS-CoV-2 spread model; Cristancho Fajardo et al. [2021] for a theoretical metapopulation model. It is important to note that if we evaluate a stochastic model multiple times for the same values of uncertain input parameters, we obtain different values for the model output, due to the intrinsic randomness of the model. In this framework, putting in practice the aforementioned approach relies on a fine trade-off between exploration of parameter space and repetition [Mazo, 2021]. To avoid computational burdens due to the cost of model evaluations, it is possible to perform GSA on a metamodel. A metamodel/ surrogate model is a simplified model approximating the actual model at a much lower cost. In Marrel et al. [2012], the mean and variance are jointly emulated by Gaussian process regression (see, e.g., Williams and Rasmussen [2006]). In Étoré et al. [2020], the sensitivity analysis of certain QoIs calculated from a stochastic differential equation, for example a hitting time, is based on metamodeling by polynomial chaos expansion. The aforementioned approach, based on deterministic QoIs to summarize the results of a stochastic model, suffers from two main drawbacks. The first is that certain parameters may prove influential on a QoI but not on others, because the GSA is carried out separately on the different QOIs. The second is that this approach provides no information on how intrinsic noise influences the model output.

More recently, a different point of view was adopted (see, e.g., Fort et al. [2021], Veiga [2021]). Stochastic models are interpreted as deterministic models whose output for a given set of model parameter values is a probability distribution. Then it is possible to emulate the output probability distribution (see, e.g., Zhu and Sudret [2021, 2023]) and/or to define sensitivity indices that measure the sensitivity of the output probability distribution to variations of input parameters.

In the framework we consider in the present paper, it is possible to control intrinsic randomness (e.g., by fixing the seed in a code). It is the framework considered in Hart et al. [2017] where Sobol' indices are computed for each realization of the variable controlling the internal noise. It results in random Sobol' indices, whose randomness is inherited from the intrinsic randomness of the model. Still in this framework, a different approach is adopted in Le Maître and Knio [2015], Jimenez et al. [2017], Navarro Jimenez et al. [2016]. It is the one we focus on in our work. The stochastic algorithm used to simulate the model output is reinterpreted as a deterministic one with an augmented set of inputs comprising both the uncertain parameters and the latent variables controlling intrinsic randomness. This allows to apply standard GSA tools to the algorithm defined on the augmented input space, and, thereby, quantify not only the effects of the uncertain parameters (the epidemic parameters in our setting), but also those of intrinsic randomness and the interaction effects between the uncertain parameters and intrinsic randomness. We shall call this the complete GSA method.

In Le Maître and Knio [2015], Jimenez et al. [2017], Navarro Jimenez et al. [2016], the Modified Next Reaction Method algorithm, which is well-known in the field of simulation of chemical reaction networks, was used to simulate the Markov
at exploring this issue, which to the best of our knowledge has not been tackled
in the literature yet. Once this issue has been discussed in the general framework
of sensitivity analysis of stochastic models, we investigate which representations
can be introduced for CTMC compartmental models. We discuss the choice of
representations, which depends on the type of results we are interested in.

The paper is organized as follows. We first introduce in Section 2.1 the class
of compartmental models we are interested in, whose mathematical formulation
is a CTMC. Then in Section 2.2 we recall the definition of variance-based Sobol'
indices [Sobol', 1993] for sensitivity analysis of deterministic models. Then we
explain in Section 3 how it is possible to perform a complete GSA of a stochastic
model if we can represent it as a deterministic function of controlled latent vari-
ables (modeling the intrinsic randomness) and uncertain parameters. We also ex-
hibit in Section 3.2 a toy example showing that such a representation is not unique
and that GSA results depend on the choice of representations. Thus GSA results
have to be interpreted with caution. In Section 4, we propose three deterministic
representations for CTMC stochastic compartmental models based on three differ-
ent simulation algorithms. We illustrate our approach in Section 5 by considering a
parsimonious SARS-CoV-2 spread model as a case study. We compare and discuss
the GSA results obtained with the different representations presented in Section 4.
Finally the main conclusions of our study are recalled in Section 6.

## 2 Preliminaries

We first recall in Section 2.1 the definition of CTMC stochastic compartmen-
tal models we are interested in. Then in Section 2.2 we recall the definition of
variance-based Sobol' indices [Sobol', 1993] used for GSA in the framework of
deterministic models.

2.1 CTMC stochastic compartmental models

Consider a finite, closed (i.e. of constant size over time) population in which each
individual has a health status (susceptible, infectious, and so on) evolving over
time. The set of all possible health statuses is denoted by **V**. Since those health
 statuses induce a partition of the whole population at any given time, the elements
of **V** are also called compartments. Every time an individual changes compart-
ments, we say that a transition occurs. Only certain types of transitions can occur.
Let **E** denote the set of all possible types of transitions. By definition, for α, β ∈ **V**,
an individual can move from α to β if (α, β) ∈ **E**. The pair (V, **E**) can be iden-
tified with a directed graph where **V** is the set of nodes and **E** is the set of arrows
connecting the compartments between which the individuals can move. Assuming
an ordering of the compartments has been chosen, let us identify the set **E** with
a subset of ℝ^|V| as follows: with each (α, β) ∈ **E**, associate the vector of length
|**V**| with components equal to zero except at the positions corresponding to α and
chain. However, in epidemiology, it is customary to use Gillespie algorithms. One
may then ask whether it is possible to apply the complete GSA method based
on these algorithms, and wether this would lead to similar results for sensitivity
analysis as intrinsic randomness is modeled differently from one algorithm to the
other. To address this problem, we first review three exact stochastic simulation
algorithms—namely, Gillespie Direct Method [Gillespie, 1976], Gillespie First Re-
action Method [Gillespie, 1976] and Modified Next Reaction Method [Anderson,
2007]—and then apply the complete GSA method on each of these algorithms.
Our contribution is to show, both mathematically and numerically, that the results
of the variance-based sensitivity analysis depend on the chosen algorithm. To the
best of our knowledge, this point has never been discussed in the literature so far.

The paper is organized as follows. We provide in Section 2.1 a description
of the class of compartmental models we are interested in, whose mathematical
formulation is a CTMC. In Section 2.2 we recall the definition of variance-based
Sobol' indices [Sobol', 1993] for sensitivity analysis of deterministic models. In
Section 3 we review the methodology of global sensitivity analysis of a stochastic
model consisting in representing this model as a deterministic model with input
space the set of epidemic parameters augmented with the set of latent variables
modeling intrinsic randomness. In Section 3.2, we exhibit a toy example showing
that such a representation is not unique and that GSA results depend on it. Thus
GSA results have to be interpreted with caution. In Section 4, we review the most
common exact simulation algorithms for CTMC stochastic compartmental models
and for each of them we provide the corresponding deterministic algorithm used for
a complete GSA method. We illustrate our approach in Section 5 by considering a
parsimonious SARS-CoV-2 spread model as a case study. We compare and discuss
the GSA results obtained with the different representations presented in Section 4.
Finally the main conclusions of our study are recalled in Section 6.

## 2 Preliminaries

We first recall in Section 2.1 the definition of CTMC stochastic compartmen-
tal models we are interested in. Then in Section 2.2 we recall the definition of
variance-based Sobol' indices [Sobol', 1993] used for GSA in the framework of
deterministic models.

2.1 CTMC stochastic compartmental models

Consider a finite, closed (i.e. of constant size over time) population in which each
individual has a health status (susceptible, infectious, and so on) evolving over
time. The set of all possible health statuses is denoted by $V$. Since those health
 statuses induce a partition of the whole population at any given time, the elements
of $V$ are also called compartments. Every time an individual changes compart-
ments, we say that a transition occurs. Only certain types of transitions can occur.
$\beta$, where the components are -1 and 1, respectively. The elements of $\mathbf{E}$ seen as a subset of $\mathbb{R}^{|\mathbf{V}|}$ are called transition vectors.

For every $\alpha \in \mathbf{V}$ and a vector of epidemic parameters $\theta \in \Theta \subset \mathbb{R}^d$, let $W_{\alpha}^{\theta}(t)$ be the number of individuals in compartment $\alpha$ at time $t$. Since the population is closed, we have that $\sum_{\alpha \in \mathbf{V}} W_{\alpha}^{\theta}(t)$ is constant over time. Denote by $W^{\theta} = (W_{\alpha}^{\theta})_{\alpha \in \mathbf{V}}$ the stochastic process that describes the whole population over time. It is assumed that $W^{\theta}$ is a continuous-time Markov chain with state space $\mathcal{E} \subset \mathbb{N}^{|\mathbf{V}|}$ and positive rate functions $g_{\mathbf{u}}(\theta, \xi)$, $\mathbf{u} \in \mathbf{E}$, $\xi \in \mathcal{E}$, given by

$$ \mathrm{Pr}(W^{\theta}(t+s) = \xi + \mathbf{u} | W^{\theta}(t) = \xi) = g_{\mathbf{u}}(\theta, \xi)s + o(s) \text{ as } s \to 0. $$

The initial state $W^{\theta}(0) =: \xi_0 \in \mathcal{E}$ is supposed to be fixed. A description of the CTMC of the classical SIR model is given in Example 1.

**Example 1.** The classical SIR model is described as follows. There are three compartments $\mathbf{V} = \{S, I, R\}$ and two types of transitions: infection $(S, I)$ and removal $(I, R)$ so that $\mathbf{E} = \{(S, I), (I, R)\} = \{(-1, 1, 0), (0, -1, 1)\}$. Infection is characterized by the transition vector $\mathbf{u}_{S,I} = (-1, +1, 0)$ and the rate function $g_{S,I} = \frac{\beta}{N}W_I W_S$, where $\beta$ is some parameter and $N$ is the total size of the population. Removal has transition vector $\mathbf{u}_{I,R} = (0, -1, +1)$ and rate function $g_{I,R} = \gamma W_I$, where $\gamma$ is some parameter. The vector of parameters is $\theta = (\beta, \gamma)$. The graph of the SIR model is given below:

<figure><img src="image_2.png" /></figure>

## 2.2 Global Sensitivity Analysis

In this section, we first recall the definition of variance-based Sobol' sensitivity indices introduced in Sobol' [1993] for deterministic models with scalar outputs. Following the paradigm of GSA, we model uncertain inputs by a random vector of independent components $\mathbf{X} = (X_1, \dots, X_m)$. Let $E_1, \dots, E_m$ be subsets of $\mathbb{R}$ and $f_1 : E_1 \times \dots \times E_m \to \mathbb{R}$ be some function such that $\mathbb{E} (f_1(\mathbf{X})^2) < +\infty$. Then first-order and total Sobol' indices (see, e.g., Sobol' [1993], Homma and Saltelli [1996]) of the output $Y = f_1(\mathbf{X})$ associated with input $X_j, j = 1, \dots, m$, are respectively defined as:

$$ S_{X_j} = \frac{\mathrm{Var}(\mathbb{E}[Y | X_j])}{\mathrm{Var}(Y)}, \qquad (1) $$

$$ S_{X_j}^{\mathrm{tot}} = 1 - \frac{\mathrm{Var}(\mathbb{E}[Y | X_1, \dots, X_{j-1}, X_{j+1}, \dots, X_m])}{\mathrm{Var}(Y)} \\ =: 1 - \frac{\mathrm{Var}(\mathbb{E}[Y | X_{-j}])}{\mathrm{Var}(Y)}. \qquad (2) $$
Let E denote the set of all possible types of transitions. By definition, for $\alpha, \beta \in V$, an individual can move from $\alpha$ to $\beta$ if $(\alpha, \beta) \in E$. The pair $(V, E)$ can be identified with a directed graph where $V$ is the set of nodes and $E$ is the set of arrows connecting the compartments between which the individuals can move. Assuming an ordering of the compartments has been chosen, let us identify the set $E$ with a subset of $\mathbb{R}^{|V|}$ as follows: with each $(\alpha, \beta) \in E$, associate the vector of length $|V|$ with components equal to zero except at the positions corresponding to $\alpha$ and $\beta$, where the components are -1 and 1, respectively. The elements of $E$ seen as a subset of $\mathbb{R}^{|V|}$ are called transition vectors.

For every $\alpha \in V$ and a vector of epidemic parameters $\theta \in \Theta \subset \mathbb{R}^d$, let $W_{\alpha}^{\theta}(t)$ be the number of individuals in compartment $\alpha$ at time $t$. Since the population is closed, we have that $\sum_{\alpha \in V} W_{\alpha}^{\theta}(t)$ is constant over time. Denote by $W^{\theta} = (W_{\alpha}^{\theta})_{\alpha \in V}$ the stochastic process that describes the whole population over time. It is assumed that $W^{\theta}$ is a continuous-time Markov chain with state space $\mathcal{E} \subset \mathbb{N}^{|V|}$ and positive rate functions $g_u(\theta, \xi)$, $u \in E$, $\xi \in \mathcal{E}$, given by

$$ \Pr(W^{\theta}(t+s) = \xi + u|W^{\theta}(t) = \xi) = g_u(\theta, \xi)s + o(s) \text{ as } s \to 0. $$

The initial state $W^{\theta}(0) =: \xi_0 \in \mathcal{E}$ is supposed to be fixed. A description of the CTMC of the classical SIR model is given in Example 1.

**Example 1.** The classical SIR model is described as follows. There are three compartments $V = \{S, I, R\}$ and two types of transitions: infection $(S, I)$ and removal $(I, R)$ so that $E = \{(S, I), (I, R)\} = \{(-1, 1, 0), (0, -1, 1)\}$. Infection is characterized by the transition vector $u_{S,I} = (-1, +1, 0)$ and the rate function $g_{S,I} = \frac{\beta}{N}W_I W_S$, where $\beta$ is some parameter and $N$ is the total size of the population. Removal has transition vector $u_{I,R} = (0, -1, +1)$ and rate function $g_{I,R} = \gamma W_I$, where $\gamma$ is some parameter. The vector of parameters is $\theta = (\beta, \gamma)$. The graph of the SIR model is given below:

<figure><img src="image_2.png" /></figure>

## 2.2 Global Sensitivity Analysis

In this section, we first recall the definition of variance-based Sobol' sensitivity indices introduced in Sobol' [1993] for deterministic models with scalar output. Following the paradigm of GSA, we model uncertain inputs by a random vector of independent components $X = (X_1, \dots, X_m)$. Let $E_1, \dots, E_m$ be subsets of $\mathbb{R}$ and $f_1 : E_1 \times \dots \times E_m \to \mathbb{R}$ be some function such that $\mathbb{E} (f_1(X)^2) < +\infty$. Then first-order and total Sobol' indices (see, e.g., Sobol' [1993], Homma and Saltelli [1996]) of the output $Y = f_1(X)$ associated with input $X_j, j = 1, \dots, m$, are
The definition of first-order and total Sobol' indices can be extended to models with vectorial or functional outputs (see, e.g., Lamboni et al. [2011], Gamboa et al. [2014]). Let $(Y_1, \dots, Y_p) := f_2(X_1, \dots, X_m)$ be a vectorial output where $f_2: E_1 \times \dots \times E_m \to \mathbb{R}^p$ is some function and $\mathbb{E} (\|Y\|^2) < +\infty$, with $\| \cdot \|$ denoting the Euclidean norm on $\mathbb{R}^p$. Aggregated first-order and total Sobol' indices are defined as:

$$S_{X_j} = \frac{\sum_{k=1}^{p} \mathrm{Var}(Y_k) S_{X_j,k}}{\sum_{k=1}^{p} \mathrm{Var}(Y_k)}, \quad (3)$$

$$S_{X_j}^{\mathrm{tot}} = \frac{\sum_{k=1}^{p} \mathrm{Var}(Y_k) S_{X_j,k}^{\mathrm{tot}}}{\sum_{k=1}^{p} \mathrm{Var}(Y_k)}, \quad (4)$$

where $S_{X_j,k}$ and $S_{X_j,k}^{\mathrm{tot}}$ are the first-order and total Sobol' indices of the scalar output $Y_k$ associated with the input $X_j$, for $k = 1, \dots, p$ and $j = 1, \dots, m$. If the output of the model of interest is a function of time, it can be reduced to a vectorial output through discretization of time. Then aggregated first-order and total Sobol' indices can be computed using (3) and (4), where the output components $Y_k$ would be the values of the function at the time points of the discretization. Also first-order and total indices defined in Equations (1) and (2) can be computed at each time point of the discretization in order to obtain dynamics of Sobol' indices.

# 3 Complete GSA for stochastic models

In this section we aim at quantifying partial variances due on one hand to uncertain parameters or groups of uncertain parameters, on the other hand to intrinsic randomness, and finally to the interaction between arbitrary groupings of uncertain parameters and intrinsic randomness. We call this a complete GSA. A strategy to achieve this aim requires controlling the latent variables modeling intrinsic randomness.

## 3.1 Deterministic representations of stochastic models

Let $Y^\theta$ denote the random output of some stochastic model with parameters $\theta \in \Theta \subset \mathbb{R}^d$. For instance, $Y^\theta$ might be the stochastic process $W^\theta$ introduced in Section 2.1, or any scalar (or vectorial) quantity of interest defined as a functional of the process $W^\theta$. Note that distinct values of the parameters encoded in $\theta$ correspond to distinct epidemiological patterns. We thus consider the collection $\{Y^\theta, \theta \in \Theta\}$ and assume mutual independence between its members.

As in Section 2.2, uncertain parameters are modeled by a random vector $X = (X_1, \dots, X_d)$ with independent components. In addition, it is assumed that $X$ is independent of $\{Y^\theta, \theta \in \Theta\}$. The pair $(X, Y^X)$ represents the input/output pair that an external observer would see should they draw uncertain parameters at random. The object $Y^\theta$ is then the output observed conditionally on $X = \theta$.
respectively defined as:

$$S_{X_j} = \frac{\mathrm{Var}(\mathbb{E}[Y | X_j])}{\mathrm{Var}(Y)}, \qquad (1)$$

$$S_{X_j}^{\mathrm{tot}} = 1 - \frac{\mathrm{Var}(\mathbb{E}[Y | X_1, \dots, X_{j-1}, X_{j+1}, \dots, X_m])}{\mathrm{Var}(Y)} \\ =: 1 - \frac{\mathrm{Var}(\mathbb{E}[Y | \mathbf{X}_{-j}])}{\mathrm{Var}(Y)}. \qquad (2)$$

The definition of first-order and total Sobol' indices can be extended to models with vectorial or functional output (see, e.g., Lamboni et al. [2011], Gamboa et al. [2014]). Let $Y = (Y_1, \dots, Y_p) := f_2(X_1, \dots, X_m)$ be a vectorial output where $f_2: E_1 \times \dots \times E_m \to \mathbb{R}^p$ is some function and $\mathbb{E}(\|Y\|^2) < +\infty$, with $\|\cdot\|$ denoting the Euclidean norm on $\mathbb{R}^p$. Aggregated first-order and total Sobol' indices are defined as:

$$S_{X_j} = \frac{\sum_{k=1}^{p} \mathrm{Var}(Y_k) S_{X_j, k}}{\sum_{k=1}^{p} \mathrm{Var}(Y_k)}, \qquad (3)$$

$$S_{X_j}^{\mathrm{tot}} = \frac{\sum_{k=1}^{p} \mathrm{Var}(Y_k) S_{X_j, k}^{\mathrm{tot}}}{\sum_{k=1}^{p} \mathrm{Var}(Y_k)}, \qquad (4)$$

where $S_{X_j,k}$ and $S_{X_j,k}^{\mathrm{tot}}$ are the first-order and total Sobol' indices of the scalar output $Y_k$ associated with the input $X_j$, for $k = 1, \dots, p$ and $j = 1, \dots, m$. If the output of the model of interest is a function of time, it can be reduced to a vectorial output through discretization of time. Then aggregated first-order and total Sobol' indices can be computed using (3) and (4), where the output components $Y_k$ would be the values of the function at the time points of the discretization. Also first-order and total indices defined in Equations (1) and (2) can be computed at each time point of the discretization in order to obtain dynamics of Sobol' indices.

# 3 Complete GSA for stochastic models

In this section, we review the methodology of global sensitivity analysis of a stochastic model consisting in representing this model as a deterministic model with input space the set of uncertain parameters augmented with the set of latent variables modeling intrinsic randomness. This methodology permits a quantification of the sensitivity of model output to a variation of the uncertain parameters but also of latent variables modeling intrinsic randomness, and finally of the interaction between both. In the following, we call this a complete GSA. A strategy to achieve this aim requires controlling the latent variables modeling intrinsic randomness.

## 3.1 Deterministic representations of stochastic models

Let $Y^\theta$ denote the random output of some stochastic model with parameters $\theta \in \Theta \subset \mathbb{R}^d$. For instance, $Y^\theta$ might be the stochastic process $W^\theta$ introduced in
It is important to note that in $Y^X$ there are two sources of variability (and hence uncertainty): the one coming from the uncertainty of the parameters (that is, modeled by vector $\mathbf{X}$), and the one coming from the intrinsic randomness of the stochastic model (that is, for a fixed $\theta$ the variability in $Y^\theta$). A complete GSA aims at separating these two sources of uncertainty and quantifying interactions between both.

To achieve this aim, it is necessary to control the latent variables modeling intrinsic randomness. More precisely, we aim at finding a function $f$ and a latent or a set of latent variables $Z$, independent of $\mathbf{X}$, such that the probability distributions of $Y^\theta$ and $f(\theta, Z)$ coincide. Since $\mathbf{X}$ is independent of $Z$ and $Y^\theta$, we immediately have that the input/output pairs $(\mathbf{X}, Y^X)$ and $(\mathbf{X}, f(\mathbf{X}, Z))$ are equal in distribution. The pair $(f, Z)$ is called a *deterministic representation* (or simply a representation) of the stochastic model $Y^\theta$. Often, the function $f$ is the function induced by a (deterministic) algorithm which, if the inputs of that algorithm were drawn from the right distribution, would produce an output statistically equal to the given stochastic model.

From the viewpoint of GSA, one advantage of constructing a deterministic representation of a stochastic model is that standard methods of GSA for deterministic models can be applied straightforwardly. For instance, we can compute first-order and total Sobol' indices by letting $m = d + 1$ and $X_{d+1} = Z$ in (1) and (2) (or in (3) and (4) if the output is vectorial or functional).

In general, there is no unique deterministic representation of a stochastic model. The set of latent variables modeling intrinsic randomness $Z$ and the function $f$ may vary from one representation to the other. More precisely, if $Y^\theta$ is a stochastic model and $(f, Z)$ a representation of it, there may exist another representation $(\tilde{f}, \tilde{Z})$ of $Y^\theta$ such that the laws of $(\mathbf{X}, \tilde{f}(\mathbf{X}, \tilde{Z}))$ and $(\mathbf{X}, f(\mathbf{X}, Z))$ coincide. (Here the probability distribution of $\tilde{Z}$ may be different from that of $Z$.) An example of a toy stochastic model with two different representations is provided in Example 2 below.

**Example 2.** Let $Z \sim \mathcal{U}([0, 1])$ independent of $(\mathbf{X}, Z_1, Z_2) \sim \mathcal{N}(0_\mathbb{R}^3, \text{Id}_3)$. Consider the stochastic model $Y^\theta \sim \mathcal{N}(\theta, 1)$. This model can be represented by using $f(\theta, Z_1, Z_2) = \theta + \frac{1}{\sqrt{2}}(Z_1 + Z_2)$ or $\tilde{f}(\theta, Z) = \theta + \Phi^{-1}(Z)$, where $\Phi$ is the cumulative distribution function of the standard normal distribution.

## 3.2 How does a complete GSA depend on the chosen representation?

As different representations can be exhibited for a same stochastic model, we can wonder how the choice of representations affects GSA results. Let us consider $(f, Z)$ and $(\tilde{f}, \tilde{Z})$ two distinct representations of a same stochastic model with uncertain parameters $\mathbf{X} = (X_1, \cdots, X_d)$ and output $Y^X$. We say that an index $S_{X_j}$ is *representation free* if $S_{X_j}(f, Z) = S_{X_j}(\tilde{f}, \tilde{Z})$, where here $S_{X_j}(f, Z)$ and $S_{X_j}(\tilde{f}, \tilde{Z})$ denote the values of the index $S_{X_j}$ based on the representations $(f, Z)$ and $(\tilde{f}, \tilde{Z})$, respectively.
Section 2.1, or any scalar (or vectorial) quantity of interest defined as a func-
tional of the process $W^θ$. Note that distinct values of the parameters encoded in
$θ$ correspond to distinct epidemiological patterns. We thus consider the collection
{$Y^θ, θ ∈ Θ$} and assume mutual independence between its members.

As in Section 2.2, uncertain parameters are modeled by a random vector $\mathbf{X} = (X_1, \dots, X_d)$ with independent components. In addition, it is assumed that $\mathbf{X}$ is independent of $\{Y^\theta, \theta \in \Theta\}$. The pair $(\mathbf{X}, Y^\mathbf{X})$ represents the input/output pair that an external observer would see should they draw uncertain parameters at random. The object $Y^\theta$ is then the output observed conditionally on $\mathbf{X} = \theta$.

It is important to note that in $Y^\mathbf{X}$ there are two sources of variability (and hence uncertainty): the one coming from the uncertainty of the parameters (that is, modeled by vector $\mathbf{X}$), and the one coming from the intrinsic randomness of the stochastic model (that is, for a fixed $\theta$ the variability in $Y^\theta$). A complete GSA aims at separating these two sources of uncertainty and quantifying interactions between both.

To achieve this aim, it is necessary to control the latent variables modeling in-
trinsic randomness. More precisely, we aim at finding a function $f$ and a latent
or a set of latent variables $Z$, independent of $\mathbf{X}$, such that the probability distri-
butions of $Y^\theta$ and $f(\theta, Z)$ coincide. Since $\mathbf{X}$ is independent of $Z$ and $Y^\theta$, we
immediately have that the input/output pairs $(\mathbf{X}, Y^\mathbf{X})$ and $(\mathbf{X}, f(\mathbf{X}, Z))$ are equal
in distribution. The pair $(f, Z)$ is called a *deterministic representation* (or simply
a representation) of the stochastic model $Y^\theta$. Often, the function $f$ is the function
induced by a (deterministic) algorithm which, if the inputs of that algorithm were
drawn from the right distribution, would produce an output statistically equal to the
given stochastic model.

From the viewpoint of GSA, one advantage of constructing a deterministic rep-
resentation of a stochastic model is that standard methods of GSA for deterministic
models can be applied straightforwardly. For instance, we can compute first-order
and total Sobol' indices by letting $m = d + 1$ and $X_{d+1} = Z$ in (1) and (2) (or
in (3) and (4) if the output is vectorial or functional).

In general, there is no unique deterministic representation of a stochastic model.
The set of latent variables modeling intrinsic randomness $Z$ and the function $f$
may vary from one representation to the other. More precisely, if $Y^\theta$ is a stochas-
tic model and $(f, Z)$ a representation of it, there may exist another representation
$(\tilde{f}, \tilde{Z})$ of $Y^\theta$ such that the laws of $(\mathbf{X}, \tilde{f}(\mathbf{X}, \tilde{Z}))$ and $(\mathbf{X}, f(\mathbf{X}, Z))$ coincide. (Here
the probability distribution of $\tilde{Z}$ may be different from that of $Z$.) An example of
a toy stochastic model with two different representations is provided in Example 2
below.

**Example 2.** Let $Z \sim \mathcal{U}([0, 1])$ independent of $(\mathbf{X}, Z_1, Z_2) \sim \mathcal{N}(0_{\mathbb{R}^3}, \text{Id}_3)$. Consider the stochastic model $Y^\theta \sim \mathcal{N}(\theta, 1)$. This model can be represented by using $f(\theta, Z_1, Z_2) = \theta + \frac{1}{\sqrt{2}}(Z_1 + Z_2)$ or $\tilde{f}(\theta, Z) = \theta + \Phi^{-1}(Z)$, where $\Phi$ is the cumulative distribution function of the standard normal distribution.
**Proposition 1.** First-order Sobol' indices associated with uncertain parameters are representation free.

*Proof.* From Section 3 we know that

$$ (\mathbf{X}, Y^{\mathbf{X}}) \sim (\mathbf{X}, f(\mathbf{X}, Z)) \sim (\mathbf{X}, \tilde{f}(\mathbf{X}, \tilde{Z})). \qquad (5) $$

Thus, for $j = 1, \dots, d$, we have almost surely the following equality:

$$ \mathbb{E}[f(\mathbf{X}, Z) | X_j] = \mathbb{E}[\tilde{f}(\mathbf{X}, \tilde{Z}) | X_j] $$

from which we deduce, using (1), that $S_{X_j}(f, Z) = S_{X_j}(\tilde{f}, \tilde{Z})$, where here $S_{X_j}(f, Z)$ and $S_{X_j}(\tilde{f}, \tilde{Z})$ denote the first-order Sobol' indices associated with $X_j$ based on $(f, Z)$ and $(\tilde{f}, \tilde{Z})$, respectively. $\square$

**Proposition 2.** Total Sobol' indices associated with intrinsic randomness are representation free.

*Proof.* Put $X'_j = X_j$, $j = 1, \dots, d$, $X'_{d+1} = Z$ and $m = d + 1$ so that $\mathbf{X}' := (X'_1, \dots, X'_m) = (\mathbf{X}, Z)$. From (5), we deduce the following almost sure equality:

$$ \mathbb{E}[f(\mathbf{X}, Z) | \mathbf{X}'_{-m}] = \mathbb{E}[f(\mathbf{X}, Z) | \mathbf{X}] = \mathbb{E}[\tilde{f}(\mathbf{X}, \tilde{Z}) | \mathbf{X}] = \mathbb{E}[\tilde{f}(\mathbf{X}, \tilde{Z}) | \mathbf{X}'_{-m}], $$

from which we deduce, using (2), that the total Sobol' indices $S_{\tilde{Z}}^{\text{tot}}(f, Z)$ and $S_{\tilde{Z}}^{\text{tot}}(\tilde{f}, \tilde{Z})$ associated with intrinsic randomness and based respectively on $(f, Z)$ and $(\tilde{f}, \tilde{Z})$ are equal. $\square$

**Proposition 3.** First-order Sobol' indices associated with intrinsic randomness and total Sobol' indices associated with uncertain parameters depend on the choice of representations in general.

To show that Proposition 3 is true, it suffices to exhibit an example where two distinct representations lead to distinct first-order Sobol' indices associated with intrinsic randomness and distinct total Sobol' indices associated with uncertain parameters. Before giving the example, let us give some intuition behind Proposition 3. Note that the random variables $\mathbb{E}[f(\mathbf{X}, Z) | (\mathbf{X}_{-j}, Z)]$ and $\mathbb{E}[\tilde{f}(\mathbf{X}, \tilde{Z}) | (\mathbf{X}_{-j}, \tilde{Z})]$ have different probability distributions in general. Indeed, since $(f, Z) \neq (\tilde{f}, \tilde{Z})$, the way each function $f$ or $\tilde{f}$ combines its (set of) latent variable(s) $Z$ with input $X_j$ to generate outputs may differ. Thus there is no reason for total Sobol' indices associated with uncertain parameters to be representation free. Also, there is no reason for the random variables $\mathbb{E}[f(\mathbf{X}, Z) | Z]$ and $\mathbb{E}[\tilde{f}(\mathbf{X}, \tilde{Z}) | \tilde{Z}]$ to have the same probability distributions and hence the first-order Sobol' index associated with intrinsic randomness to be representation free. This is illustrated on the toy Example 3 below.
## 3.2 How does a complete GSA depend on the chosen representation?

As different representations can be exhibited for a same stochastic model, we can wonder how the choice of representations affects GSA results. Let us consider $(f, Z)$ and $(\tilde{f}, \tilde{Z})$ two distinct representations of a same stochastic model with uncertain parameters $\mathbf{X} = (X_1, \dots, X_d)$ and output $Y^\mathbf{X}$. We say that an index $S_{X_j}$ is representation free if $S_{X_j}(f, Z) = S_{X_j}(\tilde{f}, \tilde{Z})$, where here $S_{X_j}(f, Z)$ and $S_{X_j}(\tilde{f}, \tilde{Z})$ denote the values of the index $S_{X_j}$ based on the representations $(f, Z)$ and $(\tilde{f}, \tilde{Z})$, respectively.

**Proposition 1.** First-order Sobol' indices associated with uncertain parameters are representation free.

*Proof.* From Section 3 we know that

$$ (\mathbf{X}, Y^\mathbf{X}) \sim (\mathbf{X}, f(\mathbf{X}, Z)) \sim (\mathbf{X}, \tilde{f}(\mathbf{X}, \tilde{Z})). \quad (5) $$

Thus, for $j = 1, \dots, d$, we have almost surely the following equality:

$$ \mathbb{E}[f(\mathbf{X}, Z) | X_j] = \mathbb{E}[\tilde{f}(\mathbf{X}, \tilde{Z}) | X_j] $$

from which we deduce, using (1), that $S_{X_j}(f, Z) = S_{X_j}(\tilde{f}, \tilde{Z})$, where here $S_{X_j}(f, Z)$ and $S_{X_j}(\tilde{f}, \tilde{Z})$ denote the first-order Sobol' indices associated with $X_j$ based on $(f, Z)$ and $(\tilde{f}, \tilde{Z})$, respectively. $\square$

**Proposition 2.** Total Sobol' indices associated with intrinsic randomness are representation free.

*Proof.* Put $X'_j = X_j$, $j = 1, \dots, d$, $X'_{d+1} = Z$ and $m = d + 1$ so that $\mathbf{X}' := (X'_1, \dots, X'_m) = (\mathbf{X}, Z)$. From (5), we deduce the following almost sure equality:

$$ \mathbb{E}[f(\mathbf{X}, Z) | \mathbf{X}'_{-m}] = \mathbb{E}[f(\mathbf{X}, Z) | \mathbf{X}] = \mathbb{E}[\tilde{f}(\mathbf{X}, \tilde{Z}) | \mathbf{X}] = \mathbb{E}[\tilde{f}(\mathbf{X}, \tilde{Z}) | \mathbf{X}'_{-m}], $$

from which we deduce, using (2), that the total Sobol' indices $S_Z^{\text{tot}}(f, Z)$ and $S_Z^{\text{tot}}(\tilde{f}, \tilde{Z})$ associated with intrinsic randomness and based respectively on $(f, Z)$ and $(\tilde{f}, \tilde{Z})$ are equal. $\square$

**Proposition 3.** First-order Sobol' indices associated with intrinsic randomness and total Sobol' indices associated with uncertain parameters depend on the choice of representations in general.

To show that Proposition 3 is true, it suffices to exhibit an example where two distinct representations lead to distinct first-order Sobol' indices associated with intrinsic randomness and distinct total Sobol' indices associated with uncertain parameters. Before giving the example, let us give some intuition behind Proposition 3. Note that the random variables $\mathbb{E}[f(\mathbf{X}, Z) | (\mathbf{X}_{-j}, Z)]$ and
**Example 3.** Let $X$ be a random variable independent of $Z$ and $\tilde{Z}$ where $Z$ and $\tilde{Z}$ are i.i.d. under $N(0, 1)$. Define two functions: $f(X, Z) = XZ$ and $\tilde{f}(X, \tilde{Z}) = X^2\tilde{Z}$. If $X$ is distributed such that $\mathbb{P}(X = -1) = \mathbb{P}(X = 1) = \frac{1}{2}$ then $(X, f(X, Z)) \sim (X, \tilde{f}(X, \tilde{Z}))$. Thus, $(f, Z)$ and $(\tilde{f}, \tilde{Z})$ represent the same stochastic model but: $\mathbb{E}[f(X, Z) | Z] = 0$ while $\mathbb{E}[\tilde{f}(X, \tilde{Z}) | \tilde{Z}] = \tilde{Z}$. Simple calculations then lead to $S_Z(f, Z) = 0$ while $S_{\tilde{Z}}(\tilde{f}, \tilde{Z}) = 1$. Total Sobol' indices associated to $X$ can easily be deduced: $S_X^{\text{tot}}(f, Z) = 1 - S_Z(f, Z) = 1$ and $S_X^{\text{tot}}(\tilde{f}, \tilde{Z}) = 1 - S_{\tilde{Z}}(\tilde{f}, \tilde{Z}) = 0$.

We conclude from this section that intrinsic randomness can be modeled in different manners. GSA results naturally depend on the modeling choice. Different modelings bring different insights. Thus the choice of representations must be made with caution, in consultation with epidemiologists. We discuss this point further on a SARS-CoV-2 spread model in Section 5. In the following section, we exhibit different meaningful representations for CTMC stochastic compartmental models, based on different simulation algorithms.

# 4 Deterministic representations for CTMC stochastic compartmental models

Following Section 2.1, let $W^\theta$ be a CTMC stochastic compartmental model with uncertain parameters $\theta$. As explained in Section 3, we wish to rewrite the trajectories of $W^\theta$ as a function of the uncertain parameters and some set of latent variables so as to perform a complete GSA. Based on three different exact simulation algorithms, we propose in this section three different deterministic representations for the same generic CTMC stochastic compartmental model.

One of the first and most basic procedure to simulate the trajectory of a CTMC is as follows. Given that the chain is at some state $W^\theta(s) = \xi$ at time $s$, the holding time until the next jump is distributed as an exponential random variable with parameter $\sum_{u \in E} g_u(\theta, \xi)$ and then the chain moves to state $\xi + u$ with probability $g_u(\theta, \xi) / \sum_{u \in E} g_u(\theta, \xi)$. See, e.g. Karlin and Taylor [1981] for more details. In epidemiology, this procedure is known as Gillespie Direct Method [Gillespie, 1976].

Algorithm 1 is a slightly modification of Gillespie Direct Method based on pseudo-random number generators $RG_j, j = 1, 2$. Here each pseudo-random number generator is seen as an infinite sequence of pseudo-random numbers determined by the seed of the generator. Thus, to get a deterministic representation of the stochastic model $W^\theta$, we let $Z$ represent the value of the seed and $f(\theta, z)$ represent the output of Algorithm 1 for some specific seed value $Z = z$. In practice, the seed values are drawn uniformly over a large number of possible seed values.

We remark that the time-marching method of Algorithm 1 does not require the storage of holding times. Instead, the sequence of holding times is determined on the fly using the pseudo random number generator $RG_1$. In the same way the sequence of types of transition is determined using the random number generator
$\mathbb{E}[\tilde{f}(\mathbf{X}, \tilde{Z}) | (\mathbf{X}_{-j}, \tilde{Z})]$ have different probability distributions in general. Indeed, since $(f, Z) \neq (\tilde{f}, \tilde{Z})$, the way each function $f$ or $\tilde{f}$ combines its (set of) latent variable(s) $Z$ with input $X_j$ to generate the output may differ. Thus there is no reason for total Sobol' indices associated with uncertain parameters to be representation free. Also, there is no reason for the random variables $\mathbb{E}[f(\mathbf{X}, Z) | Z]$ and $\mathbb{E}[\tilde{f}(\mathbf{X}, \tilde{Z}) | \tilde{Z}]$ to have the same probability distributions and hence the first-order Sobol' index associated with intrinsic randomness to be representation free. This is illustrated on the toy Example 3 below.

**Example 3.** Let $X$ be a random variable independent of $Z$ and $\tilde{Z}$ where $Z$ and $\tilde{Z}$ are i.i.d. under $N(0, 1)$. Define two functions: $f(X, Z) = XZ$ and $\tilde{f}(X, \tilde{Z}) = X^2\tilde{Z}$. If $X$ is distributed such that $\mathbb{P}(X = -1) = \mathbb{P}(X = 1) = \frac{1}{2}$ then $(X, f(X, Z)) \sim (X, \tilde{f}(X, \tilde{Z}))$. Thus, $(f, Z)$ and $(\tilde{f}, \tilde{Z})$ represent the same stochastic model but: $\mathbb{E}[f(X, Z) | Z] = 0$ while $\mathbb{E}[\tilde{f}(X, \tilde{Z}) | \tilde{Z}] = \tilde{Z}$. Simple calculations then lead to $S_Z(f, Z) = 0$ while $S_{\tilde{Z}}(\tilde{f}, \tilde{Z}) = 1$. Total Sobol' indices associated to $X$ can easily be deduced: $S_X^{\text{tot}}(f, Z) = 1 - S_Z(f, Z) = 1$ and $S_X^{\text{tot}}(\tilde{f}, \tilde{Z}) = 1 - S_{\tilde{Z}}(\tilde{f}, \tilde{Z}) = 0$.

We conclude from this section that intrinsic randomness can be modeled in different manners. GSA results naturally depend on the modeling choice. Different modelings bring different insights. We discuss this point further on a SARS-CoV-2 spread model in Section 5. In the following section, we exhibit different meaningful representations for CTMC stochastic compartmental models, based on different simulation algorithms.

# 4 Deterministic representations for CTMC stochastic compartmental models

Following Section 2.1, let $W^\theta$ be a CTMC stochastic compartmental model with uncertain parameters $\theta$. As explained in Section 3, we wish to rewrite the trajectories of $W^\theta$ as a function of the uncertain parameters and some set of latent variables so as to perform a complete GSA. Based on three different exact simulation algorithms from the state of the art, we propose in this section three different deterministic representations for the same generic CTMC stochastic compartmental model.

One of the first and most basic procedure to simulate trajectories of a CTMC is as follows. Given that the chain is at some state $W^\theta(s) = \xi$ at time $s$, the holding time until the next jump is distributed as an exponential random variable with parameter $\sum_{u \in E} g_u(\theta, \xi)$ and then the chain moves to state $\xi + u$ with probability $g_u(\theta, \xi) / \sum_{u \in E} g_u(\theta, \xi)$. See, e.g. Karlin and Taylor [1981] for more details. In epidemiology, this procedure is known as Gillespie Direct Method [Gillespie, 1976].

Algorithm 1 is a slight modification of Gillespie Direct Method based on two pseudo-random number generators $RG_1, RG_2$. The first is used to find when the
$RG_2$. In practice, the sequences of holding times and types of transition are entirely determined by the seeds initializing the pseudo-random number generators $RG_j$, $j = 1, 2$. Identical sequences of holding times and types of transition are produced if the same seeds are used. Compared to other algorithms for simulating CTMCs, the main advantage of Algorithm 1 is its low simulation cost.

Algorithm 1 combines all types of transitions to determine the next holding time and then the associated type of transition. However, interpretation of intrinsic randomness in this case is difficult. To overcome this issue, we introduce Algorithm 2 below, which is a slight modification of Gillespie First Reaction Method using pseudo-random number generators. Contrarily to Algorithm 1, Gillespie First Reaction Method analyzes the various types of transitions separately.

Although Gillespie Direct Method and Gillespie First Reaction Method are the most common simulation algorithms in epidemiology, we will compare their results to those obtained from a last simulation algorithm based on the so-called random-time change representation introduced in Kurtz [1982] (see also Ethier and Kurtz [1986]). More precisely, the random state $W^\theta(t)$ can be expressed, for every $t \ge 0$, through

$$W^\theta(t) = W^\theta(0) + \sum_{u \in E} Y_u \left( \int_0^t g_u(\theta, W^\theta(s)) ds \right) u,$$

where $\{Y_u(t), t \ge 0\}$ are independent unit-rate Poisson processes associated with the types of transition $u \in E$. Introduced in Navarro Jimenez et al. [2016] as a tool for doing complete GSAs of chemical reaction network models [Le Maître et al., 2015], Algorithm 3 below is a slight modification of the Modified Next Reaction Method proposed by Anderson [2007].
next transition occurs and the second is used to determine which type of transitions occurs at that time. Each pseudo-random number generator is seen as an infinite sequence of pseudo-random numbers determined by a positive integer called the seed of the generator. For $z = (z_1, z_2)$ a realization of the random seed $Z = (Z_1, Z_2)$, we denote by $f(\theta, z)$ the output of Algorithm 1. Then $(f, Z)$ is a deterministic representation of $W^\theta$ in the sense of Section 3.1. In practice, random seeds $Z_1$ and $Z_2$ are drawn independently from a uniform distribution over a large set of positive integers.

<caption>**Algorithm 1: Gillespie Direct Method**</caption>

**Inputs**: $t_{end}, \theta, Z := (RG_1, RG_2)$
**Data**: $\xi_0, E, \{g_u, u \in E\}$
**Output**: $\{W^\theta(s), s \in [0, t_{end}]\}$
Initialization: $s \leftarrow 0, W^\theta(s) \leftarrow \xi_0$;
**while** $s < t_{end}$ **do**
    $\Sigma \leftarrow \sum_{u \in E} g_u(\theta, W^\theta(s))$;
    Take $r_1$ from $RG_1$;
    $\Delta \leftarrow -\log(r_1)/\Sigma$;
    **for** $u \in E$ **do**
        $p_u \leftarrow g_u(\theta, W^\theta(s)) / \Sigma$;
    **end**
    Divide the interval $(0, 1)$ into $|E|$ sub-intervals of length $p_u$, $u \in E$;
    Take $r_2$ from $RG_2$ and let $\bar{u}$ such that $r_2$ lies within the sub-interval of
    length $p_{\bar{u}}$;
    $W^\theta(s + \Delta) \leftarrow W^\theta(s) + \bar{u}$;
    $s \leftarrow s + \Delta$;

**end**

Using Algorithm 1, it possible to analyze the sensitivity of model output to $Z$, the random vector controlling the intrinsic randomness. However, with this algorithm we cannot conduct a finer sensitivity analysis, by quantifying separately the impact of intrinsic randomness associated to each transition type. To fix this drawback, we introduce Algorithm 2 below, which is a slight modification of Gillespie First Reaction Method using pseudo-random number generators. Algorithm 2 uses as many random numbers as the number of transition types per step. It is thus possible from this algorithm to analyze separately the sensitivity of model output to intrinsic randomness associated to each transition type.
**Algorithm 1: Gillespie Direct Method**

**Inputs**: $t_{\text{end}}, \theta, Z := (RG_1, RG_2)$

**Data**: $\xi_0, \mathbf{E}, \{g_u, u \in \mathbf{E}\}$

**Output**: $\{W^\theta(s), s \in [0, t_{\text{end}}]\}$

Initialization: $s \leftarrow 0, W^\theta(s) \leftarrow \xi_0$;

**while** $s < t_{\text{end}}$ **do**

$\Sigma \leftarrow \sum_{u \in \mathbf{E}} g_u(\theta, W^\theta(s));$

Take $r_1$ from $RG_1$;

$\Delta \leftarrow -\log(r_1)/\Sigma$;

**for** $u \in \mathbf{E}$ **do**

| $p_u \leftarrow g_u(\theta, W^\theta(s)) / \Sigma$;

**end**

Divide the interval $(0, 1)$ into $|\mathbf{E}|$ sub-intervals of length $p_u$, $u \in \mathbf{E}$;

Take $r_2$ from $RG_2$ and let $\bar{u}$ such that $r_2$ lies within the sub-interval of length $p_{\bar{u}}$;

$W^\theta(s + \Delta) \leftarrow W^\theta(s) + \bar{u}$;

$s \leftarrow s + \Delta$;

**end**

**Algorithm 2: Gillespie First Reaction Method**

**Inputs**: $t_{\text{end}}, \theta, Z := \{RG_u, u \in \mathbf{E}\}$

**Data**: $\xi_0, \mathbf{E}, \{g_u, u \in \mathbf{E}\}$

**Output**: $\{W^\theta(s), s \in [0, t_{\text{end}}]\}$

Initialization: $s \leftarrow 0, W^\theta(s) \leftarrow \xi_0$;

**while** $s < t_{\text{end}}$ **do**

**for** $u \in \mathbf{E}$ **do**

Take $r_u$ from $RG_u$;

$a_u \leftarrow g_u(\theta, W^\theta(s));$

$\Delta_u \leftarrow \frac{-\log(r_u)}{a_u}$

**end**

$\bar{u} \leftarrow \operatorname{argmin}_u \Delta_u$;

$\Delta \leftarrow \Delta_{\bar{u}}$;

$W(s + \Delta) \leftarrow W(s) + \bar{u}$;

$s \leftarrow s + \Delta$;

**end**
<caption>**Algorithm 2: Gillespie First Reaction Method**</caption>

**Inputs**: $t_{end}$, $\theta$, $Z := \{RG_u, u \in E\}$
**Data**: $\xi_0$, $E$, $\{g_u, u \in E\}$
**Output**: $\{W^\theta(s), s \in [0, t_{end}]\}$
Initialization: $s \leftarrow 0, W^\theta(s) \leftarrow \xi_0$;
**while** $s < t_{end}$ **do**
| **for** $u \in E$ **do**
| | Take $r_u$ from $RG_u$;
| | $a_u \leftarrow g_u(\theta, W^\theta(s))$;
| | $\Delta_u \leftarrow \frac{-\log(r_u)}{a_u}$
| **end**
| $\bar{u} \leftarrow \text{argmin}_u \Delta_u$;
| $\Delta \leftarrow \Delta_{\bar{u}}$;
| $W(s + \Delta) \leftarrow W(s) + \bar{u}$;
| $s \leftarrow s + \Delta$;
| **end**

As an alternative to Gillespie First Reaction Method, we review the simulation algorithm introduced in Kurtz [1982] (see also Ethier and Kurtz [1986]), based on the so-called random-time change representation. More precisely, the random state $W^\theta(t)$ can be expressed, for every $t \ge 0$, through

$$W^\theta(t) = W^\theta(0) + \sum_{u \in E} Y_u \left( \int_0^t g_u(\theta, W^\theta(s)) ds \right) u,$$

where $\{Y_u(t), t \ge 0\}$, $u \in E$, are independent unit-rate Poisson processes associated with transition types $u \in E$. The next reaction method generates exact sample paths while only needing one random number per step. Introduced in Navarro Jimenez et al. [2016] as a tool for doing complete GSA of chemical reaction network models [Le Maître et al., 2015], Algorithm 3 below is a slight modification of the Modified Next Reaction Method proposed by Anderson [2007].
<caption>**Algorithm 3: Modified Next Reaction Method**</caption>

**Inputs**: $t_{\text{end}}, \theta, Z := \{RG_u, u \in E\}$
**Data**: $\xi_0, E, \{g_u, u \in E\}$
**Output**: $\{W^\theta(s), s \in [0, t_{\text{end}}]\}$
Initialization:
**for** $u \in E$ **do**
> Take $r_u$ from $RG_u$;
> $t_u \leftarrow 0$, $t_u^+ \leftarrow -\log(r_u)$;

**end**

$s \leftarrow 0, W^\theta(s) \leftarrow \xi_0$;

**while** $s < t_{\text{end}}$ **do**

**for** $u \in E$ **do**

| $a_u \leftarrow g_u(\theta, W(s)); \Delta_u \leftarrow \frac{t_u^+ - t_u}{a_u}$

**end**

$\bar{u} \leftarrow \operatorname{argmin}_u \Delta_u$;

$\Delta \leftarrow \Delta_{\bar{u}}$;

$W(s + \Delta) \leftarrow W(s) + \bar{u}$;

$s \leftarrow s + \Delta$;

**for** $u \in E$ **do**

| $t_u \leftarrow t_u + a_u \Delta$

**end**

Take $r_{\bar{u}}$ from $RG_{\bar{u}}$;

$t_{\bar{u}}^{+} \leftarrow t_{\bar{u}}^{+} - \log(r_{\bar{u}})$;

**end**

In summary, in this section we introduced three deterministic representations of compartmental models based on CTMCs using three different stochastic simulation algorithms: Gillespie Direct Method (Algorithm 1), Gillespie First Reaction Method (Algorithm 2) and Modified Next Reaction Method (Algorithm 3). As they analyze the various types of transitions separately, Algorithms 2 and 3 are more insightful, but they are also more computationally expensive than Algorithm 1.

The aim of Section 5 will be to implement a complete GSA using each of the above representations on a SARS-CoV-2 spread model. From the theoretical results of Section 3.2, we expect GSA results to differ from one representation to the other, with the take home message that the representation must be chosen with caution and in consultation with epidemiologists.

# 5 Application to a SARS-CoV-2 spread model

Although our approach for a complete GSA of stochastic compartmental models is generic, we propose in this section, as a case study, to do a sensitivity analysis of a parsimonious SARS-CoV-2 spread model, which is a simplified but still resalistic version of the model introduced in Cazelles et al. [2021]. In this section, we do not
<caption>**Algorithm 3: Modified Next Reaction Method**</caption>

**Inputs**: $t_{\text{end}}, \theta, Z := \{RG_u, u \in E\}$
**Data**: $\xi_0, E, \{g_u, u \in E\}$
**Output**: $\{W^\theta(s), s \in [0, t_{\text{end}}]\}$
**Initialization**:
**for** $u \in E$ **do**
| Take $r_u$ from $RG_u$;
| $t_u \leftarrow 0$, $t_u^+ \leftarrow -\log(r_u)$;

**end**

$s \leftarrow 0, W^\theta(s) \leftarrow \xi_0;$

**while** $s < t_{\text{end}}$ **do**

**for** $u \in E$ **do**

| $a_u \leftarrow g_u(\theta, W(s)); \Delta_u \leftarrow \frac{t_u^+ - t_u}{a_u}$

**end**

$\bar{u} \leftarrow \operatorname{argmin}_u \Delta_u;$

$\Delta \leftarrow \Delta_{\bar{u}};$

$W(s + \Delta) \leftarrow W(s) + \bar{u};$

$s \leftarrow s + \Delta;$

**for** $u \in E$ **do**

| $t_u \leftarrow t_u + a_u \Delta$

**end**

Take $r_{\bar{u}}$ from $RG_{\bar{u}}$;

$t_{\bar{u}}^+ \leftarrow t_{\bar{u}}^+ - \log(r_{\bar{u}});$

**end**

We introduced in this section three deterministic representations of compart-
mental models, each one based on a slight modification of an exact simulation al-
gorithm: Gillespie Direct Method (Algorithm 1), Gillespie First Reaction Method
(Algorithm 2) and Modified Next Reaction Method (Algorithm 3). Now, the aim
of Section 5 is to implement a complete GSA using each of the above representa-
tions on a SARS-CoV-2 spread model. Depending on the number of compartments
in our model, the computational cost may vary from one algorithm to the other.
However for this specific SARS-CoV-2 model, simulation times were comparable.
From the theoretical results of Section 3.2, we expect GSA results to differ from
one representation to the other.

5 Application to a SARS-CoV-2 spread model

We propose in this section, as a case study, to apply the methodology of global sen-
sitivity analysis reviewed in the previous sections to a parsimonious SARS-CoV-2
spread model, which is a simplified but still realistic version of the model intro-
duced in Cazelles et al. [2021]. We do not pretend to provide the most suitable
model for the propagation of SARS-CoV-2, we rather aim at demonstrating the ef-
pretend to provide the most suitable model for the propagation of SARS-CoV-2, we
rather aim at demonstrating the effectiveness of the approach presented in Section
3 for a complete GSA of stochastic compartmental models. In order to illustrate the
statement in Section 3.2 that sensitivity analysis results depend on the determinis-
tic representation chosen for its implementation, we compare the results by using
each of Modified Gillespie Direct Method (Algorithm 1), Modified Gillespie First
Reaction Method (Algorithm 2) and Modified Next Reaction Method (Algorithm
3) for simulations. Recall that these algorithms have been presented in Section
4. In Section 5.1 we describe the considered SARS-CoV-2 model. Then in Sec-
tion 5.2 we introduce the quantities of interest and detail our numerical setting for
sensitivity analysis. Finally in Section 5.3 we present the results of the sensitivity
analyses obtained from the different simulation algorithms of Section 4.

**5.1 A SARS-CoV-2 spread model**

Recall from Section 2 that each process $W_{\alpha}^{\theta}, \alpha \in \mathbf{V}$, counts the number of individuals in compartment $\alpha$ over time. In this section we let $\mathbf{V} = \{S, E, A, I, H, R, D\}$, where the seven compartments represent seven possible health statuses: an individual can be susceptible (S), exposed (E) (i.e. infected but not yet infectious), asymptomatic infectious (A), symptomatic infectious (I), hospitalized (H), recovered (R) or dead (D). There are nine possible types of transition between these compartments, see Figure 1. Note that infection is neglected within hospitals so that hospitalized individuals cannot infect. Moreover, it is assumed that recovered individuals get perfectly immunized so they cannot be susceptible after recovering. The vector of uncertain parameters is given by $\theta = (\beta, \gamma_E, \gamma_A, \gamma_I, \gamma_H, p_{(E,A)}, p_C, p_{D|C}, p_{(H,D)})$. The different types of transition and their characteristics (transition vector $\mathbf{u}$ and associated rate function $g_{\mathbf{u}}$) are described in Table 1.
fectiveness of the approach presented in Section 3 for a complete GSA of stochastic compartmental models. In order to illustrate the statement in Section 3.2 that sensitivity analysis results depend on the deterministic representation chosen for its implementation, we compare the results by using each of Modified Gillespie Direct Method (Algorithm 1), Modified Gillespie First Reaction Method (Algorithm 2) and Modified Next Reaction Method (Algorithm 3) for simulations. Recall that these algorithms have been reviewed in Section 4. In Section 5.1 we describe the considered SARS-CoV-2 model. Then in Section 5.2 we introduce the quantities of interest and detail our numerical setting for sensitivity analysis. Finally in Section 5.3 we present the results of the sensitivity analyses obtained from the different simulation algorithms of Section 4. The code developed to perform the numerical experiments is available at https://hal.inrae.fr/MATHNUM/hal-03565729.

## 5.1 A SARS-CoV-2 spread model

Recall from Section 2 that each process $W_{\alpha}^{\theta}, \alpha \in V$, counts the number of individuals in compartment $\alpha$ over time. In this section we let $V = \{S, E, A, I, H, R, D\}$, where the seven compartments represent seven possible health statuses: an individual can be susceptible (S), exposed (E) (i.e. infected but not yet infectious), asymptomatic infectious (A), symptomatic infectious (I), hospitalized (H), recovered (R) or dead (D). There are nine possible types of transition between these compartments, see Figure 1. Note that infection is neglected within hospitals so that hospitalized individuals cannot infect. Moreover, it is assumed that recovered individuals get perfectly immunized so they cannot be susceptible after recovering. The vector of uncertain parameters is given by $\theta = (\beta, \gamma_E, \gamma_A, \gamma_I, \gamma_H, p_{(E,A)}, p_C, p_{D|C}, p_{(H,D)})$. The different types of transition and their characteristics (transition vector $\mathbf{u}$ and associated rate function $g_{\mathbf{u}}$) are described in Table 1.
<figure><img src="image_3.png" /><figcaption>Figure 1: Compartmental model of the spread of SARS-CoV-2. The nodes are the possible health statuses and the arrows connecting them are the possible types of transition. The labels above the arrows are the corresponding rate functions.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 1: Compartmental model of the spread of SARS-CoV-2. The nodes are the possible health statuses and the arrows connecting them are the possible types of transition. The labels above the arrows are the corresponding rate functions.</figcaption></figure>

The list of uncertain epidemic parameters is given in Table 2. Instead of pa-
rameterizing the transitions from *I* to *R*, *H* or *D* with the probabilities *p*<sub>IR</sub>, *p*<sub>IH</sub>
and *p*<sub>ID</sub>, that correspond to the probability for an individual in *I* to recover, to be
hospitalized or to die, that are clearly not independent as *p*<sub>IH</sub> + *p*<sub>ID</sub> + *p*<sub>IR</sub> = 1,
we chose to use instead parameters *p*<sub>C</sub> corresponding to the probability for an in-
dividual in compartment *I* of being in a critical state and *p*<sub>D|C</sub> corresponding to
the probability to die without being hospitalized conditionally on being in a critical
state. These last two parameters are considered as independent from each other.
This re-parameterization was inspired, e.g., from what is done in [Da Veiga et al.,
2021, Chapter 7, page 191]. More generally all the parameters listed in Table
2 are assumed to be independent. The nominal value and range of variation for
each parameter has been chosen in agreement with current knowledge and rep-
resent at least plausible values; see, e.g., [Knock et al., 2021, Table S2 on page
15 of the Supplementary Material] and, specifically for parameters *p*<sub>C</sub> and *p*<sub>D|C</sub>,
[Da Veiga et al., 2021, Chapter 7].
<table><caption>Table 1: Description of the model transitions between states {S, E, A, I, H, R, D}.</caption><thead><tr><th>Transition</th><th>Type</th><th>Transition vector <b>u</b></th><th>Rate function g<sub>u</sub></th></tr></thead><tbody><tr><td>(S, E)</td><td>infection</td><td>(-1, 1, 0, 0, 0, 0)</td><td><sup>β</sup>/<sub>N</sub> &middot; W<sub>S</sub> &middot; (W<sub>A</sub> + W<sub>I</sub>)</td></tr><tr><td>(E, A)</td><td>asymptomatic<br>infectiousness activation</td><td>(0, -1, +1, 0, 0, 0)</td><td>γ<sub>E</sub> &middot; p<sub>(E,A)</sub> &middot; W<sub>E</sub></td></tr><tr><td>(E, I)</td><td>symptomatic<br>infectiousness activation</td><td>(0, -1, 0, +1, 0, 0)</td><td>γ<sub>E</sub> &middot; (1 - p<sub>(E,A)</sub>) &middot; W<sub>E</sub></td></tr><tr><td>(A, R)</td><td>recovery of an<br>asymptomatic</td><td>(0, 0, -1, 0, 0, +1)</td><td>γ<sub>A</sub> &middot; W<sub>A</sub></td></tr><tr><td>(I, R)</td><td>recovery of a<br>symptomatic</td><td>(0, 0, 0, -1, 0, +1)</td><td>γ<sub>I</sub> &middot; (1 - p<sub>C</sub>) &middot; W<sub>I</sub></td></tr><tr><td>(I, H)</td><td>hospitalization<br>of a symptomatic</td><td>(0, 0, 0, -1, +1, 0)</td><td>γ<sub>I</sub> &middot; p<sub>C</sub> &middot; (1 - p<sub>D|C</sub>) &middot; W<sub>I</sub></td></tr><tr><td>(I, D)</td><td>death of a<br>symptomatic</td><td>(0, 0, 0, -1, 0, 0, +1)</td><td>γ<sub>I</sub> &middot; p<sub>C</sub> &middot; p<sub>D|C</sub> &middot; W<sub>I</sub></td></tr><tr><td>(H, R)</td><td>recovery of a<br>hospitalized</td><td>(0, 0, 0, 0, -1, +1, 0)</td><td>γ<sub>H</sub> &middot; (1 - p<sub>(H,D)</sub>) &middot; W<sub>H</sub></td></tr><tr><td>(H, D)</td><td>death of a<br>hospitalized</td><td>(0, 0, 0, 0, -1, 0, +1)</td><td>γ<sub>H</sub> &middot; p<sub>(H,D)</sub> &middot; W<sub>H</sub></td></tr></tbody></table>

<table><caption>Table 2: Model parameter nominal values and their range of variation in the sensitivity analysis.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Nominal value</th><th>Range of variation</th></tr></thead><tbody><tr><td>β</td><td>transmission rate</td><td>2.175</td><td>(0.35, 4)</td></tr><tr><td>1/γ<sub>E</sub></td><td>mean sojourn duration in E</td><td>4.5 days</td><td>(2, 7)</td></tr><tr><td>1/γ<sub>A</sub></td><td>mean sojourn duration in A</td><td>2 days</td><td>(1, 3)</td></tr><tr><td>1/γ<sub>I</sub></td><td>mean sojourn duration in I</td><td>4 days</td><td>(3, 5)</td></tr><tr><td>1/γ<sub>H</sub></td><td>mean sojourn duration in H</td><td>9.5 days</td><td>(7, 12)</td></tr><tr><td>p<sub>(E,A)</sub></td><td>probability for an exposed to become asymptomatic</td><td>0.5</td><td>(0.3, 0.7)</td></tr><tr><td>p<sub>C</sub></td><td>probability for an individual in compartment I of being in a critical state</td><td>0.175</td><td>(0.15, 0.2)</td></tr><tr><td>p<sub>D|C</sub></td><td>probability to die without being hospitalized knowing that the individual is in a critical state</td><td>0.175</td><td>(0.15, 0.2)</td></tr><tr><td>p<sub>(H,D)</sub></td><td>probability for a hospitalized to die</td><td>0.0505</td><td>(0.001, 0.1)</td></tr></tbody></table>

## 5.2 Setting for sensitivity analysis

We consider a population of $N = 2005$ individuals including five exposed individuals at the start of the epidemic $t = 0$, so that the process $W^θ$ has the initial state
<table><caption>Table 1: Description of the model transitions between states {S, E, A, I, H, R, D}.</caption><thead><tr><th>Transition</th><th>Type</th><th>Transition vector <b>u</b></th><th>Rate function g<sub>u</sub></th></tr></thead><tbody><tr><td>(S, E)</td><td>infection</td><td>(-1, 1, 0, 0, 0, 0)</td><td><sup>β</sup>/<sub>N</sub> &middot; W<sub>S</sub> &middot; (W<sub>A</sub> + W<sub>I</sub>)</td></tr><tr><td>(E, A)</td><td>asymptomatic<br>infectiousness activation</td><td>(0, -1, +1, 0, 0, 0)</td><td>γ<sub>E</sub> &middot; p<sub>(E,A)</sub> &middot; W<sub>E</sub></td></tr><tr><td>(E, I)</td><td>symptomatic<br>infectiousness activation</td><td>(0, -1, 0, +1, 0, 0)</td><td>γ<sub>E</sub> &middot; (1 - p<sub>(E,A)</sub>) &middot; W<sub>E</sub></td></tr><tr><td>(A, R)</td><td>recovery of an<br>asymptomatic</td><td>(0, 0, -1, 0, 0, +1)</td><td>γ<sub>A</sub> &middot; W<sub>A</sub></td></tr><tr><td>(I, R)</td><td>recovery of a<br>symptomatic</td><td>(0, 0, 0, -1, 0, +1)</td><td>γ<sub>I</sub> &middot; (1 - p<sub>C</sub>) &middot; W<sub>I</sub></td></tr><tr><td>(I, H)</td><td>hospitalization<br>of a symptomatic</td><td>(0, 0, 0, -1, +1, 0)</td><td>γ<sub>I</sub> &middot; p<sub>C</sub> &middot; (1 - p<sub>D|C</sub>) &middot; W<sub>I</sub></td></tr><tr><td>(I, D)</td><td>death of a<br>symptomatic</td><td>(0, 0, 0, -1, 0, 0, +1)</td><td>γ<sub>I</sub> &middot; p<sub>C</sub> &middot; p<sub>D|C</sub> &middot; W<sub>I</sub></td></tr><tr><td>(H, R)</td><td>recovery of a<br>hospitalized</td><td>(0, 0, 0, 0, -1, +1, 0)</td><td>γ<sub>H</sub> &middot; (1 - p<sub>(H,D)</sub>) &middot; W<sub>H</sub></td></tr><tr><td>(H, D)</td><td>death of a<br>hospitalized</td><td>(0, 0, 0, 0, -1, 0, +1)</td><td>γ<sub>H</sub> &middot; p<sub>(H,D)</sub> &middot; W<sub>H</sub></td></tr></tbody></table>

<table><caption>Table 2: Model parameter nominal values and their range of variation in the sensitivity analysis.</caption><thead><tr><th>Parameter</th><th>Description</th><th>Nominal value</th><th>Range of variation</th></tr></thead><tbody><tr><td>β</td><td>transmission rate</td><td>2.175</td><td>(0.35, 4)</td></tr><tr><td>1/γ<sub>E</sub></td><td>mean sojourn duration in E</td><td>4.5 days</td><td>(2, 7)</td></tr><tr><td>1/γ<sub>A</sub></td><td>mean sojourn duration in A</td><td>2 days</td><td>(1, 3)</td></tr><tr><td>1/γ<sub>I</sub></td><td>mean sojourn duration in I</td><td>4 days</td><td>(3, 5)</td></tr><tr><td>1/γ<sub>H</sub></td><td>mean sojourn duration in H</td><td>9.5 days</td><td>(7, 12)</td></tr><tr><td>p<sub>(E,A)</sub></td><td>probability for an exposed to become asymptomatic</td><td>0.5</td><td>(0.3, 0.7)</td></tr><tr><td>p<sub>C</sub></td><td>probability for an individual in compartment I of being in a critical state</td><td>0.175</td><td>(0.15, 0.2)</td></tr><tr><td>p<sub>D|C</sub></td><td>probability to die without being hospitalized knowing that the individual is in a critical state</td><td>0.175</td><td>(0.15, 0.2)</td></tr><tr><td>p<sub>(H,D)</sub></td><td>probability for a hospitalized to die</td><td>0.0505</td><td>(0.001, 0.1)</td></tr></tbody></table>

## 5.2 Setting for sensitivity analysis

We consider a population of $N = 2005$ individuals including five exposed individuals at the start of the epidemic $t = 0$, so that the process $W^θ$ has the initial state
$$ \xi_0 = (W_S(0) = 2000, W_E(0) = 5, W_A(0) = W_I(0) = W_H(0) = W_R(0) = W_D(0) = 0). $$

We focus on two quantities of interest (QoIs). First we consider a scalar QoI, namely the extinction time $Y_{ext}^{\theta}$ of the epidemic, defined as the first instant at which there are no exposed (E) nor infectious (A or I) individuals anymore:

$$ Y_{ext}^{\theta} = \inf\{t \ge 0 : W_{E}^{\theta}(t) + W_{A}^{\theta}(t) + W_{I}^{\theta}(t) = 0\}. $$

Note that for all $\theta \in \Theta$, $Y_{ext}^{\theta}$ is well-defined, i.e. $Y_{ext}^{\theta} < +\infty$. Indeed, by considering the compartmental model described in Figure 1, after a finite number of transitions, the stochastic process will necessarily reach an absorbing state with empty compartments E, A and I. We display in Figure 2 two hundred independent realizations of $Y_{ext}^{\theta}$ for each of the simulation algorithms of Section 4. The uncertain parameters $\theta$ were set to the nominal values given in Table 2. These values agree with current knowledge and are considered as plausible (see, e.g., Knock et al. [2021], Da Veiga et al. [2021]). The three boxplots are similar, which was expected since the distributions of the processes returned by each of the three algorithms are the same.

<figure><img src="image_4.png" /><figcaption>Figure 2: Boxplot of 200 simulations of $Y_{ext}^{\theta}$ performed with (left, red) Gillespie Direct Method (middle, green) Gillespie First Reaction (right, blue) Modified Next Reaction algorithm with uncertain parameters $\theta$ set to the nominal values given in Table 2.</figcaption></figure>

The second QoI we are considering is the dynamic of the number of symptomatic infectious individuals:

$$ Y_{I}^{\theta} = \{W_{I}^{\theta}(t), t \in [0, t_{end}]\}, $$

where $t_{end}$ was set to 60. (The process dies out at around that time, see Figure 3.) We display in Figure 3 twenty independent realizations of the process $Y_{I}^{\theta}$ for each of the simulation algorithms of Section 4. The input parameter vector $\theta$ was set to
$$ \xi_0 = (W_S(0) = 2000, W_E(0) = 5, \\ W_A(0) = W_I(0) = W_H(0) = W_R(0) = W_D(0) = 0). $$

We focus on two quantities of interest (QoIs). First we consider a scalar QoI, namely the extinction time $Y_{ext}^{\theta}$ of the epidemic, defined as the first instant at which there are no exposed (E) nor infectious (A or I) individuals anymore:

$$ Y_{ext}^{\theta} = \inf\{t \ge 0 : W_{E}^{\theta}(t) + W_{A}^{\theta}(t) + W_{I}^{\theta}(t) = 0\}. $$

Note that for all $\theta \in \Theta$, $Y_{ext}^{\theta}$ is well-defined, i.e. $Y_{ext}^{\theta} < +\infty$. Indeed, by considering the compartmental model described in Figure 1, after a finite number of transitions, the stochastic process will necessarily reach an absorbing state with empty compartments E, A and I. We display in Figure 2 two hundred independent realizations of $Y_{ext}^{\theta}$ for each of the simulation algorithms of Section 4. The uncertain parameters $\theta$ were set to the nominal values given in Table 2. The three boxplots are similar, which was expected since the distributions of the processes returned by each of the three algorithms are the same.

<figure><img src="image_4.png" /><figcaption>Figure 2: Boxplot of 200 simulations of $Y_{ext}^{\theta}$ performed with (left, red) Gillespie Direct Method (middle, green) Gillespie First Reaction (right, blue) Modified Next Reaction algorithm with uncertain parameters $\theta$ set to the nominal values given in Table 2.</figcaption></figure>

The second QoI we are considering is the dynamic of the number of symptomatic infectious individuals:

$$ Y_{I}^{\theta} = \{W_{I}^{\theta}(t), t \in [0, t_{end}]\}, $$

where $t_{end}$ was set to 60. (The process dies out at around that time, see Figure 3.) We display in Figure 3 twenty independent realizations of the process $Y_{I}^{\theta}$ for each of the simulation algorithms of Section 4. The input parameter vector $\theta$ was set to the nominal values given in Table 2. The three charts in Figure 3 display similar sample paths for $Y_{I}^{\theta}$, which was expected since the distribution of the processes returned by each of the three algorithms is the same.
the nominal values given in Table 2. The three charts in Figure 3 display similar sample paths for $Y_t^\theta$, which was expected since the distribution of the processes returned by each of the three algorithms is the same.

<figure><img src="image_5.png" /><figcaption>Figure 3: 20 independent realizations of $t \rightarrow W_t^\theta(t)$, with $\theta$ set to nominal values given in Table 2, from (left) Gillespie Direct Method (middle) Gillespie First Reaction (right) Modified Next Reaction algorithm.</figcaption></figure>

As expected, for $\theta$ fixed to its nominal value, we do not observe on Figure 3 a significant difference in the simulations run with Algorithm 1, 2 or 3.

In practice, simulations are carried out using the R Statistical Software [R Core Team, 2021]. Sensitivity indices are estimated by using the R package sensitivity [Iooss et al., 2021]. The function soboljansen() is used for total Sobol' index estimation while sobol2007() is used for first-order Sobol' index estimation. A priori distributions for model parameters are uniform distributions as described in Table 2 and the a priori for intrinsic randomness is modeled by seeds uniformly distributed in $\{1, \dots, 10^9\}$. Sensitivity indices are estimated from two independent designs of $n = 2000$ input-output samples, where for each sample a trajectory of $W^\theta$ is simulated through either Algorithm 1, Algorithm 2 or Algorithm 3. As the dimension of the input space is large (at least 7 model parameters plus inputs modeling intrinsic randomness whose number depends on the simulation algorithm), we use Latin Hypercube Sampling (see, e.g., Lin and Tang [2015]). Latin Hypercube Samples are generated by using the R package DiceDesign [Dupuy et al., 2015].

## 5.3 Sensitivity analysis results

This section is devoted to the presentation and comparison of sensitivity analysis results obtained for the algorithms presented in Section 4, namely Gillespie Direct Method, Gillespie First Reaction Method and Modified Next Reaction Method. In Section 5.3.1 we present the results for the scalar output of interest, namely the extinction time of the epidemic $Y_{ext}^\theta$. Then in Section 5.3.2 we present the sensitivity analysis results for the functional output corresponding to the dynamic of the number of symptomatic infectious individuals $Y_t^\theta$. Finally in Section 5.4 we discuss the choice of algorithms, depending on the practitioner's objectives. Recall that in
<figure><img src="image_5.png" /><figcaption>Figure 3: First row: 20 independent realizations of $t \rightarrow W_I^{\theta}(t)$. Second row: empirical quantile of level 0.1, \dots, and 0.9 as a function of time, computations are based on 1000 independent sample paths of $t \rightarrow W_I^{\theta}(t)$. For each row, the components of the parameter vector $\theta$ were set to their nominal value (see Table 2) and simulations were obtained from Gillespie Direct Method (left), Gillespie First Reaction Method (middle), Modified Next Reaction Method (right).</figcaption></figure>

In the top row of Figure 3, we show 20 independent trajectories of the process $t \mapsto W_I^{\theta}(t)$, obtained from Gillespie Direct Method (left), Gillespie First Reaction Method (middle) and Modified Next Reaction Method (right) with the components of the parameter vector $\theta$ fixed to their nominal value (see Table 2). On the bottom row of the same figure, we show the evolution over time of quantiles of different order (from 0.1 to 0.9) calculated from 1000 independent trajectories of each algorithm (from left to right). These plots are in coherence with the fact that Algorithms 1, 2 and 3 are all exact simulation algorithms of the same stochastic process.

In practice, simulations are carried out using the R Statistical Software [R Core Team, 2021]. Sensitivity indices are estimated by using the R package sensitivity [Iooss et al., 2021]. The function soboljansen() is used for total Sobol' index estimation while sobol2007() is used for first-order Sobol' index estimation. A priori distributions for model parameters are uniform distributions as described in Table 2 and the a priori for intrinsic randomness is modeled by seeds uniformly distributed in $\{1, \dots, 10^9\}$. Sensitivity indices are estimated from two independent designs of $n = 2000$ input-output samples, where for each sample a trajectory of $W^\theta$ is sim-
all the results presented in this section, sensitivity indices were estimated from two independent designs of $n = 2000$ input-output samples, and the estimation was repeated independently 50 times for the different boxplots.

### 5.3.1 Sensitivity analysis results for $Y_{ext}^{\theta}$

We display on Figure 4 boxplots of first-order and total Sobol' index estimates.

<figure><img src="image_6.png" /><figcaption>Figure 4: Boxplots of 50 independent estimates of the (top) first-order and (bottom) total Sobol' indices for $Y_{ext}^{\theta}$, each computed with two independent designs of $n = 2000$ input-output samples, for each simulation algorithm: (red, left) Gillespie Direct Method, (green, middle) Gillespie First Reaction, (blue, right) Modified Next Reaction.</figcaption></figure>

In accordance with the results stated in Section 3.2, we observe on the top of Figure 4 that there are no significant differences between the three algorithms for
ulated through either Algorithm 1, Algorithm 2 or Algorithm 3. As the dimension
of the input space is large (at least 7 model parameters plus inputs modeling intrin-
sic randomness whose number depends on the simulation algorithm), we use Latin
Hypercube Sampling (see, e.g., Lin and Tang [2015]). Latin Hypercube Samples
are generated by using the R package DiceDesign [Dupuy et al., 2015].

**5.3 Sensitivity analysis results**

This section is devoted to the presentation and comparison of sensitivity analysis
results obtained for the algorithms presented in Section 4, namely Gillespie Direct
Method, Gillespie First Reaction Method and Modified Next Reaction Method. In
Section 5.3.1 we present the results for the scalar output of interest, namely the ex-
tinction time of the epidemy $Y_{ext}^{\theta}$. Then in Section 5.3.2 we present the sensitivity
analysis results for the functional output corresponding to the dynamic of the num-
ber of symptomatic infectious individuals $Y_{I}^{\theta}$. Finally in Section 5.4 we discuss
the choice of algorithms, depending on the practitioner’s objectives. Recall that in
all the results presented in this section, sensitivity indices were estimated from two
independent designs of $n = 2000$ input-output samples, and the estimation was
repeated independently 50 times for the different boxplots.

**5.3.1 Sensitivity analysis results for $Y_{ext}^{\theta}$**

We display on Figure 4 boxplots of first-order and total Sobol' index estimates.
the first-order Sobol' index estimates. The only input parameters with a significant first-order effect are $\beta$ and $\gamma_E$. The sum of the first-order index estimates is far below 1 which means that interactions are not negligible. We observe on the bottom of Figure 4 that almost all inputs have a total effect significantly greater than zero. The interaction strength varies from one simulation algorithm to the other. This is due to the fact that the modeling of intrinsic randomness depends on each simulation algorithm. In particular, we observe that the total Sobol' index estimates corresponding to the Modified Next Reaction algorithm are never less than their counterpart computed from the First Reaction Method algorithm. This reflects a stronger interaction with intrinsic random noise for Modified Next Reaction algorithm. Finally, as expected from the theoretical results in Section 3.2, the total index estimates associated with intrinsic randomness do not depend on the chosen simulation algorithm.

### 5.3.2 Sensitivity analysis results for $Y_i^\theta$

Since $Y_i^\theta$ is a dynamical process, we can consider the sensitivity of the whole trajectory or the sensitivity time by time. In the numerical experiments, $Y_i^\theta$ is discretized over a regular grid of size 1000 of the interval $[0, t_{end}]$. The sensitivity of the whole trajectory consists of computing estimates of the aggregated sensitivity indices of Section 2.2. These provide a scalar summary for the dynamical evolution of first-order and total Sobol' indices. They are displayed in Figure 5. While the three algorithms show similar first-order Sobol index estimates, they show some significant differences for the total index estimates. We observe that the total index estimates for the uncertain parameters $\gamma_A, \gamma_I, \gamma_H, p_{E,A}, p_C, p_{D|C}$ and $p_{H,D}$ are significantly higher for Modified Next Reaction Method, indicating that each of those parameters interacts more with the variable $Z$. Then, using Algorithms 2 or 3, it is possible to decompose $Z$ into components that correspond to the different types of transition. On Figure 6, we plotted first-order and total Sobol' index estimates associated with each of those components, for both algorithms. While there were no difference between the three algorithms for total index estimates associated with the intrinsic noise $Z$ as a whole (see the bottom of Figure 5), the analysis by type of transition reveals that the total sensitivity estimates of its components significantly differ from one algorithm to the other (see the bottom of Figure 6).
<figure><img src="image_6.png" /><figcaption>Figure 4: Boxplots of 50 independent estimates of the (top) first-order and (bottom) total Sobol' indices for $Y_{ext}^\theta$, each computed with two independent designs of $n = 2000$ input-output samples, for each simulation algorithm: (red, left) Gillespie Direct Method, (green, middle) Gillespie First Reaction, (blue, right) Modified Next Reaction.</figcaption></figure>

In accordance with the results stated in Section 3.2, we observe on the top of Figure 4 that there are no significant differences between the three algorithms for the first-order Sobol' index estimates. The only input parameters with a significant first-order effect are $\beta$ and $\gamma_E$. The sum of the first-order index estimates is far below 1 which means that interactions are not negligible. We observe on the bottom of Figure 4 that almost all inputs have a total effect significantly greater than zero. The interaction strength varies from one simulation algorithm to the other. This is due to the fact that the modeling of intrinsic randomness depends on each simulation algorithm. In particular, we observe that the total Sobol' index estimates corresponding to the Modified Next Reaction algorithm are never less
<figure><img src="image_7.png" /><figcaption>Figure 5: Aggregated (top) first-order and (bottom) total Sobol' index estimates for $Y_I^\theta$. For each input parameter (x-axis), boxplots are displayed for each simulation algorithm: (red/left) Gillespie Direct Method, (green,middle) Gillespie First Reaction, (blue,right) Modified Next Reaction. Each boxplot represents 50 independent index estimates, each of them computed with two independent designs of $n = 2000$ input-output samples.</figcaption></figure>
than their counterpart computed from the First Reaction Method algorithm. This reflects a stronger interaction with intrinsic random noise for Modified Next Reaction algorithm. Finally, as expected from the theoretical results in Section 3.2, the total index estimates associated with intrinsic randomness do not depend on the chosen simulation algorithm.

### 5.3.2 Sensitivity analysis results for $Y_I^\theta$

Since $Y_I^\theta$ is a dynamical process, we can consider the sensitivity of the whole trajectory or the sensitivity time by time. In the numerical experiments, $Y_I^\theta$ is discretized over a regular grid of size 1000 of the interval $[0, t_{end}]$. The sensitivity of the whole trajectory consists of computing estimates of the aggregated sensitivity indices of Section 2.2. These provide a scalar summary for the dynamical evolution of first-order and total Sobol' indices. They are displayed in Figure 5. While the three algorithms show similar first-order Sobol index estimates, they show some significant differences for the total index estimates. We observe that the total index estimates for the uncertain parameters $\gamma_A, \gamma_I, \gamma_H, p_{E,A}, p_C, p_{D|C}$ and $p_{H,D}$ are significantly higher for Modified Next Reaction Method, indicating that each of those parameters interacts more with the variable $Z$. Then, using Algorithms 2 or 3, it is possible to decompose $Z$ into components that correspond to the different types of transition. On Figure 6, we plotted first-order and total Sobol' index estimates associated with each of those components, for both algorithms. While there were no difference between the three algorithms for total index estimates associated with the intrinsic noise $Z$ as a whole (see the bottom of Figure 5), the analysis by type of transition reveals that the total sensitivity estimates of its components significantly differ from one algorithm to the other (see the bottom of Figure 6).
<figure><img src="image_8.png" /><figcaption>Figure 6: Aggregated (top) first-order and (bottom) total (b) Sobol' indices associated to each component of Z (each type of transition described in Table 1). For each type of transition (x-axis), boxplots are displayed for (red/left) Gillespie First Reaction and (blue right) Modified Next Reaction. Each boxplot represents 50 independent index estimates, each of them computed with two independent designs of n = 2000 input-output samples.</figcaption></figure>

The mean dynamical evolution of first-order and total Sobol' index estimates is displayed on Figure 7. Each mean is computed from 50 independent repetitions. At the beginning of the epidemic, the number of infected individuals is mostly sensitive to Z—that is, to random fluctuations inherent to the model. This confirms that intrinsic randomness rules the dynamics in the emergence phase of an epidemic disease. While the epidemic evolves, the main effect of Z quickly drops and some uncertain parameters—namely, $β$, $γ_E$ and to a less extent $p_{EA}$—gain more influence. The uncertain parameter $β$, in particular, becomes much more important than any other input and remains so until the end. Notice that, except $β$, the main
<figure><img src="image_7.png" /><figcaption>Figure 5: Aggregated (top) first-order and (bottom) total Sobol' index estimates for $Y_I^\theta$. For each input parameter (x-axis), boxplots are displayed for each simulation algorithm: (red/left) Gillespie Direct Method, (green,middle) Gillespie First Reaction, (blue,right) Modified Next Reaction. Each boxplot represents 50 independent index estimates, each of them computed with two independent designs of $n = 2000$ input-output samples.</figcaption></figure>
effect of every input (both the uncertain parameters and intrinsic noise) approaches
zero as the epidemic goes to its end, while the opposite is true for total effects. This
indicates that interactions become more prevalent near the end of the epidemic.

Although the most salient features of the performed sensitivity analyses are
shared between the three algorithms, we do observe some differences in the mean
dynamics across the three algorithms. These differences seem to be significant:
see Figure 9, where the sampling variability of the dynamics of first-order and total
Sobol' index estimates associated to p<sub>EA</sub> are displayed with functional boxplots,
namely highest density region (HDR) boxplots, obtained by using the R package
rainbow developed by Hyndman and Shang [2010]. The HDR boxplot is a vizual-
ization tool for functional data based on kernel density estimation of the scores as-
sociated to the two first principal components of the functional data (see Hyndman
[1996] for further details). The picture clearly indicates that the differences in the
mean dynamics obtained from the three different algorithms cannot be attributed to
sample variability alone. As another example, a zoom in the time t = 60 (see Fig-
ure 8) shows significant differences for the total index estimates of the parameters
γ<sub>I</sub>, γ<sub>H</sub>, p<sub>C</sub>, p<sub>DC</sub> and p<sub>HD</sub>, and to a less extent γ<sub>1</sub> and p<sub>EA</sub>.
<figure><img src="image_8.png" /><figcaption>Figure 6: Aggregated (top) first-order and (bottom) total (b) Sobol' indices associated to each component of Z (each type of transition described in Table 1). For each type of transition (x-axis), boxplots are displayed for (red/left) Gillespie First Reaction and (blue/right) Modified Next Reaction. Each boxplot represents 50 independent index estimates, each of them computed with two independent designs of n = 2000 input-output samples.</figcaption></figure>

The mean dynamical evolution of first-order and total Sobol' index estimates is displayed on Figure 7. Each mean is computed from 50 independent repetitions. At the beginning of the epidemic, the number of infected individuals is mostly sensitive to Z—that is, to random fluctuations inherent to the model. This confirms that intrinsic randomness rules the dynamics in the emergence phase of an epidemic disease. While the epidemic evolves, the main effect of Z quickly drops and some uncertain parameters—namely, $β$, $γ_E$ and to a less extent $p_{EA}$—gain more influence. The uncertain parameter $β$, in particular, becomes much more important than any other input and remains so until the end. Notice that, except $β$, the main
<figure><img src="image_9.png" /><figcaption>Figure 7: Mean dynamical evolution of (left) first-order and (right) total Sobol' indices for $Y_I$ with respect to (Figures (a) and (b)) Gillespie Direct Method, (Figure (c) and (d)) Gillespie First Reaction, (Figures (e) and (f)) Modified Next Reaction algorithm. The mean is computed from 50 independent repetitions of the estimation procedure performed with two independent designs of $n = 2000$ input-output samples.</figcaption></figure>
effect of every input (both the uncertain parameters and intrinsic noise) approaches
zero as the epidemic goes to its end, while the opposite is true for total effects. This
indicates that interactions become more prevalent near the end of the epidemic.

Although the most salient features of the performed sensitivity analyses are
shared between the three algorithms, we do observe some differences in the mean
dynamics across the three algorithms. These differences seem to be significant:
see Figure 9, where the sampling variability of the dynamics of first-order and total
Sobol' index estimates associated to p<sub>EA</sub> are displayed with functional boxplots,
namely highest density region (HDR) boxplots, obtained by using the R package
rainbow developed by Hyndman and Shang [2010]. The HDR boxplot is a vizual-
ization tool for functional data based on kernel density estimation of the scores as-
sociated to the two first principal components of the functional data (see Hyndman
[1996] for further details). The picture clearly indicates that the differences in the
mean dynamics obtained from the three different algorithms cannot be attributed to
sample variability alone. As another example, a zoom in the time t = 60 (see Fig-
ure 8) shows significant differences for the total index estimates of the parameters
γ<sub>I</sub>, γ<sub>H</sub>, p<sub>C</sub>, p<sub>DC</sub> and p<sub>HD</sub>, and to a less extent γ<sub>1</sub> and p<sub>EA</sub>.
<figure><img src="image_10.png" /><figcaption>Figure 8: Zoom on total Sobol' indices at time point $t = 60$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 9: Functional HDR boxplots of dynamical total indices associated to parameter $p_{E,A}$, obtained from 50 independent repetitions of the estimation procedure. (left) Gillespie Direct Method, (middle) Gillespie First Reaction Method, (right) Modified Next Reaction Method. The 50% HDR is plotted in dark gray, the 100% HDR in light gray and the modal curve, that is the curve in the sample with the highest density is represented by a black solid line.</figcaption></figure>

## 5.4 Some thoughts about the choice of representations

The numerical experiments confirm that the sensitivity analysis results depend on the choice of the simulation algorithm. An interesting conclusion is that Gillespie algorithms are less prone to interactions between uncertain parameters and intrinsic randomness. It implies that simulations with Gillespie algorithms are more robust to a local perturbation of uncertain input parameters as we can see below by
<figure><img src="image_9.png" /><figcaption>Figure 7: Mean dynamical evolution of (left) first-order and (right) total Sobol' indices for $Y_I$ with respect to (Figures (a) and (b)) Gillespie Direct Method, (Figure (c) and (d)) Gillespie First Reaction, (Figures (e) and (f)) Modified Next Reaction algorithm. The mean is computed from 50 independent repetitions of the estimation procedure performed with two independent designs of $n = 2000$ input-output samples.</figcaption></figure>
perturbing parameter $\beta$.

To plot Figure 10, we first simulate, for each simulation algorithm (Gillespie Direct Method, Gillespie First Reaction, Modified Next Reaction), 2000 trajectories (corresponding to 2000 different seeds) of the difference between the number of symptomatic infectious individuals computed with all uncertain parameters fixed to their nominal value and the number of symptomatic infectious individuals computed by perturbing only parameter $\beta$ by 5% from its nominal value. Then in Figure 10 are plotted highest density region (HDR) boxplots.

<figure><img src="image_12.png" /></figure>

Figure 10: Functional HDR boxplots of differences of the dynamical number of symptomatic infectious individuals computed with uncertain parameters fixed to their nominal value and by perturbing parameter $\beta$ by 5% from its nominal value: (left) Gillespie Direct Method, (middle) Gillespie First Reaction Method, (right) Modified Next Reaction Method. The 50% HDR is plotted in dark gray, the 100% HDR in light gray and the modal curve, that is the curve in the sample with the highest density is represented by a black solid line. Functional HDR boxplots were drawn from 2000 independent realizations.

# 6 Conclusion

In this work, we proposed a methodology for global sensitivity analysis of stochastic compartmental models described by continuous-time Markov chains. This approach consists in leveraging exact simulation algorithms of the continuous-time Markov chain to propose a representation of the stochastic compartmental model as a deterministic function of the uncertain parameters and controlled latent variables modeling intrinsic randomness. Our study takehome message is that sensitivity analysis results depend on the choice of the representation, which, therefore, must be made with caution and guided by the practitioners’ expertise. We exhibited three different representations for CTMC compartmental models, based respectively on Gillespie Direct Algorithm, Gillespie First Reaction Method and Modified Next
<figure><img src="image_10.png" /><figcaption>Figure 8: Zoom on total Sobol' indices at time point $t = 60$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 9: Functional HDR boxplots of dynamical total indices associated to parameter $p_{E,A}$, obtained from 50 independent repetitions of the estimation procedure. (left) Gillespie Direct Method, (middle) Gillespie First Reaction Method, (right) Modified Next Reaction Method. The 50% HDR is plotted in dark gray, the 100% HDR in light gray and the modal curve, that is the curve in the sample with the highest density is represented by a black solid line.</figcaption></figure>

## 5.4 Some thoughts about the choice of representations

The numerical experiments confirm that the sensitivity analysis results depend on the choice of the simulation algorithm. An interesting conclusion is that Gillespie algorithms are less prone to interactions between uncertain parameters and intrinsic randomness. It implies that simulations with Gillespie algorithms are more robust to a local perturbation of uncertain input parameters as we can see below by
Reaction Method. Gillespie Direct Algorithm is computationally advantageous but cannot provide a sensitivity analysis by type of transition, contrarily to Gillespie First Reaction and Modified Next Reaction algorithms. We found that Gillespie Direct algorithm and Gillespie First Reaction algorithm are more robust to local perturbations of the uncertain parameters. We applied our approach to a stochastic compartmental model of SARS-CoV-2 spread.

In the present paper, we considered Markovian models only. However an in-
teresting follow-up would be to extend our results to non-Markovian stochastic
processes by using Sellke's construction (Sellke [1983]).

Declaration

Declarations of interest: none.

References

D. F. Anderson. A modified next reaction method for simulating chemical systems with time dependent propensities and delays. *The Journal of Chemical Physics*, 127(21):214107, 2007. doi: 10.1063/1.2799998.

P. Bittihn and R. Golestanian. Stochastic effects on the dynamics of an epidemic due to population subdivision. *Chaos: An Interdisciplinary Journal of Nonlinear Science*, 30(10):101–102, Oct. 2020. ISSN 1054-1500. doi: 10.1063/5.0028972.

F. Brauer. *Compartmental Models in Epidemiology*, pages 19–79. Springer Berlin Heidelberg, Berlin, Heidelberg, 2008. ISBN 978-3-540-78911-6. doi: 10.1007/978-3-540-78911-6_2.

T. Britton. Stochastic epidemic models: a survey. *Mathematical biosciences*, 225 1:24–35, 2009.

B. Cazelles, C. Champagne, B. Nguyen-Van-Yen, C. Comiskey, E. Vergu, and B. Roche. A mechanistic and data-driven reconstruction of the time-varying reproduction number: Application to the covid-19 epidemic. *PLOS Computational Biology*, 17(7):1–20, 07 2021. doi: 10.1371/journal.pcbi.1009211.

A. Courcoul, H. Monod, M. Nielen, D. Klinkenberg, L. Hogerwerf, F. Beaudeau, and E. Vergu. Modelling the effect of heterogeneity of shedding on the within herd coxiella burnetii spread and identification of key parameters by sensitivity analysis. *Journal of Theoretical Biology*, 284(1):130–141, 2011. ISSN 0022-5193. doi: https://doi.org/10.1016/j.jtbi.2011.06.017.

L. Cristancho Fajardo, P. Ezanno, and E. Vergu. Accounting for farmers’ control decisions in a model of pathogen spread through animal trade. *Scientific Reports*, 11(1):9581, May 2021. ISSN 2045-2322. doi: 10.1038/s41598-021-88471-6.
perturbing parameter $\beta$.

To plot Figure 10, we first simulate, for each simulation algorithm (Gillespie Direct Method, Gillespie First Reaction, Modified Next Reaction), 2000 trajectories (corresponding to 2000 different seeds) of the difference between the number of symptomatic infectious individuals computed with all uncertain parameters fixed to their nominal value and the number of symptomatic infectious individuals computed by perturbing only parameter $\beta$ by 5% from its nominal value. Then in Figure 10 are plotted highest density region (HDR) boxplots. It is clear on these plots that the small perturbation applied to parameter $\beta$ has a much stronger impact when Modified Next Reaction Method is used for simulations. We thus expect that quantities calculated from a Monte-Carlo sampling scheme are less robust to an inaccurate estimate of parameter $\beta$ when using simulations based on Modified Next Reaction Algorithm.

<figure><img src="image_12.png" /><figcaption>Figure 10: Functional HDR boxplots of differences of the dynamical number of symptomatic infectious individuals computed with uncertain parameters fixed to their nominal value and by perturbing parameter $\beta$ by 5% from its nominal value: (left) Gillespie Direct Method, (middle) Gillespie First Reaction Method, (right) Modified Next Reaction Method. The 50% HDR is plotted in dark gray, the 100% HDR in light gray and the modal curve, that is the curve in the sample with the highest density is represented by a black solid line. Functional HDR boxplots are plotted from 2000 independent realizations.</figcaption></figure>

# 6 Conclusion

In this work, we leveraged three different exact simulation algorithms for continuous-time Markov chains from the state of the art which we combined with common tools from variance-based sensitivity analysis to perform a global sensitivity analysis of stochastic compartmental models. Such a methodology was introduced by Navarro Jimenez et al. [2016] in the framework of chemical reaction networks, and
S. Da Veiga, F. Gamboa, B. Iooss, and C. Prieur. *Basics and trends in sensitivity analysis: Theory and practice in R*. SIAM, 2021.

D. Dupuy, C. Helbert, and J. Franco. DiceDesign and DiceEval: Two R packages for design and analysis of computer experiments. *Journal of Statistical Software*, 65(11):1–38, 2015.

S. N. Ethier and T. G. Kurtz. *Markov processes – characterization and convergence*, chapter 4, 6. Wiley Series in Probability and Mathematical Statistics: Probability and Mathematical Statistics. John Wiley & Sons Inc., New York, 1986. ISBN 0-471-08186-8.

P. Étoré, C. Prieur, D. K. Pham, and L. Li. Global sensitivity analysis for models described by stochastic differential equations. *Methodology and Computing in Applied Probability*, 22(2):803–831, June 2020. ISSN 1387-5841, 1573-7713. doi: 10.1007/s11009-019-09732-6.

J.-C. Fort, T. Klein, and A. Lagnoux. Global sensitivity analysis and Wasserstein spaces. *SIAM/ASA Journal on Uncertainty Quantification*, 9(2):880–921, 2021. doi: 10.1137/20M1354957.

F. Gamboa, A. Janon, T. Klein, and A. Lagnoux. Sensitivity analysis for multidimensional and functional outputs. *Electronic Journal of Statistics*, 8(1):575–603, 2014. doi: 10.1214/14-EJS895.

D. T. Gillespie. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. *Journal of Computational Physics*, 22(4):403–434, 1976. ISSN 0021-9991. doi: https://doi.org/10.1016/0021-9991(76)90041-3.

S. Goel, S. K. Bhatia, J. P. Tripathi, S. Bugalia, M. Rana, and V. P. Bajiya. SIRC epidemic model with cross-immunity and multiple time delays. *Journal of Mathematical Biology*, 87(3):42, 2023.

K. Hanthanan Arachchilage, M. Y. Hussaini, N. Cogan, and M. H. Cortez. Exploring how ecological and epidemiological processes shape multi-host disease dynamics using global sensitivity analysis. *Journal of Mathematical Biology*, 86 (5):83, 2023.

J. L. Hart, A. Alexanderian, and P. A. Gremaud. Efficient computation of Sobol' indices for stochastic models. *SIAM Journal on Scientific Computing*, 39(4): A1514–A1530, 2017. doi: 10.1137/16M106193X.

W. Hoeffding. A Class of Statistics with Asymptotically Normal Distribution. *The Annals of Mathematical Statistics*, 19(3):293 – 325, 1948. doi: 10.1214/aoms/1177730196.
using simulations from Modified Next Reaction Algorithm. In this paper, we discussed for the first time the impact of the choice of the algorithm used for model simulations on the result of global sensitivity analysis. We implemented and compared three sensitivity analyses based on simulations obtained from different exact simulation algorithms of a SARS-CoV-2 epidemic model. We observed that the different simulation algorithms are not equivalent in terms of robustness with respect to the uncertainty on epidemic parameters. Indeed with Figure 10 we exhibited that variations in the value of parameter $β$ have a stronger influence on the variability of the simulated number of infectious individuals by using simulations produced with Modified Next Reaction algorithm.

In the present paper, we considered Markovian models only. However an interesting follow-up would be to extend our results to non-Markovian stochastic processes by using Sellke's construction (Sellke [1983]).

## Acknowledgment

We thank the Associate Editor and two anonymous reviewers for their thorough reading and constructive feedback that led to an improved version of this manuscript.

## Declaration

Declarations of interest: none.

## References

D. F. Anderson. A modified next reaction method for simulating chemical systems with time dependent propensities and delays. *The Journal of Chemical Physics*, 127(21):214107, 2007. doi: 10.1063/1.2799998.

P. Bittihn and R. Golestanian. Stochastic effects on the dynamics of an epidemic due to population subdivision. *Chaos: An Interdisciplinary Journal of Nonlinear Science*, 30(10):101–102, Oct. 2020. ISSN 1054-1500. doi: 10.1063/5.0028972.

F. Brauer. *Compartmental Models in Epidemiology*, pages 19–79. Springer Berlin Heidelberg, Berlin, Heidelberg, 2008. ISBN 978-3-540-78911-6. doi: 10.1007/978-3-540-78911-6_2.

T. Britton. Stochastic epidemic models: a survey. *Mathematical biosciences*, 225:1–24–35, 2009.

B. Cazelles, C. Champagne, B. Nguyen-Van-Yen, C. Comiskey, E. Vergu, and B. Roche. A mechanistic and data-driven reconstruction of the time-varying reproduction number: Application to the covid-19 epidemic. *PLOS Computational Biology*, 17(7):1–20, 07 2021. doi: 10.1371/journal.pcbi.1009211.
T. Homma and A. Saltelli. Importance measures in global sensitivity analysis of nonlinear models. *Reliability Engineering and System Safety*, 52(1):1–17, 1996. ISSN 0951-8320. doi: https://doi.org/10.1016/0951-8320(96)00002-6.

R. J. Hyndman. Computing and graphing highest density regions. *The American Statistician*, 50(2):120–126, 1996.

R. J. Hyndman and H. L. Shang. Rainbow plots, bagplots, and boxplots for functional data. *Journal of Computational and Graphical Statistics*, 19(1):29–45, 2010.

B. Iooss, S. Da Veiga, A. Janon, and G. Pujol. sensitivity: *Global Sensitivity Analysis of Model Outputs*, 2021. URL https://CRAN.R-project.org/package=sensitivity. R package version 1.24.0.

M. N. Jimenez, O. P. Le Maître, and O. M. Knio. Nonintrusive polynomial chaos expansions for sensitivity analysis in stochastic differential equations. *SIAM/ASA Journal on Uncertainty Quantification*, 5(1):378–402, Jan. 2017. ISSN 2166-2525. doi: 10.1137/16M1061989.

S. Karlin and H. M. Taylor. *A second course on stochastic processes*. Academic Press, 1981.

E. S. Knock, L. K. Whittles, J. A. Lees, P. N. Perez-Guzman, R. Verity, R. G. FitzJohn, K. A. M. Gaythorpe, N. Imai, W. Hinsley, L. C. Okell, A. Rosello, N. Kantas, C. E. Walters, S. Bhatia, O. J. Watson, C. Whittaker, L. Cattarino, A. Boonyasiri, B. A. Djaafara, K. Fraser, H. Fu, H. Wang, X. Xi, C. A. Donnelly, E. Jauneikaite, D. J. Laydon, P. J. White, A. C. Ghani, N. M. Ferguson, A. Cori, and M. Baguelin. Key epidemiological drivers and impact of interventions in the 2020 sars-cov-2 epidemic in england. *Science Translational Medicine*, 13(602): 42–62, 2021. doi: 10.1126/scitranslmed.abg4262.

T. G. Kurtz. Representation and approximation of counting processes. In W. H. Fleming and L. G. Gorostiza, editors, *Advances in Filtering and Optimal Stochastic Control*, pages 177–191, Berlin, Heidelberg, 1982. Springer Berlin Heidelberg. ISBN 978-3-540-39517-1.

M. Lamboni, M. Hervé, and M. David. Multivariate sensitivity analysis to measure global contribution of input factors in dynamic models. *Reliability Engineering and System Safety*, 96(4):450–459, 2011. ISSN 0951-8320. doi: https://doi.org/10.1016/j.ress.2010.12.002.

O. Le Maître and O. Knio. PC analysis of stochastic differential equations driven by wiener noise. *Reliability Engineering & System Safety*, 135:107–124, 2015. ISSN 0951-8320. doi: https://doi.org/10.1016/j.ress.2014.11.002.
A. Courcoul, H. Monod, M. Nielen, D. Klinkenberg, L. Hogerwerf, F. Beauveau, and E. Vergu. Modelling the effect of heterogeneity of shedding on the within herd coxiella burnetii spread and identification of key parameters by sensitivity analysis. *Journal of Theoretical Biology*, 284(1):130-141, 2011. ISSN 0022-5193. doi: https://doi.org/10.1016/j.jtbi.2011.06.017.

L. Cristancho Fajardo, P. Ezanno, and E. Vergu. Accounting for farmers’ control decisions in a model of pathogen spread through animal trade. *Scientific Reports*, 11(1):9581, May 2021. ISSN 2045-2322. doi: 10.1038/s41598-021-88471-6.

S. Da Veiga, F. Gamboa, B. Iooss, and C. Prieur. *Basics and trends in sensitivity analysis: Theory and practice in R*. SIAM, 2021.

D. Dupuy, C. Helbert, and J. Franco. DiceDesign and DiceEval: Two R packages for design and analysis of computer experiments. *Journal of Statistical Software*, 65(11):1-38, 2015.

S. N. Ethier and T. G. Kurtz. *Markov processes – characterization and convergence*, chapter 4, 6. Wiley Series in Probability and Mathematical Statistics: Probability and Mathematical Statistics. John Wiley & Sons Inc., New York, 1986. ISBN 0-471-08186-8.

P. Étoré, C. Prieur, D. K. Pham, and L. Li. Global sensitivity analysis for models described by stochastic differential equations. *Methodology and Computing in Applied Probability*, 22(2):803-831, June 2020. ISSN 1387-5841, 1573-7713. doi: 10.1007/s11009-019-09732-6.

J.-C. Fort, T. Klein, and A. Lagnoux. Global sensitivity analysis and Wasserstein spaces. *SIAM/ASA Journal on Uncertainty Quantification*, 9(2):880-921, 2021. doi: 10.1137/20M1354957.

F. Gamboa, A. Janon, T. Klein, and A. Lagnoux. Sensitivity analysis for multidimensional and functional outputs. *Electronic Journal of Statistics*, 8(1):575-603, 2014. doi: 10.1214/14-EJS895.

D. T. Gillespie. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. *Journal of Computational Physics*, 22(4):403-434, 1976. ISSN 0021-9991. doi: https://doi.org/10.1016/0021-9991(76)90041-3.

S. Goel, S. K. Bhatia, J. P. Tripathi, S. Bugalia, M. Rana, and V. P. Bajiya. SIRC epidemic model with cross-immunity and multiple time delays. *Journal of Mathematical Biology*, 87(3):42, 2023.

K. Hanthanan Arachchilage, M. Y. Hussaini, N. Cogan, and M. H. Cortez. Exploring how ecological and epidemiological processes shape multi-host disease dynamics using global sensitivity analysis. *Journal of Mathematical Biology*, 86 (5):83, 2023.
O. P. Le Maître, O. M. Knio, and A. Moraes. Variance decomposition in stochastic simulators. *The Journal of Chemical Physics*, 142(24):244115, 2015. doi: 10.1063/1.4922922.

C. Lin and B. Tang. Latin hypercubes and space-filling designs. In A. Dean, M. Morris, J. Stufken, and D. Bingham, editors, *Handbook of design and analysis of experiments*, Handbooks of modern statistical methods. Chapman & Hall/CRC, 2015.

S. Marino, I. B. Hogue, C. J. Ray, and D. E. Kirschner. A methodology for performing global uncertainty and sensitivity analysis in systems biology. *Journal of Theoretical Biology*, 254(1):178–196, 2008. ISSN 0022-5193. doi: https://doi.org/10.1016/j.jtbi.2008.04.011. URL https://www.sciencedirect.com/science/article/pii/S0022519308001896.

A. Marrel, B. Iooss, S. Da Veiga, and M. Ribatet. Global sensitivity analysis of stochastic computer models with joint metamodels. *Statistics and Computing*, 22(3):833–847, May 2012. ISSN 1573-1375. doi: 10.1007/s11222-011-9274-8.

M. Massard, R. Eftimie, A. Perasso, and B. Saussereau. A multi-strain epidemic model for covid-19 with infected and asymptomatic cases: Application to french data. *Journal of Theoretical Biology*, 545:111117, 2022. ISSN 0022-5193. doi: https://doi.org/10.1016/j.jtbi.2022.111117. URL https://www.sciencedirect.com/science/article/pii/S0022519322001151.

G. Mazo. A trade-off between explorations and repetitions for estimators of two global sensitivity indices in stochastic models induced by probability measures. *SIAM/ASA Journal on Uncertainty Quantification*, 9(4):1673–1713, 2021. doi: 10.1137/19M1272706.

M. Navarro Jimenez, O. P. Le Maître, and O. M. Knio. Global sensitivity analysis in stochastic simulators of uncertain reaction networks. *The Journal of Chemical Physics*, 145(24):244106, 2016. doi: 10.1063/1.4971797.

C. N. Ngonghala, M. I. Teboh-Ewungkem, and G. A. Ngwa. Persistent oscillations and backward bifurcation in a malaria model with varying human and mosquito populations: implications for control. *Journal of mathematical biology*, 70(7): 1581–1622, 2015.

R Core Team. *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing, Vienna, Austria, 2021. URL https://www.R-project.org/.

Q. Richard, S. Alizon, M. Choisy, M. T. Sofonea, and R. Djidjou-Demasse. Age-structured non-pharmaceutical interventions for optimal control of covid-19 epidemic. *PLOS Computational Biology*, 17(3):1–25, 03 2021. doi: 10.1371/journal.pcbi.1008776.
J. L. Hart, A. Alexanderian, and P. A. Gremaud. Efficient computation of Sobol' indices for stochastic models. *SIAM Journal on Scientific Computing*, 39(4): A1514–A1530, 2017. doi: 10.1137/16M106193X.

T. Homma and A. Saltelli. Importance measures in global sensitivity analysis of nonlinear models. *Reliability Engineering and System Safety*, 52(1):1–17, 1996. ISSN 0951-8320. doi: https://doi.org/10.1016/0951-8320(96)00002-6.

R. J. Hyndman. Computing and graphing highest density regions. *The American Statistician*, 50(2):120–126, 1996.

R. J. Hyndman and H. L. Shang. Rainbow plots, bagplots, and boxplots for functional data. *Journal of Computational and Graphical Statistics*, 19(1):29–45, 2010.

B. Iooss, S. Da Veiga, A. Janon, and G. Pujol. sensitivity: *Global Sensitivity Analysis of Model Outputs*, 2021. URL https://CRAN.R-project.org/package=sensitivity. R package version 1.24.0.

M. N. Jimenez, O. P. Le Maître, and O. M. Knio. Nonintrusive polynomial chaos expansions for sensitivity analysis in stochastic differential equations. *SIAM/ASA Journal on Uncertainty Quantification*, 5(1):378–402, Jan. 2017. ISSN 2166-2525. doi: 10.1137/16M1061989.

S. Karlin and H. M. Taylor. *A second course on stochastic processes*. Academic Press, 1981.

E. S. Knock, L. K. Whittles, J. A. Lees, P. N. Perez-Guzman, R. Verity, R. G. FitzJohn, K. A. M. Gaythorpe, N. Imai, W. Hinsley, L. C. Okell, A. Rosello, N. Kantas, C. E. Walters, S. Bhatia, O. J. Watson, C. Whittaker, L. Cattarino, A. Boonyasiri, B. A. Djaafara, K. Fraser, H. Fu, H. Wang, X. Xi, C. A. Donnelly, E. Jauneikaite, D. J. Laydon, P. J. White, A. C. Ghani, N. M. Ferguson, A. Cori, and M. Baguelin. Key epidemiological drivers and impact of interventions in the 2020 sars-cov-2 epidemic in england. *Science Translational Medicine*, 13(602): 42–62, 2021. doi: 10.1126/scitranslmed.abg4262.

T. G. Kurtz. Representation and approximation of counting processes. In W. H. Fleming and L. G. Gorostiza, editors, *Advances in Filtering and Optimal Stochastic Control*, pages 177–191, Berlin, Heidelberg, 1982. Springer Berlin Heidelberg. ISBN 978-3-540-39517-1.

M. Lamboni, M. Hervé, and M. David. Multivariate sensitivity analysis to measure global contribution of input factors in dynamic models. *Reliability Engineering and System Safety*, 96(4):450–459, 2011. ISSN 0951-8320. doi: https://doi.org/10.1016/j.ress.2010.12.002.
L. Rimbaud, C. Bruchou, S. Dallot, D. R. J. Pleydell, E. Jacquot, S. Soubeyrand, and G. Thébaud. Using sensitivity analysis to identify key factors for the propagation of a plant epidemic. *Royal Society Open Science*, 5(1):171435, 2018. doi: 10.1098/rsos.171435.

A. Saltelli, K. Chan, and E. M. Scott. *Sensitivity Analysis*. John Wiley & Sons, 2000.

T. Sellke. On the asymptotic distribution of the size of a stochastic epidemic. *Journal of Applied Probability*, 20(2):390–394, 1983. doi: 10.2307/3213811.

I. M. Sobol'. Sensitivity analysis for non-linear mathematical models. *Mathematical Modelling and Computational Experiment*, 1:407–414, 1993.

A. J. Torii, R. Begnini, H. M. Kroetz, O. M. I. Matar, R. H. Lopez, and L. F. F. Miguel. Global sensitivity analysis for mathematical models comparison. *Computational and Applied Mathematics*, 42(8):345, 2023.

S. D. Veiga. Kernel-based ANOVA decomposition and shapley effects – application to global sensitivity analysis, 2021. arXiv:2101.05487.

J. Yang, Y. Chen, and F. Xu. Effect of infection age on an SIS epidemic model on complex networks. *Journal of mathematical biology*, 73:1227–1249, 2016.

X. Zhu and B. Sudret. Global sensitivity analysis for stochastic simulators based on generalized lambda surrogate models. *Reliability Engineering & System Safety*, 214:107815, 2021. ISSN 0951-8320. doi: https://doi.org/10.1016/j.ress.2021.107815.

X. Zhu and B. Sudret. Stochastic polynomial chaos expansions to emulate stochastic simulators. *International Journal for Uncertainty Quantification*, 13(2):31–52, 2023. ISSN 2152-5080.
O. Le Maître and O. Knio. PC analysis of stochastic differential equations driven by wiener noise. *Reliability Engineering & System Safety*, 135:107-124, 2015. ISSN 0951-8320. doi: https://doi.org/10.1016/j.ress.2014.11.002.

O. P. Le Maître, O. M. Knio, and A. Moraes. Variance decomposition in stochastic simulators. *The Journal of Chemical Physics*, 142(24):244115, 2015. doi: 10.1063/1.4922922.

C. Lin and B. Tang. Latin hypercubes and space-filling designs. In A. Dean, M. Morris, J. Stufken, and D. Bingham, editors, *Handbook of design and analysis of experiments*, Handbooks of modern statistical methods. Chapman & Hall/CRC, 2015.

S. Marino, I. B. Hogue, C. J. Ray, and D. E. Kirschner. A methodology for performing global uncertainty and sensitivity analysis in systems biology. *Journal of Theoretical Biology*, 254(1):178-196, 2008. ISSN 0022-5193. doi: https://doi.org/10.1016/j.jtbi.2008.04.011. URL https://www.sciencedirect.com/science/article/pii/S0022519308001896.

A. Marrel, B. Iooss, S. Da Veiga, and M. Ribatet. Global sensitivity analysis of stochastic computer models with joint metamodels. *Statistics and Computing*, 22(3):833-847, May 2012. ISSN 1573-1375. doi: 10.1007/s11222-011-9274-8.

M. Massard, R. Eftimie, A. Perasso, and B. Saussereau. A multi-strain epidemic model for covid-19 with infected and asymptomatic cases: Application to french data. *Journal of Theoretical Biology*, 545:111117, 2022. ISSN 0022-5193. doi: https://doi.org/10.1016/j.jtbi.2022.111117. URL https://www.sciencedirect.com/science/article/pii/S0022519322001151.

G. Mazo. A trade-off between explorations and repetitions for estimators of two global sensitivity indices in stochastic models induced by probability measures. *SIAM/ASA Journal on Uncertainty Quantification*, 9(4):1673-1713, 2021. doi: 10.1137/19M1272706.

M. Navarro Jimenez, O. P. Le Maître, and O. M. Knio. Global sensitivity analysis in stochastic simulators of uncertain reaction networks. *The Journal of Chemical Physics*, 145(24):244106, 2016. doi: 10.1063/1.4971797.

C. N. Ngonghala, M. I. Teboh-Ewungkem, and G. A. Ngwa. Persistent oscillations and backward bifurcation in a malaria model with varying human and mosquito populations: implications for control. *Journal of mathematical biology*, 70(7): 1581-1622, 2015.

R Core Team. *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing, Vienna, Austria, 2021. URL https://www.R-project.org/.
Q. Richard, S. Alizon, M. Choisy, M. T. Sofonea, and R. Djidjou-Demasse. Age-structured non-pharmaceutical interventions for optimal control of covid-19 epidemic. *PLOS Computational Biology*, 17(3):1–25, 03 2021. doi: 10.1371/journal.pcbi.1008776.

L. Rimbaud, C. Bruchou, S. Dallot, D. R. J. Pleydell, E. Jacquot, S. Soubeyrand, and G. Thébaud. Using sensitivity analysis to identify key factors for the propagation of a plant epidemic. *Royal Society Open Science*, 5(1):171435, 2018. doi: 10.1098/rsos.171435.

A. Saltelli, K. Chan, and E. M. Scott. *Sensitivity Analysis*. John Wiley & Sons, 2000.

T. Sellke. On the asymptotic distribution of the size of a stochastic epidemic. *Journal of Applied Probability*, 20(2):390–394, 1983. doi: 10.2307/3213811.

I. M. Sobol'. Sensitivity analysis for non-linear mathematical models. *Mathematical Modelling and Computational Experiment*, 1:407–414, 1993.

A. J. Torii, R. Begnini, H. M. Kroetz, O. M. I. Matar, R. H. Lopez, and L. F. F. Miguel. Global sensitivity analysis for mathematical models comparison. *Computational and Applied Mathematics*, 42(8):345, 2023.

S. D. Veiga. Kernel-based ANOVA decomposition and shapley effects – application to global sensitivity analysis, 2021. arXiv:2101.05487.

C. K. Williams and C. E. Rasmussen. *Gaussian processes for machine learning*, volume 2. MIT press Cambridge, MA, 2006.

J. Yang, Y. Chen, and F. Xu. Effect of infection age on an SIS epidemic model on complex networks. *Journal of mathematical biology*, 73:1227–1249, 2016.

X. Zhu and B. Sudret. Global sensitivity analysis for stochastic simulators based on generalized lambda surrogate models. *Reliability Engineering & System Safety*, 214:107815, 2021. ISSN 0951-8320. doi: https://doi.org/10.1016/j.ress.2021.107815.

X. Zhu and B. Sudret. Stochastic polynomial chaos expansions to emulate stochastic simulators. *International Journal for Uncertainty Quantification*, 13(2):31–52, 2023. ISSN 2152-5080.