key: cord-0118201-et7j76f4
authors: Pang, Guodong; Pardoux, Etienne
title: Functional central limit theorems for epidemic models with varying infectivity
date: 2020-09-24
journal: nan
DOI: nan
sha: 3baad3e85db9a941d61e3374514798805ec0050a
doc_id: 118201
cord_uid: et7j76f4

In this paper, we prove functional central limit theorems (FCLTs) for a stochastic epidemic model with varying infectivity and general infectious periods recently introduced in Forien, Pang and Pardoux (2020).The infectivity process (total force of infection at each time) is composed of the independent infectivity random functions of each infectious individual, which starts at the time of infection. These infectivity random functions induce the infectious periods (as well as exposed, recovered or immune periods in full generality), whose probability distributions can be very general. The epidemic model includes the generalized non-Markovian SIR, SEIR, SIS, SIRS models with infection-age dependent infectivity. In the FCLTs for the generalized SIR and SEIR models, the limits of the diffusion-scaled fluctuations of the infectivity and susceptible processes are a unique solution to a two-dimensional Gaussian-driven stochastic Volterra integral equations, and then given these solutions, the limits for the infected (exposed/infectious) and recovered processes are Gaussian processes expressed in terms of the solutions to those stochastic Volterra integral equations. We also present the FCLTs for the generalized SIS and SIRS models.

It has been observed in recent studies of the Covid-19 pandemic (see e.g. [11] ) that the infectivity of infectious individuals decreases from the epoch of symptom first appearing to full recovery. The varying infectivity characteristics also appears in many other epidemic diseases [14, 6] . We have presented a stochastic epidemic model with varying infectivity in [9] , where the various individuals have i.i.d. infectivity random functions, and the total force of infection at each time is the aggregate infectivity of all the individuals that are currently infectious. The i.i.d. assumption on the infectivity functions can be justified since we model the spread of the same disease and the infectivity is largely determined by viral load progression. We have proved a functional law of large numbers (FLLN) for the epidemic dynamics which results in a deterministic epidemic model, which is the model described as an "age-of-infection epidemic model" in [14, 5, 6] . For the early stage, we have studied the stochastic model directly starting with a small number of infectious individuals, and proved that the epidemic grows at an exponential rate on the event of non-extinction using an approximation by a non Markovian branching process. In addition, we have deduced the initial basic reproduction number R 0 from the limit process and computed its value for the case of the early phase of the Covid-19 epidemic in France. We have concluded a decreased value of R 0 induced by the decrease of the infectivity when the age of infection increases.

In this paper, we study the stochastic fluctuations of the evolution dynamics around the deterministic limits for the stochastic epidemic models with varying infectivity. The infectivity random functions can take a very general form (see Assumption 2.3) and start with a value zero for a period of time, which then in turn determines the durations of the exposed and infectious periods. The joint distribution of these two periods is determined by the law of the random function, and can be very general. Clearly, this also includes the special case with only infectious periods, where the random functions do not start from zero. We start with the generalized SIR model, where "I" represents the "infected", including either the case of only infectious periods, or the case of both exposed and infectious periods. We then study the generalized SEIR model, where "E" and "I" represent the "exposed" and "infectious" periods, which can be regarded as a more detailed model in the second above scenario by considering the separate compartments. Although it may appear this is a special case of the generalized SIR model, in many settings, one would be interested in the dynamics of the detailed "E" and "I" compartments, for example, the infectious population in isolation or hospitalization. Moreover, mathematically, for the FCLTs, one can deduce the driving Gaussian processes for the limits of the SIR model from those of the SEIR model, but not the other way around (see also Remarks 2.2 and 2.3). Therefore, for the sake of readability, we first describe the generalized SIR model and state the limit theorems, and then do the same for the generalized SEIR model, while we only provide the detailed proofs for the generalized SEIR model. Would we provide the proofs for the SIR case only, that would not simplify significantly the arguments.

For both models, in the FLLNs, we study the mean total infectivity process jointly with the proportions of the compartment counting processes. In particular, the infectivity and susceptible functions in the limit are uniquely determined by a two-dimensional Volterra integral equation, and given these two functions, the other compartment limits are given by Volterra integral formulas. In the FCLTs, we prove the joint convergence of the diffusion-scaled infectivity process with the compartment counting processes. In the limit, the diffusion-scaled infectivity and susceptible processes is determined by a two-dimensional Gaussian-driven linear stochastic Volterra integral equation. In particular, the diffusion-scaled instantaneous infectivity rate process has a limit that is a linear functional of the susceptible and infectivity limit processes. Given these, the limits of the other compartment counting processes are expressed in terms of the solution of the above stochastic Volterra integral equation. These results extend the FCLT for the classical non-Markovian SIR and SEIR models in [17] .

The main challenge in the proof of the FCLT lies in the convergence of the aggregate infectivity process. We allow the individual infectivity random functions to be piecewise continuous with a finite number of discontinuities as stated in Assumption 2.3, which covers many practical settings, see Examples 2.1-2.2. We use Poisson random measures (PRMs) induced by the laws of these random functions in the functional space D, and take advantage of some useful properties of stochastic integrals with respect to the corresponding compensated PRMs. We first give a useful decomposition of this process, and construct two auxiliary processes by replacing the random instantaneous infectivity rate process by its deterministic limit function in the FLLN. For these auxiliary processes, we employ the moment method to prove their tightness, using the criterion for tightness in [4, Theorem 13.5] . This condition is more precisely (13.14) from [4] , and is (6.3) in the Appendix. This, together with the convergence of finite dimensional distributions, proves their weak convergence. The martingale approach employed in [17] (see the proof of Lemma 4.2) could not be used for the aggregate infectivity process. It seems hard to construct appropriate martingales for the model considered in the present paper.

A major difficulty in the proof of tightness lies in the estimation of the expectation of the supremum of the square of the diffusion-scaled processes (see Lemma 3.5 and equation (3. 36)), since we cannot use a semimartingale decomposition and exploit Doob's maximal inequality. In order to establish this, we have developed a new criterion to interchange sup t and the expectation of the square of a stochastic process, provided that the process satisfies the well-known tightness criterion using the increments at three time points: see equation (3. 2) as in Theorem 6.1 (which is a variant of Theorem 13.5 of [4] ), as well as the sufficient moment criterion in (3.4) ). This result is stated in Theorem 3.1. It should be of independent interest and we expect that it will be useful for other purposes. For our model, we employ the moment formulas (specifically the mixed second moment formula in (3.7)) for the stochastic integrals with respect to the PRMs to establish the tightness moment bound for the increments of the infectivity process (see Propositions 3.1 and 3.2).

In addition, another challenge is to prove the joint convergence of the aggregate infectivity processes and the compartment counting processes. We have developed an approach using a common PRM to represent the subprocesses associated with the newly infected individuals in the expressions of the infectivity and compartment processes, and then carried out calculations with characteristic functions of these processes by making use of the properties of PRMs (See the proof of Lemma 4.2) .

We also state the FCLTs without proofs for the generalized SIS and SIRS models with varying infectivity (the proofs follow with slight modifications). For the SIS model, the epidemic dynamics is determined by the aggregate infectivity process and the infectious process, whose limits are given by a two-dimensional Gaussian-driven linear stochastic Volterra integral equation (Theorem 5.1) . For the SIRS model, the epidemic dynamics is determined by the aggregate infectivity process and the infectious and recovered/immune processes, whose limits are given by a three-dimensional Gaussian-driven linear stochastic Volterra integral equation (Theorem 5.2).

This work contributes to the literature of stochastic epidemic models in the aspects of infectionage dependent infectivity, and general infectious periods. It establishes a useful result describing the fluctuations of the stochastic individual based model around its law of large numbers limit. The existing work on epidemic models with infection-age dependent (varying) infectivity has all been about the deterministic models, pioneered by Kermack and McKendrick in 1927 [14] and more recently proposed by Brauer [5] (see also [6, Chapter 4.5] ), and the PDE models (see, e.g., [12, 20, 13, 15] ). These were not considered as the FLLN limits of a well specified stochastic model either, except our work in [9] . For non-Markovian epidemic models with general infectious periods, although some deterministic models (including Volterra integral equations) appeared in the literature (see, e.g., [6, Chapter 4.5] and references therein), the works that rigorously establish them as a FLLN from a stochastic model are limited. Both FLLNs and FCLTs are proved in [21, 22, 23] for population models, including an epidemic SIR model which has an infection rate dependent on the number of infectious individuals and general infectious periods for newly infected individuals, and which also has the remaining infectious periods of the initially infected individuals dependent on the elapsed infectious durations. A FLLN for the associated measure-valued processes was established for the SIR model with the infection rate dependent on time and the number of infected individuals using Stein's method in [18] . For the general SIS, SIR, SEIR and SIRS models with a constant infection rate and the associated multi-patch models, both FLLNs and FCLTs are recently established in [17, 16] . The results in [17] are similar to those in [21, 22, 23] for the SIR model with general infectious periods, since the limits are given as Volterra integral equations, deterministic in the FLLNs and stochastic driven by Gaussian processes in the FCLTs. However, the proofs are rather different, and the FCLTs in [17] do not require any condition on the c.d.f. of the infectious periods, while those in [21, 22, 23 ] assume a C 1 condition. As mentioned at the beginning, the model with varying infectivity studied in this paper is much more general than that in [17] . It is also quite different from the models in [21, 22, 23, 18] since an infection rate dependent on the number of infectious individuals captures neither the randomness of infectivity as the i.i.d. random infectivity function of each individual, nor the age of infection dependent infectivity of each individual. Distinct from [17] , this work focuses on the convergence of the aggregate infectivity process and its joint convergence with the compartment counting processes, which is much more involved as discussed above. For other aspects of epidemic models with general infectious periods, including the Sellke construction and the final size of an epidemic, we refer the readers to [19, 1, 2, 3] and the recent survey [7] .

The paper is organized as follows. In Section 2.1, we first describe the generalized SIR model with varying infectivity and state the limit theorems. Here the infected periods possibly include the combined exposed and infectious durations. We then describe a generalized SEIR model with separate exposed and infectious compartments in Section 2.2, and state the limit theorems for this model. We provide a detailed proof for the FCLT of the generalized SEIR model in Sections 3 and 4. We first state several technical preliminaries in Subsection 3.1, including the criterion for moment estimate of the supremum of stochastic processes in D, the formulas of the Laplace functional and moments of stochastic integrals with respect to PRMs, and an estimate of the increments of the infectivity function. We then proceed to the proofs of the convergence of the number of susceptible individuals and the aggregate infectivity process in Subsections 3.2-3.6, with a proof roadmap given in Subsection 3.2. The joint convergence of the exposed, infectious and recovered processes together with the two processes mentioned above are given in Section 4. Finally, the FCLTs for the generalized SIR and SIRS models are stated in Sections 5.1 and 5.2, respectively.

1.1. Notation. Throughout the paper, N denotes the set of natural numbers, and R k (R k + ) denotes the space of k-dimensional vectors with real (nonnegative) coordinates, with R(R + ) for k = 1. For x, y ∈ R, denote x ∧ y = min{x, y} and x ∨ y = max{x, y}. Let D = D([0, +∞), R) denote the space of R-valued càdlàg functions defined on [0, +∞). Throughout the paper, convergence in D means convergence in the Skorohod J 1 topology, see chapter 3 of [4] . Also, D k stands for the k-fold product equipped with the product topology. We use C to denote be the subset of D consisting of continuous functions and C 1 the subset of differentiable functions whose derivative is continuous.

All random variables and processes are defined on a common complete probability space (Ω, F, P). The notation ⇒ means convergence in distribution. The notation x y means that there exists some constant c such that x ≤ cy for x, y ∈ R. We use 1(·) or 1 {·} for indicator function.

We write F (t) = t 0 F (ds) for a cumulative distribution function (c.d.f.) F on R + . For any measure µ on R and f a measurable and µ-integrable function, the integral b a f (t)µ(dt) will mean (a,b] f (t)µ(dt). We use the following conventional notation in the paper: for the scaled processes and quantities (all indexed by the population size N ) and their limits in the FLLN, we use an upper bar, and for those in the FCLT, we use a hat. This means that for any process X N , we letX N = N −1 X N , andX N := √ N (X N −X). For the processes and quantities associated with the initial quantities, we use a super index '0'. The letters S, E, I, R with different indices (e.g., index '0' for the initial quantities) represent the Susceptible, Exposed, Infectious and Recovered compartments. In the case of the SIR model, I stands for Infected.

2.1. Generalized SIR model with varying infectivity. In this subsection, S stands for the compartment of susceptible individuals, I for infected (either exposed, and not yet infectious, or infectious) and R for Recovered. Those who are recovered have a permanent immunity, and do not play a role in the epidemic anymore. Those who died from the disease are placed in that R compartment.

Let the population size be N , and S N (t), I N (t), R N (t) be the numbers of susceptible, infected and recovered individuals at each time t, respectively. We have the balance equation N = S N (t) + I N (t) + R N (t) for t ≥ 0. We assume that R N (0) = 0, S N (0) > 0 and I N (0) > 0. If there are recovered individuals at time 0, we can suppress them from the population, since they play no role in the epidemic.

An individual who gets infected (i.e., jumps from the S to the I compartment) at time t > 0 will have at time s > t an infection age infectivity λ(s − t). The function λ is random, in the sense that each individual experiences a specific realization of that random function (due to the specificity of each individual's immune system reaction to the illness). We assume that the realizations associated to the various individuals are i.i.d. Let χ = inf{s > 0, λ(r) = 0, for all r > s}.

(2.1)

The time interval [t, t + χ) is the time interval during which the individual is infected (including either exposed or infectious). At time t+χ the individual jumps from the I into the R compartment. Let {λ i (t), t ≥ 0} denote the infectivity function of the i-th infected individual after time 0. We assume that the collection of random functions {λ i } i≥1 is i.i.d., and hence, the corresponding infected periods {χ i } i≥1 are also i.i.d. The I N (0) individuals who are infected at time 0 are thought of as having been infected in the past (at some time t < 0). Therefore their random infectivity functions {λ 0 j (t), t ≥ 0} 1≤j≤I N (0) , starting at time t = 0, are i.i.d., and their common law differs for that of λ i . We also assume that the random functions λ 0 1 , . . . , λ 0 I N (0) , λ 1 , λ 2 , . . . are independent. Also, let us define, for j = 1, . . . , I N (0), Remark 2.1. One possible choice for {λ 0 j , 1 ≤ j ≤ I N (0)} is as follows. Suppose we are given a sequence of i.i.d. random functions {λ i } i∈Z . As above, λ i (t), for i ≥ 1, describes the infectivity function of the i-th infected individual after time 0. To each 1 ≤ j ≤ I N (0), we associate the random function λ −j (t), and a positive random variable κ j , which represents the infection age at time 0 of the initially infected individual j, which we assume to satisfy κ j < χ −j . For 1 ≤ j ≤ I N (0), we let λ 0 j (t) = λ −j (κ j + t). For t > 0, let A N (t) denote the number of individuals who have been infected on the time interval (0, t]. This process can be expressed as

is the instantaneous infection rate function at time t, and Q is a standard Poisson random measure (PRM) 1 on R 2 + . In the formula for Υ N (t), I N (t) is the total force of infection which is exerted on the susceptibles at time t. If {τ N i } i≥1 denote the successive times of infection (i.e., the jump times of the counting process A N (t)),

(2.5) We can now precise the evolution of the numbers of infected and recovered individuals. Now the numbers of infected and recovered individuals at time t, I N (t) and R N (t), are given by 

We also make the following assumptions on the random infectivity functions. See the discussions on the assumptions in Section 2.3.

Assumption 2.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, satisfying the following properties.

(i) There exists a constant λ * < ∞ such that sup t∈[0,T ] max{λ 0 (t), λ(t)} ≤ λ * almost surely.

(ii) There exist a given number k ≥ 1, a random sequence 0 = ξ 0 < ξ 1 < · · · < ξ k = +∞ and random functions λ j ∈ C, 1 ≤ j ≤ k such that

In addition, for any T > 0, there exists a deterministic nondecreasing function

In [9] , the following functional law of large numbers (FLLN) is established for the renormalized quantities

in probability, whereS

(s)Ī(s)ds,

(t − s)S(s)Ī(s)ds,

The purpose of the present paper is to establish a functional central limit theorem (FCLT) for the fluctuations of the stochastic sequence around its deterministic limit. We define the renormalized differences

We make the following assumption on the initial conditions. Before we state the FCLT for the generalized SIR model, let us state the following definition. Definition 2.1. Let (Ŵ S ,Ŵ I ,Ŵ I ,Ŵ R ) is a four-dimensional centered Gaussian process, independent ofÎ(0), with covariance functions: for t, t ′ ≥ 0,

We shall establish below the following FCLT result, for which we impose the following further conditions on the random function λ(t). Assumption 2.3. In addition to the conditions in Assumption 2.1, the following conditions hold.

(i) There exist nondecreasing functions φ and ψ in C and α > 1/2 and β > 1 such that for all 0 ≤ r ≤ s ≤ t, denotingλ 0 (t) = λ 0 (t) −λ 0 (t),

(ii) The function ϕ T from Assumption 2.1 (ii) satisfies

Also, if F j denotes the c.d.f. of the r.v. ξ j , there exist C ′ and ρ > 0 such that for any 12) and in addition, for any 1 ≤ j ≤ k − 1, r > 0, 

where the process Ŝ ,Î is the unique solution of the SDÊ

(s)ds,

(2.14) Moreover, the pair Î N ,R N is given aŝ

Generalized SEIR model with varying infectivity. In this section, we consider a more detailed model by studying the separate compartments for the exposed and infectious individuals, which were combined into the infected compartment in the previous section. Changing our notation, we use now I N (t) to denote the number of infectious individuals at time t in this model. And we add E N (t) to denote the number of exposed individuals at time t. The processes S N (t) and R N (t) again represent the numbers of susceptible and recovered individuals at time t, respectively. We have the balance equation

In this model, the infected period χ in (2.1) is split into the exposed and infectious periods: for each newly infected individual i,

That is,

However, for the initially infected individuals, we encounter two groups: initially exposed individuals (going through the remaining exposed period and then the infectious period) and initially infectious individuals (only going through the remaining infectious period). Unlike the newly infected individuals as described above in section 2.1, we will use two separate collections of variables and quantities to describe the dynamics of these two groups, which is necessary.

Let {λ 0 j (·)} and {λ 0,I k (·)} be the infectivity processes associated with each initially exposed and initially infectious individual, respectively. Assume that the sequence {λ 0 j (·)} is i.i.d., and so are {λ 0,I k (·)}. The remaining exposed period and the infectious period (ζ 0 j , η 0 j ) of an initially exposed individual are given as: 16) and the remaining infectious period η 0,I k of an initially infectious individual (note that inf{t > 0, λ 0,I k (t) > 0} = 0) is:

(2.17) Under the i.i.d. assumptions of the corresponding infectivity processes, the random vectors

, and so is the sequence {η 0,I k : k ∈ N}. In addition, these three sequences are mutually independent.

Let H(du, dv) denote the law of (ζ, η), H 0 (du, dv) that of (ζ 0 , η 0 ) and F 0,I the c.d.f. of η 0,I . Define

(Note that Φ(t) and Φ 0 (t) were used for the infected periods in the previous section for the generalized SIR model; if the infected periods have both exposed and infectious periods, then they have the same meanings as defined above in the generalized SEIR model. See Remark 2.3 for further discussions.

is the conditional law of η 0 , given that ζ 0 = u. In the case of independent exposed and infectious periods, it is reasonable that the infectious periods of the initially exposed individuals have the same distribution as the newly exposed ones, that is, F 0 = F . Note that in the independent case,

We again use A N (t) to denote the number of individuals that are exposed in (0, t], and τ N i to denote the time at which the i th individual gets exposed. Then A N (t) has the same expression as in (2.2). However, by taking into account the detailed initial conditions with initially exposed and infectious individuals, the total force of infection I N (t) at time t is expressed as

The epidemic dynamics of the model can be described by

Remark 2.2. We remark that in addition to the extra care in the initial conditions, since we only count the infectious individuals in I N (t), the third term in the expression of I N (t) differs from the second term in the expression of I N (t) in (2.6) in the generalized SIR model.

We remark that the model in the previous section includes the special case without any exposed periods. Specifically, the random infectivity function λ i (t) does not equal to zero immediately after time 0, that is, an infected individual is immediately infectious, so ζ = ζ 0 = 0 a.s., and there are no exposed individuals, E N (t) = 0 for all t ≥ 0. In that case, the FLLN and FCLT results in Theorems 2.1 and 2.2 hold with Φ being equal to F , the c.d.f. of the infectious period in this section. On the other hand, the results for the SIR model without exposed periods can be also obtained by setting G to have mass 1 at 0 so that Φ(t) = F (t) andĒ(0) = 0 in the following two theorems.

In [9, Theorem 2.1], the following FLLN is proved.

in probability. The limit (S,Ī) is the unique solution of the following system of integral equations:

and the limits (Ē,Ī,R) are given by the following formulas:

We also haveῩ N →Ῡ in D in probability as N → ∞, wherē

We make the following assumption on the initial quantities for the FCLT.

and Cov(Ŵ S (t),Ŵ S (t ′ )), Cov(Ŵ S (t),Ŵ I (t ′ )), and Cov(Ŵ S (t),Ŵ R (t ′ )) are the same as those in Definition 2.1. 

The limit process (Ŝ,Î) is the unique solution to the following system of stochastic integral equations:

andS(t) andĪ(t) are given by the unique solutions to the integral equations (2.20) and (2.21). The limit processesÊ,Î andR are given by the expressions:

S has continuous paths, and ifλ 0 andλ 0,I are in C, which implies that G 0 and F 0,I are continuous, thenÎ,Ê,Î andR are continuous.

On the assumptions on λ(t). We remark that the conditions in Assumption 2.3 are not required to establish the FLLN [9] . It is not surprising that the FCLT requires additional assumptions, compared with the FLLN.

These additional conditions are required to establish the moment criterion for tightness of the aggregate infectivity process (see Propositions 3.1 and 3.2 for the moment bounds for the increments of the associated processes). In fact, condition (2.13) on {ξ j } is used only once in the end of the proof of Proposition 3.1, while (2.12) is used twice in the proofs of both Proposition 3.1 and 3.2. The conditions in (iii) allow fairly general random infectivity functions λ(t), which can have càdlàg paths. They evidently impose a restriction on the distribution functions of the exposed and infectious periods in the case where λ jumps at the end of those periods, which in particular implies that they are continuous. Recall that no condition is required on the distribution functions of the exposed and infectious periods in establishing the FLLNs and FCLTs of the standard SIR and SEIR models in [17] . The restriction on these distributions as a result of (2.12) and (2.13) is a compromise of allowing the random infectivity functions to have discontinuities.

On the other hand, as a special case with k = 1, λ(·) is continuous satisfying (2.11) , no conditions of the type (2.12) and (2.13) has to be imposed. In that case, one can allow the durations of the exposed and infectious periods to take discrete values, and as a consequence, the same results of the FLLNs and FCLTs hold without requiring any condition on the distribution functions of those exposed and infectious periods.

We now describe examples of λ's which satisfy the conditions in Assumption 2.3. We also refer the reader to [9, Section 2.5] for an example of λ adapted to the Covid-19 epidemic.

Example 2.1. Suppose we are given a random element (ζ, η) of (0, +∞) 2 , and let λ be a random function with trajectories in C which is such that

and which is bounded and satisfies (2.11). Then such a λ satisfies all the conditions in Assumption 2.3, and the law of the pair (ζ, η) is completely arbitrary. This corresponds to the case k = 1 in Assumption 2.3 (ii).

Example 2.2. Suppose we are given again a random element (ζ, η) of (0, +∞) 2 , and λ a random function with now trajectories in D of the form

where λ 1 ∈ C is bounded and satisfies (2.11). Note that unlike in the above example, we allow λ 1 (ζ) > 0 and λ 1 (ζ + η) > 0, hence λ(ζ) > λ(ζ − ) = 0 and λ((ζ + η) − ) > λ(ζ + η) = 0. In this case, Assumption 2.3 (ii) is satisfied with k = 3, ξ 1 = ζ and ξ 2 = ζ + η, provided the distribution functions of both ζ and ζ + η satisfy (2.12), and (2.13) holds. We can of course add more jumps on the time interval (ξ 1 , ξ 2 ) = (ζ, ζ + η). In the particular case λ 1 (ζ) = 0 but λ 1 (ζ + η) > 0, (resp. λ 1 (ζ) > 0 but λ 1 (ζ + η) = 0), we are in the case of Assumption 2.3 (ii) with k = 2 and ξ 1 = ζ + η (resp. ξ 1 = ζ), and condition (2.12) needs to be imposed to that r.v. only, while condition (2.13) is irrelevant in this case. Finally if λ 1 (ζ) = λ 1 (ζ + η) = 0, we are back in the previous example.

3. Proofs for the Convergence of (Ŝ N ,Î N )

In this section we prove the convergence of (Ŝ N ,Î N ) → (Ŝ,Î) in D 2 in Theorem 2.4. The joint convergence with (Ê N ,Î N ,R N ) will be studied in the next section. We first provide the following technical preliminaries used in the proofs.

3.1. Technical Preliminaries. In this subsection, we gather a few technical results which will be used later in our proofs. Recall the moment criterion (6.3) (which is (13.14) of [4] . We will first deduce from that condition a uniform moment estimate.

Suppose now that a sequence {X N } satisfies the condition (6.2) (which is (13.13) in [4] ). Under an additional mild condition, it implies that {X N } is tight in D, hence in particular that for any T > 0, sup 0≤t≤T |X N (t)| is a tight sequence of R + -valued r.v.'s. We now show how we can deduce a bound on the second moment of that supremum.

. Suppose that the two following conditions are satisfied

and that for some α > 1/2, β > 1/2, and for all 0

where G is a nondecreasing continuous function satisfying G(0) = 0. Then

As noted in [4] , a sufficient condition for (3.2) is that

We shall use this result with β = 1. Condition (3.2) with any β ≥ 0 plus a minor condition, much weaker than the following reinforced version of (3.1):

implies tightness of the sequence X N . Note that going from (3.5) to (3.3) is usually easy to achieve if X N is a semimartingale, using Doob's inequality for the martingale part. In our non Markovian setup, we do not have enough martingales at our disposal, and the above theorem will be useful to us.

Theorem 10.3 in [4] says that, provided that the process X N (t) is right continuous, (3.2) implies that there exists a constant K α,β such that for all N ≥ 1, x > 0,

where we have used (3.6) for the second inequality. The result now follows from (3.1) and the fact that 2β > 1.

We shall also use the following technical lemma.

Proof. The result follows from a combination of several results from [4] . First the Corollary of Theorem 7.4 says that our assumption implies that for all T > 0,

Combined with the inequality (12.7), that last fact implies that the second condition of Theorem 13.2 applies. Moreover since our assumption implies that the jumps of X N tend to 0 as N → ∞, and X N (0) = 0, condition (i") from the Corollary of Theorem 13.2 applies, which replaces the first condition of the Theorem, hence that Theorem applies, and the sequence X N is tight in D([0, T ]) for all T > 0, hence in D.

Let us recall well-known formulas for the exponential moment and moments of an integral with respect to a compensated PRM. These follow, e.g., rather easily from Theorem VI.2.9 and Exercise VI.2.21 in [8] .

We will need the following bounds on the increments of the infectivity function λ(·), which are used below in the proofs of Propositions 3.1-3.2 and Lemmas 3.8-3.9, Lemma 3.3. For t ≥ s ≥ 0, with α > 1/2 as in Assumption 2.3 (ii),

Moreover,

Thus the first statement follows from Assumption 2.3 (ii). The other statements follow readily from the first one.

We shall also use several times the following straightforward inequality: for

(3.8)

3.2. Representation of (Ŝ N ,Î N ) and proof roadmap. By the expressions of I N in (2.18) and I in (2.21), we obtain

For the process A N (t), we have the decomposition

with Q(ds, du) := Q(ds, du) − dsdu being the compensated PRM. The process {M N A (t) : t ≥ 0} is a square-integrable martingale with respect to the filtration {F N t : t ≥ 0} defined by

It has the quadratic variation (see e.g. [8, Chapter VI] By (3.12), we haveÂ

It then follows that

Roadmap of the proofs ahead. In what follows, we will first prove the joint convergence of Î N 0,1 ,Î N 0,2 in Lemma 3.4, which is analogous to the proof of CLT for i.i.d. random processes, but with a variation of a random summation term I N (0) or E N (0). We then prove the joint convergences of (Î N 1 ,Î N 2 ). We observe that both processesÎ N 1 andÎ N 2 can be written as stochastic integrals with respect to the associated PRMs:Î N 1 uses a PRM on R 2 + × D, whose mean measure is the product of the Lebesgue mesure on R 2 + with the law of λ, whileÎ N 2 uses a standard PRM on R 2 + . We employ the moment criterion (6.3) to establish the tightness of these two processes which requires calculations of moments of the process increments. In order to apply the formulas in Lemma 3.2, we introduce two associated processes I N 1 and I N 2 by replacing the stochastic intensity by a deterministic function, and provide the required moment estimates in Section 3.4. The next crucial step is to prove the finiteness of the second moment of the supremum of the processesŜ N andÎ N , which results in the same property ofΥ N (Lemma 3.5). For that, we employ the new criterion in Theorem 3.1 to first establish the bounds for the supremum of the processes I 

where (Î 0,1 ,Î 0,2 ) is a centered 2-dimensional Gaussian process with the following covariance functions: for t, t ′ ≥ 0, Cov(Î 0,1 (t),Î 0,2 (t ′ )) = 0, and

By the CLT for the random elements in D (see Theorem 2 in [10] , whose conditions (i) and (ii) are satisfied thanks to Assumption 2.3 (i) (a) and (b), respectively) and by the independence of the sequences {λ 0 j } j≥1 and {λ 0,I k } k≥1 , we obtain

It then suffices to show that

in probability. We focus on I N 0,2 −Î N 0,2 ⇒ 0. It is clear from the definition in (3.19) and the i.i.d. property of λ 0 j (·) that for each t ≥ 0, E I N 0,2 (t) −Î N 0,2 (t) = 0, and

where v 0 (t) < ∞, and the convergence follows from Assumption 2.4 and the dominated convergence theorem. It then remains to show that { I N 0,2 −Î N 0,2 : N ∈ N} is tight in D. We have

We use the moment criterion (6.3), and consider the moment: for t ′ ≤ t ≤ t ′′ ,

(3.23) Recallλ 0 j (t) = λ 0 j (t) −λ 0 (t), and we drop the subscript j for the generic variableλ 0 (t). Then by the i.i.d. and mean zero properties ofλ 0 j (t), and by the independence betweenĒ N (0) andλ 0 j (t), we obtain that the moment above is equal to 3.4. Moment estimates associated with the processesÎ N 1 ,Î N 2 andΥ N . We first consider I N 1 . Let us introduce the canonical probability space associated to the random function λ. If we equip D with its Borel σ-field, and the probability measure which is the law of the random function λ, then the identity mapping λ → λ from D into itself defines our random function λ. So that below λ will be both the same object as defined above, and the generic element in D. As such, it will also be a dummy variable in the next integral. We introduce a PRMQ on R + × D × R + , which to the point τ N i associates the copy λ i of the random function λ, so that the mean measure of the PRM is ds × Law of λ × du . With the notationλ := λ −λ,Î N 1 can be written

We note that if we replace in the aboveQ by its mean measure, then the resulting integral vanishes. Consequently we also havê

where Q is the compensated PRM ofQ. Hence E[Î N 1 (t)] = 0 and

In the next proposition, we prove the moment bound for the increments of I N 1 , which will be used in the proof of its convergence. Note that by (3.13) and (3.14), we haveῩ(t) ≤ λ * for any t ≥ 0. Proposition 3.1. There exist C > 0 and β > 0 such that for any r < t < v,

Proof. We deduce from (3.7) the following identity

(3.30)

We now bound each of the four terms on the right of (3.30). It follows from Lemma 3.3 and (2.11) and (2.12) in Assumption 2.3 that for some constant C, the first term on the right of (3.30) is bounded by

We next consider the third term. The absolute value of the first summand in the square is clearly bounded by (λ * ) 3 (t − r). The second summand is bounded by

where we have used twice (3.8) . Finally the third term is bounded by a constant times

By similar arguments, the first factor in the last term is bounded by a constant times

Since the second factor can be estimated in the same way, with however t − r replaced by v − t, we conclude that the last term is again bounded by a constant times (v − r) 2 + (v − r) 4α . It remains to consider the second term, which is a bit more delicate. We disregard the factor 1/N . s) . We disregard the factor 1/N . The integrand in the integral from s = 0 to s = r is bounded from above by a constant times

The first term on the right hand side of the last inequality is bounded by (v − r) 4α . The expectation of the second term is bounded by (t − r) 2α P(A i (s)), whose integral from s = 0 to s = r is bounded by (t − r)(v − r) 2α , and the next term is treated exactly in the same way. We finally treat the last term. We note that (1 A i (s) − P(A i (s))) 2 (1 B j (s) − P(B j (s))) 2 is the square of a real number which is between −1 and +1, hence it is bounded by its absolute value. Consequently

It remains to bound each of those two last terms. Note that P(A i (s)∩B j (s)) = 0 only if i < j ≤ k−1, and in that case

thanks to (2.13) in Assumption 2.3. Since by (3.8), r 0 P(A i (s))ds ≤ t − r, the integral from s = 0 to s = r of the first term has the right bound. The second term is bounded by the same expression, using this time (2.12) instead of (2.13).

Concerning the processÎ N 2 (t), we can represent it aŝ

In the next proposition, we prove the moment bound for the increments of I N 2 . Proposition 3.2. There exist C > 0 and β > 0 such that for any r < t < v,

Proof. As in the proof of Proposition 3.2, we exploit (3.7). We obtain

For the first term on the right hand side, by Lemma 3.3 as well as the assumptions in (2.11) and (2.12), we obtain the upper bound: for some constant C > 0,

For the second term, we obtain the upper bound

where we have used (2.12) and (3.8) . For the third term, we obtain

Finally, we get the following bound for the last term,

Therefore we have obtain the desired upper bound, with β = (4α − 1) ∧ ρ ∧ 1 ∧ 2α. Proof. We first show that for any T > 0, 

We shall use (3.15 ) and the two integral representations in (3.9) and (3.16) . We first obtain the following estimates. It is clear that sup N sup t∈[0,T ] E |M N A (t)| 2 ≤ λ * T and from Assumption 2.4, there exists a constant C > 0 such that for all N ,

It is also clear that 

Combining (3.9), (3.15) and (3.16) with the simple boundsS(t) ≤ 1 andĪ N (t) ≤ λ * (Ī N (0) + A N (t)) ≤ 2λ * , and Gronwall's inequality, we obtain the claims in (3.37).

We next prove (3.35). By Doob's maximal inequality, we have sup ForÎ N (t), we first consider the two processesÎ N 0,1 (t) andÎ N 0,2 (t), which can be treated in the same way, so we focus onÎ N 0,2 (t) as in the proof of Lemma 3.4. RecallĨ N 0,2 (t) in (3.19) . Under Assumption 2.3 (ii), we deduce from a computation similar to the one leading to (3.24) and Theorem 3.1 that

(3.39)

We next considerÎ N 0,2 (t) −Ĩ N 0,2 (t) which is given in (3.22) . By (3.23) and (3.24), applying again Theorem 3.1, we obtain that the moment property in (3.39) holds forÎ N 0,2 (t) −Ĩ N 0,2 (t). Combining these two moment estimates, we deduce that the moment property in (3.39) holds forÎ N 0,2 (t).

For the processesÎ N 1 (t), we write it asÎ N 1 (t) = I N 1 (t)+(Î N 1 (t)− I N 1 (t)) and treat the decomposed terms separately. By Proposition 3.1, applying Theorem 3.1, we obtain that the moment property in (3.39) holds for I N 1 (t). Now for the differenceÎ N 1 (t) − I N 1 (t), we havê

Combining the above arguments, we obtain that the moment property in (3.39) holds forÎ N 1 (t). Similarly, we writeÎ N 2 (t) = I N 2 (t) + (Î N 2 (t) − I N 2 (t)) and treat the decomposed terms separately. We obtain the moment property in (3.39) for I N 2 (t) by Proposition 3.2 and Theorem 3.1. For the differenceÎ N 2 (t) − I N 2 (t), we havê

Both terms on the right hand side are increasing in t. Thus,

Here the first integral converges to zero as N → ∞, and the second term is bounded by (3.38). Thus combining these, we have shown that the moment property in (3.39) holds forÎ N 2 (t). Finally, combining (3.9) with the above estimates yields thatÎ N satisfies (3.35).

In this subsection, we will show the following result. Lemma 3.6. Under Assumptions 2.3 (ii) and 2.4,

where (Î 1 ,Î 2 ) is a centered Gaussian process with covariance functions: for t, t ′ ≥ 0, Cov(Î 1 (t),Î 2 (t ′ )) = 0 and

We shall first show that (

Given the results in Propositions 3.1 and 3.2, Theorem 6.1 tells us that it remains to prove that the finite dimensional distributions of ( I N 1 , I N 2 ) converge to those (Î 1 ,Î 2 ), which we establish in the first Lemma which follows. It will then remain to prove thatÎ N 1 − I N 1 → 0 andÎ N 2 − I N 2 → 0 in D in probability, as N → ∞ which will be done in the next two Lemmas.

Lemma 3.7. For any k ≥ 1, 0 < t 1 < t 2 < · · · < t k , as N → ∞,

Proof. Recall (3.28) and (3.32) . Note that the process I N 2 (t) can be equivalently written as follows using Q:

In order to simplify the notations, we start with the case k = 2. For any θ 1 , θ 2 , θ ′ 1 , θ ′ 2 ∈ R, t, t ′ > 0, we compute the limit as N → ∞ of

We apply Lemma 3.2 to the particular case E = R + × D × R + , Q =Q and

We have proved that the two dimensional distributions of ( I N 1 , I N 2 ) converge to those of a centered Gaussian process, whose covariances are the ones of (Î 1 ,Î 2 ) as given in Theorem 2.4. Note that it is immediate that the k dimensional distributions also converge to those of a centered Gaussian process, and the covariances are fully determined by the above computation, hence the result. 

Here the convergence follows from |λ(t − s)|Q(ds, dλ, du). 

We have

N (Ῡ N (s)∧Ῡ(s)) Q(ds, du)

where the second inequality follows from Lemma 3.3. It is clear that the above upper bound is increasing in v. Thus, we obtain that for any ǫ > 0,

Q(ds, du) > ǫ/4

N (Ῡ N (s)∧Ῡ(s)) Q(ds, du) > ǫ/4

The first term is bounded by 16

where the first inequality follows from Q(ds, du) = Q(ds, du) + ds × du. We note that the first term on the right of (3.48) tends to 0 as N → ∞, while the lim sup N of the second term multiplied by δ −1 tends to 0, as δ → 0, by the moment property ofΥ N in Lemma 3.5, which is exactly what we want.

The second term is bounded by 16

By (2.19) , the first term converges to zero as N → ∞, while, thanks to the fact that 2α > 1 and Lemma 3.5, lim sup N of the second term multiplied by δ −1 tends to 0, as δ → 0, which again is exactly what we want. The third term on the right hand side of (3.47) is bounded by 16

Here the first term is equal to twice

which converges to zero as N → ∞ by (2.19) . The second term satisfies, thanks to (3.8) and Lemma 3.5,

The fourth and last term on the right hand side of (3.47) is bounded by 16

Here the first term converges to zero as N → ∞ by (3.45). The second term also satisfies (3.49). It is then clear that (3.46) holds for Ξ N . This completes the proof. 

By Lemma 3.1, it suffices to show that for ℓ = 1, 2,

For the process Ξ N 1 (t), we have

Q(ds, du)

We already know how to treat the first term, see (3.48). By Lemma 3.3, the second term on the right hand side is bounded by

which is nondecreasing in v. Thus, we obtain that for any ǫ > 0,

Q(ds, du) > ǫ/3

N (Ῡ N (s)∧Ῡ(s)) Q(ds, du) > ǫ/3

The first term is bounded as in (3.48). The second term is bounded by 9

This upper bound satisfies the proper bound (3.50), by the same argument as already used in the proof of Lemma 3.8. The third term on the right hand side of (3.51) is bounded by 9

Here the first term converges to zero as N → ∞ by (2.19) . The second term satisfies (3.49). Next for the process Ξ N 2 (t), we have

The first term is bounded from above by

while the second term on the right hand side can be bounded by

Those two upper bounds are nondecreasing in v, and by already used arguments, we easily establish that Ξ N 2 satisfies (3.50). This completes the proof of the lemma. 3.6. Completing the proof of the convergence (Ŝ N ,Î N ) → (Ŝ,Î) in D 2 . We are now ready to complete the proof of the convergence of (Ŝ N ,Î N ) → (Ŝ,Î), stated in Theorem 2.4.

Proof of the convergence (Ŝ N ,Î N ) → (Ŝ,Î) . We first prove the joint convergence

in R 2 × D 5 as N → ∞. By the independence of the variables associated with the initially and newly infected individuals, it suffices to show the joint convergences Recall the process I N 1 (t) defined in (3.28), whereQ can be replaced by Q. Also, recall the process I N 2 (t) in (3.43) using Q. A minor extension of the computation done in Lemma 3.7 yields the joint convergence of the finite dimensional distributions of the processes M N A , I N 1 , I N 2 . Note that M N A being a martingale, its tightness is easily established. Hence we deduce the joint convergence

Concerning the pair (M A ,Î 2 ), we also observe thatM A is a non-standard Brownian motion, and I 2 (t) = The mapping which to those data associates the solution is continuous in the Skorohod J 1 topology, see Lemma 8.1 in [17] . Thus, by the joint convergence in (3.52), we apply the continuous mapping theorem to conclude (2.23) .

It is clear from (3.9) that the limit termŴ I (t) in the expression ofÎ(t) in Theorem 2.4 is equal toŴ I (t) =Î 0,1 (t) +Î 0,2 (t) +Î 1 (t) +Î 2 (t). Given the above results, it is straightforward to derive the covariance functions ofŴ I (t). Specifically, Cov(Ŵ

Finally we show that the limit processesÎ 1 andÎ 2 have a continuous version in C, given their consistent finite dimensional distributions. Since its increments of the processes are Gaussian and centered, it suffices to show the continuity of the covariance functions. It is easy to calculate that

Then we can verify the continuity property under the conditions in Assumption 2.3 (iii).

In this section we complete the proof of Theorem 2.4 by establishing the convergence of (Ê N ,Î N ,R N ) and their joint convergence with of (Ŝ N ,Î N ). The former follows essentially the same arguments as that of Theorem 3.2 in [17] , for which we only highlight the differences. However, the joint convergence is nontrivial since the exposed and infectious periods are induced by the random infectivity functions.

We have the following representations for the processes (Ê N ,Î N ,R N ):

,

, 

jointly with the convergence Î N 0,1 ,Î N 0,2 ⇒ Î 0,1 ,Î N 0,2 in (3.17). (Î 0,1 ,Î 0,1 ,R 0,1 ) and (Î 0,2 ,Ê 0 ,Î 0,2 ,R 0,2 ) are independent three-dimensional centered Gaussian processes, independent of (Ê(0),Î (0)), driven by the common random source λ 0,I (·) and λ 0 (·), respectively. (Î 0,1 ,Î 0,1 ,R 0,1 ) has covariance functions: for t, t ′ ≥ 0,

). (Î 0,2 ,Ê 0 ,Î 0,2 ,R 0,2 ) has covariance functions: for t, t ′ ≥ 0,

Proof. By the independence of the sequences {λ 0 j } j≥1 and {λ 0,I k } k≥1 , it suffices to prove the joint convergence of Î N 0,1 ,Î N 0,1 ,R N 0,1 and Î N 0,2 ,Ê N 0 ,Î N 0,2 ,R N 0,2 separately. Recall the processes I N 0,1 and I N 0,2 defined in (3.18) and (3.19) , respectively. Similarly we define 

→ 0 in probability in D 7 , as N → ∞. The convergence for I N 0,1 −Î N 0,1 , I N 0,2 −Î N 0,2 → 0 in probability in D 2 is shown in (3.21) . For the other process, the convergence follows similar arguments as already used above. See also the proofs of Lemmas 6.1 and 8.1 in [17] . This completes the proof.

Remark 4.1. We remark that the limit (Î 0,1 ,R 0,1 ) can be written aŝ

where W 0,1 is a Gaussian white noise on R + satisfying

The limits (Ê 0 ,Î 0,2 ,R 0,2 ) can be written aŝ

where W 0,2 is a Gaussian white noise on R 2 + , independent of W 0,1 , such that

for 0 ≤ s ≤ t and 0 ≤ s ′ ≤ t ′ .

jointly with the convergence of Î N 1 ,Î N 2 ⇒ Î 1 ,Î 2 , where Î 1 ,Î 2 is given in Lemma 3.6. (Î 1 ,Î 2 ,Ê 1 ,Î 1 ,R 1 ) is a four-dimensional centered Gaussian process, independent of (Ê(0),Î(0)), (Î 0,1 ,Î 0,1 ,R 0,1 ) and (Î 0,2 ,Ê 0 ,Î 0,2 ,R 0,2 ), with the covariance functions: for t, t ′ ≥ 0,

follows from the same argument as in the proofs of Lemmas 8.3 and 8.4 in [17] . Let us sketch the argument, and in particular the coupling between Ê N 1 ,Î N 1 ,R N 1 and Î N 1 ,Î N 2 , referring the reader to [17] for the technical details. We have shown the joint convergence of Î N 1 ,Î N 2 in (3.53) using the processes I N 1 , I N 2 which are defined via the PRMQ(ds, dλ, du). Recall the definition of the pair (ζ, η) as a function of λ in (2.15 ). This defines a map φ : D → R 2 + such that (ζ, η) = φ(λ). Recall that ζ is the duration of the exposed period, and η the duration of the infectious period. We now define Λ : R + × D × R + → R 4 + by Λ(s, λ, u) = (s, φ (s) (λ), u), where φ (s) (λ) = (s + φ 1 (λ), s + φ 1 (λ) + φ 2 (λ)) .

The intuition behind this definition is that if s is the time of infection, φ (s) 1 stands for the end of the exposed period (transition from exposed to infectious), and φ (s) 2 the end of the infectious period (transition from infectious to recovered). We finally define a PRM Q 1 on R 4 + , which is the image ofQ by the mapping Λ, i.e., for any Borel subset A ⊂ R 4 + , Q 1 (A) =Q(Λ −1 (A)) .

Let Q 1 denote its associated compensated measure. We havê

and define the auxiliary procesŝ

It is not hard to show that the three processes M N A (t), L N 1 (t) and R N 1 (t) are martingales (with respect to three different filtrations), whose tightness is easy to establish. Moreover,

Hence the tightness of E N 1 and I N 1 . In order to establish the joint convergence I N 1 , I N 2 , E N 1 , I N 1 , R N 1 ⇒ Î 1 ,Î 2 ,Ê 1 ,Î 1 ,R 1 in D 5 as N → ∞, it remains to prove the joint convergence of the finitedimensional distributions, that is, for any k ≥ 1, 0 < t 1 < t 2 < · · · < t k , as N → ∞, ( I N 1 (t 1 ), I N 2 (t 1 ), E N 1 (t 1 ), I N 1 (t 1 ) converges to (Î 1 (t 1 ),Î 2 (t 1 ),Ê 1 (t 1 ),Î 1 (t 1 ),R 1 (t 1 )), . . . , (Î 1 (t k ),Î 2 (t k ),Ê 1 (t k ),Î 1 (t k ),R 1 (t k )) in distribution in R 5k . Noting that an integral w.r.t. the compensated measure Q 1 can be considered in fact as an integral w.r.t. the compensated measure Q, it is clear that the computation done in the proof of Lemma 3.7 can be extended to this situation, yielding the announced result. We just show one piece of this computation, namely, for ϑ, θ, ϑ ′ , θ ′ ∈ R and t, t ′ > 0, we get

→ 0 in D 3 in probability. We will show thatR N 1 − R N 1 → 0 and the other two follow from similar arguments. It is clear that

We have

Denote these two terms as R N 1 (t) and R N 2 (t). We show tightness of these two processes. Observe that both are increasing in t. Thus, we only need to verify the following (see the Corollary on page 83 in [4] ): for any ǫ > 0, and i = 1, 2,

By some direct calculations (see also the similar derivations in the proofs of Lemmas 6.3 and 8.3 in [17] for the models with constant infectivity rate), this condition for both processes R N 1 (t) and R N 2 (t) reduces to show that lim sup

as δ → 0. It is clear the expectation in (4.8) is bounded by δ 2 sup 0≤s≤T E |Υ N (s)| 2 , so the claim follows by (3.38) . For (4.9), we first observe that

Thus, we have where the last inequality follows from applying (3.8) (noting that the second integral requires using interchange of the order of integration first). Thus by (3.36), we obtain (4.9). This completes the proof.

Remark 4.2. We remark that the limit (Ê 1 ,Î 1 ,R 1 ) can be written aŝ

where W 1 is a Gaussian white noise on R 3 + , independent of the pair (W 0,1 , W 0,2 ), such that

Completing the proof of Theorem 2.4. The representations of E N (t), I N (t) and R N (t) in (4.1), (4.2) and (4.3), give a natural integral mapping from (Ê N (0),Î N (0)), Ê N 0 ,Î N 0,1 ,Î N 0,2 ,R N 0,1 ,R N 0,2 , Ê N 1 ,Î N 1 ,R N 1 , Ŝ N ,Ĵ N and the fixed functions G 0 (t), F 0,I (t), Ψ 0 (t),λ 0 (t),λ 0,I (t),λ(t),S(t) and I(t) to the processes Ê N ,Î N ,R N . The mapping is continuous in the Skorohod J 1 topology (by a slight modification of Lemma 8.1 in [17] ). We can then apply the continuous mapping theorem to conclude the convergence Ê N ,Î N ,R N ⇒ Ê ,Î,R in D 3 . Given their joint convergence with Î N 0,1 ,Î N 0,2 in Lemma 4.1 and (Î N 1 ,Î N 2 ) in Lemma 4.2, we can also conclude the joint convergence of the processes Ê N ,Î N ,R N with (Ŝ N ,Î N ).

Finally, for the covariance functions of (Ŵ S ,Ŵ I ,Ŵ E ,Ŵ I ,Ŵ R ), we have observed in the previous section from (3.9) thatŴ I =Î 0,1 (t)+Î 0,2 (t)+Î 1 (t)+Î 2 (t), and we also observe from (4.1), (4.2) and (4.3) thatŴ E =Ê 0 (t)+Ê 1 (t),Ŵ I =Î 0,1 (t)+Î 0,2 (t)+Î 1 (t),Ŵ R =R 0,1 (t)+R 0,2 (t)+R 1 (t). Then the covariance functions can be easily computed from the the covariance functions of (Î 0,1 ,Î 0,1 ,R 0,1 ), (Î 0,2 ,Ê 0 ,Î 0,2 ,R 0,2 ) and (Î 1 ,Î 2 ,Ê 1 ,Î 1 ,R 1 ) in Lemmas 4.1 and 4.2. This completes the proof of Theorem 2.4.

In this section, we discuss how the results can be generalized to the SIS and SIRS models. We state the FCLTs for these models without proofs, since they can be done analogously. 5.1. Generalized SIS models with varying infectivity. In the SIS model, individuals become susceptible immediately after going through the infectious periods. Since S N (t) = N − I N (t), the epidemic dynamics is determined by the two dimensional processes I N (t) and I N (t). As stated in Remark 2.3 of [9] , the FLLN limit (Ī(t),Ī(t)) is determined by the two-dimensional integral equations:Ī (t) =Ī(0)λ 0 (t) + 

Generalized SIRS models with varying infectivity. In the SIRS model, individuals experience the infectious and recovered/immune periods, and then become susceptible. We use I N (t) and R N (t) to denote the numbers of infectious and recovered/immune individuals, respectively, corresponding to the processes E N (t) and I N (t) in the SEIR model. Similarly, we use {λ 0 j } j≥1 and {λ i } i≥1 to denote the infectivity processes of the initially and newly infectious individuals, respectively, and also use {(ξ 0 j , η 0 j )} j≥1 and {(ξ i , η i )} i≥1 for the infectious and recovered/immune periods for the initially and newly infected individuals, respectively. Denote the remaining immune time of the initially recovered/immune individuals by η 0,R k (changing notation η 0,I to η 0,R accordingly). Of course, {λ 0 j } j≥1 and {λ i } i≥1 only take positive values over the intervals [0, ζ 0 j ) and [0, ζ i ), respectively. The definitions of the variables (ξ i , η i ), (ξ 0 j , η 0 j ) and η 0,R k in (2.15), (2.16) and (2.17) also need to change accordingly in a natural manner. The c.d.f.'s G 0 , G denote the distributions of infectious periods of the initially and newly infectious individuals, and the c.d.f.'s F 0,R and F denote the distributions of the recovered/immune periods of the initially and newly recovered individuals. Similarly for the notation Ψ 0 , Ψ, Φ 0 , Φ.

Since S N (t) = N − I N (t) − R N (t), the epidemic dynamics is described by the three processes (I N , I N , R N ). As stated in Remark 2.3 in [9] , the FLLN limit Ī ,Ī,R is determined by the three-dimensional integral equations: We next state the FCLT for these processes. (Ŵ I ,Ŵ I ,Ŵ R ) is a three-dimensional centered Gaussian process, independent ofÎ(0) andR(0), which has the covariance functions as the processes (Ŵ I ,Ŵ E ,Ŵ I ) in the SEIR model withS = 1 −Ī −R in the corresponding formulas. Ifλ 0 (t), G 0 and F 0,R are continuous, the limitsÎ,Î and R have continuous paths.

We recall Theorem 13.5 from [4] . In the cited reference, the result states conditions which imply convergence in D([0, 1]; R), while we want to obtain convergence in D := D([0, +∞); R). This implies that in our case condition (13.12) in [4] is not required, and the above quoted theorem becomes the following. Theorem 6.1. Let {X N , N ≥ 1} be a sequence of random processes with trajectories in D. X being a process with trajectories in D, suppose that the two following conditions are satisfied.

For any k ≥ 1, t 1 , t 2 , . . . , t k ≥ 0, (X N t 1 , X N t 2 , . . . , X N t k ) ⇒ (X t 1 , X t 2 , . . . , X t k ), (6.1)

and there exists a continuous nondecreasing function F from R + into itself, α > 1/2 and β ≥ 0 such that for any 0 ≤ r ≤ s ≤ t, N ≥ 1 and λ > 0,

Then X N ⇒ X in D.

As explained in [4] , the following moment condition is stronger than (6.2)

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

The duration of the closed stochastic epidemic

Convergence of probability measures

Age-of-infection and the final size relation

Stochastic epidemic models with inference

Probability and Stochastics

Epidemic models with varying infectivity

Temporal dynamics in viral shedding and transmissibility of COVID-19

An age dependent epidemic model

A mathematical model for chagas disease with infection-age-dependent infectivity

A contribution to the mathematical theory of epidemics

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

Multi-patch epidemic models with general exposed and infectious periods

Functional limit theorems for non-Markovian epidemic models

The asymptotic evolution of the general stochastic epidemic

On the asymptotic distribution of the size of a stochastic epidemic

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

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