While the 2019-nCoV is spreading crazily around the globe, I’d like to investigate into the classical epidemic models that academia uses to simulate and analyze these diseases.$\newcommand{P}{\text{P}}\newcommand{d}{\text{d}}$ The most famous ones are referred to as compartmental models where the total population is divided into several compartments with assumptions upon transition rates between these Markov states. The model was first introduced in the early 20th century in literature like Kermack and McKendrick (1927).

Background Knowledge

Before we start walking through the models, let’s first introduce some basic yet necessary knowledge.

Notations

The compartmental model, as the name indicates, involve four compartments:

• $S_t$: susceptible people (i.e. healthy people)
• $E_t$: exposed people (who’re infected but not yet infectious)
• $I_t$: infectious people
• $R_t$: recovered people (who’re immune to the disease)
• $N_t$: total population

We have

$$ N_t = S_t + I_t + R_t $$

and $N_t\equiv N$ within a short period of time (assuming no vital population dynamics).

In the meantime, we have parameters in these models as below:

• $\beta$: infectious rate, which represents the probability of transmitting disease between a susceptible and an infectious individual.
• $\gamma=1/D$: recovery rate, which is determined by the average duration $D$ of infection.
• $\sigma$: incubation rate, which is the rate of latent individuals becoming infectious (average duration of incubation is $1/\sigma$).
• $\xi$: this represents the rate which recovered individuals return to the susceptible statue due to loss of immunity
• $\mu$: birth rate
• $\nu$: death rate (among infectious people)

Solution to a Special ODE

Claim. The solution to ODE

$$ \frac{\d x_t}{\d t} + a_1 x_t + a_2 x_t^2 = 0\quad (a_1\neq 0, a_2\neq 0) $$

with initial condition $x_0$ given, is

$$ x_t = \frac{\exp(-a_1 t)}{1/x_0 - [\exp(-a_1 t) - 1]a_2/a_1}. $$

Proof. We have

$$ \begin{align*} &\frac{\d x_t}{\d t} + a_1 x_t + a_2 x_t^2 = 0\ \Rightarrow& -\frac{\d}{\d t}\left(\frac{1}{x_t}\right) + a_1 \left(\frac{1}{x_t}\right) + a_2 = 0 \end{align*} $$

which, by letting $y_t:=1/x_t$, is

$$ \begin{align*} &\frac{\d y_t}{\d t} - a_1 y_t = a_2\ \Rightarrow& \frac{\d }{\d t}[\exp(-a_1 t) y] = a_2 \exp(-a_1 t). \end{align*} $$

Integrate both sides and

$$ \exp(-a_1 t) y_t - y_0 = [1-\exp(-a_1 t)]a_2/a_1. $$

Therefore, by substituting $x_t=1/y_t$ and $x_0=1/y_0$ we have

$$ x_t = \frac{\exp(-a_1 t)}{1/x_0 + [1 - \exp(-a_1 t)]a_2/a_1} $$

which is exactly the formula given in the claim.Q.E.D.

SI Model

SI model is the simplest among its siblings: there’re only susceptible and infectious people in our problem and they together form a closed-loop Markov chain. $\newcommand{P}{\text{P}}\newcommand{d}{\text{d}}$When there are no vital population dynamics, we have

SI w/o Vital Dynamics

$$ \frac{\d I_t}{\d t} = -\frac{\d S_t}{\d t} = \frac{\beta S_tI_t}{N_t} $$

where $N_t = S_t + I_t\equiv N$ remains a constant. We have closed-form solution as follows

$$ \begin{cases} \displaystyle{I_t = \frac{N_tI_0}{I_0 + S_0\exp(-\beta t)}}\\ S_t = N_t - I_t \end{cases} $$

which indicates $S_t\to 0$ with time, i.e. under this setting the total population becomes infected eventually.

SI w/ Vital Dynamics

With significant birth or death happening, we need to adjust our model s.t. the two rates, $\mu$ and $\nu$, are taken into consideration. The ODE then becomes

$$ \begin{cases} \displaystyle{\frac{\d S_t}{\d t}=\mu N_t - \nu S_t - \frac{\beta S_t I_t}{N_t}}\\ \displaystyle{\frac{\d I_t}{\d t}=\frac{\beta S_t I_t}{N_t} - \nu I_t} \end{cases} $$

where $N_t=S_t+I_t$ is not a constant any more. Instead of a closed-form solution, here we’re more interest in the stationary distribution. When

$$ \mu=\nu=\frac{\beta S_t}{N_t} = \beta\left(1 - \frac{I_t}{N_t}\right) $$

we know the final epidemic size is (in percentage)

$$ \lim_{t\to\infty}\frac{I_t}{N_t} = \left(1 - \frac{\nu}{\beta}\right)_{!!+}. $$

SIS Model

Notice in the SI model we didn’t consider the possibility that a infectious person can recover and become susceptible again. By assuming so and taking the recovery rate $\gamma$ into our model, we have

SIS w/o Vital Dynamics

$$ \frac{\d I_t}{\d t} = -\frac{\d S_t}{\d t} = \frac{\beta S_tI_t}{N_t} - \gamma I_t $$

where $N_t=S_t+I_t\equiv N$ is a constant. The solution is

$$ \begin{cases} \displaystyle{I_t = \frac{\varphi N_tI_0}{[\varphi N_t - I_0\beta]\exp(-\varphi t)+1}}\\ S_t = N_t - I_t \end{cases} $$

where $\varphi=\beta-\gamma$. The asymptotic epidemic size is

$$ \lim_{t\to\infty} \frac{I_t}{N_t} = \left(1 - \frac{\gamma}{\beta}\right)_{!!+}. $$

SIS w/ Vital Dynamics

With birth and death, the SIS model becomes

$$ \begin{cases} \displaystyle{\frac{\d S_t}{\d t}=\mu N_t - \nu S_t + \gamma I_t - \frac{\beta S_t I_t}{N_t}}\\ \displaystyle{\frac{\d I_t}{\d t}=\frac{\beta S_t I_t}{N_t} - \nu I_t - \gamma I_t} \end{cases} $$

and by solving $\d I_t/\d t=0$ we have the asymptotic epidemic size

$$ \lim_{t\to\infty}\frac{I_t}{N_t} = \left(1 - \frac{\gamma+\nu}{\beta}\right)_{!!+}. $$

SIR Model

Now we consider a third state of people: recovered. These people have recovered from the disease and will never be susceptible again thanks to immunity.

SIR w/o Vital Dynamics

First we assume no vital dynamics. The model, where $N_t=S_t+I_t+R_t\equiv N$, now becomes

$$ \begin{cases} \displaystyle{\frac{\d S_t}{\d t} = -\frac{\beta S_t I_t}{N_t}}\\ \displaystyle{\frac{\d I_t}{\d t} = \frac{\beta S_t I_t}{N_t} - \gamma I_t}\\ \displaystyle{\frac{\d R_t}{\d t} = \gamma I_t} \end{cases} $$

which have no closed-form solution but still we can find its asymptotics. It can be proved that $I_t\to 0$ with time (i.e. all infectious people will recover eventually) and the epidemic size attains its maximum at $(1-\gamma/\beta)_+$, which coincides with the limit from the SIS model.

SIR w/ Vital Dynamics

When including potential vital dynamics into the SIR model, the ODE becomes

$$ \begin{cases} \displaystyle{\frac{\d S_t}{\d t} = \mu N_t - \nu S_t -\frac{\beta S_t I_t}{N_t}}\\ \displaystyle{\frac{\d I_t}{\d t} = \frac{\beta S_t I_t}{N_t} - \nu I_t - \gamma I_t}\\ \displaystyle{\frac{\d R_t}{\d t} = \gamma I_t - \nu R_t} \end{cases} $$

which, just like without the dynamics, has no known closed-form solution. Moreover, since new births can bring new susceptible individuals to the population, the disease doesn’t die out like in the previous case. Instead, it will attain its maximum at $[1-(\gamma+\nu)/\beta]_+$ and eventually reach a stationary state at $[\nu/(\nu+\gamma)-\nu/\beta]_+$.

SIRS Model

The SIR model assumes people carry lifelong immunity to a disease upon recovery; this is the case for a variety of diseases. For another class of airborne diseases, for example seasonal influenza, an individual’s immunity may wane over time. In this case, the SIRS model is used to allow recovered individuals to return to a susceptible state.

SIRS w/o Vital Dynamics

The model is

$$ \begin{cases} \displaystyle{\frac{\d S_t}{\d t} = - \frac{\beta S_t I_t}{N_t} + \xi R_t}\\ \displaystyle{\frac{\d I_t}{\d t} = \frac{\beta S_t I_t }{N_t} - \gamma I_t}\\ \displaystyle{\frac{\d R_t}{\d t} = \gamma I_t - \xi R_t} \end{cases} $$

where $N_t = S_t+ I_t + R_t \equiv N$ is a constant. If there is sufficient influx to the susceptible population, at equilibrium the dynamics will be in an endemic state with damped oscillation due to people losing immunity after a period of time and becoming susceptible again.

SIRS w/ Vital Dynamics

We can also add vital dynamics into the SIRS model. The ODE becomes

$$ \begin{cases} \displaystyle{\frac{\d S_t}{\d t} = \mu N_t - \nu S_t - \frac{\beta S_t I_t}{N_t} + \xi R_t}\\ \displaystyle{\frac{\d I_t}{\d t} = \frac{\beta S_t I_t }{N_t} - \nu I_t - \gamma I_t}\\ \displaystyle{\frac{\d R_t}{\d t} = \gamma I_t - \xi R_t - \nu R_t} \end{cases} $$

where $N_t=S_t+I_t+R_t$.

The stationary distributions of the SIRS models, despite the oscillation, can be derived in the exactly same way as we discussed above.

SEIR Model

Finally let’s introduce the last epidemic state, i.e. being exposed. In this category of epidemic models, individuals experience a long incubation period when they’re infected but not yet infectious to others.

SEIR w/o Vital Dynamics

Without considering birth or death, we may assume model as follows

$$ \begin{cases} \displaystyle{\frac{\d S_t}{\d t} = -\frac{\beta S_t I_t}{N_t}}\\ \displaystyle{\frac{\d E_t}{\d t} = \frac{\beta S_t I_t}{N_t}}-\sigma E_t\\ \displaystyle{\frac{\d I_t}{\d t} = \sigma E_t - \gamma I_t}\\ \displaystyle{\frac{\d R_t}{\d t} = \gamma I_t} \end{cases} $$

where $N_t = S_t + E_t + I_t + R_t\equiv N$ is a constant.

SEIR w/ Vital Dynamics

Similar to SIR with vital dynamics, the SEIR model gives an oscillated convergence towards a steady state should we allow for birth and death to be significant:

$$ \begin{cases} \displaystyle{\frac{\d S_t}{\d t} = \mu N_t - \nu S_t -\frac{\beta S_t I_t}{N_t}}\\ \displaystyle{\frac{\d E_t}{\d t} = \frac{\beta S_t I_t}{N_t}}-\nu E_t\\ \displaystyle{\frac{\d I_t}{\d t} = \sigma E_t - \gamma I_t - \nu I_t}\\ \displaystyle{\frac{\d R_t}{\d t} = \gamma I_t - \nu R_t} \end{cases} $$

where $N_t=S_t + E_t + I_t + R_t$.

SEIRS Model

In the SIER model we assume people carry lifelong immunity to a disease, which in many cases is not true. By introducing an outflow from R to S we have SEIRS model.

SEIRS w/o Vital Dynamics

The ODE is

$$ \begin{cases} \displaystyle{\frac{\d S_t}{\d t} = \xi R_t-\frac{\beta S_t I_t}{N_t}}\\ \displaystyle{\frac{\d E_t}{\d t} = \frac{\beta S_t I_t}{N_t}}-\sigma E_t\\ \displaystyle{\frac{\d I_t}{\d t} = \sigma E_t - \gamma I_t}\\ \displaystyle{\frac{\d R_t}{\d t} = \gamma I_t-\xi R_t} \end{cases} $$

where $N_t = S_t + E_t + I_t + R_t\equiv N$ is a constant.

SEIRS w/ Vital Dynamics

The ODE is

$$ \begin{cases} \displaystyle{\frac{\d S_t}{\d t} = \mu N_t - \nu S_t + \xi R_t-\frac{\beta S_t I_t}{N_t}}\\ \displaystyle{\frac{\d E_t}{\d t} = \frac{\beta S_t I_t}{N_t}}-\nu E_t\\ \displaystyle{\frac{\d I_t}{\d t} = \sigma E_t - \gamma I_t - \nu I_t}\\ \displaystyle{\frac{\d R_t}{\d t} = \gamma I_t - \nu R_t - \xi R_t} \end{cases} $$

where $N_t=S_t + E_t + I_t + R_t$.

Visualization

Below is a simple visualization of the SEIRS model which you can interactive with. Specifically, here I also include the fatality rate $v$ due to the underlying infectious disease. By doing so, the final model we simulate here becomes a so-called deadly SEIRS model:

$$ \begin{cases} \displaystyle{\frac{\d S_t}{\d t} = \mu N_t - \nu S_t + \xi R_t-\frac{\beta S_t I_t}{N_t}}\\ \displaystyle{\frac{\d E_t}{\d t} = \frac{\beta S_t I_t}{N_t}}-\nu E_t\\ \displaystyle{\frac{\d I_t}{\d t} = \sigma E_t - \gamma I_t - \nu I_t - v I_t}\\ \displaystyle{\frac{\d R_t}{\d t} = \gamma I_t - \nu R_t - \xi R_t} \end{cases} $$

with $N_t=S_t + E_t + I_t + R_t$. With tiny perturbation in certain parameters below, you may find dramatic difference in the epidemic course within the population. If you need to reset the plot together with all sliders, simply double-click anywhere on the figure.

# import packages
import pyperclip
import numpy as np
from bokeh.embed import components
from bokeh.events import DoubleTap
from bokeh.io import show, output_notebook
from bokeh.layouts import gridplot, column
from bokeh.plotting import figure, output_file, show, ColumnDataSource
from bokeh.models import CustomJS, Slider, Legend, Range1d, NumeralTickFormatter


# plot settings
output_format = 'script'  # 'script' or 'plot'

# model settings
s0, e0, i0, r0 = .999, .000, .001, .000
beta = .10
gamma = .025
sigma = .10
xi = .01
v = .0017
mu = .0124 / 360
nu = .0071 / 360
T = 1000

# bokeh plotting
t = np.arange(T + 1)


def iterate(beta, gamma, sigma, xi, v, mu, nu):
    global s, e, i, r
    n = s + e + i + r
    ds = n * mu + xi * r - s * (nu + beta * i)
    de = s * beta * i - e * (nu + sigma)
    di = e * sigma - i * (nu + gamma + v)
    dr = i * gamma - r * (nu + xi)
    s = max(s + ds, 0)
    e = max(e + de, 0)
    i = max(i + di, 0)
    r = max(r + dr, 0)
    return s, e, i, r


def simulate(initial_states, parameters):
    global s, e, i, r
    s, e, i, r = initial_states
    beta, gamma, sigma, xi, v, mu, nu = parameters
    state0 = (s0, e0, i0, r0)
    return zip(*([state0] + [iterate(beta, gamma, sigma, xi, v, mu, nu) for _ in t[1:]]))

s, e, i, r = simulate([s0, e0, i0, r0], [beta, gamma, sigma, xi, v, mu, nu])

output_notebook(hide_banner=True)
source = ColumnDataSource(data=dict(t=t, s=s, e=e, i=i, r=r))

p = figure(
    x_range=Range1d(0, T + 20, bounds=(0, T)),
    plot_width=400, plot_height=300, tools=''
)
handlers = p.varea_stack(
    ['s', 'e', 'i', 'r'], x='t', source=source,
    color=['#97D52D', '#FFB108', '#850606', '#69951d'])
handler_s, handler_e, handler_i, handler_r = handlers

p.y_range.start = 0
p.yaxis.formatter = NumeralTickFormatter(format='0%')
p.yaxis.major_label_orientation = 'vertical'
p.outline_line_alpha = 0
p.xgrid.grid_line_color = None
p.ygrid.grid_line_color = None
p.xaxis.major_label_text_font = 'Palatino'
p.yaxis.major_label_text_font = 'Palatino'
p.xaxis.major_label_text_font_size = '10pt'
p.yaxis.major_label_text_font_size = '10pt'
p.xaxis.axis_label_text_font = 'Palatino'
p.yaxis.axis_label_text_font = 'Palatino'
p.xaxis.axis_label_text_font_size = '12pt'
p.yaxis.axis_label_text_font_size = '12pt'
p.xaxis.axis_label = 'Time (days)'
p.yaxis.axis_label = 'Population'
p.title.text = 'Deadly SEIRS Model'
p.title.text_font = 'Palatino'
p.title.text_font_size = '15pt'
p.title.align = 'center'

legend = Legend(items=[
    ('Susceptible', [handler_s]),
    ('Exposed', [handler_e]),
    ('Infectious', [handler_i]),
    ('Recovered', [handler_r])
], location='center', orientation='horizontal')
legend.label_text_font = 'Palatino'
legend.label_text_font_size = '10pt'
legend.background_fill_alpha = 0
legend.border_line_alpha = 0
p.add_layout(legend, 'below')

e0_slider = Slider(start=0, end=.1, value=e0, step=.001,
                   title='Initial Exposed', format='0.0%')
i0_slider = Slider(start=0, end=.1, value=i0, step=.001,
                   title='Initial Infectious', format='0.0%')
r0_slider = Slider(start=0, end=.1, value=r0, step=.001,
                   title='Initial Recovered', format='0.0%')
beta_slider = Slider(start=0, end=1, value=beta, step=.01,
                     title='Infectious Rate', format='0%')
gamma_slider = Slider(start=0, end=.1, value=gamma, step=.001, format='0.0%',
                      title='Recovery Rate (1 / Days of Treatement)')
sigma_slider = Slider(start=0, end=1, value=sigma, step=.01, format='0%',
                      title='Incubation Rate (1 / Days before Symptoms)')
xi_slider = Slider(start=0, end=1, value=xi, step=.01, format='0%',
                   title='Resusceptibility Rate (1 / Days of Immunity)')
v_slider = Slider(start=0, end=.005, value=v, step=.0001,
                   title='Disease Fatality Rate', format='0.00%')
mu_slider = Slider(start=0, end=.2, value=mu * 360, step=.0001,
                   title='Annual Birth Rate', format='0.00%')
nu_slider = Slider(start=0, end=.2, value=nu * 360, step=.0001,
                   title='Annual Death Rate', format='0.00%')

func_js = '''
    const data = source.data;
    var t = data['t'];
    var s = data['s']; var _s = s0;
    var e = data['e']; var _e = e0;
    var i = data['i']; var _i = i0;
    var r = data['r']; var _r = r0;
    
    var ds, de, di, dr;
    for (var j = 0; j < t.length; j++) {
        s[j] = Math.max(_s, 0);
        e[j] = Math.max(_e, 0);
        i[j] = Math.max(_i, 0);
        r[j] = Math.max(_r, 0);

        nj = _s + _e + _i + _r;
        ds = nj * mu + xi * _r - _s * (nu + beta * _i);
        de = _s * beta * _i - _e * (nu + sigma);
        di = _e * sigma - _i * (nu + gamma + v);
        dr = _i * gamma - _r * (nu + xi);

        _s += ds;
        _e += de;
        _i += di;
        _r += dr;
    }
    source.change.emit();
    p.reset.emit();
'''

callback = CustomJS(
    args=dict(
        p=p, source=source,
        e0_slider=e0_slider, i0_slider=i0_slider, r0_slider=r0_slider,
        beta_slider=beta_slider, gamma_slider=gamma_slider, sigma_slider=sigma_slider,
        xi_slider=xi_slider, v_slider=v_slider, mu_slider=mu_slider, nu_slider=nu_slider
    ),
    code='''
    var e0 = e0_slider.value;
    var i0 = i0_slider.value;
    var r0 = r0_slider.value;
    var s0 = 1 - e0 - i0 - r0;
    var beta = beta_slider.value;
    var gamma = gamma_slider.value;
    var sigma = sigma_slider.value;
    var xi = xi_slider.value;
    var v = v_slider.value;
    var mu = mu_slider.value / 360;
    var nu = nu_slider.value / 360;
    ''' + func_js
)

e0_slider.js_on_change('value', callback)
i0_slider.js_on_change('value', callback)
r0_slider.js_on_change('value', callback)
beta_slider.js_on_change('value', callback)
gamma_slider.js_on_change('value', callback)
sigma_slider.js_on_change('value', callback)
xi_slider.js_on_change('value', callback)
v_slider.js_on_change('value', callback)
mu_slider.js_on_change('value', callback)
nu_slider.js_on_change('value', callback)


p.js_on_event(DoubleTap, CustomJS(
    args=dict(
        p=p, source=source,
        s0=s0, e0=e0, i0=i0, r0=r0,
        beta=beta, gamma=gamma, sigma=sigma, xi=xi, v=v, mu=mu, nu=nu,
        e0_slider=e0_slider, i0_slider=i0_slider, r0_slider=r0_slider,
        beta_slider=beta_slider, gamma_slider=gamma_slider, sigma_slider=sigma_slider,
        xi_slider=xi_slider, v_slider=v_slider, mu_slider=mu_slider, nu_slider=nu_slider
    ),
    code='''
    e0_slider.value = e0;
    i0_slider.value = i0;
    r0_slider.value = r0;
    beta_slider.value = beta;
    gamma_slider.value = gamma;
    sigma_slider.value = sigma;
    xi_slider.value = xi;
    v_slider.value = v;
    mu_slider.value = mu * 360;
    nu_slider.value = nu * 360;
    ''' + func_js
))

c = column(
    i0_slider, e0_slider, r0_slider,
    beta_slider, gamma_slider, sigma_slider, xi_slider, v_slider, mu_slider, nu_slider
)

layout = gridplot(
    [[p], [c]],
    sizing_mode='scale_width', toolbar_location=None, toolbar_options=dict(logo=None)
)


# output in the desired format
if output_format == 'plot':
    show(layout)
elif output_format == 'script':
    script, div = components(layout)
    script = '<script src="https://cdnjs.cloudflare.com/ajax/libs/bokeh/1.3.4/bokeh.min.js"></script>' + \
             '<script src="https://cdnjs.cloudflare.com/ajax/libs/bokeh/1.3.4/bokeh-widgets.min.js"></script>' + \
             '<script src="https://cdnjs.cloudflare.com/ajax/libs/bokeh/1.3.4/bokeh-tables.min.js"></script>' + \
             script + '\n' + div
    pyperclip.copy(script)
else:
    raise ValueError('unknown output format')