Stochastic differential equations: Application to epidemiology

A Dynamical System for Epidemiology

By Nadhir Hassen in Julia ODE Dynamical Systems

June 1, 2021

Dynamical Systems with Machine Learning in Engineering

Introduction

A SIR model is an epidemiological model that calculates the theoretical number of people infected with a contagious disease in a closed population over time. The name of this class of models derives from the fact that they involve coupled equations connecting the number of susceptible people S (t), the number of infected people I (t) and the number of people who have recovered R (t). One of the simpler SIR models is the Kermack-McKendrick model .

The Kermack-McKendrick model is a SIR model for the number of people infected with a contagious disease in a closed population over time. It has been proposed to explain the rapid increase and decrease in the number of infected patients observed in epidemics such as plague (London 1665-1666, Bombay 1906) and cholera (London 1865). It assumes that the population size is fixed (i.e. No birth, death from disease, or death from natural causes), the incubation period of the infectious agent is instantaneous, and the incubation period of the infectious agent is instantaneous. duration of infectivity is the same as duration of illness. It also assumes a completely homogeneous population with no age, spatial or social structure.

Deterministic SIR model

The model consists of a system of the following ordinary differential equations

dSdt=−βNSI dIdt=βNSI−γI dRdt=γI,

where β and γ is the infection birth and death rate, in other words we can interpret β as the transmission rate and γ the recovery rate of the disease, in addition we denote the effective reproduction factor (the rate of contamination) R=βS(0)γN.
Assume that S(0),I(0),R(0)≥0 and S(0)+I(0)+R(0)=N. R represents the average number of secondary infections produced by an infected individual (susceptible to infection) during the incubation period. At t=0, we assume that β=γ. We can notice that R=R0S(0)/N. We can obtain that R≤1, then I(t) is decreasing monotonically and therefore the epidemic is non-existent. If R>1, then I(t) starts by increasing and ends up decreasing until reaching zero: we are therefore in the presence of an epidemic. In addition, an important notion to take into consideration is the severity of the epidemic, this is the total number of cases or the final size of the epidemic noted as R(∞). If we assume that R(0)=0 and I(0)=1 we can therefore obtain the value of R(∞) from the differential equations dIdt and dSdt

dIdS=−if(β/N)−γIif(β/N) =γNβS−1.

The solution can be obtained by the method of separations of variables and by using the initial conditions S(0)=N−1 and I(0)=1, we obtain I(t)+S(t)=γβNlnS(t)+−γβNln(N−1)+N,
when t→∞ and I(∞)=0 implies S(∞)=N(γβln(S(∞)N−1)+1).
This result gives us R(∞)=N−S(∞).

Stochastic SIR model

Let S(t),I(t) and R(t) be the random variables which correspond to the number of individuals likely to be infected, individuals infected and individuals immunized respectively. We therefore have S(t)+I(t)+R(t)=N. In this stochastic model there is no latent period, ie infected people are also infectious. As S(t) and I(t) are independent then the transition probabilities are

pΔS(t)=i,ΔI(t)=j|S(t),I(t)={βNS(t)I(t)Δt+o(Δt),if (i,j)=(−1,1), γI(t)Δt+o(Δt),if (i,j)=(0,−1), 1−re,if (i,j)=(0,0), o(Δt)otherwise

where re=(βNS(t)I(t)+γI(t))Δt+o(Δt). In the case when ΔI(t)=−1 then ΔR(t)=1. Consider that the initial conditions are S(0)=s0≥0 and I(0)=i0≥0 we have pi,j(t)=p(S(t)=i),I(t)=j. To find the evolution of the following Markov chain we use the Kolmogorov equations such that dp(i,j)(t)dt=βN(i+1)(j−1)p(i+1,j−1)(t)+γ(j+1)p(i,j+1)(t)−(βNij+γj)p(i,j)(t),

where i=0,1,2,…,N, j=0,1,2,…,N−1 and i+j≤N. So when S(0)=Nj≈N and I(0)=j small enough, we have R=βγS(0)N≈βγ=R0. We are therefore in the presence of a process of pure birth and death. Indeed, the death rate corresponds to the reimission rate μ=γ and the birth rate to the rate of new infections λ≈β. At the start of the pandemic i.e. when I(0)≈j is small enough, the likelihood of an epidemic ending very quickly (or avoiding an epidemic) could be considered a process. of pure birth and death. So the probability of not having an epidemic can be summarized as {1if R≤1 (1R0)j if R0>1.

The SEIR model- a variant of the original model

To make the model more realistic, we modify the model, the states are now S→E→E→I→R. We always assume that the size of the population remains fixed S(t)+E(t)+I(t)+R(t)=N for all t. We redefine the transition rates between states as being

  • β(t) the transmission rate or effective contact rate (the rate at which individuals go from S to E).
  • σ(t) the infection rate (the rate at which individuals go from exposed to infected) \item γ the recovery rate (the rate at which individuals go from I to R and stay in R assuming they die or create immunity)

As in the previous section, we represent the model by a system of the following ordinary differential equations dSdt=−βNSI dEdt=βNSI−σE dIdt=σE−γI dRdt=γI.

As transmissions depend on time, we define R0:=β(t)/γ, this implies that

  • The transitions from the state I to R are done at a Poisson rate γ, the average time spent in the infected state I is therefore 1/γ.
  • Prolonged interactions occur at the rate β, so that a new individual entering the infected state will potentially transmit the virus to other individuals with an average of R0=βγ.

To solve this problem we can make the following reparametrization β(t)=γR0(t), moreover we define s:=S/N as being the proportion of individuals in each state, by dividing each equation by N we obtain the following system dSdt=−γR0sidedt=γR0si−σedidt=σe−γidRdt=γi.

We can therefore express the proportion of individuals restored by r=1−s−e−i due to the fact that the states form a partition. By solving this system of differential equations, we can visualize the dynamics of the model in the first figure by considering the case of province of Quebec on November 20,2020 as well as the proportions corresponding to this population in second.

2021 PSB Poster (letter)

The SEIR model with confinement policy

So far we have assumed that R0(t) depends only on the infection rate β(t) and recovery rate γ(t). A containment measure can disturb the rate of contamination R0(t). We therefore introduce an additional constraint dR0dt=ν(R0−R0), the parameter ν controls the speed of dispersion from R0 to the reference value R0 following the containment policy (social distancing, closure of shops for example). Another element so far considered constant is the mortality rate, releasing this assumption to better adapt the model to reality, then the variation in the mortality rate is expressed by ddtd(t)=δγi. Thus by adapting the present model to reality, we fix σ=1/5 the incubation rate of the disease, γ=1/14 which corresponds to 14 days according to the ministry of public health, the reproduction rate of the reference disease R0=R0=1.5 and δ=0.004 the mortality rate in Quebec. In the figure of cumulative cases confirms that after a containment policy, the rate of contamination delays the cases of infections (the blue line). 2021 PSB Poster (letter)

2021 PSB Poster (letter)

The SEIR model with random shocks

We now want to incorporate several random effects that affect the rate of transmission R0(t) from shocs through the process β(t), this can be explained by the fact of have different behavior of the population coming from different sources of information (social networks, presidential elections, protests for individual freedom, etc.), but also shocks can come from different regional policies (red zone, orange zone, green zone) which are not subject to the same confinement restrictions.

To reflect this randomness, we can incorporate into the variation of the contamination rate R0(t) a stochastic volatility σ√R0. To ensure that the process is a stationary process in the weak sense we have to add √R0. Then we define the instantaneous variation of the contamination rate such that

dRt=ν(ˉR0(t)−R0(t))dt+σ√R0(t)dWt

where W is a standard Brownian motion (or Weiner process). To get even closer to reality, we must also incorporate a random effect in the death rate, indeed not all individuals react in the same way to infection, this can be caused by age, regional area, state of health, etc. Let δ(t) be the mortality rate, we make some assumptions below

  • δt represents a mortality reference, which would be the average of the process.
  • a diffusion term with α√δ(1−δ) volatility, with 0≤δ≤1.

From what follows, we define the death rate as a stochastic process according to dδt=θ(δ−δt)dt+α√δt(1−δt)dWt

with θ is the parameter which controls the dispersion between the current mortality rate at time t and the reference mortality rate to be maintained. The general form of the stochastic differential equation becomes dSdt=−γR0sididt=γR0si−γidrdt=(1−δ)γidddt=δγidR0dt=ν(R0(t)−R0)dδdt=θ(δ−δ).
2021 PSB Poster (letter)

The solution of this stochastic differential equation leads to several possible realizations due to the random effect of shocks (more precisely to the white noise processes incorporated in our modeling). These two realization of the solution which represents two possible scenarios to control the pandemic. In the trajectory figure we can clearly see that with a containment measure, the mortality rate reaches a peak and stabilizes, the same observation is made for the recovery rate. Another important element to consider is the effect of containment on controlling the rate of contamination R0 which greatly depends on the rate of infections.