Tutorial Seven — Exercises

Epidemic Spreading

The following helper function implements a discrete-time SI model on a NetworkX graph. You may use and adapt it throughout this tutorial.

import numpy as np

def run_si(G, beta, seed_node, n_steps):
    """
    Run a discrete-time SI model on graph G.

    Parameters
    ----------
    G         : NetworkX graph
    beta      : infection probability per edge per step
    seed_node : node to place in Infected state at t=0
    n_steps   : number of time steps to simulate

    Returns
    -------
    S, I : arrays of length n_steps+1 with fraction of nodes in each state
    """
    nodes = list(G.nodes())
    n = len(nodes)
    state = {v: 'S' for v in nodes}
    state[seed_node] = 'I'

    S, I = [1 - 1/n], [1/n]

    for _ in range(n_steps):
        new_state = state.copy()
        for v in nodes:
            if state[v] == 'I':
                for u in G.neighbors(v):
                    if state[u] == 'S' and np.random.rand() < beta:
                        new_state[u] = 'I'
        state = new_state
        s = sum(state[v] == 'S' for v in nodes) / n
        S.append(s)
        I.append(1 - s)

    return np.array(S), np.array(I)

1. Load the network at http://www.networkatlas.eu/exercises/20/1/data.txt. Run the SI model 10 times for each of β ∈ {0.05, 0.1, 0.2}, seeding a randomly chosen node at each run.

(a) For each value of β, plot the mean fraction of infected nodes as a function of time step, and report the average time step at which 80% of the network is infected.

(b) Examine the shape of the infection curves. What does the early-time growth tell you about the dynamics? What does the curve converge to, and why?

2. Load the network at http://www.networkatlas.eu/exercises/20/2/data.txt. Run the same SI model with β ∈ {0.05, 0.1, 0.2} and the same seeding procedure as in Question 1.

(a) One of the two networks is an Erdős–Rényi random graph; the other has a power-law degree distribution. Based on the infection curves and the time to 80% infection, can you identify which is which?

(b) For each network, compute ⟨k²⟩/⟨k⟩. How does this quantity help explain the difference in spreading speed you observed? (Hint: consider what this ratio tells you about the epidemic threshold.)

3. Extend the SI model to an SIS model, in which infected nodes recover and become susceptible again at rate μ at each time step. Run the SIS model for 100 steps on both networks, with β = 0.2 and μ ∈ {0.05, 0.1, 0.2}. Seed a single randomly chosen node as the initial infected node. Plot the fraction of infected nodes over time for each combination.

(a) For which combinations of network and μ do you observe an endemic state (a sustained non-zero fraction of infected nodes)?

(b) The epidemic threshold for the SIS model on a network is:

\[\frac{\beta}{\mu} > \frac{1}{\lambda_1}\]

where λ₁ is the largest eigenvalue of the adjacency matrix. Compute λ₁ for each network and determine the theoretical threshold value of β/μ. Does it agree with your simulations? (Note: randomness may affect individual runs — repeat the experiment several times.)

4. Extend the SI model to an SIR model, in which recovered nodes gain permanent immunity and move to a Removed state at rate μ. The following function implements this model with an explicit seed node:

def run_sir(G, beta, mu, seed_node, n_steps):
    """
    Run a discrete-time SIR model on graph G.

    Parameters
    ----------
    G         : NetworkX graph
    beta      : infection probability per edge per step
    mu        : recovery probability per step
    seed_node : node to place in Infected state at t=0
    n_steps   : number of time steps to simulate

    Returns
    -------
    S, I, R : arrays of length n_steps+1 with fraction of nodes in each state
    """
    nodes = list(G.nodes())
    n = len(nodes)
    state = {v: 'S' for v in nodes}
    state[seed_node] = 'I'

    S, I, R = [1 - 1/n], [1/n], [0.0]

    for _ in range(n_steps):
        new_state = state.copy()
        for v in nodes:
            if state[v] == 'I':
                if np.random.rand() < mu:
                    new_state[v] = 'R'
                else:
                    for u in G.neighbors(v):
                        if state[u] == 'S' and np.random.rand() < beta:
                            new_state[u] = 'I'
        state = new_state
        s = sum(state[v] == 'S' for v in nodes) / n
        i = sum(state[v] == 'I' for v in nodes) / n
        r = 1 - s - i
        S.append(s); I.append(i); R.append(r)

    return np.array(S), np.array(I), np.array(R)

(a) Run the SIR model for 400 steps on both networks, with β = 0.2 and μ ∈ {0.01, 0.02, 0.04}. Plot the fraction of nodes in the Removed state over time. How quickly does the epidemic run its course for each parameter combination?

(b) Does it matter which node you seed? For β = 0.2, μ = 0.02, compare SIR dynamics when seeding: (i) the highest-degree node, (ii) the highest-PageRank node, and (iii) a randomly chosen node. Run each condition 10 times and plot the mean fraction removed over time with 25th–75th percentile bands. Do the three seeding strategies lead to measurably different outcomes?

Discussion question: Would you expect the advantage of targeted seeding to be larger or smaller on the power-law network than on the Erdős–Rényi network? Use your findings from Questions 2 and 3 to support your answer.