Fast prototyping of epidemiological models based on reaction equations. Analyze the ODEs analytically or numerically, or run stochastic simulations on networks/well-mixed systems.

## Project description

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)

## Idea

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.

## Install

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.

## Dependencies

`epipack` directly depends on the following packages which will be
installed by `pip` during the installation process

`numpy>=1.17``scipy>=1.3``sympy==1.6``pyglet<1.6``matplotlib>=3.0.0``ipython>=7.14.0``ipywidgets>=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).

## Documentation

The full documentation is available at epipack.benmaier.org.

## Changelog

Changes are logged in a separate file.

## License

This project is licensed under the MIT License. Note that this excludes any images/pictures/figures shown here or in the documentation.

## Contributing

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.

## Examples

Let’s define an SIRS model with infection rate `eta`, recovery rate
`rho`, and waning immunity rate `omega` and analyze the system

### Pure Numeric Models

#### Basic Definition (EpiModel)

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])

#### Functional Rates

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

#### Basic Definition

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), ])

#### Analytical Evaluations

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}

#### Numerical Evaluations

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)

#### Temporally Varying Rates

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))

#### Interactive Analyses

`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))

### Pure Stochastic Models

#### On a Network

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)

#### Visualize

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})

#### On a Lattice

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, } )

### Reaction-Diffusion Models

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) ), ])

## Dev notes

Fork this repository, clone it, and install it in dev mode.

git clone git@github.com:YOURUSERNAME/epipack.git make

If you want to upload to PyPI, first convert the new `README.md` to
`README.rst`

make readme

It will give you warnings about bad `.rst`-syntax. Fix those errors in
`README.rst`. Then wrap the whole thing

make pypi

It will probably give you more warnings about `.rst`-syntax. Fix those
until the warnings disappear. Then do

make upload

## Project details

## Download files

Download the file for your platform. If you're not sure which to choose, learn more about installing packages.

Filename, size | File type | Python version | Upload date | Hashes |
---|---|---|---|---|

Filename, size epipack-0.1.5.tar.gz (93.4 kB) | File type Source | Python version None | Upload date | Hashes View |