The SIR epidemic model from a PDE point of view
aa r X i v : . [ q - b i o . P E ] D ec The SIR epidemic model from a PDE point of view
Fabio A. C. C. Chalub
Departamento de Matem´atica and Centro de Matem´atica e Aplica¸c˜oes, UniversidadeNova de Lisboa, Quinta da Torre, 2829-516, Caparica, Portugal.
Max O. Souza
Departamento de Matem´atica Aplicada, Universidade Federal Fluminense, R. M´arioSantos Braga, s/n, 22240-920, Niter´oi, RJ, Brasil.
Abstract
We present a derivation of the classical SIR model through a mean-fieldapproximation from a discrete version of SIR. We then obtain a hyperbolicforward Kolmogorov equation, and show that its projected characteristicsrecover the standard SIR model. Moreover, we show that the long time limitof the evolution will be a Dirac measure. The exact position will depend onthe well-know R parameter, and it will be supported on the correspondingstable SIR equilibrium.
1. Introduction
A very fruitful modeling paradigm in epidemiology is the so-called com-partmental models, with dynamics governed by mass-action laws. Most clas-sical epidemiological models are of this type, and this has led to a numberof both quantitative and qualitative predictions in the disease dynamics [2].More recently, there is a growing interest in discrete, agent-based, models [7].See also [6] for a comparison between different models. In many cases, thesemodels are thought to be more realistic, and able to capture important dy-namical features that are not present in the continuous models.
Email addresses: [email protected] (Fabio A. C. C. Chalub), [email protected] (Max O. Souza)
November 9, 2018 ere, we follow the ideas in [3], to study the the large-population regimeof the discrete dynamics. In this way, we obtain a Hyperbolic Forward Kol-mogorov equation for the probability density evolution. It is well known,through the method of characteristics, that there is a strong linkage betweensolutions to first-order Partial Differential Equations (PDEs) and systems ofOrdinary Differential Equations (ODEs). Thus, its not entirely surprisingthat projected characteristics of this PDE will be related to the classical SIR(Susceptible-Infected-Removed) ODE model.This modeling through PDEs has some advantages, in particular, it al-lows the introduction of higher-order effects, like, for example, stochasticity,adding a second-order term to the equation. On the other hand, if we con-sider a discrete model in population dynamics and consider its limit of largepopulation (under suitable condition) we naturally obtain a PDE for p ( t, x ),the probability density to find the population at state x at time t . The ODEcan then be obtained as the hyperbolic limit of the PDE, or, alternativelyas the initial dynamics of the PDE. In short, this means that the dynam-ics of the discrete population can be approximated for short times and largepopulation by a certain ODE — the derivation of this ODE requires an intro-duction of intermediate models, a stochastic differential equation or a partialdifferential one.The ODE approach can be seen in [1], while the PDE modeling was thesubject of a previous work from the authors, where the replicator equationwas obtained as the limit of the finite-population discrete Moran process [3].The resulting equation is of singular type and required a specific analysis ofits behavior [4].In this work, we will study in a certain level of detail the SIR epidemicmodel from the PDE point of view . This is one of the most elementary andwell studied model in mathematical epidemiology. See [5, 2].In particular, we shall prove that the solution of the SIR-hyperbolic-PDE obtained as a first order expansion in the inverse of the populationsize from the discrete-SIR converges when t → ∞ to a Dirac-delta measuresupported at the unique stable equilibrium of the SIR-ODE. Moreover (andthis will be proved in a forthcoming work) the solution of the SIR-hyperbolic The expression SIR appears in the literature in two different context: one as a discreteevolutionary system, used in general in computer simulations; the second as an ODEsystem. We hope that all these different meanings are clear from the context.
2s approximated by the SIR-ODE in all time scales and approximate thediscrete-SIR for short time scales. This provides a framework to unify allthese descriptions. In a forthcoming work, we will also introduce the SIR-parabolic-PDE which has the inverse behavior (approximate the discrete-SIRfor all time scales and the ODE-SIR for short times).
2. Discrete and continuous SIR models
Consider a discrete SIR model, i.e., consider a fixed size population of N individuals, each one in one of the three states: n individuals S usceptible(i.e, individuals with no imunity), m individuals I nfected (i.e., individualcurrenctly infected by a given infectious disease and able to transmit it tothe susceptibles) and N − n − m R emoved (after infection idividuals have atemporary imunity and then are removed from the dynamics; after certaintime they become susceptible again). This is a very simple model for nonlethal diseases transmitted by contact, e.g., normal influenza. At each timestep of size ∆ t > • If it is S , then it changes to I with probability proportional to thefraction of I in the remainder, αm/ ( N − • If it is I , then it changes to R with constant probability β ; • If it is R , then it changes to S with constant probability γ .This can be summarized in the following diagram: S + I α −→ I + I , I β −→ R , R γ −→ S . This model is also called SIRS, and when γ = 0 we recover the classicalSIR model. For simplicity, however, we will call it the SIR model in thiswork and all results presented here include the case γ = 0.Constants α , β and γ depend, in principle, in N and ∆ t . As we areinterested in the limit behavior when N → ∞ , ∆ t → scaling relationslim N →∞ , ∆ t → αN ∆ t = a , lim N →∞ , ∆ t → βN ∆ t = b , lim N →∞ , ∆ t → γN ∆ t = c , with ab = 0. For further informations on scalings, see [3].3et P ( N, ∆ t ) ( t, n, m ) be the probability that at time t we have n susceptible, m infected and N − n − m removed, where the total population N is constantand the time step is given by ∆ t >
0. Therefore P ( N, ∆ t ) ( t + ∆ t, n, m ) = α ( n + 1)( m − N ( N − P ( N, ∆ t ) ( t, n + 1 , m − β m + 1 N P ( N, ∆ t ) ( t, n, m + 1) + γ N − n − m + 1 N P ( N, ∆ t ) ( t, n − , m )+ (cid:20) nN (cid:18) − α mN − (cid:19) + mN (1 − β ) + N − n − mN (1 − γ ) (cid:21) P ( N, ∆ t ) ( t, n, m ) . Now, define x = n/N , y = m/N and p ( t, x, y ) = P ( t, xN, yN ; N ). Then,using p ( t, x, y ) = p and keeping terms until order 1 /N : p ( t + ∆ t, x, y ) = α (cid:0) x + N (cid:1) (cid:0) y − N (cid:1)(cid:0) − N (cid:1) p (cid:18) t, x + 1 N , y − N (cid:19) + β (cid:18) y + 1 N (cid:19) p (cid:18) t, x, y + 1 N (cid:19) + γ (cid:18) − x − y + 1 N (cid:19) p (cid:18) t, x − N , y (cid:19) + (cid:18) x (cid:18) − αy − N (cid:19) + y (1 − β ) + (1 − x − y )(1 − γ ) (cid:19) p ( t, x, y ) ≈ p + 1 N [( α ( y − x ) + β + γ ) p + ( αxy − γ (1 − x − y )) ∂ x p + ( βy − αxy ) ∂ y p ]= p + 1 N [ ∂ x (( αxy − γ (1 − x − y )) p ) + ∂ y (( β − αx ) yp )]Finally, ∂ t p = ∂ x (( axy − c (1 − x − y )) p ) + ∂ y (( b − ax ) yp ) . (1)subject to probability conservation, i.e,dd t Z p ( t, x )d x = 0 . (2)
3. The SIR model as a transport problem
Let x = ( x, y ), and let Φ t ( x ) be the flow map associated to the SIRsystem ˙ X = c (1 − X − Y ) − aXY , ˙ Y = ( aX − b ) Y , p ( x ) is C ( R ) function, and p ( t, x, y ) = e Q ( x ) − Q (Φ − t ( x )) p (Φ − t ( x )),with x = ( x, y ), and Q satisfies F · ∇ Q = −∇ · F, where F denotes the right hand side of the SIR system. Let S denote theunit simplex in R .Fix x ∈ S . Then, we have thate − Q ( φ t ( x )) p ( t, φ t ( x )) = e Q ( x ) p ( x ) . Hence,0 = e Q ( φ t ( x )) dd t (cid:2) e − Q ( φ t ( x )) p ( t, φ t ( x )) (cid:3) = − F ( φ t ( x )) · ∇ Q ( φ t ( x ) p ( t, φ t ( x )) + ∂ t p ( t, φ t ( x ) + F ( φ t ( x )) ∇ p ( t, φ t ( x ))= ∇ · F p + F ∇ p + ∂ t p = ∂ t p + ∇ · ( pF ) . Thus, the SIR system are the characteristics of (1), and the probabilitydensity should be transported along them. We now make this calculationmore precise.
Definition 1.
Let F : U ⊂ R n → R n , be a Lipschitz vector field, where U is an open set, and let Ω ⊂ U be compact. We say that Ω is regularlyattracting for F , if there is an open set V with a piecewise smooth boundaryand Ω ⊂ V ⊂ ¯ V ⊂ U , such that, if Φ t denotes the flow by F restricted to V ,then we have that ω ( V ) ⊂ ¯Ω . Theorem 1.
Let Ω be a domain with piecewise smooth boundary ∂ Ω . Let F : U ⊂ R N → R N be Lipschitz, with Ω ⊂ U being a regularly attracting setfor F . Let p ∈ L (Ω) be nonnegative. Then the equation ∂ t p + ∇ · ( pF ) = 0 , p (0 , x ) = p ( x ) (3) has a unique solution p ( t, x ) = e Q ( x ) − Q (Φ − t ( x )) p (Φ − t ( x )) Moreover, p is nonnegative, and supp( p ( t, · )) ⊂ Ω .In addition, if supp( p ) ⊂ Ω , then dd t Z Ω p ( t, x )d x = 0 . (4)5 roof. Existence can be shown as follows by considering the weak formula-tion, Let W = [0 , ∞ ) × Ω, and let ψ ∈ C c ( W ). Z Z W p ( t, x ) ∂ t ψ ( t, x )d x d t + Z Z W p ( t, x ) ∇ ψ ( t, x )d x d t ++ Z Ω p (0 , x ) ψ (0 , x )d x d t = 0 . (5)Choose ψ ∈ C c ((0 , ∞ ) × Ω). Then (5) becomes
Z Z W p ( t, x ) ∂ t ψ ( t, x )d x d t + Z Z W p ( t, x ) ∇ ψ ( t, x )d x d t = 0 . Let y = φ t ( x ) and η t ( x ) = det( ∂ x Φ t ( x )). Also let W t = Φ t ( W ). Then, thefirst integral becomes Z Z W t e Q (Φ t ( y )) − Q ( y ) p ( y ) ∂ t ψ ( t, Φ t ( y )) η t ( y )d y d t. The second integral becomes
Z Z W t e Q (Φ t ( y )) − Q ( y ) p ( y ) ∇ · ψ ( t, Φ t ( y )) η t ( y )d y d t. Combining both integrals, we can write
Z Z W t e Q (Φ t ( y )) − Q ( y ) p ( y ) dd t ψ ( t, Φ t ( y )) η t ( y )d y d t. On integrating by parts, we have that
Z Z W t e − Q ( y ) p ( y ) ψ ( t, Φ t ( y ))e Q (Φ t ( y )) [ F · ∇ Q − ∇ · F ] d y d t = 0 . We have that p is clearly nonnegative. Moreover, Let V = supp( u ( t, · )) andlet V t = (Φ − t ( W )). Let ψ be a vanishing function in Ω. Then0 = Z Φ t (Ω) u ( t, x ) ψ ( x )d x = Z Ω e Q (Φ t ( y ) − Q ( y ) u ( y ) ψ ( y )div( F )d y Hence supp( u ( t, · )) ⊂ Φ t (Ω). Since Ω is regularly attracting for F , we havethat Φ t (Ω) ⊂ Ω. If supp( p ( y )) ⊂ Ω, then we extend p by defining it to be6ero in R N − U . Then we have thatdd t Z Ω u ( t, · )d x = dd t Z V u ( t, · )d x = − Z V div( u ( t, · ) ˆ F )d x = 0 . Let F S ( x, y ) = ( c (1 − x − y ) − αxy, αxy − by ) . We define the reflected SIR field by F RS ( x, y ) = ( c (1 − x − y ) − axy, ( ax − b ) | y | )We immediately have Lemma 1.
Let Ω be the simplex in the nonnegative orthant of R . Then F RS is a C vector field in R , and Ω is regularly attracting for F RS . Thus for smooth solutions, we have
Theorem 2.
Consider the Cauchy problem for (1) , with a L non-negativeinitial condition p . Then there exists a unique solution satisfying (2) . More-over, p ( t, x ) ≥ .
4. Asymptotic Behavior and measure solutions
We recall that the dynamics of SIR are controlled by the parameter R := a/b (see [2]), i.e, Proposition 1.
Let x = (1 , , and x = (cid:16) R , cc + b (cid:16) − R (cid:17)(cid:17) be the twoequilibria of SIR. x is referred to as the disease free equilibrium and x asthe endemic equilibrium. If R ≤ , then any solution that starts in thenonnegative orthant of R approaches x for large time. If R > , then x is the limiting point. With this point of view we have 7 heorem 3.
Let p be a solution of (1) , satisfying (2) . Then, in the Wasser-stein metric, we have that lim t →∞ p ( t, x ) = (cid:26) δ x , R ≤ δ x , R > . Proof.
We deal with the case R ≤
1; the case R > R ≤
1, we have that x is the globally asymptotic stable equilibrium. Thengiven, δ >
0, we can find
T >
0, such that, for t > T , we have thatΦ t ( S ) ⊂ B δ ( x ) . Let ψ ( x ) be a continuous function. Then, for t > T , we have that Z S p ( t, x ) ψ ( x )d x = Z B δ ( x ) p ( t, x ) ψ ( x )d x . But, let ǫ > ψ is continuous, we have δ >
0, such that ψ ( x ) − ǫ ≤ Z B δ ( x ) p ( t, x ) ψ ( x )d x ≤ ψ ( x ) + ǫ, and this proves the claim, since we have that Z B ǫ ( x ) p ( t, Φ t ( y )d y = Z B ǫ ( x ) e Q (Φ t ( y )) − Q ( y )) p ( y ) ψ (Φ t ( y )) η t ( y )d y Acknowledgments
FACCC is partially supported by FCT/Portugal, grants PTDC/MAT/68615/2006and PTDC/MAT/66426/2006, PTDC/FIS/7093/2006. MOS is partiallysupported by FAPERJ grants 170.382/2006 and 110.174/2009. Both FACCand MOS are partially supported by the bilateral agreement Brazil-Portugal(CAPES-FCT). 8 eferences [1] Linda J. S. Allen. An introduction to stochastic epidemic models. In
Mathematical epidemiology , volume 1945 of
Lecture Notes in Math. , pages81–130. Springer, Berlin, 2008.[2] Nicholas F. Britton.
Essential mathematical biology . Springer Undergrad-uate Mathematics Series. Springer-Verlag London Ltd., London, 2003.[3] F.A.C.C Chalub and M. O. Souza. From discrete to continuous evolutionmodels: A unifying approach to drift-diffusion and replicator dynamics.
Theoretical Population Biology , 76(4):268–277, 2009.[4] F.A.C.C Chalub and M. O. Souza. A non-standard evolution problemarising in population genetics.
Commun. Math. Sci. , 7(2):489–502, 2009.[5] J. D. Murray.
Mathematical biology. I , volume 17 of
InterdisciplinaryApplied Mathematics . Springer-Verlag, New York, third edition, 2002.An introduction.[6] G. Schneckenreithera, N. Poppera, G. Zaunera, and F. Breiteneckera.Modelling sir-type epidemics by ODEs, PDEs, difference equations andcellular automata a comparative study.
Simulation Modelling Practiceand Theory , 16(8):1014–1023, 2008.[7] J. Verdasca, M. M. Telo da Gama, A Nunes, N. R. Bernardino, J. M.Pacheco, and M. C. Gomes. Recurrent epidemics in small world networks.