SIR epidemics driven by Feller processes

Abstract

We consider a stochastic SIR (susceptible $\rightarrow$ infective $\rightarrow$ removed) model in which the infectious periods are modulated by a collection of independent and identically distributed Feller processes. Each infected individual is associated with one of these processes, the trajectories of which determine the duration of his infectious period, his contamination rate, and his type of removal (e.g. death or immunization). We use a martingale approach to derive the distribution of the final epidemic size and severity for this model and provide some general examples. Next, we focus on a single infected individual facing a given number of susceptibles, and we determine the distribution of his outcome (number of contaminations, severity, type of removal). Using a discrete-time formulation of the model, we show that this distribution also provides us with an alternative, more stable method to compute the final epidemic outcome distribution.

1. Introduction

Compartmental epidemic models are a popular way to describe the spread of an infectious disease in a population. Among them, one of the most popular is the SIR (susceptible $\rightarrow$ infective $\rightarrow$ removed) model, in which the infected individuals remain contagious for a specific period of time (the infectious period) before being permanently removed from the set of potential vectors of the infection. This typical structure is widely used in the literature, for both deterministic and stochastic models. We refer to the books by Daley and Gani [11] and Andersson and Britton [1] for a good introduction to the mathematical analysis of epidemic processes.

In particular, the so-called general epidemic model has received much attention in the literature. In this classical stochastic SIR process, the state of the population is given by the vector (S(t), I(t), R(t)) of the number of susceptibles, infected, and removed individuals present at time t. The individuals behave independently of each other and the infectious periods are exponentially distributed. Moreover, there can be only one contamination at a time, at a rate $\beta(S(t),I(t)) = \beta S(t) I(t)$ for some $\beta >0$ . This model has been studied by means of different methods including the use of martingale arguments (see e.g. Picard [20, 21] and Picard and Lefèvre [22]) and probabilistic reasoning (see e.g. Ball [3] and Gani and Purdue [13]). Various approximation methods have also been developed, for instance in Ball and Donnelly [5] and Barbour and Reinert [7].

The general epidemic has been extended to various more flexible and realistic models by means of different methods. A first possibility of extension is to make the contamination rate $\beta(s,i)$ a more general function of the number of individuals present in each class, for instance in Neuts and Li [18], Clancy [9, 10], and O’Neill [19]. Another possible extension is to relax the assumption of exponentially distributed infectious periods, by allowing transitions between different states during the infection. This was done in Black and Ross [8], Picard and Lefèvre [22], and Simon [24], for example. Some authors also used block-structured Markov chains to describe the epidemic, in order to allow for random environments and state-dependent contamination rates (see e.g. Artalejo et al. [2] and Lefèvre and Simon [16]). Finally, the homogeneity assumption can also be relaxed. It was done for instance in Ball et al. [6], where the authors considered two levels of mixing (local and global) in the population.

In this paper we consider a stochastic SIR model in which the infectious periods are driven by independent, identically distributed Feller processes. These processes are used to represent the evolution of the infectivity and removal rates of the individuals throughout their infectious period. More precisely, one of these processes is activated when an individual gets infected. His contamination rate varies over the infectious period: it is of the form $\beta(s,x) = \beta_s \gamma(x)$ , where $\beta_s$ is an arbitrary function of the number of available susceptibles and x is the state occupied by his associated Feller process. Moreover, the removal of an individual can take two different forms (e.g. death or immunization) and the removal rates are also a function of x, so the duration of the infectious periods and the type of removal depend on the trajectory of the process. The wide class of Feller processes makes this model quite general and gives great flexibility to include different characteristics of the disease and of the infection mechanism in the model, while the general form $\beta_s$ makes it possible to include a large choice of dependences between the number of susceptibles and the contamination rate of the infected individuals. This flexibility is illustrated by different general examples at the end of Section 3.

In Section 3 we build a family of martingales involving different statistics related to the final epidemic outcome: the final number of susceptibles, the total severity, and the final numbers of removals of both types. We then show how to exploit these martingales to determine the joint distribution of the statistics and we provide a few examples.

Section 4 is devoted to the presentation of an alternative method for the computation of the final epidemic outcome distribution. This is motivated by the fact that the numerical procedure obtained through this martingale approach often becomes unstable when the population size is too big, which can happen for relatively small populations (see e.g. Demiris and O’Neill [12] and House et al. [15]). A first step towards this goal is taken in Section 4.1. Here, we focus on a single, generic infected individual facing a given number of susceptibles and we determine the distribution of his personal outcome: the distribution of his own number of contaminations, his contribution to the final severity, and his type of removal. This study also provides a better understanding of the transmission mechanism of the disease and can be useful, e.g. for parameter estimation purposes.

Next, we establish a link between the total epidemic outcome and the personal outcome of a generic infected individual. This is done by building a discrete-time epidemic model, the final state of which is distributed as the final state of the original model. This construction has already been used in the setting of the general epidemic with linear contamination rates $\beta_s=\beta s$ (see e.g. Lefèvre and Utev [17] and Ball [3]). In Section 4.2 we show that such a link can still be established when the contamination rates $\beta_s$ are arbitrary and when the infectious periods are represented by Feller processes. Finally, we show the practical interest of this link by using the equivalent discrete-time model to obtain an alternative, more stable procedure for the computation of the final outcome distribution. We conclude with a small illustration showing that this procedure is numerically stable for much larger population sizes.

2. Model

Consider a closed population of $n+m$ individuals ( $n,m \in \mathbb{N}_0$ ) and a Feller process $X = \{X(u) \mid u \in \mathbb{R}^+\}$ taking values in a metric space $\mathcal{E}$ . Let $\{\Phi_u\}$ denote the semigroup of transition operators associated with X, such that

\begin{equation*}\Phi_u\,f(x) = \mathbb{E}_{x}[ f(X(u)) ] \equiv \mathbb{E}[ f(X(u))\mid X(0)=x ]\end{equation*}

for all f in the set $B(\mathcal{E})$ of bounded and measurable functions from $\mathcal{E}$ to $\mathbb{R}$ . The dynamics of the process X is characterized by its (extended) generator $(\mathcal{L},\mathcal{D}_\mathcal{L})$ , where

\begin{equation*}\mathcal{L}\;:\; \mathcal{D}_\mathcal{L} \subseteq B(\mathcal{E}) \to B(\mathcal{E})\end{equation*}

is given by

\begin{equation*}\mathcal{L} f = \lim_{u \to 0^+} \dfrac{1}{u}(\Phi_u\,f -f),\end{equation*}

where the limit is in the sense of the supremum norm. The set $\mathcal{D}_\mathcal{L}$ contains all the functions in $B(\mathcal{E})$ for which this limit exists. Finally, the initial state X(0) is a random variable with distribution given by the probability measure $\alpha(\!\cdot\!)$ on $\mathcal{E}$ . Throughout the paper we refer to X as the ‘infectivity process’.

The SIR epidemics considered in this paper have the following dynamics. We assume that when an individual (labelled i, say) gets infected, his infectious period is governed by a process $X_i$ distributed as X. As long as he is not removed, the infected individual can contaminate the susceptibles one by one. His contamination rate depends on the number s of available susceptibles and of the state x occupied by $X_i$ . It is of the form $\beta_s \gamma(x)$ , where $\beta_s >0$ , $\beta_s \neq \beta_u$ for $s \neq u$ and $\gamma\;:\; \mathcal{E} \to \mathbb{R}^+$ is a càdlàg function. Most traditional models assume that $\beta_s = \beta s$ for some constant $\beta>0$ . The form $\beta_s$ gives more flexibility in the dependence between the number of susceptibles and the contamination rates. For instance, taking $\beta_s = \beta s^{\lambda}$ for some $\lambda >0$ allows us to incorporate more group effects in the model, as the infection rate no longer varies linearly with the number of susceptibles.

To assess the severity and the impact of the epidemic, it can be useful to distinguish between different kinds of elimination, for instance to count separately the individuals who died from the disease and those who recovered. For that reason, we assume that each infected individual can become a removed of either type 1 or type 2 at the end of his infectious period. The removal rate and the type of removal are also functions of his associated infectivity process:

\begin{equation*} \begin{aligned}& \mathbb{P}( \tau_{\star} \in \!]u, u+{\textrm{d}} u[,\, \mathbb{I}_{1}=1 \mid X_i(v), 0\le v \le u, \tau_{\star}>u ) = h_1(X_i(u))\, {\textrm{d}} u + {\textrm{o}}({\textrm{d}} u), \\[5pt] & \mathbb{P}( \tau_{\star} \in \!]u, u+{\textrm{d}} u[, \, \mathbb{I}_{2}=1 \mid X_i(v), 0\le v \le u, \tau_{\star}>u ) = h_2(X_i(u))\, {\textrm{d}} u + {\textrm{o}}({\textrm{d}} u),\end{aligned}\end{equation*}

where $\tau_{\star}$ is the removal time (when $u=0$ corresponds to the beginning of the infectious period), $\mathbb{I}_{1}, \mathbb{I}_{2}$ are the indicators that he becomes a removed of types 1 and 2, respectively, and $h_1$ , $h_2$ are two positive and càdlàg functions. We assume that $h_1, h_2$ are chosen so that all infectious periods have an almost surely finite duration.

Finally, we assume that the infected individuals behave independently of each other, that is, the various versions $X_i$ that activate throughout the epidemic are independent of each other.

Let $S(t)$ , $I(t)$ , $R_1(t)$ , and $R_2(t)$ denote the number of susceptibles, infected individuals, and removed cases of types 1 and 2 at time t. The epidemic starts at time $t=0$ with $S(0)=n$ and $I(0)=m$ . The number of susceptibles can only decrease over time, so $\{S(t)\}$ takes its values in $\{0,1,\ldots,n\}$ . The processes $\{I(t)\}$ , $\{R_1(t)\}$ , and $\{R_2(t)\}$ take their values in $\{0,1,\ldots,n+m\}$ . As the population is closed, $S(t) + I(t) + R_1(t) + R_2(t) = n + m$ for all $t \ge 0$ . The epidemic terminates at the first time

\begin{equation*}T = \inf\{t\ge 0 \mid I(t)=0\}\end{equation*}

when there are no more infected individuals. Note that $T<\infty$ almost surely since the population size is finite and all infectious periods are almost surely finite.

Below, we are interested in the joint distribution of four statistics related to the final epidemic outcome for this class of SIR epidemics. The first three statistics are the numbers $S(T)$ , $R_1(T)$ , and $R_2(T)$ of susceptibles and removed individuals at time T. The last statistic is the final severity $A(T)$ , defined as the total cost of the epidemic when an infected individual whose infectivity process is in state x has an instantaneous cost a(x) for some fixed càdlàg function $a\;:\; \mathcal{E} \to \mathbb{R}^+$ . More precisely, the severity up to time t is given by

(2.1) \begin{equation}A(t) = \int_0^t \sum_{i=1}^{I(u)} a(Y_i(u)) \, {\textrm{d}} u,\end{equation}

where $Y_i(t)$ denotes the state occupied at time t by the infectivity process associated to the ith infected individual present at time t: if $\xi_i$ is the infection time of this individual, then $Y_i(\xi_i) = X_i(0) \sim \alpha(\!\cdot\!)$ and

\begin{equation*}Y_i(t) = X_i(t-\xi_i).\end{equation*}

Note that only individuals truly contagious at time t are considered in the sum in (2.1). As an example, if $a(x)=1$ for all $x \in \mathcal{E}$ , then $A(T)$ gives the cumulative duration of all infectious periods up to the ending time of the epidemic.

3. Final epidemic outcome

To determine the joint distribution of the four statistics above, we build a family of martingales involving $S(t)$ , $I(t)$ , $R_1(t)$ , $R_2(t)$ , and $A(t)$ , and we make use of the optional stopping theorem. Let $(\mathcal{F}_t)_{t \ge 0}$ denote the filtration generated by the history of the epidemic process up to time t. We start with the following result.

Proposition 3.1. Fix $\theta \ge 0$ and $z_1,z_2 \in\!]0,1]$ . If $f \in \mathcal{D}_\mathcal{L}$ and $c\;:\; \{0,1,\ldots,n\} \to \mathbb{R}^+$ satisfy

(3.1) \begin{align}& c(s) \mathcal{L} f(x) + [ c(s-1) \beta_s \gamma(x) - c(s) (\theta a(x)+ \beta_s \gamma(x) + h_1(x) + h_2(x) ) ]\, f(x) \notag \\[5pt] & \quad + c(s) (z_1h_1(x) + z_2h_2(x)) = 0\end{align}

for all $x \in \mathcal{E}$ and $s \in \{0,1,\ldots,n\}$ , then the process

\begin{equation*} M(t) = c(S(t)) \biggl(\int_{\mathcal{E}} f(y) \, {\textrm{d}}\alpha(y)\biggr)^{S(t)}\, {\textrm{e}}^{-\theta A(t)} \, z_1^{R_1(t)} \, z_2^{R_2(t)}\, \prod_{i=1}^{I(t)} f(Y_i(t))\end{equation*}

is a martingale with respect to $(\mathcal{F}_t)_{t \ge 0}$ .

Proof. We look for some coefficients b(s) and some bounded and measurable function $f\;:\; \mathcal{E} \to ]0,1]$ such that the process

\begin{equation*}M(t) = b(S(t)) \, {\textrm{e}}^{-\theta A(t)} \, z_1^{R_1(t)} \, z_2^{R_2(t)}\, \prod_{i=1}^{I(t)} f(Y_i(t))\end{equation*}

is a martingale with respect to $(\mathcal{F}_t)_{t \ge 0}$ . For that, we fix $t_0\ge 0$ and define for all $t \ge t_0$

\begin{equation*}m(t) = \mathbb{E}\bigl[ M(t) \mid \mathcal{F}_{t_0} \bigr].\end{equation*}

If $b(\!\cdot\!)$ and $f(\!\cdot\!)$ are chosen in such a way that the right derivative $m_r^{\prime}(t)$ of m(t) equals zero for all $t \ge t_0$ , then M is a martingale. To differentiate m(t), we first use the fact that $\mathcal{F}_{t_0} \subseteq \mathcal{F}_{t}$ and the tower property to write

(3.2) \begin{equation}m(t+{\textrm{d}} t) = \mathbb{E}\bigl[ M(t+{\textrm{d}} t) \mid \mathcal{F}_{t_0} \bigr]= \mathbb{E}\bigl[ \mathbb{E}[ M(t+{\textrm{d}} t) \mid \mathcal{F}_{t} ] \mid \mathcal{F}_{t_0} \bigr]. \end{equation}

Next, we examine the possible events on the time interval $[t, t+{\textrm{d}} t[$ : there may be a contamination during this interval, or the elimination of an infected individual, or neither of the two. The events involving the contamination/removal of more than one individual on $[t, t+{\textrm{d}} t[$ occur with probability ${\textrm{o}}({\textrm{d}} t)$ , so we obtain

\begin{align*}& \mathbb{E}[ M(t+{\textrm{d}} t) \mid \mathcal{F}_{t} ] \\[5pt] &\quad = b(S(t) - 1) \, z_1^{R_1(t)} \, z_2^{R_2(t)} \, \mathbb{E}\Biggl[{\textrm{e}}^{-\theta A(t+{\textrm{d}} t)}\prod_{i=1}^{I(t)} f(Y_i(t+{\textrm{d}} t)) \biggm| \mathcal{F}_{t} \Biggr]\biggl(\int_{\mathcal{E}} f(y) \, {\textrm{d}}\alpha(y)\biggr)\\[5pt] &\quad\quad\quad \times \Biggl(\sum_{i=1}^{I(t)} \beta_{S(t)}\gamma(Y_i(t)) \, {\textrm{d}} t + {\textrm{o}}({\textrm{d}} t) \Biggr) \\[5pt] &\quad\quad + \sum_{i=1}^{I(t)} b(S(t)) \, z_1^{R_1(t)+1} \, z_2^{R_2(t)} \,\mathbb{E}\!\left[\rule{0pt}{24pt}\right.\! {\textrm{e}}^{-\theta A(t+{\textrm{d}} t)}\prod_{\substack{v=1 \\ v\neq i}}^{I(t)}f(Y_v(t+{\textrm{d}} t)) \Biggm| \mathcal{F}_{t} \!\left.\rule{0pt}{24pt}\right]\! ( h_1(Y_i(t)) \, {\textrm{d}} t + {\textrm{o}}({\textrm{d}} t) ) \\[5pt] &\quad\quad + \sum_{i=1}^{I(t)} b(S(t)) \, z_1^{R_1(t)} \, z_2^{R_2(t)+1} \,\mathbb{E} \!\left[\rule{0pt}{24pt}\right.\! {\textrm{e}}^{-\theta A(t+{\textrm{d}} t)}\prod_{\substack{v=1 \\ v\neq i}}^{I(t)}f(Y_v(t+{\textrm{d}} t)) \Biggm| \mathcal{F}_{t} \!\left.\rule{0pt}{24pt}\right]\! ( h_2(Y_i(t)) \, {\textrm{d}} t + {\textrm{o}}({\textrm{d}} t) ) \\[5pt] &\quad\quad + b(S(t)) \, z_1^{R_1(t)} \, z_2^{R_2(t)} \,\mathbb{E}\Biggl[ {\textrm{e}}^{-\theta A(t+{\textrm{d}} t)}\prod_{i=1}^{I(t)}f(Y_v(t+{\textrm{d}} t)) \biggm| \mathcal{F}_{t} \Biggr] \\[5pt] &\quad\quad \quad \times \Biggl(1-\sum_{i=1}^{I(t)}(\beta_{S(t)}\gamma(Y_i(t)) + h_1(Y_i(t)) + h_2(Y_i(t))) \,{\textrm{d}} t + {\textrm{o}}({\textrm{d}} t) \Biggr) \\[5pt] &\quad\quad + {\textrm{o}}({\textrm{d}} t).\end{align*}

Taking $b(\!\cdot\!)$ of the form $b(s) = c(s) \lambda^s$ with $\lambda = \int_{\mathcal{E}} f(y) \, {\textrm{d}}\alpha(y)$ and using the notation

\[P(t,{\textrm{d}} t)= {\textrm{e}}^{-\theta A(t+{\textrm{d}} t)}\prod_{i=1}^{I(t)} f(Y_i(t+{\textrm{d}} t)),\]

this equality becomes

\begin{align*}& \mathbb{E}[ M(t+{\textrm{d}} t) \mid \mathcal{F}_{t} ] \\[5pt] &\quad = \sum_{i=1}^{I(t)} c(S(t) - 1) \, \lambda^{S(t)}\, z_1^{R_1(t)} \, z_2^{R_2(t)} \, \mathbb{E} [P(t,{\textrm{d}} t) \mid \mathcal{F}_{t} ]\beta_{S(t)}\gamma(Y_i(t)) \, {\textrm{d}} t \\[5pt] &\quad\quad + \sum_{i=1}^{I(t)} c(S(t)) \, \lambda^{S(t)} \, z_1^{R_1(t)} \, z_2^{R_2(t)} \,\mathbb{E} [P(t,{\textrm{d}} t) \mid \mathcal{F}_{t} ]\, \dfrac{z_1 h_1(Y_i(t)) + z_2 h_2(Y_i(t))}{f(Y_i(t+{\textrm{d}} t))} \, {\textrm{d}} t \\[5pt] &\quad\quad - \sum_{i=1}^{I(t)} c(S(t)) \, \lambda^{S(t)}\, z_1^{R_1(t)} \, z_2^{R_2(t)} \,\mathbb{E} [P(t,{\textrm{d}} t) \mid \mathcal{F}_{t} ](\beta_{S(t)}\gamma(Y_i(t)) + h_1(Y_i(t)) + h_2(Y_i(t))) \,{\textrm{d}} t \\[5pt] &\quad\quad + c(S(t)) \, \lambda^{S(t)}\, z_1^{R_1(t)} \, z_2^{R_2(t)} \,\mathbb{E}[P(t,{\textrm{d}} t) \mid \mathcal{F}_{t} ] + {\textrm{o}}({\textrm{d}} t).\end{align*}

Now, subtract m(t) from both sides of (3.2), divide each side by ${\textrm{d}} t$ , and then take ${\textrm{d}} t \to 0$ . Using the last equality and the dominated convergence theorem, this yields

(3.3) \begin{align}m_r^{\prime}(t)&= \mathbb{E}\Biggl[\sum_{i=1}^{I(t)} \Pi(t) \biggl( c(S(t) - 1) \,\beta_{S(t)}\gamma(Y_i(t))+ c(S(t)) \, \dfrac{z_1 h_1(Y_i(t)) + z_2 h_2(Y_i(t))}{f(Y_i(t))} \nonumber\\[5pt] &\quad \quad - c(S(t)) (\beta_{S(t)}\gamma(Y_i(t)) + h_1(Y_i(t)) + h_2(Y_i(t))) \biggr) \nonumber \\[5pt] & \quad \quad + c(S(t)) \, \lambda^{S(t)}\, z_1^{R_1(t)} \, z_2^{R_2(t)} \,\lim_{{\textrm{d}} t\to 0^+} \dfrac{1}{{\textrm{d}} t}\mathbb{E} [P(t,{\textrm{d}} t)-P(t)\mid \mathcal{F}_{t} ] \biggm| \mathcal{F}_{t_0}\Biggr], \end{align}

where

\begin{align*} P(t)= {\textrm{e}}^{-\theta A(t)}\prod_{i=1}^{I(t)} f(Y_i(t)),\quad\Pi(t) = \lambda^{S(t)} z_1^{R_1(t)} z_2^{R_2(t)} P(t).\end{align*}

To determine the limit in (3.3), we note from (2.1) that

\begin{align*}{\textrm{e}}^{-\theta A(t+{\textrm{d}} t)} &= {\textrm{e}}^{-\theta A(t) -\theta\sum_{i=1}^{I(t)} a(Y_i(t)) \, {\textrm{d}} t + \omega(t,{\textrm{d}} t)} \\[5pt] &= {\textrm{e}}^{-\theta A(t)} \Biggl(1-\theta \sum_{i=1}^{I(t)} a(Y_i(t)) \, {\textrm{d}} t + \omega(t,{\textrm{d}} t)\Biggr),\end{align*}

where $\omega(t,{\textrm{d}} t)$ denotes any random variable such that $\lim_{{\textrm{d}} t \to 0}\mathbb{E}[ \omega(t,{\textrm{d}} t)\mid \mathcal{F}_{t}]/{\textrm{d}} t=0$ almost surely. Thus

\begin{align*}& \lim_{{\textrm{d}} t\to 0^+} \dfrac{1}{{\textrm{d}} t}\mathbb{E}[P(t,{\textrm{d}} t)-P(t)\mid \mathcal{F}_{t} ] \\[5pt] &\quad = {\textrm{e}}^{-\theta A(t)} \lim_{{\textrm{d}} t\to 0^+} \dfrac{1}{{\textrm{d}} t}\mathbb{E}\Biggl[\prod_{i=1}^{I(t)} f(Y_i(t+{\textrm{d}} t)) - \prod_{i=1}^{I(t)} f(Y_i(t)) \biggm| \mathcal{F}_{t} \Biggr] - P(t) \sum_{i=1}^{I(t)} \theta a(Y_i(t)).\end{align*}

Using the fact that the processes $Y_i$ are independent of each other and their generator is $\mathcal{L}$ , we obtain

\begin{align*}\lim_{{\textrm{d}} t\to 0^+} \dfrac{1}{{\textrm{d}} t}\mathbb{E}[P(t,{\textrm{d}} t)-P(t)\mid \mathcal{F}_{t} ]= P(t) \sum_{i=1}^{I(t)}\dfrac{\mathcal{L} f(Y_i(t))}{f(Y_i(t))} - P(t) \sum_{i=1}^{I(t)} \theta a(Y_i(t)).\end{align*}

Substituting the last equality in (3.3) yields

(3.4) \begin{align} m_r^{\prime}(t)& = \mathbb{E}\Biggl[\Pi(t) \sum_{i=1}^{I(t)} \biggl[ c(S(t)- 1) \,\beta_{S(t)}\gamma(Y_i(t)) + c(S(t)) \, \dfrac{z_1 h_1(Y_i(t)) + z_2 h_2(Y_i(t))}{f(Y_i(t))} \nonumber\\[5pt] &\quad + c(S(t))\biggl(\dfrac{\mathcal{L} f(Y_i(t))}{f(Y_i(t))} - \beta_{S(t)}\gamma(Y_i(t)) - \theta a(Y_i(t)) - h_1(Y_i(t)) -h_2(Y_i(t)) \biggr)\biggr] \biggm| \mathcal{F}_{t_0}\Biggr] .\end{align}

If $c(\!\cdot\!)$ and $f(\!\cdot\!)$ satisfy (3.1), then each term inside the sum in (3.4) is equal to zero for all possible values of $S(t)$ and $Y_i(t)$ . It implies that $m_r^{\prime}(t) = 0$ for all $ t \ge t_0$ , and therefore that M is a martingale.

We now derive different solutions of equation (3.1), which provide us with a family of martingales. Henceforth and by convention, a product $ \prod_{j=r}^{r-k} c_j$ is equal to one if $k=1$ and is equal to zero if $k>1$ .

Theorem 3.1. For all $k \in \{0,1,2,\ldots,n\}$ , $\theta \ge 0$ , and $z_1,z_2 \in\!]0,1]$ , the process

(3.5) \begin{equation}\Biggl\{ \Biggl(\prod_{l=k+1}^{S(t)}\dfrac{\beta_l}{\beta_l-\beta_k} \Biggr) q_{\alpha}(k,\theta,z_1,z_2)^{S(t)} \, {\textrm{e}}^{-\theta A(t)} \, z_1^{R_1(t)} \, z_2^{R_2(t)}\, \prod_{i=1}^{I(t)} q(k,\theta,z_1,z_2,Y_i(t)) \Biggr\} \end{equation}

is a martingale with respect to $(\mathcal{F}_t)_{t \ge 0}$ when $q(k,\theta,z_1,z_2,x)$ is the solution in $\mathcal{D}_\mathcal{L}$ of the equation

(3.6) \begin{equation}\mathcal{L} f(x) = (\theta a(x)+ \beta_k \gamma(x) + h_1(x) +h_2(x) )\,f(x) - (z_1h_1(x)+z_2h_2(x)),\end{equation}

and when

(3.7) \begin{equation}q_{\alpha}(k,\theta,z_1,z_2) = \int_{\mathcal{E}} q(k,\theta,z_1,z_2,y) \, {\textrm{d}}\alpha(y).\end{equation}

Proof. Equation (3.1) can be rewritten as

(3.8) \begin{align}& c(s) [\mathcal{L} f(x) - (\theta a(x)+ \beta_s \gamma(x) + h_1(x) + h_2(x) )\,f(x) + z_1h_1(x) + z_2h_2(x)] \notag \\[5pt] & \quad = - c(s-1) \beta_s \gamma(x) f(x) \quad \text{{for all $x \in \mathcal{E}$, $ s \in \{0,1,\ldots,n\}$.}}\end{align}

Fix $k \in \{0,1,2,\ldots,n\}$ and choose the function $f\equiv f_k$ so that the bracket on the left-hand side of (3.8) is always equal to zero for $s=k$ , that is, $f_k$ is a solution of (3.6). With this choice, for $s=k$ , equation (3.8) becomes

\begin{equation*}0 = c(k-1) \beta_k \gamma(x) f_k(x) \quad \text{{for all $x \in \mathcal{E}$,}}\end{equation*}

which implies that $c(k-1)=0$ and therefore $c(s)=0$ for all $s<k$ . For $s>k$ , (3.8) can be rewritten as

\begin{align*}& c(s) [(\mathcal{L} f_k(x) - (\theta a(x)+ \beta_k \gamma(x) + h_1(x) + h_2(x))\,f_k(x) + z_1h_1(x) + z_2h_2(x) )] \\[5pt] & \quad = c(s) (\beta_s-\beta_k) \gamma(x) f_k(x) - c(s-1) \beta_s \gamma(x) f_k(x) \quad \text{{for all $x \in \mathcal{E}$, $ s \in \{0,1,\ldots,n\}$.}}\end{align*}

Since $f_k$ satisfies (3.6), the left-hand side of this equality is equal to zero and we obtain

\begin{equation*}c(s) (\beta_s-\beta_k) = c(s-1) \beta_s, \quad {s>k.}\end{equation*}

Choosing $c(k)=1$ , this recurrence is easily solved and yields

\begin{equation*} c(s)=\prod_{l=k+1}^{s}\beta_l/(\beta_l-\beta_k). \end{equation*}

A solution of (3.6) can always be obtained as follows: let $\tau_{\star}$ be the duration of the infectious period for an individual with infectivity process X. Then we have the following result.

Proposition 3.2. The function

(3.9) \begin{equation}q(k,\theta,z_1,z_2,x) = \mathbb{E}_{x}\Bigl[ z_1^{\mathbb{I}_{1}} \, z_2^{\mathbb{I}_{2}} \,{\textrm{e}}^{-\int_0^{\tau_{\star}} (\theta a(X(u)) + \beta_k \gamma(X(u)) ) \, {\textrm{d}} u} \Bigr] \end{equation}

is a solution of (3.6) for all $k \in \{0,1,2,\ldots,n\}$ , $\theta \ge 0$ , and $z_1,z_2 \in\!]0,1]$ .

Proof. Start from the expectation in (3.9) (where we assume that the process X starts at time 0 in state x) and condition on whether or not the removal occurs on the time interval [0, t]. For small t, denoting $r(x) = \theta a(x) + \beta_k \gamma(x)$ , we have

\begin{align*}q(k,\theta,z_1,z_2,x) &= \mathbb{E}_x\Bigl[{z_1 \,{\textrm{e}}^{-\int_0^{\tau_{\star}} r(X(u) \, {\textrm{d}} u} \mid \tau_{\star} \le t, \mathbb{I}_{1}=1}\Bigr] (h_1(x)t + {\textrm{o}}(t)) \\[3pt] &\quad + \mathbb{E}_x\Bigl[{z_2 \,{\textrm{e}}^{-\int_0^{\tau_{\star}} r(X(u) \, {\textrm{d}} u} \mid \tau_{\star} \le t, \mathbb{I}_{2}=1}\Bigr] (h_2(x)t + {\textrm{o}}(t)) \\[3pt] &\quad + \mathbb{E}_x\Bigl[{z_1^{\mathbb{I}_{1}} \, z_2^{\mathbb{I}_{2}} \,{\textrm{e}}^{-\int_0^{\tau_{\star}} r(X(u) \, {\textrm{d}} u} \mid \tau_{\star} >t}\Bigr] (1-h_1(x)t -h_2(x)t + {\textrm{o}}(t))\\[3pt] &= (1-r(x)t+{\textrm{o}}(t))(z_1h_1(x)t + z_2h_2(x)t +{\textrm{o}}(t)) \\[3pt] &\quad + \mathbb{E}_x\Bigl[{z_1^{\mathbb{I}_{1}} \, z_2^{\mathbb{I}_{2}} \,{\textrm{e}}^{-\int_0^{\tau_{\star}} r(X(u) \, {\textrm{d}} u} \mid \tau_{\star} >t}\Bigr] (1-h_1(x)t -h_2(x)t + {\textrm{o}}(t)).\end{align*}

Letting $(\mathcal{F}^X_t)_{t \ge 0}$ denote the filtration generated by the history of X up to time t, and using the tower property and the Markov property,

\begin{align*}\mathbb{E}_x\Bigl[{z_1^{\mathbb{I}_{1}} \, z_2^{\mathbb{I}_{2}} \,{\textrm{e}}^{-\int_0^{\tau_{\star}} r(X(u) \, {\textrm{d}} u} \mid \tau_{\star} >t}\Bigr]& = \mathbb{E}_x\Bigl[ \mathbb{E}\Bigr[{z_1^{\mathbb{I}_{1}} \, z_2^{\mathbb{I}_{2}} \,{\textrm{e}}^{-\int_0^{\tau_{\star}} r(X(u) \, {\textrm{d}} u} \mid \mathcal{F}^X_t, \tau_{\star} >t}\Bigr] \Bigr] \\[3pt] & = \mathbb{E}_x\Bigl[ {\textrm{e}}^{-\int_0^{t} r(X(u) \, {\textrm{d}} u} \mathbb{E}_{X(t)}\Bigl[{z_1^{\mathbb{I}_{1}} \, z_2^{\mathbb{I}_{2}} \,{\textrm{e}}^{-\int_0^{\tau_{\star}} r(X(u) \, {\textrm{d}} u} }\Bigr] \Bigr] \\[3pt] & = \mathbb{E}_x [ (1-r(x)t+{\textrm{o}}(t)) \, q(k,\theta,z_1,z_2,X(t)) ].\end{align*}

Therefore

\begin{align*}q(k,\theta,z_1,z_2,x) &= (z_1h_1(x) + z_2h_2(x))t + \mathbb{E}_{x}[ q(k,\theta,z_1,z_2,X(t)) ] \\[5pt] &\quad -\mathbb{E}_{x}[ q(k,\theta,z_1,z_2,X(t)) ] (r(x)+h_1(x)+h_2(x))t + {\textrm{o}}(t).\end{align*}

Dividing both sides of the last equality by t, after some rearranging we obtain

\begin{align*}& \dfrac{1}{t} (\mathbb{E}_{x}[ q(k,\theta,z_1,z_2,X(t)) ] - q(k,\theta,z_1,z_2,x))\\[5pt] &\quad = (r(x)+h_1(x)+h_2(x))\,\mathbb{E}_{x}[ q(k,\theta,z_1,z_2,X(t)) ] - (z_1h_1(x) + z_2h_2(x))+ \dfrac{{\textrm{o}}(t)}{t}.\end{align*}

It then suffices to take the limit $t\to 0^+$ to see that $q(k,\theta,z_1,z_2,x)$ satisfies equation (3.6).

Applying the optional stopping theorem to the martingales (3.5) with T as the stopping time, we obtain a transform of the joint distribution of the final number of susceptibles $S(T)$ , the final severity $A(T)$ , and the final number of removals.

Corollary 3.1. For all $k \in \{0,1,2,\ldots,n\}$ , $\theta \ge 0$ , and $z_1,z_2 \in\!]0,1]$ ,

(3.10) \begin{align}&\mathbb{E}\Biggl[ \Biggl(\prod_{l=k+1}^{S(T)}\dfrac{\beta_l}{\beta_l-\beta_k} \Biggr) q_{\alpha}(k,\theta,z_1,z_2)^{S(T)} \, {\textrm{e}}^{-\theta A(T)} \,z_1^{R_1(T)} \,z_2^{R_2(T)} \Biggr] \notag \\[5pt] & \quad = q_{\alpha}(k,\theta,z_1,z_2)^{n+m} \prod_{l=k+1}^{n}\dfrac{\beta_l}{\beta_l-\beta_k}.\end{align}

From (3.10), we can compute the joint distribution of $S(T)$ , $A(T)$ , $R_1(T)$ , and $R_2(T)$ .

Corollary 3.2. For all $\theta \ge 0$ , the numbers

\begin{equation*}y_{(s,r)}(\theta) = \mathbb{E}\bigl[ {\textrm{e}}^{-\theta A(T)} \, \unicode{x1d7d9}_{S(T)=s, R_1(T)=r} \bigr], \quad {0 \le s \le n,\ 0 \le r \le n+m-s,}\end{equation*}

constitute the solution of the linear system of equations

(3.11) \begin{align}&\sum_{s=k}^{n} \Biggl(\prod_{l=k+1}^{s}\dfrac{\beta_l}{\beta_l-\beta_k} \Biggr) \sum_{j=0}^{r} \left(\substack{s\\[4pt] j}\right) \, q_{\alpha}(k,\theta,1,0)^j \, q_{\alpha}(k,\theta,0,1)^{s-j} \, y_{(s,r-j)}(\theta) \nonumber \\[5pt] & \quad = \left(\substack{n+m\\[4pt] r}\right) \Biggl(\prod_{l=k+1}^{n}\dfrac{\beta_l}{\beta_l-\beta_k}\Biggr) q_{\alpha}(k,\theta,1,0)^{r} \, q_{\alpha}(k,\theta,0,1)^{n+m-r}\end{align}

for $0\le k \le n$ , $0 \le r \le n+m-k$ .

Proof. Taking $z_1\equiv z$ and $z_2=1$ , we can rewrite (3.10) as

\begin{align*}& \sum_{s=k}^{n} \Biggl(\prod_{l=k+1}^{s}\dfrac{\beta_l}{\beta_l-\beta_k} \Biggr) q_{\alpha}(k,\theta,z,1)^{s} \, \mathbb{E}\bigl[ {\textrm{e}}^{-\theta A(T)} \,z^{R_1(T)} \,\unicode{x1d7d9}_{S(T)=s} \bigr]\\[5pt] & \quad = q_{\alpha}(k,\theta,z,1)^{n+m} \prod_{l=k+1}^{n}\dfrac{\beta_l}{\beta_l-\beta_k}.\end{align*}

Differentiating r times with respect to z and taking $z\to 0^+$ , we obtain (3.11) after having noticed that

\begin{align*} &\dfrac{{\textrm{d}}^j}{{\textrm{d}} z^j}\mathbb{E}\bigl[ {\textrm{e}}^{-\theta A(T)} \,z^{R_1(T)} \,\unicode{x1d7d9}_{S(T)=s} \bigr]\big|_{z=0} = j!\,y_{(s,j)}(\theta), \\[5pt] & \dfrac{{\textrm{d}}^j}{{\textrm{d}} z^j}q_{\alpha}(k,\theta,z,1)^{s}\big|_{z=0^+} = \dfrac{s!}{(s-j)!} q_{\alpha}(k,\theta,1,0)^j \, q_{\alpha}(k,\theta,0,1)^{s-j},\end{align*}

the last equality following from (3.6) or (3.9).

For each value of $\theta \ge 0$ , the equalities in (3.11) constitute a system of

\[\sum_{k=0}^n (n + m - k + 1) = (n + 1)\biggl(\dfrac{n}{2} + m + 1\biggr)\]

linear equations to determine the numbers $y_{(s,r)}(\theta)$ . In particular, the probabilities $\mathbb{P}( S(T)=s$ , $R_1(T)=r )$ are obtained by solving this system when $\theta=0$ .

Alternatively, we can obtain the final epidemic outcome from Corollary 3.1, through the transform

\begin{equation*}u_s(\theta,z_1,z_2) = \mathbb{E}\bigl[ {\textrm{e}}^{-\theta A(T)} \,z_1^{R_1(T)}\,z_2^{R_2(T)} \,\unicode{x1d7d9}_{S(T)=s} \bigr].\end{equation*}

For that, it suffices to rewrite (3.10) as

(3.12) \begin{equation}\sum_{s=k}^{n} \Biggl(\prod_{l=k+1}^{s}\dfrac{\beta_l}{\beta_l-\beta_k} \Biggr) q_{\alpha}(k,\theta,z_1,z_2)^{s} \, u_s(\theta,z_1,z_2) = q_{\alpha}(k,\theta,z_1,z_2)^{n+m} \prod_{l=k+1}^{n}\dfrac{\beta_l}{\beta_l-\beta_k}\end{equation}

for all $k=0,1,\ldots,n$ . For each chosen value of $\theta$ , $z_1$ , and $z_2$ , (3.12) constitutes a triangular system of $n+1$ linear equations that determine $u_s(\theta,z_1,z_2)$ . The probabilities $\mathbb{P}( S(T)=s )$ are obtained by solving this system when $\theta=0$ and $z_1=z_2=1$ . The Laplace transform of $A(T)$ and the probability generating function of $(R_1(T),R_2(T))$ are given by

\begin{align*}\mathbb{E}\bigl[ {\textrm{e}}^{-\theta A(T)} \bigr] &= \sum_{s=0}^n u_s(\theta,1,1), \\[5pt] \mathbb{E}\bigl[ z_1^{R_1(T)}\,z_2^{R_2(T)} \bigr] &= \sum_{s=0}^n u_s(0,z_1,z_2).\end{align*}

Finally, differentiating with respect to $\theta$ , $z_1$ , or $z_2$ in (3.10), we can also obtain various relations between the moments of $A(T)$ , $R_1(T)$ , $R_2(T)$ , and those of $S(T)$ . In particular, the expectations of $A(T)$ and $R_1(T)$ have the simple expressions

\begin{align*}\mathbb{E}[ A(T) ] &= (n+m-\mathbb{E}[ S(T) ]) \dfrac{{\textrm{d}}}{{\textrm{d}}\theta} q_{\alpha}(0,\theta,1,1)\big|_{\theta = 0}, \\[5pt] \mathbb{E}[ R_1(T) ] &= (n+m-\mathbb{E}[ S(T) ]) \dfrac{{\textrm{d}}}{{\textrm{d}} z} q_{\alpha}(0,0,z,1)\big|_{z = 1},\end{align*}

which are intuitive since the derivative $({\textrm{d}}/{\textrm{d}}\theta)q_{\alpha}(0,\theta,1,1)|_{\theta = 0}$ gives the mean severity due to an infected individual whose contamination process starts with initial distribution $\alpha(\!\cdot\!)$ , and $({\textrm{d}}/{\textrm{d}} z)q_{\alpha}(0,0,z,1)|_{z = 1}$ is the probability that he becomes a removed of type 1 at the end of his infectious period.

Note that in order to solve the systems (3.11) and (3.12), the coefficients $q_{\alpha}(k,\theta,z_1,z_2)$ must first be computed. This can be done in several ways, the most convenient option depending on the nature of the process X: these coefficients can be computed from the representation (3.9) either by direct calculation or by simulating the process X. Alternatively, they can be obtained by solving equation (3.6).

Remark 3.1. In the case where the infection rates are linear with the number of susceptibles ( $\beta_s = \beta s$ ), another way to determine the final epidemic outcome is to adapt the approach developed in Ball [4] to the case of Feller infectivity processes.

To close this section, we emphasize the flexibility of the model by presenting a few interesting classes of Feller processes that can be used to represent the infectious periods. We comment on the possible interpretations and applications of each class, and we show how to derive the corresponding expressions for the function $q(k,\theta,z_1,z_2,x)$ , whose knowledge allows us to solve the systems (3.11) and (3.12).

Example 3.1. (General epidemic with arbitrary infectious period distribution.) Let us assume as in Clancy [10] that each infected individual remains contagious for a duration distributed as a random variable $Z>0$ with cumulative distribution function $F(\!\cdot\!)$ . During this period, he makes contaminations at rate $\beta_{S(t)}$ , where $\beta_s$ is an arbitrary positive function of s. We consider here an extension of this model where there are two removal classes and where the type of removal of each individual depends on the duration of his infectious period. This can add interesting information in the model, for instance in situations where infected people are more likely to die when their illness lasts longer.

To describe this model in our setting, we take the infectivity process X as a deterministic process on $\mathcal{E}=\mathbb{R}^+$ such that $X(u)=u$ . Then $X=x$ simply means that the infectious period started x units of time ago. The contamination and severity functions are constant: $\gamma(x)=a(x)=1$ for all $x \in \mathbb{R}^+$ . When $X=x$ , the instantaneous removal rate is the failure rate associated to $F(\!\cdot\!)$ , i.e. ${\textrm{d}} F(x)/(1-F(x))$ . Of course, each infected individual starts from state 0, so $\alpha$ is the probability measure on $\mathbb{R}^+$ with all the mass concentrated on $x=0$ .

Let us assume that if the infectious period of an individual lasts y units of time, the probability that he becomes a removed of type 1 (resp. type 2) is p(y) (resp. $1-p(y)$ ). The removal rate functions $h_1$ and $h_2$ are then given by

\begin{equation*}h_1(x) = p(x) \dfrac{{\textrm{d}} F(x)}{1-F(x)}, \quad h_2(x) = (1-p(x)) \dfrac{{\textrm{d}} F(x)}{1-F(x)}.\end{equation*}

From (3.9), we then obtain

\begin{align*}q(k,\theta,z_1,z_2,x) &= \mathbb{E}_{x}\bigl[ z_1^{\mathbb{I}_{1}} \, z_2^{\mathbb{I}_{2}} \,{\textrm{e}}^{-(\theta + \beta_k)\tau_{\star}} \bigr] \\[5pt] &= \mathbb{E}\bigl[ z_1^{\mathbb{I}_{1}} \, z_2^{\mathbb{I}_{2}} \,{\textrm{e}}^{-(\theta + \beta_k)(Z-x)} \mid Z>x \bigr] \\[5pt] &= \mathbb{E}\bigl[ \mathbb{E}\bigl[ z_1^{\mathbb{I}_{1}} \, z_2^{\mathbb{I}_{2}}\mid Z \bigr] \,{\textrm{e}}^{-(\theta + \beta_k)(Z-x)} \mid Z>x \bigr],\end{align*}

and therefore

\begin{align*}q(k,\theta,z_1,z_2,x) &= \mathbb{E}\bigl[ (z_1 p(Z) + z_2(1-p(Z)))\, {\textrm{e}}^{-(\theta + \beta_k)(Z-x)} \mid Z>x \bigr]\\[5pt] &= \dfrac{1}{1-F(x)} \int_x^{\infty} (z_1 p(y) + z_2 (1-p(y))) \, {\textrm{e}}^{-(\theta + \beta_k)(y-x)} \, {\textrm{d}} F(y).\end{align*}

Alternatively, this expression for $q(k,\theta,z_1,z_2,x)$ can be obtained by using the fact that the generator of X is the differential operator $\mathcal{L} f = f'$ . So (3.6) gives the differential equation

\begin{equation*}f'(x) = \biggl(\theta + \beta_k + \dfrac{{\textrm{d}} F(x)}{1-F(x)} \biggr)\,f(x) - (z_1p(x)+z_2 (1-p(x)))\dfrac{{\textrm{d}} F(x)}{1-F(x)},\end{equation*}

which can be easily solved under the condition that f is bounded. Note that we only need to determine the coefficient $q_{\alpha}(k,\theta,z_1,z_2) \equiv q(k,\theta,z_1,z_2,0)$ in order to apply the transform (3.10).

Example 3.2. (SIR models driven by Markov chains.) Now consider the case where the infectivity process X is a Markov chain defined on a finite or discrete state space $\mathcal{E}$ . This process is characterized by the $|\mathcal{E}|\times |\mathcal{E}|$ intensity matrix Q (such that $Q_{jl}$ is the instantaneous transition rate from state j to state l when $l\neq j$ and $\sum_{l \in \mathcal{E}}Q_{jl} = 0 $ for all j) and by the $1\times |\mathcal{E}|$ initial vector $\boldsymbol{\alpha}$ with components $\alpha_j = \mathbb{P}( X(0)=j )$ . Assume that an individual whose infectivity process is in state $j \in \mathcal{E}$ has a cost a(j) per unit of time and makes contaminations at rate $\beta_{s} \gamma(j)$ if there are s available susceptibles. Moreover, when in state j, an individual becomes removed of type 1 at rate $h_1(j)$ and of type 2 at rate $h_2(j)$ .

Here, the generator of X is the operator $\mathcal{L}$ given by

\begin{equation*}(\mathcal{L} f)(j) = \sum_{l \in \mathcal{E}} Q_{jl}\, f(l).\end{equation*}

From (3.6), the coefficients $q(k,\theta,z_1,z_2,j)$ satisfy the recursion

(3.13) \begin{equation} \sum_{l \in \mathcal{E}} Q_{jl} f(l) = (\theta a(j)+ \beta_k \gamma(j) + h_1(j) +h_2(j) ) f(j) - z_1h_1(j) - z_2h_2(j)\end{equation}

for all $j \in \mathcal{E}$ . Below, we take a closer look at two particular cases in this class of models.

(1) In the case where $\mathcal{E} = \{1,2,\ldots,L\}$ , the states of X can be seen as a finite number of possible degrees of the illness a contaminated individual can go through during his infectious period, each of them with specific contamination and removal rates. These different degrees can be used to distinguish, for instance, periods of improvement and deterioration in the patient’s condition, or the possible reactions following a treatment against the disease. They can also simply reflect the fact that the parameters are subject to a more general random environment specific to each individual.

Identifying the functions $g\;:\; \mathcal{E} \to \mathbb{R}$ to the column vectors $\boldsymbol{g} = (g(1),g(2),\ldots,g(L))'$ and writing $\Delta_g = \text{diag}(\boldsymbol{g})$ , equation (3.13) can be rewritten in matrix form as

\begin{equation*}Q \boldsymbol{f} = (\theta \Delta_a+ \beta_k \Delta_{\gamma} + \Delta_{h_1} + \Delta_{h_2} ) \boldsymbol{f} - z_1 \boldsymbol{h}_1 - z_2 \boldsymbol{h}_2.\end{equation*}

The unique solution of this equation is

\begin{equation*}\boldsymbol{q}(k,\theta,z_1,z_2) = (\theta \Delta_a+ \beta_k \Delta_{\gamma} + \Delta_{h_1} + \Delta_{h_2} - Q )^{-1} (z_1 \boldsymbol{h}_1 + z_2 \boldsymbol{h}_2),\end{equation*}

and the coefficient $q_{\alpha}(k,\theta,z_1,z_2) = \boldsymbol{\alpha}\boldsymbol{q}(k,\theta,z_1,z_2)$ is easily computed.

(2) Now consider the case where X is an absorbing birth-and-death process on $\mathcal{E} = \{0,1,2,\ldots\}$ with birth rates $\lambda_j$ and death rates $\mu_j$ , where $\lambda_j, \mu_j>0$ for $j>0$ and $\lambda_0=\mu_0=0$ . Here, the generator Q is the infinite matrix of which the only non-zero components are

\begin{equation*}Q_{j,j-1} = \mu_j, \quad Q_{j,j+1} = \lambda_j, \quad Q_{j,j} = -\lambda_j-\mu_j, \quad {j\ge 1.}\end{equation*}

This kind of model can be used, for instance, to represent situations where the infection is caused by a population of diseased cells which reproduce inside the organism of the infected individual, the birth-and-death process keeping track of the evolution of their number.

Assume that when his infectivity process is in state $j>0$ , an individual transmits the disease at rate $\beta_s \gamma_j$ (when there are s available susceptibles) and has a cost $a_j$ per time unit. He can die from the illness and become a removed case of type 1 at rate $h_1(j)$ , but he can also recover and become a removed case of type 2 if his process reaches the absorbing state 0 before succumbing to the illness. Then the function $q(k,\theta,z_1,z_2,j)$ is the solution of

(3.14) \begin{equation}\mu_j\, f(j-1) - ( \lambda_j + \mu_j + \theta a(j)+ \beta_k \gamma(j) + h_1(j) ) f(j) + \lambda_j f(j+1) = -z_1 h_1(j) \end{equation}

for $j\ge 1$ , subject to the conditions that f is bounded and $f(0)=z_2$ . In the particular case of a linear birth-and-death process ( $\lambda_j = \lambda j$ , $\mu_j = \mu j$ , $a_j = aj$ , $\gamma_j = j$ , $h_1(j) = \kappa j$ ), the recursion (3.14) is easily solved and yields

\begin{equation*}q(k,\theta,z_1,z_2,j) = C + (z_2-C) \biggl(\dfrac{\lambda + \mu + \beta_k + a\theta + \kappa - \sqrt{(\lambda + \mu + \beta_k + a\theta + \kappa)^2 - 4 \lambda \mu}}{2\lambda}\biggr)^j,\end{equation*}

where $C = z_1 \kappa /(\beta_k + a\theta + \kappa)$ .

Example 3.3. (SIR models driven by diffusions.) In this last example, X is a diffusion process taking values in $\mathcal{E} \subseteq \mathbb{R}$ , i.e. X satisfies a stochastic differential equation of the form

\begin{equation*}{\textrm{d}} X(u) = \lambda(X(u)) \, {\textrm{d}} u + \sigma(X(u))\, {\textrm{d}} W(u),\end{equation*}

where W is a standard Brownian motion and the functions $\lambda\;:\; \mathcal{E} \to \mathbb{R}$ and $\sigma\;:\; \mathcal{E} \to \mathbb{R}^+$ satisfy the usual Lipschitz continuity condition. More generally, X could take its values in a subset of $\mathbb{R}^n$ . The corresponding SIR models can be regarded as Brownian perturbations of more classical models (like the one in the first example), where the Brownian component represents a noise resulting from errors in the parameter estimation, for example. As in the last example, the diffusion can also represent the density of infected cells inside the organism, which are responsible for the (transmission of the) disease.

Here, the generator of X is of the form

\begin{equation*}\mathcal{L} f(x) = a(x) f'(x) + \dfrac{1}{2}b(x)^2 f{''}(x),\end{equation*}

and a necessary condition for f to be part of $\mathcal{D}_\mathcal{L}$ is that f is bounded and twice differentiable. One of the simplest examples of diffusion is Brownian motion with zero drift, i.e. the diffusion with $\lambda(x)=0$ and $\sigma(x) = \sigma >0$ . For illustration, let us consider this example when the Brownian motion is reflected at level zero (so that $\mathcal{E} = \mathbb{R}^+$ ) and the infection rate $\gamma(x)=x$ is linear with the level of X. Let us also fix $a(x)=1$ so that the severity is the duration of the infectious periods. An infected individual recovers and becomes a removed of type 1 after a period of time $\textrm{Exp}(\mu_1)$ and thus $h_1(x) = \mu_1$ . He can also die from the illness at a rate $h_2(x)=\mu_2x$ , which increases with the level of the infectivity process. Then he becomes a removed of type 2. From (3.6), the function $q(k,\theta,z_1,z_2,x)$ satisfies the equation

(3.15) \begin{equation}\dfrac{1}{2}\sigma^2 f''(x) = [\theta + \mu_1 + (\beta_k + \mu_2)x ]\, f(x) - z_1\mu_1 - z_2\mu_2x,\end{equation}

with the constraint $f'(0)=0$ to account for the reflection of X at level zero. Defining the function g such that

\begin{equation*} f(x) = (\beta_k + \mu_2)^{-1} \pi(\omega z_1 \mu_1 - \eta z_2 \mu_2) g(\omega x + \eta) + (\beta_k + \mu_2)^{-1} z_2 \mu_2 ,\end{equation*}

with $\omega = 2^{{{1}/{3}}} (\beta_k+\mu_2) \sigma^{-{{2}/{3}}} (\beta_k + \mu_2)^{-{{2}/{3}}}$ and $\eta = 2^{{{1}/{3}}} (\theta+\mu_1) \sigma^{-{{2}/{3}}} (\beta_k + \mu_2)^{-{{2}/{3}}}$ , (3.15) turns into the standard inhomogeneous Airy equation

(3.16) \begin{equation} g''(y) - y g(y) + \dfrac{1}{\pi} = 0, \end{equation}

the solutions of which are well known (see e.g. Gil et al. [14]). Solving (3.16) with the constraints that f is bounded and that $f'(0)=0$ , we obtain

\begin{equation*}q(k,\theta,z_1,z_2,x) = \dfrac{z_2\mu_2}{\beta_k + \mu_2} + \dfrac{\omega z_1 \mu_1 - \eta z_2 \mu_2}{(\beta_k + \mu_2) \textit{Ai}'(\omega)} (\textit{Ai}'(\omega)\textit{Gi}(\eta+\omega x)-\textit{Gi}'(\omega)\textit{Ai}(\eta + \omega x) ),\end{equation*}

where the Airy function $\textit{Ai}(x)$ and the Scorer function $\textit{Gi}(x)$ are given by

\begin{equation*}\textit{Ai}(x) = \dfrac{1}{\pi} \int_0^{\infty} \cos\biggl(\dfrac{y^3}{3} + xy\biggr) \,{\textrm{d}} y, \quad\textit{Gi}(x) = \dfrac{1}{\pi} \int_0^{\infty} \sin\biggl(\dfrac{y^3}{3} + xy\biggr) \,{\textrm{d}} y.\end{equation*}

4. Contaminations per infected and final epidemic outcome computation

In Section 4.1 we analyse the severity, removal type, and number of contaminations made by a unique infected individual facing a given number of susceptibles. In Section 4.2 we show that this analysis provides an alternative procedure to compute the final epidemic outcome distributions obtained in Section 3.

4.1. Outcome per infected individual

To assess the risk of an epidemic spreading in a population, an interesting measure is the effect of a single infected individual on a fixed group of s susceptibles. This measure is simpler to obtain than the final epidemic outcome distribution and provides many useful estimators of the threat caused by the disease (such as the mean number of contacts per individual, the cost per individual, the probability of death).

To measure this effect, we are going to consider a simpler epidemic model in which only one infected individual is present and faces a fixed number of susceptibles who are unable to spread the disease (even if they are contaminated). This model can be seen as a generalization of the carrier-borne epidemic of Weiss [25]. More precisely, we consider a (unique) generic infected individual whose infectivity process is $\{X(t)\mid t \ge 0\}$ with X(0) given by the probability measure $\alpha(\!\cdot\!)$ . Here, $t=0$ corresponds to the time when the infection period starts and $t=\tau_{\star}$ is the time when it ends. The individual is facing an arbitrary but fixed number s of susceptibles, and we assume that those of them who are contaminated are directly removed (i.e. they play no role in the contamination of the remaining susceptibles). We are going to determine the joint distribution of several statistics. The first one is the total number $\mathcal{N}_s$ of contaminations made by the individual among the s susceptibles. The second one is the total severity $D(\tau_{\star})$ due to the individual, with

\begin{equation*} D(t) = \int_{0}^{t} a(X(u)) \, {\textrm{d}} u.\end{equation*}

Finally, $\mathbb{I}_{1}$ (resp. $\mathbb{I}_{2}$ ) are the indicators that he becomes a removed of type 1 (resp. type 2) at the end of his infectious period. We start with the following preliminary result, which will be used both here and in Section 4.2.

Lemma 4.1. For each $k \in \{0,1,2,\ldots,n\}$ , $\theta \ge 0$ , and $z_1,z_2 \in\!]0,1]$ ,

(4.1) \begin{equation}\mathbb{E}\Biggl[ \Biggl(\prod_{l=k+1}^{s-\mathcal{N}_s}\dfrac{\beta_l}{\beta_l-\beta_k} \Biggr) \,{\textrm{e}}^{-\theta D(\tau_{\star})} \, z_1^{\mathbb{I}_{1}} \, z_2^{\mathbb{I}_{2}} \Biggr] = \Biggl(\prod_{l=k+1}^{s}\dfrac{\beta_l}{\beta_l-\beta_k} \Biggr) q_{\alpha}(k,\theta,z_1,z_2), \end{equation}

where $q_{\alpha}(k,\theta,z_1,z_2)$ is given in (3.7).

Proof. Let $\mathcal{N}_s(t)$ be the number of contaminations made by the infected individual after t units of time and let $\mathbb{I}_{j}(t)$ be the indicator of the event $[t\ge \tau_{\star}$ , $ \mathbb{I}_{j}=1]$ (for $j=1,2$ ). To show (4.1), the reasoning is the same as that used to obtain (3.10). We look for a function $b(\!\cdot\!)$ and a bounded and measurable $f\;:\; \mathcal{E} \cup \{\star\} \to ]0,1]$ such that the process

\begin{equation*} M(t) = b(\mathcal{N}_s(t)) \, {\textrm{e}}^{-\theta D(t)} \, z_1^{\mathbb{I}_{1}(t)} \, z_2^{\mathbb{I}_{2}(t)} \, f(X(t))\end{equation*}

is a martingale with respect to the filtration $\mathcal{G} = (\mathcal{G}_t)_{t \ge 0}$ containing the history of X(t), $D(t)$ , and $\mathcal{N}_s(t)$ up to time t. Here, $\{\star\}$ denotes an artificial, absorbing state in which the process X is sent as soon as the infected individual is removed. We proceed as in the proof of Proposition 3.1: we fix $t_0>0$ and look for some $b(\!\cdot\!)$ and $f(\!\cdot\!)$ such that the right derivative $m_r^{\prime}(t)$ of $m(t) = \mathbb{E}[ M(t) \mid \mathcal{G}_{t_0} ]$ equals zero for all $t \ge t_0$ . The tower property allows us to write

(4.2) \begin{equation}m(t+{\textrm{d}} t) = \mathbb{E}[ M(t+{\textrm{d}} t) \mid \mathcal{G}_{t_0} ]= \mathbb{E}\bigl[ \mathbb{E}[ M(t+{\textrm{d}} t) \mid \mathcal{G}_{t} ] \mid \mathcal{G}_{t_0} \bigr]. \end{equation}

By examining the possible events on the time interval $[t, t+{\textrm{d}} t[$ , we find that the following equality holds up to ${\textrm{o}}({\textrm{d}} t)$ :

\begin{align*} & \mathbb{E}[ M(t+{\textrm{d}} t) \mid \mathcal{G}_{t} ] \\[5pt] &\quad = b(\mathcal{N}_s(t) + 1) \, \mathbb{E}\bigl[{\textrm{e}}^{-\theta D(t+{\textrm{d}} t)} f(X(t+{\textrm{d}} t)) \mid \mathcal{G}_{t} \bigr] (\beta_{s-\mathcal{N}_s(t)}\gamma(X(t)) \,{\textrm{d}} t + {\textrm{o}}({\textrm{d}} t) ) \\[5pt] &\quad\quad + b(\mathcal{N}_s(t)) \,\mathbb{E}\bigl[ {\textrm{e}}^{-\theta D(t+{\textrm{d}} t)} \mid \mathcal{G}_{t} \bigr]\, f(\!\star\!)(z_1 h_1(X(t)) \,{\textrm{d}} t + z_2 h_2(X(t)) \,{\textrm{d}} t + {\textrm{o}}({\textrm{d}} t) ) \\[5pt] & \quad \quad + b(\mathcal{N}_s(t)) \, \mathbb{E}\bigl[ {\textrm{e}}^{-\theta D(t+{\textrm{d}} t)}f(X(t+{\textrm{d}} t)) \mid \mathcal{G}_{t} \bigr] \\[5pt] &\quad\quad \times (1- (\beta_{s-\mathcal{N}_s(t)}\gamma(X(t)) + h_1(X(t)) + h_2(X(t))) \,{\textrm{d}} t + {\textrm{o}}({\textrm{d}} t) ).\end{align*}

Subtracting m(t) from both sides of (4.2), dividing each side by ${\textrm{d}} t$ , and taking ${\textrm{d}} t \to 0$ , we obtain as in the proof of Proposition 3.1 that a sufficient condition for $m_r^{\prime}(t)$ to be equal to zero is that $b(\!\cdot\!)$ and $f(\!\cdot\!)$ satisfy

(4.3) \begin{align}& b(l) [\mathcal{L} f(x) - (\beta_{s-l}\gamma(x) + \theta a(x) + h_1(x) + h_2(x) )\,f(x) + (z_1h_1(x) + z_2h_2(x)) f(\!\star\!) ] \nonumber \\[5pt] & \quad = - \beta_{s-l}\, b(l+1) \, \gamma(x) \,f(x) \quad \text{{for all $x \in \mathcal{E}$, $ l \in \{0,1,\ldots,s\}$.}}\end{align}

Next, we proceed as in the proof of Theorem 3.1: for each fixed $k \in \{0,1,2,\ldots,n\}$ , a solution of (4.3) is obtained by choosing $f \equiv f_k$ such that $f_k(\!\star\!)=1$ and $f_k(x) = q(k,\theta,z_1,z_2,x)$ for $x \in \mathcal{E}$ . In that way, $f_k$ satisfies (3.6), and (4.3) reduces to the equation

\begin{align*} b(l) & = 0 \quad \text{for }l>s-k, \\[5pt] \beta_{s-l-k} b(l) &= \beta_s b(l+1) \quad \text{for }l\le s-k,\end{align*}

which is easily solved and yields $b(l)=\prod_{l=k+1}^{s-l}(\beta_l/(\beta_l-\beta_k))$ . As a consequence, the process

\begin{equation*}\Biggl\{ \Biggl(\prod_{l=k+1}^{s-\mathcal{N}_s(t)}\dfrac{\beta_l}{\beta_l-\beta_k} \Biggr) \,{\textrm{e}}^{-\theta D(t)} \, z_1^{\mathbb{I}_{1}(t)} \, z_2^{\mathbb{I}_{2}(t)} \, q(k,\theta,z_1,z_2,X(t))^{\unicode{x1d7d9}_{t<\tau_{\star}}} \Biggr\}\end{equation*}

is a martingale with respect to $\mathcal{G}_t$ , and (4.1) follows by applying the optional stopping theorem with $\tau_{\star}$ as the stopping time.

Henceforth, we use the product notation

(4.4) \begin{equation}\prod_{\substack{j=l \\ j\neq v}}^{L} c_j \equiv \Biggl( \prod_{j=l}^{v-1} c_j \Biggr) \Biggl( \prod_{j=v+1}^{L} c_j \Biggr), \end{equation}

with the usual convention that $ \prod_{j=r}^{r-k} c_j$ is equal to one if $k=1$ and equal to zero if $k>1$ . Note that with this convention, (4.4) is equal to zero if $v \notin \{l, l+1,\ldots,L\}$ .

From Lemma 4.1, we can obtain an explicit expression (in terms of the coefficients $q_{\alpha}(k,\theta,z_1,z_2)$ ) for the joint distribution of $\mathcal{N}_s$ , $D(\tau_{\star})$ , $\mathbb{I}_{1}$ , and $\mathbb{I}_{2}$ .

Proposition 4.1. For all $\theta \ge 0$ and $x,z_1,z_2 \in\!]0,1]$ ,

(4.5) \begin{align}\mathbb{E}\bigl[ x^{s-\mathcal{N}_s}\,{\textrm{e}}^{-\theta D(\tau_{\star})} \, z_1^{\mathbb{I}_{1}} \, z_2^{\mathbb{I}_{2}} \bigr] = \sum_{l=0}^{s} \Biggl(\sum_{v=l}^{s} \, Q_s(l,v) \, q_{\alpha}(v,\theta,z_1,z_2)\Biggr) x^l, \end{align}

where

\begin{equation*} Q_s(l,v) = \Biggl( \prod_{j=l+1}^{s} \beta_j \Biggr) \!\left(\rule{0pt}{24pt}\right.\! \prod_{\substack{j=l \\ j\neq v}}^{s}\dfrac{1}{\beta_j-\beta_v} \!\left.\rule{0pt}{24pt}\right)\!.\end{equation*}

Proof. Let us first show that for $0 \le l$ , $L \le s$ ,

(4.6) \begin{equation}\Biggl( \prod_{j=l+1}^{L} \beta_j \Biggr) \sum_{v=l}^s \!\left(\rule{0pt}{24pt}\right.\! \prod_{\substack{j=l \\ j\neq v}}^{L}\dfrac{1}{\beta_j-\beta_v} \!\left.\rule{0pt}{24pt}\right)\! = \delta_{l,L}, \end{equation}

where $\delta$ is the Kronecker delta. Note that the left-hand side of (4.6) can be rewritten as

(4.7) \begin{equation}H(l,L) = \Biggl( \prod_{j=l+1}^{L} \beta_j \Biggr) \sum_{v=l}^L \!\left(\rule{0pt}{24pt}\right.\! \prod_{\substack{j=l \\ j\neq v}}^{L}\dfrac{1}{\beta_j-\beta_v}\!\left.\rule{0pt}{24pt}\right)\!. \end{equation}

We clearly have $H(L,L) = 1$ . For $l>L$ , the first product appearing in (4.7) equals zero so that $H(l,L) = 0$ . For $0<l<L$ , we show that $H(l,L) = 0$ by writing the polynomial $P(x) = x$ in the Lagrange basis $\{p_0(x), p_1(x),\ldots, p_{L-l}(x)\}$ , where

\begin{equation*}p_v(x) = \prod_{\substack{j=0 \\ j\neq v}}^{L-l}\dfrac{\beta_{l+j} - x}{\beta_{l+j}-\beta_{l+v}}.\end{equation*}

This yields

\begin{equation*}x = \sum_{v=0}^{L-l} \beta_{v+l} \, p_v(x),\end{equation*}

and this equality with $x=0$ is equivalent to $H(l,L)=0$ . Finally, we show that $H(0,L)=0$ if $L>0$ in a similar way, writing

\begin{equation*}H(0,L) = 1 - \sum_{v=1}^L \!\left(\rule{0pt}{24pt}\right.\! \prod_{\substack{j=1 \\ j\neq v}}^{L}\dfrac{\beta_j}{\beta_j-\beta_v}\!\left.\rule{0pt}{24pt}\right)\! ,\end{equation*}

and using the fact that the sum in the last equality is the constant polynomial $P(x)=1$ written in the Lagrange basis $\{p_1(x), p_1(x),\ldots, p_{L}(x)\}$ .

Now, from (4.6), where L is replaced by the random variable $s-\mathcal{N}_s$ , we obtain

\begin{align*}&\mathbb{E}\left[ x^{s-\mathcal{N}_s}\,{\textrm{e}}^{-\theta D(\tau_{\star})} \, z_1^{\mathbb{I}_{1}} \, z_2^{\mathbb{I}_{2}} \right] \\[5pt] & \quad = \mathbb{E}\!\left[\rule{0pt}{24pt}\right.\! \sum_{l=0}^{s} \Biggl(\prod_{j=l+1}^{s-\mathcal{N}_s} \beta_j \Biggr)\sum_{v=l}^{s} \!\left(\rule{0pt}{24pt}\right.\!\prod_{\substack{j=l \\ j\neq v}}^{s-\mathcal{N}_s}\dfrac{1}{\beta_j-\beta_v}\!\left.\rule{0pt}{24pt}\right)\! x^l \,{\textrm{e}}^{-\theta D(\tau_{\star})} \, z_1^{\mathbb{I}_{1}} \, z_2^{\mathbb{I}_{2}}\!\left.\rule{0pt}{24pt}\right]\! \\[5pt] & \quad =\sum_{v=0}^{s} \sum_{l=0}^{v} \mathbb{E}\Biggl[\Biggl(\prod_{j=v+1}^{s-\mathcal{N}_s}\dfrac{\beta_j}{\beta_j-\beta_v} \Biggr) \,{\textrm{e}}^{-\theta D(\tau_{\star})} \, z_1^{\mathbb{I}_{1}} \, z_2^{\mathbb{I}_{2}}\Biggr] \Biggl(\prod_{j=l+1}^{v} \beta_j \Biggr) \Biggl(\prod_{j=l}^{v-1}\dfrac{1}{\beta_j-\beta_v}\Biggr) x^l.\end{align*}

Using Lemma 4.1, we then have

\begin{align*}&\mathbb{E}\bigl[ x^{s-\mathcal{N}_s}\,{\textrm{e}}^{-\theta D(\tau_{\star})} \, z_1^{\mathbb{I}_{1}} \, z_2^{\mathbb{I}_{2}} \bigr] \\[5pt] & \quad = \sum_{v=0}^{s} \sum_{l=0}^{v} \Biggl(\prod_{j=v+1}^{s}\dfrac{\beta_j}{\beta_j-\beta_v} \Biggr) q_{\alpha}(v,\theta,z_1,z_2) \Biggl(\prod_{j=l+1}^{v} \beta_j \Biggr) \Biggl(\prod_{j=l}^{v-1}\dfrac{1}{\beta_j-\beta_v}\Biggr) x^l \\[5pt] &\quad = \sum_{l=0}^{s} \sum_{v=l}^{s} \Biggl( \prod_{j=l+1}^{s} \beta_j \Biggr) \!\left(\rule{0pt}{24pt}\right.\!\prod_{\substack{j=l \\ j\neq v}}^{s}\dfrac{1}{\beta_j-\beta_v}\!\left.\rule{0pt}{24pt}\right)\! q_{\alpha}(v,\theta,z_1,z_2) \, x^l,\end{align*}

hence the stated expression.

From Proposition 4.1, we can easily determine the distribution of various quantities of interest related to the behaviour of a typical infected individual. By differentiating both sides of (4.5) $s-r$ times with respect to x and taking $x=0$ , we obtain that for all $\theta \ge 0$ , $z_1,z_2 \in\!]0,1]$ , and $r \in \{0,1,\ldots,s\}$ ,

(4.8) \begin{align}\mathbb{E}\bigl[ {\textrm{e}}^{-\theta D(\tau_{\star})} \, z_1^{\mathbb{I}_{1}} \, z_2^{\mathbb{I}_{2}} \, \unicode{x1d7d9}_{\mathcal{N}_s=r} \bigr] = \sum_{v=s-r}^{s} \, Q_s(s-r,v) \, q_{\alpha}(v,\theta,z_1,z_2). \end{align}

In particular, the distribution of $\mathcal{N}_s$ together with the type of elimination is given by

(4.9) \begin{equation}\mathbb{P}( \mathcal{N}_s=s-r, \mathbb{I}_{1} = 0 ) = \sum_{v=r}^{s} \, Q_s(r,v) \, q_{\alpha}(v,0,0,1), \end{equation}

and $\mathbb{P}( \mathcal{N}_s=s-r )$ is given by the same expression except that $q_{\alpha}(k,0,0,1)$ is replaced by $q_{\alpha}(k,0,1,1)$ . The factorial moments of $s-\mathcal{N}_s$ are obtained by differentiating both sides of (4.5) k times and then taking $\theta=0$ and $x =z_1=z_2=1$ :

\begin{align*} \mathbb{E}[ (s-\mathcal{N}_s)_{[k]} ] = \sum_{l=k}^{s} \dfrac{l!}{(l-k)!} \Biggl(\sum_{v=l}^{s} \, Q_s(l,v) \, q_{\alpha}(v,0,1,1)\Biggr).\end{align*}

Special case. When $\beta_l=l$ for all l, the above expressions become much simpler. We have

\begin{equation*}Q_s(l,v) = (\!-\!1)^{v-l} \, \dfrac{s!}{l! \,(v-l)! \, (s-v)!},\end{equation*}

and using the binomial theorem, (4.5) can be rewritten as

\begin{align*}\mathbb{E}\bigl[ x^{s-\mathcal{N}_s}\,{\textrm{e}}^{-\theta D(\tau_{\star})} \, z_1^{\mathbb{I}_{1}} \, z_2^{\mathbb{I}_{2}} \bigr]&= \sum_{v=0}^{s} \binom{s}{v} q_{\alpha}(v,\theta,z_1,z_2) \sum_{l=0}^{v} \binom{v}{l} (\!-\!1)^{v-l} x^l \\[5pt] &= \sum_{v=0}^{s} \binom{s}{v} \, q_{\alpha}(v,\theta,z_1,z_2)\, (x-1)^v.\end{align*}

In particular, the probability that an infected individual facing s susceptibles makes k contaminations and then becomes a removed of type 1 is

(4.10) \begin{equation}\mathbb{P}( \mathcal{N}_s=k, \mathbb{I}_{1} = 0 ) = \binom{s}{k} \sum_{l=0}^{k} \binom{k}{l} (\!-\!1)^{k-l} \, q_{\alpha}(s-l,0,0,1), \end{equation}

and the factorial moments of $s-\mathcal{N}_s$ are given by the simple formula

\begin{equation*} \mathbb{E}\bigl[ (s-\mathcal{N}_s)_{[k]} \bigr] = \dfrac{s!}{(s-k)!} \, q_{\alpha}(k,0,1,1).\end{equation*}

Note that from its definition (3.7) (with $\beta_l = l$ ), the coefficients $q_{\alpha}(k,0,z_1,z_2)$ can be written here as $q_{\alpha}(k,0,z_1,z_2) = \mathbb{E}_{\alpha}[ z_1^{\mathbb{I}_{1}} \, z_2^{\mathbb{I}_{2}} \, W^k ]$ , with

\begin{equation*}V = \int_0^{\tau_{\star}} \gamma(X(u)) \,{\textrm{d}} u, \quad W = {\textrm{e}}^{-V}.\end{equation*}

Using this representation in (4.10), we have

\begin{align*} \mathbb{P}( \mathcal{N}_s=k ) &= \mathbb{E}_{\alpha}\Biggl[\binom{s}{k} (\!-\!1)^k \, W^{s-k} \sum_{l=0}^{k} \binom{k}{l} (\!-\!1)^{l} \, W^{k-l}\Biggr] = \mathbb{E}_{\alpha}\biggl[\binom{s}{k} W^{s-k} \, (1-W)^{k}\biggr].\end{align*}

So, in this special case, the variable $\mathcal{N}_s$ follows a mixed binomial distribution with parameters s and the (random) probability $1-W$ . Note that this result is intuitive: given that we know the cumulative infection rate generated by the infective over his entire infectious period (i.e. the random variable V), each of the s susceptibles has probability $1-{\textrm{e}}^{-V}$ of escaping infection. The result then follows from the fact that the susceptibles behave independently of each other.

In the general case where the rates $\beta_s$ are arbitrary, this intuition no longer holds. The probabilities $\mathbb{P}( \mathcal{N}_s=k )$ must then be computed using (4.9). Alternatively, they can be computed by simulating the random variable V, noting that from (4.9) we have the representation

\begin{equation*}\mathbb{P}( \mathcal{N}_s \le k ) = \mathbb{E}_{\alpha}[\mathbb{P}( T_{k} > V ) ]\end{equation*}

for $k<s$ , where $T_{k}$ is independent of V and follows a hypoexponential distribution with parameters $\beta_s, \beta_{s-1}, \ldots, \beta_{s-k}$ .

4.2. Artificial time and final outcome computation

In Section 3 we showed that the final epidemic outcome distribution can be computed by solving the linear system (3.11) or (3.12). However, these systems are often very unstable in practice, even for small values of n and m (see e.g. Demiris and O’Neill [12], House et al. [15], and the numerical illustration below). In this section we link our epidemic process to another, discrete-time process with the same outcome distribution. Such discrete-time equivalent models have already been used to obtain a system like (3.12) in the setting of the general epidemic (see e.g. Lefèvre and Utev [17] and Ball [3]), and to facilitate the study of various discrete-time epidemics (see e.g. Scalia-Tomba [23]). Here, we extend the approach to our setting and we take the opposite direction: we use the discrete-time model to provide an alternative, more stable method to compute the final epidemic outcome of the original model.

The idea is to make an artificial time change to follow the infections caused by every individual taken one after another. Informally stated, the dynamics of the original epidemic process is modified in the following way. At the beginning of the epidemic, all the initial infected individuals are sent to a queue outside the system. Then they are activated one after another. When an individual (labelled i, say) is activated, his infectious period starts and is governed by an independent copy $\{X_i(u)\}$ of $\{X(u)\}$ . When s susceptibles are present, he makes contaminations at rate $\beta_s \gamma(X_i(u))$ and he is removed in one of the two removal classes at rates $h_1(X_i(u))$ and $h_2(X_i(u))$ . When a susceptible is contaminated, he is sent to the queue outside the system. At the end of his infectious period, the infective is removed in one of the two types of elimination, and the next individual in the queue is activated.

In this modified system, there is only one infective in the population at a time. It leads to a discrete-time representation $\{(\mathcal{S}_{\tau},\mathcal{I}_{\tau},\mathcal{R}^{(1)}_{\tau},\mathcal{R}^{(2)}_{\tau},\mathcal{A}_{\tau}) \mid \tau= 0, 1, 2, \ldots\}$ , where $\mathcal{S}_{\tau}$ , $\mathcal{I}_{\tau}$ , $\mathcal{R}^{(1)}_{\tau}$ , and $\mathcal{R}^{(2)}_{\tau}$ , respectively, are the number of remaining susceptibles, infected individuals, and removed cases after the first $\tau$ removals. The variable $\mathcal{A}_{\tau}$ is the severity due to the $\tau$ first activated infectives. The initial state of the process is similar to that of the continuous-time model, i.e. $\mathcal{S}_0 = n$ , $\mathcal{I}_0 = m$ , and $\mathcal{A}_0=\mathcal{R}_0^{(1)}=\mathcal{R}_0^{(2)}=0$ . The epidemic terminates at the first time $\mathcal{T}$ when $\mathcal{I}_{\mathcal{T}}=0$ , i.e. when the active infective is removed and there are no more individuals waiting in the queue. It can be expressed as

(4.11) \begin{equation}\mathcal{T} = \inf \{\tau \ge 0 \mid \tau + \mathcal{S}_{\tau} = n+m \}. \end{equation}

Considering this discrete-time model instead of the original one, we lose all information about the transient behaviour of the epidemic. However, Lemma 4.1 can be used to show that the final epidemic outcome is the same in distribution for both models.

Proposition 4.2. For all $\theta \ge 0$ and $x, z_1,z_2 \in\!]0,1]$ ,

\begin{equation*} \mathbb{E}\bigl[ x^{S(T)} \, {\textrm{e}}^{-\theta A(T)} \, z_1^{R_1(T)} \, z_2^{R_2(T)} \bigr] = \mathbb{E}\Bigl[ x^{\mathcal{S}_{\mathcal{T}}} \, {\textrm{e}}^{-\theta \mathcal{A}_{\mathcal{T}}} \, z_1^{\mathcal{R}_{\mathcal{T}}^{(1)}} \, z_2^{\mathcal{R}_{\mathcal{T}}^{(2)}} \Bigr].\end{equation*}

Proof. We work with the discrete-time model introduced in this section. Let $D_{\tau}=_d D$ be the total severity due to the $\tau$ th activated infected individual, $\mathcal{N}_{\mathcal{S}_{\tau}}$ his number of contaminations and $\mathbb{I}_{1,\tau} =_d \mathbb{I}_{1}$ , $\mathbb{I}_{2,\tau} =_d \mathbb{I}_{2}$ the indicators of the events that he becomes a removed of type 1 or 2. The following equalities hold:

\begin{align*}& \mathcal{S}_{\tau+1} = \mathcal{S}_{\tau} - \mathcal{N}_{\mathcal{S}_{\tau}}, \\[5pt] & \mathcal{A}_{\tau+1} = \mathcal{A}_{\tau} + D_{\tau+1}, \\[5pt] & \mathcal{R}_{\tau+1}^{(r)} = \mathcal{R}_{\tau}^{(r)} + \mathbb{I}_{r,\tau+1}, \quad {r=1,2.}\end{align*}

Let $(\mathcal{H}_{\tau})_{\tau \ge 0}$ be the filtration $\sigma\{\mathcal{S}_{s}, \mathcal{A}_{s}, \mathcal{R}_{s}^{(1)},\mathcal{R}_{s}^{(2)} \mid s \le \tau\}$ and fix $k \in \{0,1,\ldots,n\}$ . Using the three relations above, we can write

\begin{align*}& \mathbb{E}\Biggl[ \Biggl(\prod_{l=k+1}^{\mathcal{S}_{\tau+1}}\dfrac{\beta_l}{\beta_l-\beta_k} \Biggr) \,{\textrm{e}}^{-\theta \mathcal{A}_{\tau+1}} \, {z_1}^{\mathcal{R}_{\tau+1}^{(1)}} \, {z_2}^{\mathcal{R}_{\tau+1}^{(2)}} \biggm| \mathcal{H}_{\tau} \Biggr] \\[5pt] & \quad = {\textrm{e}}^{-\theta \mathcal{A}_\tau} \, {z_1}^{\mathcal{R}_{\tau}^{(1)}} \, {z_2}^{\mathcal{R}_{\tau}^{(2)}} \, \mathbb{E}\Biggl[ \Biggl(\prod_{l=k+1}^{\mathcal{S}_{\tau}-\mathcal{N}_{\mathcal{S}_{\tau}}}\dfrac{\beta_l}{\beta_l-\beta_k} \Biggr) \,{\textrm{e}}^{-\theta D_{\tau+1}} \, {z_1}^{\mathbb{I}_{1,\tau+1}} \, {z_2}^{\mathbb{I}_{2,\tau+1}} \biggm| \mathcal{S}_{\tau} \Biggr] \\[5pt] & \quad = {\textrm{e}}^{-\theta \mathcal{A}_{\tau}} \, {z_1}^{\mathcal{R}_{\tau}^{(1)}} \, {z_2}^{\mathcal{R}_{\tau}^{(2)}} \, \Biggl(\prod_{l=k+1}^{\mathcal{S}_{\tau}}\dfrac{\beta_l}{\beta_l-\beta_k} \Biggr) q_{\alpha}(k,\theta,z_1,z_2),\end{align*}

where the last equality follows directly from (4.1). Dividing both sides of the last equality by $q_{\alpha}(k,\theta,z_1,z_2)^{\tau+1}$ , we see that the process

\begin{equation*}\Biggl\{\Biggl(\prod_{l=k+1}^{\mathcal{S}_{\tau}}\dfrac{\beta_l}{\beta_l-\beta_k} \Biggr) \,{\textrm{e}}^{-\theta \mathcal{A}_{\tau}} \, {z_1}^{\mathcal{R}_{\tau}^{(1)}} \, {z_2}^{\mathcal{R}_{\tau}^{(2)}}\, q_{\alpha}(k,\theta,z_1,z_2)^{-\tau}\Biggr\}\end{equation*}

is a martingale with respect to $(\mathcal{H}_{\tau})_{\tau \ge 0}$ . We can therefore apply the optional stopping theorem with $\mathcal{T}$ as the stopping time, and obtain

\begin{equation*}\mathbb{E}\Biggl[ \Biggl(\prod_{l=k+1}^{\mathcal{S}_{\mathcal{T}}}\dfrac{\beta_l}{\beta_l-\beta_k} \Biggr) \,{\textrm{e}}^{-\theta \mathcal{A}_\mathcal{T}} \, {z_1}^{\mathcal{R}_{\mathcal{T}}^{(1)}} \, {z_2}^{\mathcal{R}_{\mathcal{T}}^{(2)}}\, q_{\alpha}(k,\theta,z_1,z_2)^{-\mathcal{T}} \Biggr] = \prod_{l=k+1}^{n}\dfrac{\beta_l}{\beta_l-\beta_k}.\end{equation*}

From (4.11), we have the relation $\mathcal{T} = n+m-\mathcal{S}_{\mathcal{T}}$ and thus

(4.12) \begin{align} &\mathbb{E}\Biggl[ \Biggl(\prod_{l=k+1}^{\mathcal{S}_{\mathcal{T}}}\dfrac{\beta_l}{\beta_l-\beta_k} \Biggr) q_{\alpha}(k,\theta,z_1,z_2)^{\mathcal{S}_{\mathcal{T}}} \, {\textrm{e}}^{-\theta \mathcal{A}_\mathcal{T}} \, {z_1}^{\mathcal{R}_{\mathcal{T}}^{(1)}} \, {z_2}^{\mathcal{R}_{\mathcal{T}}^{(2)}} \Biggr] \notag \\[5pt] &\quad = q_{\alpha}(k,\theta,z_1,z_2)^{n+m}\prod_{l=k+1}^{n}\dfrac{\beta_l}{\beta_l-\beta_k}.\end{align}

The result follows by comparing (4.12) and (3.10), since their left-hand sides take the same value for all $k \in \{0,1,\ldots,n\}$ , $\theta \ge 0$ , and $z_1,z_2 \in\!]0,1]$ .

Using Proposition 4.2, we easily obtain a recursive procedure to compute the final outcome of the original continuous-time epidemic model of Section 2.

Proposition 4.3. For all $\theta \ge 0$ , $z_1,z_2 \in\!]0,1]$ , and $g\;:\; \{0,1,\ldots,n\}\to \mathbb{R}$ ,

\begin{equation*} \mathbb{E}\bigl[ g(S(T)) \, {\textrm{e}}^{-\theta A(T)} \, z_1^{R_1(T)} \, z_2^{R_2(T)} \mid S(0)=s, I(0)=i \bigr] = G_{s,i}(\theta,z_1,z_2),\end{equation*}

where the $G_{s,i}(\theta,z_1,z_2)$ satisfy the recursion

(4.13) \begin{equation} \begin{aligned}G_{s,i}(\theta,z_1,z_2) & = \sum_{r=0}^s L(s,r)\, G_{s-r,i+r-1}(\theta,z_1,z_2), \quad i >0,\\[5pt] G_{s,0}(\theta,z_1,z_2) & = g(s),\end{aligned}\end{equation}

where L(s, r) denotes the expectation (4.8).

Figure 1. Distribution of S(T) in Example 3.1, obtained from the system (3.12) (a) or from the recursion (4.13) (b). The parameters are $n=47$ , $m=3$ , $\beta_s = 2s/n$ , and $Z\sim \text{Beta}(2,3)$ .

Figure 2. Distribution of S(T) in Example 3.1, obtained from the system (3.12) (a) or from the recursion (4.13) (b). The parameters are $n=95$ , $m=5$ , $\beta_s = 2s/n$ , and $Z\sim \text{Beta}(2,3)$ .

Figure 3. Distribution of S(T) in Example 3.1, obtained from the recursion (4.13) when $n=490, m=10$ (a) or $n=980, m=20$ (b). The other parameters are $\beta_s = 2s/n$ and $Z\sim \text{Beta}(2,3)$ .

Proof. It is obvious that $G_{s,0}(\theta,z_1,z_2) = g(s)$ . When $i>0$ , we use Proposition 4.2 to write

\begin{align*}G_{s,i}(\theta,z_1,z_2) &= \mathbb{E}\bigl[ g(S(T)) \, {\textrm{e}}^{-\theta A(T)} \, z_1^{R_1(T)} \, z_2^{R_2(T)} \mid S(0)=s, I(0)=i \bigr] \\[5pt] & = \mathbb{E}\bigl[ g(\mathcal{S}_{\mathcal{T}}) \, {\textrm{e}}^{-\theta \mathcal{A}_{\mathcal{T}}} \, z_1^{\mathcal{R}^{(1)}_{\mathcal{T}}} \, z_2^{\mathcal{R}^{(2)}_{\mathcal{T}}}\mid \mathcal{S}(0)=s, \mathcal{I}(0)=i \bigr].\end{align*}

Working in the artificial discrete-time version, we rewrite the last equality by conditioning on the effect of the first infected individual who is activated. Letting $D$ denote his severity, $\mathcal{N}_s$ his number of contaminations, and $\mathbb{I}_{1}$ , $\mathbb{I}_{2}$ the indicators of his removal type, it suffices to use the tower property, the strong Markov property, and the independence between infected individuals to obtain

\begin{align*}G_{s,i}(\theta,z_1,z_2) &= \mathbb{E}\Bigl[ \mathbb{E}\Bigl[ g(\mathcal{S}_{\mathcal{T}}) \, {\textrm{e}}^{-\theta \mathcal{A}_{\mathcal{T}}} \, z_1^{\mathcal{R}^{(1)}_{\mathcal{T}}} \, z_2^{\mathcal{R}^{(2)}_{\mathcal{T}}}\mid D,\mathbb{I}_{1},\mathbb{I}_{2} \Bigr]\mid \mathcal{S}(0)=s,\mathcal{I}(0)=i \Bigr] \\[5pt] &= \sum_{r=0}^s \mathbb{E}\bigl[ {\textrm{e}}^{-\theta D} \, z_1^{\mathbb{I}_{1}} \, z_2^{\mathbb{I}_{2}} \, \unicode{x1d7d9}_{\mathcal{N}_s=r} \bigr] G_{s-r,i+r-1}(\theta,z_1,z_2),\end{align*}

which is (4.13).

We conclude with a small numerical illustration. In Figures 1, 2, and 3, we consider the distribution of S(T) for the ‘general epidemic model with arbitrary infectious periods’ described in the first example of Section 3, when the duration of the infectious periods is distributed as $Z\sim \text{Beta}(2,3)$ and the contamination rate per infected individual is $\beta_s = 2s/n$ . We compare the results obtained by solving the system (3.12) with those obtained by applying the recursion (4.13), and illustrate the fact that the latter are more stable than the former even for the most classical SIR models. In Figure 1 the population size is small ( $n=47$ , $m=3$ ), and we see that both methods deliver the same result. In Figure 2 the population size is higher ( $n=95$ , $m=5$ ). We see that the system (3.12) has become unstable and yields an aberrant distribution, while the recursion (4.13) delivers a sensible result. Finally, in Figure 3 we show that the recursion (4.13) keeps working sensibly for bigger population sizes ( $n=490$ , $m=10$ and $n=980$ , $m=20$ ).