key: cord-0584857-h6ivpmmm authors: Gairat, Alexander; Shcherbakov, Vadim title: Discrete SIR model on a homogeneous tree and its continuous limit date: 2020-05-10 journal: nan DOI: nan sha: 47e0cd057c39dde10e4d0c304ae253374fd4e5de doc_id: 584857 cord_uid: h6ivpmmm This paper concerns a stochastic Susceptible-Infected-Recovered (SIR) model for the spread of infectious disease on a homogeneous tree. Under general assumptions on infection rates and recovery times we derive the distribution of the time it takes for a susceptible individual to get infected in terms of a solution of a non-linear integral equation. Then we study the discrete model in the limit when the vertex degree of the tree tends to infinity, and the infection rates are appropriately scaled. We show that in this limit our result for the discrete model implies an equation for the susceptible population compartment. This is a master equation in the sense that both the infectious and the recovered compartments can be explicitly expressed in terms of its solution. In other words, the master equation implies a continuous SIR model for the joint time evolution of all three population compartments. The obtained continuous SIR model can be naturally interpreted in terms of interaction between compartments and provides a flexible technical framework for modeling the spread of infectious disease in a homogeneous population. Studying the spread of infectious disease has been of great interest for a long time and has motivated a lot of mathematical models. Susceptible-Infected-Recovered (SIR) type models are among those that are commonly studied. A SIR model is a compartmental model, in which a population of individuals is divided into three distinct groups (compartments). The first compartment consists of individuals that are susceptible to the disease, but are not yet infected. The second compartment is represented by infected individuals. Finally, the remaining third compartment is a group of individuals, who have been infected and recovered from the disease. A recovered individual is immune. In this paper we study a continuous time discrete stochastic SIR model and its continuous limit (to be explained). In a discrete model SIR a population is modeled by vertices of a graph, and infection is transmitted from infected vertices to susceptible ones via edges of the graph. Such a model is usually studied under additional assumptions on the graph, infection rates and recovery times (e.g., see [1] , [3] , [6] , [11] , [12] , [14] , [16] , and references therein). For example, if the underlying graph is complete, infection rates are constant and the recovery times are exponentially distributed, then the model is an immediate discrete version of the classic SIR model of A.G. McKendrick and W.O. Kermack ( [9] ). We consider the discrete SIR model for the spread of infectious disease in an infinite closed homogeneous population, in which all individuals have the same number of social contacts. Such a population is modeled by a homogeneous (regular) tree. The latter is an infinite connected constant vertex degree graph without cycles. The constant vertex degree means that each vertex has the same number of adjacent vertices (neighbors). We assume that an infected vertex emits germs towards a susceptible adjacent vertex (a susceptible neighbor) according to a Poisson process with a time-varying rate. A susceptible vertex can be also infected by itself, according to another Poisson process. An infected vertex recovers in a certain period of time (the recovery time) given by a random variable. The recovery times are assumed to be independent and identically distributed. All Poisson processes under consideration are independent of each other, and they are also independent of the recovery times. Our main result for the discrete model concerns the distribution of the time it takes for a susceptible vertex of a homogeneous tree to get infected (the time to infection). Under general assumptions for infection rates and recovery times we obtain a simple analytical expression for this distribution in terms of a solution of a non-linear integral equation. In some special cases, the integral equation is equivalent to the Bernoulli type differential equation, which can be solved analytically. The structure of a homogeneous tree plays the essential role in our analysis of the discrete model. The key observation is that a susceptible vertex splits the homogeneous tree into a finite number of identical subgraphs, and infection processes on these subgraphs are independent and identically distributed. It should be noted that despite the fact that discrete SIR models on graphs have received considerable attention over the years, the SIR model on trees has been overlooked (at least, to the best of our knowledge). In the second part of the paper we study the discrete model in the limit, as the vertex degree of the tree tends to infinity, and the infection rates decrease proportionally. We show that in this limit our results for the discrete model imply an equation for the susceptible compartment. We call it the master equation, as both the infectious and the recovered compartments can be explicitly expressed in terms of its solution. In other words, the master equation implies a continuous SIR model for the joint time evolution of the three population compartments. In a sense, this generalizes the result of [8] for the Kermack-McKendrick model to the case of SIR models with time-dependent infection rates and general recovery times. Briefly, in [8] they derived a non-linear differential equation equivalent to the system of equations of the classic SIR model (see Remarks 11 and 12 below for more detail), which allowed them to obtain an analytical solution for the model equations in an exact parametric form. It should be noted that in [8] the master equation was obtained analytically from equations of the SIR model. In contrast, the master equation in the present paper is obtained from the discrete stochastic model. This is in the spirit of [2] , where a master equation for the infectious compartment was obtained from an underlying stochastic process. In fact, our master equation implies a family of continuous SIR models. A specific type of the SIR model depends on the structure and interpretation of the parameters of the master equation. This provides a flexible technical framework for modeling memory effects and arbitrarily distributed recovery times. The rest of the paper is organized as follows. In Section 2 we consider the discrete SIR model. The model is formally defined in Section 2.1. The main result concerning the distribution of the time to infection is stated and proved in Section 2.2. In Section 3.1 we use this result to derive the master equation for the susceptible compartment in the limit, as the tree vertex degree tends to infinity. Continuous SIR models implied by the master equation are considered in Section 3.2. We briefly comment of the relationship of our model with fractional SIR models in Section 4. Finally, in Section 5 we consider some special cases of the discrete model, which might be of interest in their own right. 2 The discrete SIR model In this section we formally define the discrete SIR model on a general graph, state and prove the main result in the case when the underlying graph is given by a homogeneous tree. We start with defining the discrete continuous time SIR model on a general graph (since a particular structure of the graph is not important at the definition). Let T be a connected (and possibly infinite) graph. With some abuse of notation, we will associate a graph with the set of its vertices. Given vertices x, y ∈ T we write x ∼ y, if these vertices are connected by an edge, in which case we call them neighbors. Given a vertex its vertex degree is defined as the number of its neighbors. A vertex can be either susceptible, or infected, or recovered (and immune). Initially, i.e. at time t = 0, each vertex is either susceptible, or infected (in which case we assume that it gets infected at time t = 0). When a vertex becomes infected, it starts emitting infectious germs towards a given neighbor and continues to do so until the moment of its recovery from the disease. Germs are emitted according to a Poisson process with a time dependent rate. Namely, if a vertex y ∈ T becomes infected at time t y and recovers in time H y , then, at time t ∈ [t y , t y + H y ] it infects a susceptible neighbor with the rate ε t−ty , where (ε t , t ≥ 0) is a non-negative deterministic function. In the general case the recovery times {H y , y ∈ T } are given by i.i.d. random variables with an arbitrary common distribution. A susceptible vertex can be also infected by itself, according to a Poisson process with the time-dependent rate λ t , where (λ t , t ≥ 0) is a non-negative deterministic function, so that at time t a susceptible vertex x is infected with the total infection rate We assume that all Poisson processes are independent of each other, and they are also independent of the recovery times. Proposition 1 below follows from the model definition. Proposition 1. Let ϕ t be the conditional probability that a susceptible vertex is not infected until time t by its given infectious neighbor infected at time t = 0, and let f t be the probability that a susceptible vertex is not infected until time t by itself. Then where H is a random variable, which has the same distribution as the recovery times. Remark 1. If ε t ≡ const, the recovery time is exponentially distributed, the underlying graph is complete (i.e. any two vertices are neighbors) and λ t ≡ 0, then the corresponding discrete SIR model is a discrete version of the classic continuous SIR model. Note that the case, when the recovery time is deterministic and equal to a given constant H, can be modeled by assuming that ε t = 0 for t ≥ H. Setting formally H = ∞ gives the model, in which an infected individual never recovers and stays infectious forever, which, of course, is not entirely realistic. However, by considering a function ε, which decays to zero sufficiently fast, one can model a situation, when an infected individual becomes less contagious to others over time. Remark 3. Our main result (Theorem 1 below) concerns the discrete SIR model on a homogeneous tree with the vertex degree N , i.e. where each vertex has N neighbors. However, we also consider the model on a rooted homogeneous tree with the vertex degree N ≥ 2 (in the proof of Theorem 1 below). A rooted tree is a tree, where one of the vertices, called the root, has N − 1 neighbors, while any other vertex has N neighbors. We assume throughout that all random variables are realized on a certain probability space (Ω, F, P), and the expectation with respect to the probability P is denoted by E. We also assume that all integrals and derivatives under consideration exist. In this section we state and prove the main result (Theorem 1 below) for the discrete SIR model on a homogeneous tree. In this theorem we obtain the distribution of the time it takes for a susceptible vertex to get infected in terms of a solution of an integral equation. Theorem 1. Let T be a homogeneous tree with the vertex degree n + 1, where n ≥ 1. Assume that at time t = 0 a vertex x ∈ T is infected with probability p and is susceptible with probability 1 − p independently of other vertices. Let τ be the time it takes for a susceptible vertex to get infected. Then where the function s t satisfies the integral equation where functions ϕ and f are defined in (2) and ϕ is the derivative of ϕ. Proof of Theorem 1. Consider a susceptible vertex x ∈ T and define the probability which does not depend on a neighbor due to homogeneity of both the tree and the initial condition. Recall that infectious neighbors infect the vertex x independently of each other, and also independently on self-infection. Therefore, we have that where the factor f t is defined in (2). Next, consider the probability s t in more detail. First, observe that removing the vertex x and all edges connecting x to its neighbors generates n + 1 subgraphs given by rooted trees with roots y 1 , ..., y n+1 , that are neighbors of x. Further, given a neighbor y ∼ x consider an auxiliary SIR model on the rooted tree T y with the root y. Assume that the auxiliary SIR model is specified by the same parameters as the original SIR model on the tree T . In addition, assume that at time t = 0 any vertex in the auxiliary model is either infected with probability p or susceptible with probability 1 − p, independently of other vertices. Then, it follows then from the law of total probability and the definition of the function ϕ, that where ν is the distribution of the time to infection of the root vertex y in the auxiliary SIR model on the graph T y , and ϕ is the function defined in (2) . Using similarity of the rooted trees and independence, we get, analogously to (6), the following equation Integrating by parts in (7) and using that dν(t) = −(1 − p) (f t s n t ) dt, we obtain that Recalling that ϕ t−u = 1 for u > t, s u → 0, as u → ∞, and f 0 = s 0 = 1, we get that and, hence, the integral equation (4), as claimed. Remark 5. Consider the SIR model on a homogeneous tree, in which recovery times are given i.i.d. random variables (including the degenerate case of deterministic recovery time). Recall the function (ϕ t , t ≥ 0) defined in (2) and definẽ It is easy to see that, as long as one is interested in the distribution of the time to infection, they can consider an equivalent model with just susceptible and infected compartments, in which the recovery mechanism is somehow embedded into the new infection rate given by the functionε t . For example, consider a SIR model with the infection rate (ε t t ≥ 0) such that ε t > 0 for all t ≥, and the deterministic recovery time given by a constant H > 0, theñ Trivially, if H = ∞, thenε t = ε t for all t ≥ 0. However, the functionε can differ significantly from the original function ε in the case of the random recovery time (e.g., see Corollary 3 in Section 5). 3 Continuous limit of the discrete model In this section we analyze integral equation (4) in the limit, as the tree vertex degree goes to infinity. Specifically, we show that in this limit the integral equation implies an equation for the susceptible population proportion. Consider the discrete SIR model on the homogeneous tree T with the vertex degree n + 1. Suppose that an infected vertex infects a susceptible neighbor with the rate 1 n+1 ε t after time t of being infected, where (ε t , t ≥ 0) is a non-negative function. In addition, suppose that the other model parameters (i.e. the recovery times and the rate of self-infection) do not depend on n. Let S t,n be the expected susceptible population in this SIR model. Then, any limit point S t of the sequence of functions (S t,n , t ≥ 0), n ≥ 1, must satisfy the following equation where Proof. Note first, that the corresponding ϕ-function (see equation (2)) is given by where S 0 = S 0,n = 1 − p, and the function s t,n satisfies the equation Observe that for sufficiently large n. Therefore, It is easy to see that for sufficiently large n, which implies that where the function γ is defined in (12) . In addition, note that (S 0 f u ) 1 n+1 ∼ 1 and (S u,n ) n n+1 ∼ S u,n , as n → ∞. Therefore, we can rewrite (15) as follows Recalling that log(f t ) = − t 0 λ u du, we obtain that which implies equation (11) for any limit point S t , as claimed. Differentiating (11) gives the master equation in the differential form which, integrating by parts, can be rewritten as follows Remark 6. The existence and the uniqueness of solution of equation (11) follows from general results for integral equations with delay ( [10] ). It can be shown that, under mild assumptions, the sequence of functions (S t,n , t ≥ 0), n ≥ 1, is equicontinuous. Therefore, there exists a subsequence that uniformly converges to the solution of (11). We skip the technical details. Remark 7. It follows from equation 11 that stationary value S ∞ satisfies the following equation Remark 8. Note that equation (11) (or its differential equivalent (20)) is a standalone equation for the susceptible population S t , namely that this equation does not involve neither the infected, nor the recovered populations. where ε t is the rate of infection in the discrete model. In particular, if H = ∞, then γ t = ε t . In general, these two functions can be significantly different (see Example 2 below). In this section we show that the master equation (11) argument is reinforced. Note that for simplicity of exposition and without loss of generality we assume throughout this section that λ t ≡ 0, i.e. there is no self-infection. The basic continuous SIR model implied by the master equation is a two-compartmental model, in which the population is divided into two compartments, namely, the compartment of susceptible individuals, described by the variable S t , and the compartment of infected ones, described by the variable I t , so that 1 = S t + I t and for t ≥ 0. Then S t = −I t , which allows to rewrite equation (20) as follows Equation (24) Example 1 (The model with latent period). Suppose that an infected individual is latent for a non-random period of time of length L > 0. In addition, suppose that the rate of infection is constant. Then γ t = ε1 {t≥L} , where ε is the rate of infection, and equation (24) becomes as follows S t = −εS t I t−L and I t = −S t . Suppose that the function γ t is of the following form where ε > 0 is a given constant and β t is a non-increasing positive function, such that β 0 = 1 and β t → 0, as t → ∞. Then, the master equation implies the continuous SIR model with the three standard compartments, in which an infected individual recovers in a time given by random variable ξ with the tail distribution P(ξ > t) = β t , and during its infectious period it infects any susceptible one with the constant rate ε. Indeed, under these assumptions, equation (20) is as follows (recall that (22)) Define It is easy to see that Moreover, one can show that both I t ≥ 0 and R t ≥ 0 (we skip the details). Therefore, variables I t and R t can be interpreted as the population proportions of infected and recovered individuals respectively in the continuous SIR model with the constant rate of infection ε and the random recovery time with the tail distribution given by the function β. Indeed, in this model the infected compartment at time t consists of 1) those who were infected at time 0 and did not recover before time t (which gives the first term I 0 β t in (27)), and 2) those, who were infected at time u ∈ (0, t] and did not recover before time t (integrating over time gives the integral term in (27)). A similar argument gives (40), which also follows from (27) and (29) combined with the initial condition S 0 + I 0 = 1. Differentiating both (27) and (28), and combining them with (26), we get the following system of equations Remark 10. The system of equations (30)-(32) is similar to the system of equations of the delay model proposed in [5] for modeling the spread of Covid-19 in Italy. Example 2 (The classic SIR model). Consider the discrete SIR model on a homogeneous tree, in which an infected individual infects its susceptible neighbor with the constant rate ε and recovers with the constant rate µ. In addition, assume that there is no self-infection, i.e. λ t = 0 for t ≥ 0. The corresponding function α is given by where H is a random variable exponentially distributed with parameter µ. A direct computation gives that Since λ t ≡ 0 and β t = −µe −µt for t ≥ 0, the system of equations (30)-(32) becomes as follows which is the system of equations of the classic SIR model (with the infection rate ε and the recovery rate µ), where S t , I t and R t stand for population proportions of susceptible, infected and recovered individuals respectively. Remark 11. Note that in Example 2 it is probably more transparent to start with equation (11) , which in this case is Then, differentiating (33) gives that Combining (34) with (33) we obtain the following differential form of the master equation in the classic case Then, setting I t = 1 − S t + µ ε log St S 0 , one can proceed as in Example 2. Remark 12. It should be noted that our equation (35) is equivalent to equation (26) in [8] . In [8] they analytically obtained the equation from equations of the classic SIR model and used it to derive an exact analytical solution of the SIR model in a parametric form. Remark 13. Note that equating the time derivative to zero in equation (35) , i.e. S t = 0, gives the known equation for the stationary population proportion of susceptible individuals S := S ∞ in the SIR model (e.g. see equation (7) in [4] and references therein). In particular, this equation shows that the limit value S depends only on the ratio µ/ε = 1/R 0 , where R 0 is the basic reproduction number. Note also that equation (36) is just a special case of more general equation (21) for the stationary susceptible state. Example 3 (Constant rate of infection and deterministic recovery). Consider a model, in which the infection rate is given by a constant ε > 0, and the recovery time is given by a deterministic constant H > 0. This model can be obtained by setting γ t = ε for t ∈ [0, H] and γ t = 0 for t > H. This gives the following model equations where S t = I t = 0 for t < 0. In this section we generalize SIR models considered in Sections 3.2.1 and 3.2.2. Suppose that the function γ is of the following form where (w t , t ≥ 0) is a non-negative function and the function (β t , t ≥ 0) is the tail distribution of some positive random variable (i.e. similarly to what we assumed in Section 3.2.2). Arguing as in Section 3.2.2, we obtain the continuous SIR model described by the following equations where, S t , I t and R t are population proportions of susceptible, infected and recovered individuals respectively, so that 1 = S t + I t + R t for t ≥ 0. As before, we assumed that R 0 = 0. In this model an infected individual recovers in a random time given by a random variable ξ with the tail distribution P(ξ > t) = β t and, if it is infected at time u, then it infects any susceptible one with the rate w t at the time u + t. Recall, that we also assume (22). In the differential form the model equations are as follows Remark 14. By choosing appropriate functions (w t , t ≥ 0) and (β t , t ≥ 0) one can model various infection rates and recovery distributions. For example, using the power law functions allows to model memory effects observed in real data (e.g. see [2] and references therein). In the rest of this section we use the idea from [2] in order to rewrite equations (41)-(43) in terms of a certain kernel. The idea is based on the fact that these equations contain convolutions, which makes it possible to apply the Laplace transform. Let (L{g} t , t ∈ R + ) be the Laplace transform of a function (g t , t ∈ R + ). It follows from (42) that and, hence, Thus, for any appropriate function (g t , t ∈ R + ) we have that Since L −1 (L{a}L{b}) is equal to the convolution a * b, we can rewrite (44) as follows where K is a kernel defined by Finally, using (45) with g t = −γ t in (41), and with g t = −β t in (42) and (43) gives the system of the model equations in the kernel form Remark 15. Note that equation (46) is an analogue of equation (16) in [2] . One of the recognized drawbacks of the classic SIR model is that both the infection rate and the recovery rate do not depend on the state of the system, i.e. the model is memoryless. A popular approach to modeling memory effects consists in using fractional SIR models (e.g., see [15] and references therein). Some of these models are obtained by formal replacement of ordinary derivatives by fractional derivatives of a certain type. This gives a system of fractional differential equations that is equivalent to a system of integro-differential equations with a power-law kernel. For example, replacing ordinary derivatives in the classic SIR model by Caputo fractional derivatives gives a system of fractional differential equations, which are equivalent to the following system of integro-differential equations with the kernel K y = y α−2 Γ(α−1) , where α ∈ (0, 1] and Γ is the Gamma-function. However, we are not quite clear what physical/biological process is described by equations (50)-(52). In contrast, equations (41)-(43) and their equivalents in the kernel form, i.e. equations (47)-(49), can be naturally interpreted in terms of the interaction between compartments. For example, in the absence of any external factors and self-infection, the flux into the infected compartment is equal to the flux out of the susceptible compartment. The susceptible compartment decreases at the rate proportional to its current value S t . The value K(γ) t−u (in (47)-(49)) describes the impact made on the susceptible compartment at time t by those individuals, who were infected earlier and is still infectious. The coefficient of proportionality, i.e. the integral term t 0 I u K(γ) t−u du, measures the total impact of the infected compartment on the susceptible one over the time period [0, t]. This generalizes the interaction between susceptible and infected compartments in the classic SIR model, where only the current value I t is taken into account. It should be also noted that the continuous SIR model in the present paper is obtained by passing to the limit in the discrete stochastic SIR model. This is in line with SIR models in the kernel form that are derived from stochastic processes based on natural biological assumptions (e.g., see [2] , [5] and references therein). In this section, we consider some special cases of the discrete SIR model on the homogeneous tree with the vertex degree n + 1. In these cases the integral equation (4) can be rewritten in an equivalent differential form, which is of interest on its own right. For simplicity of notations we assume that all vertices are initially susceptible (i.e. p = 0 in Theorem 1). s t = ε(n − 1) + λ ε(n − 1)e −λt + λe ε(n−1)t so that P(τ > t) = e −λt e − 2ε λ (e −λt −1+λt) , if n = 1; (54) P(τ > t) = e −λt ε(n − 1) + λ ε(n − 1)e −λt + λe ε(n−1)t n+1 n−1 , if n ≥ 2. (55) Proof of Corollary 1. Since ε t = ε1 {t≥0} and λ t = λ1 {t≥0} , we have that ϕ t = 1, t < 0, e −εt , t ≥ 0, and f t = 1, t < 0, and ϕ t = −εe −εt for t ≥ 0. The integral equation (4) in this case is as follows Differentiating (57) and simplifying gives the following differential equation Proof of Corollary 2. In this case we have that f t = e −λt and ϕ t = e −ε min(t,H) for t ≥ 0, and equation (4) becomes as follows s t = ϕ t + ε t 0 e −λu s n t ϕ t−u du for t < H, ϕ t + ε t t−H e −λu s n u ϕ t−u u for t ≥ H. A direct computation gives that We are going to consider an example of the model with a random recovery time. Corollary 3 (Exponential recovery time). Suppose that the recovery time is exponentially distributed with the parameter µ and functions (ε t , t ∈ R + ) and (λ t , t ∈ R + ) are as in Corollaries 1 and 2. Then integral equation (4) is equivalent to the following differential equation Recall Remark 5 and definẽ ε t := − (log(ϕ t )) = ε(µ + ε) µe (µ+ε)t + ε 1 {t≥0} for t ≥ 0, so that ϕ t = −ε t ϕ t for t ≥ 0. A direct computation gives that ε t −ε 2 t = −(µ + ε) γ t for t > 0. Using the equation in the preceding display and differentiating equation (4) gives that s t = −ε t ϕ(t) + εf t s n t + t 0 f (u)s n u ϕ t−u ε t−u − ε 2 t−u du = −ε t ϕ t + εf t s n t − (µ + ε) t 0 f u s n u ϕ t−uεt−u du = −ε t ϕ(t) + εf t s n t − (µ + ε)(s t − ϕ t ) = −(µ + ε)s t + εf t s n t + (µ + ε −ε t )ϕ t . Noting that gives the equation s t = −(µ+ε)s t +εf t s n t +µ, which is the Bernoulli equation with the additional term µ, as claimed. Stochastic epidemic models and their statistical analysis A Fractional Order Recovery SIR Model from a Stochastic Process Stochastic Epidemic Models with Inference Accurate closed-form solution of the SIR epidemic model Solvable delay model for epidemic spreading: the case of Covid-19 in Italy Exploring the threshold of epidemic spreading for a stochastic SIR model with local and global contacts Contact process without revival on tree Exact analytical solutions of the Susceptible-Infected-Recovered (SIR) epidemic model and of the SIR model with equal death and birth rates A contribution to the mathematical theory of epidemics Global existence and uniqueness of solutions of integral equations with delay: progressive contractions A stochastic SIR model on a graph with epidemiological and population dynamics occurring over the same time scale Epidemic outbreaks in complex heterogeneous networks Who solved the Bernoulli differential equation and how did they do it? Exact solution of a stochastic susceptible-infectious-recovered model Review of fractional epidemic models Fractal scale-free networks resistant to disease spread It is easy to verify that the function defined in (53) is a solution of equation 58.Remark 16. The equation (58) is the well known Bernoulli equation (e.g., see [13] ). Under assumptions of Corollary 1 the model was originally considered in [7] .Remark 17. It follows from equation (55) thatfor sufficiently large n. Further, if ε = ε n , where ε n n → c > 0, as n → ∞, thenIn other words, the probability P(τ ≤ t) converges to a linear transformation of the logistic curve 1 1 + c λ e −(c+λ)t . Corollary 2. Suppose that the recovery time is given by a deterministic constant H > 0, λ t = λ1 {t≥0} and ε t = ε1 {0≤t≤H} , where ε > 0 and λ > 0 are given constants. Then integral equation (4) is equivalent to the following differential equation