key: cord-0229801-wzon0u1g
authors: Forien, Raphael; Pang, Guodong; Pardoux, Etienne
title: Recent Advances in Epidemic Modeling: Non-Markov Stochastic Models and their Scaling Limits
date: 2021-06-15
journal: nan
DOI: nan
sha: fa829d171662a54a2dbb77083383ab275180a688
doc_id: 229801
cord_uid: wzon0u1g

In this survey paper, we review the recent advances in individual based non--Markovian epidemic models. They include epidemic models with a constant infectivity rate, varying infectivity rate or infection-age dependent infectivity, infection-age recovery rate (or equivalently, general law of infectious period), as well as varying susceptibility/immunity. We focus on the scaling limits with a large population, functional law of large numbers (FLLN) and functional central limit theorems (FCLT), while the large and moderate deviations for some Markovian epidemic models are also reviewed. In the FLLN, the limits are a set of Volterra integral equations, and in the FCLT, the limits are stochastic Volterra integral equations driven by Gaussian processes. We relate our deterministic limits to the results in the seminal papers by Kermack and McKendrick published in 1927, 1932 and 1933, where the varying infectivity and susceptibility/immunity were already considered. We also discuss some extensions, including models with heterogeneous population, spatial models and control problems, as well as open problems.

Had this paper been written in 2019 or before, we would have needed, in order to drive the interest of the reader, to recall the major pandemics of the previous centuries: the black plague pandemic which killed between 30% and 50% of Europe's population in the 14th century, the plague epidemic which killed almost half of the population of Marseille and a quarter of the population of Provence in 1720, the so-called Spanish flu which killed between 50 and 100 million people, not forgetting HIV/AIDS, malaria and tuberculosis, which together killed more than 3 million humans in 2011. This long list is not without a silver lining, thanks to the huge success of vaccination, which in particular has reduced the number of deaths due to measles by 94% and has permitted the eradication of smallpox in the late 20th century.

However, in 2021, everyone has heard about infectious diseases, the basic reproduction number, herd immunity and the importance of vaccination. This is a result of the Covid-19 pandemic, which started at the end of 2019 in China, and by spring 2020 had hit Europe and North America, filling intensive care units in every country. During the spring of 2020, many countries implemented drastic lockdown measures, which were decided after the leaders of those countries had learned the predictions of mathematical models about the number of deaths that the epidemic was likely to cause, if no such measures were taken. One year later, in the spring of 2021, most wealthy countries are striving to vaccinate a large proportion of their population, in the hope that they can get rid of the epidemic, or at least lower the pressure on hospitals to a manageable level.

The aim of this survey paper is not primarily to present all those notions, or to review all the efforts aimed at modelling this particular pandemic, but rather to highlight some recent progress in epidemic modelling to which the authors of this paper have contributed.

The use of mathematics and mathematical models as a tool to understand and control the propagation of infectious diseases has a long history. Around 1760, Daniel Bernoulli, a member of a famous family of mathematicians, who had also been trained as a physician, exploited a mathematical model in order to convince his contemporaries of the advantage of inoculation (the ancestor of vaccination) against smallpox, discussing already the balance between the benefit of inoculation and the associated risk, a question which is much debated these days. The foundations of modern epidemic modelling was mainly the result of efforts of physicians, rather than mathematicians. During the second half of the 19-th century, the Russian physician P. D. En'ko was probably the first scientist who created a chain binomial stochastic model of epidemics, similar to the better known Reed-Frost model, which was formulated in 1928, but published only in the 1950's.

Most modern mathematical epidemic models are formulated as deterministic compartmental models. Around 1910, Ross introduced the concept of the basic reproduction number R 0 , and argued that malaria would stop if the proportion of mosquitos to humans were maintained under a certain threshold. His arguments, which were based on the understanding of the large time behaviour of dynamical systems, were not accepted by many of his contemporaries, who claimed that malaria would continue as long as mosquitos would be present.

The description of the transmission of communicable diseases via compartmental models was pioneered by Kermack and McKendrick in a series of three papers, published in 1927, 1932 and 1933, see [51, 52, 53] . Their models were very refined. In their 1927 paper, they consider infection-age dependent infectivity (i.e., the infectivity of an infectious individual depends upon the time elapsed since he/she was infected), as well as infection-age dependent recovery rate which, as we shall explain below, corresponds to the fact that the duration of the infectious periods can be very general (any absolutely continuous distribution). In one section of the paper, they consider the case of constant rates, i.e., constant infectivity, and constant recovery rate, the latter imposing that the law of the duration of the infectious period be an exponential distribution. In this particular case, the deterministic model is a system of ordinary differential equations (ODEs), instead of the more complex system of Volterra type integro-differential equations in the general case. It is rather clear that the general model considered in [51] can be made much more realistic by a proper choice of the coefficients adapted to each particular situation (both to a specific illness and to a specific society with its interactions) than the particular case of constant rates. However, almost all epidemic models which were considered since 1927 treat the special case of ODE models, and when the seminal paper [51] is quoted, most of the time reference is made only to the special case of constant rates.

Furthermore in their second and third papers [52, 53] , Kermack and McKendrick considered the loss of immunity (and study the endemic situations), again with a very realistic point of view: they consider that loss of immunity is not sudden, but progressive. It is hard to find a recent work which adopts that point of view (one exception is [46] ). In their 1932 paper [52] , they also pioneer the use of PDE models for the description of epidemic models with infection-age dependent infectivity and recovery rate and recovery-age dependent level of immunity.

The goal of the present survey paper is twofold. First, we want to draw the attention of the readers to the complex models of Kermack and McKendrick, the more classical ODE models being in our opinion a rather unrealistic approximation of the former, which should be used only when both have a similar behaviour. We shall discuss that point below. Second we want to derive the deterministic models (whether ODEs, integral equations or PDEs) as law of large numbers limits, in the asymptotic of a large population, of individual based stochastic models. This can be seen as an analogue of many recent works which establish certain equations of physics as limits of stochastic particle systems, as the number of particles tends to infinity (see in particular the book by Kipnis and Landim [54] , and the references therein). We shall also discuss the difference between the stochastic and the deterministic models, via the central limit theorem, moderate and large deviation principles.

1.1. Literature review. The Markov models and their limiting ODE models have long been the standard tools to study epidemics, see the recent survey [23] and monographs [5, 7, 60, 16, 21] . We refer the readers to the above for the related literature. The functional law of large numbers (FLLN) and functional central limit theorem (FCLT) results for Markov models can be found in [23] , which use the standard Poisson random measure representations and martingale convergence arguments in [30] . Most relevant results concerning the large deviation and moderate deviation principles in Section 2 can be found in the recent works [70, 56, 71, 23, 69] .

Since the seminal work by Kermack and McKendrick [51, 52, 53] , Volterra integral equations have been developed, without proving an FLLN rigorously, in various epidemic models with general infectious periods, see, e.g., [20, 25, 27, 43, 79, 31, 21] . Also, PDE models have been developed, again without proving an FLLN (with one exception [36] ), for various epidemic models, Markovian or non-Markovain, see, e.g., [77, 47, 59, 44, 58, 89, 50, 39] . This survey does not focus on these various deterministic Volterra integral equations and PDE models, but on individual based infectious models for which such equations arise as scaling limits, in particular, our recent works in [64, 33, 32, 65, 66, 67, 34] .

Let us also provide a brief literature review on some other existing methods and results on these general models. Sellke [75] developed an approach, the so-called "Sellke construction", to find the distribution of the number of remaining uninfected individuals in an SIR/SEIR epidemic model with a large population. Ball [9] developed a unified approach using a Wald's identity to find the distribution of the total size and the total area under the trajectory of infectious individuals. Ball [10] further developed this approach to study multi-type epidemic models. Barbour [14] proved limit theorems for the distribution of the duration from the first infection to the last removal in a closed epidemic. The LLN and CLT results concerning the final number of infected individuals are presented in [23] . Non-Markovian SIR epidemic models were studied as piecewise Markov deterministic processes using the associated martingales in [24, 40] to analyze the distribution of the number of survivors of the epidemic and the population transmission number and the infection probability of a given susceptible individual.

We now review the existing literature on functional limit theorems for the non-Markovian epidemic models preceding our recent works. Wang [82, 83, 84] proved FLLN and FCLT for some age and density dependent stochastic population models, including the SIR model which has an infection rate depending on the number of infectious individuals and allows an initial condition with infection-age dependent infectivity. The strategy of the proof in [82, 83, 84] is different from ours, does not make use of Poisson random measures, and also assumes a C 1 condition on the distribution of the infectious period for the FCLT. Some asymptotic properties of the limiting deterministic integral equations were studied [85] . Reinert [74] used Stein's method with a generalized Sellke construction to prove a deterministic limit (LLN) for the empirical measure describing the system dynamics of the generalized SIR model with the infection rate dependent upon time and state of infection; however, no FCLTs have been establish using her approach. In developing a PDE model to study the Covid-19 pandemic, [36] establishes the PDE model as a a law of large numbers limit of stochastic individual based models.

The models and main results in our recent works are reviewed in this article, which are briefly described in the next subsection.

Organization of the paper. The organization of the paper is as follows. In Section 2, we consider the "special case" of constant rates, which leads to a Markov individual-based stochastic model, an ODE model in the FLLN and a diffusion model in the FCLT. We also discuss the large deviation and moderate deviation results in the Markov models, which draws upon [69] . In Section 3, we consider the case of a constant infectivity, but a general law of the infectious period, or equivalently an infection age dependent recovery rate, in which case the stochastic model is non-Markov, and the deterministic LLN limiting model is a system of integro-differential equations, i.e., a system of equations with memory. The stochastic limiting processes in the FCLT are Gaussian-driven Volterra integral equations. This section draws on our first paper in the series in [64] . We also describe the discrete spatial model, i.e., a multi-patch non-Markov model with constant infectivity, where individuals may migrate from one patch to another and individuals may be infected locally within each patch or from some distance (which can be thought of as the result of infections taking place during short stays of individuals outside of their current patch). This part draws upon [65] . In Section 4, we add the infection-age infectivity. At the level of the stochastic model, we assume that the infectivity functions of the various individuals are i.i.d. copies of a given càdlàg random function. It turns out that only the mean of this random function appears in the limiting LLN deterministic model, which, as we shall see, is precisely the model introduced in [51] . We then present the limiting Gaussian-driven stochastic integral equations obtained in the FCLT. We also discuss the use of these models to model the Covid-19 epidemic. These results are proved in [33, 32, 66] . We shall then discuss the PDE point of view of the same model in Section 5. In addition to the total infectivity process, we use a stochastic process that tracks the number of infected individuals at each time that have been infected for a certain amount of time. The LLN limit is again a system of integro-differential equations, while the density of the more detailed process with respect to the elapsed infectious time have a PDE representation that coincides with the equations of [52] if the distribution of the infectious period is absolutely continuous. In addition, we obtain a PDE for the models with a deterministic infectious period. These PDE results are established in [67] . In Section 6 we discuss the case where the infectivity depends upon the age of infection, the duration of the infectious period has a general distribution, and the loss of immunity is a random function of the time elapsed since recovery. Again, those varying immunities of the various individuals are i.i.d. copies of a given random function. It turns out that in this case, the limiting deterministic model involves the whole distribution of this varying immunity function (and not just its expectation), and in general our LLN model is very different from the model introduced in [52] , unless the loss of immunity is described by the same deterministic function of the time elapsed since the random recovery time for all individuals. This result is proved in [34] . Finally, in Section 7, we discuss various extensions, including models with heterogeneous population, spatial models and control problems. We then discuss open problems.

Starting with Section 3, this paper describes recent results obtained by the authors, as the output of a research effort which started in January 2020. Needless to say, our program has not yet been completed. In particular, the moderate and large deviations results have so far been obtained only in the Markov case, and little has been done until now on spatial model outside the Markov case. We nevertheless believe that it is now a good time to put together a series of results which both relate stochastic and deterministic results, and insist upon the rich and complex models which Kermack and McKendrick introduced almost a century ago, and have been unfortunately largely neglected and/or forgotten.

1.3. Basic vocabulary of epidemic models. Compartments. In a compartmental model, each individual belongs to one of the following compartments:

-S denotes the compartment of susceptible individuals: those who are not infected, but are susceptible to the disease, which means that they might get infected if they meet an infectious individual. -E denotes the compartment of exposed individuals, who are infected, but not yet infectious.

-I denotes the compartment of infectious individuals, who are infected, and able to transmit the disease to susceptible individuals. Note that when considering the models with varying infectivity, we shall denote by I the compartment of infected individuals, whether they are exposed or infectious. Their infectivity is ≥ 0, they are infectious when it is > 0. -R denotes the compartment of removed or recovered individuals. Those who have been infected, and have recovered from the disease. They are neither infected, nor susceptible, and are immune to the disease. Often one includes in that compartment those who died from the disease. In the model with varying immunity/susceptibility which we shall consider in Section 6, we will merge the S and R compartments into the S compartment, each individual having a varying susceptibility. Whenever that susceptibility is 0, he/she cannot be infected, as when he/she is in the R compartment. The above compartments are the most commonly used ones, but certain models consider other compartments, e.g., V for vaccinated. Also, in particular concerning the COVID-19, some authors have decomposed the I compartment into several ones, distinguishing, e.g., the symptomatic and the asymptomatic infectious individuals, creating a compartment for those who are hospitalized, and another one for those in Intensive Care Units. We shall see below that our formalism allows one to reformulate very general models with many compartments as non-Markov models with only a few compartments.

Various types of models. The most classical model is the SIR model. In this model, when a susceptible individual is infected, he/she leaves the S compartment and enters the I compartment. While in that compartment, he/she is infectious and can infect susceptible individuals. After some time, the infectious individual recovers, moves to the R compartment, and stays there for ever. This means that he/she is immune and cannot be infected a second time.

A variant of the SIR model is the SEIR model, in which infected individual first enter the E compartment. While in that compartment, the individual is not infectious. He/she becomes infectious when entering the I compartment. As in the SIR model, the individual eventually recovers when entering the R compartment.

Next we have the models where the recovered individuals lose their immunity after some time, and become susceptible again. The simplest of those is the SIS (or the SEIS) model, where upon recovery the individual becomes susceptible again, i.e., there is no period of immunity. Another class of such models is the SIRS (or the SEIRS) model, where upon recovery the individual first stays for some time in the R compartment, where he/she is immune and cannot be infected. Later he/she loses immunity, and becomes susceptible again.

In the simplest models, the size of the population is fixed. This makes sense if the epidemic is considered over a time interval during which there are not many births and deaths, and the deaths due to the epidemic are possibly included in the R compartment. Quite a few models however include the demography, and the total size of the population is allowed to fluctuate. The latter are called S(E)IR(S) models with demography.

As already explained, contrary to the traditional approach, if we consider an infection-age dependent infectivity, we need not distinguish the compartments E and I (while in E, the infectivity is zero), and if we consider a recovery age dependent susceptibility, we need not distinguish the compartments R and S (while in R, the susceptibility is zero). We shall also see that the non-Markov models / models with memory allow to have a precise description of the propagation of the disease, without increasing the number of compartments, as is commonly done with Markov / ODE models. R 0 . A fundamental concept of infectious disease modeling is the so-called basic reproduction number, denoted R 0 , which is the mean number of susceptible individuals whom an infectious individual infects during its infectious period, at the beginning of the epidemic (i.e., while essentially all members of the population are susceptible).

Note that when a significant fraction of the population has been hit by the disease, the mean number of susceptible individuals whom an infectious individual infects during its infectious period will be different, and is sometimes called the effective reproduction number, which depends upon time, since it depends upon the evolution of the epidemic. More precisely, ifS(t) denotes the proportion of susceptible individuals in the population, R ef f (t) = R 0 ×S(t). If R ef f (t) < 1, the epidemic regresses and eventually goes extinct. For that to occur, we need to haveS(t) < R −1 0 , i.e., the proportion of immune individuals should be greater than 1 − R −1 0 . In such a situation, one says that herd immunity has been achieved.

In this section, we shall discuss stochastic Markov models, and their limiting ODE models. In this case, the proof of the LLN is rather easy, and has been known for a long time. Also, the fluctuations of the stochastic model around its law of large numbers limit has been fully studied. Not only do we have a Central Limit Theorem, but also moderate and large deviations have been studied in this simpler case.

2.1. The SIR Markov model and its LLN limit. We shall follow the recent presentation in Section 2 of [23] . Let us first describe the stochastic individual based model. We shall use the notions of a Poisson process and a Poisson Random Measure, two notions which are introduced in Subsection 8.1 of the Appendix below.

Suppose that we have a population of fixed size N , which is distributed in the three compartments S, I and R. Let S N (t) (resp. I N (t), resp. R N (t)) denote the number of susceptible (resp. infectious, resp. recovered) individuals at time t. We have S N (t) + I N (t) + R N (t) = N . We assume that each infectious individual meets others at rate β. If the encountered individual is susceptible, which at time t happens with probability S N (t)/N (since we make the homogeneity assumption that the individual who is met is chosen uniformly in the population), then the encounter results in a new infection with probability p. Hence, if we use the notation λ = β × p, the rate at which one particular infectious individual infects susceptibles is λS N (t)/N . Hence the total rate of new infections in the population at time t is

This means that the number of new infections on the time interval [0, t] is

where P inf (t) is a standard Poisson process. The fact that this process at time t involves P inf up to a random time creates sometimes technical problem. Therefore we shall use in further sections an alternative description which we now introduce. Given a standard Poisson random measure Q on R 2 + , an alternative equivalent description of the counting process of infections is

The crucial point now is that in the Markov model considered in this section, we assume that the duration of the infection period is Exp(γ) (the durations for various individuals are independent, and independent of the rest of the population), where γ −1 is the mean duration of this infectious period. The end of the infectious period of a given individual is thus the first jump time of a rate γ Poisson process. Since the sum of mutually independent Poisson processes is a Poisson process with rate the sum of the rates, the counting process of the number of recoveries on the interval [0, t] (i.e., the number of jumps from the I to the R compartment) is

where P rec (t) is a standard Poisson process, independent of P inf (t). Finally, we obtain the following system of stochastic differential equations for the evolution of the numbers of susceptible, infectious and recovered individuals:

(2.1) Remark 2.1. Value of R 0 It is easy to compute R 0 in the present model. "At the beginning of the epidemic", means "while N −1 S N (t) = 1". In that case, each infectious individual infects a mean λ susceptible individuals per time unit. The mean duration of the infectious period is 1/γ.

. We need to formulate one assumption.

Assumption 2.1. We assume that as N → ∞,

Remark 2.2. If eitherS(0) = 0 orĪ(0) = 0, there would be no epidemic in the limiting LLN deterministic model. Typically an epidemic starts with a small number of initially infectious individuals, which is not of the order of N . The description of the first phase of the epidemic, until the number of infectious reaches a "positive fraction of N " must be done by a stochastic model, the deterministic model becomes valid once a significant fraction of the total population is infectious. The stochastic model at the start of the epidemic can be well approximated by a branching process. We shall explain this below in Section 4.

In the next LLN result, we shall denote by D the Skorokhod space of càdlàg real-valued functions on R + , endowed with the usual topology. The reader is referred to Subsection 8.3 in the Appendix below for its definition and properties. 

with the initial condition specified by Assumption 2.1.

Note that (S(t),Ī(t),R(t)) are, respectively, the proportions of susceptible, infectious and recovered individuals in the limit as N → ∞.

Proof. Let us consider the proportions in the three compartments, i.e., we divide equation (2.1) by N , and define M inf (t) := P inf (t) − t, M rec (t) := P rec (t) − t. We obtain, with the notation

Indeed, the pointwise convergence follows directly from the classical law of large numbers, and then the uniform convergence from the fact that the function t → N −1 P (N t) is increasing, and converges to the continuous function t → t, thanks to the second Dini Theorem. This, combined with Assumption 2.1, allows one to take the limit in (2.3), and deduce (2.2). More details can be found in Section 2.2 of [23] .

2.2. The Markov SIR model : Central Limit Theorem. Let us now rescale the differences between the proportions in the N model and the limiting proportions. We define

In order to obtain a limit of the above processes, we need to formulate an assumption concerning (Ŝ N (0),Î N (0),R N (0)).

There exists a random vector (Ŝ(0),Î(0),R(0)) such that

We have the following FCLT. 

where (Ŝ(t),Î(t),R(t)) is the unique solution of the following linear SDE: 

where B inf (t) and B rec (t) are two mutually independent standard Brownian motions, which are globally independent of (Ŝ(0),Î(0),R(0)). If (Ŝ(0),Î(0),R(0)) is Gaussian, then (Ŝ(t),Î(t),R(t)) is a Gaussian process.

Note that the notion of a Brownian motion is defined in Section 8.2 in the Appendix below.

Proof. We take the difference between (2.3) and (2.2), and multiply by √ N . Hence

The result now follows from Theorem 2.1, the next Lemma and rather standard arguments, see for the details Section 2.3 in [23] . 

where B(t) is a standard Brownian motion.

Proof. Note that M(t) := N −1/2 M (N t) is a square integrable martingale, whose associated increasing process is given as M t = t. Hence tightness in D follows readily by the criterion from Proposition 8.1 in the Appendix below. It thus suffices to show that for any n ≥ 1, any 0 = t 0 < t 1 < · · · < t n ,

By independence of the increments of both M and B, it suffices to show that for any t > 0,

. This is easily verified by a characteristic function computation. Indeed, for any u ∈ R,

2.3. Markovian SIS and SIRS models, and SIR model with demography. In the SIR model, the number of susceptibles who can be infected is limited, and therefore, the epidemic goes soon or later to an end. However, there are several models where there is a constant flux of susceptibles, which allow the establishment of an endemic disease. Let us describe three such models.

2.3.1. The SIS model. In this model, contrary to the SIR model, when an infectious individual recovers, he/she becomes susceptible again. There is no immunity. The stochastic model reads:

and the LLN limiting deterministic model reads:

Note that exploiting the identities S N (t) + I N (t) = N ,S(t) +Ī(t) = 1, we can write in fact equations for I N (t) andĪ(t) only, which read:

Again in this model R 0 = λ/γ. If R 0 ≤ 1, the last equation has the unique equilibrium I * = 0, while if R 0 > 1, this disease-free equilibrium is unstable, and there is a stable endemic equilibrium

The SIRS model. In this model, an individual is first removed (i.e., immune) when he/she recovers, but he/she loses its immunity at a given rate ρ. This gives the following stochastic model (which we write for the two quantities S N (t) and I N (t):

where "loim" is an abbreviation for "loss of immunity", and the following deterministic model:

Again R 0 = λ/γ, if R 0 ≤ 1, the only equilibrium is (1, 0), while if R 0 > 1, we have an endemic equilibrium (γ/λ, (1 − R −1 0 )(γ + ρ) −1 ρ). 2.3.3. The SIR model with demography. In this model, the recovered individuals do not lose their immunity, but births produce a constant flux of susceptible individuals. The stochastic model reads:

Remark 2.3. We model the birth as a constant flux at rate µN , instead of µ times the number of individuals in the population, in order to avoid the pitfall of critical branching processes, which go extinct in finite time a.s. As a result, if the individuals in the R compartment die at rate µ as well, the total population in the stochastic model remains close to N . Therefore we approximate the proportion of susceptibles in the population by S N (t)/N .

This time, R 0 = λ γ+µ . If R 0 > 1, the endemic equilibrium reads ((γ + µ)/λ, (1 − R −1 0 )(γ + µ) −1 µ).

Deviations from the law of large numbers and extinction of an endemic disease.

In the above three models, the endemic equilibrium is stable whenever R 0 > 1. This means that the LLN deterministic model, starting from a positive value ofĪ(0), will never go extinct. However, if we consider the stochastic model, it is easily seen that for any N ≥ 1, disease free states are accessible with positive probability. Since moreover they are absorbing, in the stochastic model the epidemic stops soon or later. We would like to know how long we have to wait for this to happen. One approach is to evaluate the time needed for the stochastic system to diverge enough from its deterministic limit, so that I N (t) = 0. For that purpose, we shall exploit three tools that Probability theory gives us, in order to estimate the difference between a stochastic process and its law of large numbers limit, namely the Central Limit Theorem, Moderate Deviations and Large Deviations. Let us first see what the CLT tells us. Consider first the SIR model with demography. In that model, at the endemic equilibrium,

It is clear that γ is much larger than µ (in inverse of years, compare 52 to 1/75, since γ −1 is of the order of 1 week, and µ −1 is of the order of 75 years). Consequently we can consider that R 0 ∼ λ γ and ε := µ γ+µ ∼ µ γ . Now the CLT tells us that I N (t) is approximately Gaussian, with mean N ε(1 − R −1 0 ) and standard deviation N/R 0 (the asymptotic variance ofÎ(t) is close to R −1 0 , see [23] page 62). If N is such that the standard deviation is at least the mean divided by 3, then it is likely that I N (t) will hit zero in time of order 1. This leads to the idea of a critical population size N c given by

In the case of measles, R 0 = 15. With the above approximations for µ and γ, we arrive at N c in the order of a few million. This confirms the empirical observation that, prior to vaccination, measles was continuously endemic in countries like UK, and died out quickly in Iceland (and was then reintroduced by infected visitors).

Remark 2.4. Would we consider the SIS model instead of the SIR model with demography, then we would find a much smaller N c . Indeed, the asymptotic variance in the CLT is about the same, butĪ * is much larger. Indeed, if everyone who recovers becomes susceptible again, we are likely to have a much larger proportion of infectious at equilibrium. In the SIR model with demography, contrary to the situation in the SIS model, those who get infected are infected only once in their life. The ratio between the twoĪ * 's is the above ε.

As explained above, sooner or later the stochastic process I N (t) will hit zero (and then stay there for ever). The CLT allows us to guess for which population sizes extinction is likely to happen in time of order 1. We will now discuss what large deviations tell us on this problem. In the three above examples, we have an R d -valued ODE of the form

which is the LLN limit of a sequence of SDEs of the form

and we have b(z) = k j=1 β j (z)h j . Recall that P 1 , . . . , P k are mutually independent standard Poisson processes. Note that if we rewrite the SDE (2.6) in the form (with M j (t) = P j (t) − t)

from the fact that the last term on the right hand side of this equation tends to 0 as N → ∞, we may regard the SDE (2.6) as a "random perturbation of the dynamic system (2.5)", Freidlin and Wentzell have studied such perturbations as an application of large deviations theory, see [38] . Let us first state what kind of information large deviations give us, concerning the convergence of Z N toward z. We shall state the results without giving the precise technical conditions under which one can establish them, referring the reader to Section 4.2 of [23] for all the technical details and proofs. Consider the function :

We let

One essential property of this functional is that I T ≥ 0 and I T = 0 iff φ solves the ODE (2.5). One may think of I T (φ) as a sort of measure of how much φ differs from being a solution to (2.5) .

Large deviations theory gives us both

Note that, as is to be expected, those results give us information only in case O (resp. F ) does not contain the solution of (2.5) starting from x. Following the ideas of Freidlin and Wentzell, one can deduce from those two statements a rather precise statement about the time taken by the random perturbations to drive the process Z N to a disease free situation. Denote by A the subset of points of R d + which are accessible by our system. Suppose that I N is the first component of Z N . We are interested in the time needed for Z N (t) to reach the subset of R d + where its first component is 0. Let us denote by z * ∈ A the endemic equilibrium, i.e., the point in A such that b(z * ) = 0 and z * 1 > 0, which we assume to be unique. We now define the "quasi-potential" (D T,A stand for D([0, T ]; A)):

For z ∈ A, we define the extinction time of the process Z N as 

Note that V is the value function of an optimal control problem. It can be computed explicitly in case of the SIS model, in which case V = log R 0 − 1 + R −1 0 . Unless V is quite small, we expect exp(V N ) to be very large.

In between the CLT and Large Deviations, we have the theory of "moderate deviations". Let us explain what we can learn from this theory, concerning our problem. We will only give a brief sketch of the ideas, referring the reader to [69] for the details.

Large Deviations discusses the probability of observing deviations from the LLN of the order of 1, as well as the time we have to wait for observing such deviations. The CLT predicts deviations from the LLN of the order of N −1/2 . Moderate Deviations discusses deviations of the order of N −α , for some 0 < α < 1/2: both the probability of observing such deviations, and the time we have to wait to see such deviations. This should allow us to predict extinction in less time that what Large Deviations predicts, with a critical population size larger than that associated to the CLT.

Since we want to discuss deviations from the endemic equilibrium z * of the order of N −α , let us consider the process Z N z (t), starting from Z N z (0) = z * + N −α z, where z ∈ R d is arbitrary. We want to study the Moderate Deviations of that process, which amounts to studying the Large Deviations

The case α = 1/2 is covered by the CLT, the case α = 0 by Large Deviations. Moderate Deviations, fills the gap between those two regimes. If the population size N is such that z * 1 is of the order of N −α , for some 0 < α < 1/2, Moderate Deviations will predict extinction in time of the order of exp[N 1−2α V a ], with a = N α z * 1 .

3. Non-Markov and integro-differential models 3.1. The SIR model. In this section, we will still assume that the infectivity is constant, but we shall let the infectious period have a general probability distribution. The way Kermack and McKendrick assumed a general distribution for the infectious period in their 1927 paper [51] was to choose an infection age dependent rate of recovery. Let us first show that this formulation covers all absolutely continuous distributions for the duration of the infectious period. Indeed, let X be an R + -valued random variable, F (t) = P(X ≤ t) its distribution function, and F c (t) = 1 − F (t). If F has a density f (t) (i.e., f (t) = F (t)), then we define its hazard function as the quantity γ(t) := f (t)/F c (t). Let P (t) be a standard Poisson process. Then the law of X coincides with that of the first jump of the counting process P t 0 γ(s)ds , which follows from the following computations.

In other words, by choosing an infection age dependent recovery rate, Kermack and McKendrick allowed a general absolutely continuous distribution for the duration of the infectious period. We shall use a different formulation, and allow a completely arbitrary distribution for the infectious period.

Concerning the infection process, we have the same infection rate as in Section 2.1, namely

Let A N (t) denote the cumulative counting process of newly infected individuals on the time interval (0, t]. We have, as above,

Clearly the following balance equations hold

To each newly infected individual i ∈ N, we associate a random variable η i to represent the infectious duration. We assume that η i 's are i.i.d. with a cumulative distribution function (c.d.f.) F , and let F c := 1 − F . For each initially infectious individual j = 1, . . . , I N (0), let η 0 j be the remaining infectious period. We also assume that η 0

1. An initially infected individual is thought of as having been infected at some time τ 0 < 0. Assuming that the law of the duration of the infectious period of this individual is F , the probability that he/she is still infectious at some time t > 0, given the time of infection, equals F c (t − τ 0 )/F c (−τ 0 ), the conditional law of still being infectious at time t, given that he/she was still infectious at time 0. In the case where F is exponential, this is exactly F c (t), hence there is no reason to choose an F 0 different from F . But in all cases where F is not the exponential distribution, it is quite natural to choose F 0 = F .

Denoting by τ N i , i ≥ 1 the successive jumps times of A N , the dynamics of I N (t) can be described by

Define the fluid-scaled processX N := N −1 X N for any process X N . We make the following assumption on the initial quantities. 

in probability as N → ∞, where the limit process (S,Ī,R) is the unique solution to the system of deterministic Volterra integral equations

The limits (S,Ī) are uniquely determined by the two equations (3.4) and (3.5) . Existence and uniqueness for such a system is well-known. Given the solution (S,Ī), the limitR is given by (3.6) . Assuming that F 0 and F have densities f 0 and f , by taking derivatives in (3.4), (3.5) and (3.6), we obtainS

Note that, as expected,S (t) +Ī (t) +R (t) = 0.

Let us now sketch the proof of Theorem 3.1.

Proof. DefineQ(ds, du) = Q(ds, du) − dsdu the compensated measure. We havē

Since 0 ≤Ῡ N (s) ≤ λ, the first term on the right of (3.7) is an increasing process which is Lipchitz continuous, with a Lipschitz constant bounded by λ. Hence that sequence is equi-continuous, hence is tight in C, hence also in D. The second term is a martingale, which satisfies

hence from Doob's maximal inequality, it converges to 0 in mean square, locally uniformly in t. As a consequence, along a subsequence,Ā N ⇒Ā in D, where for 0 ≤ s < t, 0 ≤Ā(t) −Ā(s) ≤ λ(t − s), and along the same subsequence,

It is not too hard to deduce from the law of large numbers, see Theorem 14.3 in [15] that as N → ∞,

It is not hard to deduce from a slight extension of the Portmanteau theorem (see e.g. Lemma 4.4 in [32] ) that along a subsequence along whichĀ N ⇒Ā, for each t > 0,

From a tightness argument, we deduce that this convergence holds in fact in D. Finally we consider

With some additional effort, one can show that in fact V N (t) → 0 in probability, locally uniformly in t, see the proof of Lemma 5.2 in [64] . Moreover one can show by similar arguments thatR N ⇒R, and that the limiting equations have a unique deterministic solution, hence the whole sequence converges, and the convergence is in probability.

We next turn to the central limit theorem. For that sake, we need to state an appropriate assumption concerning the initial quantities.

In addition, we assume that

We can now state the following result. 

where the limit (Ŝ,Î,R) is the unique solution to the following set of linear stochastic Volterra integral equations driven by Gaussian processes:

(3.8)

withS(t) andĪ(t) given in Theorem 3.1. Here (Î 0 ,R 0 ), independent of (Ŝ(0),Î(0),R(0)), is a mean-zero two-dimensional Gaussian process with the covariance functions: for t, t ≥ 0,

, is a continuous three-dimensional Gaussian process, independent of (Ŝ(0),Î(0),R(0),Î 0 ,R 0 ), and has the representation

where W F is a Gaussian white noise process on R 2 + with mean zero and

for 0 ≤ a ≤ b and 0 ≤ c ≤ d. The limit processŜ has continuous sample paths andÎ andR have càdlàg sample paths. If the c.d.f. F 0 is continuous, thenÎ andR have continuous sample paths. If (Ŝ(0),Î(0),R(0)) is a Gaussian random vector, then (Ŝ,Î,R) is a Gaussian process.

Note that the notion of white noise in defined in Section 8.2 in the Appendix below. From the representation of the limit processes (M A ,Î 1 ,R 1 ) using the white noise W F , we easily obtain their covariance functions: for t, t ≥ 0,

In the FCLT, the limits (Ŝ,Î) are the unique solution of the system of stochastic Volterra integral equations (3.8) and (3.9) . OnceŜ andÎ are specified,R is given by the formula (3.10) .

For the proof of Theorem 3.2, we refer the reader to Section 6 of [64] .

When the infectious periods are deterministic, that is, η i is equal to a positive constant η with probability one, the dynamics of I N (t) can be written as

In the FCLT, we havê

, is a continuous mean-zero Gaussian process with the covariance function

andÎ 1 (t), t ≥ 0, is a continuous mean-zero Gaussian process with the covariance function

Note that the effect of the initial quantities vanish after time η, that is, in the stochastic integral equation (3.12) ofÎ(t), the componentsÎ 0 (t) andÎ(0)(1 − t/η) + vanish after η.

An alternative initial condition. In the model above, we have assumed the remaining infectious periods have a different distribution F 0 . That modeling approach is mostly due to lack of information of infections for these initially infected individuals. An alternative modeling approach is to assume that the infection times of the initially infected individuals are known. We assume the laws of the infectious durations of all individuals are the same, given by the c.d.f. F . This is reasonable since the model is for the same disease. Of course, there is no difference in the two modeling approaches in the Markovian setting due to the lack of memory property of exponential distributions. Suppose that the initially infected individuals are infected at times τ N j,0 , j = 1, . . . , I N (0). Theñ τ N j,0 = −τ N j,0 , j = 1, . . . , I N (0), represent the amount of time that an initially infected individual has been infected by time 0, that is, the age of infection at time 0. WLOG, we can assume that 0 > τ N 1,0 > · · · > τ N I N (0),0 . We use the same notation η 0 j to denote the remaining infectious duration for j = 1, . . . , I N (0). Then the distribution of η 0 j will naturally depend on the elapsed infectious timeτ N j,0 . In particular, the conditional distribution of η 0 j given thatτ N j,0 = s > 0 is given by

Note that the η 0 j 's are independent but not identically distributed. Setτ N 0,0 = 0. Let I N (0, x) = max{j ≥ 0 :τ N j,0 ≤ x}. Assume that there existsx > 0 such that I N (0) = I N (0,x).

We will have the same description of the dynamics of I N (t) in (3.3), however, the variables η 0 j implicitly depend onτ N j,0 . One can explicitly write I N (t) as

Instead of Assumption 3.1, we assume the following holds.

Then it can be shown that the FLLN in Theorem 3.1 holds with the sameS(t) in (3.4) , and

This recovers the result in [82] , specialized to the SIR model. In addition to the FLLN limits, one can also establish the FCLT and obtain the Gaussian limits. Instead of Assumption 3.2, we assume the following holds.

Assumption 3.4. There exist a deterministic functionĪ(0, ·) ∈ (0, 1) and a stochastic procesŝ I(0, ·), and a constant vector (S(0),R(0)) ∈ [0, 1] 2 and a random vector

Then it can be shown that the FCLT in Theorem 3.2 holds with the sameŜ(t) in (3.8), and

withĪ(t) is given in (3.15) , and the limitsÎ 0 (t) andR 0 (t) are as given in Theorem 3.2. However, the limits (Î 0 ,R 0 ) are continuous Gaussian processes with covariance functions: for t, t ≥ 0,

This recovers the result in [84] , specialized to the SIR model.

3.3. The SEIR model. In the SEIR model, the population is split into four groups of individuals: Susceptible, Exposed, Infectious and Recovered/Immune. Let S N (t), E N (t), I N (t) and R N (t) be the numbers of susceptible, exposed and infectious and removed/immune individuals. Infections happen in the same way as the SIR model, that is, through contacts of infectious individuals and susceptible ones according to a Poisson process with rate λ. Thus, the instantaneous infection rate Υ N (t) is given as in (3.1) . Then the cumulative process of newly infected individuals in (0, t], A N (t), has the same expression as in (3.2) . Let L N (t) be the number of individuals that have become infectious after being exposed by time t. We have the following balance equations: for each t ≥ 0,

Each newly infected individual i is associated with the time epoch of being exposed τ N i , exposing duration ξ i and infectious duration η i . Each initially infectious individual j = 1, . . . , I N (0), is associated with the remaining infectious period η 0 j . Each initially exposed individual k = 1, . . . , E N (0), is associated with the remaining exposing time ξ 0 k and the infectious duration η 0 k . Then, we can represent the dynamics of (S N , E N , I N , R N ) as follows: for t ≥ 0,

Assume that (ξ i , η i )'s are i.i.d. bivariate random vectors with a joint distribution H(du, dv), which has marginal c.d.f.'s G and F for ξ i and η i , respectively, and a conditional c.d.f. of η i , F (·|u) given that ξ i = u. Assume that (ξ 0 j , η j )'s are i.i.d. bivariate random vectors with a joint distribution H 0 (du, dv), which has marginal c.d.f.'s G 0 and F for ξ 0 j and η j , respectively, and a conditional c.d.f. of η j , F 0 (·|u) given that ξ 0 j = u. (Note that the pair (ξ 0 j , η j ) is the remaining exposing time and the subsequent infectious period for the i th individual initially being exposed.) In addition, we assume that the sequences

Note that in the case of independent ξ i and η i , e have F (dv) = F (dv|u), and

Similarly, with independent ξ 0 j and η j , we have F 0 (dv) = F 0 (dv|u) = F (dv), and

Define the fluid-scaled processes as in the SIR model. 

The limitS is in C andĒ,Ī andR are in D. If G 0 and F 0 are continuous, then they are in C.

For the proof of Theorem 3.3, as well as for the associated Functional Central Limit Theorem, we refer the reader to [64] . An alternative initial condition. For the initially exposed individuals, let τ N j,0 , j = 1, . . . , E N (0) be the time being exposed, that is,τ N j,0 = −τ N j,0 is the corresponding elapsed duration of being exposed at time zero. For the initially infectious individual k = 1, . . . , I N (0), let ς N k,0 be the time when an initially infectious individual becomes infectious, and thusς N k,0 = −ς N k,0 is the elapsed time at time 0 since being infectious. WLOG, assume that 0 > τ N 1,0 > · · · > τ N E N (0),0 (equivalently, 0 <τ N 1,0 < · · · <τ N E N (0),0 ). Setτ N 0,0 = 0. Then we can write E N (0, x) = max{j ≥ 0 :τ N j,0 ≤ x}. Similarly, assume that 0 > ς N 1,0 > · · · > ς N I N (0),0 (equivalently, 0 >ς N 1,0 > · · · >ς N I N (0),0 ) and set ς N 0,0 = 0. Then we can write I N (0, x) = max{k ≥ 0 :ς N k,0 ≤ x}. We also assume that there exist constants 0 ≤x e < ∞ and 0 ≤x < ∞ such that E N (0) = E N (0,x e ) and I N (0) = I N (0,x) a.s.

For the initially infectious individuals, given their elapsed infection timesς N k,0 , k = 1, . . . , I N (0), we assume that their remaining infectious times are conditional independent and have distributions dependent on their own infection ages, that is, given that theς N k,0 = s,

For the initially exposed individuals, the pairs {(ξ 0 j , η E,0 j )} are assumed to be conditionally independent given the elapsed exposed times {τ N j,0 }, and the distribution of (ξ 0 j , η E,0 j ) depends on the τ N j,0 . In particular, given thatτ N j,0 = s, the joint distribution of (ξ 0 j , η E,0 j ) is given by H 0 (du, dv|s). Assume that the marginal distribution of ξ 0 j given thatτ N j,0 = s is given by

and the conditional distribution of η E,0 j givenτ N j,0 = s and ξ 0 j = u, and given by F 0 (·|u, s) = F 0 (·|s+u). Thus, H 0 (du, dv|s) = G(du|s)F 0 (dv|u, s) = G(du|s)F 0 (dv|u + s).

Let

Similarly, with independent ξ 0 j and η E,0 j given thatτ N j,0 = s, we have

If, in addition, F 0 (·|s) = F (·), then 

An FCLT can be similarly proved, which we omit for brevity.

Multipatch SIR model. Unlike in the case of Markov models, the LLN and CLT for non-Marvovian model in discrete space is not an immediate Corollary of the general result. Multi-patch non Marvovian SEIR models have recently been studied in [65] . We now present the SIR version of that model. The patches may refer to populations in different locations, for example, a densely populated city and a less populated rural area. Individuals in each patch are infected locally and from distance, namely, infectious individuals in each patch may infect any susceptible inside the patch (locally), and susceptible individuals in each patch may get infected from "short" visits to other patches, and are counted as infected in its own patch (from distance). The rate of infection is different in the patches (because of the differences in the density of population or the use of public transportations), while the law of the infectious period is the same (same illness).

Let N be the total population size and L be the number of patches. For each i = 1, . . . , L, let S N i (t), I N i (t) and R N i (t) denote the numbers of individuals in patch i that are susceptible, infectious and recovered at time t, respectively. Then we have the balance equation:

. . , L. Let λ i be the infection rate of patch i, i = 1, . . . , L. The instantaneous infection rate process of patch i is given by,

. . , L , where κ ii = 1 and 0 ≤ κ ij < 1 for i = j represent the infectivity from distance, and 0 ≤ γ ≤ 1. In the case γ = 1, the rate of encounters of individuals in patch i by an infectious (either of the same patch, or from another patch) is fixed, equal to λ i , whatever the total population in patch i at time t may be. In case γ = 0, the same rate is proportional to the total population of patch i at time t. In the intermediate cases, the rate lies between those two extremes. The FLLN is proved for for any value of γ ∈ [0, 1], and the FCLT is only for γ ∈ [0, 1) in the general case, and for all γ ∈ [0, 1] in the case that infections are only local, i.e., κ ij = 0 for i = j. By convention, we shall assume that

, and the same in case γ = 1, i.e., 0 0 = 0 (it is of course 0 if patch i is empty).

Let A N i (t) be the cumulative counting process of individuals in patch i that become infectious during (0, t]. Then we can give a representation of the process A N i (t) via the standard Poisson random measure Q i on R 2 + (with mean measure dsdu):

Equivalently, we can write

where A i, * is a unit-rate Poisson process, and independent from each other for i = 1, . . . , L. We let {τ N j,i , j ≥ 1} denote the successive jump times of the process A N i , for i = 1, . . . , L. Susceptible (resp. infectious, resp. removed) individuals migrate from patch i to patch j at rate ν S,i,j (resp. at rate ν I,i,j , resp. at rate ν R,i,j )). Let X(t) denote the location (i.e., the patch) at time t ≥ 0 of an infected individual. It is clear that X(t) is a Markov process which alternates between states 1, . . . , L. Define p i,j (t) = P(X(t) = j|X(0) = i) for i, j = 1, . . . , L and t ≥ 0. (Note that the probability is the same for any starting time, for example, p i,j (t) = P(X(r + t) = j|X(r) = i) for any r, t ≥ 0.)

We will use X 0,k i and X k i to indicate the associated process for individual k in patch i, for the initially and newly infected ones, respectively. Note that they are all mutually independent and have the same law as described for the process X(t) above. Also note that the processes X k i start from the time becoming infected τ N k,i while the processes X 0,k i start from time 0. We now provide a representation of the epidemic evolution dynamics:

where P S,i, , P R,i, , i, = 1, . . . , L, are all unit-rate Poisson processes, mutually independent, and also independent of P A,i . Here, the first term in I N i (t) represents the number of initially infected individuals from patch = 1, . . . , L that remain infected and are in patch i at time t, and the second term represents the number of newly infected individuals from patch = 1, . . . , L that remain infected and are in patch i at time t. The first term in R N i (t) represents the number of initially infected individuals from patch = 1, . . . , L that have recovered by time t and were in patch i at the time of recovery, and the second term represents the number of newly infected individuals from patch = 1, . . . , L that have recovered by time t, and were in patch i at the time of recovery.

It is not easy to take the limit as N → ∞ in the formula (3.22) . We give another representation of the process I N i (t). Lemma 3.1. We have

where P I,i,j , i, j = 1, . . . , L, are all unit-rate Poisson processes, mutually independent, and also independent of P A,i , P S,i,j and P R,i,j .

Let us comment on this formula. The first term counts the number of initially infectious individuals in patch i, and the second term adds the number of those who get infected in patch i on the time interval (0, t]. The last two terms describe the movements of infectious individuals out of i, and into i. The third terms subtracts the number of individuals initially infected in any patch, who have recovered before time t in patch i, and the fourth term subtracts the number of individuals infected on the time interval (0, t] in any patch, who have recovered before time t in patch i.

Define a PRM Q (ds, du, dv, dθ) on R 3 + × {1, . . . , L}, which is the sum of the Dirac masses at the points (τ N j, , U N j, , η j,i , X j (η j, )) with mean measure ds × du × F (dv) × µ (v, dθ), where for each v > 0, µ (v, { }) = p , (v), and an infection occurs at time τ N j, in case U N j, ≤ Υ N (τ N j, ). We can then write for , = 1, . . . , L,

Remark 3.4. If we set some migration rates to zero, then the corresponding patches could be considered as sub-groups (like age groups) of the population, which interact and infect one another.

. . , L) in probability in R 3L as N → ∞. In addition, assume that F 0 is continuous. 

in probability, locally uniformly on [0, T ], where (S i (t),Ī i (t),R i (t), i = 1, . . . , L) ∈ C 3L is the unique solution to the following set of deterministic integral equations:

For any process Z N , letẐ N := √ N (Z N −Z) be the diffusion-scaled process whereZ N is the fluid-scaled process andZ is its limit. where the limits are the unique solution to the following set of stochastic Volterra integral equations driven by Gaussian processes:

Here, with the notationĪ (i) (t) = L j=1 κ ijĪj (t),

with B A,i , B S,i,j , B I,i,j , B R,i,j being mutually independent standard Brownian motions, and with the deterministic functionsS i ,Ī i ,R i being the limits in Theorem 3.4. The processesÎ 0 i,j andÎ i,j are continuous Gaussian processes with mean zero and covariance functions:

otherwise.

In addition,Î 0 i,j andÎ i,j are independent, and also independent of the Brownian terms. Remark 3.5. The analysis can be easily extended to the multi-patch SIS model, where the population in each patch has susceptible and infectious groups, and when infectious individuals recover, they become susceptible immediately. The epidemic evolution dynamics is described as

, for i = 1, . . . , L. Thus, in the FLLN, we obtain the same limitĪ i in (3.24) as in the multi-patch SIR model, and the limitS i (t):

Similarly in the FCLT, we obtain the same limitÎ i as in (3.25) for the multi-patch SIR model, and the limitŜ i (t):

4. Models with varying infectivity and limiting integro-differential models 4.1. Stochastic model with varying infectivity, LLN and CLT. In this section, we shall consider the same model as in the original work of Kermack and McKendrick [51] , except that we shall formulate a continuous time stochastic individual based model, which as the size N of the population tends to ∞, converges to their model (but our model is slightly more general, since we do not assume that the law of the infectious period is absolutely continuous).

As usual, the population consists of three groups of individuals, susceptible, infected and recovered. Let N be the population size, and S N (t), I N (t), R N (t) denote the sizes of the three groups, respectively. We have the balance equation N = S N (t) + I N (t) + R N (t) for t ≥ 0. Assume that S N (0) > 0, I N (0) > 0 and R N (0) are such that S N (0) + I N (0) + R N (0) = N . Infections occur through interactions of infected individuals with the susceptibles, as in the standard models.

Each initially infected individual is associated with an infectivity process λ 0 j (t), j = 1, . . . , I N (0), which are assumed to be i.i.d. Each newly infected individual is associated with an infectivity process λ i (t), i ∈ N, which are also assumed to be i.i.d. We assume moreover that (S N (0), I N (0), R N (0)), {λ 0 j } j≥1 and {λ i } i≥1 are mutually independent. Assume that λ 0 j (0) = 0 and λ i (0) = 0 with probability one. These processes are only taking effect during the infectious periods. Define 

Thus, the instantaneous infectivity rate function at time t is

Observe that in comparison with the Υ N (t) in (3.1) of the standard model, we have replaced λI N (t) by the total force of infection F N (t) in the generalized model. It is clear that the standard SIR model is the particular case of the present result, where λ(t) = λ1 t<η , η being the random duration of the infectious period. The cumulative infection process A N (t) is expressed exactly as in (3.2) , using the instantaneous infectivity rate function Υ N (t) in (4.2). The epidemic dynamics of the model can be described in the same way as the standard SIR models, that is,

Remark 4.1. The SEIR model. Suppose that λ i (t) = 0 for t ∈ [0, ξ i ), where ξ i < η i , and denote I as the compartment of infected (not necessarily infectious) individuals. An individual who gets infected at time τ N i is first exposed during the time interval [τ N i , τ N i + ξ i ), and then infectious during the time interval

All what follows covers perfectly this situation. In other words, our model acomodates perfectly an exposed period before the infectious period, which is important for many infectious diseases, including the Covid-19. However, we distinguish only three compartments, S for susceptible, I for infected (either exposed or infectious), R for recovered. Note that we could also describe the evolution of the numbers of individuals in the four compartments S, E, I and R as it is done in [32] .

We make the following assumptions on λ 0 and λ.

Assumption 4.1. The random functions λ(t) (resp. λ 0 (t)), of which λ 1 (t), λ 2 (t), . . . (resp. λ 0 1 (t), λ 0 2 (t), . . .) are i.i.d. copies, satisfy the following assumptions. There exists a constant λ * < ∞ such that sup t∈[0,T ] max{λ 0 (t), λ(t)} ≤ λ * almost surely, and in addition there exist a given number k ≥ 1, a random sequence 0 = ξ 0 ≤ ξ 1 ≤ · · · ≤ ξ k = η and random functions λ j ∈ C, 1 ≤ j ≤ k such that

for t ≥ 0. Also, let v 0 (t) = Var(λ 0 (t)) and v(t) = Var(λ(t)) for t ≥ 0. 5) and the limitsĪ andR are given by the following integral equations: Suppose now that F and F 0 have densities f and f 0 . Then we can rewrite the above deterministic model as follows. For the convenience of the comparison with the model in [51] , we replace the set of variables (S(t),F(t),Ī(t),R(t)) by the set of variables (S(t),Ῡ(t),Ī(t),R(t)), whereῩ(t) =S(t)F(t).

If we assume that F 0 ≡ F , then the last system of equations is exactly the system of equations (12) , (13) , (15) and (14) on page 704 of [51] . Indeed, it follows from the computation at the start of Section 3.1 that the function B t of [51] is our F c (t), while their C t is our f (t). Moreover their A t is ourλ(t) (indeed, one can think of our λ(t) as being the product of a deterministic function of t (their φ t ) multiplied by 1 η>t , so thatλ(t) = φ t F c (t).

We now sketch the proof of Theorem 4.1.

Proof. Thanks to Assumption 4.1, it is clear thatῩ N (t) ≤ λ * . Hence the first step of the proof of Theorem 3.1 remains valid here, i.e., we have the convergence, along a subsequence, of (S N ,Ā N ). We next consider the sequencē

Concerning the first termF N 0 , as in the proof of Theorem 3.1, we first consideȓ

which converges thanks to a LLN for random elements in D, see Theorem 1 in [72] . The differencē

is treated as in the proof of Theorem 3.1. Concerning the termF N 1 , we first consideȓ

The argument for the weak convergence of that sequence towardsF(t) = t 0λ (t − s)dĀ(s), along any subsequence along whichĀ N ⇒Ā is similar to a similar result in the proof of Theorem 3.1, with slightly more tricky arguments. For the details, as well as for the proof of the fact that F N 1 (t) −F N 1 (t) → 0, we refer to Section 4 of [32] . It remains to prove that (Ī N (t),R N (t)) ⇒ (Ī(t),R(t)), which requires similar arguments as in the first steps of the proof. Finally, one can show that the limiting equation has a unique deterministic solution, hence the whole sequence converges, and the convergence is in probability.

For the FCLT, we need the following additional conditions on the random infectivity functions.

In addition to the conditions in Assumption 4.1, the random functions λ(t) (resp. λ 0 (t)) satisfy the following conditions.

(i) There exist nondecreasing functions φ and ψ in C and α > 1/2 and β > 1 such that for all

(ii) Either λ ∈ C and satisfies (4.6)-(4.7) below, or else it satisfies (4.3) and the additional conditions below. There exists a nondecreasing function ϕ ∈ C satisfying ϕ(r) ≤ Cr α , with α > 1/2 and C > 0 arbitrary,

such that |λ j (t) − λ j (s)| ≤ ϕ(|t − s|), a.s., (4.7) for all t, s ≥ 0, 1 ≤ j ≤ k. Also, if F j denotes the c.d.f. of the r.v. ξ j , then the exist C and ρ > 0 such that for any 0 ≤ j ≤ k, 0 ≤ s < t,

and in addition, for any 1 ≤ j ≤ k, r > 0, The limit process (Ŝ,F) is the unique solution of the following system of stochastic integral equations:

(s)ds, (4.8)

andS(t) andF(t) are given by the unique solutions to the integral equations (4.4) and (4.5),M A , F 0 ,F 1 andF 2 are centered Gaussian processes which are globally independent of (Ê(0),Î(0)). Moreover, the processes in (F 0 ,F 1 , (F 2 ,M A )) are independent, and the covariances of each of those four processes (the last one being 2-dimensional) are given as follows :

Concerning the pair (M A ,F 2 ),M A is a non-standard Brownian motion, andF 2 (t) = t 0λ (t − s)M A (ds).Ŝ has continuous paths, and ifλ 0 andλ 0,I are in C, thenF is also continuous.

The limits (Î,R) are given bŷ

where (Î 0 ,R 0 ), independent ofÎ(0), is a mean-zero two-dimensional Gaussian process with the covariance functions as given in (3.11) , and the limits (Î 1 ,R 1 ), independent of (F 0 ,Î 0 ,R 0 ) andÎ(0), are a continuous two-dimensional Gaussian process with mean zero and covariance functions, for t, t ≥ 0,

In addition, for t, t ≥ 0,

An alternative initial condition. In the above formulation, we have assumed that for the initially infected individuals, their infectivity functions λ 0 j (·) are i.i.d., and may follow a different law from those of the newly infected individuals λ i (·), hence, the distribution F 0 of the remaining infected periods η 0 j generated from λ 0 j (·), is different from F of the infected periods η i generated from λ i (·). However, we can assume that the random infectivity functions of all individuals, {λ 0 j (·)} j and {λ i (·)} i are all i.i.d., while for the initially infected individuals, the time epochs of them becoming infected before time 0 are known, τ N j,0 , j = 1, . . . , I N (0). Thenτ N j,0 = −τ N j,0 , is the elapsed time at time 0 since infection. Setτ N 0,0 = 0, and let I N (0, x) = max{j ≥ 0 :τ N j,0 ≤ x}. Assume that there existsx ∈ R + , such that I N (0) = I N (0,x).

for t ≥ 0. The remaining infected period is given by η 0 j = inf{t > 0 : λ 0 j (τ N j,0 + r) = 0, ∀r ≥ t}. It depends on the elapsed infection timeτ N j,0 , and independent from the remaining infected durations of the other individuals due to the i.i.d. assumption of {λ 0 j (·)} j . Given thatτ N j,0 = s > 0, the distribution of η 0 j is given as in (3.13) , that is,

Instead of (4.1), the total force of infectivity at time t can be written as

All the other processes have the same representations. Recall that the process I N (t) is given as in (3.14) with the variables η 0 j implicitly depending onτ N j,0 . Under Assumption 3.3, we can show that the FLLN holds withS(t) in (4.4) , and the limitF(t) is given asF

(y + t)Ī(0, dy) + t 0λ (t − s)S(s)F(s)ds , (4.13) and the limitsĪ andR are given bȳ

Under Assumption 3.4, we can show that the FCLT holds with the limitŜ(t) in (4.8), the limit G(t) given byF

whereΥ(t) is given in (4.10), andS(t) andF(t) are given by the unique solutions to the integral equations (4.4) and (4.13),F 0 is a continuous Gaussian process with mean zero and covariance function: for t, t ≥ 0,

Cov(λ(y + t), λ(y + t ))Ī(0, dy), and the other limitsM A ,F 1 andF 2 are centered Gaussian processes as given in Theorem 4.2.

The limits (Î,R) are given bŷ

where (Î 0 ,R 0 ) are continuous Gaussian processes with covariance functions given in (3.17)- (3.19) .

In addition,F 0 andÎ 0 (t),R 0 (t) have covariance functions: for t, t ≥ 0,

The early phase of the epidemic. In this subsection we follow again [32] , to which we refer the reader for the proofs. Theorem 4.1 shows that the deterministic system of equations (4.4)-(5.7) accurately describes the evolution of the stochastic process defined in the previous subsection when the initial number of infectious individuals is of the order of N . But epidemics typically start with only a handful of infectious individuals, and it takes some time before the epidemic enters the regime of Theorem 4.1. Exactly how long this takes depends on the population size N and on the growth rate of the epidemic. To determine this growth rate, we study the behavior of the stochastic process when the initial number of infectious individuals is kept fixed as N → ∞.

Recall that R 0 = ∞ 0λ (t)dt, and let ρ ∈ R be the unique solution of ∞ 0 λ(t)e −ρt dt = 1. (4.14)

If R 0 ≤ 1, the total number of infected individuals remains small as N → ∞, while if R 0 > 1 (which we assume in what follows), with positive probability a major outbreak takes place, i.e., a positive fraction of the N individuals is infected at some point during the course of the epidemic. It is well-known, see e.g. section 1.2 in [23] , that during its early stage, an epidemic can be well approximated by a continuous time branching process. Indeed, each infectious infects individuals in the population, and as long as almost all the individuals in the population are susceptible, the probability that two distinct infectious individuals try to infect the same susceptible is close to 0. As a result the "progenies" of the various infectious individuals are essentially independent, thus the branching property. Of course, that approximation breaks down as soon as a significant number of individuals have been hit by the disease. Using an approximation of the early phase by a (in our case non-Markov) branching process, it has been shown in [32] that on the event that a major outbreak takes place, for any ε < 1 − R −1 0 , if T N ε denotes the first time at which the proportion of infected individuals is at least ε, as N → ∞, T N ε = 1 ρ log(N ) + O(1), which means an exponential growth with rate ρ.

Next one can show that, still at the start of the epidemic, our LLN deterministic model also grows at the same rate ρ. More precisely, if we assume that we can replaceS(t) by 1, the LLN model becomes (after remultiplication by N ) the following linear system:

In the next statement, ρ is specified by (4.14). Suppose first that R 0 > 1, hence ρ > 0. Then, if λ 0 = λ ρ and F 0 = F ρ , the above linear system admits the following solution F(t) = ρ e ρt , I(t) = i e ρt , R(t) = r e ρt t ≥ 0.

If, however, R 0 < 1, hence ρ < 0, then the above linear system (with λ 0 = λ ρ and F 0 = F ρ ) admits the following solution

Let us suppose that λ is only known up to a constant factor µ > 0, i.e.,

where µ is unknown but g is known (for example from medical data on viral shedding). We now assume w.l.o.g. that g has been normalized in such a way that ∞ 0 g(t)dt = 1. We can then estimate µ (and R 0 ) from the growth rate ρ, which can be measured easily at the beginning of the epidemic (ρ = log(2)/d, where d is the doubling time of the daily number of newly infected individuals), using the relation (4.14). The following is thus a corollary of Theorem 4.3. (4.15)

Note that R 0 can be thought of as the growth rate of the epidemic from one generation to the next, while ρ is the growth rate of the epidemic in real time. The formula (4.15) is formula (2.7) in [81] .

Application to the Covid-19 epidemic in France. We now explain how the type of model described in this section can be used to model the Covid-19 epidemic. As we have seen, the increase in realism with respect to the classical "Markovian" models (where the infectivity is constant and fixed across the population, and the Exposed and Infectious periods follow an exponential distribution) is paid by replacing a system of ODEs by a system of Volterra integral equations. However, we have a small benefit in that the flexibility induced by the fact that the law of λ is arbitrary allows us to reduce the number of compartments in the model, so that we can replace a system of ODEs by a system of Volterra type equations of smaller dimension. Figure 1 . Flow chart of the SEIRU model of [57] and of our SIR model. We are able to replace the six compartments of the SEIRU model with only three compartments by using the equations described in Theorem 4.1.

To be more specific, let us describe the SEIRU model of [57] . An individual who is infected is first "Exposed" E, then "Infectious" I. Soon after, the infectious individual either develops significant symptoms, and then will be soon "Reported" R, and isolated so that he/she does not infect any more; while the alternative is that this infectious individual is asymptomatic: he/she develops no or very mild symptoms, so remains "Unreported" U, and continues to infect susceptible individuals for a longer period. Both unreported and reported cases eventually enter the "Removed" (Rem) compartment. In this model, there are 6 compartments: S like susceptible, E like exposed, I like infectious, R like reported, U like unreported, and Rem like removed.

Our approach allows us to have a more realistic version of this model with only 3 compartments (see Figure 1 ): S like susceptible, I like infected (first exposed, then infectious), R like removed (which includes the Reported individuals, since they do not infect any more, and will recover soon or later). As already explained, we do not need to distinguish between the exposed and infectious, since the function λ is allowed to remain equal to zero during a certain time interval starting from the time of infection. More importantly, since the law of λ is allowed to be bimodal, we can accommodate in the same compartment I individuals who remain infectious for a short duration of time, and others who will remain infectious much longer (but probably with a lower infectivity). Moreover, since we know, see [42] , that the infectivity decreases after a maximum which in the case of symptomatic individuals, seems to take place shortly before symptom onset, our varying infectivity model allows us to use a model corresponding to what the medical science tells us about this illness. Note that our version of the SEIRU model from [57] is the same as the one which we have already used in [33] (except that there we had to distinguish the E and the I compartments). However, the main novelty here is that the infectivity decreases after a maximum near the beginning of the infectious period. More precisely, we consider that t → g(t) increases linearly on the time interval [ζ, ζ + η/5], from 0 to 1 for reported individuals, and from 0 to α for unreported individuals, and that it then decreases linearly to 0 on the interval [ζ + η/5, ζ + η], as shown on Figure 2 . We then take (X 1 , X 2 ) a pair of independent Beta random variables with parameters (2, 2) and we assume that ζ = 2 + 2X 1 , η = 3 + X 2 for reported individuals, 8 + 4X 2 for unreported individuals.

This joint law of (ζ, η) is the one that was used in [33] to study the Covid-19 epidemic in France (where the infectivity was assumed to be constant and uniform among individuals in this work), and these values are compatible with the results described in [42] .

Kermack and McKendrick pioneered the introduction of PDE models to describe infection age dependent infectivity and recovery age dependent susceptibility in their 1932 paper [52] . Here we shall describe the PDE model for an infection age dependent infectivity, in the framework of a SIR/SEIR, i.e., where the recovered individuals do not loose their immunity. As in the rest of this paper, we shall obtain the deterministic (here PDE/integral equation) model as a LLN limit of individual based stochastic models. The underlying assumptions are the same as in the previous section, except for a slight modification concerning the initially infected individuals. But we shall give a different description of the model, as we shall see now.

We have the same compartments as in Section 4.1, S N (t), I N (t) and R N (t) are as above, and again S N (t) + I N (t) + R N (t) ≡ N . Let now I N (t, x) be the number of infected individuals at time t that have been infected for a duration less than or equal to x. Note that for each t, I N (t, x) is nondecreasing in x, which is the distribution of I N (t) over the infection-ages. Let A N (t) be the cumulative number of newly infected individuals in (0, t], with the infection times {τ N i : i ∈ N}. Each individual who has been infected after time 0 has an infectivity process λ i (·), and we assume that these random functions are i.i.d.. Let η i = inf{t > 0 : λ i (r) = 0, ∀r ≥ t} be the infected period corresponding to the individual that gets infected at time τ N i . The η i 's are i.i.d., with a cumulative distribution function (c.d.f.) F . Let F c = 1 − F .

Let {τ N j,0 , j = 1, . . . , I N (0)} be the times at which the initially infected individuals at time 0 were infected. Thenτ N j,0 = −τ N j,0 , j = 1, . . . , I N (0), represents the age of infection of individual j at time 0. W.l.o.g., we assume that 0 > τ N 1,0 > τ N 2,0 > · · · > τ N I N (0),0 (or equivalently 0 <τ N 1,0 <τ N 2,0 < · · · < τ N I N (0),0 ). Setτ N 0,0 = 0. We define I N (0, x) = max{j ≥ 0 :τ N j,0 ≤ x}, the number of initially infected individuals that have been infected for a duration less than or equal to x at time 0. Assume that there exists 0 ≤x < ∞ such that I N (0) = I N (0,x) a.s.

To each initially infected individual j = 1, . . . , I N (0), is associated an infectivity process λ 0 j (·), and we assume that they are also i.i.d., with the same law as λ i (·). For each j, let η 0 j = inf{t > 0 : λ 0 j (τ N j,0 + r) = 0, ∀r ≥ t} be the remaining infectious period, which depends on the elapsed infection timeτ N j,0 , but is independent of the elapsed infection times of other initially infected individuals. In particular, the conditional distribution of η 0 j given thatτ N j,0 = s > 0 is given by

Note that the η 0 j 's are independent but not identically distributed. For an initially infected individual j = 1, . . . , I N (0), the infection age is given byτ N j,0 + t for 0 ≤ t ≤ η 0 j , during the remaining infectious period. For a newly infected individual i, the infection age is given by t − τ N i , for τ N i ≤ t ≤ τ N i + η i during the infectious period. Note that λ i (·) and λ 0 j (·) are equal to zero on R − .

The aggregate force of infection at time t is given by

We have again (4.2) and (3.2) . Moreover the total number of individuals infected at time t that have been infected for a duration which is less than or equal to x:

x) is the number of initially infected individuals who have been infected for a duration less than or equal to x at time t, which is given as

and I N 1 (t, x) is the number of newly infected individuals who have been infected for a duration less than or equal to x at time t, which equals

Note that for each t, I N 0 (t, ·) has support over [0, t +x] and I N 1 (t, ·) has support over [0, t]. Thus I N (t) = I N 0 (t, t +x) + I N 1 (t, t) = I N (t, ∞), t ≥ 0. The sample paths of I N (t, x) belong to the space D D , denoting D(R + ; D(R + ; R)), the D-valued D space.

Define the fluid-scaled processesX N = N −1 X N for any processes X N . We make the following assumptions on the initial quantities. 

The functionĪ(t, x) is nondecreasing in x for each t, the integrals w.r.t. d xĪ (0, y) and d xĪ (t, x) are Lebesgue-Stieltjes integrals with respect to the measure which coincides with the distributional derivative ∂ xĪ (0, ·) =Ī x (0, ·) (resp. ∂ xĪ (t, ·) =Ī x (t, ·)). As a consequence,Ī N →Ī in D in probability as N → ∞ wherē

The proof of Theorem 5.1 is similar to the proofs of the FLLNs in the previous sections, with the additional complications that we have one function of two parameters. We refer the reader to Section 4 of [67] for that proof.

We now turn to deriving a Partial Differential Equation for the derivative with respect to x of I(t, x), when it exists. In the next result, we shall denote by γ(x) = f (x)/F c (x) the hazard function of the r.v. η, where f (x) = F (x). Proposition 5.1. Suppose that F is absolutely continuous, with the density f , and thatĪ(0, x) is differentiable with respect to x, with the density functionī(0, x). Then for t > 0, the increasing functionĪ(t, ·) is absolutely continuous, andī(t, x) := ∂ xĪ (t, x) satisfies (t, x) a.e. in (0, +∞) 2 ,

with the boundary conditionsī(0, x) =Ī x (0, x) for x ∈ [0,x], and {ī(t, 0) , t ≥ 0} is the unique non-negative solution of the following Volterra equation

with the initial conditionī(0, 0) = ∂ xĪ (0, x)| x=0 . In addition,

Moreover, the PDE (5.10) has a unique solution which is given as follows. For x ≥ t,

Remark 5.2. We remark that the PDE in [52] resembles (5.10), see equations (28)- (29) , see also equation (2.2) in [46] . In particular, the function γ(x) is interpreted as the recovery rate at infection age x. Equivalently, it is the hazard function of the infectious duration.

Proof. By the fact that F has a density, we see that the two partial derivatives ofĪ exist (t, x) a.e. From (5.6), the sum of the two partial derivatives satisfies for t > 0 and x > 0,

By taking the derivative on both sides of the last identity with respect to x (possibly in the distributional sense for each term on the left), we obtain the expression Thus we obtain the expression (5.11). We next prove that equation (5.11) has a unique non-negative solution.

Observe that x(t) =ī(t, 0) is also a solution to

and any non-negative solution of (5.11) solves (5.16). First, note that for any t ≥ 0, 0 ≤ x 0λ (y + t)ī(0, y)dy ≤ λ * Ī (0) and 0 ≤λ(t) ≤ λ * , from which we conclude that t 0ī (s, 0)ds ≤S(0). Indeed, if that were not the case, there would exist a time TS (0) < t such that

0ī (s, 0)ds =S(0), hence t 0ī (s, 0)ds ≥S(0) and from (5.16), we would have x(t) = 0 for any t ≥ TS (0) .

With those bunds, it is not hard to prove that the Volterra equation (5.16) has a unique non negative solution.

We next derive (5.13) and (5.14). We first deduce from (5.15 ) that for x ≥ 0, where

for t, x ≥ 0. The density functionī(t, x) = ∂Ī(t,x) ∂x , if it exists, satisfies again (5.10). The same calculations as in the case of the SIR model lead to (5.13), (5.14) and

However, the formula forS(t) is different in the case of the SIS model. We havē

Thus, the Volterra equation on the boundary reads

whose form is similar to the one for the SIR model. Recall that the standard SIS model has a nontrivial equilibrium pointĪ * = 1 − µ/λ if µ < λ, where λ is the infection rate (the bar over λ is dropped for convenience), and 1/µ is the mean of the infectious periods. See Section 4.3 in [64] for the account of the SIS model with general infectious periods. Here we consider the model in the generality of infection-age dependent infectivity.

Proposition 5.2. If R 0 = ∞ 0λ (y)dy ≤ 1, the only equilibrium isĪ * = 0 (the disease free equilibrium). In the complementary case, R 0 = ∞ 0λ (y)dy > 1, there is a unique endemic equilibrium with a proportion of infectious individuals equal tō

The density functionī(t, x) has an equilibriumī * (x) in the age of infection x, given bȳ

dt is the expectation of the duration of the infectious period. If F has a density f , then the equilibrium densityī * (x) satisfies

Proof. Assume that the equilibriumĪ * (x) :=Ī(∞, x) exists. Then it must satisfȳ

where F e (x) = µ

x 0 F c (s)ds, the equilibrium (stationary excess) distribution. Letting x → ∞ in this formula, we deduceĪ

Combining the last two equations, we obtain

Plugging this formula in the previous identity, we deduce that We believe that the SIR models considered in Sections 4 and 5 can be made, by a proper choice of the parameters, much more realistic that the Markov / ODE models of Section 2. One may however ask the question whether those "refined" models make a real difference, as compared to Markov / ODE models. It is not clear that the large time behaviors of the two kind of models are significantly different. Indeed, we note that the formula for the endemic equilibrium in the SIS model (5.17) in Proposition 5.2 is exactly the same, when expressed in terms of the basic reproduction number R 0 , as the formula obtained in Section 2.3.1.

On the other hand, the transitory behavior can be drastically different in the two types of models. Indeed, we can implement in our "refined" model the memory of recent situations of the epidemic, so that when the rate of contact between individuals changes drastically (e.g. when a government enforces a lockdown), in the "refined" model the number of daily infections will take more time an ODE model obtained as the limit of a Markov stochastic SEIR model and a model with memory obtained as a limit of a non-Markov SEIR model. Both models have the same mean exposed and infectious period, and are chosen so that they have the same initial growth rate. After 28 days, the contact rate is reduced instantly, in such a way that both models should have the same rate of decay of newly infected individuals. We see that the epidemic in the model with memory "slows down" less rapidly than in the ODE model, due to its greater inertia, which causes a larger number of individuals to be infected.

to go down than in the ODE model, see Figure 3 . Some authors who use ODE models correct that behavior by making the contact rate change gradually after the time of lockdown, which does not correspond to the behavior of the population. As observed several times and in various places during the Covid-19 pandemic, daily infections, as well as hospital admissions and hospital deaths, continue to grow for a relatively long time (up to several weeks) after strict preventive measures are taken, a pattern that arises naturally in models with memory but is much harder to reproduce using ODE models. This behavior has consequences both for inference methods and for decision making, since ODE models can underestimate the inertia of the epidemic on the short term.

As in the previous sections, we start with a population of fixed size N , and we enumerate the individuals in the population with the parameter k, 1 ≤ k ≤ N .

The model in this section is a sort of SIRS model, except that we do not really distinguish between the states R and S. We shall only consider the compartments S and I, and formally our model is a SIS model, although individuals experience immunity after recovery, before being susceptible again. Each individual who is infected first draws a random infectivity function, as in the previous section, and may infect other individuals as before. At the end of the infectious period, the individual is first immune (i.e., its susceptibility is equal to zero), but this acquired immunity then wanes with time, and we assume that the susceptibility of the individual increases gradually according to some random function. When such a partially susceptible individual is the target of an infectious contact, the probability that this individual becomes reinfected is given by its susceptibility (so this quantity evolves between 0 and 1). Whenever an individual in the compartment S has susceptibility 0, he/she is in fact immune.

Since each individual might get infected an arbitrary number of times, and the infectivity and susceptibility functions of the infection age are a priori different after each new infection, we attach to each individual a countable family of infectivity and susceptibility functions. More precisely, we consider an independent family of elements of D 2 , (λ k,i , γ k,i ) 1≤k≤N, i≥0 , which are such that they are i.i.d. for 1 ≤ k ≤ N and i ≥ 1, but may have a different law with the (λ k,0 , γ k,0 ) 1≤k≤N . Those last quantities represent the infectivity and susceptibility starting at time t = 0, while for i ≥ 1, (λ k,i , γ k,i ) represents the infectivity and susceptibility of the k-th individual after his/her i-th infection (not counting a possible infection before time 0). Note that all λ k,i take values in [0, λ * ] and all γ k,i take values in [0, 1].

At time 0, individual k can be susceptible (or "naive"). In that case, λ k,0 (t) ≡ 0 and γ k,0 (t) ≡ 1. A second possibility is that individual k is infected. In that case λ k,0 ≥ 0 and γ k,0 (0) = 0. A third possibility is that individual k has recovered at time 0 from a past infection. In that case, λ k,0 (t) ≡ 0 and the function γ k,0 is arbitrary.

In addition to what has been explained above, we assume that

We introduce the following notations:

Let us now describe our individual based stochastic model. Contrary to what we did in the previous sections, we do not just count the number of infections in the population on the time interval (0, t]. We shall denote by A N k (t) the number of times that the individual k has been infected on the time interval (0, t]. Let σ N k (t) denote the age of infection of the individual k at time t, i.e., σ N k (t) := t − sup{s ∈ [0, t], A N k (s) = A N k (s − ) + 1} ∨ 0 , where as usual the sup is 0 if the concerned set is empty. At time t, the infectivity of the individual k is λ k,A N k (t) (σ N k (t)), and its susceptibility is γ k,A N k (t) (σ N k (t)). Note that in the case where A N k (t) = 0, we recover the above description of the situation prior to the first (re)infection.

The total force of infection in the population at time t is

According to the above description, we expect that the rate at which the individual k gets infected is

Let now {Q k , 1 ≤ k ≤ N } be a collection of mutually independent standard Poisson random measures on R 2 + . We assume that the number of infections endured by the individual k on the interval (0, t] is given by

. We finally define the average susceptibility in the population as

We will show that the pair (S N (t),F N (t)) converges to a deterministic pair (S(t),F(t)), locally uniformly in time. Before we state this convergence result, let us study the limiting equation. Let (x, y) ∈ D 2 be a solution to the following set of equations:

(t − s)x(s)y(s)ds . Proof. We first prove uniqueness. We need an a priori bound on the solutions. Suppose (x, y) is a non negative solution of (6.1), i.e., a solution satisfying x(t) ≥ 0, y(t) ≥ 0 for all t ≥ 0. Since γ 0 (t) ≤ 1 and γ(t) ≤ 1, we deduce from the first equation that

If we multiply the first equation in (6.1) by y(t), we obtain an identity which shows that the derivative with respect to t of the right hand side of the above inequality is zero, hence that upper bound equals its value at time t = 0, which is 1. We have proved that x(t) ≤ 1. Next from the second equation and Gronwall's Lemma we deduce that y(t) ≤ exp(λ * t). With the help of those bounds, it is not very hard to show that (6.1) has at most one negative solution. Existence can be shown using a Picard iteration procedure, thanks to the estimates which are used for uniqueness. Note that the solution starts with x(0) > 0 and y(0) > 0, and it is not hard to see that neither x nor y can hit 0 in finite time.

We can now state the main result of this section. In other words, the pair (S,F) solves the system of integral equations

(t − s)S(s)F(s)ds . In that case, the first equation in (6.3) reduces tō

which is the solution of (4.4). So our new result is consistent with Theorem 4.1.

We will sketch the proof of this Theorem, refering the reader to [34] for the details. The main idea of this proof is to replace the collection {A N k , 1 ≤ k ≤ N } by an i.i.d. sequence of random processes, which will be close to that collection in an appropriate sense. The idea of the construction of that sequence is the following. A N k depends upon N only throughF N , which is a mean field interaction. Suppose we try to replaceF N (t) by a deterministic function m(t), in such a way that m(t) = E λ A (m) k (t) (σ (m) (t)) . Then we would already have constructed the wished limit. Let Q be a standard Poisson Random Measure on R 2 + , and (λ i , γ i ) i≥1 an i.i.d. sequence, each one having the law of (λ 1,1 , γ 1,1 ), and which is globally independent of Q. Also take (λ 0 , γ 0 ) independent of the previous sequence and distributed as (λ 1,0 , γ 1,0 ). To each deterministic m ∈ D, we associate the solution A (m) (t) of the following SDE:

where

The next Lemma is crucial for our proof. 

Let (x, y) denote the unique solution of (6.1). Let us choose m = y. Then comparing the last identity to the first equation in (6.1), we deduce that Θ (y) = x, and comparing the previous identity to the second equation of (6.1), we deduce that Ψ (y) = y. Conversely, if Ψ (m) = m, we have that (Θ (m) , m) solves (6.1), hence the result.

We now exploit again the sequence of independent PRMs {Q k , k ≥ 1}, and for each k ≥ 1, we let A k be the A (m * ) associated to Q k . More explicitly, we define for each k ≥ 1,

Note that it follows from Lemma 6.1 that for each k ≥ 1,

The next step in the proof is the following Lemma.

Lemma 6.2. For any k ≥ 1 and T > 0, we have

Proof. The first inequality is rather obvious. We now establish the second inequality. We will use repeatedly the fact that the r.v. sup 0≤r≤t |A N k (r) − A k (r)| is either 0 or else ≥ 1. First note that

The first term on the right is bounded by

By a standard computation, the second term on the right hand side is bounded by λ * / √ N , and the first term by λ * P sup 0≤r≤t |A N k (r) − A k (r)| . If we define δ N (t) := P sup 0≤r≤t |A N k (r) − A k (r)| , combining the above computations yields

The result now follows from Gronwall's Lemma.

Completing the proof of Theorem 6.1. The remainder of the proof of Theorem 6.1 can be sketched as follows. The r.v. whose expectation is close to 0 by Lemma 6.2 is 0 with probability close to 1 for large N . It is then not hard to deduce that both

tend to 0 in probability. Since moreover the collection (A N k , σ N k ) 1≤k≤N is exchangeable, and the sequence (A k , σ k ) k≥1 is i.i.d., the result essentially follows from Proposition I.2.2 in [76] . Remark 6.2. While the random infectivity appears in the limiting LLN deterministic equations only through its mean functionλ(t), a complicated mixed moment-exponential moment of the trajectory of the random susceptibility γ(t) appears in the deterministic version of our varying infectivity / varying susceptibility model. This is due to the possibility of reinfection of the individuals who are experiencing a graduate loss of their immunity / gain of their susceptibility.

Using the same techniques, we can also obtain the limiting equations for the proportion of susceptible and infectious individuals. As before, we let

Then define

i.e., the number of infectious individuals at time t (recall that η k,0 = 0 if the k-th individual is initially susceptible). Also set

Then, settingĪ N (t) = 1 N I N (t) andS N (t) = 1 N S N (t), we have the following convergence. where (η 0 , γ 0 ) is distributed as (η 1,0 , γ 1,0 ) and (η, γ) as (η 1,1 , γ 1,1 ).

Note that we can recover the fact thatĪ(t) +S(t) = 1 for all t ≥ 0 from the fact that γ(s) = 0 for all s < η (resp. γ 0 (s) = 0 for all s < η 0 ) and the fact that the right hand side of (6.2) equals 1.

Proof of Corollary 6.1. This convergence follows from Lemma 6.2 in much the same way as the convergence ofS andF. Lemma 6.2 implies that, for any fixed k,

tends to zero in probability as N → ∞, and the convergence ofĪ N follows again from the exchangeability of the (A N k , σ N k ) 1≤k≤N and Proposition I.2.2 in [76] . The convergence ofS N then follows from the fact thatS N (t) = 1 −Ī N (t) andS(t) = 1 −Ī(t).

Let us now show how our equations take a more explicit form in a particular case, which allows a comparison with the equations in Kermack and McKendrick [52] . Suppose that for i ≥ 1, our function γ(t) = g(t − η), where g is deterministic and equals 0 on R − , and η is the duration of the infectious period. We assume that he r.v. η has the density f (t).

We now need to specify γ 0 (t). Suppose we have a random variable χ, which is independent of all other random sources, and takes the value S with probabilityS(0), the value I with probabilityĪ(0) and the value R with probabilityR(0), whereS(0) +Ī(0) +R(0) = 1. We assume that

where η 0 has the density f 0 (t) and ξ has the density h(t).

Let us now rewrite the first equation of (6.3) in the just specified particular case. We get: (6.6) Remark 6.3. It can be shown that (6.6) together with the second equation in (6.3) are equivalent to the system of non local PDEs from [52] , see [34] for the details of the calculations.

We thus see that (6.3) yields a more general deterministic model than the one proposed in [52] , in the case where the susceptibility is a random function of the time elapsed since the individual has recovered.

7.1. Non homogeneous models. All the models presented so far are homogeneous, in the sense that whenever one infectious individual meets someone else, anyone in the population has the same chance to be met, and to be possibly infected if he/she was susceptible (with the exception of the previous section, where the susceptibility of the various individuals plays a role in that choice, but it does not contradict the homogeneity). There are many reasons why this is not realistic, and we shall indicate several attempts to correct the homogeneous model, and make the epidemic models more realistic. However, the reader should realize that too complicated models may not be really useful, among other reasons because they involve too many parameters, which might not be easy to estimate.

A first complexification of the above models is to distribute the population into age groups. This is quite reasonable concerning the Covid -19 epidemic, since the proportion of severe cases and deaths among those who catch the disease depends very much upon the age. This is not too difficult to implement, provided one can exploit informations about the contact rates between those age groups, which might be found in the sociology literature. Also, in several countries data concerning the numbers of hospitalized patients, those in intensive care units, and those who die, are available by age class. Note however that the health condition of the patients (and their weight) is almost as important as their age.

It has been recently pointed out, see [22] , that the heterogeneity in levels of social activity between individuals has an impact on the herd immunity. Indeed, those who have a higher rate of social contacts have more chances to get infected towards the beginning of the epidemic, and they will infect more people than those with a lower rate of social contacts. Once most of the socially very active individuals have been infected and become immune, one may think that the progress of the epidemic might slow down. For that reason, an epidemic which infects a certain percentage of the population will be more efficient towards building herd immunity than vaccinating the same proportion of individuals. This raises also questions concerning a vaccination campaign. Should one vaccinate first those more at risk, but who have few social contacts and do not contribute much to the propagation of the epidemic, or rather those who have much social contacts?

The reason why reality is not homogeneous has to do with the fact that each individual has frequent contact with those in his close environment (those who share the same household and workplace), and much less frequent contact with people who are met in public transportation, shops, various social activities. There has been a lot of effort to adapt the epidemic models to such situation with various levels of contact rates. [11] contains a recent review of the works in that direction. See also [35] for a mean field model approach to household epidemic models. 7.2. Spatial models. A good reason for non homogeneity is the spatial dispersion of the population. There is of course a strong motivation for studying epidemics models for population distributed in discrete or continuous space. There is a quite significant body of literature on deterministic models in those two situation, see in particular [3] , [4] , [73] and reviews in Chapter 15 of [60] and Chapter 14 of [21] . Recently, some authors have proved law of large numbers and central limit theorems for Markov epidemic models in continuous space, see in particular [63] and [19] . Concerning non-Markov models, we have recently studied such a model in discrete space, see Section 3.4.

There is also an extensive literature on epidemic models on random graphs, including various limit theorems and asymptotic results with large population, large graphs, and various network topologies, see, e.g., [13, 6, 26, 12, 48, 61, 49, 37, 80] , and also the recent survey [78] . As far as we know, most of these works are for Markovian models. It would be interesting to investigate non-Markovian epidemic transmissions on random graphs. 7.3. Control problems. Optimal control problems in the Markovian and limiting ODE epidemic models have been studied extensively in the literature. Isolation, vaccination and immunization strategies have been developed to minimize the epidemic size and costs associated with the implementation of them. For example, optimal isolation strategy to minimize the total number of infected individuals [1] or to minimize the total infectious burden over an outbreak together with a cost for implementing the control in [87, 62] , optimal vaccination strategy [2, 87, 62] , optimal immunization strategy [86] , optimal combined isolation-vaccination strategy with resource constraints [41] , and optimal control to minimize the total number of infectious and the time needed for the infection to go extinct [18, 17] . To cope with Covid-19 pandemic, various lockdown, social distancing, testing and vaccination strategies have been implemented by governments. Some studies have been conducted for their effect, see, e.g., [28, 29, 88] . Incentives for individuals for participate in the mitigation process are also studied from the game theory perspective, see, e.g., [45, 8, 55] . It would be interesting to study how robust these control strategies with respect to the Markovian assumption, in particular, when the infectious periods are assumed to have a general distribution rather than exponential. 7.4. Open problems. There is clearly a need for more work on non-Markovian spatial epidemic models, both in discrete and continuous space. Also, endemic situations should be studied in the varying infectivity / varying susceptibility situation which we have exposed in section 6. Then the study of large and moderate deviations from the limiting deterministic LLN model opens new questions. We expect to address these questions in future work.

A standard Poisson process P (t) is a counting process (a process which counts a number of events which has happened during the interval [0, t]), which is such that P (0) = 0, P has independent increments 1 and for any 0 ≤ s < t, the law of P (t) − P (s) is Poi(t − s). Equivalently, for any t ≥ 0, the time after t until the next event is independent of what happened before t and its law is Exp(1). If P (t) is a standard Poisson process and λ > 0, P (λt) is a rate λ Poisson process (i.e., for s < t, P inf (t) − P inf (s) Poi(λ(t − s)), the waiting time until the next event after t is Exp(λ). More generally, for a deterministic function λ(t), P t 0 λ(s)ds is a rate λ(t) Poisson process. A Poisson Random Measure (abbreviated PRM) Q on a measurable set E with mean measure µ is a sum of Dirac measures at random points, which is such that the number of those points in disjoint subsets are independent, and for any measurable set A, Q(A) Poi(µ(A)). A PRM Q on a subset of R d will be called standard if its mean measure is the Lebesgue measure. Note that what we have called above a standard Poisson process is the distribution function of a standard PRM on R + .

It is not very hard to show that if λ(t) is a measurable locally bounded R + -valued function, then the two processes where P is a standard Poisson process and Q a standard Poisson random measure on R 2 + , have the same law (i.e., the same finite dimensional distributions).

Brownian motion and space-time white noise. A standard Brownian motion {B(t), t ≥ 0} is a Gaussian process with continuous paths and independent increments, and such that for any t ≥ 0, B(t) ∼ N (0, t), i.e., B(t) is a Gaussian r.v. with mean 0 and variance t. A non standard Brownian motion could have a non zero mean, and a different variance.

We use in the statement of Theorem 3.2 the notion of a white noise on R 2 . A standard white noise W on R 2 is a generalized Gaussian process {W (f ), f ∈ L 2 (R 2 )} whose law is specified by the fact that f → W (f ) is linear, and W (f ) ∼ N (0, f 2 L 2 (R 2 ) . Equivalently, for any Borel subset A ⊂ R 2 with finite Lebesgue measure, W (A) := W (1 A ) ∼ N (0, Leb(A)). A non standard white noise on R 2 is associated with a measure µ on R 2 + , such that W (A) ∼ N (0, µ(A)). The law of W (A) is specified, provided µ(A) < ∞.

8.3. The space D. In this paper, we denote by D := D([0, +∞)) the space of functions from [0, +∞) into R which are right continuous and possess a left limit at any time t > 0. Such a function is said to be càlàg, an acronym for continuà droite et limitéà gauche. If x ∈ D, whenever t n → t, with t n ≥ t for any n ≥ 1, x(t n ) → x(t), and we shall write x(t − ) for the value of lim n x(t n ), whenever t n < t for all n ≥ 1. It is not convenient to equip D with the supnorm topology, since we want that after a small modification of the time of a jump, the resulting function be close to the original one.

A sequence converges in D iff it converges in D([0, T ]) for all T > 0. It then suffices to discuss the convergence in D([0, T ]). A distance on D([0, T ]) can be defined as follows. Let Λ T denote the set of continuous strictly increasing functions from [0, T ] into itself, which map 0 into 0 and T into T . If · T denotes the supnorm on [0, T ] and I the identity mapping, a possible choice for the distance is d(x, y) = inf

The associated topology is sometimes called the Skorokhod J 1 topology. That distance makes D([0, T ]) separable. If we want D([0, T ]) to be complete, we better use a slightly different distance, whose definition is given by replacing λ − I T by sup 0≤s<t≤T log λ(t)−λ(s) t−s . It is crucial for us to have conditions under which a sequence of stochastic processes with trajectories in D is tight, which implies that such a sequence has a subsequence which converges weakly for the topology of D. Let us formulate the celebrated Aldous tightness criterion. A sufficient condition for a sequence X n of random elements of D to be tight is that the two following conditions are satisfied:

(i) For any T > 0, lim sup n P ( X n T ≥ a) → 0, as a → ∞.

(ii) For any ε, η, T > 0, there exists δ 0 > 0 and n 0 such that if δ ≤ δ 0 and n ≥ n 0 , for any discrete X n -stopping time τ ≤ T , P (|X n (τ + δ) − X n (τ )| ≥ ε) ≤ η. A proof of this criterion can be found e.g. on pages 178-179 of [15] .

In the case a-of a semi-martingale, we have a very simple criterion to verify Aldous's condition. The following is Proposition 37 in [68] : where M n (t) is a martingale, and M n t its associated predictable increasing process (i.e., M n t is predictable and |M n (t)| 2 − M n t is a martingale). If both {X n 0 }, and {sup 0≤t≤T (|ϕ n (t)| + ψ n (t))} are tight for all T > 0, then X n is tight in D. Tightness of semimartingales is usually not too hard to establish. With Markov processes are always associated martingales. However, with our non-Markov processes, we do not necessary have martingales to help us. This is why more delicate techniques are involved in the proofs for the non-Markov processes.

An optimal isolation policy for an epidemic

Optimal immunisation policies for epidemics

Asymptotic profiles of the steady states for an SIS epidemic patch model

Asymptotic profiles of the steady states for an SIS epidemic reaction-diffusion model. Discrete and Continuous Dynamical Systems

Infectious diseases of humans: dynamics and control

Limit theorems for a random graph epidemic model

Stochastic epidemic models and their statistical analysis

Optimal incentives to mitigate epidemics: a Stackelberg mean field game approach

A unified approach to the distribution of total size and total area under the trajectory of infectives in epidemic models

The final size and severity of a generalised stochastic multitype epidemic model

Stochastic SIR epidemics in structured populations

Epidemics on random intersection graphs

Epidemics and random graphs

The duration of the closed stochastic epidemic

Convergence of probability measures

Optimal control of epidemic size and duration with limited resources

React or wait: which optimal culling strategy to control infectious diseases in wildlife

A spatial stochastic epidemic model: Law of large numbers and central limit theorem. submitted

On a nonlinear integral equation for population growth problems

A mathematical model reveals the influence of population heterogeneity on herd immunity to SARS-CoV-2

Stochastic epidemics in a homogeneous community

Sir epidemic models with general infectious period distribution

An epidemic equation with immigration

Large graph limit for an SIR process in random network with heterogeneous connectivity

Limiting behaviour in an epidemic model

Optimal COVID-19 epidemic control until vaccine deployment. medRxiv

Contact rate epidemic control of COVID-19: an equilibrium view

Markov processes: characterization and convergence

Epidemiological models with non-exponentially distributed disease stages and applications to disease control

Epidemic models with varying infectivity

Estimating the state of the Covid-19 epidemic in France using a model with memory

Epidemic models with varying infectivity and immunity/susceptibility. in preparation

Household epidemic models and McKean-Vlasov Poisson driven SDEs

From individual-based epidemic models to McKendrick-von Foerster PDEs: A guide to modeling and inferring COVID-19 dynamics

SIR epidemics and vaccination on random graphs with clustering

Random perturbations of dynamical systems

Understanding and monitoring the evolution of the Covid-19 epidemic from medical emergency calls: the example of the Paris area. medRxiv preprint

On SIR epidemic models with generally distributed infectious periods: Number of secondary cases and probability of infection

Optimal control of epidemics with limited resources

Temporal dynamics in viral shedding and transmissibility of COVID-19

An SIS epidemic model with variable population size and a delay

An age dependent epidemic model

Incentives, lockdown, and testing: from Thucydides's analysis to the COVID-19 pandemic

Kermack and McKendrick revisited: the variable susceptibility model for infectious diseases

A mathematical model for Chagas disease with infectionage-dependent infectivity

Law of large numbers for the SIR epidemic on a random graph with given degrees

Near-critical SIR epidemic on a random graph with given degrees

OM Forum: COVID-19 scratch models to support local decisions. Manufacturing and Service Operations Management

A contribution to the mathematical theory of epidemics

A contribution to the mathematical theory of epidemics II. the problem of endemicity

Contributions to the mathematical theory of epidemics

Scaling Limits of Interacting Particle Systems

Dynamic games of social distancing during an epidemic: Analysis of asymmetric solutions

Large deviations for infectious diseases models

A COVID-19 epidemic model with latency period

Lyapunov functional and global asymptotic stability for an infection-age model

Two-group infection age model including an application to nosocomial infection

An introduction to mathematical epidemiology

Epidemics on networks with large initial conditions or changing structure

On the optimal control of a deterministic epidemic

A SIR model on a refining spatial grid i -Law of Large Numbers

Guodong Pang andÉtienne Pardoux. Functional limit theorems for non-Markovian epidemic models

Multi-patch epidemic models with general exposed and infectious periods

Functional central limit theorems for epidemic models with varying infectivity

Functional law of large numbers and PDEs for epidemic models with infection-age dependent infectivity

Probabilistic models of population evolution

Moderate deviations and extinction of an epidemic

Large deviation principle for epidemic models

Large deviation principle for reflected Poisson driven SDEs in epidemic models

The law of large numbers for D[0, 1]-valued random variables

Spatial deterministic epidemics

The asymptotic evolution of the general stochastic epidemic

On the asymptotic distribution of the size of a stochastic epidemic

Topics in propagation of chaos

How may infection-age-dependent infectivity affect the dynamics of HIV/AIDS?

Stochastic epidemics in a heterogeneous community

A simple SIS epidemic model with a backward bifurcation

Critical epidemics, random graphs, and Brownian motion with a parabolic drift

How generation intervals shape the relationship between growth rates and reproductive numbers

Limit theorems for age and density dependent stochastic population models

A central limit theorem for age-and density-dependent population processes

Gaussian approximation of some closed stochastic epidemic models

Asymptotic behavior of some deterministic epidemic models

Optimal immunization rules for an epidemic with recovery

Optimal isolation policies for deterministic and stochastic epidemics

Control strategies for COVID-19 epidemic with vaccination, shield immunity and quarantine: A metric temporal logic approach

A sirs epidemic model with infection-age dependence