Epidemic Models on Networks, in depth: SIR and SIS on graphs, the excess-degree threshold, percolation and event-driven simulation

SIR and SIS on a graph

Each node is in one compartment. In SIR a node moves Susceptible to Infectious to Removed and never returns; removed means recovered with immunity, isolated or dead. In SIS an infectious node recovers straight back to susceptible, which models infections without lasting immunity, such as many bacterial infections or a reinstallable piece of malware. Variants add an Exposed state (SEIR) or waning immunity (SIRS), but the network effects show up already in SIR and SIS.

In continuous time the network model has two rates. Each edge between an infectious node and a susceptible node transmits at rate beta, so a susceptible node with three infectious neighbours is infected at total rate 3 beta. Each infectious node recovers at rate gamma, so its infectious period is exponential with mean 1/gamma. Note the units: beta here is per edge, not per person as in the well-mixed ODE, which is a common source of confusion when people port parameters between the two worlds.

The key derived quantity for SIR is the transmissibility T, the probability that an infected node passes the infection along a given edge before it recovers. With exponential transmission and recovery times the two clocks race, and T = beta / (beta + gamma). With a fixed infectious period tau instead, T = 1 - exp(-beta tau). The difference matters: two models with the same mean infectious period but different period distributions give different T and therefore different outbreaks.

Why the threshold depends on the second moment

Follow an infection down a chain. A node reached by an edge did not get there at random: it was picked by following an edge, and a node with degree k is k times more likely to be at the end of a random edge than a node of degree 1. Once infected, it can pass the infection to k - 1 new neighbours, because one edge leads back to its infector. The average number of onward edges, the mean excess degree, is therefore (⟨k²⟩ - ⟨k⟩) / ⟨k⟩, where ⟨k⟩ is the mean degree and ⟨k²⟩ the mean squared degree.

Each onward edge fires with probability T, so on a locally tree-like random graph the basic reproduction number is R0 = T (⟨k²⟩ - ⟨k⟩) / ⟨k⟩. A large outbreak is possible only when R0 is above 1, which gives the critical transmissibility T_c = ⟨k⟩ / (⟨k²⟩ - ⟨k⟩). This is the bond percolation threshold of the configuration model. With a fixed infectious period the SIR outbreak is exactly the cluster of the starting node after every edge is kept independently with probability T. With a variable period, such as the exponential one, a node's outgoing edges share one recovery clock and are correlated: the mapping still gives the threshold and the size of a large outbreak, but not the probability that one happens. That mapping is the most useful fact in the subject, because percolation results about cluster sizes become epidemic results about outbreak sizes.

For SIS the heterogeneous mean-field approximation, which groups nodes by degree, gives a threshold beta/gamma above ⟨k⟩ / ⟨k²⟩. A finer approximation that keeps the actual adjacency matrix, the quenched mean field, gives a threshold of 1 / Λ_max, where Λ_max is the largest eigenvalue of the adjacency matrix. Both say the same thing: heterogeneity lowers the threshold. For a scale-free degree distribution with exponent at or below 3, ⟨k²⟩ grows without bound as the network grows, so the threshold tends to zero. Pastor-Satorras and Vespignani pointed this out in 2001; it explains why computer viruses with low transmissibility persisted on the internet.

Worked example: same mean degree, different epidemic

Take two networks of 10,000 nodes, both with mean degree ⟨k⟩ = 4. Network A is an Erdős–Rényi random graph, so degrees are Poisson and ⟨k²⟩ = ⟨k⟩² + ⟨k⟩ = 20. Network B has a heavy tail with the same mean but ⟨k²⟩ = 60. Let beta = 0.1 per day per edge and gamma = 0.2 per day, so the mean infectious period is 5 days and T = 0.1 / 0.3 ≈ 0.333.

For network A the final outbreak size has a closed form. Bond-percolating a Poisson graph with probability T leaves a Poisson graph with mean ⟨k⟩T = 1.33, and the giant cluster fraction S solves S = 1 - exp(-1.33 S). Iterating from S = 0.5 converges to S ≈ 0.454: if a large outbreak happens, about 45% of nodes are eventually infected. The chance that a single seed sparks a large outbreak is lower. With a fixed infectious period it would also be 0.454, but with exponential periods some nodes recover almost at once and others stay infectious for a long time, and solving the branching process with that variation gives about 0.32. So roughly two thirds of introductions die out by chance even though R0 is above 1. That stochastic die-out is invisible in ODE models and it is why one simulation run tells you very little.

The vaccination row shows the policy cost of heterogeneity. Vaccinating a fraction v at random removes nodes uniformly, scaling R0 by (1 - v), so the critical coverage is 1 - 1/R0. On network A, 25% coverage brings R0 to exactly 1. On network B the same pathogen needs about 79%. Same mean contact count, very different public health problem, and the difference sits entirely in ⟨k²⟩.

Exact event-driven simulation

Analytical thresholds assume a locally tree-like graph and large size. Real networks have clustering, finite size and correlations, so you simulate. The exact way to simulate Markovian SIR is event-driven: when a node becomes infectious, draw its recovery time, then draw a transmission time for each neighbour and keep it only if it comes before recovery and before any earlier scheduled infection of that neighbour. A priority queue processes events in time order. There is no time step, so there is no discretisation error, and the cost is about O(E log N).

import heapq
import random

def sir_event_driven(adj, beta, gamma, seeds, t_max=float("inf"), rng=random):
    # adj: dict node -> list of neighbours. Returns infection and recovery times.
    infected_at, recovered_at = {}, {}
    pending = {}                       # node -> earliest scheduled infection time
    queue = [(0.0, "inf", s) for s in seeds]
    for s in seeds:
        pending[s] = 0.0
    heapq.heapify(queue)
    while queue:
        t, kind, u = heapq.heappop(queue)
        if t > t_max:
            break
        if kind == "rec":
            recovered_at[u] = t
            continue
        if u in infected_at or pending.get(u) != t:
            continue                   # stale event, already infected earlier
        infected_at[u] = t
        t_rec = t + rng.expovariate(gamma)
        heapq.heappush(queue, (t_rec, "rec", u))
        for v in adj[u]:
            if v in infected_at:
                continue
            t_inf = t + rng.expovariate(beta)
            if t_inf < t_rec and t_inf < pending.get(v, float("inf")):
                pending[v] = t_inf
                heapq.heappush(queue, (t_inf, "inf", v))
    return infected_at, recovered_at

def outbreak_stats(adj, beta, gamma, runs=500, threshold=0.05, rng=random):
    n, sizes = len(adj), []
    for _ in range(runs):
        seed = rng.choice(list(adj))
        inf, _ = sir_event_driven(adj, beta, gamma, [seed], rng=rng)
        sizes.append(len(inf) / n)
    big = [s for s in sizes if s >= threshold]
    return len(big) / runs, (sum(big) / len(big) if big else 0.0)

The second function reports two numbers, and you should always report both: the probability of a large outbreak and the mean size conditional on one occurring. Averaging all runs together mixes a bimodal distribution into a meaningless middle value. On network A from the worked example, expect a large-outbreak probability near 0.32 and a conditional size near 0.45, with finite-size scatter; a mismatch in the second number points to a bug. For SIS use the Gillespie algorithm instead: keep the total rate gamma times the infected count plus beta times the number of S-I edges, draw the next event time from an exponential with that total, and pick the event in proportion to its rate.

Interventions that use structure

Because hubs drive R0 through ⟨k²⟩, removing hubs is far more efficient than removing random nodes. Targeted vaccination by degree can stop an outbreak on a heavy-tailed network with a small fraction of the coverage random vaccination needs, but it requires knowing the whole graph, which public health agencies almost never do. Acquaintance immunization, proposed by Cohen, Havlin and ben-Avraham in 2003, gets most of the benefit with only local information: pick random people, and vaccinate a random contact of each. A random contact is reached by following an edge, so it is biased toward high degree, which is the same friendship-paradox argument that produced the excess degree.

def acquaintance_targets(adj, fraction, rng=random):
    # Vaccinate random neighbours of random nodes until the budget is spent.
    budget, chosen, nodes = int(fraction * len(adj)), set(), list(adj)
    while len(chosen) < budget:
        u = rng.choice(nodes)
        if adj[u]:
            chosen.add(rng.choice(adj[u]))
    return chosen

def remove_nodes(adj, removed):
    return {u: [v for v in nbrs if v not in removed]
            for u, nbrs in adj.items() if u not in removed}

Evaluate strategies by running the simulator on the graph with vaccinated nodes removed and comparing large-outbreak probability and size at equal budgets. Centrality measures beyond degree can do better still on clustered graphs, since they account for position, not just edge count; the same ranking problem appears in reverse in influence maximization, where you choose seeds to maximise spread rather than blockers to minimise it.

Failure modes

• Porting a well-mixed beta. An ODE beta is a per-person contact rate times transmission probability. Dropping it into a per-edge network model inflates or deflates R0 by a factor of the mean degree. Recalibrate from an observed growth rate or a measured secondary attack rate.
• Fixed time steps. A discrete-time simulation with a large step lets a node be infected and transmit in the same step, overestimating speed. Use event-driven or Gillespie simulation, or shrink the step until results stop moving.
• Trusting tree-like formulas on clustered graphs. In a household of four, many edges point back to already infected members, so the excess-degree formula overestimates R0. Use simulation or message-passing methods that handle clustering.
• Sampled networks. A contact survey truncated at, say, ten contacts per person cuts off the tail, and ⟨k²⟩ is dominated by the tail. Under-sampled hubs make the threshold look higher than it is.
• Static graphs for temporal contacts. An edge that exists only on Tuesdays cannot carry infection on Wednesday. Aggregating a temporal network into a static graph creates paths that do not respect time ordering and overstates spread.
• Reporting a single run or the overall mean. Stochastic die-out makes outcomes bimodal. Report outbreak probability, conditional size and their spread over hundreds of runs.

Trade-offs between methods

What to do next

1. Write down your model's units: is beta per edge or per person, and is the infectious period exponential or fixed? Compute T accordingly.
2. Estimate ⟨k⟩ and ⟨k²⟩ from your contact data, check how the tail was sampled, and compute the tree-like R0 and T_c as a baseline.
3. Implement the event-driven simulator above, validate it on an Erdős–Rényi graph against S = 1 - exp(-⟨k⟩T S), then run it on your real graph.
4. Report outbreak probability and conditional size over at least several hundred runs, with confidence intervals.
5. Compare random, degree-targeted and acquaintance vaccination at equal budgets, and add a centrality-based strategy if your graph is clustered.
6. If contacts change over time, rebuild the simulation on a time-ordered edge list before trusting any intervention ranking.

Keep learning: influence maximization turns the same spreading model into a seed-selection problem, eigenvector centrality and PageRank give rankings for targeting, and Louvain community detection finds the clusters that slow spread down.