SEIR-Type models¶

Introduction¶

In epidemiology the most widespread type of simple model is the Lotka-Volterra-like set of coupled ordinary differential equations of Kermack and McKendrick (1927), and its variants (SIR, SEIR, SECIR, SEIRD etc etc). The letters SEIR stand for the “Susceptible”, “Exposed”, “Infected” and “Recovered” portions of the population. These models are called “compartimental models” because the S, E, I and R portions of the population are called compartiments, and individuals move from the S \(\rightarrow\) E \(\rightarrow\) I \(\rightarrow\) R compartiments. We do not include birth and natural death in the model.

Additional compartiments such as “Deceased”, “Carrier” or a variety of other compartiments can be added to these models, giving rise to a wide variety of compartimental model types.

The EpiDemo package includes integrators for the SIR and SEIR models, but only for the purpose of comparison to the preferred model: the lookback integrator. Because we implement these models only for comparison purpose, we do not provide the corresponding GUIs for them. But you can look how the GUIs for the lookback integrator were built, and build your own GUI for the SIR or SEIR models, if you want.

The SIR model¶

The simplest of the compartimental models is the SIR model with the “Susceptible”, “Infected” and “Recovered” compartiments. Susceptible means that an individual can be infected (is not immune). Infected means, an individual is infectuous. Recovered means the individual is no longer infectuous. We define \(S(t)\), \(I(t)\) and \(R(t)\) as normalized functions:

\[\begin{split}\begin{split} S(t) &= N_S(t)/N_{\mathrm{pop}} \\ I(t) &= N_I(t)/N_{\mathrm{pop}} \\ R(t) &= N_R(t)/N_{\mathrm{pop}} \end{split}\end{split}\]

The \(N_S\), \(N_I\) and , \(N_R\) are the number of people in each of the compartiments. We have \(N_S+N_I+N_R=N_{\mathrm{pop}}\), and likewise \(S+I+R=1\).

The equations are:

\[\begin{split}\begin{split} \frac{dS}{dt} &= - \beta\, S\, I\\ \frac{dI}{dt} &= \beta\, S\, I -\gamma\, I\\ \frac{dR}{dt} &= \gamma\, I \end{split}\end{split}\]

where (at least in EpiDemo) we define time \(t\) in units of days. The coefficient \(\beta\) is the rate by which susceptible individuals are being infected. It represents the transition from the S compartiment to the I compartiment. The coefficient \(\gamma\) is the rate by which an infectuous individual recovers. It is associated by the mean duration of illness (or better: of the infectuousness) \(t_{\mathrm{ifc}}\) through

\[\gamma = \frac{1}{t_{\mathrm{ifc}}}\]

The basic reproduction number \(R_0\) is defined as the mean number of new individuals a single infectuous individual can infect in a fully susceptible population. On average the duration of infectuousness is \(t_{\mathrm{ifc}}\) and the rate is \(\beta\), so by definition:

\[R_0 = \frac{\beta}{\gamma}\]

From epidemiological data you typically do not know the \(t_{\mathrm{ifc}}\). Instead, during the initial phases of an epidemic, you measure the doubling time scale \(t_2\) or the ten-folding time scale \(t_{10}\) or the e-folding time scale \(t_e\):

\[t_e = t_2/\ln(2) = t_{10}/\ln(10) = \frac{t_1-t_0}{\ln(I(t_1))-\ln(I(t_0))}\]

From \(R_0\) and \(t_e\) you can compute the value of \(\gamma\):

\[\gamma = \frac{1}{(R_0-1)t_e}\]

and then you can compute \(\beta\) from:

\[\beta = \gamma R_0\]

The fact that \(\gamma = 1/(R_0-1)t_e\) can be derived by finding an exponential solution for the case \(S=1\), i.e. the very beginning of the epidemic.

The \(R_0\) value is usually computed from the rise of case numbers by defining the generation time \(t_g\). For instance, the Robert Koch Institute (RKI) computes \(R_0\) under the assumption that \(t_g=4\) days. The \(R_0\) is then defined as:

\[R_0 = \frac{I(t+t_g)}{I(t)} = \frac{\delta N_{\mathrm{inf}}(t+t_e)}{\delta N_{\mathrm{inf}}(t)}\]

during the exponetial rise stage at the beginning, when \(S=1\). In fact, given that \(R_0\) is already well defined, this equation is actually the definition of the generation time \(t_g\). This can be expressed more explicitly:

\[t_g = t_e\,\ln R_0\]

So if you get, for instance from the reports of the RKI the value for \(R_0\), and knowing that the RKI assumes \(t_g=4\) days, you can compute \(t_e\), then \(\gamma\), and finally \(\beta\). The duration of the illness and infectuousness, \(t_{\mathrm{ifc}}=1/\gamma\) is then an outcome of this calculation. If this duration is not quite exact, it will not affect the exponential rise, but it might affect somewhat the expected herd immunity.

Conversely, if you are sure you know the value of \(t_{\mathrm{ifc}}\), and you measure \(t_e\) from the development of the case numbers, you can compute \(R_0\) by solving \(\gamma = 1/(R_0-1)t_e\), and then compute \(\beta\) from \(\beta=\gamma R_0\), and you can compute the generation time from \(t_g = t_e\,\ln R_0\).

The SIR model has two exponential stages: an exponential rise, followed by an exponential decay after herd immunity.

Note that there are several conceptual differences between the SIR model and the lookback model:

The \(N_S\), \(N_I\) and \(N_R\) are the number of individuals having the current status “susceptible”, “infectuous” or “recovered”. In the lookback model, on the other hand, \(\delta N_{\mathrm{inf}}\) is the number of new people being infected that day and \(N_{\mathrm{inf}}\) is the cumulative version of that. In the lookback model there are no compartiments by default: you have to compute these numbers aposteriori.
In the compartimental models the equations are coupled sets of ordinary differential equations. Numerically integrating these equations is most accurate for the smallest time steps, but good numerical integrators can get accurate enough results for time steps of 1 day. The lookback integrator is an autoregression equation. The equations are designed to be integrated by time steps of 1 day (in the formulation used in EpiDemo).
In the SIR model (and all compartimental models) the individuals within a compartiment have no memory of how long they have been in that compartiment. There is no well-defined time after becoming infected that the patient will recover. Instead, each patient (member of the I compartiment) has a fixed chance each day (\(\gamma\)) of recovering. It is like rolling a die each day. Or to the physicist it is like radioactive decay. Whether this is a realistic description of the real progression of a disease is a question I can not answer here. The lookback model can also treat illness in this way, by setting the weights of the autoregression accordingly. But is can also handle more deterministic progressions of a disease. In the simplest version of the lookback model the duration of infectuousness is exactly defined as a time window. As a result, during the post-herd-immunity phase, the simple autoregression model with a constant infection rate in a fixed time window decays faster than the SIR and SEIR models.

In EpiDemo you can run SIR models as follows:

import epidemo
Npop = 83000000
model = epidemo.SIRModel(Npop,1.3/4,1/4,'2020-08-16','2021-06-01')
model.reset_init_infected(300)
model.run()

import matplotlib.pyplot as plt
import matplotlib.dates as mdates
fig, axes = plt.subplots()
plt.semilogy(model.time,model.Npop*model.S,label='S')
plt.semilogy(model.time,model.Npop*model.I,label='I')
plt.semilogy(model.time,model.Npop*model.R,label='R')
axes.xaxis.set_major_locator(mdates.MonthLocator([1,3,5,7,9,11]))
axes.xaxis.set_major_formatter(mdates.DateFormatter('%d %b'))
plt.ylim(ymin=1e2,ymax=1e8)
plt.legend()
plt.show()

The SEIR model¶

In the SEIR model an extra compartiment E is added, to account for a incubation time between a person becoming infected (called E for “Exposed”), and becoming infectuous (called I for “Infectuous”). The transition from E to I is a rate \(\sigma\). The mean time in the E compartiment is then \(t_{\mathrm{xp}}=1/\sigma\) (where xp is used instead of exp, to avoid confusion with the exponential growth time scale \(t_e\)). The equations are then:

\[\begin{split}\begin{split} \frac{dS}{dt} &= - \beta\, S\, I\\ \frac{dE}{dt} &= \beta\, S\, I -\sigma\, E\\ \frac{dI}{dt} &= \sigma\, E -\gamma\, I\\ \frac{dR}{dt} &= \gamma I \end{split}\end{split}\]

Several things remain the same as in the SIR model:

\[\begin{split}\begin{split} R_0 &= \frac{\beta}{\gamma}\\ t_g &= t_e\,\ln R_0 \end{split}\end{split}\]

But some things chance with respect to the SIR model. For instance the exponential growth time scale \(t_e\) and the new parameter \(\sigma\) are now related via:

\[\left(\frac{1}{\sigma t_e}+1\right)\left(\frac{1}{t_e}+\gamma\right) = \beta\]

By setting \(\beta=\gamma R_0\), you can solve for \(\gamma\):

\[\gamma = \left[\frac{R_0}{\left(\frac{1}{\sigma t_e}+1\right)}-1\right]^{-1}\frac{1}{t_e}\]

and then compute \(\beta\). Conversely, if you know \(\beta\), \(\sigma\), and \(\gamma\), you can compute the exponential growth time: by solving the above quadratic equations, obtaining:

\[t_e = \frac{\sigma+\gamma}{2\sigma(\beta-\gamma)} \left(1+\sqrt{1+\frac{4\sigma(\beta-\gamma)}{(\sigma+\gamma)^2}}\right)\]

It shows that in the SEIR model the relation between the measured parameters and the coefficients of the equations can be a bit tricky. In practice one would simply tune the parameters to the data using, for instance, an automatic fitting procedure (such as Markov-Chain-Monte-Carlo).

In addition to the differences already discussed for the SIR model, the SEIR model has an additional difference with the lookback model:

The incubation time is, as with all compartiments, merely a mean. The transition from E into I is a rate, and while the mean incubation time \(t_{\mathrm{inc}}\) is a parameter of the model (\(\sigma=1/t_{\mathrm{inc}}\)), the most likely incubation time is still 0. In contast, in the lookback model there is typically a time period after being infected with zero chance of already being infectuous, although the “decay-like” transition from E to I can be mimicked in the lookback model by using appropriate weighting values.

In EpiDemo you can run a SEIR model as follows:

import epidemo
Npop = 83000000
model = epidemo.SEIRModel(Npop,1.3/3,1/1,1/3,'2020-08-16','2021-06-01')
model.reset_init_infected(300)
model.run()

import matplotlib.pyplot as plt
import matplotlib.dates as mdates
fig, axes = plt.subplots()
plt.semilogy(model.time,model.Npop*model.S,label='S')
plt.semilogy(model.time,model.Npop*model.E,label='E')
plt.semilogy(model.time,model.Npop*model.I,label='I')
plt.semilogy(model.time,model.Npop*model.R,label='R')
axes.xaxis.set_major_locator(mdates.MonthLocator([1,3,5,7,9,11]))
axes.xaxis.set_major_formatter(mdates.DateFormatter('%d %b'))
plt.ylim(ymin=1e2,ymax=1e8)
plt.legend()
plt.show()