key: cord-0131556-cs6ce6n5
authors: Pang, Guodong; Pardoux, Etienne
title: Multi-patch epidemic models with general exposed and infectious periods
date: 2020-06-25
journal: nan
DOI: nan
sha: aee9127021206b0d404b84c3473e87245e904e71
doc_id: 131556
cord_uid: cs6ce6n5

We study multi-patch epidemic models where individuals may migrate from one patch to another in either of the susceptible, exposed/latent, infectious and recovered states. We assume that infections occur both locally with a rate that depends on the patch as well as"from distance"from all the other patches. The exposed and infectious periods have general distributions, and are not affected by the possible migrations of the individuals. The migration processes in either of the three states are assumed to be Markovian, and independent of the exposed and infectious periods. We establish a functional law of large number (FLLN) and a function central limit theorem (FCLT) for the susceptible, exposed/latent, infectious and recovered processes. In the FLLN, the limit is determined by a set of Volterra integral equations. In the special case of deterministic exposed and infectious periods, the limit becomes a system of ODEs with delays. In the FCLT, the limit is given by a set of stochastic Volterra integral equations driven by a sum of independent Brownian motions and continuous Gaussian processes with an explicit covariance structure.

Multi-patch epidemic models have been used to study various infectious diseases, for example, nosocomial infection [24] , vector-borne diseases [23] , HIV/AIDS transmission [22] , SARS epidemic [26] , and so on. Patches refer to different locations, for example, a densely populated city and a less populated rural area, and thus such models capture geographic heterogeneity. It also helps to study the effect of migrations or lock-down measures among different population groups or locations. In the Covid-19 pandemic, it has been observed that the infectivity in different regions may vary and is impacted by various social-distance and lock-down measures [32, 34] .

Compartment ordinary differential equation (ODE) models are often used to study the dynamics of such multi-patch epidemic models. It is well known that the ODE dynamics arises from the Markovian assumptions in the stochastic multi-patch epidemic model, that is, the infection process is Poisson, the infectious (and/or exposed/latent) periods are exponentially distributed and the migration processes are also Markovian [2, 10, 25, 26, 5, 23] .

In this paper, we study multi-patch SEIR models, in which each individual may experience successively Susceptible (S), Exposed (E), Infectious (I) and Recovered (R) periods, and the exposed/latent and infectious periods have a general joint distribution (possibly correlated), while the migration processes are Markovian. SEIR models are widely used to study the propagation of various epidemics in that they capture both exposed/latent and infectious stages, see, e.g., [2, 10] . The infection in the multi-patch model is assumed to be both local, and from distance. That is, the infection rate in a given patch depends on the susceptible population in that patch, and on the infectious population in all the patches. Individuals may migrate from one patch to another in each of the Susceptible, Exposed (Latent), Infected and Recovered stages. The reason for infection at distance is twofold. First, if there were no migration among a subset of the patches (i.e., the migration rates among them are zero), we could consider some of the patches as substructures of the population, like age classes, which infect each other. Second, some of the movements of the population should not be considered as migrations, but visits from one patch to another, during a week-end or holiday, with a return at home at the end of a short period. Such movements may produce infection of a susceptible individual in patch i by an infectious individual from patch i ′ = i. Those displacements would be very complicated to model as such. We think that infection at distance is a reasonable way to model infections due to such movements. Mathematically, this requires a much more sophisticated model for the instantaneous infection rate, see (2.2) . The formulation is new in the literature, even in the Markovian setting. Moreover, it brings new challenges in the proofs of the limit theorems (see Lemmas 4.10-4.12 and Lemmas 5. 3-5.4 ), compared to those of the one-patch non-Markovian models in [31] .

We describe the evolution dynamics by tracking the time epochs of becoming exposed and/or infectious and the location of an individual at these event times. Specifically, in the multi-patch SEIR model, each individual tracks the time epochs of becoming exposed, infected and recovered, and is associated with two Markov chains that are used to track their movement starting when the individual becomes exposed and infectious, respectively. For the initially exposed and/or infected individuals, we also assume that their remaining exposing and/or infectious periods have general distributions, which may be different from those of the newly exposed/infected individuals. For these initially exposed/infectious individuals, we also track their movement among the patches using Markov chains while being exposed/infectious. Although the idea of tracking time epochs of each individual is analogous to the one-patch SIR/SEIR models in [31] , the migrations among different patches makes the system dynamics much more challenging to describe (see the equations (2.5)-(2.8)) and analyze, despite the independence of the Markovian migration processes from the exposed/infectious durations. The formulation and proofs constitute non-trivial generalizations of the one-patch SEIR model.

Given the representations with these time epochs and location processes, we show a functional law of large numbers (FLLN) and a functional central limit theorem (FCLT). More precisely, we divide by N (which is the total population size) the equations for the evolution of the numbers of individuals in each patch and compartment, thus obtaining equations for the proportions of individuals in each patch and compartment. Those proportions are shown to converge in probability, locally uniformly in t, towards the solution of a system of integral equations. This is our FLLN. Next, if we multiply by √ N the difference between the proportions in the N stochastic model and the limiting proportion, that renormalized difference is shown to converge in distribution to the solution of a set of linear integral equations driven by Gaussian noise. This is the FCLT. Needless to say, the deterministic system of equations obtained for the FLLN limit is much simpler than the N stochastic model. It can be used for predictions of an epidemic which affects a population of reasonably large size, at least away form the very early and very final stages, when the number of infected individuals is not of the order of N . Note that the FCLT tells us that the error made on the proportions by using the deterministic model is of the order of N −1/2 .

The FLLN limits are determined by a set of Volterra integral equations. When the infectious (and exposed/latent) periods are deterministic, we can write the fluid integral equations as a set of ODEs with delay (Remark 3.1). The limit processes in the FCLT are determined by a set of stochastic Volterra integral equations, driven by a sum of independent Brownian motions and continuous Gaussian processes with a certain covariance structure. When the infectious (and exposed/latent) periods are deterministic, the limits become stochastic differential equations with linear drifts and delay (Remark 3.3). In both FLLN and FCLT limits, the effects of migrations are exhibited through the transition probability functions and transition rates of the migration Markov processes. We discuss how the results simplify in the SIR model as a special case, and also how the approach and results can be extended to study multi-patch SIS and SIRS models (Section 3.3).

In the proofs of these results, we employ Poisson random measures (PRMs) that are constructed as the sums of the Dirac masses at the time epochs of becoming exposed and infectious, the infectious and exposing periods and the Markov chain starting from the location of each individual at those epochs. While to the Markov migration process have naturally associated martingales, which are easily proved to be tight in the appropriate path space, the non-Markovian nature of the epidemic process does not produce obvious martingales. However, as in our previous work [31] , we are able to find various martingales attached to various and non-standard filtrations, which are constructed from our representation of the epidemic by integrals with respect to PRMs. We use the martingale properties and convergence theorems as critical tools in the proofs. For the single-patch SIR and SEIR models with general infectious and exposing periods, an approach using PRMs that are constructed at the time epochs of becoming infectious (and/or exposed), was developed in [31] . The approach is further developed in this paper for multi-patch SIR and SEIR models, to track the locations of each individual at each event epoch. Incorporating infection from distance in addition to local infections in the model also brings in new technical challenges in the proofs of both the FLLN and FCLT.

This paper contributes to the limited literature on stochastic epidemic models with general exposed/infectious periods. We refer the readers to the overview in Chapter 3.4 of [10] on the common approaches to study non-Markovian epidemic models and the limit theorems for the final sizes of the epidemic; see also the recent method using piecewise Markov deterministic processes in [12] and [18] for the SIR model. FLLNs and FCLTs are proved for some age and density dependent population models in [36, 37, 38] , which includes the SIR model with the infection rate depending on the number of infectious individuals, and general infectious period as a special case. Reinert [33] proves a FLLN for the empirical measure of the SIR epidemic dynamics using Stein's method, while no FCLT has been proved with that approach. In [31] , both FLLN and FCLT were established for the epidemic dynamics in the classical models (SIR, SIS, SEIR, SIRS) where the PRM representations of the dynamics plays a fundamental role in the proofs. The FCLT limit for the SIR model in [31] is similar to that in [36, 37, 38] , however, the proof approaches are completely different; in addition, the distribution function of the infectious periods is assumed to be continuously differentiable in [36, 37, 38] while no condition is imposed in [31] . We highlight that the distribution functions of the exposed and infectious periods in this paper are general without requiring any conditions. The integral equations for the SEIR model in [31] are also used to estimate the state of the Covid-19 epidemic in [17] . For SIR and SEIR models with varying infectivity, where each individual is associated with an i.i.d. random infectivity, which is a function of the time elapsed since infection, FLLN and FCLT have recently been established in [16, 29] . Although Volterra integral equations were used to describe the proportion of infectious population in the SIS, SIR or SEIR model without proving an FLLN (see [8, 9, 13, 14, 21, 35] ), as far as we know no Volterra integral equations have been proposed so far for multi-patch epidemic models with general infectious (and/or exposed) periods. Our work shows both FLLN and FCLT for non-Markovian multi-patch models, and identify (stochastic) multidimensional Volterra integral equations as their limits.

It is also worth mentioning the multi-type epidemic models where the population splits up into multiple groups of individuals and each group may infect any other group in addition to itself (no migration), see Chapters 6.1 and 6.2 in [2] and [3, 4] . The special case of our model with zero migration rates covers that situation. In those models, proportionate mixing taking into account control measures like social distance or lockdowns can also be incorporated.

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}. For x ∈ R, let x + = x ∨ 0 and x − = −(x ∧ 0). Let D = D([0, T ], R) denote the space of R-valued càdlàg functions defined on [0, T ]. Throughout the paper, convergence in D means convergence in the Skorohod J 1 topology, see chapter 3 of [6] . Also, D k stands for the k-fold product equipped with the product topology. Let C be the subset of D consisting of continuous functions. Let C 1 consist of differentiable functions whose derivative is continuous. For any function x ∈ D, we use the notation x T = sup t∈[0,T ] |x(t)|. For two functions x, y ∈ D, we use x•y(t) = x(y(t)) denote their composition. All random variables and processes are defined on a common complete probability space (Ω, F, P). The notation ⇒ means convergence in distribution. We use 1 {·} for indicator function, and occasionally we shall write 1{.} in case the first notation is not readable enough. We use small-o notation for real-valued functions f and non-zero g: f (x) = o(g(x)) if lim sup x→∞ |f (x)/g(x)| = 0. We useî to denote the unit imaginary number. 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 consider a multi-patch epidemic model, where individuals in each patch experience the Susceptible-Exposed (Latent)-Infectious-Recovered (SEIR) process. The patches may refer to populations in different locations, for example, a densely populated city and a less populated rural area. As explained in the introduction, susceptible individuals in each patch are infected both locally, by infectious individuals located in the same patch, and at distance, by infectious individuals from other patches. The rate of infection is different in each patch (because of the differences in the density of population or in the type of available public transportations), while the law of the infectious period is the same (due to the same illness).

Let N be the total population size and L be the number of patches. The set of patches will be denoted L = {1, . . . , L}. (We use indices i, i ′ , ℓ, ℓ ′ for elements in L, and occasionally i ′′ , ℓ ′′ .)

count the numbers of individuals that are susceptible, exposed (latent), infectious and recovered in patch i at time t, respectively. We have the balance equation:

It is straightforward to allow R N i (0) to be nonzero, however the initially recovered individuals can be removed from the population under consideration from the beginning, since they do not interact with the propagation of the epidemic.

2.1. The infection process. Let λ i be the infection rate of patch i ∈ L. Define the following processes, for some 0 ≤ γ ≤ 1,

where κ ii = 1 and κ iℓ ≥ 0 for i = ℓ represent the infectivity from distance. Letκ i := L ℓ=1 κ iℓ and κ := max i∈Lκi . The rate of new infections in patch i at time t is λ i Υ N i (t). Let us explain the role of the parameter γ.

In the homogeneous population model, where L = 1, (2.1) tells us that in the unique patch, S N (t) + E N (t) + I N (t) + R N (t) = N , hence Υ N (t) is the same, irrespective of the value of γ. The rationale of this form of the infection rate is as follows. Each infectious individual meets others at rate β. Since we assume that the individual who is met is chosen uniformly at random in the whole population, he/she is susceptible with probability S N (t)/N . In that case, the encounter results in a new infection with probability p. If we let λ = β × p, we find the above formula λ i Υ N i (t) for the rate of new infections in case L = 1. Now, consider the case L > 1. We do not factorize λ into β × p anymore, or equivalently do as if p = 1. In order to make the role of the parameter γ transparent, let us define

, the total population in the patch i, and rewrite

In the case γ = 1, the rate of encounters of individuals in patch i by a given infectious is given as λ i for an infectious of the same patch, and equal to λ i κ ii ′ for an infectious from patch i ′ , whatever the total population in patch i at time t may be. In case γ = 0, the same rate is proportional to B N i (t), the total population of patch i at time t. In the intermediate cases, the rate lies between those two extremes. The case γ = 1 seems to be used in most spatial epidemics models. The values of λ i 's can correct for the different densities of population of the various patches, resulting in more or less encounters. Indeed, we believe that the rate of encounters by any individual is very different in a densely populated area, from what it is in a desert. However, especially in the stochastic model, the population size in each patch may fluctuate significantly, which we believe is a good motivation for using a model with γ < 1.

We shall prove the FLLN for any value of γ ∈ [0, 1], and the FCLT only for γ ∈ [0, 1) in the general case, and for all γ ∈ [0, 1] in the case that infections are only local, i.e., κ iℓ = 0 for i = ℓ. The reason for this restriction is that in the case γ = 1 and ℓ =i κ iℓ > 0, we are not able to establish the estimate (5.18) in Lemma 5.3 below.

By convention, we shall assume that the same holds in case γ = 1, i.e., 0

Let A N i (t) be the cumulative counting process of individuals in patch i that get infected on the time interval (0, t]. Then we can give a representation of the process A N i (t) via the standard Poisson random measure (PRM) Q i on R 2 + (with mean measure dsda), the various {Q i , i ∈ L} being mutually independent,

where P A,i is a unit-rate Poisson process, and independent from each other for i ∈ L. But the first description will be more useful for us. We let {τ N j,i , j ≥ 1} denote the successive jump times of the process A N i , for i ∈ L. (Note that all the analysis and results can be easily extended to a deterministic time-dependent rate function λ i (t). For example, in the expression above, we have an integral t 0 λ i (s)Υ N i (s)ds instead.) 2.2. On the exposed and infectious periods. The E N i (0) initially exposed individuals experience the exposed and infectious periods before recovery. Let {η 0 k,i : k = 1, . . . , E N i (0)} be the remaining exposed periods of the initially exposed individuals in patch i. After the exposed period, let {ζ −k,i : k = 1, . . . , E N i (0)} be the durations of their infectious periods. The I N i (0) initially infectious individuals experience a remaining infectious period before recovery, and let ζ 0 k,i , k = 1, . . . , I N i (0), denote their remaining infectious periods. The A N i (t) newly infected individuals in patch i experience the exposed and infectious periods. Let {η j,i : j ∈ N} and {ζ j,i : j ∈ N} be the associated exposing and infectious periods. Assume that {ζ 0 k,i }, {(η 0 k,i , ζ −k,i )} and {(η j,i , ζ j,i )} are all i.i.d. sequences of random variables having distribution functions F 0 , H 0 (du, dv) and H(du, dv), respectively, and they are also mutually independent. Note that ζ j,i is defined for j ∈ Z and i ∈ L (those with j < 0 code the infectious periods of the initially exposed individuals, while those with j > 0 code the infectious periods of the newly exposed individuals). Let G 0 and F be the marginals of H 0 for η 0 k,i and ζ −k,i , and G and F be the marginals of H for η k,i and ζ k,i , respectively. (It is reasonable to assume that the marginal distributions of ζ −k,i and ζ k,i are the same.) Also let F 0 (·|u) and F (·|u) be the conditional c.d.f.'s of ζ −k,i and ζ k,i , given that η 0 k,i = u and η k,i = u, respectively. Let

The migration processes. Individuals may migrate from patch ℓ to ℓ ′ in any of the four epidemic stages, with rates ν S,ℓ,ℓ ′ , ν E,ℓ,ℓ ′ , ν I,ℓ,ℓ ′ and ν R,ℓ,ℓ ′ for the susceptible, exposed, infectious and recovered ones, respectively. For each individual, the times between migrations in each of the stages are exponentially distributed.

In order to track the location/patch of the j-th individual who got exposed in patch ℓ at time τ N j,ℓ , we first use the Markov process X j ℓ , taking values in L, associated to the rates ν E,ℓ,ℓ ′ for ℓ ′ ∈ L. It takes effect from the time τ N j,ℓ of becoming exposed, until the time τ N j,ℓ + η j,ℓ when this individual becomes infectious. Given that this individual has migrated to patch ℓ ′ at the end of the exposed period (she/he may have done several migrations to other patches during the exposed period), that is X j ℓ (η j,ℓ ) = ℓ ′ , we then use another Markov process Y j,ℓ ℓ ′ to track the location/patch of the individual during the infectious period ζ j,ℓ , starting from ℓ ′ and associated to the rates ν I,ℓ ′ ,ℓ ′′ for ℓ ′′ ∈ L. This process Y j,ℓ ℓ ′ only takes effect from the time of becoming infectious τ N j,ℓ + η j,ℓ , until the time of recovery τ N j,ℓ + η j,ℓ + ζ j,ℓ . Suppose that the individual has migrated to patch i at the end of the infectious period, that is, Y j,ℓ ℓ ′ (ζ j,ℓ ) = i. The individual will then belong to the compartment of recovered individuals, and will migrate among patches according to the rates ν R,i,i ′ for i ′ ∈ L. Similarly we use X 0,k ℓ and Y −k,ℓ ℓ ′ for the initially exposed individuals k = 1, . . . , E N ℓ (0) that have been exposed at time 0 in patch ℓ. They are again associated with the rates ν E,ℓ,ℓ ′ for ℓ ′ ∈ L, and ν I,ℓ ′ ,ℓ ′′ for ℓ ′′ ∈ L, respectively. In addition, we also use Y 0,k ℓ for the initially infectious individuals k = 1, . . . , I N ℓ (0) that have been infectious at time 0 in patch ℓ. It is again associated with the rates ν I,ℓ,ℓ ′ for ℓ ′ ∈ L. We assume that for each j, X j ℓ and Y j,ℓ ℓ ′ are independent for ℓ, ℓ ′ ∈ L, and for each k, X 0,k ℓ and Y −k,ℓ ℓ ′ are independent for ℓ, ℓ ′ ∈ L. We also assume that all these Markov

for ℓ, ℓ ′ , ℓ ′′ ∈ L, j ≥ 1 and t ≥ 0. For each ℓ, the processes {X 0,k ℓ } k have the same transition function (p ℓ,ℓ ′ (·)) ℓ ′ ∈L as {X j ℓ } j , and the process {Y −k,ℓ ℓ ′ } k has the same transition function (q ℓ ′ ,ℓ ′′ (·)) ℓ ′ ,ℓ ′′ ∈L as {Y j,ℓ ℓ ′ } j , and the process {Y 0,k ℓ } k also has the transition function (q ℓ,ℓ ′ (·)) ℓ,ℓ ′ ∈L .

2.4. Epidemic evolution dynamics. The multi-patch SEIR epidemic evolution dynamics can be described as follows:

where P S,i,ℓ , P R,i,ℓ , i, ℓ ∈ L, are mutually independent unit-rate Poisson processes, which are globally independent of the Q i 's. The dynamics of S N i (t) is straightforward since it is simply equal to the number of initially susceptible individuals minus the number of exposed ones in patch i and then take into account the migrations. For the dynamics of E N i (t), the first term represents the number of initially exposed individuals from patch ℓ that remain exposed and are in patch i at time t, and the second term represents the number of newly exposed individuals from patch ℓ that remain exposed and are in patch i at time t. In the expression for I N i (t), the first term counts the number of initially infectious individuals from all the patches that remain infectious and are in patch i at time t, and the second term counts the numbers of initially exposed individuals from all the patches that have become infectious and are in patch i at time t (for tracking purposes, the location at the epochs of becoming infectious is recorded). Also note that we use the Markov process Y 0,k ℓ to indicate that these are for the initially exposed individuals. The third term counts the number of newly exposed individuals at all patches that have become infectious and are in patch i at time t, and we also track the patch in which each individual has become infectious. In the expression for R N i (t), the first term represents the number of initially infectious individuals from patch ℓ that have recovered by time t and were in patch i at the time of recovery, the second term represents the number of initially exposed individuals from patch ℓ that have recovered by time t, and were in patch i at the time of recovery, while becoming infectious in patch ℓ ′ , the third term represents the number of newly exposed individuals from patch ℓ that have recovered by time t, and were in patch i at the time of recovery while becoming infectious in patch ℓ ′ .

It is not easy to take the limit as N → ∞ in the formulas of E N i and I N i above. We now derive the following representations, which will be very helpful in the proofs of our results.

where P E,i,ℓ and P I,i,ℓ , i, ℓ ∈ L, are all unit-rate Poisson processes, mutually independent, and also independent of P A,i , P S,i,ℓ and P R,i,ℓ .

Before turning to the proof, let us comment on these formulas. In the expression of E N i (t), the first and third term count the number of initially exposed individuals in patch i, and the number of those whose became exposed in patch i on the time interval [0, t]. The second and fourth terms subtract the numbers of initially and newly exposed individuals in any patch who have become infectious before time t in patch i. Finally the last two term count the numbers of migrations of exposed individuals to and from patch i. The expression of I N i (t) is similar, except that the second term subtracts the number of initially infectious individuals in any patch who have recovered before time t in patch i, and the fourth and sixth terms subtract the numbers of initially and newly exposed individuals in any patch who have recovered before time t in patch i, where the patch in which they became infectious is also tracked.

Proof. In the representation of E N i (t), we observe that

and for ℓ = i,

k,i )=ℓ is the number of initially exposed individuals from patch i that are in patch ℓ at the time t ∧ η 0 k,i for k = 1, . . . , E N i (0). We also observe that

and for ℓ = i,

j,i +η j,i )=ℓ denotes the number of individuals who became exposed at time τ N j,ℓ ∈ (0, t) in patch i, and are in patch ℓ at time t ∧ (τ N j,i + η j,i ) for j = 1, . . . , A N i (t).

It is clear that

Thus, using the above identities, we obtain the expression in (2.9). A similar argument gives the expression in (2.10).

2.5. Using PRMs in order to represent some terms of the model. We end this section of presentation of our model with the description of the representations of some of the key components in the dynamics of the above model via PRMs. Those will play an important role in the proofs below. The infection process A N ℓ has the representation (2.3), which makes use of the PRM Q ℓ . Define a PRMQ ℓ (ds, da, du, dθ) on R 3 + × L, which is the sum of the Dirac masses at the points

, and again an infection occurs at time τ N j,ℓ if and only if A N j,ℓ ≤ λ ℓ Υ N ((τ N j,ℓ ) − ). We can then write for ℓ, i ∈ L,

The following FLLN is proved in Section 4.

Theorem 3.1. Under Assumption 3.1, and assume that F 0 and G 0 are continuous,

is the unique solution of the following system of deterministic integral equations:

and

Note that if the exposed and infectious periods are independent for each individual, we have

and

Remark 3.1. Suppose the exposed and infectious periods are deterministic, taking values of t e > 0 and t o > 0. Also, assume that the remaining exposed and infectious periods of the initially exposed and infectious are uniformly distributed over (0, t e ) and (0, t o ), respectively. These are the corresponding equilibrium distributions of the deterministic ones. Recall that for any c.d.

It is easy to see that we obtain a set of ODEs with delay after taking derivative.

be the diffusionscaled process whereZ N := N −1 Z N andZ is its limit.

where the limit is the unique solution of the following system of stochastic Volterra integral equations driven by continuous Gaussian processes: 

are continuous Gaussian processes, independent of the Brownian motions above, with mean zero and covariance functions:

, and Ê ℓ,i (t),Î ℓ,i (t) are independent from each other, and

Remark 3.2. The continuity of F 0 and G 0 will be important in our proofs of the FLLN and FCLT, in particular for the convergence of the processes associated with the initially exposed and infectious individuals in Lemmas 4.4 and 5.2. Note that this assumption is not really restrictive, in the sense that even when F or G is a Dirac measure (deterministic duration), the time in the past when the initially exposed (or infectious) individuals have been infected (or have become infectious) would most naturally be assumed to follow a uniform distribution on some interval dictated by F or G as discussed in Remark 3.1. 

On the multi-patch SIR, SIS and SIRS models.

Multi-patch SIR model. The multi-patch SEIR model includes the multi-patch SIR model as a special case, without the exposed periods and the associated Markov chain X. The infectious process I N i (t) becomes

which can be also expressed as

In the FCLT, we obtain

i,ℓ being mutually independent standard Brownian motions, and with the deterministic functionsS i ,Ī i ,R i given above. The processesÎ 0 ℓ,i andÎ ℓ,i are continuous Gaussian processes with mean zero and covariance functions:

otherwise.

In addition,Î 0 ℓ,i andÎ ℓ,i are independent, and also independent of the Brownian terms. 3.3.2. Multi-patch SIS model. The analysis of the multi-patch SIR model 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

Thus, in the FLLN, we obtain the same limitĪ i in (3.20) 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.24) for the multi-patch SIR model, and the limitŜ i (t):

Multi-patch SIRS model. The analysis for the multi-patch SEIR model can be easily extended to multi-patch SIRS model, where in each patch, the population is grouped into susceptible, infectious, and recovered individuals and individuals become susceptible after experiencing a recovery period. In this model, the infectious and recovered processes I N i , R N i correspond to the exposed and infectious processes E N i , I N i in the SEIR model. In the description of the epidemic dynamics, we need to change the dynamics of S N i in (2.5) by adding the individuals that have become susceptible after recovery, i.e., the first three terms in R N i in (2.8) . This is similar to the susceptible process S N i in (3.27) for the SIS model. Then it is straightforward to write down the limit processes in the FLLN and FCLT for the processes (

In this section we prove Theorem 3.1. Here we extend the approach in [31] . Specifically, we first establish that the process {(Ā N 1 , . . . ,Ā N L ) : N ≥ 1} is tight and then along its convergent subsequence with a given limit, we prove the convergence of

We then identify the limit of {(Ā N 1 , . . . ,Ā N L )}, which allows us to show that the above limits satisfy the system of equations (3.2)-(3.6). Finally we show that this system of equations (3.2)-(3.6) has at most one solution, which implies that the whole sequence converges. Due to the complications from the migration processes, the proofs of tightness for some key component processes become much more involved, see Lemmas 4.4 and 4.8, and in addition, the formula of Υ N i (t) for the infection also brings some new challenges, see Lemma 4.11. Let us first rewrite the representation (2.3) of the processes A N i (t) which uses the PRM Q i (ds, da).

It has the predictable quadratic variation: 

and

it follows readily from Doob's inequality that, as N → ∞, M N A,i → 0 in probability, locally uniformly in t . In the following we consider a convergent subsequence of {(Ā N 1 , . . . ,Ā N L )}. Before we establish the convergence of the S N i 's, let us establish a simple technical Lemma, which will be useful below. Lemma 4.2. Let A be a d × d matrix. For any F ∈ D(R + ; R d ), the equation

has a unique solution X ∈ D(R + ; R d ), and the mapping A defined by U := A(F ) is continuous from D(R + ; R d ) equipped with the Skorohod J 1 topology, into itself.

Proof. Let V t := U t − F t . It is easy to verify that U solves (4.6) iff V solves the linear ODE

This ODE has a unique solution, which is given by a well-known explicit formula, from which we deduce that

The continuity of the mapping A is now clear. 

Proof. Starting from (2.5), we can rewrite the processesS N i as

The processesM N S,ℓ,i are square integrable martingales with respect to the filtration {F N S (t) : t ≥ 0}, defined by

They have the predictable quadratic variation: For ℓ, i ∈ L, let

We first treat the components associated with the initial quantities. 

in probability, uniformly in t, where for ℓ, i ∈ L and t ≥ 0, We first focus onĨ N,0,1 ℓ,i (t). Note that, since ζ 0 k,ℓ and Y 0,k ℓ are independent,

where the expectation is taken under the condition that Y 0,k ℓ (0) = ℓ. Note that the pairs (ζ 0 k,ℓ , Y 0,k ℓ (·)) are independent over k, and have the same distributions. Thus, by the LLN of i.i.d. random variables, we obtain that for each t ≥ 0, as N → ∞,

ℓ,i (t) in probability. In order to establish locally uniform convergence in t, is suffices to establish tightness in D, which (see the Corollary of Theorem 7.4 in [6] ) will follow from the fact that lim sup

(4.11)

By the independence of the pairs

The first term is bounded by

as N → ∞, so the right hand side converges as N → ∞ to 1 t+δ t q ℓ,i (u)F 0 (du) > ǫ/2Ī ℓ (0) which vanishes for δ > 0 small enough, uniformly w.r.t. t ∈ [0, T ]. Hence, (4.11) follows. The convergence of the other processesĨ N,0,1 ℓ,i follows similarly.

We next sketch the proof for the convergence ofĨ N,0,2 ℓ,i since it follows similar steps as that of

which implies by the LLN of i.i.d. variables that for each t ≥ 0,Ī N,0,2 ℓ,i (t) ⇒Ī 0,2 ℓ,i (t) as N → ∞. The convergence of finite dimensional distribution is a straightforward extension. For tightness we use the same approach as above. We start with

We next note that each of the two terms on the right hand side is increasing in s, so that

The two probability terms on the right hand side are bounded by

respectively. Thus, as N → ∞, the right hand side converges to the sum of the two indicator terms. And the two terms in the limit both vanish for small enough δ. Thus, for any ǫ > 0,

Therefore, we can conclude thatĨ N,0,2 ℓ,i →Ī 0,2 ℓ,i in probability, uniformly in t, as N → ∞. Then, since G 0 is continuous, we can verify the continuity of the covariance function, and thus the continuity of the limit processesĪ 0,2 ℓ,i . This completes the proof. In the next proof, we will make use of the following result, which is Lemma 4.4 in [16] . In the next statement, D ↑ (R + ) (resp. C ↑ (R + )) denotes the set of real-valued nondecreasing function on R + , which belong to D(R + ) (resp. C(R + )). Lemma 4.5. Let f ∈ D(R + ) and {g N } N ≥1 be a sequence of elements of D ↑ (R + ) which is such that g N → g locally uniformly as N → ∞, where g ∈ C ↑ (R + ). Then, for any t > 0, as N → ∞,

We next prove the convergence ofĒ N ℓ,i andĪ N ℓ,i for ℓ, i ∈ L. Define the auxiliary processes: for

, t ≥ 0. We first prove these processes converge to the desired limits in the following lemma, and then show that these processes are asymptotically equivalent toĒ N ℓ,i andĪ N ℓ,i , ℓ, i ∈ L. 

in probability, where for ℓ, i ∈ L and t ≥ 0,

with Φ ℓ,i defined in (3.8).

Proof. Observe that for ℓ, i ∈ L, 

and P sup

Proof. We first considerĒ N ℓ,i . We havē

Then it is clear that for

. And by the independence of the pairs η j,ℓ , X j ℓ (·) and η j ′ ,ℓ , X j ′ ℓ (·) , we have

which, together with the upper bound in (4.3) for E Ā N ℓ (t) , implies that for any t > 0 and ǫ > 0, Observe that the three terms on the right hand side are nondecreasing in u. Thus we obtain P sup

Using the PRMQ ℓ (ds, da, du, dθ) and its compensated PRM Q ℓ (ds, da, du, dθ), we have

The first term converges to zero as N → ∞, and the second term satisfies

The second term on the right hand side of (4.21) can be treated similarly as the second term right above. Now for the third term, using (4.1), 

Therefore, we have proved (4.17) We next show (4.18), which follows from similar steps as above. We highlight the main differences below.

which implies that for any ǫ > 0,

Next, for t, s > 0, we have

Observe that each of the five terms on the right hand side is increasing in s. Thus, we have P sup

The last term is treated in the same way as the last term in (4.20) using the bound in (4.23). For the first two terms, we use the PRM representation in (2.12). We have

The first term is equal to

The second term is bounded by

Without the constant 2(λ ℓκℓ ) 2 , it satisfies

where the last step follows from the same argument as in (4.22) .

Similarly, for the second term on the right hand side of (4.24), we have

Here it is clear that the first term converges to zero as N → ∞, and the second term without the constant 2(λ ℓκℓ ) 2 satisfies 

The first terms on the right hand sides of both converge to zero as N → ∞, while the second terms are bounded by a constant times δ 2 . Thus we have shown that

Therefore, we have shown that (4.18) holds.

By the above two lemmas we have shown the following result. 

whereĒ ℓ,i andĪ ℓ,i are given in (4.13) and (4.14), respectively.

We are now ready to prove the convergence of Ē 

where the limits are the unique solution to the systems of ODEs: for i ∈ L,

withĒ 0 ℓ,i ,Ī 0,1 ℓ,i andĪ 0,2 ℓ,i being given in (4.9) and (4.10), and withĒ ℓ,i andĪ ℓ,i being defined in (4.13) and (4.14), respectively.

Proof. The proof for the process E N i (t) is similar to that of I N i (t), so we focus on I N i (t). By the representations of I N i (t) in (2.10), we havē 

in probability, locally uniformly in t. Hence, it follows from We next prove the convergence of (R N 1 , . . . ,R N L ). Similar to (2.10), we obtain the following representations for the process R N i (t):

Thus, we can represent the processesR N i (t) bȳ

Arguments very similar to those in the above proof allow us to conclude.

From the above arguments, since we have the joint convergence in Lemmas 4.4 and 4.8, we can conclude the joint convergence of (S N i ,Ē N i ,Ī N i ,R N i , i ∈ L). However, we have not yet quite explicited the limiting equations, since we have not expressedĀ i (t) in terms of (S i (t),Ē i (t),Ī i (t),R i (t), I ℓ (t), ℓ = i). It seems easy to do that since we know that

However, there is a difficulty in the case where both γ = 1 and ℓ =i κ iℓ = 0. Define the function ψ(s, e, i, r, u) = s(i+u) (s+e+i+r) γ on [0, 1] 4 ×[0,κ]. If either 0 ≤ γ < 1 or sup i ℓ =i κ iℓ = 0, ψ is continuous. However, if both γ = 1 and sup i j =i κ ij > 0, then ψ is not continuous at any point of the form (0, 0, 0, 0, u), with u > 0. Hence, if we want to include that case in our model, we need to prove that for any i ∈ L and T > 0, inf 0≤t≤T (S i (t) +Ē i (t) +Ī i (t) +R i (t)) > 0. Fortunately, we can prove such an estimate, although we do not have yet established the exact system of equations of the (S i (t),Ē i (t),Ī i (t),R i (t)), i ∈ L.

. For any i ∈ L and T > 0, there exists a constant C i,T > 0 which is such that for 0 ≤ t ≤ T ,

) is a solution of (4.7), (4.27), (4.28) and (4.29). Then we haveŪ

Differentiating equation (4.37) and exploiting the inequality W i (t) ≥ 0 for all t ≥ 0, we deduce that

for any 0 ≤ t ≤ T .

We can now explicit the processesĀ i (t).

. Then for any i ∈ L and 0 ≤ t ≤ T ,

Proof. Let D 4 denote the set of càdlàg functions from R + into [0, 1] 4 × 0,κ . For any γ ∈ [0, 1], the function ψ(s, e, i, r, u) = s(i+u) (s+e+i+r) γ is continuous for the Skorohod topology, on the subset of D 4 which is such that for any T > 0, inf 0≤t≤T {s(t) + e(t) + i(t) + r(t)} > 0. Thus, we deduce from the joint convergence

and Lemma 4.10 that

Consequently, given Lemma 4.1,

Therefore, any limit satisfies the system of integral equations given in Theorem 3.1.

The uniqueness of solutions to the set of integral equations in Theorem 3.1 follows from the next Lemma, from which we deduce that the whole sequence converges, and since the limit is deterministic, the convergence is in probability. This completes the proof of Theorem 3.1. Proof. In view of Lemma 4.10, if we take the difference between two solutions, any convex combination of those two solutions satisfies the lower bound (4.36). Sincē

and at each time t ∈ [0, T ], the derivatives of ψ with respect to each of its variables is bounded in absolute value by the supremum of 1 andκ iŪ −γ i (t) ≤κ i C −γ i,T , we can now apply a standard argument based upon Gronwall's Lemma in order to deduce uniqueness.

In this section, we prove Theorem 3.2. We thus generalize the approach in [31] to the multipatch model, by employing the standard technique of convergence of finite dimensional distributions and tightness as exposed in [6] . The migration processes require subtle care in proving the finite dimensional distribution convergence as shown in Lemmas 5.2 and 5.6 for the key components. The tightness proofs require a moment estimate for the supremum of the processes, which is challenging due to the formula Υ N (t) for infection, see Lemmas 5.3 and 5.4.

We first give the following representations of the diffusion-scaled processes. The process N i (t) can be decomposed as:

For the processŜ N i (t), we havê

For the processÊ N i (t), by the representation in (2.9), using the definitions of E N,0 ℓ,i (t) and E N ℓ,i (t) in (4.8), we obtain

where for ℓ, i ∈ L,

For the processÎ N i (t), we obtain

whereÊ N,0 ℓ,i (t) andÊ N ℓ,i (t) are defined in (5.6) and (5.7), respectively, and For the processR N i (t), we havê

where for ℓ, i ∈ L, and ℓ = i,

(5.14)

Let us consider the convergence ofĨ N,0,1 1,1 . Observe that the pairs ζ k,1 , Y 0,k 1 (·) and ζ k ′ ,1 , Y 0,k ′ 1 (·) are independent and have the same law. Thus, its proof follows in a similar approach for empirical processes, see, e.g., Theorem 14.3 in [6] . There are some differences due to the Markov process Y 0,k 1 , which we highlight below. So, we apply Theorem 13.5 in [6] . For each t > 0 and α ∈ R, we have

Similarly, it can be also shown that for any 0 < s < t < r and α 1 , α 2 ∈ R,

Thus, for the convergence of finite dimensional distributions, with t 1 < t 2 < · · · < t k and α ℓ , ℓ = 1, . . . , k, we can write k ℓ=1î α ℓĨ N,0,1 1,1 (t ℓ ) using the incrementsĨ N,0,1 1,1 (t ℓ ) −Ĩ N,0,1 1,1 (t ℓ−1 ), and carry out the calculations as shown above.

Next, to prove tightness, we employ Theorem 13.5 and verify condition (13.14) in [6] . We show that for r ≤ s ≤ t and for N ≥ 1,

for some constant C and φ(t) = t 0 q 1,1 (u)F 0 (du) which is a nonnegative, nondecreasing and continuous function. Recall that F 0 is assumed to be continuous. This will enforce condition (13.14) in [6] , which according to Theorem 13.5 implies tightness in D. Let Note that E[∆I k r,s ] = 0, E[∆I k s,t ] = 0,

By direct calculations, following similar steps in the proof of (14.9) in [6] , we obtain , i ∈ L , we obtain tightness from that of each process as established above, so it suffices to show the joint convergence of their finite dimensional distributions. Take ℓ = 1, i = 1, 2 as an example. For 0 < t 1 < t 2 and α 1 , α 2 ∈ R, E exp îα 1Ĩ ⇒ 0 in D. We have

It is clear that by Assumption 3.2,

N is tight, by Assumption 3.2 and (5.15), it suffices to prove the tightness of the first term on the right hand side of (5.15), which we denote as ∆ N,0,1 1,1 (t). By the Corollary of Theorem 7.4 in [6] , it suffices to show that for all ǫ > 0,

Since ∆ N,0,1 1,1 (t + u) is increasing in u, we only need to consider

whose lim sup N is equal to zero under Assumption 3.2, and thus we conclude that (5.16) holds. Therefore we have shownÎ N,0,1 1,1 −Ĩ N,0,1

→ 0 in D in probability, and conclude the joint convergence Î N,0,1

To prove that the limit processes Î 0,1 ℓ,i , ℓ, i ∈ L are continuous when F 0 is continuous, since they are Gaussian, it suffices to show continuity in the quadratic mean [19] , that is, for all t > 0, lim s→t E Î 0,1 ℓ,i (t) −Î 0,1 ℓ,i (s) 2 = 0. This is easily checked from the continuity of the covariance functions. We next focus on the processes Î N,0,2 (5.10) . We first prove the joint convergence Ĩ N,0,2 ℓ,i , ℓ, i ∈ L ⇒ Î 0,2 ℓ,i , ℓ, i ∈ L in D L 2 as N → ∞. We again apply Theorem 13.5 in [6] . By direct calculations, we obtain for t ≥ 0,

and for t ′ ≤ t ≤ t ′′ and α 1 , α 2 ∈ R,

Hence, we can establish the convergence of finite dimensional distributions ofÎ N,0,2 ℓ,i similarly as that ofĨ N,0,1 1,1 (t) above. For tightness, we obtain for t ′ ≤ t ≤ t ′′ and for N ≥ 1, 

it suffices, applying the Corollary of Theorem 7.4 in [6] to the first term, denoted by ∆ N,0,2

Since it is increasing in t, we consider for δ > 0, , ℓ, i ∈ L and moreover, since they are tight individually, it suffices to

Therefore we have shown the joint convergence of Ê N,0 ℓ,i ,Î N,0,2 ℓ,i , ℓ, i ∈ L . Finally for the continuity of the limit processes, it suffices to show the continuity in the quadratic mean [19] , which follows from the continuity of the covariance functions. This completes the proof of the lemma.

For the next lemma on the moment estimates, we shall need the following technical result. Lemma 5.3. In the two cases γ ∈ [0, 1) and ℓ =i κ iℓ = 0, there exists a constant C such that for any 0 ≤ t ≤ T , 1 ≤ i ≤ L, 

Moreover, if we define for 0 ≤ a ≤ 1, g(a) = ψ(s + a(s ′ − s), e + a(e ′ − e), i + a(i ′ − i), r + a(r ′ − r), u + a(u ′ − u)), we have

We havē

Suppose first that γ < 1. Clearly, the result will follow from the last formulas, if we prove that there exists C > 0 such that for all t ∈ [0, T ], N ≥ 1,

We have

for any 0 ≤ t ≤ T , where we have used Lemma 4.10 for the last inequality.

It is easy to check that in the case where the variable u disappear from the above formulas, the derivatives of ψ are bounded on [0, 1] 3 , and the result holds in the case γ = 1 as well.

We will now prove the following estimate. 

Proof. We first show that for each i ∈ L,

We shall use (5.18) . In the representations ofŜ N i (t),Ê N i (t),Î N i (t) andR N i (t) in (5.4), (5.5), (5.8) and (5.13) , respectively, the following hold: there exists a constant C > 0 such that for each i, ℓ ∈ L,

with the corresponding ν Z = ν S,i,ℓ , ν E,i,ℓ , ν I,i,ℓ , ν R,i,ℓ . Moreover, it is easy to check that

Thus, by taking squares of the processesŜ N i (t),Ê N i (t),Î N i (t) andR N i (t) in (5.4), (5.5), (5.8) and (5.13), we can apply Cauchy-Schwartz inequality and Gronwall's inequality to conclude that (t) is driven by the sequence of two dimensional r.v.'s (ζ 0 k,ℓ , Y 0,k ℓ (ζ 0 k,ℓ )), k ≥ 1. We add a sequence of i.i.d. r.v.'s, globally independent of the above sequence, U k , k ≥ 1, which all have the uniform distribution on the interval [T, T + 1]. We define a sequence of r.v.'s ζ 0 k,ℓ , k ≥ 1, as follows (noting thatζ 0 k,ℓ depends on i, which we omit for brevity)

By writing

we apply the Dvoretsky-Kiefer-Wolfowitz inequality (with Massart's optimal constant [27] ) and

On the other hand,Î

Combining the above, we have shown that the property in (5.19) 

. We also add a sequence of i.i.d. r.v.'s globally independent of the previous sequence, U k , k ≥ 1, uniformly distributed on [T, T + 1]. Define

(Note that ς 0 k,ℓ,ℓ ′ depends on i, which we omit for brevity.) Then we have

Note that Φ 0 ℓ,i (t) is the sum of the right hand side of the above equation over ℓ ′ ∈ L, as given in (3.7). Then we can writẽ

We can apply the Dvoretsky-Kiefer-Wolfowitz inequality to obtain the desired estimate forĨ N,0,2 ℓ,ℓ ′ ,i (t) for each fixed ℓ, ℓ ′ , that is,

Hence,

is again easy to treat. So we obtain the estimate forÎ N,0,2 ℓ,i (t).

We next consider the processes associated with the newly infected individualsÊ N ℓ,i (t) andÎ N ℓ,i (t). Recall the expression ofÊ N ℓ,i (t) andÎ N ℓ,i (t) in (5.7) and (5.11), respectively. Recall the expressions in (2.11) using the PRMQ ℓ (ds, du, dy, dθ) and (2.12) using the PRMQ ℓ (ds, du, dy, dz, dϑ, dθ). Also recall that Q ℓ and Q ℓ are the corresponding compensated PRMs. Thus we can writê

whereῩ ℓ (t) is given in (3.22) and is deterministic. It is not hard to check thatẼ N ℓ,i (t) (resp.Ĩ N ℓ,i (t)) is a martingale w.r.t. the filtration F E t (resp. F I t ), where F E t is generated by the restriction ofQ to the set of (s, a, u, θ) which are such that s + u ≤ t, and F I t is generated by the restriction ofQ to the set of (s, a, u, v, θ, ϑ) which are such that s + u + v ≤ t. Those martingales have the quadratic variations

Then by Doob's inequality, we obtain

We next show that

Let us establish the first estimate in (5.29) . The second can be obtained by the exact same argument. We will use below the identity

We havê

Thus we obtain

hence the first part of (5.29), thanks to (5.20) . Plugging the above estimates in (5.4), (5.5), (5.8) and (5.13), using (5.18) and Gronwall's Lemma we finally establish (5.19 ).

In the next Lemma, we shall need the following well-known result on integrals with respect to a PRM, which follows rather easily from Theorem VI.2.9 in [11] .

Lemma 5.5. Let Q be a PRM on some measurable space (E, E), with mean measure ν, and Q the associated compensated measure. Let f : E → C be measurable and such that

We are now ready to prove the convergence of the components associated with the newly exposed individuals.

where the limits are as given in Theorem 3.2.

Proof. Recall the processesẼ N ℓ,i andĨ N ℓ,i defined in (5.26) and (5.27) using the compensated PRMs Q ℓ and Q ℓ , respectively. Each of these processes being a martingale, they are easily shown to be tight. We now establish their joint final dimensional convergence. By their definitions of the two PRMs, we can regardQ ℓ (resp. Q ℓ ) as the image ofQ ℓ (resp. Q ℓ ) by the projection Π from R 4 + ×L 2 onto R 3 + × L, defined by Π(s, a, u, v, θ, ϑ) = (s, a, u, θ). In other words, we can write, together with (5.26) 

Consequently, for any α E , α I , α ′ E , α ′ I ∈ R, and for each ℓ, i, i ′ ∈ L and t, t ′ > 0,

By Lemma 5.5, we have We can decompose sign(Ῡ N ℓ (s) −Ῡ ℓ (s)) = 1ῩN ℓ (s)−Ῡ ℓ (s)>0 − 1ῩN ℓ (s)−Ῡ ℓ (s)<0 , and writeΥ N ℓ (s) = Υ N ℓ (s) + −Υ N ℓ (s) − , such that each of these will induce a process that is nondecreasing in t. It is also clear that tightness of these processes will be implied by the tightness of the following processes: Since these two processes are nondecreasing in t, by the Corollary on page 83 in [6] , it suffices to show that for any ǫ > 0, and ι = 1, 2, lim sup N →∞ 1 δ P Ξ N ι (t + δ) − Ξ N ι (t) > ǫ → 0 as δ → 0. It is clear that E Ῡ N ℓ (s) −Ῡ ℓ (s) → 0 as N → ∞ by the convergenceῩ N ℓ ⇒Ῡ ℓ and the dominated convergence theorem. Thus, the first and third terms converge to zero as N → ∞. Thanks to (5.19) , δ −1 times the second term converges to zero as δ → 0, which is exactly what we wish. Thus, in order to prove (5.31) Hence we obtain (5.33) by the same argument as in (4.22) , and the bound in Lemma 5.4. For the process Ξ N 2 (t), we have The argument for these two terms follow from that for the second and fourth terms above for Ξ N 1 (t). We next prove that as N → ∞,Î N ℓ,i (t) −Ĩ N ℓ,i (t) → 0 in probability, locally uniformly in t for each ℓ, i ∈ L. It follows a similar argument so we only highlight differences below. It is clear that

where the convergence holds by Theorem 3.1 and the dominated convergence theorem. To show that the sequence {Î N ℓ,i −Ĩ N ℓ,i } is tight, we writê 

By the monotone property of these two processes in t, we then show that for any ǫ > 0, and ι = 1, 2, lim sup N →∞ 1 δ P I N ι (t + δ) − I N ι (t) > ǫ → 0 as δ → 0. (5.34) Similar to the derivation in (5.32), we obtain

(5.35)

The first and third terms converge to zero as N → ∞ by the convergence E Ῡ N ℓ (s) −Ῡ ℓ (s) → 0 and applying the dominated convergence theorem. For the last term in (5.35), we have 

, Similarly for ψ N e,a (t), ψ N i,a (t), ψ N r,a (t) and ψ N u,a (t). Also write ψ e,0 (t), ψ i,0 (t), ψ r,0 (t), ψ u,0 (t), ψ u,0 (t) when a = 0 (noting that they no longer depend on N in this case). We then havê We can use the bounds in the proof of Lemma 5.3 to show that lim sup N sup t∈[0,T ] 1 0 ψ N s,a (t)da < ∞ and the same holds for ψ N e,a (t), ψ N i,a (t), ψ N r,a (t) and ψ N u,a (t). It is also clear that 1 0 ψ N s,a (t)da−ψ s,0 (t) → 0 as N → ∞, uniformly in t, and similarly for the others. In addition, similarly as Lemma 5.3, we can also show that sup N sup t∈[0,T ] E S N i (t) 2 < ∞ and the same for E N i (t), I N i (t) and R N i (t).

Then by Gronwall's inequality, we conclude that ( S N i −Ŝ N i , I N i −Î N i , R N i −R N i , i ∈ L) ⇒ 0. This completes the proof of the theorem.

We have extended the approach in [31] to study multi-patch SEIR models with Markov migrations among patches. In this generalization we have introduced a new formulation for the infection process, and further developed the methodology to tackle the challenges arising from that and the migration processes. However we have assumed a constant infectivity rate for each individual. In [16, 29] , in a model with homogeneous population, each individual is associated with a random infectivity function, for which FLLNs and FCLTs have been established. It would be interesting to study multi-patch models with varying infectivity. In addition, with an infection age dependent infectivity, FLLNs and PDEs have been established for the one-patch models in [30] . PDEs for the multi-patch models with an infection age dependent infectivity can be also derived. Control strategies such vaccination and isolation have been developed for epidemic models [1, 39, 20, 7] . For multi-patch models, one may also consider control strategies restricting migrations in different patches.

i (t) +îα IÎℓ,i ′ (t) +îα ′ EÊℓ,i (t ′ ) +îα ′ IÎℓ,i ′ (t ′ )

An optimal isolation policy for an epidemic

Stochastic epidemic models and their statistical analysis

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

Multi-patch and multi-group epidemic models: a new framework

Convergence of probability measures

Optimal control of epidemic size and duration with limited resources

On a nonlinear integral equation for population growth problems

Stochastic epidemics in a homogeneous community

Probability and Stochastics

SIR epidemic models with general infectious period distribution

An epidemic equation with immigration

Limiting behaviour in an epidemic model

Markov processes: characterization and convergence

Epidemic models with varying infectivity

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

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

Optimal control of epidemics with limited resources

An SIS epidemic model with variable population size and a delay

Stability and bifurcation for a multiple-group model for the dynamics of hiv/aids transmission

On the dynamics of a class of multi-group models for vector-borne diseases

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

Final size of an epidemic for a two-group SIR model

Final size of a multi-group SIR epidemic model: Irreducible and non-irreducible modes of transmission

The tight constant in the Dvoretzky-Kiefer-Wolfowitz inequality

Nolinear Volterra Integral Equations

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

Functional limit theorems for non-Markovian epidemic models

Population modeling of early COVID-19 epidemic dynamics in French regions and estimation of the lockdown impact on infection rate

The asymptotic evolution of the general stochastic epidemic

Modelling COVID-19 transmission in the United States through interstate and foreign travels and evaluating impact of governmental public health interventions

A simple SIS epidemic model with a backward bifurcation

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

Optimal isolation policies for deterministic and stochastic epidemics

Acknowledgement. The authors want to thank the anonymous reviewers for their helpful comments that have led to substantial improvements in the paper. G. Pang was supported in part by the US National Science Foundation grants DMS-1715875 and DMS-2108683 and Army Research Office grant W911NF-17-1-0019.

show the convergence of their joint finite dimensional distributions. Note thatẼ N,0 ℓ,i (t) andĨ N,0,2 ℓ ′ ,i (t ′ ) are independent for ℓ = ℓ ′ . We calculate that for α, α ′ ∈ R and t, t ′ > 0, lim N →∞ E exp îαẼ N,0 ℓ,i (t) +îα ′Ĩ N,0,2 ℓ,i ′ (t ′ ) = E exp îαẼ 0 ℓ,i (t) +îα ′Ĩ 0,2 ℓ,i ′ (t ′ )This can be extended easily to finite dimensional distributions of Ẽ N,0This requires calculations of the covariances of the cross terms for the increments ofẼ N,0 ℓ,i (t) andĨ N,0,2 ℓ,i ′ (t). For instance, we have for α, α ′ ∈ R and t 1 < t 2 ,Would we consider more distinct times, we would clearly deduce that the whole vector converges to a Gaussian random vector. The only point which requests a detailed computation is the determination of the covariances, which can be deduced from the above formula, and obvious similar formulas.In particular, it is easily seen for ℓ = ℓ ′ , that the covariance ofÊ ℓ,i (t) and The formulas for the covariances of the pair (Ê,Î) in the statement of Theorem 3.2 are easy to deduce from the above computation.It then remains to show that, for each ℓ, i ∈ L, as N → ∞,We focus on the processÊ N ℓ,i −Ẽ N ℓ,i . It is clear thatwhere the convergence holds by Lemma 4.11 and the dominated convergence theorem. We next show that the sequenceWe define a 4m-dimensional integral mapping ̥: given a i , b i , c i , d i , φ i , ψ i , ϕ i , χ i ∈ D, some constants α i , β i , γ i , κ i > 0 and functions F ℓ,i , G ℓ,i for ℓ, i = 1, . . . , m, let x i , y i , z i , w i be the solutions to the following integral mapping:The existence and uniqueness of its solution and the continuity property are stated in the following lemma.Lemma 7.1. Assume that F ℓ,i and G ℓ,i , ℓ, i = 1, . . . , m are measurable, bounded and continuous functions satisfying F ℓ,i (0) = 0 and G ℓ,i (0) = 0, and let the constants α i , β i , γ i , κ i > 0 and the functions φ i , ψ i , ϕ i , χ i be given. There exists a unique solution (x i , y i , z i , w i , i = 1, . . . , m) ∈ D 4m to the set of integrable equations defining the mapping̥. The mapping is continuous in the SkorohodIn addition, if φ i , ψ i , ϕ i , χ i are continuous, then (x i , y i , z i , w i i = 1, . . . , m) ∈ C 4m and the mapping ̥ is continuous uniformly on compact sets in [0, T ].Proof. For the existence and uniqueness of solutions, we can apply the Schauder-Tychonoff fixed point theorem, and modify the proofs in Theorems 1.2 and 2.3 in Chapter II of [28] (where these results are shown for Volterra integral equations with continuous functions). The continuity can be proved similarly as Lemma 9.1 in [31] .