EoN (Epidemics on Networks): a fast, flexible Python package for simulation, analytic approximation, and analysis of epidemics on networks

We provide a description of the Epidemics on Networks (EoN) python package designed for studying disease spread in static networks. The package consists of over $100$ methods available for users to perform stochastic simulation of a range of different processes including SIS and SIR disease, and generic simple or comlex contagions.


Introduction
EoN (EpidemicsOnNetworks) is a pure-python package designed to assist studies of infectious processes spreading through networks. It originally rose out of the book Mathematics of Epidemics on Networks [10], and now consists of over 100 user-functions.
EoN provides a set of tools for In this paper we provide brief descriptions with examples of a few of EoN's tools. The examples shown are intended to demonstrate the ability of the tools. The online documentation gives more detail about how to use them.
We model spreading processes on a contact network. In this context, many mathematicians and physicists are accustomed to thinking of individuals as nodes with their potentially infectious partnerships as edges. However, for those who come from other backgrounds this abstraction may be less familiar.
Therefore, we will describe a contact network along which an infections process spreads as consisting of "individuals" and "partnerships" rather than "nodes" and "edges". This has an additional benefit because in the simple contagion algorithm (described later), we need to define some other networks whose nodes represent possible statuses and whose edges represent transitions that can occur. Referring to "individuals" and "partnerships" when discussing the process spreading on the contact network makes it easier to avoid confusion between the different networks.
We start this paper by describing the tools for studying SIS and SIR disease through stochastic simulation and differential equations models. Then we discuss the simple and complex contagions, including examples showing how the simple contagion can be used to capture a range of standard disease models. Finally we demonstrate the visualization tools.

SIR and SIS disease Stochastic simulation
The stochastic SIR and SIS simulation tools allow the user to investigate standard SIS and SIR dynamics (SEIR/SIRS and other processes are addressed within the simple contagion model): • Markovian SIS and SIR simulations (fast SIS, Gillespie SIS, fast SIR, and Gillespie SIR).
• discrete time SIS and SIR simulations where infections last a single time step (basic discrete SIS, basic discrete SIR, and discrete SIR).
For both Markovian and non-Markovian methods it is possible for the transition rates to depend on intrinsic properties of individuals and of partnerships.
The continuous-time stochastic simulations have two different implementations: a Gillespie implementation [5,4] and an Event-driven implementation. Both approaches are efficient. They have similar speed if the dynamics are Markovian (depending on the network and disease parameters either may be faster than the other), but the event-driven implementation can also handle non-Markovian dynamics. In earlier versions, the event-driven simulations were consistently faster than the Gillespie simulations, and thus they are named fast SIR and fast SIS. The Gillespie simulations have since reached comparable speed using ideas from [8] and [2].
The algorithms can typically handle an SIR epidemic spreading on hundreds of thousands of individuals in well under a minute on a laptop. The SIS versions are slower because the number of events that happen is often much larger in an SIS simulation.

Examples
To demonstrate these, we begin with SIR simulations on an Erdős-Rényi network having a million individuals (in an Erdős-Rényi network each individual has identical probability of independently partnering with any other individual in the population).
import networkx as nx import EoN import matplotlib.pyplot as plt N = 10**6 #number of individuals kave = 5 #expected number of partners print('generating graph G with {} nodes'.format(N)) G = nx.fast_gnp_random_graph(N, kave/(N-1)) #Erdo''s-Re'nyi graph rho = 0.005 #initial fraction infected tau = 0.3 #transmission rate gamma = 1.0 #recovery rate print('doing event-based simulation') t1, S1, I1, R1 = EoN.fast_SIR(G, tau, gamma, rho=rho) #instead of rho, we could specify a list of nodes as initial_infecteds, or #specify neither and a single random node would be chosen as the index case. The run-times of fast SIR and Gillespie SIR are both comparable to the time taken to generate the million individual network G. The epidemics affect around 28 percent of the population. The differences between the simulations are entirely due to stochasticity.
We can perform similar simulations with an SIS epidemic. Because SIS epidemics take longer to simulate, we use a smaller network and specify the optional tmax argument defining the maximum stop time (by default tmax=100). We now consider an SIR disease spreading with non-Markovian dynamics. We assume that the infection duration is gamma distributed, but the transmission rate is constant (yielding an exponential distribution of time to transmission).
import networkx as nx import EoN import matplotlib.pyplot as plt import numpy as np #To reduce file size and make plotting faster, we'll just plot 1000 #data points. It's not really needed here, but this demonstrates #one of the available tools in EoN.

Differential Equations Models
EoN also provides a set of tools to numerically solve approximately 20 differential equations models for SIS or SIR disease spread in networks. The various models use different information about the network to make deterministic predictions about the number infected at different times. These use the Scipy integration tools. The derivations of the models and explanations of their simplifying assumptions are described in [10]. Depending on the model, we need different information about the network structure. The algorithms allow us to provide the information as inputs. However, there is also a version of each model which takes a network as an input instead and then measures the network properties.

Examples
We demonstrate an SIS pairwise model and an SIR edge-based compartmental model.
Our first example uses an SIS homogeneous pairwise model (section 4.3.3 of [10]). This model uses the average degree of the population and then attempts to track the number of [SI] and [SS] pairs. We assume a network with an average degree of 20. The initial condition is that a fraction ρ (rho) of the population is infected at random. import networkx as nx import EoN import matplotlib.pyplot as plt N=10000 gamma = 1 rho = 0.05 kave = 20

This produces
For ease of comparison with simulation and consistency with existing literature, the output of the model should be interpreted in terms of an expected number of individuals in each status, which requires that our values scale with N. So all of the initial conditions have a factor of N. If we are interested in proportion, we could arbitrarily set N=1, and then our solutions would give us the proportion of the population in each status.
Our second example uses an Edge-based compartmental model for an SIR disease (section 6.5.2 of [10] and also [17,14]). This model incorporates information about the degree distribution (i.e., how the number of partners is distributed), but assumes that the partnerships are selected as randomly as possible given this distribution. The model requires we define a "generating function" ψ(x) which is equal to the sum ∞ k=0 S k (0)x k where S k (0) is the proportion of all individuals in the population who both have k partners and are susceptible at t = 0. It also requires the derivative ψ (x) as well as φ S (0), the probability an edge from a susceptible node connects to another susceptible node at time 0. By default, it assumes there are no recovered individuals at time 0.
If the population has a Poisson degree distribution with mean kave and the infection is introduced by randomly infecting a proportion ρ of the population at time 0, then ψ( To be consistent with the other differential equations models, this EBCM implementation returns the expected number in each status, rather than the expected proportion. Most of the literature on the EBCM approach [17] focuses on expected proportion. By setting N=1, we have found the proporitons of the population in each status.

Simple and Complex Contagions
There are other contagious processes in networks which have received attention. Many of these fall into one of two types, "simple contagions" and "complex contagions".
In a "simple contagion" an individual u may be induced to change status by an interaction with its partner v. This status change occurs with the same rate regardless of the statuses of other partners of u (although the other partners may cause u to change to another status first). SIS and SIR diseases are special cases of simple contagions.
In a "complex contagion" however, we permit the rate at which u changes from one status to another to depend on the statuses of others in some more complicated way. Two infected individuals may cause a susceptible individual to become infected at some higher rate than would result from them acting independently. This is frequently thought to model social contagions where an individual may only believe something if multiple partners believe it [1].
The simple and complex contagions are currently implemented only in a Gillespie setting, and so they require Markovian assumptions. Although they are reasonably fast, it would typically be feasible to make a bespoke algorithm that runs significantly faster.

Simple contagions
EoN provides a function Gillespie simple contagion which allows a user to specify the rules governing an arbitrary simple contagion.
Examples are provided in the online documentation, including • SEIR disease (there is an exposed state before becoming infectious) • SIRS disease (recovered individuals eventually become susceptible again) • SIRV disease (individuals may get vaccinated) • Competing SIR diseases (there is cross immunity) • Cooperative SIR diseases (infection with one disease helps spread the other) The implementation requires the user to separate out two distinct ways that transitions occur: those that are intrinsic to an individual's current state and those that are induced by a partner. To help demonstrate, consider an "SEIR" epidemic, where individuals begin susceptible, but when they interact with infectious partners they may enter an exposed state. They remain in that exposed state for some period of time before transitioning into the infectious state independently of the status of any partner. They remain infectious and eventually transition into the recovered state, again independently of the status of any partner. Here the "E" to "I" and "I" to "R" transitions are intrinsic to the individual's state, while the "S" to "E" transition is induced by a partner.
To formalize this, we identify two broad types of transitions: • Spontaneous Transitions: Sometimes individuals change status without influence from any other individual. For example, an infected individual may recover, or an exposed individual may move into the infectious class. These transitions between statuses can be represented by a directed graph H where the nodes are not the original individuals of the contact network G, but rather the potential statuses individuals can take. The edges represent transitions that can occur, and we weight the edges by the rate. In the SEIR case we would need the graph H to have edges 'E'→'I' and 'I'→'R'. The edges would be weighted by the transition rates. Note H need not have a node 'S' because susceptible nodes do not change status on their own.

Examples
We first demonstrate a stochastic simulation of a simple contagion with an SEIR example.
To demonstrate additional flexibility we allow some individuals to have a higher rate of transitioning from 'E' to 'I' and some partnerships to have a higher transmission rate. This is done by adding weights to the contact network G which scale the rates for those individuals or partnerships. The documentation discusses other ways we can allow for heterogeneity in transition rates. Note that this process is guaranteed to terminate, so we can set tmax to be infinite. Processes which may not terminate will require a finite value. The default is 100.
import EoN import networkx as nx from collections import defaultdict import matplotlib.pyplot as plt import random N = 100000 print('generating graph G with {} nodes'.format(N)) G = nx.fast_gnp_random_graph(N, 5./(N-1)) #We add random variation in the rate of leaving exposed class #and in the partnership transmission rate. #There is no variation in recovery rate. node_attribute_dict = {node: 0.5+random.random() for node in G.nodes()} edge_attribute_dict = {edge: 0.5+random.random() for edge in G.edges()} nx.set_node_attributes(G, values=node_attribute_dict, name='expose2infect_weight') nx.set_edge_attributes(G, values=edge_attribute_dict, name='transmission_weight') # #These individual and partnership attributes will be used to scale #the transition rates. When we define \texttt{H} and \texttt{J}, we provide the name #of these attributes. Interaction between diseases can lead to interesting effects [6,3,11,13]. Now we show two cooperative SIR diseases. In isolation, each of these diseases would fail to start an epidemic. However, together they can, and sometimes they exhibit interesting osillatory behavior. To help stimulate the oscillations, we start with an asymmetric initial condition, though oscillations can be induced purely by stochastic effects for smaller initial conditions. To the best of our knowledge, this oscillatory behavior has not been studied previously.

Complex contagions
Complex contagions are implemented through Gillespie complex contagion which allows a user to specify the rules governing a relatively arbitrary complex contagion. The one criteria we note is that there is no memory -an individual will change from one status to another based on the current statuses of its neighbors, and not based on previous interactions with some neighbors who may have since changed status.
In the Gillespie implementation, we need a user-defined function which calculates the rate at which u will change status (given knowledge about the current state of the system) and another user-defined function which chooses the new status of u given that it is changing status. We finally need a user-defined function that will determine which other nodes have their rate change due to u's transition. By knowing the rates of all nodes the Gillespie algorithm can choose the time of the next transition and which node transitions. Then it finds the new state, and finally it calculates the new rates for all nodes affected by the change.
Once these functions are defined, the Gillespie algorithm is able to perform the complex contagion simulation.

Example
Previous work [16] considered a dynamic version of the Watts Threshold Model [23] spreading through clustered and unclustered networks. The Watts Threshold Model is like an SI model, except that nodes have a threshold and must have more than some threshold number of infected partners before becoming infected. The dynamic model in [16] assumed that nodes transmit independently of one another, and a recipient accumulates transmissions until reaching a threshold and then switches status. An individual v can only transmit once to a partner u. Because individuals cannot transmit to the same partner more than once it becomes nontrivial to implement this in a way consistent with the memoryless criterion.
Here we use another dynamic model that yields the same final state. Once a node has reached its threshold number of infected partners, it transitions at rate 1 to the infected state. The dynamics are different, but it can be proven that the final states in both models are identical and follow deterministically from the initial condition. The following will produce the equivalent of Fig. 2a of [16] for our new dynamic model. In that Figure, the threshold was 2.
import networkx as nx import EoN import numpy as np import matplotlib.pyplot as plt from collections import defaultdict def transition_rate(G, node, status, parameters): '''This function needs to return the rate at which \texttt{node} changes status. return 'I' def get_influence_set(G, node, status, parameters): '''this function needs to return a set containing all nodes whose rates might change because \texttt{node} has just changed status. That is, which nodes might \texttt{node} influence?
For our models the only nodes a node might affect are the susceptible neighbors. G will be a random clustered network [12,19], with the same degree distribution as before. If we use a different range of values of rho, such as for rho in np.linspace (1./80, 5./80, 8): this will produce a figure similar to Fig. 8a of [16]. Note that the threshold initial size required to trigger a cascade is smaller in this clustered network.

Visualization & Analysis
By default the simulations return numpy arrays providing the number of individuals with each state at each time. However if we set a flag return full data=True, then the simulations return a Simulation Investigation object. With the Simulation Investigation object, there are methods which allow us to reconstruct all details of the simulation. We can know the exact status of each individual at each time, as well as who infected whom.
There are also methods provided to produce output from the Simulation Investigation object. These allow us to produce a snapshot of the network at a given time. By default the visualization also includes the time series (e.g., S, I, and R) plotted beside the network snapshot. These time series plots can be removed, or replaced by other time series, for example we could plot multiple time series in the same axis, or time series generated by one of the differential equations models. With appropriate additional packages needed for matplotlib's animation tools, the software can produce animations as well.
For SIR outbreaks, the Simulation Investigation object includes a transmission tree. For SIS and simple contagions, it includes a directed multigraph showing the transmissions that occurred (this may not be a tree). However for complex contagions, we cannot determine who is responsible for inducing a transition, so the implementation does not provide a transmission tree. The transmission tree is useful for constructing synthetic phylogenies as in [18].
Example -a snapshot of dynamics and a transmission tree for SIR disease.
Using the tools provided, it is possible to produce a snapshot of the spreading process at a given time as well as an animation of the spread. We consider SIR disease spreading in the Karate Club network [24].
Example -Visualizing dynamics of SIR disease with vaccination.
We finally consider an SIRV disease, that is an SIR disease with vaccination. As the disease spreads susceptible individuals get vaccinated randomly, without regard for the status of their neighbors.
To implement this with EoN, we must use Gillespie simple contagion.
We provide an animation showing the spread. To make it easier to visualize, we use a lattice network.
import networkx as nx import EoN import matplotlib.pyplot as plt from collections import defaultdict print('generating graph G') G = nx.grid_2d_graph(100,100) #each node is (u,v) where 0<=u,v<=99 #we'll initially infect those near the middle This will create an mp4 file animating previous display over all calculated times. Depending on the computer installation, extra args will need to be modified.

Discussion
EoN provides a number of tools for studying infectious processes spreading in contact networks. The examples given here are intended to demonstrate the range of EoN, but they represent only a fraction of the possibilities.

Related Packages
There are several alternative software packages that allow for simulation of epidemics on networks. Here we briefly review some of these.

epydemic
Epydemic is a python package that can simulate SIS and SIR epidemics in networks. It is also built on networkx. It can handle both discrete-time simulations or continuous-time Markovian simulations for which it uses a Gillespie-style algorithm. It can handle more processes than just SIS or SIR disease. In fact it can handle any model which can be simulated using the EoN.simple contagion. The documentation is available at https://pyepydemic.readthedocs.io/en/latest/.

Graph-tool
Graph-tool [20] is a python package that serves as an alternative to networkx. Many of its underlying processes are written in C++, so it is often much faster than networkx.
Graph-tool has a number of built-in dynamic models, including the SIS, SIR, and SIRS models. The disease models are currently available only in discrete-time versions.
The documentation for these disease models is available at https://graph-tool.skewed.de/static/doc/dynamics.html.

EpiModel
EpiModel [9] is an R package that can handle SI, SIS, and SIR disease spread. It is possible to extend EpiModel to other models. EpiModel is built around the StatNet package. More details about EpiModel are available at https://www.epimodel.org/.