# An extension of Sellke construction and uncertainty quantification for non-Markovian epidemic models

Henri Mermoz Kouye, Clémentine Prieur, Elisabeta Vergu

► To cite this version:

Henri Mermoz Kouye, Clémentine Prieur, Elisabeta Vergu. An extension of Sellke construction and uncertainty quantification for non-Markovian epidemic models. Mathematical Modelling of Natural Phenomena, 2025, pp.1-31. 10.1051/mmnp/2025027. hal-04719348v2

HAL Id: hal-04719348

https://hal.science/hal-04719348v2

Submitted on 25 Sep 2025

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
# An extension of Sellke construction and uncertainty quantification for non-Markovian epidemic models

Henri Mermoz KOUYE$^{(1,2)}$, 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

September 25, 2025

## Abstract

Several major epidemic events over the past two decades have highlighted the importance of developing and studying non-Markovian compartmental models. Selke [1983] introduced an ingenious construction for the SIR epidemic process to study the final size of epidemics. In this paper, we extend this construction to the *SEI<sub>1</sub>I<sub>2</sub>RS* model. This model is chosen for its compactness, while including parallel infectious stages (*I*<sub>1</sub> and *I*<sub>2</sub>) and cycles (aka loops) due to reinfection. Our methodology easily generalizes to a general class of stochastic compartmental models in closed populations, including SIR-like models (a series of compartments in one row), SEIAR-like models (parallel compartments), but also models with cycles. Our construction inherits from Sellke construction its ability to handle both Markovian and non-Markovian frameworks. Also, it naturally leads to a representation of the epidemic process under the form of a deterministic function of uncertain parameters (such as epidemic parameters) and variables modeling internal noise. Based on this representation, we propose a global sensitivity analysis of the *SEI<sub>1</sub>I<sub>2</sub>RS* model. With our methodology we are able to quantify epistemic uncertainty due to the lack of knowledge on epidemic parameters and statistical uncertainty induced by stochasticity of the model. Finally we provide numerical experiments in both Markovian and non-Markovian frameworks.

**Keywords** Sellke construction, compartmental models, non-Markovian epidemic process, global sensitivity analysis.

## 1 Introduction

The COVID-19 pandemic has underlined the importance of mathematical modeling of epidemics. There are a wide variety of mathematical models that can be used to predict the spread of an epidemic or to guide public health decision-making.
Models designed to capture the spread of an epidemic in a given population include stochastic compartmental models. These models are related to the seminal model introduced in Kermack and McKendrick [1927] known as the SIR model, where the population is divided into three compartments: susceptible (S), infected (I), and removed (R). The model then follows the change in the proportion of the population belonging to each compartment through time, by reproducing numerically the transition rules from one compartment to the other. Advanced models can be designed with a larger number of compartments to better reflect the characteristics of an outbreak, as well as to match data on which dynamic inference can be performed. Resulting processes can be either Markovian or non-Markovian, depending on the modeling framework. Although the Markovian setting is quite restrictive, as it is characterized by memorylessness [Nowzari et al., 2015, Großmann et al., 2021, Sofonea et al., 2021, Saeedian et al., 2017], it has been widely studied in the literature as it is simpler to analyze from a mathematical point of view. Memory in non-Markovian processes is induced by sampling non-exponentially distributed holding times between consecutive state transitions [Van Mieghem and van de Bovenkamp, 2013, Streftaris and Gibson, 2012], or by considering stochastic systems with time delay [see Brett and Galla, 2013, e.g.].

Both Markovian and non-Markovian models are subject to uncertainties. On the one hand, epistemic uncertainty due to a lack of knowledge of the parameters. On the other hand, statistical uncertainty induced by the stochasticity of the model. In compartmental epidemic models, epistemic uncertainty is characterized by imperfect knowledge of the characteristics of the pathogen and the macroscopic behavior of the population, while statistical uncertainty (or internal noise) models the specificity of individuals with respect to infection and behavior. Most risk assessment frameworks for infectious diseases only take into account epistemic uncertainty. However, internal noise has a major impact on the evolution of epidemics that needs to be quantified [Penn et al., 2023]. One way of approaching this issue is to perform a global sensitivity analysis to quantify how the uncertainty in the output of a stochastic model is related to the uncertainty in its parameters, but also to the internal noise. The approach we consider in this paper is based on a representation of the stochastic model under the form of a deterministic function of a set of epidemic parameters and a collection of variables modeling internal noise. This approach has been successfully used in the framework of Markovian models. In Le Maître et al. [2015], the authors used the random time-change representation introduced in Ethier and Kurtz [1986], while the authors in Kouye et al. [2024] used Gillespie algorithms [Gillespie, 1976].

The objective of our work is twofold. We first aim at extending the construction introduced in Sellke [1983] and discussed hereafter to a general class of compartmental models in closed population. We then aim at performing a global sensitivity analysis based on this extended construction, which can handle both Markovian and non-Markovian dynamics.

Sellke [1983] introduced an ingenious construction used to sample the final size distribution of a SIR epidemic in a finite population. It relies on an individual-
based approach and it has the advantage to enable the simulation of a SIR epidemic
in both Markovian and non-Markovian frameworks. Since then, it has been studied
and extended in many ways. In Reinert [1995], non independent tolerance thresh-
olds of susceptible individuals to infection and non independent infectious periods
were considered. Andersson and Britton [2000] studied SIR-multitype epidemic
models. In Streftaris and Gibson [2012], the unit-rate exponential distribution of
thresholds was changed to a Weibull distribution. House [2014] generalized the
construction to the case of finite heterogeneous population, while Di Lauro et al.
[2022] focused on the age-dependent SIR model. An extension to the classical
SEIR (Susceptible-Exposed-Infectious-Recovered) model was introduced in Brit-
ton and Pardoux [2019b]. This extension is straightforward because, on the one
hand, the jumps in compartment *E* follow the same transition mechanism as those
in compartment *I* and, on the other hand, the individuals in compartment *E* are not
infectious. However, in view of all these works, it appears that so far, the general-
izations of Sellke construction mainly focused on SIR-like model. For instance, to
the best of our knowledge, none of these works consider models with reinfection
and more generally compartmental models that include parallel infectious stages
(like SEIAR: Susceptible-Exposed-Infectious-Asymptomatic-Recovered) or such
that their corresponding graph includes cycles (or loops) like the classical SIS
(Susceptible-Infectious-Susceptible).

In our work, we extend Sellke construction to a more general class of compart-
mental models in closed populations. Throughout the paper, and for sake of clarity
in the presentation, we focus on the *SEI*₁*I*₂*RS* model. However, the methodol-
ogy easily generalizes to a general class of stochastic compartmental models in
closed population, including SIR-like models (a series of compartments in one
row), SEIAR-like models (the existence of compartments in parallel), but also
models with cycles. The *SEI*₁*I*₂*RS* model was chosen for its compactness, while
including multiple parallel infectious stages (*I*₁ and *I*₂) and cycles due to reinfec-
tion. The construction we propose inherits from the original Sellke construction
the ability to cover both Markovian and non-Markovian framework, depending on
the choice of probability distributions for holding times between consecutive state
transitions. One of the advantages of the construction we propose, and also of
the one of Sellke, is that it permits a clear identification of internal noise through
the behavior and characteristics of individuals. Using the extended construction,
we identify a representation of the epidemic process in the form of a determin-
istic function of uncertain epidemic parameters and variables modelling internal
noise. Based on this representation, we perform global sensitivity analysis of the
*SEI*₁*I*₂*RS* model both in Markovian and non-Markovian frameworks.

The paper is organized as follows. In Section 2 we review Sellke original con-
struction (Section 2.1) from which we identify a representation algorithm (Section
2.2) for the SIR model. Then in Section 3 we generalize Sellke construction to
a model which includes both loops due to reinfection and two parallel infectious
stages, the *SEI1I2RS* model from which we derive a representation algorithm.
This algorithm is used in Section 4 for simulations. Section 5 presents the method-
ology and the numerical experiments of global sensitivity analysis based on the representation algorithm. In particular we identify key factors (including epidemic parameters but also variables modeling intrinsic randomness) or key parameter interactions in both Markovian and non-Markovian frameworks. Finally the paper finishes with concluding remarks (Section 6).

## 2 Sellke original construction, application to the identification of a deterministic representation of stochastic (non)-Markovian SIR model

In Section 2.1 we recall the original Sellke construction from which we identify in Section 2.2 a deterministic representation for the non-Markovian stochastic SIR model.

### 2.1 Sellke construction

Consider the Susceptible-Infected-Recovered (SIR) model (see Figure 1) and assume that the population is closed and comprises $N$ individuals. We also assume that at each time $t \ge 0$, at most one transition can occur. In the following, we denote by $\{W^\theta(t) = (W_S^\theta(t), W_I^\theta(t), W_R^\theta(t)), t \ge 0\}$ the jump process counting at each time $t \ge 0$ the number of individuals in each compartment, with $\theta = (\beta, \gamma_I)$, $\beta$ the transmission rate and $1/\gamma_I$ the mean sojourn time in $I$.

<figure><img src="image_2.png" /><figcaption>Figure 1: The classical SIR model.</figcaption></figure>

Sellke [1983] proposed an individual-based construction for the SIR model. At the start of the epidemic ($t=0$), each individual of the population of size $N$ is labeled ($i=1, \dots, N$) and depending on his initial health status, he is given a set of variables that characterize his behavior towards infection. More precisely, if the individual labeled $i$ is initially susceptible, he is given a threshold $Q_i > 0$ and a sojourn duration $L_i > 0$ in $I$. If this individual is initially infectious, then he is only given a sojourn duration $L_i$ in compartment $I$.

Two types of events are responsible for the change of health status: infection or recovery. Occurrence of infection events is linked to a function called the infectious pressure defined as $P(t) = \frac{\beta}{N} \int_0^t W_I^\theta(s) ds$ so that if the individual labeled $i$ is initially susceptible, he remains susceptible as long as $Q_i > P(t)$ for $t \ge 0$. This function accounts for the pressure exerted by the mass of infectious on susceptible individuals. At time $\tau_i^S = \inf\{s \ge 0 \mid P(s) \ge Q_i\}$, the individual labeled $i$ gets
infected and moves to compartment *I* (see Figure 2). Individuals with thresholds which remain greater than the value of the infectious pressure escape infection. In Figure 1 below, the labels of individuals have been chosen so that the sequence of tolerance thresholds ($Q_i$)$_{1\le i\le N}$ is increasing.

<figure><img src="image_3.png" /><figcaption>Figure 2: An illustration of the evolution of infectious pressure $t \mapsto P(t)$. Marks ◆ denote infection events. $0 < Q_1 < Q_2 < Q_3 < Q_4$ are tolerance thresholds. $\tau_1^S, \dots, \tau_4^S$ denote transition times from susceptible to infected associated with thresholds $Q_1, \dots, Q_4$.</figcaption></figure>

**Remark 1.** If the tolerance thresholds $Q_i$, $i = 1, \dots, N$ and the sojourn durations $L_i$, $i = 1, \dots, N$ are independent from each other and independently distributed under the standard exponential distribution, then the jump process $\{W^\theta(t), t \ge 0\}$ is a continuous-time Markov chain. Moreover, for any time $t \ge 0$, the infection rate is given by $\frac{\beta}{N} W_S^\theta(t) W_I^\theta(t)$ and the recovery rate is given by $\gamma_I W_I^\theta(t)$. This is a classical result which can be found, e.g., in Britton and Pardoux [2019a]. We briefly recall the main lines for the computation of the infection and recovery rates. Let us first focus on the computation of the infection rate. For any $t \ge 0$, the probability that an infection occurs at time $t < s \le t + \Delta$, conditionally to $W^\theta(t)$, is equal to
$$ \sum_{i: X_i(t)=S} \mathbb{P} (Q_i \le P(t+\Delta) \mid Q_i > P(t)) = \sum_{i: X_i(t)=S} \left(1 - \exp\left(-\frac{\beta}{N} W_I^\theta(t)\Delta\right)\right) $$
as the tolerance thresholds $Q_i$ are sampled
from the standard exponential distribution. Finally, this last expression is equiv-
alent to $\frac{\beta}{N} W_S^\theta(t) W_I^\theta(t) \Delta$ as $\Delta \to 0$. Regarding now the computation of the
recovery rate, suppose the individual labeled $i$ is initially infected ($t=0$) or enters
compartment $I$ at time $t > 0$. Then, the probability that a recovery occurs $t <$
$s \le t + \Delta$, conditionally to $W^\theta(t)$, is equal to $\sum_{i: X_i(t)=I} P(L_i \le t + \Delta | L_i > t) =$
$W_I^\theta(t) \times (1 - \exp(-\gamma_I \Delta))$, which is equivalent to $\gamma_I W_I^\theta(t)$ as $\Delta \to 0$.

## 2.2 A deterministic representation of the SIR model

Recall that the population is supposed closed and that individual initial health sta-
tuses, i.e. $X_1(0), \dots, X_N(0)$, are non-random with values in $\{S, I\}$. Then, it
appears that at time $t \ge 0$, the number of susceptible individuals in the population
is equal to the number of initially susceptible individuals (i.e. $X_i(0) = S$) whose
tolerance threshold remains greater than the value of the infectious pressure at time
$t$, i.e. $P(t)$. Thus

$$
W_S^\theta(t) = \sum_{i=1}^{N} 1_{Q_i > P(t), X_i(0)=S} \quad (1)
$$

with the convention $Q_i = 0$ if $X_i(0) = I$. For any $t > 0$, let $\mathcal{F}_{t-}^\theta$ denote the
$\sigma$-algebra generated by $\{W(s), s < t\}$. Note that $P(t)$ is deterministic condi-
tionally on $\mathcal{F}_{t-}^\theta$. Thus, proceeding sequentially along time, it implies that $W_S^\theta(t)$
can be written as a deterministic function of the sequence of tolerance thresholds
$Q_1, \dots, Q_N$.

In order to explicit the number of infected individuals, notice that the com-
partment *I* includes at most two types of individuals: either initially infected in-
dividuals (infected at *t* = 0) who have not recovered yet, or initially susceptible
individuals who are infected but not yet recovered. Thus

$$
W_I^\theta(t) = \sum_{i=1}^{N} 1_{L_i > t, X_i(0)=I} + \sum_{i=1}^{N} 1_{\tau_i^S + L_i > t, Q_i \le P(t), X_i(0)=S} \quad (2)
$$

Recall that $\tau_i^S$ is defined as $\tau_i^S = \inf\{s \ge 0 \mid P(s) \ge Q_i\} = P^{-1}(Q_i)$ with $P^{-1}$ the generalized inverse of $P$. This implies, together with Equation (2), that it is possible to write $W_I^\theta(t)$ as a deterministic function of the sequences of tolerance thresholds and sojourn durations $(Q_1, L_1, \dots, Q_N, L_N)$. Relying on the assumption that the population is closed, it holds that for any $t \ge 0$, $W_S^\theta(t) + W_I^\theta(t) + W_R^\theta(t) = N$. Thus,

$$
W_R^\theta(t) = N - W_S^\theta(t) - W_I^\theta(t) \quad \forall t \ge 0, \quad (3)
$$

so that for all $t > 0$, $W_R^\theta(t)$ can be written as a deterministic function of the se-
quences of tolerance thresholds and sojourn durations $(Q_1, L_1, \dots, Q_N, L_N)$. Let
$\tilde{Z} = (\tilde{Z}_1, \dots, \tilde{Z}_{2N}) := (Q_1, L_1, \dots, Q_N, L_N)$ and $\theta := (\beta, \gamma_I)$. Let $F_k^{-1}$ denote
the quantile function of $\tilde{Z}_k$, $1 \le k \le 2N$. Then for $Z = (Z_1, \dots, Z_{2N})$ uniformly distributed on $[0, 1]^{2N}$, we know that $(F_1^{-1}(Z_1), \dots, F_{2N}^{-1}(Z_{2N}))$ is distributed as $\tilde{Z}$. Note that $Z$ is parameter free. Only the quantile functions $F_k^{-1}$, $k = 1, \dots, 2N$ depend on $\theta$. We thus deduce from Equations (1), (2) and (3) that there exists a deterministic function $f$ such that, for all $t > 0$, $W^\theta(t) := (W_S^\theta(t), W_F^\theta(t), W_R^\theta(t)) \stackrel{D}{=} f(t, \theta, Z)$. If now uncertain parameter $\theta$ is modelled by a random vector $\mathbf{X}$ independent of $Z$, then $(f, Z)$ defines a deterministic representation of $(\mathbf{X}, W^\mathbf{X})$ in the sense that:

$$ (\mathbf{X}, W^{\mathbf{X}}(\cdot)) \stackrel{D}{=} (\mathbf{X}, f(\cdot, \mathbf{X}, Z)). \qquad (4) $$

Such a representation is built sequentially on time, and will be a key stone for performing global sensitivity analysis (see Section 5 for more details).

## 3 Extending Sellke construction

The construction of Sellke [1983] was originally introduced to study the distribution of the final size of epidemics modeled by the SIR model. A natural question is whether such a construction can be extended to more complex compartmental models. The answer is positive. In this section we focus on the extension of Sellke construction to the $SEI_1I_2RS$ model. This model is interesting as it includes loops (individuals may be reinfected) and a branching compartment, i.e. a compartment which leads to more than one compartment, like compartment $E$ in the $SEI_1I_2RS$ model.

### 3.1 Description of the SEI1I2RS model

The $SEI_1I_2RS$ model presented in Figure 3 is a compartmental model with health statuses $\{S, E, I_1, I_2, R\}$ and six different types of transitions $\{(S, E), (E, I_1), (E, I_2), (I_1, R), (I_2, R), (R, S)\}$. The different types of transitions between states are described in Table 1.

<figure><img src="image_4.png" /><figcaption>Figure 3: The $SEI_1I_2RS$ model</figcaption></figure>
This model enables to describe the propagation of an epidemic within a population with the following characteristics: existence of an incubation period for the infected individuals, the presence of two categories of infectious individuals, and the absence of permanent immunity (possible reinfection). For instance, it could be used to model the spread of SARS-CoV-2 by considering individuals in compartment $I_1$ as symptomatic and those in compartment $I_2$ as asymptomatic. To each type of transition we associate a vector $u = (u_1, u_2, u_3, u_4, u_5)$, with $u_i \in \{-1, 0, 1\}$ (see last column of Table 1).

<table><caption>Table 1: Description of the model transitions between states.</caption><thead><tr><th>Type of transition</th><th>Type</th><th>Transition vector</th></tr></thead><tbody><tr><td>(S, E)</td><td>infection</td><td>u<sub>(S,E)</sub> := (-1, 1, 0, 0, 0)</td></tr><tr><td>(E, I<sub>1</sub>)</td><td>first type<br>infection activation</td><td>u<sub>(E,I<sub>1</sub>)</sub> := (0, -1, 1, 0, 0)</td></tr><tr><td>(E, I<sub>2</sub>)</td><td>second type<br>infection activation</td><td>u<sub>(E,I<sub>2</sub>)</sub> := (0, -1, 0, 1, 0)</td></tr><tr><td>(I<sub>1</sub>, R)</td><td>first type<br>recovery</td><td>u<sub>(I<sub>1</sub>,R)</sub> := (0, 0, -1, 0, 1)</td></tr><tr><td>(I<sub>2</sub>, R)</td><td>second type<br>recovery</td><td>u<sub>(I<sub>2</sub>,R)</sub> := (0, 0, 0, -1, 1)</td></tr><tr><td>(R, S)</td><td>reinfection</td><td>u<sub>(R,S)</sub> := (1, 0, 0, 0, -1)</td></tr></tbody></table>

The succession over time of these transitions completely defines the states of the epidemic process. The dynamics depend on a certain number of parameters related to the epidemic characteristics. In the current modeling context, the SEI<sub>1</sub>I<sub>2</sub>RS model depends on six epidemic parameters listed in Table 2.

<table><caption>Table 2: Epidemic parameters of the Markovian SEI<sub>1</sub>I<sub>2</sub>RS model.</caption><thead><tr><th>Parameter Name</th><th>Parameter Role</th></tr></thead><tbody><tr><td>β</td><td>transmission rate</td></tr><tr><td>1/μ<sub>E</sub></td><td>mean sojourn duration in E</td></tr><tr><td>p</td><td>probability for an exposed<br>individual to move to I<sub>1</sub></td></tr><tr><td>1/μ<sub>1</sub></td><td>mean sojourn time in I<sub>1</sub></td></tr><tr><td>1/μ<sub>2</sub></td><td>mean sojourn time in I<sub>2</sub></td></tr><tr><td>1/δ</td><td>mean sojourn time in R</td></tr></tbody></table>

In the following, $\theta := (\beta, \mu_E, p, \mu_1, \mu_2, \delta)$ is defined as the vector of epidemic parameters. The set of values for $\theta$ is denoted by $\mathcal{X} \subseteq \mathbb{R}^6$. We denote by $W(t) := (W_S(t), W_E(t), W_{I_1}(t), W_{I_2}(t), W_R(t))$ so that $W = \{W(t), t \ge 0\}$ represents
the epidemic process corresponding to the SEI₁I₂RS model.

## 3.2 Extension of Sellke construction for the SEI1I2RS model

The construction introduced in Sellke [1983] for the SIR model cannot be directly applied to the SEI₁I₂RS model as for this last model, there exists a branching compartment (compartment E) and the possibility of reinfection (R → S). The existence of the branching compartment E induces a bifurcation thus, preliminary to any transition from E, one has to activate a selection mechanism to choose either E → I₁ or E → I₂. Reinfection introduces loops in the structure of the model. Indeed, it is possible for an individual to exit a compartment and to return back to it after a finite number of transitions. The infection mechanism described in Section 2 has thus to be generalized in a suitable way. This is the purpose of what follows.

As in the original Sellke construction, assume that the N individuals of the closed population are labeled i = 1, ..., N. Recall that the health status of the individual labeled i at time t ≥ 0 is the random variable Xᵢ(t) taking values in V = {S, E, I₁, I₂, R}. As previously, we define for any t > 0 the σ-algebra Fᵢ⁻ generated by {Wθ(s), s < t}. In the following, when there is no ambiguity, we omit superscript θ.

### Infection mechanism: transitions S → E

As in Section 2, the infection mechanism still depends on a pressure function. There are two groups of infectious individuals which exert pressure over susceptible individuals within the population. We define the pressure function as:

$$P: t \mapsto \frac{\beta}{N} \int_{0}^{t} W_{I}(s) ds \qquad (5)$$

with for any s ≥ 0, Wᵢ(s) = Wᵢ₁(s) + Wᵢ₂(s). With this formulation, we assume that individuals in I₁ and I₂ share a common transmission rate β. One may assume instead that each infectious group is characterized by its own transmission rate, say β₁ and β₂. In that case, the pressure function would be given by

$$t \mapsto \int_{0}^{t} \left( \frac{\beta_{1}}{N} W_{I_{1}}(s) + \frac{\beta_{2}}{N} W_{I_{2}}(s) \right) ds.$$

We now describe the infection mechanism. If the individual labeled i is susceptible at initial time t = 0, he is assigned an initial tolerance pressure Qᵢ,₀ > 0 and he gets infected as soon as t ≥ Qᵢ,₀. Otherwise Qᵢ,₀ = 0. Then, each time the individual labeled i becomes susceptible, he is assigned a new tolerance threshold Qᵢ,ⱼ₊₁ > 0 and he remains susceptible as long as the excess pressure remains below that threshold. More precisely, if we denote by ηᵢⱼ,S the entrance times of individual labeled i in compartment S, then for t ≥ ηᵢⱼ,S, the individual labeled
$i$ remains susceptible as long as $P(t) - P(\eta_{i,j}^S) < Q_{i,j+1}$. He gets infected at the first time the pressure function reaches and exceeds $P(\eta_{i,j}^S) + Q_{i,j+1}$, i.e. $\tau_{i,j+1}^S = \inf\{t \ge \eta_{i,j}^S \mid P(t) \ge P(\eta_{i,j}^S) + Q_{i,j+1}\}$. In the following we denote by $Q_i$ the sequence $(Q_{i,0}, Q_{i,1}, \dots)$. We also remark that the pressure function computed at time $t > 0$, $P(t)$, is deterministic conditionally on $\mathcal{F}_{t-}$.

## Handling transitions from compartments *E*, *I*₁, *I*₂ and *R*

Transitions from compartments *E*, *I*₁, *I*₂ and *R* are handled as transition from *I* in the original Sellke construction. We define an intuitive partial order on **V** := {*S*, *E*, *I*₁, *I*₂, *R*}: *S* < *E* < *I*₁ < *R* and *S* < *E* < *I*₂ < *R*. By a slight abuse of notation, we define compartment *I* as the union *I*₁∪*I*₂, leading to *S* < *E* < *I* < *R*. As soon as an individual enters one of the compartments in {*S*, *E*, *I*, *R*}, he is assigned a sojourn duration. At the end of this period, the individual leaves the compartment. In the following, we denote by $L_i^\alpha = (L_{i,0}^\alpha, L_{i,1}^\alpha, \dots)$ the successive sojourn durations of individual labeled *i* in compartment $\alpha \in \{E, I, R\}$, with $L_{i,0}^\alpha = 0 \forall \alpha < X_i(0)$. For an individual in compartment *E*, the transition leads to compartment *I*₁ with a probability *p* in (0, 1) and to compartment *I*₂ with probability $1-p$. The selection mechanism is modeled with a multinomial probability distribution (binomial in the present setting with only two infectious compartments), denoted by $M(\{E \to I_1, E \to I_2\}, \mathbf{p} = (p, 1-p))$. Then, each time an individual leaves compartment *S*, a random variable is sampled from this multinomial distribution to select either $E \to I_1$ or $E \to I_2$. Successive samplings for individual labeled *i* are denoted by $M_i = (M_{i,0}, M_{i,1}, \dots)$, with $M_{i,0} = 0$ if $I_k < X_i(0)$, $k = 1, 2$, $M_{i,0} = E \to I_1$ if $X_i(0) = I_1$ and $M_{i,0} = E \to I_2$ if $X_i(0) = I_2$.

## Retrieving the Markovian case

Similarly to what happens for the construction introduced by Sellke for the SIR model, the transition mechanisms we introduced above lead to a continuous-time Markov chain if one samples tolerance thresholds and sojourn durations adequately. The precise result is stated in Proposition 1 below.

**Proposition 1.** *Assume that*

* the sequences $Q_i$, $L_i^\alpha$, $\alpha \in \{E, I, R\}$ and $M_i$, $1 \le i \le N$ are mutually independent,

* the variables in each of the sequences $Q_i$, $L_i^\alpha$, $\alpha \in \{E, R\}$ and $M_i$, $1 \le i \le N$ are independent and identically distributed,

* for each $1 \le i \le N$, the variables $L_{i,j}^I$ are independent and identically distributed conditionally on $M_{i,j}$,

* the tolerance thresholds are sampled from the standard exponential distribution,
* the sojourn durations in compartment E (resp. R) are sampled from the exponential distribution with mean $1/\mu_E$ (resp. $1/\delta$),

* for each $1 \le i \le N$, the sojourn durations $L_{i,j}^I$ in compartment I are sampled conditionally on $M_{i,j}$ from the exponential distribution with mean $1/\mu_k$ if $M_{i,j} = E \to I_k$, $k = 1, 2$.

Then, the transition rates of the continuous-time homogeneous Markov process $\{W(t), t \ge 0\}$ are given in Table 3 below.

<table><caption>Table 3: Transition rates for the Markovian model.</caption><thead><tr><th>Type of transition</th><th>Rates</th></tr></thead><tbody><tr><td>(S, E)</td><td><sup>β</sup>/<sub>N</sub> &middot; W<sub>S</sub> &middot; (W<sub>I1</sub> + W<sub>I2</sub>)</td></tr><tr><td>(E, I<sub>1</sub>)</td><td>p &middot; μ<sub>E</sub> &middot; W<sub>E</sub></td></tr><tr><td>(E, I<sub>2</sub>)</td><td>(1 &minus; p) &middot; μ<sub>E</sub> &middot; W<sub>E</sub></td></tr><tr><td>(I<sub>1</sub>, R)</td><td>μ<sub>1</sub> &middot; W<sub>I1</sub></td></tr><tr><td>(I<sub>2</sub>, R)</td><td>μ<sub>2</sub> &middot; W<sub>I2</sub></td></tr><tr><td>(R, S)</td><td>δ &middot; W<sub>R</sub></td></tr></tbody></table>

The proof of Proposition 1 is based on classical arguments. It is postponed to Appendix A.

### 3.3 A deterministic representation of the SEI1I2RS model

As in the original Sellke construction, assume that the $N$ individuals of the closed population are labeled $i = 1, \dots, N$. Recall that the health status of the individual labeled $i$ at time $t \ge 0$ is the random variable $X_i(t)$ taking values in $V = \{S, E, I_1, I_2, R\}$. As already mentioned, such an individual is characterized by a vector of sequences $(Q_i, L_i^E, L_i^I, L_i^R, M_i)$. Moreover these sequences are infinite as, due to the presence of loops in the model, an individual may visit a compartment infinitely many times. Let us now define $\tau_{i,j}^S$, $j \ge 1$, as the $j$th exit time of individual labeled $i$ from compartment $S$. We set $\tau_{i,0}^S = 0$ except if individual labeled $i$ is in compartment $S$ at time $t = 0$, in which case we set $\tau_{i,0}^S = \inf\{t > 0 \text{ such that } P(t) \ge Q_{i,0}\}$. For $j \ge 0$, the following recurrence relation holds:

$$ \tau_{i,j+1}^S = \inf\{t > \tau_{i,j}^S + \Delta_{i,j+1}^S \text{ such that } P(t) \ge P(\tau_{i,j}^S + \Delta_{i,j+1}^S) + Q_{i,j+1}\} \quad (6) $$

with $\Delta_{i,j+1}^S$ the cumulated sojourn duration of individual number $i$ in compartments other than $S$ after his $j$th exit from $S$ and before his $(j+1)$th arrival in $S$. It is given by:

$$ \Delta_{i,j+1}^S = L_{i,j}^E + L_{i,j}^I + L_{i,j}^R. \quad (7) $$

Note that $\tau_{i,j}^S + \Delta_{i,j+1}^S$ is the time when the individual labeled $i$ becomes susceptible for the $j+1$-th time. Thus, as long as the infection pressure is below
$P(\tau_{i,j}^S + \Delta_{i,j+1}^S) + Q_{i,j+1}$, this individual remains in $S$. Overall, it turns out that the number of susceptible individuals can be reformulated as:

$$W_S(t) = \sum_{i=1}^{N} \sum_{j \ge 0} \mathbb{1}_{\tau_{i,j}^S + \Delta_{i,j+1}^S \le t < \tau_{i,j+1}^S}. \quad (8)$$

Now, by Equation (6) and System (7), and as $t \mapsto P(t)$ is non decreasing, the sequence $\{\tau_{i,j}^S, j \ge 0\}$ can be written as a deterministic function of $\{Q_i, L_i^E, L_i^I, L_i^R, M_i\}$. We then deduce from Equation (8) that $W_S := \{W_S(t), t \ge 0\}$ can be written as a deterministic function of $\{Q_i, L_i^E, L_i^I, L_i^R\}, i = 1, \dots, N$.

Similarly, exit times of individuals from other compartments obey explicit recurrence formulas. Indeed, for $\alpha \in \{E, I, R\}$, it holds that:

$$\tau_{i,j+1}^{\alpha} = \tau_{i,j}^{\alpha} + \Delta_{i,j+1}^{\alpha} \quad (9)$$

with

$$\Delta_{i,j+1}^{\alpha} = \tau_{i,j+1}^{S} - \tau_{i,j}^{S} + \sum_{\gamma \le \alpha, \gamma \ne S} (L_{i,j+1}^{\gamma} - L_{i,j}^{\gamma}), \quad (10)$$

with $\tau_{i,0}^\alpha = L_{i,0}^\alpha$. Details leading to (10) are postponed to Appendix B. The quantity $\Delta_{i,j+1}^\alpha$ corresponds to the time spent by the individual labeled $i$ outside the $\alpha$ compartment before his $(j+1)$-th visit to $\alpha$, added to the time he spends in $\alpha$ during this new stay. Thus, $\tau_{i,j}^\alpha + \Delta_{i,j+1}^\alpha - L_{i,j+1}^\alpha$ corresponds to the time at which this individual enters the compartment $\alpha$ so that for $t \in [\tau_{i,j}^\alpha + \Delta_{i,j+1}^\alpha - L_{i,j+1}^\alpha, \tau_{i,j+1}^\alpha)$, $X_i(t) = \alpha$. Thus:

$$W_{\alpha}(t) = \sum_{i=1}^{N} \sum_{j \ge 0} \mathbb{1}_{\tau_{i,j}^{\alpha} + \Delta_{i,j+1}^{\alpha} - L_{i,j+1}^{\alpha} \le t < \tau_{i,j+1}^{\alpha}}. \quad (11)$$

Now, by Equation (9) and System (10), and as $t \mapsto P(t)$ is non decreasing, each of both sequences $\{\tau_{i,j}^\alpha, j \ge 0\}$ and $\{\Delta_{i,j}^\alpha, j \ge 1\}$ can be written as a deterministic function of $\{Q_i, L_i^E, L_i^I, L_i^R, M_i\}$. We then deduce from Equation (11) that $W_\alpha := \{W_\alpha(t), t \ge 0\}$ can be written as a deterministic function of $\{Q_i, L_i^E, L_i^I, L_i^R, i = 1, \dots, N\}$.

The description of the different exit times through recurrence formulas enables to build an algorithm that describes any trajectory of the epidemic process. Such an algorithm is presented in Algorithm 1 below. As in Table 1, for any $(\alpha, \alpha') \in V \times V$, the transition vector $u_{\alpha,\alpha'}$ is a 5-dimensional vector.

In the following, for each $1 \le i \le N$, we sample the sequence $L_i^I$ in the following manner: we first sample $M_{i,j}$, then we sample $L_{i,j}^{I_1}$ and $L_{i,j}^{I_2}$ independently from each other and independently from $M_{i,j}$, finally we compute $L_{i,j}^I = L_{i,j}^{I_1} \mathbb{1}_{M_{i,j}=E \to I_1} + L_{i,j}^{I_2} \mathbb{1}_{M_{i,j}=E \to I_2}$. Let $Z_Q = (Z_{Q,k})_{k \ge 1}$, $Z_\alpha = (Z_{\alpha,k})_{k \ge 1}$, $\alpha \in \{E, I_1, I_2, R\}$ and $Z_M = (Z_{M,k})_{k \ge 1}$ be independent sequences of i.i.d. random variables uniformly distributed on $[0, 1]$. Let $\tilde{Z} := (Q, L^E, L^{I_1}, L^{I_2}, L^R, M)$.
Arguing as in Section 2.2, it is possible to exhibit a deterministic map $G_{\theta}$ (defined from quantile functions of independent marginals of $\tilde{Z}$), depending on the vector of parameters $\theta$, such that: $\tilde{Z} \stackrel{D}{=} G_{\theta}(Z_{Q}, Z_{E}, Z_{I_{1}}, Z_{I_{2}}, Z_{R}, Z_{M})$. Then in Algorithm 1, $Z$ is identified with the sequences of seeds $RG_{Q}, RG_{\alpha}, \alpha \in \{E, I_{1}, I_{2}, R\}, RG_{M}$ used to initialize the pseudorandom number generators of $Z_{Q}, Z_{\alpha}, \alpha \in \{E, I_{1}, I_{2}, R\}$ and $M$. Algorithm 1 defines a deterministic function $f$ such that, for all $t \ge 0$, $W^{\theta}(t) \stackrel{D}{=} f(t, \theta, Z)$ with $\theta = (\beta, \mu_{E}, p, \mu_{1}, \mu_{2}, \delta)$. Also, if $\theta$ is sampled from a random vector $\mathbf{X}$ with $\mathbf{X}$ and $Z$ mutually independent, then $(f, Z)$ defines a deterministic representation of $(\mathbf{X}, W^{\mathbf{X}})$ in the sense that:

$$ (\mathbf{X}, W^{\mathbf{X}}(\cdot)) \stackrel{D}{=} (\mathbf{X}, f(\cdot, \mathbf{X}, Z)). $$

As already mentioned, such a representation is a key stone for global sensitivity analysis (see Section 5 for more details).

**Algorithm 1:**

**Data:** $\xi_{0}(=w(0)), T$

**inputs:** $\theta = (\beta, \mu_{E}, p, \mu_{1}, \mu_{2}, \delta)$,

$$ Z = \{RG_{Q}, RG_{E}, RG_{I_{1}}, RG_{I_{2}}, RG_{R}, RG_{M}\} $$

**output:** $w = \{w(t); t \in [0, T]\}$

/* Initialization */

1

$t \leftarrow 0, w(t) \leftarrow \xi_{0}, P_{S}(t) \leftarrow 0, j \leftarrow 0$

2 $N \leftarrow w_{S}(t) + w_{E}(t) + w_{I_{1}}(t) + w_{I_{2}}(t) + w_{R}(t)$

/* Draw thresholds for initially susceptible individuals */

3

4 **for** $i$ **such that** $X_{i}(0) = S$ **do**

5 **Draw** $Q_{i,j}$ **using the seed** $RG_{Q}$

6 **end**

7 **for** $\alpha \in \{E, I_{1}, I_{2}, R\}$ **do**

/* Draw sojourn variables for individuals initially in */

$E, I_{1}, I_{2}, R$

8

9 **for** $i$ **such that** $X_{i}(0) = \alpha$ **do**

10 **Draw** $L_{i,j}^{\alpha}$ **using the seed** $RG_{\alpha}$

11 **end**

12 **end**

12 **end**

12 **end**
while $t < T$ do

/* Pick the smallest threshold */

13 $\tilde{Q}_{min} \leftarrow \min\{Q_{i,j} : i \text{ such that } X_i(t) = S\}$

14 $P_{deriv} \leftarrow (\beta/N) \times (w_{I_1}(t) + w_{I_2}(t))$

/* Compute putative next infection time */

15 $\tau^S \leftarrow t + (\tilde{Q}_{min} - P_S(t)) / P_{deriv}$

/* Compute putative next exit time from other compartments */

16 for $\alpha \in \{E, I_1, I_2, R\}$ do

17 $\lfloor \tau^\alpha \leftarrow \min\{\tau^\alpha_{i,j} : i \text{ such that } X_i(t) = \alpha\} \rfloor$

/* Compute next event time */

18 $\bar{\tau} \leftarrow \min\{\tau^\alpha : \alpha \in \{S, E, I_1, I_2, R\}\}$

/* Find next exit compartment */

19 $\alpha^* \leftarrow \operatorname{argmin}\{\tau^\alpha \in \{S, E, I_1, I_2, R\}\}$

/* Compute pressure value at next exit time */

20 $P_S(\bar{\tau}) \leftarrow P_S(t) + P_{deriv} \times (\bar{\tau} - t)$

/* Identify the individual that makes the transition */

21 if $\alpha^* = S$ then

22 $i^* \leftarrow \operatorname{argmin}\{Q_{i,j} : X_i(t) = S\}$

23 Delete $Q_{i^*,j}$

24 else if $\alpha^* \in \{E, I_1, I_2, R\}$ then

25 $i^* \leftarrow \operatorname{argmin}\{\tau^\alpha_{i,j} : X_i(t) = \alpha^*\}$

26 Delete $\tau^\alpha_{i^*,j}$

/* Pick randomly destination compartment for a transition from E */

27 if $\alpha^* = E$ then

28 $\lfloor \text{Using the seed } RG_M, \text{ draw } (E, I_1) \text{ or } (E, I_2) \text{ with probability } p \\ \text{or } 1-p \text{ respectively.} \rfloor$

29 else

/* Destination compartment for transitions from other compartments */

30 $\gamma^* \leftarrow E$ if $\alpha^* = S$

31 $\gamma^* \leftarrow R$ if $\alpha^* = I_1$ or $\alpha^* = I_2$

32 $\gamma^* \leftarrow S$ if $\alpha^* = R$

/* Updates */

33 $w(\bar{\tau}) \leftarrow w(t) + \mathbf{u}_{\alpha^*, \gamma^*}$

34 $t \leftarrow \bar{\tau}, j \leftarrow j+1$

35 if $\gamma^* = S$ then

36 $\lfloor \text{Draw } Q_{i^*,j} \text{ using the seed } RG_Q \rfloor$

37 $\tilde{Q}_{i^*,j} \leftarrow P(t) + Q_{i^*,j}$

38 if $\gamma^* \in \{E, I_1, I_2, R\}$ then

39 $\lfloor \text{Draw } L^\gamma_{i^*,j} \text{ using the seed } RG_\gamma \rfloor$

40 $\tau^\gamma_{i^*,j} \leftarrow t + L^\gamma_{i^*,j}$

14
In the framework of Proposition 1, the sequences of independent and identi-
cally distributed random variables $Q_i, L_i^E, L_i^{I_1}, L_i^{I_2}$ and $L_i^R$ are mutually indepen-
dent and all sampled from an exponential distribution, and the resulting process
$\{W(t), t \ge 0\}$ is a homogeneous continuous-time Markov chain. Now, if the se-
quences of independent and identically distributed random variables $Q_i, L_i^E, L_i^{I_1}$,
$L_i^{I_2}$ and $L_i^R$ are mutually independent and if at least one of them is sampled from
any non negative continuous probability distribution, different from the exponen-
tial distribution, and the other ones are sampled from an exponential distribution,
then the corresponding process is non-Markovian.

4 Simulations

In this section, we present simulations of the SEI1I2RS model, obtained from Algorithm 1. We highlight differences between outputs depending on whether the model is Markovian or not. For this purpose, we consider the following epidemic scenario inspired from COVID-19 pandemic. The population is assumed to be composed of 2505 individuals with 5 initially exposed individuals (i.e. $W_E(0) = 5$). Population size is chosen arbitrary but not too large in order to avoid unnecessary computation burdens for our study which is methodological and does not claim to be a real case study. Regarding the epidemic characteristics, we assume that the incubation period $1/\mu_E$ fluctuates uniformly over [4.5, 5.8] with mean 5.1 [Lauer et al., 2020]. Infectious periods $1/\mu_1$ and $1/\mu_2$ are respectively fixed to $1/\mu_1 = 5$ and $1/\mu_2 = 5$ according to Davies et al. [2020] (see Supplementary Table 1 in their paper). We then set the transmission rate to $\beta = 0.442$, so that for $1/\mu_1 = 5$ and $1/\mu_2 = 5$, the basic reproduction number $R_0 := \beta / (p \times \mu_1 + (1-p) \times \mu_2) = 2.21$, which seems reasonable [Locatelli et al., 2021]. In addition, we set the proportion $p$ to 18, 1% (see Prague et al. [2020]). In order to set the average period of loss of immunity $1/\delta$, we rely on Schuler et al. [2021]. In this paper, authors conclude that reinfection occurs from 90 to 180 days after recovery. Therefore, we set $1/\delta = 135$ corresponding to the mean of a uniform distribution over [90, 180].

For these parameter values, two scenarios are considered for simulations, which correspond to the Markovian framework and a non-Markovian alternative. In both cases, sojourn times in a given compartment are independently and identically distributed (i.i.d.) from Gamma distributions. For a Gamma distribution, length of memory can be read in the shape parameter. Indeed, a Gamma distribution with shape parameter equal to one coincides with an exponential distribution; thus it is used to model the Markovian setting characterized by memoryless. To model scenarios with memory, we use in the following Gamma distributions with shape parameter greater than one. More precisely, we compare two scenarios. The Markovian scenario corresponds to exponential sampling of sojourn times in E, $I_1$, $I_2$ and R, with parameter chosen to fit respectively mean duration in E, $I_1$, $I_2$ and R. We also consider a non-Markovian scenario, with sojourn times in E, $I_1$ and $I_2$ sampled from a Gamma distribution with shape parameter equal to four (see [Davies
et al., 2020, Supplementary Table 1)) and scale parameter chosen to fit respectively mean duration in E, $I_1$ and $I_2$. The shape parameter for sojourn times in R is chosen equal to 1. We assume that the immune system of individuals completely loses immunity from one infection to the other. Tolerance thresholds are drawn in both Markovian and non-Markovian experiments from an exponential distribution with parameter 1.

<figure><img src="image_5.png" /><figcaption>Figure 4: Evolution over time of quantiles of order 0.1 to 0.9 of the number of individuals in $I_2$ (top) and the one of recovered individuals (bottom). Quantiles are computed from 1000 independent trajectories in both Markovian (left) and non-Markovian (right) frameworks.</figcaption></figure>

Using Algorithm 1, we simulate independently 50 sample paths or trajectories of the epidemic process over time period [0, 500], for both Markovian and non-Markovian scenarios described above. Internal noise in Algorithm 1 is represented by $Z = (Q, L^E, L^{I_1}, L^{I_2}, L^R, M)$. Each realization $q^{(i)}$ (respectively $\ell^{E(i)}, \ell^{I_1(i)}, \ell^{I_2(i)}, \ell^{R(i)}$ and $m^{(i)}$) of Q (respectively $L^E, L^{I_1}, L^{I_2}, L^R$ and $M$) is obtained from a pseudorandom number generator initialized with a random seed $z_1^{(i)}$ (respectively $z_2^{(i)}, z_3^{(i)}, z_4^{(i)}, z_5^{(i)}$ and $z_6^{(i)}$). Within each scenario, two trajectories differ only through the realization of Z. These simulations are provided in Figure
C.1 in Appendix C. Differences can be observed from one scenario to the other. However, to limit the sampling impact of any conclusion, we simulate 1000 independent trajectories in both Markovian and non-Markovian frameworks, and plot in Figure 4 the evolution over time of quantiles of different order (ranging from 0.1 to 0.9) for the number of individuals in $I_2$ (top) and for the number of recovered individuals (bottom). On Figure 4, we read that out of our thousand simulations, at least 50% show a trend towards extinction for the non-Markovian scenario (top right), whereas this proportion falls to 20% for the Markovian scenario (top left). We also remark a significant difference in the peak for the dynamics of individuals in $I_2$ between the two scenarios (top line). The peak is significantly higher in the non-Markovian setting (top right) due to memory effect induced by the Gamma sampling of sojourn times.

With these simulations, we highlight the importance of modeling memory effects as they truly influence model outputs. To better account for the impact of memory effects, we perform in the next section a global sensitivity analysis.

# 5 Global sensitivity analysis of the SEI₁I₂RS model

This section is devoted to global sensitivity analysis of the SEI₁I₂RS model. We first review variance based sensitivity analysis in Section 5.1. Then in Section 5.2 we explain how to perform a global sensitivity analysis for a stochastic model based on a deterministic representation of the model, i.e. a deterministic function of the epidemic parameters and the variables modelling the internal noise. Finally in Section 5.3, we present numerical experiments for global sensitivity analysis of the SEI₁I₂RS model.

## 5.1 Sobol' indices

Let $\mathbf{X} = (X_1, \dots, X_p)$ with $p \in \mathbb{N}^*$ be a random vector with known probability distribution. Let $f$ be a real function taking $p$ real arguments. Assume $\mathrm{Var}(f(\mathbf{X})) < +\infty$. Moreover, suppose that $X_1, \dots, X_p$ are mutually independent so that Sobol'-Hoeffding decomposition [Sobol', 1993, Hoeffding, 1948] yields:

$$ \mathrm{Var}(f(\mathbf{X})) = \sum_{u \subseteq \{1, \dots, p\}, u \neq \emptyset} \mathrm{Var}(f_u(\mathbf{X}_u)), \quad (12) $$

where, for each $u \subseteq \{1, \dots, p\}$, $f_u$ is a function of $X_u =: \{X_j, j \in u\}$ such that $\mathbb{E}[f_u(\mathbf{X}_u) | \mathbf{X}_v] = 0$ for every $v \not\subseteq u$. Based on (12), Sobol' [1993] introduced variance-based sensitivity indices, nowadays known as Sobol' indices. In this paper we focus on first-order and total Sobol' indices.

**First-order Sobol' indices** Let $u$ be a nonempty subset of $\{1, \dots, p\}$. The first-order Sobol' index of $\mathbf{X}_u$ is defined as:

$$ S_{\mathbf{X}_u} = \frac{\mathrm{Var}(\mathbb{E}[f(\mathbf{X}) | \mathbf{X}_u])}{\mathrm{Var}(f(\mathbf{X}))}. \quad (13) $$
Such an index is a ratio of the variance of the part of the output due only to the variation of $X_u$ and the global variance. It is the so-called main effect of $X_u$.

**Total Sobol' indices** The total Sobol index of $X_u$ is defined as

$$ST_{X_u} = 1 - \frac{\mathrm{Var}(\mathbb{E}[f(\mathbf{X}) | \mathbf{X}_{\sim u}])}{\mathrm{Var}(f(\mathbf{X}))}, \quad (14)$$

where $X_{\sim u} = \{X_j, j \notin u\}$. This index summarizes the main effect of $X_u$ as well as its interactions with all other inputs of the model. Indeed, for any $j \in \{1, \dots, p\}$, we have:

$$ST_j = \sum_{u \subseteq \{1, \dots, p\} \text{ such that } j \in u} S_{X_u}.$$

**Remark 2.** The definition of first-order and total Sobol' indices can be extended to models with vectorial or functional output. For a vectorial output $(Y_1, \dots, Y_p)$, and for any input parameter $X_i$, it is possible to compute the so-called aggregated first-order (or total) Sobol' index by computing the weighted sum of first-order (or total) Sobol' indices for each scalar output $Y_j$ weighted by the variance of $Y_j$. We refer to Lamboni et al. [2011], Gamboa et al. [2014] for the explicit definition. 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. As an alternative, first-order and total indices defined in Equations (13) and (14) can be computed at each discretized time.

## 5.2 Deterministic representation of stochastic models

Let $\mathcal{X}$, $\mathcal{Y}$ and $\mathcal{Z}$ be three nonempty measurable spaces. A stochastic model $g$ is defined as follows. For each set of input parameters $x \in \mathcal{X}$, the model output $g(x, \cdot)$ is a random variable with values on $\mathcal{Y}$, so that a realization of such an output is under the form $g(x, \omega)$, where $\omega$ belongs to some $\sigma$-algebra. In the paradigm of global sensitivity analysis, uncertain parameters are modeled by a random vector $\mathbf{X}$. In the following, we make the usual assumption that for every $x \in \mathcal{X}$, $\mathbf{X}$ and $g(x, \cdot)$ are mutually independent. In other words, uncertain parameters are independent from internal noise.

**Definition 1.** [Kouye et al., 2024, Section 2.2] Let $f: \mathcal{X} \times \mathcal{Z} \to \mathcal{Y}$ be a deterministic function. Let $Z$ be a random variable with values in $\mathcal{Z}$, whose probability distribution is known, with $Z$ and $\mathbf{X}$ mutually independent. Then $(f, Z)$ is said to be a deterministic representation of the stochastic model $g$ if $(\mathbf{X}, g(\mathbf{X}, \cdot)) \stackrel{\mathcal{D}}{=} (\mathbf{X}, f(\mathbf{X}, Z))$, where $\stackrel{\mathcal{D}}{=} \text{means equality in distribution. The random variable } Z$ controls the internal noise of the stochastic model.

**Remark 3.** In particular, if for every set of input parameters $x \in \mathcal{X}$, $g(x, \cdot) \stackrel{\mathcal{D}}{=} f(x, Z)$, then $(\mathbf{X}, g(\mathbf{X}, \cdot)) \stackrel{\mathcal{D}}{=} (\mathbf{X}, f(\mathbf{X}, Z))$.
A strategy to perform global sensitivity analysis of a stochastic model *g* is to first identify a deterministic representation (*f*, *Z*) of *g* and then to compute Sobol' indices as in Section 5.1 by considering the deterministic model *f* with augmented input vector (X, Z). Let *j* ∈ {1, ..., *p*}. Then the impact of the different sources of uncertainty on model output are measured by the computation of Sobol' indices: *S*<sub>X<sub>j</sub></sub> quantifies the impact of *X*<sub>j</sub>, *S*<sub>Z</sub> measures the impact of *Z*, the variable that controls the stochasticity of the model. We also calculate the total Sobol' index *ST*<sub>X<sub>j</sub></sub>, in order to quantify the impact of *X*<sub>j</sub> alone or in interaction with other epidemic parameters *X*<sub>j'</sub>, *j' ≠ j* or with *Z*, and the total Sobol' index *ST*<sub>Z</sub>.

In the next section, we apply this methodology to perform a global sensitiv-
ity analysis of simulations of the SEI<sub>1</sub>I<sub>2</sub>RS model obtained from Algorithm 1
presented in Section 3.3.

**5.3 Sensitivity analysis**

In this section, we perform a global sensitivity analysis for two different quantities of interest in the SEI₁I₂RS model. We compare the results obtained in a Markovian framework and in a non-Markovian framework. As the results are different, this highlights the importance of taking into account the memory (if any) effect when modeling an epidemic. Throughout this section, we can think about I₁ and I₂ as the compartments for asymptomatic and symptomatic infectious individuals. Each epidemic parameter is sampled from a uniform distribution whose support is set according to the literature on COVID-19. Indeed, even if we do not pretend to consider a real case study, we aim to consider a realistic setting.

Two quantities of interest are under study. The first quantity is the dynamics
over time of infected individuals, namely W<sub>E</sub>(t), t ∈ [0, 500]. The second quan-
tity is the moment the peak of the number individuals in I<sub>2</sub> is reached, namely
argmax<sub>t∈[0,500]</sub> W<sub>I<sub>2</sub></sub>(t). In the context of COVID-19, this corresponds to the mo-
ment the peak of the number of symptomatic individuals is reached. In practice,
this quantity is of particular interest. It is linked to the overload of the medical
system, which proved to be a major problem during the COVID-19 epidemic.

**Numerical setting**

In our numerical experiments, we consider the Markovian and the non-Markovian
frameworks described in Section 4. The size of the population, as well as initial
conditions, are also chosen as in Section 4. Epidemic parameters β, p, μ<sub>E</sub>, μ<sub>1</sub>, μ<sub>2</sub>
and δ are sampled uniformly (the range of variation for each parameter is described
in Table 4).
<table><caption>Table 4: Nominal values and range of variation of uncertain epidemic parameters</caption><thead><tr><td>Parameter Name</td><td>Nominal value</td><td>Range of variation</td></tr></thead><tbody><tr><td>β</td><td>0.442</td><td>[0.23, 0.53]</td></tr><tr><td>1/μ<sub>E</sub></td><td>5.1 days</td><td>[4.5, 5.8]</td></tr><tr><td>p</td><td>0.181</td><td>[0.167, 0.192]</td></tr><tr><td>1/μ<sub>1</sub></td><td>5 days</td><td>[5, 8]</td></tr><tr><td>1/μ<sub>2</sub></td><td>5 days</td><td>[5, 8]</td></tr><tr><td>1/δ</td><td>135 days</td><td>[90, 180]</td></tr></tbody></table>

To sample the variable Z, that controls internal noise, we sample the seed of pseudo-random number generator uniformly in {1, ..., 10⁹}. We compute Sobol' indices by using the function soboljansen of the R-package sensitivity [Iooss et al., 2020]. Each evaluation of Sobol' indices uses a n-sample of (p, μ<sub>E</sub>, μ<sub>1</sub>, μ<sub>2</sub>, δ, Z) with n = 2500. The time-dependent output W<sub>E</sub>(t), t ∈ [0; 500] is evaluated on a regular grid of size 1000 denoted by t<sub>0</sub> = 0 < t<sub>1</sub> < ... < t<sub>999</sub> < t<sub>1000</sub> = 500. Then we can compute Sobol' indices at each time t<sub>i</sub> of the grid, or alternatively compute aggregated Sobol' indices (see Remark 2 for a definition). To account for variability due to input sampling, each Sobol' index is evaluated N<sub>rep</sub> = 30 times.

**Results for W<sub>E</sub>(t), t ∈ [0, 500]** Figure 5 shows the evolution over time of first-order (top) and total (bottom) Sobol' indices. For both Markovian and non-Markovian scenarios we consider, we note that, except during the first phase of the epidemic, when the impact of internal noise Z is the strongest, β is the most important input. Only slight differences can be observed on the evolution of the influence of β from one scenario to the other. Though, a significant difference can be observed for the impact of internal noise Z between both scenarios. We observe that Z has much more impact in the Markovian scenario. This is even clearer on Figure 6, where we have plotted first-order and total aggregated Sobol' indices to summarize the influence of inputs on the dynamical output W<sub>E</sub>(t), t ∈ [0, 500].
<figure><img src="image_6.png" /><figcaption>Figure 5: Time evolution of the median of a sample of $N_{rep} = 30$ realizations of first-order (top) and total (bottom) Sobol' index estimators in Markovian (left) and non-Markovian (right) frameworks. The estimators are built from an input sample of size $n = 2500$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Boxplots obtained from $N_{rep} = 30$ realizations of first-order (left) and total (right) aggregated Sobol' index estimators in our Markovian (red) and non-Markovian (blue) scenarios, for the output defined as $\{W_E(t), t \in [0, 500]\}$. The estimators are built from an input sample of size $n = 2500$.</figcaption></figure>
**Results for argmax<sub>t∈[0,500]</sub> W<sub>12</sub>(t)** On Figure 7 are plotted the results for first-order and total Sobol' index estimation.

<figure><img src="image_8.png" /><figcaption>Figure 7: Boxplots obtained from $N_{rep} = 30$ realizations of first-order (top) and total (bottom) Sobol' index estimators in our Markovian (red) and non-Markovian (blue) scenarios, for the output defined as $\operatorname{argmax}_{t \in [0,500]} W_{I_2}(t)$. The estimators are built from an input sample of size $n = 2500$.</figcaption></figure>

We note that $\beta$ is always the most important input for this quantity of interest. From the plots of first-order indices, $\beta$ explains a little bit more than 37.5% of the variance of the output in the Markovian scenario, while in the non-Markovian one, this proportion increases up to 55%. Considering the total Sobol' indices, the dominance of the influence of $\beta$ is confirmed. In addition, there are significant differences in the total effects of all other input parameters between the Markovian scenario and the non-Markovian one. The larger difference is observed in the total effect of $Z$. The impact of internal noise is greater in the Markovian scenario than in the non-Markovian one. To highlight the pure interaction effects, we compute for each input parameter the value of the total Sobol' index minus the one of the first-order Sobol' index. Corresponding boxplots are plotted in Figure 8.
<figure><img src="image_9.png" /><figcaption>Figure 8: Boxplots obtained from $N_{rep} = 30$ realizations of pure interaction effects in our Markovian (red) and non-Markovian (blue) scenarios. The estimators are built from an input sample of size $n = 2500$.</figcaption></figure>

## Impact of shape parameter uncertainty on model variability

So far, we have assumed that the shape parameters for the Gamma distributions of sojourn times in compartments *E*, *I*₁, *I*₂, and *R* are fixed. However, in practice, determining these parameters is not straightforward as they are also subject to uncertainty. In this section, we aim to quantify the impact of these uncertainties on the global variability of outputs from the *SEI*₁*I*₂*RS* model. To do so, we follow [Lauer et al., 2020, Table 2 in Appendix] and assume that the shape parameters of the different Gamma distributions vary uniformly over [3.585, 13.865]. The other uncertain parameters are still sampled uniformly and independently, with respective variation ranges given in Table 4. Sobol' indices are estimated from a *n*-sample of input parameters of size *n* = 2500 using the function *soboljansen* from the R-package *sensitivity*. As outputs, we still consider *W*_E(*t*), *t* ∈ [0, 500] and argmax*_{t∈[0,500]}* *W*_{I₂}(*t*).

**Results for W_E(t), t ∈ [0, 500]** Figure 9 displays the results of the dynamic Sobol' indices for both first-order (main) and total effects. A first remark is that the estimated dynamics in Figure 9 are less smooth than the ones in Figure 5. The main reason for that is probably that we kept the input sample size unchanged while the input space dimension increased from 7 to 11 (4 shape parameters). The analysis of Figure 9 reveals that the epidemic parameter $\beta$ has the greatest overall influence on output variability, closely followed by $\delta$ on the period $t \in [200, 500]$. Concerning the intrinsic randomness $Z$, its main effect decreases from initial date to $t = 200$, while its total effect strongly decreases until time 50 then increases untill time 200 and finally seems to stabilize around 0.70. Comparing left and right subplots in Figure 9, we remark that shape parameters are not that influent at the start of the epidemics. Their influence strongly increases around time $t = 200$ and
this influence is mainly due to interaction effects. The shape parameter associated
with the sojourn time distribution in *R* turns out to be one of the most important
parameter for the output $W_E(t)$ on the period $t \in [200, 500]$.

<figure><img src="image_10.png" /><figcaption>Figure 9: Time evolution of the median of a sample of $N_{rep} = 50$ realizations of first-order (left) and total (right) dynamical Sobol' index estimators. The estimators are built from an input sample of size $n = 2500$.</figcaption></figure>

**Results for argmax<sub>t∈[0,500]</sub> W<sub>I<sub>2</sub></sub>(t)** For the peak in *I*<sub>2</sub>, Figure 10 shows that, apart from *β*, uncertain parameters (including *Z*) have no significant main effect. However, for all uncertain parameters, the total effect is much larger than the main effect, indicating strong interaction effects. This suggests that the shape parameters are just as important as the epidemic parameters, with the exception of *β*. Overall, interactions play a predominant role in the variability of this output.
<figure><img src="image_11.png" /><figcaption>Figure 10: Boxplots obtained from $N_{rep} = 30$ realizations of first-order (green) and total (red) Sobol' index estimators for the output defined as $\operatorname{argmax}_{t \in [0,500]} W_{I_2}(t)$. The estimators are built from an input sample of size $n = 2500$.</figcaption></figure>

In summary, the full sensitivity analysis highlights the significant influence of shape parameters, mainly due to interaction effects. As a future work, it would be interesting to investigate the strength of second-order interactions.

# 6 Conclusion

In epidemiology, compartmental modeling plays a key role in the quantitative study and understanding of the spread of epidemics. Two frameworks are available to modelers in stochastic modeling: the Markovian framework and the non-Markovian framework. The Markovian framework consists in using Markovian processes such as continuous-time Markov chains, i.e. memoryless processes, to model epidemics. The advantage of this framework is that an arsenal of mathematical tools is available to facilitate theoretical analysis. On the other hand, the non-Markovian framework is more difficult to deal with theoretically, but it is more general and more realistic, as it does not rely on the absence of memory. The choice of one of these frameworks has consequences for the lessons or conclusions that can be drawn from the study of the resulting epidemic process.

In this paper, we use sensitivity analysis to highlight the differences between different scenarios, with or without memory effects. These differences appear in the influence of epidemic parameters but also in the influence of intrinsic noise on epidemic dynamics. The key stone for our sensitivity analysis is the extension of Sellke construction to more complex models. Note that, although we focused the presentation on the $SEI_1I_2RS$ model, our construction easily extends to any closed-population compartmental model with loops and parallel compart-
ments. Our extension of Sellke construction inherits its ability to simulate exact trajectories of a same model in a Markovian context (by sampling exponential distributions for sojourner times) but also in a non-Markovian context (by replacing exponential sampling by, e.g., Gamma sampling).

In future work, it would be interesting to compare simulation results obtained from Algorithm 1 based on our extension of Sellke construction to other algorithms dedicated to non-Markovian dynamics (e.g. Vestergaard and Génois [2015], Boguñá et al. [2014], Masuda and Rocha [2018]). In particular, robustness to internal noise may vary from one algorithm to the other, with an impact on global sensitivity analysis results, as studied in Kouye et al. [2024] in the Markovian context.

## References

Håkan Andersson and Tom Britton. *Stochastic Epidemic Models and Their Statistical Analysis*, volume 151 of Lecture Notes in Statistics. Springer New York, New York, NY, 2000. ISBN 978-0-387-95050-1 978-1-4612-1158-7. doi: 10.1007/978-1-4612-1158-7.

Marian Boguñá, Luis F. Lafuerza, Raúl Toral, and M. Ángeles Serrano. Simulating non-Markovian stochastic processes. *Phys. Rev. E*, 90:042108, Oct 2014. doi: 10.1103/PhysRevE.90.042108.

Tobias Brett and Tobias Galla. Stochastic Processes with Distributed Delays: Chemical Langevin Equation and Linear-Noise Approximation. *Phys. Rev. Lett.*, 110:250601, Jun 2013. doi: 10.1103/PhysRevLett.110.250601.

Tom Britton and Etienne Pardoux. *Chapter 3 General Closed Models*, pages 43–57. Springer International Publishing, Cham, 2019a. ISBN 978-3-030-30900-8. doi: 10.1007/978-3-030-30900-8_3.

Tom Britton and Etienne Pardoux, editors. *Stochastic Epidemic Models with Inference*, volume 2255 of Lecture Notes in Mathematics. Springer International Publishing, Cham, 2019b. ISBN 978-3-030-30899-5 978-3-030-30900-8. doi: 10.1007/978-3-030-30900-8.

Nicholas G. Davies, Petra Klepac, Yang Liu, Kiesha Prem, Mark Jit, CM-MID COVID-19 working group, Carl A. B. Pearson, Billy J. Quilty, Adam J. Kucharski, Hamish Gibbs, Samuel Clifford, Amy Gimma, Kevin Van Zandvoort, James D. Munday, Charlie Diamond, W. John Edmunds, Rein M. G. J. Houben, Joel Hellewell, Timothy W. Russell, Sam Abbott, Sebastian Funk, Nikos I. Bosse, Yueqian Fiona Sun, Stefan Flasche, Alicia Rosello, Christopher I. Jarvis, and Rosalind M. Eggo. Age-dependent effects in the transmission and control of COVID-19 epidemics. *Nature Medicine*, 26(8):1205–1211, August 2020. ISSN 1078-8956, 1546-170X. doi: 10.1038/s41591-020-0962-9.
Francesco Di Lauro, Wasiur R. KhudaBukhsh, István Z. Kiss, Eben Kenah, Max Jensen, and Grzegorz A. Rempła. Dynamic survival analysis for non-Markovian epidemic models. *Journal of The Royal Society Interface*, 19(191): 20220124, 2022. doi: 10.1098/rsif.2022.0124.

Stewart N. Ethier and Thomas 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.

Fabrice Gamboa, Alexandre Janon, Thierry Klein, and Agnès Lagnoux. Sensitivity analysis for multidimensional and functional outputs. *Electronic Journal of Statistics*, 8(1):575–603, 2014. doi: 10.1214/14-EJS895.

Daniel 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.

Gerrit Großmann, Michael Backenköhler, and Verena Wolf. Why ODE models for COVID-19 fail: Heterogeneity shapes epidemic dynamics. *medRxiv*, 2021. doi: 10.1101/2021.03.25.21254292.

Wassily Hoeffding. A Class of Statistics with Asymptotically Normal Distribution. *The Annals of Mathematical Statistics*, 19(3):293 – 325, 1948. doi: 10.1214/aoms/1177730196.

Thomas House. Non-markovian stochastic epidemics in extremely heterogeneous populations. *Math. Model. Nat. Phenom.*, 9(2):153–160, 2014. doi: 10.1051/mmnp/20149210.

Bertrand Iooss, Sebastien Da Veiga, Alexandre Janon, Gilles Pujol, with contributions from Baptiste Broto, Khalid Boumhaout, Thibault Delage, Reda El Amri, Jana Fruth, Laurent Gilquin, Joseph Guillaume, Loic Le Gratiet, Paul Lemaitre, Amandine Marrel, Anouar Meynaoui, Barry L. Nelson, Filippo Monari, Roelof Oomen, Oldrich Rakovec, Bernardo Ramos, Olivier Roustant, Eunhye Song, Jeremy Staum, Roman Sueur, Taieb Touati, and Frank Weber. sensitivity: Global Sensitivity Analysis of Model Outputs, 2020. R package version 1.19.0.

William Ogilvy Kermack and Anderson G McKendrick. A contribution to the mathematical theory of epidemics. *Proceedings of the royal society of london. Series A, Containing papers of a mathematical and physical character*, 115(772): 700–721, 1927.

Henri Mermoz Kouye, Gildas Mazo, Clémentine Prieur, and Elisabeta Vergu. Performing global sensitivity analysis on simulations of a continuous-time
Markov chain model motivated by epidemiology. *Computational and Applied Mathematics*, 43(7):409, 2024. ISSN 1807-0302. doi: 10.1007/s40314-024-02897-y.

Matieyendou Lamboni, Monod Hervé, and Makowski 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.

Stephen A. Lauer, Kyra H. Grantz, Qifang Bi, Forrest K. Jones, Qulu Zheng, Hannah R. Meredith, Andrew S. Azman, Nicholas G. Reich, and Justin Lessler. The Incubation Period of Coronavirus Disease 2019 (COVID-19) From Publicly Reported Confirmed Cases: Estimation and Application. *Annals of Internal Medicine*, 172(9):577–582, May 2020. ISSN 0003-4819, 1539-3s704. doi: 10.7326/M20-0504.

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.

Isabella Locatelli, Bastien Trächsel, and Valentin Rousson. Estimating the basic reproduction number for covid-19 in western europe. *PLOS ONE*, 16(3):1–9, 03 2021. doi: 10.1371/journal.pone.0248731.

Naoki Masuda and Luis E. C. Rocha. A Gillespie Algorithm for Non-Markovian Stochastic Processes. *SIAM Review*, 60(1):95–115, 2018. doi: 10.1137/16M1055876.

Cameron Nowzari, Masaki Ogura, Victor M. Preciado, and George J. Pappas. A general class of spreading processes with non-Markovian dynamics. In *2015 54th IEEE Conference on Decision and Control (CDC)*, pages 5073–5078, 2015. doi: 10.1109/CDC.2015.7403013.

Matthew J. Penn, Daniel J. Laydon, Joseph Penn, Charles Whittaker, Christian Morgenstern, Oliver Ratmann, Swapnil Mishra, Mikko S. Pakkanen, Christl A. Donnelly, and Samir Bhatt. Intrinsic randomness in epidemic modelling beyond statistical uncertainty. *Communications Physics*, 6(1):146, June 2023. ISSN 2399-3650. doi: 10.1038/s42005-023-01265-2.

Mélanie Prague, Linda Wittkop, Quentin Clairon, Dan Dutartre, Rodolphe Thiébaut, and Boris P. Hejblum. Population modeling of early covid-19 epidemic dynamics in french regions and estimation of the lockdown impact on infection rate. *medRxiv*, 2020. doi: 10.1101/2020.04.21.20073536.

Gesine Reinert. The Asymptotic Evolution of the General Stochastic Epidemic. *The Annals of Applied Probability*, 5(4):1061 – 1086, 1995. doi: 10.1214/aoap/1177004606.
M. Saeedian, M. Khalighi, N. Azimi-Tafreshi, G. R. Jafari, and M. Ausloos. Memory effects on epidemic evolution: The susceptible-infected-recovered epidemic model. *Phys. Rev. E*, 95:022409, Feb 2017. doi: 10.1103/PhysRevE.95.022409.

Charles F Schuler, Carmen Gherasim, Kelly O'Shea, David M Manthei, Jesse Chen, Cristyn Zettel, Jonathan P Troost, Andrew A Kennedy, Andrew W Tai, Donald A Giacherio, Riccardo Valdez, James L Baldwin, and James R Baker. Mild sars-cov-2 illness is not associated with reinfections and provides persistent spike, nucleocapsid, and virus-neutralizing antibodies. *Microbiology spectrum*, 9(2):e0008721, October 2021. ISSN 2165-0497. doi: 10.1128/spectrum.00087-21.

Thomas 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.

Mircea T. Sofonea, Bastien Reyné, Baptiste Elie, Ramsès Djidjou-Demasse, Christian Selinger, Yannis Michalakis, and Samuel Alizon. Memory is key in capturing COVID-19 epidemiological dynamics. *Epidemics*, 35:100459, 2021. ISSN 1755-4365. doi: https://doi.org/10.1016/j.epidem.2021.100459.

G. Streftaris and G. J. Gibson. Non-exponential tolerance to infection in epidemic systems—modeling, inference, and assessment. *Biostatistics*, 13(4):580–593, September 2012. ISSN 1465-4644, 1468-4357. doi: 10.1093/biostatistics/kxs011.

P. Van Mieghem and R. van de Bovenkamp. Non-Markovian Infection Spread Dramatically Alters the Susceptible-Infected-Susceptible Epidemic Threshold in Networks. *Phys. Rev. Lett.*, 110:108701, Mar 2013. doi: 10.1103/Phys-RevLett.110.108701.

Christian L. Vestergaard and Mathieu Génois. Temporal Gillespie Algorithm: Fast Simulation of Contagion Processes on Time-Varying Networks. *PLOS Computational Biology*, 11(10):1–28, 10 2015. doi: 10.1371/journal.pcbi.1004579.
A Proof of Proposition 1

The process {$W(t), t \ge 0$} is clearly a continuous-time homogeneous Markov
process. We will now prove that its transition rates are given by formulae provided
in Table 3. Let $t \ge 0$ and denote by $w$ the current state of $W$ with $\mathbb{P}(W(t) = w) >
0$. Let $\Delta > 0$. In the following, $\Delta$ will tend to zero.

**Infection rate:** The probability of infection between $t$ and $t + \Delta$ is:

$$
\begin{align*}
& \mathbb{P} (W(t + \Delta) - W(t) = \mathbf{u}_{(S,E)} | W(t) = w) \\
&= \sum_{i:X_i(t)=S} \mathbb{P} (Q_{i,j+1} + P(\eta_{i,j}^S) \le P(t+\Delta) | Q_{i,j+1} + P(\eta_{i,j}^S) > P(t)) \\
&= \sum_{i:X_i(t)=S} \left( 1 - \mathbb{P} \left( Q_{i,j+1} + P(\eta_{i,j}^S) > P(t) + \frac{\beta}{N} (W_{I_1}(t) + W_{I_2}(t)) \Delta + o(\Delta) \mid Q_{i,j+1} + P(\eta_{i,j}^S) > P(t) \right) \right) \\
&= \sum_{i:X_i(t)=S} \left( 1 - \exp \left( -\frac{\beta}{N} (W_{I_1}(t) + W_{I_2}(t)) \Delta + o(\Delta) \right) \right) \\
&= W_S(t) \left( 1 - \exp \left( -\frac{\beta}{N} (W_{I_1}(t) + W_{I_2}(t)) \Delta + o(\Delta) \right) \right).
\end{align*}
$$

$$
\text{As } \Delta \to 0, W_S(t) \left(1 - \exp\left(-\frac{\beta}{N}(W_{I_1}(t) + W_{I_2}(t))\Delta + o(\Delta)\right)\right) / \Delta \to \frac{\beta}{N} W_S(t)(W_{I_1}(t) + W_{I_2}(t)) = \frac{\beta}{N} W_S(t)W_I(t).
$$

**Transitions E → I₁:** The probability $\mathbb{P}(W(t + \Delta) - W(t) = \mathbf{u}_{(E,I_1)} | W(t) = w)$ computes as

$$
\begin{align*}
& \mathbb{P} (W(t + \Delta) - W(t) = \mathbf{u}_{(E,I_1)} | W(t) = w) \\
&= \sum_{i:X_i(t)=E} \mathbb{P} \left( L_{i,j+1}^E + \eta_{i,j}^E \le t + \Delta, M_{i,j+1} = (E, I_1) \mid L_{i,j+1}^E + \eta_{i,j}^E > t \right) \\
&= \sum_{i:X_i(t)=E} \mathbb{P} (M_{i,j+1} = (E, I_1)) \left( 1 - \mathbb{P} \left( L_{i,j+1}^E + \eta_{i,j}^E > t + \Delta \mid L_{i,j+1}^E + \eta_{i,j}^E > t \right) \right) \\
&= \sum_{i:X_i(t)=E} p \left( 1 - \exp(-\mu_E \Delta) \right) = pW_E(t) \left( 1 - \exp(-\mu_E \Delta) \right).
\end{align*}
$$

As $\Delta \to 0$, $pW_E(t)(1 - \exp(-\mu_E \Delta)) / \Delta \to p\mu_E W_E(t)$. Computations for transitions $\mathbf{E} \to \mathbf{I}_2$ are similar.
**Other transitions ($\alpha_1, \alpha_2$)** We have:

$$ 
\begin{aligned}
& \mathbb{P}(W(t + \Delta) - W(t) = \mathbf{u}_{(\alpha_1, \alpha_2)} | W(t) = w) \\
&= \sum_{i:X_i(t)=\alpha_1} \mathbb{P} \left( L_{i,j+1}^{\alpha_1} + \eta_{i,j}^{\alpha_1} \le t + \Delta \mid L_{i,j+1}^{\alpha_1} + \eta_{i,j}^{\alpha_1} > t \right) \\
&= \sum_{i:X_i(t)=\alpha_1} \left( 1 - \mathbb{P} \left( L_{i,j+1}^{\alpha_1} + \eta_{i,j}^{\alpha_1} > t + \Delta \mid L_{i,j+1}^{\alpha_1} + \eta_{i,j}^{\alpha_1} > t \right) \right) \\
&= \sum_{i:X_i(t)=\alpha_1} (1 - \exp(-\mu_{\alpha_1} \Delta)) = W_{\mu_1}(t) (1 - \exp(-\mu_{\alpha_1} \Delta)),
\end{aligned}
$$

$$ \text{with } \mu_{\alpha_1} = \begin{cases} \mu_1 & \text{if } \alpha_1 = I_1, \\ \mu_2 & \text{if } \alpha_1 = I_2, \\ \delta & \text{if } \alpha_1 = R. \end{cases} $$

## B Proof of (10) in Section 3.3

Let us define $L_i^S = (L_{i,0}^S, L_{i,1}^S, \dots)$ the successive sojourn durations of individual labeled $i$ in compartment $S$, with $L_{i,0}^S = 0$ except if $X_i(0) = S$. Then for $\alpha \in \{S, E, I, R\}, i = 1, \dots, N, j \ge 0$, we have:

$$ \tau_{i,j+1}^{\alpha} = \tau_{i,j}^{\alpha} + \sum_{\gamma > \alpha} L_{i,j}^{\gamma} + \sum_{\gamma \le \alpha} L_{i,j+1}^{\gamma}. \qquad (15) $$

Applying (15) to $\alpha = S$ we get:

$$ \tau_{i,j+1}^{S} - \tau_{i,j}^{S} = \sum_{\gamma > S} L_{i,j}^{\gamma} + L_{i,j+1}^{S}. \qquad (16) $$

Now, to prove (10), we combine (15) and (16).
C Plots

<figure><img src="image_12.png" /><figcaption>Figure C.1: 50 independent trajectories over time period [0, 500] of number of individuals in compartments *S*, *I*₁, *I*₂ and *R* in Markovian framework (a, b, c) and non-Markovian framework (d, e, f).</figcaption></figure>