Extending the basic SIR Model in R
Accounting for population turnover, vaccination and waning immunity
Compartmental models are a modelling approach where the population is divided into different ‘compartments’, representing their disease status and possibly other characteristics. Here, mathematical equations are used to model transitions between different compartments.
The SIR Model
The SIR model forms the cornerstone for infectious disease modelling.It’s a compartmental model where the individuals belong to one of the three compartments:
- Susceptible (S)
- Infected (I)
- Recovered (R)
It’s a population-based model in the sense that each compartment models the behavior of an average individual in that compartment.Of course, it’s an oversimplification of the real-life scenario but it is actually possible to build more complex models based on this.
What are the assumptions of an SIR model?
- Homogeneous population : Individuals in the same compartment is subject to the same hazards
- Well-mixed population : All susceptible have the same risk of getting infected
- Immunity forever : Individuals who recovered from the disease are immune forever
Th last assumption leads to the observation that SIR model is good for scenarios where population doesn’t change and the immunity is stable during the simulation period.
These assumptions give us the following set of differential equations:
where S,I and R represents the number of susceptible,infected and recovered individuals in each compartment respectively and N=S+I+R.
γ is the recovery rate . β is the infection rate i.e the average number of secondary infections per unit time.
Basic Reproduction Number(R₀)
Th basic reproduction number is the average number of secondary infections caused by a single infected case, in a fully susceptible population.
For an SIR model, it can be calculated as :
Effective Reproduction Number
If the population of susceptible varies with time, it’s preferable to use the effective reproduction number. It is the average number of secondary cases arising from an infected case at a given point in the epidemic.
Force of infection
It is the risk of infection of an individual, per unit time. In simple terms,it is the force acting on susceptible individuals to change them into infected.It changes with time.At the beginning of the epidemic, since there are not many infected people ,the force of infection will be small but as the epidemic grows, so does the force of infection.
It can be expressed as :
Now, we can substitute this expression in the equations considered above to obtain :
It’s possible to consider two cases for this model :
(a) SIR model with constant λ
(b) SIR model with varying λ
In the following cases, we assume that at the start of the epidemic, there is only one infected individual and no one has recovered yet.
Extensions of the basic model
The extension considered here are :
- Population turnover
- Vaccination
- Waning immunity
These factors are called modifiers of susceptibility in the population.
Modelling Population Turnover
In many cases, the people who gained immunity will be replaced by newborns who are fully susceptible.So, it will be not be realistic to consider the simple SIR model. It becomes necessary to account for births and deaths in the population. For simplicity, the disease induced mortality is ignored.
The assumptions are:
- Every individual in each compartment experiences the same background mortality
- b is the rate at which newborns are entering the susceptible compartment
The extension will include birth rate b and mortality rate μ which gives the following differential equations:
Modelling Vaccination Effects
Let p be the proportion of population vaccinated with a perfectly effective vaccine.This will be reflecting in the initial conditions as there will be a proportion of individuals in the recovered compartment.
The equations will be given as :
with initial conditions at t=0 as S=(1-p)(N-1), I=1 and R=p(N-1)
Critical Vaccination Threshold
Not everyone needs to be vaccinated to prevent the epidemic.This is because of the phenomenon called herd immunity. Herd immunity is said to be achieved when there is sufficient immunity in a population to interrupt transmission.
So the proportion of population that needs to be vaccinated(p) can be calculated as
Modelling Waning Immunity
It may happen that people gain immunity after recovering from the infection but the immunity doesn’t last forever. So, these individuals become susceptible again. Here, σ is the waning rate.
SIR Models in R
The deSolve package in R contains functions to solve initial value problems of a system of first-order ordinary differential equations (‘ODE’). The ggplot2 package provides functions for visualizations.
Simple SIR Model
Initially we consider a simple SIR model with varying force of infection(λ).
Assume that the population size is 1000000 and at the start of the epidemic, there is only infected individual .Also, the infectious person infects one person on average every 2 days and is infectious for 5 days.
So here, β=1/2 day⁻¹=0.5 day⁻¹ and γ=1/5 day⁻¹=0.2 day⁻¹
The time period of interest is 100 days in daily intervals.
The infection peaks at around 48 days.
SIR model extended to account for population turnover
Assuming the population size stays constant over time i.e at 10⁶,so b=μ.Background mortality rate is calculated based on average human lifespan, say 70 years. So, μ=1/70 year⁻¹ and b=1/70 year⁻¹.Also,β=0.4*365 year⁻¹ and γ=0.2*365 year⁻¹.
The prevalence plot for susceptible has peaks with long deep troughs in between.We also observe that the infection peaks when the susceptible pool reaches it’s maximum.
Impact of vaccination
Suppose that 40% of the population is vaccinated against the epidemic.Assume γ=0.1 day⁻¹ and β=0.4 day⁻¹.The model is run for a time period of 3 years.
We observe from the plot above that the epidemic still occurs infecting around 12.5% of the population even with 40% vaccinated.
The critical vaccination threshold comes out to be 0.75 in this case. Hence, 75% of the population needs to be vaccinated to prevent the epidemic.
We observe from the plot above that the epidemic does not occur.
Waning immunity
Assume γ=0.2 day⁻¹, β=0.4 day⁻¹ and waning rate, σ=1/10 year⁻¹ if average immunity period is taken as 10 years. The model is run for a time period of 50 years in daily intervals.
The prevalence plot appears similar to the situation with population turnover. While in the population turnover scenario, the source of new susceptibles were births, in this case, the source turns out to be individuals losing their immunity.
Combining all these mechanisms, it can help address the question raised by policy-makers regarding the proportion of people to be vaccinated, the frequency of vaccination and the like. More complex models are needed for these real world situations but the underlying dynamics remain the same.
In the next post, we will see how we can model the effect of interventions in an SIR model :)
References
- Soetaert, K. E., Petzoldt, T., & Setzer, R. W. (2010). Solving differential equations in R: package deSolve. Journal of statistical software, 33.
- https://www.coursera.org/specializations/infectious-disease-modelling
