In [1]:
from collections import Counter
from pprint import pprint
import pandas as pd
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import warnings
warnings.filterwarnings("ignore")
import epidemik
from epidemik import EpiModel
import watermark
%load_ext watermark
%matplotlib inline
We start by print out the versions of the libraries we're using for future reference
In [2]:
%watermark -n -v -m -g -iv
Python implementation: CPython Python version : 3.11.7 IPython version : 8.12.3 Compiler : Clang 14.0.6 OS : Darwin Release : 23.4.0 Machine : arm64 Processor : arm CPU cores : 16 Architecture: 64bit Git hash: 6f3f47aa98b5141fcddbadd07afeedaa91c4cfb7 pandas : 2.1.4 epidemik : 0.0.11 matplotlib: 3.8.0 watermark : 2.4.3 numpy : 1.26.4
Load default figure style and extract the color cycle for ease of reference
In [3]:
plt.style.use('./d4sci.mplstyle')
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
SIR Model¶
Our first step is to reimplement the SIR model from the previous post
In [4]:
beta = 0.2
mu = 0.1
SIR = EpiModel()
SIR.add_interaction('S', 'I', 'I', beta)
SIR.add_spontaneous('I', 'R', mu)
We can easily make sure that it looks correct
In [5]:
print(SIR)
Epidemic Model with 3 compartments and 2 transitions: S + I = I 0.200000 I -> R 0.100000 R0=2.00
And we integrate it using our standard parameter values
In [6]:
N = 100000
I0 = 10
SIR.integrate(365, S=N-I0, I=I0, R=0)
And a quick visualization
In [7]:
ax = SIR.plot(normed=True)
ax.legend(['Susceptible', 'Infectious', 'Recovered'])
Out[7]:
<matplotlib.legend.Legend at 0x148b6ff10>
This will be the basis that we will build up on.
Incubation Period - SEIR Model¶
One of the main factors that is missing from the simple SIR model, is that a Latent or incubation period. We can add this feature to our model by creating a new compartment E(xposed) that is Infected, but not yet infectious. The parameter $\epsilon$ determines the rate at which individuals develop the ability to infect others
In [8]:
epsilon = 0.4
SEIR = EpiModel()
SEIR.add_interaction('S', 'E', 'I', beta)
SEIR.add_spontaneous('E', 'I', epsilon)
SEIR.add_spontaneous('I', 'R', mu)
Now we have a third transition
In [9]:
print(SEIR)
Epidemic Model with 4 compartments and 3 transitions: S + I = E 0.200000 E -> I 0.400000 I -> R 0.100000 R0=2.00
Which has a direct impact on the time line of the system
In [10]:
SEIR.integrate(365, S=N-I0, E=0, I=I0, R=0)
The value of $R_0$ remains the same, but the epidemic takes place over a longer timeframe
In [11]:
ax = SEIR.plot(normed=True, color=[colors[0], colors[4], colors[1], colors[2]])
ax.legend(['Susceptible', 'Exposed', 'Infectious', 'Recovered'])
Out[11]:
<matplotlib.legend.Legend at 0x14feb6e90>
We can easily illustrate this by simply plotting the infected compartments of the SIR and SEIR models, side by side.
In [12]:
ax = (SEIR.I/N).plot(color=colors[1], linestyle='-', label='SEIR')
(SIR.I/N).plot(color=colors[1], linestyle='--', label='SIR', ax=ax)
ax.set_xlabel('Time')
ax.set_ylabel('Infectious')
ax.legend()
Out[12]:
<matplotlib.legend.Legend at 0x148b51e10>
We clearly see that the presence of an Esposed compartment causes a delay in the number of infectious individuals as well as a flattening of the curve
Temporary Immunity¶
There have been recent reports that immunity might not be permanent for those who recover from the infection. We can easily add this feature to our model by including an extra transition from the Recovered compartment back to Susceptible
In [13]:
rho = 0.3
SEIRS = EpiModel()
SEIRS.add_interaction('S', 'E', 'I', beta)
SEIRS.add_spontaneous('E', 'I', epsilon)
SEIRS.add_spontaneous('I', 'R', mu)
SEIRS.add_spontaneous('R', 'S', rho)
Naturally, the time scale of the system depend on the exact value of $\rho$ chosen, but the main features will remain the same. Our model is now:
In [14]:
print(SEIRS)
Epidemic Model with 4 compartments and 4 transitions: S + I = E 0.200000 E -> I 0.400000 I -> R 0.100000 R -> S 0.300000
Let's integrate the equations of the model and examine the result
In [15]:
SEIRS.integrate(365, S=N-I0, E=0, I=I0, R=0)
ax = SEIRS.plot(normed=True, color=[colors[0], colors[4], colors[1], colors[2]])
ax.legend(['Susceptible', 'Exposed', 'Infectious', 'Recovered'])
Out[15]:
<matplotlib.legend.Legend at 0x16a2ed290>
In [16]:
rho = 0.03
SEIRS2 = EpiModel()
SEIRS2.add_interaction('S', 'E', 'I', beta)
SEIRS2.add_spontaneous('E', 'I', epsilon)
SEIRS2.add_spontaneous('I', 'R', mu)
SEIRS2.add_spontaneous('R', 'S', rho)
In [17]:
SEIRS2.integrate(365, S=N-I0, E=0, I=I0, R=0)
In [18]:
ax = SEIRS2.plot(normed=True, color=[colors[0], colors[4], colors[1], colors[2]])
ax.legend(['Susceptible', 'Exposed', 'Infectious', 'Recovered'])
Out[18]:
<matplotlib.legend.Legend at 0x17764e350>
In [19]:
ax = (SEIRS2.E/N).plot(color=colors[4])
(SEIRS2.I/N).plot(color=colors[1], ax=ax)
ax.legend(['Exposed', 'Infectious'])
ax.set_xlabel('Time')
ax.set_ylabel('Population')
Out[19]:
Text(0, 0.5, 'Population')
Asymptomatic cases¶
In [20]:
rbeta = 0.75
pa = 0.4
R0 = 2.0
beta = R0*mu/(pa*rbeta+(1-pa))
SEIIR = EpiModel()
SEIIR.add_interaction('S', 'E', 'Ia', rbeta*beta)
SEIIR.add_interaction('S', 'E', 'Is', beta)
SEIIR.add_spontaneous('E', 'Ia', epsilon*pa)
SEIIR.add_spontaneous('E', 'Is', epsilon*(1-pa))
SEIIR.add_spontaneous('Ia', 'R', mu)
SEIIR.add_spontaneous('Is', 'R', mu)
In [21]:
print(SEIIR)
Epidemic Model with 5 compartments and 6 transitions: S + Ia = E 0.166667 S + Is = E 0.222222 E -> Ia 0.160000 E -> Is 0.240000 Ia -> R 0.100000 Is -> R 0.100000 R0=2.00
In [22]:
SEIIR.integrate(365, S=N-I0, Ia=0, Is=I0, E=0, R=0)
In [23]:
fig, ax = plt.subplots(1)
(SEIIR.S/N).plot(color=colors[0], linestyle='-', label='Susceptible', ax=ax)
(SEIIR.E/N).plot(color=colors[4], linestyle='-', label='Exposed', ax=ax)
(SEIIR.Ia/N).plot(color=colors[1], linestyle=':', label='Asymptomatic', ax=ax)
(SEIIR.Is/N).plot(color=colors[1], linestyle='--', label='Symptomatic', ax=ax)
(SEIIR.R/N).plot(color=colors[2], linestyle='-', label='Recovered', ax=ax)
ax.legend()
ax.set_xlabel('Time')
ax.set_ylabel('Population')
Out[23]:
Text(0, 0.5, 'Population')
In [24]:
fig, ax = plt.subplots(1)
(SEIIR.E/N).plot(color=colors[4], linestyle='-', label='Exposed', ax=ax)
(SEIIR.Ia/N).plot(color=colors[1], linestyle=':', label='Asymptomatic', ax=ax)
(SEIIR.Is/N).plot(color=colors[1], linestyle='--', label='Symptomatic', ax=ax)
Emax = np.argmax(SEIIR.E)
Imax = np.argmax(SEIIR.Ia)
ax.axvspan(xmax=Imax, xmin=Emax, alpha=0.3, color=colors[3], zorder=-1)
ax.legend()
ax.set_xlabel('Time')
ax.set_ylabel('Population')
Out[24]:
Text(0, 0.5, 'Population')
Mortality¶
In [25]:
pd = 0.10
SEIIRD = EpiModel()
SEIIRD.add_interaction('S', 'E', 'Ia', rbeta*beta)
SEIIRD.add_interaction('S', 'E', 'Is', beta)
SEIIRD.add_spontaneous('E', 'Ia', epsilon*pa)
SEIIRD.add_spontaneous('E', 'Is', epsilon*(1-pa))
SEIIRD.add_spontaneous('Ia', 'R', mu)
SEIIRD.add_spontaneous('Is', 'R', (1-pd)*mu)
SEIIRD.add_spontaneous('Is', 'D', pd*mu)
In [26]:
print(SEIIRD)
Epidemic Model with 6 compartments and 7 transitions: S + Ia = E 0.166667 S + Is = E 0.222222 E -> Ia 0.160000 E -> Is 0.240000 Ia -> R 0.100000 Is -> R 0.090000 Is -> D 0.010000 R0=2.00
In [27]:
SEIIRD.integrate(365, S=N-I0, Ia=0, Is=I0, E=0, R=0, D=0)
In [28]:
fig, ax = plt.subplots(1)
(SEIIRD.S/N).plot(color=colors[0], linestyle='-', label='Susceptible', ax=ax)
(SEIIRD.E/N).plot(color=colors[4], linestyle='-', label='Exposed', ax=ax)
(SEIIRD.Ia/N).plot(color=colors[1], linestyle=':', label='Asymptomatic', ax=ax)
(SEIIRD.Is/N).plot(color=colors[1], linestyle='--', label='Symptomatic', ax=ax)
(SEIIRD.R/N).plot(color=colors[2], linestyle='-', label='Recovered', ax=ax)
(SEIIRD.D/N).plot(color=colors[7], linestyle='-', label='Dead', ax=ax)
ax.legend()
ax.set_xlabel('Time')
ax.set_ylabel('Population')
Out[28]:
Text(0, 0.5, 'Population')
In [29]:
fig, ax = plt.subplots(1)
(SEIIRD.E/N).plot(color=colors[4], linestyle='-', label='Exposed', ax=ax)
(SEIIRD.Ia/N).plot(color=colors[1], linestyle=':', label='Asymptomatic', ax=ax)
(SEIIRD.Is/N).plot(color=colors[1], linestyle='--', label='Symptomatic', ax=ax)
(SEIIRD.D/N).plot(color=colors[7], linestyle='-', label='Dead', ax=ax)
ax.legend()
ax.set_xlabel('Time')
Out[29]:
Text(0.5, 0, 'Time')
