On this lab we will consider an implementation of disease spreading over a network. We will start from SIS model and slightly modify it to SIR model.
Just to recall from the lecture how it looks like:<br> \begin{equation} \begin{cases} \cfrac{ds_i(t)}{dt} = -\beta s_i(t)\sum\limits_j A_{ij}x_j(t) + \gamma x_i(t)\\ \cfrac{dx_i(t)}{dt} = \beta s_i(t)\sum\limits_j A_{ij}x_j(t) - \gamma x_i(t) \end{cases} \\ x_i(t) + s_i(t) = 1 \end{equation} where $x_i(t)$ and $s_i(t)$ are probabilities for a node $v_i$ to be infected or susceptable.
In [ ]:
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
from numpy.linalg import eig
from scipy.integrate import odeint
%matplotlib inline
In [ ]:
# Let's start from a complete graph
n = 100
G = nx.complete_graph(n)
# G = nx.barabasi_albert_graph(n, 3)
# Get adj. matrix
A = np.array( nx.adjacency_matrix(G).todense() )
In [ ]:
# Spreading\restoring coefficient
beta, gamma = 1.3, 0.2
# Time domain
t = np.arange(0,5,0.05)
# Initial state
idx = np.random.choice(range(n), 5)
i0 = np.zeros((n,))
i0[idx] = 1
# i0 = np.random.random_integers(0,1,[n,])
z0 = np.concatenate((1-i0,i0))
In [ ]:
# System of differential equations..
def sis(z, t, A, n, beta, gamma):
return np.concatenate((
-beta * z[0:n] * A.dot(z[n:2*n]) + gamma * z[n:2*n],
beta * z[0:n] * A.dot(z[n:2*n]) - gamma * z[n:2*n]))
# ..Solved
z = odeint(sis, z0, t, (A, n, beta, gamma))
In [ ]:
# It's a bit messy, so let's see what we have got here
z.shape
In [ ]:
# Plot probs for some node
nId = 6
s = z[:,0:n]
x = z[:,n:2*n]
fig, ax = plt.subplots(1,1,figsize=(14,6))
ax.plot(s,color = 'blue')
ax.plot(x,color = 'red')
ax.set_xlabel('$t$')
ax.set_ylabel('prob')
ax.set_title('Probability for all nodes', fontsize = 15)
Hope that you remember that stuff about the correspondence between largest eigenvalue and $\frac{\gamma}{\beta}$ ratio:
In [ ]:
w,v = eig(A)
print max(w), gamma/beta
In SIR model healed population gain immunity to the infection \begin{equation} \begin{cases} \cfrac{ds_i(t)}{dt} = -\beta s_i(t)\sum\limits_j A_{ij} x_j(t)\\ \cfrac{dx_i(t)}{dt} = \beta s_i(t)\sum\limits_j A_{ij} x_j(t) - \gamma x_i(t)\\ \cfrac{dr_i(t)}{dt} = \gamma x_i(t) \end{cases} \\ x_i(t) + s_i(t) + r+i(t) = 1 \end{equation}
Adapt the code above to produce SIR model
In [ ]:
# Let's start from a complete graph
n = 100
# G = nx.complete_graph(n)
G = nx.barabasi_albert_graph(n, 3)
# Get adj. matrix
A = np.array( nx.adjacency_matrix(G).todense( ) )
# Spreading\restoring coefficient
beta, gamma = 1.3, 0.2
# Time domain
t = np.arange(0,1,0.05)
# Initial state
idx = np.random.choice(range(n), 5)
i0 = np.zeros((n,), dtype = np.float )
i0[idx] = 1
# i0 = np.random.random_integers(0,1,[n,])
z0 = np.concatenate((1-i0,i0,np.zeros((n,))))
# System of differential equations..
def sir(z, t, A, n, beta, gamma):
return np.concatenate((
- beta * z[0:n] * A.dot(z[n:2*n]) ,
beta * z[0:n] * A.dot(z[n:2*n]) - gamma * z[n:2*n],
gamma * z[n:2*n]))
# ..Solved
z = odeint(sis, z0, t, (A, n, beta, gamma))
In [ ]:
# Plot probs for some node
nId = 6
s = z[:,0:n]
x = z[:,n:2*n]
r = z[:,2*n:3*n]
fig, ax = plt.subplots(1,1,figsize=(14,6))
ax.plot(s,color = 'blue')
ax.plot(x,color = 'red')
ax.plot(r,color = 'green')
ax.set_xlabel('$t$')
ax.set_ylabel('prob')
ax.set_title('Probability for all nodes', fontsize = 15)
The stuff that you will see below is kind of simulation model of infection spreading on graph. The description of the model is the following:
In [ ]:
size = 10
# Create Grid Graph
G = nx.grid_2d_graph(size,size)
# Make node relabelling
f = {}
for v in G.nodes_iter():
f[v] = v[0]*size+v[1]
G = nx.relabel_nodes(G, f)
nx.draw_spectral(G)
In [ ]:
def simulSIS(A, timePeriod, modelParams):
# init params
initInfected = modelParams.get('initInfected', None)
p = modelParams.get('probInfect', 0.5)
upd = modelParams.get('updateInfection', True)
maxRecTime = modelParams.get('t2Recover', 2)
# init output
n = A.shape[0]
states = np.zeros([n, timePeriod+1]) # 1 = infected, 0 = susceptable
recTime = np.zeros(n,)
# set initially infected nodes
if initInfected is None:
initInfected = np.random.choice(range(n), n/2)
states[initInfected,0] = 1
else:
states[initInfected,0] = 1
recTime[initInfected] = maxRecTime + 1
# Start simulation
for t in xrange(1, timePeriod+1):
recTime = np.maximum(recTime-1,0)
states[recTime>0, t] = 1
states[recTime==0, t] = 0
curInf = np.nonzero(states[:,t])[0]
states[curInf, t] = 1
for i in curInf:
#NN = np.setdiff1d(np.nonzero(A[i,])[0], curInf)
NN = np.nonzero(A[i,])[0]
infNN = NN[np.random.rand(len(NN))<p]
states[infNN, t] = 1
recTime[infNN] = maxRecTime + 1
return states
In [ ]:
# Running model
timePeriod = 11
modelParams = {}
modelParams['t2Recover'] = 3
modelParams['initInfected'] = None
modelParams['probInfect'] = 0.2
modelParams['updateInfection'] = True
A = np.array(nx.adj_matrix(G).todense())
states = simulSIS(A, timePeriod, modelParams)
In [ ]:
states
In [ ]:
# Plotting inferction spread
pos = nx.spectral_layout(G)
fig = plt.figure(figsize=(15,10))
for t in xrange(0,timePeriod+1):
plt.subplot(3,4,t+1)
nx.draw_spectral(G,
nodelist=np.nonzero(states[:,t])[0].tolist(),
node_color = 'r')
nx.draw_spectral(G,
nodelist=np.nonzero(1-states[:,t])[0].tolist(),
node_color = 'b')
plt.title('t = {0}'.format(t+1))
In simple branchin process (tree of depth $h$ with $k$ leaths per each parent and probability $p$ to infect a neigbour) we are considerring $R_0=pk$ - the basic reproductive number of the infection.
In [ ]:
## Checking of branching branching process
G=nx.balanced_tree(3,7, create_using=nx.DiGraph())
sp = nx.shortest_path_length(G,0)
timePeriod = 10
modelParams = {}
modelParams['t2Recover'] = 1
modelParams['initInfected'] = [0]
modelParams['probInfect'] = 0.2
modelParams['updateInfection'] = False
A = np.array(nx.adj_matrix(G).todense())
states = simulSIS(A, timePeriod, modelParams)
sum(states)
In [ ]: