Fast prototyping of epidemiological models based on reaction equations. Analyze the ODEs analytically or numerically, or run/animate stochastic simulations on networks/well-mixed systems.
- repository: https://github.com/benmaier/epipack/
- documentation: http://epipack.benmaier.org/
import epipack as epk
from epipack.vis import visualize
import netwulf as nw
network, _, __ = nw.load('cookbook/readme_vis/MHRN.json')
N = len(network['nodes'])
links = [ (l['source'], l['target'], 1.0) for l in network['links'] ]
S, I, R = list("SIR")
model = epk.StochasticEpiModel([S,I,R],N,links)\
.set_link_transmission_processes([ (I, S, 1.0, I, I) ])\
.set_node_transition_processes([ (I, 1.0, R) ])\
.set_random_initial_conditions({ S: N-5, I: 5 })
visualize(model, network, sampling_dt=0.1)
Simple compartmental models of infectious diseases are useful
to investigate effects of certain processes on disease dissemination.
Using pen and paper, quickly adding/removing compartments and transition processes
is easy, yet the analytical and numerical analysis or stochastic simulations
can be tedious to set up and debug—especially when the model changes (even slightly).
epipack aims at streamlining this process
such that all the analysis steps can be performed in an efficient manner,
simply by defining processes based on reaction equations. epipack provides
three main base classes to accomodate different problems.
EpiModel: Define a model based on transition, birth, death, fission, fusion, or transmission reactions, integrate the ordinary differential equations (ODEs) of the corresponding well-mixed system numerically or simulate it using Gillespie's algorithm. Process rates can be numerical functions of time and the system state.SymbolicEpiModel: Define a model based on transition, birth, death, fission, fusion, or transmission reactions. Obtain the ODEs, fixed points, Jacobian, and the Jacobian's eigenvalues at fixed points as symbolic expressions. Process rates can be symbolic expressions of time and the system state. Set numerical parameter values and integrate the ODEs numerically or simulate the stochastic systems using Gillespie's algorithm.StochasticEpiModel: Define a model based on node transition and link transmission reactions. Add conditional link transmission reactions. Simulate your model on any (un-/)directed, (un-/)weighted static/temporal network, or in a well-mixed system.
Additionally, epipack provides a visualization framework to animate
stochastic simulations on networks, lattices, well-mixed systems,
or reaction-diffusion systems based on MatrixEpiModel.
Check out the Example section for some demos.
Note that the internal simulation algorithm for network simulations is based on the following paper:
"Efficient sampling of spreading processes on complex networks using a composition and rejection algorithm", G.St-Onge, J.-G. Young, L. Hébert-Dufresne, and L. J. Dubé, Comput. Phys. Commun. 240, 30-37 (2019), http://arxiv.org/abs/1808.05859.
pip install epipack
epipack was developed and tested for
- Python 3.6
- Python 3.7
- Python 3.8
So far, the package's functionality was tested on Mac OS X and CentOS only.
epipack directly depends on the following packages which will be installed by pip during the installation process
numpy>=1.17scipy>=1.3sympy==1.6pyglet<1.6matplotlib>=3.0.0ipython>=7.14.0ipywidgets>=7.5.1
Please note that fast network simulations are only available if you install
SamplableSet==2.0(SamplableSet)
manually (pip won't do it for you).
The full documentation is available at epipack.benmaier.org.
Changes are logged in a separate file.
This project is licensed under the MIT License. Note that this excludes any images/pictures/figures shown here or in the documentation.
If you want to contribute to this project, please make sure to read the code of conduct and the contributing guidelines. In case you're wondering about what to contribute, we're always collecting ideas of what we want to implement next in the outlook notes.
Let's define an SIRS model with infection rate eta, recovery rate rho, and waning immunity rate omega and analyze the system
Define a pure numeric model with EpiModel.
Integrate the ODEs or simulate the system stochastically.
from epipack import EpiModel
import matplotlib.pyplot as plt
import numpy as np
S, I, R = list("SIR")
N = 1000
SIRS = EpiModel([S,I,R],N)\
.set_processes([
#### transmission process ####
# S + I (eta=2.5/d)-> I + I
(S, I, 2.5, I, I),
#### transition processes ####
# I (rho=1/d)-> R
# R (omega=1/14d)-> S
(I, 1, R),
(R, 1/14, S),
])\
.set_initial_conditions({S:N-10, I:10})
t = np.linspace(0,40,1000)
result_int = SIRS.integrate(t)
t_sim, result_sim = SIRS.simulate(t[-1])
for C in SIRS.compartments:
plt.plot(t, result_int[C])
plt.plot(t_sim, result_sim[C])
It's also straight-forward to define temporally varying (functional) rates.
import numpy as np
from epipack import SISModel
N = 100
recovery_rate = 1.0
def infection_rate(t, y, *args, **kwargs):
return 3 + np.sin(2*np.pi*t/100)
SIS = SISModel(
infection_rate=infection_rate,
recovery_rate=recovery_rate,
initial_population_size=N
)\
.set_initial_conditions({
'S': 90,
'I': 10,
})
t = np.arange(200)
result_int = SIS.integrate(t)
t_sim, result_sim = SIS.simulate(199)
for C in SIS.compartments:
plt.plot(t_sim, result_sim[C])
plt.plot(t, result_int[C])
Symbolic models are more powerful because they can do the same as the pure numeric models while also offering the possibility to do analytical evaluations
from epipack import SymbolicEpiModel
import sympy as sy
S, I, R, eta, rho, omega = sy.symbols("S I R eta rho omega")
SIRS = SymbolicEpiModel([S,I,R])\
.set_processes([
(S, I, eta, I, I),
(I, rho, R),
(R, omega, S),
]) Print the ODE system in a Jupyter notebook
>>> SIRS.ODEs_jupyter()
Get the Jacobian
>>> SIRS.jacobian()
Find the fixed points
>>> SIRS.find_fixed_points()
Get the eigenvalues at the disease-free state in order to find the epidemic threshold
>>> SIRS.get_eigenvalues_at_disease_free_state()
{-omega: 1, eta - rho: 1, 0: 1}Set numerical parameter values and integrate the ODEs numerically
>>> SIRS.set_parameter_values({eta: 2.5, rho: 1.0, omega:1/14})
>>> t = np.linspace(0,40,1000)
>>> result = SIRS.integrate(t)If set up as
>>> N = 10000
>>> SIRS = SymbolicEpiModel([S,I,R],N)the system can simulated directly.
>>> t_sim, result_sim = SIRS.simulate(40)Let's set up some temporally varying rates
from epipack import SymbolicEpiModel
import sympy as sy
S, I, R, eta, rho, omega, t, T = \
sy.symbols("S I R eta rho omega t T")
N = 1000
SIRS = SymbolicEpiModel([S,I,R],N)\
.set_processes([
(S, I, 2+sy.cos(2*sy.pi*t/T), I, I),
(I, rho, R),
(R, omega, S),
])
SIRS.ODEs_jupyter()
Now we can integrate the ODEs or simulate the system using Gillespie's SSA for inhomogeneous Poisson processes.
import numpy as np
SIRS.set_parameter_values({
rho : 1,
omega : 1/14,
T : 100,
})
SIRS.set_initial_conditions({S:N-100, I:100})
_t = np.linspace(0,200,1000)
result = SIRS.integrate(_t)
t_sim, result_sim = SIRS.simulate(max(_t))
epipack offers a classs called InteractiveIntegrator
that allows an interactive exploration of a system
in a Jupyter notebook.
Make sure to first run
%matplotlib widget
in a cell.
from epipack import SymbolicEpiModel
from epipack.interactive import InteractiveIntegrator, Range, LogRange
import sympy
S, I, R, R0, tau, omega = sympy.symbols("S I R R_0 tau omega")
I0 = 0.01
model = SymbolicEpiModel([S,I,R])\
.set_processes([
(S, I, R0/tau, I, I),
(I, 1/tau, R),
(R, omega, S),
])\
.set_initial_conditions({S:1-I0, I:I0})
# define a log slider, a linear slider and a constant value
parameters = {
R0: LogRange(min=0.1,max=10,step_count=1000),
tau: Range(min=0.1,max=10,value=8.0),
omega: 1/14
}
t = np.logspace(-3,2,1000)
InteractiveIntegrator(model, parameters, t, figsize=(4,4))
Let's simulate an SIRS system on a random graph (using the parameter definitions above).
from epipack import StochasticEpiModel
import networkx as nx
k0 = 50
R0 = 2.5
rho = 1
eta = R0 * rho / k0
omega = 1/14
N = int(1e4)
edges = [ (e[0], e[1], 1.0) for e in \
nx.fast_gnp_random_graph(N,k0/(N-1)).edges() ]
SIRS = StochasticEpiModel(
compartments=list('SIR'),
N=N,
edge_weight_tuples=edges
)\
.set_link_transmission_processes([
('I', 'S', eta, 'I', 'I'),
])\
.set_node_transition_processes([
('I', rho, 'R'),
('R', omega, 'S'),
])\
.set_random_initial_conditions({
'S': N-100,
'I': 100
})
t_s, result_s = SIRS.simulate(40)
Likewise, it's straight-forward to visualize this system
>>> from epipack.vis import visualize
>>> from epipack.networks import get_random_layout
>>> layouted_network = get_random_layout(N, edges)
>>> visualize(SIRS, layouted_network, sampling_dt=0.1, config={'draw_links': False})
A lattice is nothing but a network, we can use
get_grid_layout and get_2D_lattice_links
to set up a visualization.
from epipack.vis import visualize
from epipack import (
StochasticSIRModel,
get_2D_lattice_links,
get_grid_layout
)
# define links and network layout
N_side = 100
N = N_side**2
links = get_2D_lattice_links(N_side, periodic=True, diagonal_links=True)
lattice = get_grid_layout(N)
# define model
R0 = 3; recovery_rate = 1/8
model = StochasticSIRModel(N,R0,recovery_rate,
edge_weight_tuples=links)
model.set_random_initial_conditions({'I':20,'S':N-20})
sampling_dt = 1
visualize(model,lattice,sampling_dt,
config={
'draw_nodes_as_rectangles':True,
'draw_links':False,
}
)
Since reaction-diffusion systems in discrete space
can be interpreted as being based on reaction
equations, we can set those up using epipack's
framework.
Checkout the docs on Reaction-Diffusion Systems.
Every node in a network is associated
with a compartment and we're using MatrixEpiModel
because it's faster than EpiModel.
from epipack import MatrixEpiModel
N = 100
base_compartments = list("SIR")
compartments = [ (node, C) for node in range(N) for C in base_compartments ]
model = MatrixEpiModel(compartments)Now, we define both epidemiological
and movement processes on a hypothetical
list links.
infection_rate = 2
recovery_rate = 1
mobility_rate = 0.1
quadratic_processes = []
linear_processes = []
for node in range(N):
quadratic_processes.append(
( (node, "S"), (node, "I"), infection_rate, (node, "I"), (node, "I") ),
)
linear_processes.append(
( (node, "I"), recovery_rate, (node, "R") )
)
for u, v, w in links:
for C in base_compartments:
linear_processes.extend([
( (u, C), w*mobility_rate, (v, C) ),
( (v, C), w*mobility_rate, (u, C) ),
])
Fork this repository, clone it, and install it in dev mode.
git clone git@github.com:YOURUSERNAME/epipack.git
makeIf you want to upload to PyPI, first convert the new README.md to README.rst
make readmeIt will give you warnings about bad .rst-syntax. Fix those errors in README.rst. Then wrap the whole thing
make pypiIt will probably give you more warnings about .rst-syntax. Fix those until the warnings disappear. Then do
make upload