SIR, birth and death

SIR, birth and death#

SIR_Birth_Death()

Here we consider an SIR model in which individuals may be removed by death from each compartment at a uniform rate per person, \(\gamma\). The population is replenished via births into the susceptible compartment at the same rate, thus conserving the total population by design. For deterministic evolution, the population size remains constant whereas for stochastic evolution, the size fluctuates around this value. The equations are as follows:

\[\begin{split}\begin{aligned} \frac{\mathrm{d}S}{\mathrm{d}t} &= \mu N - \frac{\beta SI}{N} - \mu S \\ \frac{\mathrm{d}I}{\mathrm{d}t} &= \frac{\beta SI}{N} - \gamma I - \mu I \\ \frac{\mathrm{d}R}{\mathrm{d}t} &= \gamma I - \mu R \end{aligned}\end{split}\]

As an example, we study stochastic evolution of this system with measles-like parameters in 3 differently sized populations. This provides a demonstration of threshold population sizes in order to support endemic circulation of certain pathogens.

from pygom import common_models
import matplotlib.pyplot as plt
import numpy as np
import math

#####################
# Set up PyGOM object
#####################

# Parameters
n_pop=1e4
mu=0.01/365    # birth/death rate 1% per year
gamma=1/20 
R0=15
beta=R0*gamma

ode = common_models.SIR_Birth_Death({'beta':beta, 'gamma':gamma, 'mu':mu})

# Time range and increments
tmax=365*10                           # maximum time over which to run solver
dt=1                                  # timestep
n_timestep=math.ceil(tmax/dt)         # number of iterations
t = np.linspace(0, tmax, n_timestep)  # times at which solution will be evaluated

# Initial conditions (endemic equilibrium derived from stationary point)
def sir_bd_endemic_eq(mu, beta, gamma, n_pop):
    s0=math.floor((gamma+mu)*n_pop/beta)
    i0=math.floor(mu*(n_pop-s0)*n_pop/(beta*s0))
    r0=n_pop-(s0+i0)
    return [s0, i0, r0, n_pop]

x0=sir_bd_endemic_eq(mu, beta, gamma, n_pop)
ode.initial_values = (x0, t[0])

##########
# Simulate
##########
n_sim=10
np.random.seed(1)

solution, jumps, simT = ode.solve_stochast(t, n_sim, full_output=True)
y=np.dstack(solution)

############################
# try larger population size
############################
n_pop=1e5
ode = common_models.SIR_Birth_Death({'beta':beta, 'gamma':gamma, 'mu':mu})  # update parameter
x0=sir_bd_endemic_eq(mu, beta, gamma, n_pop)  # recalculate IC's
ode.initial_values = (x0, t[0])
solution_2, jumps_2, simT_2 = ode.solve_stochast(t, n_sim, full_output=True)  # simulate
y_2=np.dstack(solution_2)

#################################
# try even larger population size
#################################
n_pop=1e6
ode = common_models.SIR_Birth_Death({'beta':beta, 'gamma':gamma, 'mu':mu})  # update parameter
x0=sir_bd_endemic_eq(mu, beta, gamma, n_pop)  # recalculate IC's
ode.initial_values = (x0, t[0])
solution_3, jumps_3, simT_3 = ode.solve_stochast(t, n_sim, full_output=True)  # simulate
y_3=np.dstack(solution_3)

Show code cell output Hide code cell output

Illegal jump, x: [7.09700e+03 1.00000e+00 9.30670e+04 1.00165e+05], new x: [ 7.32000e+03 -1.00000e+00  9.28500e+04  1.00169e+05]
Illegal jump, x: [5.61400e+03 1.10000e+01 9.45900e+04 1.00215e+05], new x: [ 5.64400e+03 -2.00000e+00  9.45700e+04  1.00212e+05]
Illegal jump, x: [5.9840e+03 7.0000e+00 9.4239e+04 1.0023e+05], new x: [ 6.04100e+03 -1.00000e+00  9.41940e+04  1.00234e+05]
Illegal jump, x: [6.56600e+03 1.00000e+00 9.36870e+04 1.00254e+05], new x: [ 6.84900e+03 -2.00000e+00  9.34250e+04  1.00272e+05]
Illegal jump, x: [6.3660e+03 1.0000e+00 9.3599e+04 9.9966e+04], new x: [ 6.66900e+03 -1.00000e+00  9.33380e+04  1.00006e+05]
Illegal jump, x: [6.3670e+03 1.0000e+00 9.3599e+04 9.9967e+04], new x: [ 6.6540e+03 -2.0000e+00  9.3328e+04  9.9980e+04]
Illegal jump, x: [6.8920e+03 3.0000e+00 9.3063e+04 9.9958e+04], new x: [ 6.9610e+03 -1.0000e+00  9.2983e+04  9.9943e+04]
Illegal jump, x: [6.2920e+03 1.0000e+00 9.3601e+04 9.9894e+04], new x: [ 6.6070e+03 -1.0000e+00  9.3314e+04  9.9920e+04]
Illegal jump, x: [5.91700e+03 4.00000e+00 9.42530e+04 1.00174e+05], new x: [ 6.00000e+03 -2.00000e+00  9.41790e+04  1.00177e+05]
Illegal jump, x: [6.76200e+03 1.00000e+00 9.34790e+04 1.00242e+05], new x: [ 7.03200e+03 -2.00000e+00  9.31890e+04  1.00219e+05]
Illegal jump, x: [7.33400e+03 3.00000e+00 9.28890e+04 1.00226e+05], new x: [ 7.39300e+03 -1.00000e+00  9.28110e+04  1.00203e+05]
Illegal jump, x: [7.1960e+03 4.0000e+00 9.2763e+04 9.9963e+04], new x: [ 7.2500e+03 -3.0000e+00  9.2704e+04  9.9951e+04]
Illegal jump, x: [6.4420e+03 1.0000e+00 9.3554e+04 9.9997e+04], new x: [ 6.74300e+03 -1.00000e+00  9.32660e+04  1.00008e+05]

Illegal jump, x: [6.35540e+04 3.20000e+01 9.36019e+05 9.99605e+05], new x: [ 7.28980e+04 -5.20000e+01  9.26926e+05  9.99772e+05]
Illegal jump, x: [6.363800e+04 2.300000e+01 9.366510e+05 1.000312e+06], new x: [ 7.666500e+04 -1.300000e+01  9.237850e+05  1.000437e+06]
Illegal jump, x: [6.363900e+04 2.300000e+01 9.366510e+05 1.000313e+06], new x: [ 7.645400e+04 -3.400000e+01  9.239020e+05  1.000322e+06]
Illegal jump, x: [6.363900e+04 2.300000e+01 9.366500e+05 1.000312e+06], new x: [ 7.638700e+04 -7.000000e+00  9.239390e+05  1.000319e+06]

Plotting the results, we see that for populations of sizes 10,000 and 100,000, the infected population is critically close to zero, such that stochastic fluctuations eventually lead to disease extinction. This is of course signified by the infected class reaching zero, but also by the recovered and susceptible classes undergoing stable linear growth due to population turnover. When the population size is 1,000,000, we see that the infected subset, of typical size 500, is able to persist for the full 10 years of the simulation.

../../_images/71aa12813e784f8fa546851cae6dbf9f9fbd9b1fcf82bc2d4c547d67be9c51d5.png