key: cord-0679910-yxefr0aa
authors: Protin, Fr'ed'eric; Jules, Martel; Nguyen, Duc Thang; Nguyen, Hang T.T.; Piffault, Charles; Rodr'iguez, Willy; Iglesias, Susely Figueroa; To, Tat Dat; Tuschmann, Wilderich; Le, Hong Van; Yeo, Tenan; Nguyen, Tien Zung
title: Unified modelling of epidemics by coupled dynamics via Monte-Carlo Markov Chain algorithms
date: 2021-06-25
journal: nan
DOI: nan
sha: 318af1abd7570aa7025a11035998fc4370ff6fbc
doc_id: 679910
cord_uid: yxefr0aa

To forecast the time dynamics of an epidemic, we propose a discrete stochastic model that unifies and generalizes previous approaches to the subject. Viewing a given population of individuals or groups of individuals with given health state attributes as living in and moving between the nodes of a graph, we use Monte-Carlo Markov Chain techniques to simulate the movements and health state changes of the individuals according to given probabilities of stay that have been preassigned to each of the nodes. We utilize this model to either capture and predict the future geographic evolution of an epidemic in time, or the evolution of an epidemic inside a heterogeneous population which is divided into homogeneous sub-populations, or, more generally, its evolution in a combination or superposition of the previous two contexts. We also prove that when the size of the population increases and a natural hypothesis is satisfied, the stochastic process associated to our model converges to a deterministic process. Indeed, when the length of the time step used in the discrete model converges to zero, in the limit this deterministic process is driven by a differential equation yielding the evolution of the expectation value of the number of infected as a function of time. In the second part of the paper, we apply our model to study the evolution of the Covid-19 epidemic. We deduce a decomposition of the function yielding the number of infectious individuals into"wavelets", which allows to trace in time the expectation value for the number of infections inside each sub-population. Within this framework, we also discuss possible causes for the occurrence of multiple epidemiological waves.

The Covid-19 pandemic has given rise to manifold new studies which aim at understanding the dynamics of the disease spread and predicting its future evolution (see, for instance, [14, 17, 31, 25] ).

In a large population context, several models have been proposed to grasp the phenomena of multiple front epidemic waves. These multiple waves have for example been explained by a diversity of sub-populations in the same locality, with their own characteristics, spreading the disease from one to the other by interacting. The first type of model is for instance described in [13] , where a Forced-SIR model is used to capture the multiple waves, interpreted as subepidemics, underlying a country's overall incidence curve.

Other models rather consider the succession of contamination episodes between geographic areas. For instance, in the seminal work [3] , the authors consider a network model for a pandemic. Each node of the graph is a town, which is driven by a basic SIR model. The edges represent airline connections between towns. The population is mixed by a Markov process, superimposed to the SIR dynamics locally defined for each town. The authors define a pseudo-distance between two towns from the matrix of the Markov process. They show that, for the SARS epidemic of 2003, and for the 2009 H1N1 influenza pandemic, the pseudo-distance from a town to the origin of the pandemic is strongly correlated to the time of contamination of this town. In [4] this model is applied to a graph of towns in Mexico, where the order of contamination of these towns for the Covid-19 epidemic is retrieved. Some models represent in a common framework the heterogeneity of the population due to subdivisions into sub-populations or to geographical dispersion. Thus in [6] , the authors introduce an epidemic model composed of overlapping sub-epidemic waves, each representing the dynamics of groups in the population. They argue that these groups are determined by spatially clustered pockets, population mobility patterns, infections moving across different risk groups, and so on. We refer to [27] for a generalization of sub-epidemic waves using wavelets approach. In particular this latter method gives very reasonable forecasting. We adopt this point of view, and, in a context of large populations, we consider a population moving on a graph, whose nodes can represent geographical zones as towns or countries swapping travellers, and sub-populations more or less close in contacts. Modeling the spreading of an epidemic in a small population requires a different treatment, as the stochastic effects become important. Here [5] proposes a discrete SIR model to estimate the evolution of the expectation of the number of infected individuals for measle epidemic. In the same way as in in [5] , T. Yeo and his colleagues [20] give both a stochastic and deterministic model of an epidemic spreading. The authors describe a spatial model for the spread of a disease on a grid of a bounded domain. The stochastic model consists of a random Markov epidemic model as a Poisson process driven stochastic differential equation. They prove that the stochastic model converges to the corresponding deterministic patch model, as the size of the population tends to infinity. On the other hand, they show that the stochastic model converges to a diffusion SIR model as the size of the population tends to infinity and the mesh of the grid goes to zero. We propose here to bring the contexts of small and large populations together in a common framework. First we present a discrete stochastic model of an epidemic spreading on a graph. It uses a Monte-Carlo Markov Chain method for simulating the displacement of individuals according to given probabilities for the individuals to be in specific node at certain time step. Then we prove that when the length of the discrete time step converges to zero, this process is driven by differential equations giving the evolution of the expectation value of the number of infected and susceptible individuals as a function of time. We interpret the resulting differential equations as classical local SIR dynamics, coupled by a diffusion operator, associated to the diffusion of the virus by the interaction of individuals. In a large population context, nodes of the graph represent geographical entities or subpopulations. It is shown that, when the size of the population increases, under a natural hypothesis, the stochastic process associated to this model converges to a deterministic process. Thus this model can be used to grasp either the geographic evolution of an epidemic, or the evolution of an epidemic in a heterogeneous population, divided into homogeneous sub-populations, or, more generally, a mix of these two contexts. In a small population context, the differential equations system driving the expectation of the number of infected and susceptible individuals is still valid, but the diffusion operator has a slightly different interpretation. It expresses the occupation time of agents in each node interpreted as a geographical zone, such as a room in a building or a ship, according to a statistical schedule. Section 2 contains the description of the model and several theoretical derivations. After motivating the model in Subsection 2.1, and describing the Monte-Carlo Markov chain (MCMC) algorithm adapted to our context in Subsection 2.2, the model in abstract form is presented in Section 3 in the case where the graph is reduced to one node, for simplicity and first derivations. In particular, it is shown that under a natural hypothesis, when the size of the population increases, the stochastic process converges towards a deterministic process, in a mathematically precise sense (Theorems 2, 3).

We explain how to generalize these results in order to apply them to variations of our model, which, in particular, include SIR models taking into account saturation effects in the contamination process, as proposed for instance in [32] in the time-continuous case. An estimation of the rate of convergence is also given. In response to the obvious stability question whether our proposed discrete approach also does allow for an interpretation and version in the setting of a continuous time parameter, in subsection 3.3 we show that when the lengths of the time steps on which our model is based converge to zero, the dynamics of the average number of infected and susceptible individuals can actually be expressed by a differential equation (Equation (21)). In addition, we prove there that this is true as well if we allow for the more general setting of an arbitrarily chosen, though, of course, fixed, time delay that will naturally occur when infected individuals pass on to 'recovered" state.

We also propose another method of prediction which is based upon iterating calculations of conditional expectations. We anticipate that this approach can be successfully transferred to the description and study of models of yet higher order of complexity, for which a differential equation reflecting the average evolution of the epidemic would be rather difficult, if not impossible, to obtain. In Subsection 3.4 we prove that the variances of the number of infected, susceptible and recovered people divided by the size of the population converge towards 0 when this size increases and the time step decrease. A precise upper bound is given for the case of recovered people (Proposition 2), according to a method that can be adapted to the case of infected and susceptible individuals. In Section 4, the model is presented in generality for a graph with several nodes. Depending on the context, a node can represent a room in a building or a ship, a geographical zone, or a subpopulation. We explain how the results of Section 2 transfer to this case. We also provide a heuristic interpretation of how the diffusion effects slow down the spread of an epidemic. A model for several groups, making sense in a context of small population, is also presented in Subsection 4.3. Numerical experiments are presented in Section 5. In particular we present a fitting of the model by data from the Covid-19 in Singapore. We also show how our model allows to differentiate the effects of lockdown or social distancing measures from the effects of limiting movement measures on a large scale, imposed by, say, cancellation of plane flights or a ban on changing geographic regions for large numbers of individuals.

Related works: Simulations using network models have been used as an effective tool to study the properties of epidemics, see for example [2] and [28] .

Recently some authors have implemented network or agent-based models to simulate and study the spread of Covid-19 through various populations. In [24] the authors introduced a method for modeling disease transmission dynamics through a SEIR model which relies on a contact network between individuals rather than cities or large populations. In [21] the authors used a network model to study the spread of Covid-19, in particular this contact network was between cities rather than individuals. In our present work, we introduce a common framework containing both contexts of small and large populations.

We also refer [11] where the authors used a network-based model to study the efficacy of various interventions in France, including national lockdown, mask wearing, distancing measures. They also studied the impacts of these measures on ICU bed occupancy, numbers of total cases, and numbers of total deaths.

We present here a simulation intended to motivate our model in a small population context.

The MCMC method (Monte-Carlo Markov Chain) is traditionally used to simulate a given law on a very large state space. More precisely, the method is based on the simulation of a Markov chain converging to this law when the number of iterations tends to infinity. We refer for example to [23] for details. In the present study, we propose to use the MCMC algorithm for a different purpose, namely to simulate a random walk on a given graph, with prescribed limiting probabilities of being in each node.

For this simulation we describe the process as follows. Two populations of 100 individuals each, are moving between the rooms represented by nodes in a graph (Figure 2 .1a). Each individual is an object in the sense of OOP (Object-oriented programming), with the following attributes: state (S, I, or R, for Susceptible, Infectious or Recovered, respectively), duration of infection for the individuals in the state I, the number of the room where the individual is located, and his group G1 or G2, as defined hereafter.

The individuals of the group G1 move mainly between the rooms 1, 2, 3, 4, while the individuals of the group G2 move between the rooms 4, 5, 6. The two groups meet therefore in room 4. For each of the two groups, we can choose the probabilities for an individual of this group to be in each room. For the group G1, we have taken in this example the probabilities 0.28, 0.28, 0.28, 0.16 for each individual to be in room 1, 2, 3, 4 respectively, and for the group G2, we have taken probabilities 0.2, 0.4, 0.4 for each individual to be in room 4, 5, 6 respectively. Thus each individual is two times less often in room 4, where the groups meet, than in the other rooms that he frequents. The movement of the individuals is then performed independently according to a Markov chain simulated by MCMC method, having these given limiting probabilities to be in each room.

The contamination process is simulated as follows. For each susceptible individual at a given node, the probability that it goes to state I, i.e. it becomes infected, is

where nI is the number of Infected individuals at that node, who remains Infected for more than 5 steps (to represent a latency time). Note that this choice is equivalent to suppose the infection of a Susceptible with probability 1−e −0.01 by going through the Infected successively. It is also equivalent to transform a proportion of susceptible into infected according to a binomial law with parameters nS, the number of susceptible individuals, and 1 − e −0.01 * n I . Furthermore, each Infected can pass to Recovered state either with a probability of 0.05, or after the infection time exceeds 20 steps, and it applies to all nodes. The simulation for this process is shown in Figure 2 .1b. The two peaks in Figure 2 .1b are robust, in the sense that they appear in almost every simulation. This shape seems to depend more on the topology of the graph, the choice of the parameters, and the rules of attendance of the rooms by each group, than on pure luck. With the same conditions as above, now suppose that groups 1 and 2 meet at nodes 2, 3 and 4, and no longer only at node 4. For the group 1, we take probabilities 0.17, 0.17, 0.49, 0.17 for each individual to be in room 1, 2, 3, 4 respectively, and for the group 2, we take probabilities 0.14, 0.44, 0.14, 0.14, 0.14 for each individual to be in room 2, 3, 4, 5, 6 respectively. We show this simulation in Figure 2 .2. The groups are thus more mixed than previously, and we note that there is now only one peak (Figure 2 Reasoning for (1). We give here some heuristic reasons for the choice of (1) in the previous model.

Recall that, in a given room, an infectious individual contaminates a susceptible one according to the value of a random variable with a Bernoulli distribution. We wish to estimate the parameter of this distribution, which depends on the number of infectious individuals in the room. Suppose that the process evolves in continuous time. It can be expected that the time T before a given susceptible individual becomes infected by a given infectious individual follows an exponential distribution, that is, P(T < t) = 1 − e −λt , for all t > 0, for some parameter λ > 0. In this case, the number of infectious individuals for each unit of time follows a Poisson distribution with parameter λ. Now suppose that a number nI of infectious individuals can contaminate a susceptible individual S. For every infectious individual i = 1, . . . , nI , let us denote by Ti the time before this individual infects S. Then the contamination of S occurs in a time step ∆t if min i=1,..,n I Ti ≤ ∆t. Since the Tis are supposed to be independent, min i=1,..,n I Ti follows an exponential distribution with parameter nI λ.

Indeed, for all t > 0,

Thus we retrieve (1): the probability of infection in a given unit of time for each susceptible in a room is p λI = 1 − e −λn I , for some parameter λ.

Finally, it is straightforward to see that (1) is invariant by a change of time scale, which is an additional motivation for this choice.

Remark 1. We have supposed so far that the infection time T takes values in R + . It is straightforward to check that if T follows an exponential distribution with parameter λnI t, then the integer part T follows a geometric distribution with parameter 1 − e −λn I t . We then find (2) once again by reasoning in discrete time.

Thus, a first idea is to consider that the probability for a susceptible to be infected in a room before time t is

The parameter λ = λ(N ) depends on the fixed amount N , the total size of the population. For instance, in the basic differential SIR model it is inversely proportional to N (see e.g. Equation (2.1) in [10] ).

In this section we aim to explain the details of the steps in the MCMC algorithm and how to use it in our context. Let E = {1, ..., N } be the set of nodes of a graph G, and π = (π1, ..., πN ) the probability measure we want to simulate. In particular, we have πi ≥ 0 for all i ∈ E, and i∈E πi = 1. We can suppose without loss of generality that ∀i ∈ E, πi > 0. Recall that a matrix P is said to be aperiodic if, for all i ∈ E, the greatest common divisor of {n : (P n )i,i > 0} is equal to 1. This is equivalent to the condition (P n )i,i > 0 for all i ∈ E and n ∈ N large enough. The following important result is well known (see e.g. Chapter 18 in [15] ). Theorem 1. Let P be a transition matrix for which π is invariant, i.e. πP = π, and let (Xn) n∈N be a Markov chain on E with transition matrix P . If P is irreducible and aperiodic, then (Xn) n∈N converges in law towards π when n → ∞, for any initial distribution X0.

In fact, we have the following sufficient condition:

Proof. (of Proposition 1) The statement follows from the fact that summing over i ∈ E the equation πiPij = πjPji gives πP = π.

Let B be some irreducible transition matrix on E such that

Suppose also that B is not reversible for π. Then, we can construct a matrix Q from B as follows, for all i ∈ E:

Thus, it is straightforward to see that the matrix Q ij is reversible for π. We deduce the following algorithm that starts by taking B as the normalized adjacency matrix of a connected graph, such that (3) is satisfied, and returns a transition matrix Q which is irreducible, reversible for π and aperiodic. Note that B is irreducible since the graph is connected.

• Choose a connected graph on E.

• Let A be the adjacency matrix of this graph, normalized such that the sum of each row equals 1.

• In the case where A is already reversible for π, it is not necessary to use the construction in (4), and we take Q = A for the transition matrix. If A is not reversible for π, then we take Q given by (4) from A. It remains to possibly modify the matrix Q to make it aperiodic, which is the subject of the next step.

• If A is aperiodic, then return Q = Q . Otherwise, we can slightly modify Q in order to become aperiodic and remain irreducible and reversible, for example by returning Q = ε Id + (1 − ε)Q , for some small ε > 0, where Id is the identity matrix of the same order as Q .

We suppose here and in the sequel that the graph is non-oriented, in the sense that the matrix Q is symmetric. This corresponds to an unbiased diffusion of the agents in the graph. Equation (4) is then simplified. Moreover, by adding arrows at each node to itself (see the next paragraph for details), the matrix Q becomes aperiodic, so the last point in the previous algorithm is removed.

When we use the algorithm given above to construct the matrix associated to the MCMC approach from an adjacency matrix of a graph and a given probability law on its nodes, we can choose a nonzero value for the diagonal elements of the adjacency matrix. From the point of view of the graph, this comes down to authorizing the passage of a node towards itself. An agent is considered to be stationary during a time step, when he takes such a route from a room to itself. Thus, multiplying such passages by increasing the diagonal values of the previous adjacency matrix, allows us to slow down the speed of diffusion from a room to another one. More precisely, let Qn be a matrix constructed as previously indicated, that is the reversible transition matrix of a Markov process on a graph such as the frequency of occupation of each node has been prescribed. Here n is the common multiplicity of the nodes of the graph, considered s a parameter, while the elements in Qn out of the diagonal are fixed. Then for a fixed probability vector X, we verify

where

is the Laplacian matrix associated to Q1, the matrix Id being the identity matrix, and n denotes the value of the diagonal elements of the adjacency matrix of the graph, considered as a parameter. HereÕ(·) denotes a matrix of functions O(·) that does not depend on n. Indeed, we deduce from the algorithm presented in Subsection 2.2:

Thus 1 n appear in (5) as a coefficient of diffusion. Note that 1 n+1 is also the probability for each time step that an individuals goes from one node to another.

The model described in this subsection is suitable for small populations or large populations, but the probability vector π appearing in Theorem 1, from which the matrix Qn is built according to the algorithm described in this subsection, receives in these two contexts a different interpretation. In a small population, where, for example, individuals move within a building or a ship, the vector π prescribes the probability for each individual of being in any one of all of the given specific nodes. In the context of large populations, the nodes of the graph represent localities, or sub-populations. In the latter case, we can take as invariant vector, in order to build the matrix Qn, the vector whose coordinates represent the size of the population at each node of the graph. Indeed, here we are looking less at the movements of individuals than at those of the virus itself. An Agent-based simulation of the model would in this case take too much time, but we will establish a differential equation depending on the matrix Qn, which describes the evolution of the epidemic in continuous time.

We want to estimate the parameter n in order to simulate a realistic displacement. Let ∈ (0, 1] denote the frequency of movement, that is, the average number of displacements for each individual by unit of time. Denote by kstep the number of time steps occurring during a simulation of duration duration. We can estimate the probability that individual passes from a node to another during one time step by duration × kstep , and thus we can set

Let us show how the matrix Qn constructed in subsection 2.2 depends on the length of the time step h. Using for h the expression

Equation (8) gives = 1 nh . Thus (7) can be rewritten as

whereÕ(·) denotes a matrix of functions O(·) that does not depend on . Note that

In view of the dependence of h and kstep expressed in (9), it follows that

Thus the discrete Markov semigroup generated by the matrices Q h converges towards a continuous Markov semigroup when the absolute value of the length of the time step converges towards 0.

Remark 2. We way interpret the formula (11) as follows. For a fixed individual, we assume that the time before he passes from one node to an adjacent node follows an exponential law with parameter 1 n . Then for small times the number of individuals passing from a room to the next one follows a Poisson process with parameter 1 n . Hence the discrete model is embedded into a continuous time model by using P

in (5), where X(t) is a Poisson process with parameter 1 n . We then conclude that our equation (11) corresponds to (8.3.7) in [16] .

In this section, the model in its abstract form is presented. It is also discussed and analysed in the case where the graph is reduced to one node. The case of several nodes requires using the MCMC algorithm presented in Subsection 2.2, and will be presented in Section 4.

In regard to the ideas introduced in Subsection 2.1, we state here the model in an abstract manner. We first suppose that there is one node, and look at the number of infectious individuals in this node at each step. Suppose for simplicity that the incubation time is null. Let S h t and I h t be the number of susceptible and infectious individuals at time t, where h ∈ R + * is the time step. For simplicity we also do not indicate explicitly the dependence of the variables S h t and I h t on N . Let (Xi) i∈N be a sequence of random variables taking values in {0, 1}. We suppose that the conditional distribution of Xi given It is Bernoulli with parameter 1 − e −λI t h , for some h > 0. This is the formula (2) . The interpretation of the variables Xi is as follows. It may be assumed that for each time step the set of susceptible or infectious individuals is ordered. For all i ∈ N, Xi = 0 if the i-th susceptible at time t remains in this state at time t + h, and Xi = 1 if it becomes infected (and infectious). Denote by (Yi) i∈N a sequence of random variables with Bernoulli distribution, and parameter

for some fixed γ > 0. Each Yi is supposed to be independent of It for all t ≥ 0. The interpretation of the variables Yi is as follows. Yi will take the value 1 each time that the i-th Infected at time t becomes Recovered at time t + h. If h is chosen much smaller than the mean recovering time, we can write

Xi (12) Note that for all t = kh, k ∈ N, h ∈ R * + , the variables I h t and S h t depend on the time step h, but for convenience, we do not show this dependency in the notation, except when we want to underline it, in which case we will write I h t and S h t . Note that the dependence on h is specified in Corollary 5. I would emphasize here that by construction, according to (12) , the sequence I h kh , S h kh k∈N is a homogeneous Markov chain.

Notations. In the sequel, we denote by N = I h t + S h t + R h t the total number of individuals in the population, which is independent of t. Here R h t is the number of recovered individuals, recursively defined by

Yi.

The notation h → O(h) denotes a function that converges to 0 when h → 0, whose value is unimportant and may change from one location to another, even within a line. We use also the notation ON (·) when it makes sense to emphasize the dependence on N .

Remark that all these matrices are inversible. Let us also define the forecast

.

(In order to keep notation simple, we do not indicate the dependency on the parameter γ). Let us give another expression of the function forecast λ,h . Taking the conditional expectation in (12) with respect to variables S h

Computing by the same way E S h t+h |I h t , S h t , we obtain

for all k ∈ N and E(·) the usual expectation. Note that a rescaling can be done as follows:

3.2 Convergence of the process.

In this subsection we continue to consider the abstract setting in Subsection 3.1. We will establish that the process normalized by the amount of the population N converges towards a deterministic process when N grows, and certain conditions are given. Let us state these conditions. The first one reads

for some constant λ > 0. This inequality becomes an equality in the basic differential SIR model (see e.g. Equation (2.1) in [10] ).

The second condition, also satisfied by the basic SIR differential model (see e.g. Equation (2.1) in [10] ), is:

The quantity N λ(N ) converges when N → +∞.

The convergence of the process under Condition (18) is specified by Theorem 2 below. We will use the following Lemma. For r ∈ N, let us denote by F h r the σ-algebras generated by the variables S h kh and I h kh for k ∈ N, k ≤ r.

is a martingale for the filtration F h k k . (We use the convention that

Proof. Indeed, we have for k ≥ n + 1 

We can now establish:

Theorem 2. Suppose that Condition (18) is verified. Then for all t ∈ Nh the relative number of infectious individuals I h t N becomes in the following manner more and more closer to

as N becomes larger: There exists a constant K > 0 independent of N , ∀n ∈ N, ∀ > 0,

Proof. Denote by P0 : R 2 → R the projection to the first coordinate, i.e. (x, y) → x, ∀(x, y) ∈ R 2 , and let P1 : R 2 → R the function that returns the second coordinate. Fix h > 0, and denote by

Then Lemma 4 implies that there exists a constant C > 0 independent of N such that

The penultimate equality follows from (16) , and the last equality follows from (12) . 

In other words, ∀ > 0,

Here T h n denotes the σ-algebras generated by the variables S h nh and I h nh . In particular, after taking expectation,

is a martingale by Lemma 1. Then a concentration inequality for martingales (Theorem 33 in [9] , see also p.14 of this reference for the notations), together with Lemma 4, provided that the constant C > 0 has been chosen sufficiently large, allows us to conclude that for all > 0, we have

We obtain by exactly the same way, ∀ > 0,

In particular, for all t ∈ Nh the variables Proof. Define for q ∈ R the matrix

In the model described by (12) , each individual in a state S, I or R has a certain probability to change his state given by the matrix Q λN,γ,h I h t N , indexed in the same order. In the terminology of [29] , the proportion of infectious individuals I h t N is the state of the system at time t. Noticing that, for q ∈ R + being fixed, the matrix Q λN,γ,h (q) converges uniformly in R 9 when N → +∞ by (19), a Mean Field limit theorem (Theorem [MF:Thrm] in [29] ) implies the statement. Theorem 2 and Proposition 3 together imply to the following statement: Corollary 1. Suppose that Conditions (18) and (19) are verified. Then, for all t ∈ Nh, the expectation of the variables S t N and I t N converge when N → +∞. Moreover, these variables themselves converge in probability towards these respective limits when N → +∞.

Proof. It suffices to note that the initial conditions S 0 N and I 0 N in the statement of Theorem 3 can be made arbitrarily close to 0.

Generalisations. Now we shall propose a generalization of the results concerning the model (12) obtained in this subsection to the case of the following variations of this model. If the function N → λ(N ) is replaced by a time-dependent random variable (t, N, S h t , I h

Considering further the abstract setting in Subsection 3.1, in this subsection we will derive differential equations driving the dynamics of the expectation of the number of susceptible and infectious individuals given by (12) when the time step h converges to 0 (Proposition 1). Indeed, it is expected that this passage to limit in the model is needed to represent the continuity of the real time parameter. Taking the conditional expectation with respect to variables S h t , I h t , for some t ∈ Nh, gives

Remark 4. This equality can be used from an instance of time t to a future instance t + h, as long as h is much smaller than the mean recovering time. Associated with this estimation, the equation (12) give the variance

. We will actually see in Proposition 5, Appendix B, how to use (20) for forecasting in more distant future. (20) the conditional expectation with respect to the variable I h t we obtain:

Taking expectation, we have

Denote by i(t) := lim 

Dividing by h and letting h → 0 + , we obtain, thanks to Corollary 4 in Appendix B, and by proceeding in the same way from the second equation in (12) , the following result: Proposition 1. The functions t → i(t) and t → s(t) are differentiable and satisfy the differential equations of the basic SIR dynamics, that is to say

We thus retrieve in average the well-known basic continuous SIR dynamics (see e.g. Equation (2.1) in [10] ). We will see in Section 4 that when the graph underlying the model has several nodes, diffusion effects appear, that slow down the spread of the epidemics.

Introduction of a time delay. One may wish to introduce a delay between the contamination of a susceptible and its contagiousness. This delay can be used to model, for instance, the case of individuals infected with Covid-19. We can suppose that this delay is a random variable T , independent of S h t and I h t for all t ∈ R + . Then (12) can be replaced by

where PT denotes the measure PT (A) := P(T ∈ A) on R and * denotes the convolution product. For instance, if T = t0 is constant, we have

Note that this system is similar to (21) with a time delay t0 in this case.

In this subsection, we will give an upper bound on the variance at each time of the number of recovered individuals R h t N in the discrete model expressed by (12) (Proposition 2). Here t = kh for some k ∈ N. Recall that h is the time step in this model, and that N is the total amount of the population. This will imply that the variances of First let us recall a well-known relation between the variance V (X) of a random variable X with finite variance defined on a probability space A, and the conditional variance V B (X) = E B ((X − E B (X)) 2 )according to a σ-field B defined on A (see e.g. [30] p. 385-386.)

For simplicity, we denote by V X (·) the conditional variance according to the σ-field σ(X) generated by some random variable X. For further convenience, we will also use also the notation E X (·) for the conditional expectation E (· | X). We have from (13):

It. Taking expectation, thanks to Lemma 3 we obtain:

Thus

). Thus we have the upper bound:

Then Lemma 3 in Appendix A leads to:

According to the fact that V (R0) = 0, and recalling that t = kh, k ∈ N, we can rewrite (22) using a telescopic sum:

The inequality · 2 ≤ · 1 between p-norms in R k gives the upper bound

To it, we can now apply a discrete Gronwall-type inequality (Consequence 1 of Theorem 105 in [8] ). After dividing by N 2 the terms of the previous inequality, we obtain, provided that h ≤ 1:

Since t = kh, V (It) ≤ N 2 and E(It) ≤ N , we obtain by supposing that h ≤ min 1, 1 γ :

Finally, recalling that V (R0) = 0, a telescopic sum argument gives the following statement:

Proposition 2. For all h ∈ R + such that h ≤ min 1, 1 γ , and all t ∈ Nh, we have the upper bound:

In particular,

A similar method can be used to majorize the variances V 

From

Then Corollary 3 or Corollary 4 in Appendix B, together with (25) and (23), gives the following fact, proving (24):

4 The case of several nodes.

Let us now suppose that there are several nodes in the graph of the model. Here a node can represent a locus such as a room in a context of a small population, or a homogeneous sub-population or a geographical zone, in a context of a large population. Let j = 1, .., n be the nodes, and suppose that each individual moves from a node to another one according to a Markov chain with transition matrix Q. Here Q is the matrix constructed by the MCMC algorithm described in Section 2.2. When necessary to exhibit the dependency of the matrix Q on h, as expressed by (10), we will denote it by Q h as in (10) . We also suppose that the infection rate depends on the node. For instance, if a node represents a room on a ship, the infection rate in this room is depending on its surface and its degree of ventilation. Thus we can assign a weight to each room, and this in such a way that e.g. a room has double weight than another one if its power of infection is the same as in the second one for a susceptible individual staying twice longer. In the same way, if a node represents a sub-population, its infection rate is depending on numerous parameters, such that its density, its mobility, etc.

Let I h j (t) and S h j (t) be the number of infected and susceptible agents at time t in the node j. Denote byS t := (S1(t), S2(t), ..., Sn(t)) T ,Ĩ := (I1(t), I2(t), ..., In(t)) T . The model at each time t is described from the initial number of infectious and susceptible individuals by the following recursive process:

• First, the evolution of I h j (t) and S h j (t) in each node j = 1, · · · , n from time t to t + h is governed by equations (12) . We thus obtain vectorsÎ h t+h andŜ h t+h fromĨ h t andS h t , that express the new number of susceptible and infectious individuals in each node before the displacement of individuals.

• Then, the new quantities of susceptible and infectious individuals are given byS h t+h =Ŝ h t+h T Q andŜ h t+h T Q respectively. This operation expresses the diffusion of the individuals in the graph according to the matrix Q.

We obtain by the same reasoning as before

where λ := (λ1, ..., λn) T ∈ (R + ) n , γ := (γ1, ..., γn) T ∈ (R + ) n , andĨ h tS h t denotes the vector

and so on. For a function f : R → R, the notation f (x1, · · · , xn) := (f (x1), f (x2), ..., f (xn)) T was used in (26) for the exponential function, and it will also be used in the sequel.

Convergence of the process. Let us just indicate how the reasoning of Section 3 can be adapted in our context to get the same type of results. Denote by

and by forecast λ,γ,h (u1, v1, · · · , un, vn) := N λ,γ,h (u1, · · · , un) u1 . . . un v1 . . . vn Q h .

This last formula is reduced to (14) when there is a single node. A proof similar to the one of Lemma 1 shows that the following sequence is a martingale:

Thus the statements of Theorem 2 and Theorem 3 and their proofs carry over to our present situation. We conclude as in Subsection 3.1 that if each coordinate of the vector λ satisfies Conditions (18) and (19) , then for all t ∈ Nh, the variablesS t N andĨ t N converge in probability towards a constant when N → +∞.

In this subsection we shall derive a differential equation expressing the evolution of the expectations ofSt andĨt. According to (11) we can suppose that there is a semigroup of stochastic matrices t → Q t , with Q 0 = I, Q h = Q, and Q t = I − t∆ T + o(t). Here ∆ is the Laplacian matrix defined in (6) , and ≥ 0 is a parameter that expresses the scale between the time of displacement and the time of the contamination process. Denote byĩ(t) := lim E(S h t ). These quantities are well defined thanks because of Corollary 2 below. Indeed, thanks to the commutativity of the operations of taking expectation of a vector of random variables and multiplying by the matrix Q h , the steps leading to Corollary 6 in Appendix B are still valid. It can be restated as:

Corollary 2. ∀t ∈ Nh, the limits lim l→0 + E(Ĩ l t ) and lim l→0 + E(S l t ) exist, and the convergence is uniform. Moreover lim

By letting h → 0 in (26) after taking expectation, as previously in subsection 3.3, we obtain from Corollary 2 the following result:

Proposition 3. The vector valued functions t →ĩ(t) and t →s(t) are differentiable and satisfy the following differential equations, that reduce to (21) when the graph underlying the model has a single node:

where ≥ 0 is a parameter that expresses the scale between the time of displacement and the time of the contamination process.

Remark 5. By summing up the coordinates in (29), we obtain an equation that explains the slowing down of infection due to diffusion, compared to the classical SIR model with a single compartment, such as expressed in (21) . Indeed, denoting by i(t) and s(t) the total number of infected and susceptible individuals at time t respectively, we obtain:

where · | · denotes the scalar product on R n . The Cauchy-Schwarz inequality implies that the instantaneous speed of the propagation is maximal when the vector of the number of infectious individuals in each node is proportional to that of the number of susceptible individuals, which actually does not happen in reality, since the onset of the epidemic is spatially localized. This fact can be interpreted by saying that when compared to a standard SIR model, the diffusion effects slow down the contagion. The greater the value of the diffusion parameter is, the sooner the diffusion tends to mix individuals and evolves the distribution of individuals closer to homogeneity, and thus approaches the picture given by a standard SIR model. We will see in Section 5 that such diffusion effects can also explain the occurrence of multiple front waves.

Note also that, thanks to (5), for > 0 sufficiently small, the solutions of the equations (29) are well approximated by those of the following equation:

where P is the matrix obtained by the MCMC method from an adjacency matrix of a graph whose diagonal values of are 1 . We recover equations (7) of [4] , together with a statistical interpretation, that allows us to use them in both contexts, i.e. in big or small populations.

Let us illustrate equations (29) with the case where there are just two rooms, each containing initially 1500 individuals. Two individuals are initially infected the room 1. Here P = 9.99999992e − 01 8.00000000e − 09 8.00000000e − 09 9.99999992e − 01 . Remark 6. Let us make a few remarks on possible generalizations. The model of diffusion on a non-oriented graph presented here could be generalized without difficulty to directed graphs, if the need for modeling there arose. Indeed the loss of symmetry in the P matrix does not have an impact on the calculations. The addition of states other than Susceptible, Infected and Recovered, for example by structuring the population by age, is also possible. It seems also meaningful to replace the Markov diffusion process by a "wilder" process such as a Levy Flight in the sense of [18] , by modifying (5) using equations (4.4) and (4.5) in this book. Indeed, this stochastic process, under certain circumstances, models well the human walk [22] .

In the context of a small population, we can consider several groups with a schedule for each, i.e. for each group given probabilities for each individual of being in each node, as illustrated in Subsection 2.1. Specifically, we investigate now K groups instead of one, with transition matrices P 1 , . . . , P K . The same reasoning as in Subsection 3.3 gives

We have denoted byS k t ,Ĩ k t the vectors whose N components are respectively the number of susceptible and infectious individuals of group k for each rooms. (We have not indicated the dependence on the time step h for simplicity). Denote byĩ k (t) := lim h→0 E(Ĩ k t ) and bys k (t) := lim h→0 E(Ĩ k t ) for k = 1, · · · , K. These quantities are well defined thanks to Corollary 2. In continuous time we obtain as before corresponding equations:

5 Numerical experiments. Then, the epidemic spreads in each room according to (12) . The simulation is represented by the blue curve of Figure 5 .1a.

The green curve in the figure 5.1a depicts the predictions (27) and (28), with Q = 0.92 0.08 0.08 0.92 here used in (27) . The formula (28) is applied iteratively until the desired date of prediction, at the jump of the green curve. Red curve is given by the differential equation (32) . Compare the result with the differential equation given by (29) , that is representated in Figure 5 .1b. Here E is the compartment containing the exposed individuals. It is the subset of those individuals who are suffering from the illness, but since they are still in the incubation phase, they are not yet contagious. Let us compare the simulation with the simulation of a model slightly different from (12) , more suited to the study of the evolution of an epidemic in short time, on a one-day scale. Inside a stadium, for example, the contagion may take place within a period of 9 hours. We simulate the spread of an epidemic into a population moving in a structure with 2 rooms. Each room contains initially 1550 individuals. Then 50 individuals are contaminated in room 1. At each step, each individual moves to the other room with a probability of 0.008, or stays. Then, the epidemic spreads in each room according to the following SEI dynamics, where we use the same notations that (12) :

Xi.

Indeed, it is here assumed that, during one day, newly infected individuals are not yet contagious, being in an Exposed state, denoted by E. We suppose also that there is no recovering during a day.

We have used λ = (0.0025, 0.0025). The simulation is shown by the blue curve of Figure 5 .2b. Ten steps are used for drawing this curve, but the number of steps is essentially inessential.

The green curve in the same figure uses for prediction the following formula, iteratively applied, derived as (28) , and using the same notation:

Here P = 0.992 0.008 0.008 0.992 , andĨ = I0 I1

, where Ii indicates the number of infected in room i.

The notation is the same for the vectorsS,Ẽ,R of susceptible, exposed and recovered individuals, respectively. This formula is applied iteratively until the desired date of prediction, here 8h later (the jump before 8h is an artefact of visualisation and obviously does not really exist).

Red curve is the solution of the following differential equation, comparable to (32) :

Now let us illustrate how (29) can be used for forecasting. We pursue the idea of [27] , according to which the evolution of the epidemic can be decomposed into several fronts waves. This approach makes it possible to detect the rise of a wave very early and to estimate its strength. First we collect the number of daily cases in Singapore ( [7] ) for the Covid 19 pandemic. Neglecting the number of deaths, we estimate the number of active cases by:

Active Cases at day t = Cumulative Cases until day t − Cumulative Recovered until day t.

We have also supposed that only a percentage of the daily cases has been really detected. We make this assumption in order to reproduce one more time the phenomenon of asymptomatic cases present in the Covid-19 epidemic. An error function has been implemented, as the mean square of the difference between active cases at Singapore, multiplied by the rate of detection previously mentioned, and solutions of (20) or (26) . The parameters in these equations minimising this error function, including the detection rate, have been obtained by the basin-hopping method. A total population of N = 50000000 has been used, distributed between the 2 subpopulations according to an unknown proportion that is also a parameter of the error function. Figure 5 .3 shows the active cases compared to the model. Data until t = 180, to the left to the dashed green vertical line, was used for training. This method can be an alternative to [27] to detect the onset of a wave. 

Let us first illustrate the effects of multiples nodes or lockdown of the generation of several epidemic waves. Limiting movement measures on a large scale will be modeled by a decrease of the diffusion coefficient ≥ 0 in equation (29), while social distancing and lockdown measures will be modeled by a decrease of the incidence coefficients λ in the same equation. Figure 5 .4 shows left the expected number of infectious individuals given by (29) with 1 node (in this case, = 0). Here the values of the parameters are N = 100000, λ = 0.02, γ = 0.01. The right graphics use the same parameters, except that λ = 0.012 until t = 10000. At this time, it remains 68639 susceptible individuals. Then λ increases to 0.02, simulating an easing of lockdown measures. Then a second wave arises. On the other hand, in the left graphics of Figure 5 .7, measures limiting movements on a large scale are simulated before t = 1100 by dividing by 10 the diffusion coefficient . It can be seen that the height of the contamination peak will be reduced. It can again be seen that the height of the peak of contamination is reduced. We conclude that multiple front waves can occur by diffusion effects, between sub-populations or geographical zones according to the signification of the nodes. These waves can again be multiplied by the effects of distancing or lockdown measures, leading to more complex effects. 

Let us gather here some general lemmas that we have used so far. The first one results from a straightforward computation.

Lemma 2. Let B ⊂ B two σ-algebras defined on a space A, and X a random with finite variance defined on A. Then

Lemma 3. Let B a σ-algebra defined on a space A, and X ≥ 0 a random variable with finite variance defined on A. Then the function B → V B E B (X) is increasing on the set of σ-algebras containing B, and bounded from above by V B (X). Moreover,

is decreasing on the set of σ-algebras containing B, and bounded from above by V B (X).

Proof. Let three σ-algebras B ⊂ B ⊂ B be given. Equation (33) implies

This proves the first statement. The second then follows from (33), that gives

This subsection gives a justification for using (26) for long term forecasting. We refer to Subsection 3.1 for the notations.

Recall that the function forecast λ,h : R 2 → R 2 is defined by (14) . For convenience, for all n ∈ N, we denote by T h n the σ-algebra generated by I nh and S nh . Recall that, being a random variable X, the conditional expectation E X | T h k is a random variable that is a unique measurable function, a.e. defined, of the pairing (S h kh , I h kh ) (see e.g. Proposition 3 in [26] ). Let us denote a continuous representing of this function, if exists, by E X | T h k (·, ·). (In literature, the notation (x, y) → E X | S h kh = x, I h kh = y is also used). Then (16) implies, for all t = kh, k ∈ Nh and x, y ∈ R + ,

Proof. We can suppose y ≥ 1. We can rewrite (14) as

Thanks to (34), we have for all h, x, y ∈ R + and k ∈ N:

Then the first inequality of the statement follows from the fact that the eigenvalues of the matrices e −γh 1 0 e −λh are < 1. The second inequality in the statement follows from e −kγ y e −kλy x ≤ forecast (k) λ,h (x, y).

Now we can state the following result: In other words, for all t ∈ Nh and r ∈ N, one has

More generally, for t = kh and t = rh, r ≤ k, we have

Moreover for all k ∈ N * , the value E (E(I kh ), E(S kh )) has a unique inverse by the function forecast λ,h , given by:

where the Yi's are independent Bernoulli variables with law 1 − e −γh , and the Ui(x)'s are independent

Bernoulli variables with law 1 − e −λxh . It appears that (x, y) → E I h t | T k (x, y) is dependent on t − t but not on the actual value of t and t1.

Finally, (37) can be deduced from (35) by E(I kh ) E(S kh ) = N −1 λ,h (E(I kh )) forecast λ,h (E(I kh ), E(S kh )) = N −1 λ,h (E(I kh )) E S (k+1)h E I (k+1)h .

Proposition 5 can be used alternatively to (29) for forecasting. The latter method can presumably be transposed to more complex models, with memory effects for example (i.e. non Markovian models), when obtaining a differential equation would seem difficult to obtain. Let us state the following corollary: (18) is verified, the function ON (·) can be chosen such that ON (·) N 2 converges pointwise towards 0 when N → +∞.

Proof. Set t = nh for some n ∈ N and h ∈ R + * . We have to proof that E(S h t I h t ) = E S h t E I h t + ON (h). This follows from Proposition 5 that gives

It results from (18) that ON (·) can be chosen such that ON (·) N 2 converges pointwise towards 0 when N → +∞.

We have also the following corollary, saying that the expectations of (I h t ) and (S h t ) do not change much by rescaling h, as long as h is not too large. 

where the function O(·) is independent of t and k.

Proof. First let us prove the rescaling formula for 0 ≤ x, y ≤ N :

where the function O(·) is independent of k. It is obtained using the fact that l ∈ R + , forecast λ,l (x, y) = N λ,l (x) x y = x − λyxl + xO l 2 y + y(λx − γ)l + yO l 2 ,

which leads, provided that h > 0 is sufficiently small, for all k ∈ N and r ≤ k, to

≤ forecast Thus for all k ∈ N,

for some function O(·) independent of k, x and y, which reinserted into (42) give by a straightforward recurrence argument

Let us just write the main argument of this recurrence. For r ≤ k, (42) with l = h k , together with (43), leads to

Letting r = k − 1, we obtain (44). This formula by comparison with (42) implies (41). Hence, thanks to Proposition 5,

for all t ∈ Nh, that prove the first equality in the statement. The second equality in the statement comes by exactly the same way.

One consequence of this corollary is the following:

Corollary 6. ∀t ∈ Nh, the limits lim 

Hence lim sup l→0 + E(I l t ) = lim inf h→0 + E I h t .

Thus the limit superior in (45) can be replaced by a limit. Moreover, the same relations holds for E(S l t ) instead of E(I l t ).

Estimating the parameters of susceptible-infected-recovered model of COVID-19 cases in India during lockdown periods

Dynamical patterns of epidemic outbreaks in complex heterogeneous networks

The hidden geometry of complex, network-driven contagion phenomena

Epidemic model on a network: Analysis and applications to covid-19

Predictability in a highly stochastic system: final size of measles epidemics in small populations

A novel sub-epidemic modeling framework for short-term forecasting epidemic waves

Covid-19 singapore

Some gronwall type inequalities and applications

Concentration inequalities and martingale inequalities: a survey

The mathematics of infectious diseases

A stochastic agent-based model of the sars-cov-2 epidemic in france

Constructive proofs of concentration bounds

Multiple epidemic wave model of the covid-19 pandemic (preprint)

Projecting the transmission dynamics of sars-cov-2 through the postpandemic period

Probability Theory: A Comprehensive Course

Probabilistic Properties of Deterministic Systems

First-wave covid-19 transmissibility and severity in china outside hubei after control measures, and second-wave scenario planning: a modelling impact assessment

Fractional Dynamics on Networks and Lattices

Five proofs of chernoff's bound with applications

A sir model on a refining spatial grid i: Law of large numbers

Network-based prediction of the 2019-ncov epidemic outbreak in the chinese province hubei

On the levy-walk nature of human mobility

Monte Carlo Statistical Methods (Springer Texts in Statistics)

Simulationfree estimation of an individual-based seir model for evaluating nonpharmaceutical interventions with an application to covid-19 in the district of columbia

Forecasting models for coronavirus disease (covid-19): A survey of the state-of-the-art

Conditional measures and applications (pure and applied mathematics: A series of monographs and textbooks/177) marcel dekker, inc

Epidemic dynamics via wavelet theory and machine learning with applications to covid-19

Sir dynamics in random networks with heterogeneous connectivity

A mean field limit

A Course in Probability

The threshold of a stochastic sirs epidemic model with saturated incidence

We will use the following result: Proof. Using the same notations as in (12) , we can compute:Recalling the independence of the Bernoulli's variables Xi and Yi according to the σ-algebra generated by both I h t and S h t , we see that the last term is equal toThe conclusion follows from the boundedness of S h t and I h t by N from above.One consequence is the following corollary. It shows that under condition 18 the number of susceptible individuals and infectious individuals become less correlated when N is large and h is small.Proof. The law of total covariance together with Proposition 4 give, for h > 0 sufficiently small,The conclusion follows by a recurrence argument.Let us denote a partial order relation on R 2 asλ,h (x, y) the k-th iterate of the function forecast λ,h (x, y). We have: Lemma 4. For h > 0, there exist a constant C h > 0 such that for all x, y ∈ R + and k ∈ N * , (18) is verified, for h > 0, there exist a constant C h > 0 independent of N such that for all x, y ∈ R + , y ≤ N ,

Proof. Let us prove (36) for the second coordinate. For convenience of the reader, we will not indicate the indices λ, h, supposed to be fixed, of the function forecast λ,h . The proof for the first coordinate proceeds exactly in the same way. First note that for t, t0, t1 ∈ Nh with t ≥ t0 ≥ t1, we haveIndeed, by construction, according to (12) Let g λ,h (x, y) := y 1 − e −λxh . Let the recurrence hypothesis beand Qn ⇐⇒ ∀t ∈ N * h, ∀l ∈ 1, n − 1 , ∀k ∈ 0, n − l ,P1 is true by (17) . Let us show that Pn =⇒ Qn. Suppose Pn. Recall that P0 : R 2 → R and P1 : R 2 → R denote the functions that returns the first and second coordinate respectively.Now let us show that Qn =⇒ Pn+1. Suppose Pn. Then for t = kh,In summary, we have shown by recurrence that Pn is true for all n ∈ N * . In particular,It remains to prove that Pn is still true for k = 0, i.e. that ∀n ∈ N * , ∀n ≤ k, P0forecast (n) (I h 0 , S h 0 ) = E I h nh | T0 = E I h nh . (The last equality follows from the fact that I0 is constant, as initial condition, so T0 is trivial.) For this, it suffices to prove that the function (t, t1) → E I h t | T t 1 h (x) is only depending on the difference t − t1 and not on the actual value of t and t1. Then it follow, for t = nh, n ∈ N,The fact that the function x → E I h t | T t 1 h (x) is depending only on t − t1 follows from the fact that the sequence I h kh , S h kh k∈N is a homogeneous Markov chain, but let us prove it directly. Let t = kh, t = k h, with k, k ∈ N such that k > k . Then (12) implies ∀x, y ∈ N,