Discrete time. (week by week, or here abstract into step by step)
$N$ popluation size.
Start 1 infectious, $N-1$ susceptible
Every infectious indiviual infect all the rest of population with probability $p$ during each time step and then recover. So essentially every timestep it is only the newly infected from the last timestep that are infectious.
Goes on till no new infected.

Basic reproduction numner $R_0$¶

$R_0$: average number of new infections caused bu a typical infected during early phase of an outbreek.

\begin{align*} R_0 = (N-1)*p \end{align*}

If $N = 1000$ and $p = 0.0015$ then

\begin{align*} R_0 \sim 1.4985 \end{align*}

Simulate the Model 10 000 times¶

# import packages
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import plotly.express as px
import plotly.graph_objects as go
import plotly.offline as py
py.init_notebook_mode(connected=False)

# simulation
total_infected, time_steps = [],[]

for tt in range(10000):
    R = 999*0.0015
    time_step = 0
    s = 1 # number of sick person to start
    N = 1000 # population total
    del_s = s #start value of newly infected =s
    ppl_infected = [s]
    ppl_remain = N-s
    while del_s >0 and ppl_remain >0:
        time_step = time_step + 1
        nbr_trial, prob = del_s, 0.0015
        infections = np.random.binomial(nbr_trial, prob, ppl_remain)
        del_s = np.sum(infections >=1) # newly infected
        s = s + del_s # over all ever infected
        ppl_remain = N-s
    
    total_infected.append(s)
    time_steps.append(time_step)
df_1 = pd.DataFrame({'time_steps':time_steps,
                   'total_infected':total_infected
                  })

df_1.head()

Distribution of total infected¶

font_dict = dict(
        family="Courier New, monospace",
        size=16 )
fig = px.histogram(df_1, x="total_infected", #labels = {"total_infected":"nbr infected"}, 
                   nbins=100, range_y=[0,1000])
fig.update_layout(font=font_dict)

fig=px.scatter(df_1, x='time_steps',y='total_infected', opacity=0.5 )
fig.update_traces(marker=dict(size=12))
fig.update_layout(font=font_dict)

Most of the cases, the epidemic last 20-40 timesteps (weeks) and infect 500-700 people in total, given total population is 1000 and $R_0 \sim$ 1.5.

Equation and solution for total infected population fraction $\tau$¶

1. Simple model with one type of population response to infection¶

Since $1-\tau$ should be the probability of not getting infected while exposed to $N\tau$ infectious trials, $\tau$ should satisfy the following condition:

(eq. 1)

\begin{align*} 1 - \tau &= (1-p)^{N\tau} \\ &= (1-\frac{R_0}{N})^{N\tau} \\ &= (1-\frac{1}{{N}/{R_0}})^{\frac{N}{R_0}*R_0\tau} \\ &= e^{-R_0\tau} \end{align*}

2. Model with two type of population responses to infection¶

In a slightly more mixed situation, fraction $w$ of the population could be less suceptible to the virus, and having a fraction $v$ of the original suceptibility $p$. So the equaction becomes:
(eq. 2)

\begin{align*} 1 - \tau &= (1-p)^{N(1-w)\tau}(1-vp)^{Nw\tau} \\ &= (1-\frac{R_0}{N})^{N(1-w)\tau} (1-\frac{vR_0}{N})^{Nw\tau} \\ &= (1-\frac{1}{{N}/{R_0}})^{\frac{N}{R_0}*R_0(1-w)\tau}(1-\frac{1}{(1/v)(N/R_0)})^{(1/v)(N/R_0)*vR_0w\tau} \\ &= e^{-R_0(1-w)\tau}e^{-R_0vw\tau}\\ &= e^{-R_0\tau(1-w+vw)} \end{align*}

(eq.2) becomes (eq.1) when $w = 0$.

Solving numerically $\tau$, $R_0$ in range $\tau = 0 \sim 1$ given $R_0 = 0 \sim 20$ for simple case and mixed case:¶

def eq_mix(tau, R, w, v):
    return np.abs(1-tau-np.exp(-R*tau*(1-w+w*v)))

def solve_eq(eq, w, v, tau_range = np.arange(0,1,0.001), R_range = np.arange(0, 25, 0.01)):
    tt, rr = np.meshgrid(tau_range, R_range)
    rs = np.array([R_range[i] for i in eq(tt, rr, w=w, v=v).argmin(axis = 0)])
    res = np.array([eq(tt, rr, w=w, v=v) for tt, rr in zip(tau_range, rs) ])
    tau, rs =  tau_range[res<1e-3], rs[res<1e-3]
    tau, rs = np.insert(tau,len(tau), max(tau_range)), np.insert(rs, len(rs), max(R_range))
    return tau, rs

# simple case w = 0, v = 1 the same as 
tau0, rs0 = solve_eq(eq=eq_mix, w = 0, v=1)
label0 = 'all equaly susceptible'

# mixed case w = 0.5, v = 0.5 the same as 
tau1, rs1 = solve_eq(eq=eq_mix, w = 0.5, v=0.5)
label1 = 'half the ppl. are half susceptible'

# mixed case w = 0.7, v = 0 the same as 
tau2, rs2 = solve_eq(eq=eq_mix, w = 0.7, v=0)
label2 = '70% the ppl. are vaccinated'

# mixed case w = 0.9, v = 0 the same as 
tau3, rs3 = solve_eq(eq=eq_mix, w = 0.9, v=0)
label3 = '90% the ppl. are vaccinated'

# mixed case w = 0.95, v = 0 the same as 
tau4, rs4 = solve_eq(eq=eq_mix, w = 0.95, v=0)
label4 = '95% the ppl. are vaccinated'

# mixed case w = 0.99, v = 0 the same as 
tau5, rs5 = solve_eq(eq=eq_mix, w = 0.99, v=0)
label5 = '99% the ppl. are vaccinated'

data_tuples = [(tau0, rs0, label0), (tau1, rs1, label1), (tau2, rs2, label2), (tau3, rs3, label3), 
                (tau4, rs4, label4), (tau5, rs5, label5)]
disease_tuples = [(2,3, 'Corona'), (5,7, 'Polio'), (12,18, 'Measles') ]
traces=[go.Scatter(x=rs, y=tau, name=label) for tau, rs, label in data_tuples]
fig=go.Figure(data=traces)
for x0, x1, label in disease_tuples:
    fig.add_vrect(x0=x0, x1=x1, annotation_text=label, annotation_position="top",annotation_font=dict(size=14),
                  opacity=0.3, fillcolor="LightSalmon", layer="above", line_width=0)
fig.update_layout(font=font_dict, legend_title_text='Succeptibility', showlegend = True, xaxis=dict(title='R_0'),
                  yaxis=dict(title='infected fraction'))
fig.show()

Importance of vaccination¶

It is evident that with less suceptible population the less total infection we end up with. With 70% vaccination the total population can be safe for viruses with moderate $R_0$ values such as coronavirus ($R_0 = 2\sim3$ ). If 90% of the population are vaccinated., the outbreaks of infectious diseases such as polio ($R_0 = 5\sim7$) can be prevented, but not for measles ($R_0 = 12\sim18$ ). With 99% vaccination rate, the outbreaks of very infectious measles can be safely prevented too.

HIT (herd immunity threshold)¶

(ref: https://en.wikipedia.org/wiki/Herd_immunity) When a critial proportion of the population is immune, the disease cannot persist in the population. This proportion is called HIT (herd immunity threshold). The formula for HIT is:

(eq. 3)

\begin{align*} HIT = 1-\frac{1}{R_0} \end{align*}

For measles it is 92 $\sim$ 95%, for polio it is 80 $\sim$ 86%.

Peak infection as influenced by $R_0$¶

$R_0$ can be factorized:

\begin{align*} R_0 = pcl \end{align*}

where

$p$: transmission probability of a contact.
$c$: number of c "contacs" per day
$l$: duration of infectious period

Prevention measures that decreases $R_0$:

$p$: face mask, hand-washing, condom ...
$c$: quarantine, avoid crowds, avoid public transports ...
$l$: quicker diagnosis, isolation

Given 100 timesteps (weeks) get peak infections graph by simulating 1000 times.

n_simulations = 1000
n_timesteps = 100
infection_cases = np.zeros(n_timesteps)

dfs = []
for R in (0.5, 0.75, 1, 1.25, 1.5, 1.75, 2, 2.25):
    for tt in range(n_simulations):
        s = 1 # number of sick person to start
        N = 1000 # population total
        del_s = s #start value of newly infected =s
        ppl_infected = np.zeros(n_timesteps)
        ppl_infected[0] = s
        for timestep in np.arange(n_timesteps-1):
            if del_s >0 and ppl_remain >0:
                ppl_remain = N-s
                nbr_trial, prob = del_s, R/N
                infections = np.random.binomial(nbr_trial, prob, ppl_remain)
                del_s = np.sum(infections >=1) # newly infected
                ppl_infected[timestep+1] = del_s
                s = int(ppl_infected.sum())
        infection_cases = infection_cases+ppl_infected

    infection_cases = infection_cases/n_simulations
    dfs.append(pd.DataFrame({f'R = {R}':infection_cases}))
df = pd.concat(dfs, axis = 1)

px.line(df, range_x=[0,40])

With $R_0$ reduced, the peak can be much smaller and much delayed. So prevention measures such as masks, hand-washing, quarantine, isolation can reduce the load on the healthcare systems.

from IPython.display import Javascript
from nbconvert import HTMLExporter

def save_notebook():
    display(
        Javascript("IPython.notebook.save_notebook()"),
        include=['application/javascript']
    )

def output_HTML(read_file, output_file):
    import codecs
    import nbformat
    exporter = HTMLExporter()
    # read_file is '.ipynb', output_file is '.html'
    output_notebook = nbformat.read(read_file, as_version=4)
    output, resources = exporter.from_notebook_node(output_notebook)
    codecs.open(output_file, 'w', encoding='utf-8').write(output)

def save_as_html(notebook_title):    
    import time
    # save_notebook()
    time.sleep(3)
    current_file = f'{notebook_title}.ipynb'
    output_file = f'{notebook_title}.html'
    output_HTML(current_file, output_file)
save_as_html('Simple Epidemic Model')

	time_steps	total_infected
0	24	578
1	32	578
2	25	505
3	36	558
4	31	489

Table of Contents

Intro¶

Background¶

Assumptions:¶

Model¶