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, policy_function=None, pf_kwargs=None, *, rng=None, return_full_data=False, sim_kwargs=None)[source]

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

fast_simple_contagion is probably a better choice for most scenarios. In my tests, it is orders of magnitude faster than this algorithm (turning days of calculations into 15 minutes in some examples).

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 edition 1 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 (e.g.) 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.

We also allow for optional policy considerations which will provide additional scalings of the rates, potentially in response to the global system state. An example might be that things like vaccination rates may change in response to the amount of infection, or contact tracing rates may have capacity constraints. This will require a global function introduced which takes the overall “natural” rates of all of the different transitions, the status of all nodes, and the graph as inputs.

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. If G has self-loops, the current implementation will break if it encounters a self-loop. I would like to modify to allow for node to alter itself.

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 for a transition.

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. A defaultdict may be useful if most nodes start off with one state.

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)

policy_function(rates, status, G, *kwargs) function (default None)

this function will adjuts the rates in response to any policy considerations/ capacity constraints or other things that may alter the rates based on the global status of the epidemic. It returns a dictionary rates whose keys are the transitions that can happen and whose values are the global rate of that transition. This will be used to determine the time of the next event and what type of transition will happen. Once we choose which transition type happens we choose which individual transitions based on the original background rates.

rng random number generator

If None, will be set to np.random.default_rng()

return_full_data boolean

Tells whether to return a Simulation_Investigation object or not

rng random number generator

If None, will be set to np.random.default_rng()

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') #everyone defaults to susceptible
for node in range(200):       #set these nodes to infected
    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')