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
 neighborinduced  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 networkG
.  The neighborinduced transitions can be represented by a “transitions
graph” whose nodes are length2 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 networkG
.
[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 bygamma
andnbr_induced_transition_graph
:('I', 'S')
>('I', 'I')
with the edge weighted bytau
. For SIR they would be:
spontaneous_transition_graph
:'I'
>``’R’`` with weightgamma
andnbr_induced_transition_graph
:('I', 'S')
>('I', 'I')
with ratetau
. ]
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 networkG
.There are two ways we introduce individual or pairlevel 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 inG
, 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 thatH.adj['I']['R']['weight_label'] = 'recovery_weight'
. If you define the attribute'weight_label'
for an edge inH
, then it will be assumed that every node inG
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 neighborinduced transition graph. The edges of the graph
'G'
have a correspondingG[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 bygamma
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 functionrate_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 inG
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 userdefined 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)
wherespont_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 isStatus1
, the rate at whichu
transitions toStatus2
israte
if neitherweight_label
norrate_function
is defined.rate*G.nodes[u][weight_label]
ifweight_label
is defined.rate*rate_function(G, u, **spont_kwargs)
ifrate_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'
withrate
equal togamma
(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 userdefined 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 neitherweight_label
norrate_function
is defined.rate*G.adj[source][target][weight_label]
ifweight_label
is defined.rate*rate_function(G, source, target, **nbr_kwargs)
ifrate_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 betau
. 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 neighborinduced transition. If any of the neighborinduced 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 inreturn_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:
simplecontagionsection
.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./(N1)) #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')