EoN.Gillespie_simple_contagion

EoN.Gillespie_simple_contagion(G, spontaneous_transition_graph, nbr_induced_transition_graph, IC, return_statuses, tmin=0, tmax=100, spont_kwargs=None, nbr_kwargs=None, return_full_data=False, sim_kwargs=None)[source]

Performs simulations for epidemics, allowing more flexibility than SIR/SIS.

This does not handle complex contagions. It assumes that when an individual changes status either he/she has received a “transmission” from a single neighbor or he/she is changing status independently of any neighbors. So this is like SIS or SIR. Other examples would be like SEIR, SIRS, etc

There is an example below demonstrating an SEIR epidemic.

We allow for nodes to undergo two types of transitions. They can be:

  • spontaneous - so a node of status A becomes B without any neighbor’s influence
  • neighbor-induced - so an edge between status A and status B nodes suddenly becomes an edge between a status A and a status C node because the status B node changes status due to the existence of the edge. (in principle, both could change status, but the code currently only allows the second node to change status).

Both types of transitions can be represented by weighted directed graphs. We describe two weighted graphs whose nodes represent the possible statuses and whose edges represent possible transitions of statuses.

  • The spontaneous transitions can be represented by a graph whose nodes are the possible statuses and an edge from 'A' to 'B' represent that in the absence of any influence from others an indivdiual of status 'A' transitions to status 'B' with default rate given by the weight of the edge in this spontaneous transition graph. The rate may be modified by properties of the nodes in the contact network G.
  • The neighbor-induced transitions can be represented by a “transitions graph” whose nodes are length-2 tuples. The first entry represents the first individual of a partnership and the second represents the second individual. only the second individual changes status. An edge in the transitions graph from the node ('A', 'B') to the node ('A', 'C') represents that an 'AB' partnership in the contact network can cause the second individual to transition to status 'C'. The weight of the edge in the represents the default transition rate. The rate may be modified by properties of the nodes or the edge in the contact network G.

[for reference, if you look at Fig 4.3 on pg 122 of Kiss, Miller & Simon the graphs for SIS would be:

spontaneous_transition_graph: 'I'-> 'S' with the edge weighted by gamma and

nbr_induced_transition_graph: ('I', 'S') -> ('I', 'I') with the edge weighted by tau.

For SIR they would be:

spontaneous_transition_graph: 'I'->``’R’`` with weight gamma and

nbr_induced_transition_graph: ('I', 'S') -> ('I', 'I') with rate tau. ]

These graphs must be defined and then input into the algorithm.

It is possible to weight edges or nodes in the contact network G (that is, not the 2 directed networks defined above, but the original contact network) so that some of these transitions have different rates for different individuals/partnerships. These are included as attributes in the contact network. In the most general case, the transition rate depends on some function of the attributes of the nodes or edge in the contact network G.

There are two ways we introduce individual or pair-level heterogeneity in the population. The first way is through introducing weights to individuals or their contacts which simply multiply the default rate. The second is through including a function of nodes or pairs of nodes which will explicitly calculate the scaling factor, possibly including more information than we can include with the weights. For a given transition, you can use at most one of these.

  • We first describe examples of weighting nodes/edges in the population

So for the SIR case, if some people have higher recovery rate, we might define a node attribute 'recovery_weight' for each node in G, and the recovery would occur with rate G.nodes[node][‘recovery_weight’]*gamma. So a value of 1.1 would correspond to a 10% increased recovery rate. Since I don’t know what name you might choose for the weight label as I write this algorithm, in defining the spontaneous transition graph (H), the 'I'->``’R’`` edge would be given an attribute 'weight_label' so that H.adj['I']['R']['weight_label'] = 'recovery_weight'. If you define the attribute 'weight_label' for an edge in H, then it will be assumed that every node in G has a corresponding weight. If no attribute is given, then it is assumed that all transitions happen with the original rate.

We similarly define the weight_labels as edge attributes in the neighbor-induced transition graph. The edges of the graph 'G' have a corresponding G[u,v]['transmission_weight']

  • Alternately we might introduce a function.

So for the SIR case if the recovery rate depends on two attributes of a node (say, age and gender), we define a function rate_function(G,node) which will then look at G.nodes[node][‘age’] and G.nodes[node][‘gender’] and then return a factor which will be multiplied by gamma to give the recovery rate. Similarly, if we are considering a neighbor induced transition and the rate depends on properties of both nodes we define another function rate_function(G,u,v) which may use attributes of u or v or the edge to find the appropriate scaling factor.

Arguments:
G NetworkX Graph
The underlying contact network. If G is directed, we assume that “transmissions” can only go in the same direction as the edge.
spontaneous_transition_graph Directed networkx graph

The nodes of this graph are the possible statuses of a node in G. An edge in this graph is a possible transition in G that occurs without any influence from neighbors.

An edge in this directed graph is labelled with attributes

  • 'rate' [a number, the default rate of the transition]
  • 'weight_label' (optional) [a string, giving the label of
    a node attribute in the contact network G that scales the transition rate]
  • 'rate_function' (optional not combinable with 'weight_label'
    for some edge.) [a user-defined function of the contact network and node that will scale the transition rate. This cannot depend on the statuses of any nodes - we must be able to calculate it once at the beginning of the process.] It will be called as rate_function(G, node,**spont_kwargs) where spont_kwargs is described below.

Only one of 'weight_label' and 'rate_function' can be given.

In the description below, let’s use
  • rate = spontaneous_transition_graph.adj[Status1][Status2]['rate']
  • weight_label = spontaneous_transition_graph.adj[Status1][Status2]['weight_label']
  • rate_function = spontaneous_transition_graph.adj[Status1][Status2]['rate_function']

For a node u whose status is Status1, the rate at which u transitions to Status2 is

  • rate if neither weight_label nor rate_function is defined.
  • rate*G.nodes[u][weight_label] if weight_label is defined.
  • rate*rate_function(G, u, **spont_kwargs) if rate_function is
    defined.

So for example in the case of an SIR disease, this would be a graph with an isolated node 'S' and an edge from node 'I' to 'R' with rate equal to gamma (the recovery rate). It would not actually be necessary to have the node 'S'.

Note that the rate_function version is more general, and in principle we don’t need the weight_label option. It’s here for backwards compatibility and general simplicity purposes.

nbr_induced_transition_graph Directed networkx graph

The nodes of this graph are tuples with possible statuses of nodes at the end of an edge. The first node in the tuple is the node that could be affecting the second. So for example for the SIR model we would expect a node ('I', 'S') with an edge to ('I', 'I').

An edge in this directed graph is labelled with attributes

  • 'rate' [a number, the default rate of the transition]
  • 'weight_label' (optional) [a string, giving the label of
    an edge* attribute in the contact network G that scales the transition rate]
  • 'rate_function' (optional not combinable with 'weight_label'
    for some edge.) [a user-defined function of the contact network and source and target nodes that will scale the transition rate. This cannot depend on the statuses of any nodes - we must be able to calculate it once at the beginning of the process. It will be called as rate_function(G, source, target, *nbr_kwargs)]

Only one of 'weight_label' and 'rate_function' can be given.

In the description below, let’s use

  • rate = spontaneous_transition_graph.adj[Status1][Status2]['rate']
  • weight_label = spontaneous_transition_graph.adj[Status1][Status2]['weight_label']
  • rate_function = spontaneous_transition_graph.adj[Status1][Status2]['rate_function']

If the transition is (A,B) to (A,C) and we have a source node of status A joined to a target node of status B then the target node transitions from B to C with rate

  • rate if neither weight_label nor rate_function is defined.
  • rate*G.adj[source][target][weight_label] if weight_label is defined.
  • rate*rate_function(G, source, target, **nbr_kwargs) if rate_function is defined

So for example in the case of an SIR disease with transmission rate tau, this would be a graph with an edge from the node ('I','S') to ('I', 'I'). The attribute 'rate' for the edge would be tau.

IC dict
states the initial status of each node in the network.
return_statuses list or other iterable (but not a generator)
The statuses that we will return information for, in the order we will return them.
tmin number (default 0)
starting time
tmax number (default 100)
stop time
spont_kwargs dict or None (default None)
Any parameters which might be needed if the user has defined a rate function for a spontaneous transition. If any of the spontaneous transition rate functions accepts these, they all need to (even if not used). It’s easiest to define the function as def f(…, **kwargs)
nbr_kwargs dict or None (default None)
Any parameters which might be needed if the user has defined a rate function for a neighbor-induced transition. If any of the neighbor-induced transition rate functions accepts these, they all need to (even if not used). It’s easiest to define the function as def f(…, **kwargs)
return_full_data boolean
Tells whether to return a Simulation_Investigation object or not
sim_kwargs keyword arguments
Any keyword arguments to be sent to the Simulation_Investigation object Only relevant if return_full_data=True
Returns:
(times, status1, status2, …) tuple of numpy arrays
first entry is the times at which events happen. second (etc) entry is an array with the same number of entries as times giving the number of nodes of status ordered as they are in return_statuses
SAMPLE USE:

This does an SEIR epidemic. It treats the nodes and edges as weighted. So some individuals have a higher E->I transition rate than others, and some edges have a higher transmission rate than others. The recovery rate is taken to be the same for all nodes.

There are more examples in the online documentation at :ref:simple-contagion-section.

import EoN
import networkx as nx
from collections import defaultdict
import matplotlib.pyplot as plt
import random

N = 100000
G = nx.fast_gnp_random_graph(N, 5./(N-1))

#they will vary in the rate of leaving exposed class.
#and edges will vary in transition 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')


H = nx.DiGraph()
H.add_node('S')
H.add_edge('E', 'I', rate = 0.6, weight_label='expose2infect_weight')
H.add_edge('I', 'R', rate = 0.1)

J = nx.DiGraph()
J.add_edge(('I', 'S'), ('I', 'E'), rate = 0.1, weight_label='transmission_weight')
IC = defaultdict(lambda: 'S')
for node in range(200):
    IC[node] = 'I'

return_statuses = ('S', 'E', 'I', 'R')

t, S, E, I, R = EoN.Gillespie_simple_contagion(G, H, J, IC, return_statuses,
                                        tmax = float('Inf'))

plt.semilogy(t, S, label = 'Susceptible')
plt.semilogy(t, E, label = 'Exposed')
plt.semilogy(t, I, label = 'Infected')
plt.semilogy(t, R, label = 'Recovered')
plt.legend()

plt.savefig('SEIR.png')