1.2.1. The SIS model

The deterministic SIS model is the following. Let s(t) (respectively i(t)) denote the proportion of susceptible (respectively infectious) individuals in the population. Given an infection parameter λ and a recovery parameter γ, the deterministic SIS model can be written as

Since clearly s(t) + i(t) ≡ 1, the system can be reduced to a one-dimensional ordinary differential equation. If we let z(t) = i(t), we have s(t) = 1 − z(t), and we obtain the ordinary differential equation

It is easy to verify that this ordinary differential equation has a so-called “disease-free equilibrium”, which is z(t) = 0. If λ > γ, this equilibrium is unstable, and there is a stable endemic equilibrium z(t) = 1 − γ/λ.

The corresponding stochastic model is as follows. Let (respectively denote the proportion of susceptible (respectively infectious) individuals in a population of total size N.

Here P_inf (t) and P_rec(t) are two mutually independent standard (i.e. rate 1) Poisson processes. Let us give some explanations, first concerning the modeling, then concerning the mathematical formulation.

Let (respectively ) denote the number of susceptible (respectively infectious) individuals in the population. The equations for those quantities are the above equations, multiplied by N. The argument of P_inf (t) can be written as

The justification for such a rate of infections in the total population is as follows. Each infectious individual meets other individuals in the population at some rate β. The encounter results in a new infection with probability p if the partner of the encounter is susceptible, which happens with probability since we assume that each individual in the population has the same probability of being that partner, and with probability 0 if the partner is an infectious individual. Letting λ = βp and summing over the infectious at time t gives the above rate. Concerning recovery, it is assumed that each infectious recovers at rate γ, independently of the others.

REMARK 1.1.– Let us comment about the fact that we write our stochastic models in terms of Poisson processes. The fact that the infection events happen according to a Poisson process is a rather natural assumption. However, concerning the recovery from infection, our model assumes that the duration of the infectious period follows an exponential distribution. This is not realistic. We are forced to make such an assumption if we want to have a Markov model. We must confess that this assumption is done for mathematical convenience. However, we expect to extend our results to non-Markovian models in forthcoming publications.

Note that there is an equivalent, but slightly more complicated way of writing the Poisson terms, which we now present. Let and denote two mutually independent Poisson random measures on (0, +∞)², with mean measure the Lebesgue measure.

and

Again we have and satisfies

1.2.2. The SIRS model

In the SIRS model, contrary to the SIS model, an infectious who heals is first immune to the illness, he is “recovered”, and only after some time does he lose his immunity and turn susceptible. The deterministic SIRS model can be written as

while the stochastic SIRS model can be written as

These two models could be reduced to two-dimensional models for (respectively

1.2.3. The SIR model with demography

In this model, recovered individuals remain immune forever, but there is a flux of susceptibles by births at rate μN, while individuals from each of the three compartments die at rate μ. Thus, the deterministic model

whose stochastic variant can be written as

REMARK 1.2.– We may think that it would be more natural to decide that births happen at rate µ times the total population. Then the total population process would be a critical branching process, which would go extinct in finite time a.s., which we do not want. Next it might seem more natural to replace, in the infection rate, the ratio which is the actual ratio of susceptibles in the population at time t. It is easy to show that is close to N, so we choose the simplest formulation.

Again, we can reduce these models to two-dimensional models for z(t) = (i(t), s(t)) (respectively by deleting the r (respectively R^N) component.

1.3.1. The stochastic model

The three above-mentioned stochastic models are of the following form.

[1.1]

where are mutually independent standard Poisson processes, and takes its values in

In the case of the SIS model, d = 1, k = 2, h₁ = 1, β₁(z) = λz(1 − z), h₂ = −1 and β₂(z) = γz.

In the case of the SIRS model,

In the case of the SIR model with demography, we can restrict ourselves to d = 2, while

While the above formulation has the advantage of being concise, for certain purposes it is more convenient to rewrite [1.1] using the equivalent formulation already described in the case of the SIS model. Let be mutually independent Poisson random measures on with mean measure the Lebesgue measure, and let We can rewrite [1.1] in the form

[1.2]

in the sense that the joint law of {Z^N, N ≥ 1} is the same law of a sequence of random elements of the Skorohod space D ([0, T]; ), whether we use [1.1] or [1.2] for its definition.

We will now state a few results, without specifying particular assumptions. Those results are valid at least in the case of the three above examples. See Britton and Pardoux (2019b) for details of the proofs, and precise assumptions under which those results hold true.

Concerning the initial condition, we assume that for some z ∈ [0, 1]^d, z_N = [N_z]/N, where is the vector whose i-th component is the integer part of the real number Nzⁱ.

1.3.2. Law of large numbers

We have a law of large numbers

THEOREM 1.1.– Let denote the solution of the stochastic differential equation [1.1]. Assume that the β_j are locally bounded, b is locally Lipschitz and the unique solution of equation [1.3] does not explode in finite time. Then a.s. locally uniformly in t, where {z_t, t ≥ 0} is the unique solution of the ordinary differential equation

[1.3]

The main argument in the proof of the above theorem is the fact that, locally uniformly in t,

1.3.3. Central Limit Theorem

We also have a Central Limit Theorem. Let

THEOREM 1.2.– Assume in addition to the hypotheses of Theorem 1.1 that b is of class C¹. Then, as for the topology of locally uniform convergence, where {U_t, t ≥ 0} is a Gaussian process of the form

[1.4]

where are mutually independent standard Brownian motions.

1.3.4. Large deviations and extinction of an epidemic

We denote by the set of absolutely continuous functions from [0, T] into For any denote the (possibly empty) set of functions c ∈ such that c_j (t) = 0 a.e. on the set {t, β_j(φ_t) = 0} and

We define the rate function

where as usual the infimum over an empty set is +∞, and

with g(ν, ω) = ν log(ν/ω) − ν + ω. We assume in the definition of g(ν, ω) that for all ν > 0, log(ν/0) = ∞ and 0 log(0/0) = 0 log(0) = 0. It is not hard to verify that I_T(φ) =0 if and only if φ solves the ordinary differential equation [1.3]. I_t(φ) can be interpreted as an energy needed for letting φ deviate from being a solution of [1.3].

The collection Z^N obeys a large deviations principle, in the sense that

THEOREM 1.3.– For any open subset

For any closed subset

where for any z

A slight reinforcement of this theorem allows us to conclude a Wentzell–Freidlin type of result. Wentzell and Freidlin have studied small random perturbations of an ordinary differential equation like [1.3] (see Freidlin and Wentzell (2012)). One of their main results is to compute, asymptotically, the time needed for a small random perturbation of such an equation to drive the solution outside of the basin of attraction of a stable equilibrium. The theory has been originally developed for Brownian perturbations. Here we give a statement of the same type, for a Poissonian perturbation. In what follows, we assume that the first component of (respectively z(t)) is (respectively i(t)). Assume that the deterministic ordinary differential equation [1.3] has a unique stable equilibrium z* whose first component satisfies We define

Let now

We have the

THEOREM 1.4.– Given any η > 0, for any z with z₁ > 0,

Moreover, for all η > 0 and N large enough,

It is important to evaluate the quantity Note that it is the value function of an optimal control problem. In case of the SIS model, which is one-dimensional, we can solve this control problem explicitly with the help of Pontryagin’s maximum principle¹, see Pontryagin et al. (1962) or for a concise introduction adapted to this application section A.6 in Britton and Pardoux (2019b), and deduce in that case that For other models, we can compute numerically the value of for each given value of the parameters.

1.4.1. CLT and extinction of an endemic disease

Consider the SIR with demography.

We assume that λ > γ + μ, in which case there is a unique stable endemic equilibrium, namely, Following section 4.1 in Britton and Pardoux (2019b), we can study the extinction of an epidemic in the above model using the CLT. We note that the basic reproduction number R₀ (the expected number of infectious contacts by one infectious at the start of the epidemic, i.e. when and the expected relative time of a life an individual is infected, ε, are given by

[1.5]

The rate of recovery γ is much larger than the death rate μ (52 compared to 1/75 for a one week infectious period and 75 year life length) so for all practical purposes, the two expressions can be approximated by R₀ ≈ λ/γ and ε ≈ μ/γ. Denote again by the fraction of the population which is infectious in a population of size N. The law of large numbers tells us that for N and t large, is close to i*. The CLT tell us that converges to a Gaussian process, whose asymptotic variance can be shown to be well approximated by This suggests that for large t, the number of infectious in the population is approximately Gaussian, with mean Ni* and standard deviation Since we expect a Gaussian process with marginal N(0, 1) to hit −2 fairly quickly, we expect that if then the epidemic will stop rather quickly, while if it is not clear that the time of extinction will be of order 1 (as a function of N). This gives a critical population size roughly of the order of

Note that the factor 9 is rather arbitrary. This N_c is rather large since i* is relatively small. Clearly, even if everybody in the population gets ill at some point, being ill one week in a lifetime of 75 years on average gives a small fraction of infectious in the population.

Consider measles prior to vaccination. If we assume that R₀ ≈ 15 and the infectious period is 1 week (1/52 years) and life duration 75 years, implying that we arrive at N_c ≈ 9(3750)²/15 ≈ 8 · 10⁶. Therefore, if the population is at most a couple of million, we expect that the disease will go extinct quickly, whereas the disease will become endemic (for a rather long time) in a population larger than, for example, 20 million people. This confirms the empirical observation prior to vaccination that measles was continuously endemic in the UK, whereas it died out quickly in Iceland (and was later reintroduced by infectious people visiting the country), see Anderson and May (1991).

1.4.2. Moderate deviations

If the CLT allows us to predict extinction of an endemic disease for population sizes under a given threshold N_c, and large deviations gives predictions for arbitrarily large population sizes, it is fair to look at moderate deviations, which describes ranges of fluctuations between those of the CLT and those of the LD. We shall present the moderate deviations approach in the specific case of the SIS model. In other words, our model from now on is one-dimensional, which can be written in the form

We consider the case λ > γ and recall that the unique stable equilibrium of the deterministic model then is We assume that We have

It follows that

Consequently,

[1.6]

We can also rewrite the above stochastic differential equation in the form

[1.7]

A combination of [1.6] and [1.7] yields the existence of a constant C such that

[1.8]

For the bound by we first take the square in [1.6].

We now define, for 0 < α < 1/2,

and deduce from [1.7]

[1.9]

It follows from [1.9] that the map is continuous from D([0, T]) into itself. Here we equip D([0, T]) with the sup norm topology, which makes it a Hausdorff topologic vector space (equipped with the Skorohod topology, D([0, T]) is not a topologic vector space).

We are interested in the large deviations of which means moderate deviations of Z^N – z*. Note that the deviations of N^α(Z^N – z*) in case α = 1/2 are analyzed by the central limit theorem, and in case α = 0 by the large deviations. So with 0 < α < 1/2, we are clearly in a regime here which is intermediate between the CLT and LD, which is called the regime of moderate deviations.

We first note that the LD of N^α(Z^N – z*) will be deduced from those of thanks to the contraction principle, see for example, Theorem 4.2.1 in Dembo and Zeitouni (1998). So we essentially have to analyze the LD of In fact, we will proceed in three steps. In the first step, we shall analyze the large deviations of

at the speed N^2α–1, or in other words, the moderate deviations of

The second step will consist of showing that have the same behavior as regards large deviations. Finally, the third step will consist of applying the contraction principle, in order to deduce the LD of N^α(Z^N – z*).

1.4.2.1. Step 1: moderate deviations of

We shall use the notations a_N = N^2α–1 and Given a signed measure ν on [0, T], we write

for the logarithmic moment generating function of

The crucial step of our derivation is the

PROPOSITION 1.1.– For any signed measure ν on [0, T], as N → ∞,

where

PROOF.– We first rewrite in the form

here β₀(z) = λz(1 − z) and β₁(z) = γ_z and the processes {M_{j, , j} = 0, 1; 1 ≤ ≤ N} are i.i.d. compensated standard Poisson processes. We now have

as N → ∞,, where we have used the notation □

The next step consists of establishing exponential tightness of the laws of in the sense that

PROPOSITION 1.2.– For any R > 0, there exists a compact set such that

The proof of this Proposition essentially follows the lines of the proof of exponential tightness in section 4.2.4 of Britton and Pardoux (2019b).

We now define the Fenchel–Legendre transform of Λ. Recall that we equip D([0,T]) with the supnorm topology. For each

From Proposition 1.1 and Proposition 1.2 combined with an approximation of by a piecewise linear continuous process (see Pardoux (2019) for the details), we deduce from Corollary 4.6.14 from Dembo and Zeitouni (1998) the following.

THEOREM 1.5.– The sequence satisfies the Large Deviation Principle in with the convex, good rate function Λ* and with speed a_N, in the sense that for any Borel subset

Let us compute Λ*. With the notation s ∧ t := inf(s, t),

It is easily seen that Let now such that φ(0) = 0. The gradient of the map can be written as

We look for ν* such that this gradient equals 0. This implies that

From those identities, combined with φ(0) = 0, we deduce that

1.4.2.2. Step 2: moderate deviations of

We want to show in this step that satisfies exactly the same large deviations result as This will follow if we prove that satisfies Proposition 1.1 (with the same expression in the limit) and Proposition 1.2. Let us state a property that allows us to conclude that satisfies Proposition 1.1 with the correct limit.

PROPOSITION 1.3.– For any C > 0, as N →∞,

[1.10]

We first prove

COROLLARY 1.1.– Given Proposition 1.1, if Proposition 1.3 holds true, then for any signed measure v on [0, T], as N →∞,

PROOF.– For any δ > 0, we deduce from Hölder’s inequality

so that, if we combine Proposition 1.1 and Proposition 1.3, we deduce that

and letting δ − 0, we conclude that

For the inequality in the other direction, we note that, by similar arguments,

so that

hence, letting δ → 0 we conclude that

□

Before we prove Proposition 1.3, we first need to establish a technical Lemma.

LEMMA 1.1.– Let be a standard Poisson random measure on and – dtdu the associated compensated measure. If ϕ is an ₊-valued predictable process, such that has exponential moments of any order, and then there exists a constant C such that for any 0 ≤ t ≤ T,

PROOF.– Consider with b ≥ 0 the process

[1.11]

It follows from Itô’s formula that

It follows from Lemma 1.2 below that is a martingale. Hence, e^X is a martingale, if b = (e^a − 1 − a), a submartingale if we replace = by <, and a supermartingale if we replace = by >. Hence if b ≥ (e^a − 1 − a), Now, first using Doob’s L² inequality for submartingales, and later Cauchy’s inequality, we have

If 2b = e^2a − 1 − 2a, the first factor on the last right-hand side equals 1. □

In order to complete the proof of Lemma 1.1, we still need to establish

LEMMA 1.2.– The process ϕ satisfying the same assumptions as in Lemma 1.1, and X_t being given by [1.11], is a martingale.

PROOF.– It is plain that M_t is a local martingale, whose predictable quadratic variation is given as

All we need to show is that the above quantity is integrable. It is clearly a consequence of the assumption in case a < 0. In case a > 0, the second factor of the right-hand side has finite exponential moments, so is square integrable, and all we need to show is that

[1.12]

Using Itô’s formula we have

It is easy to conclude that It follows from Cauchy–Schwartz that

and the result follows from our assumption on ϕ. □

We now turn to the

PROOF OF PROPOSITION 1.3 We note that

Proposition 1.3 will follow from the fact that for any C > 0, as N →∞,

[1.13]

[1.14]

[1.15]

We shall prove [1.13] and [1.15]. The proof of [1.14] is quite similar to that of [1.13].

STEP 1: PROOF OF [1.13] It suffices to consider one of the terms in the sum over j, and we suppress the index j for simplicity. We note that

It is not hard to see that we can treat each of the two terms on the right separately, and treat only the first term, with the treatment of the second one being quite similar. We note that there exists a compensated Poisson process on ₊ M such that this first term can be rewritten as

We need to estimate If we decompose the signed measure ν as the difference of two measures as follows ν = ν₊ − ν₋, we again have two terms, and it suffices to treat one of them, say ν₊. Of course, it suffices to treat the case where Since the positive constant C is arbitrary, without loss of generality, we can assume that ν₊ is a probability measure on [0, T]. It is then clear that

We choose a new parameter and we split the expression whose expectation needs to be estimated in two terms.

[1.16]

We now estimate the first term on the right-hand side of [1.16]. In this case, we define the stopping time

and note that

Consequently, the expectation of the first term on the right of [1.16] is bounded from above by

where the first inequality follows from Lemma 1.1, and the second one exploits the Lipschitz property of β. Consider now the second term on the right-hand side of [1.16].

where the second inequality follows from Lemma 1.1 and the boundedness of β. For the second factor in the last expression, we need to consider

Both probabilities are estimated in a similar way. By an exponential estimate,

[1.17]

for N large enough. Finally, the expectation of the second term of the right-hand side of [1.16] is bounded by expjciN with c₁, c₂ > 0, and

From the inequality log(a + b) ≤ log(2) + log(sup(a, b)), for N large enough,

which establishes [1.13].

STEP 2: PROOF OF [1.15] We, in fact, must prove that for any C > 0, as N →∞,

In this proof, C will denote a constant whose value may change from line to line. We now introduce a new process, where

the event

and the stopping time

where the constant b will be chosen below. From [1.8], the fact that and Cauchy–Schwartz,

[1.18]

[1.19]

we take the limit successively in the two terms of the above right-hand side.

STEP 2a: estimate of [1.18] We have

It remains to be noted that

and

for some positive constant C, so that finally there exist two positive constants C₁ and C₂ such that, for N large enough,

STEP 2b: estimate of [1.19] Since Y^N is a martingale, it is clear that the process

is a submartingale. Consequently, from Doob’s L² submartingale inequality,

[1.20]

Consider first the first factor on the right-hand side of [1.20]. We have

and the result follows from [1.13] and [1.14].

We finally consider the second term on the right-hand side of [1.20]. We have

Hence, the second term on the right of [1.20] satisfies

[1.21]

We now write

with

We have proved that the first factor on the right-hand side of [1.21] remains bounded, as N → ∞. Next consider the second term on the right-hand side of [1.21]. We have

It follows that the second factor in [1.20] is bounded from above by

where C1 and C2 are two positive constants. This last expression is bounded say by 2, as soon as N is large enough. □

1.4.2.3. Step 3: moderate deviations of Z^N − z*

Recall that is the image of by the mapping from D([0, T]) into itself, which is continuous if we equip D([0, T]) with the supnorm topology, defined by:

Note also that the above mapping is a bijection. The following result is then a consequence of Corollary 1.1 and the contraction principle.

THEOREM 1.6.– The collection of processes satisfies a large deviation principle with the good rate function

REMARK 1.3.– While the LD rate function of a Poisson-driven stochastic differential equation is very different from the rate function of LDs for its Brownian-driven diffusion approximation, the rate function for moderate deviations of Poisson-driven stochastic differential equations is identical to that of LDs for its Brownian-driven diffusion approximation.

1.4.2.4. Wentzell–Freidlin theory and extinction of an epidemic

We want to conclude from the Wentzell–Freidlin theory an estimate of the time needed for to make a deviation of −c, i.e. to go from 0 to −c, which, for the value of N such that z* = cN^-α, means the time for to hit 0. In this case, we first compute

An application of Pontryagin’s maximum principle, see Pontryagin et al. (1962), yields

Using the same arguments as in Kratz and Pardoux (2018) and Britton and Pardoux (2019b), we then deduce from Theorem 1.6

THEOREM 1.7.– Let For any δ > 0,

Moreover,

Recall that In the CLT regime, while in the LD regime,

Let us now compute the corresponding critical population size. is the order of magnitude of the time needed for to make a deviation of size cN^-α. This is sufficient to extinguish an epidemic, provided i* is of the same order, so that the corresponding critical size is N_α ~ (1/i*)^1/α, that is roughly the CLT critical population size raised to the power 1/2α. In the case of the SIR model with demography for measles, the CLT critical population size is of the order of a few million; thus, for example, with α = 1/3, we go from 10⁶ to 10⁹, i.e. a few billion, which is the order of magnitude of the biggest countries, i.e. China and India.