Network models

NetworkProcess spreads infection over a fixed contact network. Each node is a person, and each graph edge is a potential route of transmission. Transmission is driven by a contact rate rather than a single coin flip per edge: along each edge the time from an infector becoming infectious to its next infectious contact with a neighbour is drawn from a contact-interval distribution (the kernel), and transmission happens only while the infector is still inside its infectious window. Because that hazard races the infector's recovery or isolation, shortening the infectious window genuinely curtails onward spread — something the earlier coin-flip-per-edge model could not express.

Each node can be infected once, and the graph does the job the offspring distribution does in BranchingProcess: it sets who can infect whom. Node attributes and the clinical timeline behave as they do elsewhere.

Defining a network

The process takes the contact network directly: pass an adjacency list, where adjacency[i] holds the nodes connected to node i, and the contact-interval kernel — a Distributions.jl distribution shared by every edge. The disease natural history is a progression of Transitions attached with a ModelSpec: a terminal removal transition sets how long a node stays infectious once infected, and an optional latent-period transition delays the start of that window. The kernel times contacts from the window's start state, derived from the progression (:infectious when a latent period produces it, otherwise :infection).

using EpiBranch
using EpiNetwork
using Distributions
using StableRNGs

# Twenty households of four, linked into a ring of households by one
# bridge edge between consecutive households.
function household_ring(n_households, household_size)
    n = n_households * household_size
    adj = [Int[] for _ in 1:n]
    for h in 0:(n_households - 1)
        members = (h * household_size + 1):(h * household_size + household_size)
        for i in members, j in members
            i != j && push!(adj[i], j)
        end
        # bridge to the next household
        a = h * household_size + 1
        b = (mod(h + 1, n_households)) * household_size + 1
        push!(adj[a], b)
        push!(adj[b], a)
    end
    return [sort(unique(a)) for a in adj]
end

adjacency = household_ring(20, 4)
model = ModelSpec(NetworkProcess(adjacency, Exponential(3.0));
    progression = [
        Transition(:infectious; from = :infection, delay = LogNormal(1.6, 0.5)),
        Transition(:recovered; from = :infectious, delay = 7.0, terminal = true),
    ])

ModelSpec(NetworkProcess(nodes=80, edges=140, kernel=Exponential); 0 interventions, 2 transitions)

The kernel here has a mean contact interval of three days, and each node stays infectious for seven days, so most edges have time to transmit but not all do.

A weighted adjacency matrix works too. NetworkProcess(A, kernel) reads any nonzero A[i, j] as an undirected edge; the matrix marks the graph structure only, and every edge shares the kernel.

Generating a network with Graphs.jl

Building an adjacency list by hand suits small or bespoke structures, but for realistic contact networks it is easier to use the generators in Graphs.jl, the standard Julia graph library. Load Graphs.jl and pass a graph straight to NetworkProcess: each vertex becomes a node and each vertex's neighbours become its contacts.

using Graphs

# A small-world network: mostly local contacts (high clustering) with a
# few long-range links, from the Watts–Strogatz model.
g = watts_strogatz(400, 6, 0.1; rng = StableRNG(1))

model_ws = ModelSpec(NetworkProcess(g, Exponential(3.0));
    progression = [Transition(:recovered; from = :infection, delay = 7.0, terminal = true)])
state = simulate(model_ws; n_initial = 1, rng = StableRNG(3))
println("Final size: ", state.cumulative_cases, " of ", nv(g))

Final size: 400 of 400

Any generator that returns a graph works, so the structure the outbreak spreads on is a modelling choice. A few that map onto common assumptions:

• watts_strogatz(n, k, β) — small-world: local clustering with a few long-range links.

• barabasi_albert(n, k) — scale-free: a heavy-tailed degree distribution, so a minority of highly-connected nodes drive spread.

• stochastic_block_model(...) — block structure: dense within blocks and sparse between them, a natural fit for households or communities.

• euclidean_graph(n, d; cutoff) — a random geometric graph where nodes close in space are linked, so clustering emerges from proximity (the mechanism spatial outbreak-network models use); it returns the graph and the distances, so take the first element.

The kernel, progression and attributes attach exactly as before — only the source of the graph changes. Interventions attach through the infectious window: one that removes a case from transmission — Isolation — shortens that window and curtails spread (see below). ContactTracing applies too, because a node's contacts are its graph neighbours and quarantining a traced neighbour closes that neighbour's own window; see Contact tracing on a network. An intervention whose effect is purely a per-contact competing risk against the infection event, such as leaky vaccination, has no representation on the continuous-time network path and is reported with a warning rather than applied. Graphs.jl is an optional dependency: this constructor becomes available once you load Graphs.jl, and the adjacency-list and matrix constructors need nothing extra. For a directed graph, a node's out-neighbours are the contacts it can infect.

Simulating

simulate returns a SimulationState, and linelist renders it as a one-row-per-case DataFrame.

rng = StableRNG(42)
state = simulate(model; n_initial = 1, rng = rng)

println("Final outbreak size: ", state.cumulative_cases, " of ", length(adjacency))

Final outbreak size: 76 of 80

The population is the graph, so an outbreak saturates at the number of nodes instead of growing without bound. For a batch of independent runs, simulate over a sequence of seeds:

sizes = [simulate(model; n_initial = 1, rng = StableRNG(i)).cumulative_cases
         for i in 1:200]
println("Mean size: ", round(sum(sizes) / length(sizes), digits = 1))

Mean size: 56.7

Attributes belong to the node

Each node is built once, so its attributes are drawn once and stay fixed for the run. Node properties like age or risk group are part of the network and are carried into the line list.

attrs = [
    demographics(age_distribution = Uniform(0, 80)),
    clinical_presentation(incubation_period = LogNormal(1.6, 0.5)),
]

model_attrs = ModelSpec(NetworkProcess(adjacency, Exponential(3.0));
    progression = [Transition(:recovered; from = :infection, delay = 7.0, terminal = true)],
    attributes = attrs)
state = simulate(model_attrs; n_initial = 1, rng = StableRNG(7))

infected = filter(is_infected, state.individuals)
ages = [ind.state[:age] for ind in infected]
println("Cases: ", length(infected),
    "; mean age: ", round(sum(ages) / length(ages), digits = 1))

Cases: 32; mean age: 41.3

Isolation curtails onward spread

Because transmission is a hazard racing removal, shortening a node's infectious window cuts onward transmission — the feature a coin-flip-per-edge model cannot express. Isolation is part of the clinical timeline: adding an :isolated transition closes the infectious window early, so a node that isolates soon after becoming infectious contacts fewer neighbours before it stops transmitting.

Here the same fast kernel spreads through the whole ring when nothing stops it, but isolating each case a couple of days after infection holds the outbreak back.

kernel = Exponential(1.5)

baseline = ModelSpec(NetworkProcess(adjacency, kernel);
    progression = [
        Transition(:recovered; from = :infection,
            delay = (rng, ind) -> 10.0, terminal = true)])

isolating = ModelSpec(NetworkProcess(adjacency, kernel);
    progression = [
        Transition(:recovered; from = :infection,
            delay = (rng, ind) -> 10.0, terminal = true),
        Transition(:isolated; from = :infection, delay = Exponential(2.0))])

base_sizes = [simulate(baseline; n_initial = 1, rng = StableRNG(i)).cumulative_cases
              for i in 1:200]
iso_sizes = [simulate(isolating; n_initial = 1, rng = StableRNG(i)).cumulative_cases
             for i in 1:200]

println("Mean size, no isolation:   ",
    round(sum(base_sizes) / length(base_sizes), digits = 1))
println("Mean size, with isolation: ",
    round(sum(iso_sizes) / length(iso_sizes), digits = 1))

Mean size, no isolation:   80.0
Mean size, with isolation: 9.8

Contact tracing on a network

On a branching process a case's contacts are drawn fresh from an offspring distribution, so tracing reaches people who exist only as that case's offspring. On a network they are the node's graph neighbours, which is a stronger statement: neighbourhoods overlap, so the same person can be a contact of several cases, and tracing a clustered graph repeatedly finds people who have already been found.

ContactTracing composes onto NetworkProcess unchanged. A traced neighbour is quarantined at its trace time, which closes that neighbour's own infectious window if and when it is infected, so tracing acts on the same window that isolation does.

clinical = clinical_presentation(incubation_period = LogNormal(1.0, 0.3),
    prob_asymptomatic = 0.0)
iso = Isolation(onset_to_isolation_delay = Exponential(2.0), test_sensitivity = 1.0)

ws = watts_strogatz(400, 6, 0.1; rng = StableRNG(1))
build(ivs) = ModelSpec(NetworkProcess(ws, Exponential(16.0));
    progression = [Transition(:recovered; from = :infection, delay = 7.0,
        terminal = true)],
    interventions = ivs, attributes = clinical)

meansize(ivs) = sum(simulate(build(ivs); n_initial = 1,
                        rng = StableRNG(s)).cumulative_cases for s in 1:100) / 100

println("no control:              ", round(meansize(AbstractIntervention[]), digits = 1))
println("isolation:               ", round(meansize([iso]), digits = 1))
for p in (0.5, 1.0)
    ct = ContactTracing(probability = p, isolation_to_trace_delay = Exponential(1.0))
    println("isolation + $(round(Int, 100p))% tracing:  ", round(meansize([iso, ct]), digits = 1))
end

no control:              234.8
isolation:               18.2
isolation + 50% tracing:  8.0
isolation + 100% tracing:  4.1

Tracing runs when the race settles a case, which is the first point at which that case's trace time is known. It therefore reaches the neighbours that are not yet settled themselves. Tracing backwards, to the already-settled neighbour a case was infected by, is not supported on either engine.

An intervention whose effect is a competing risk against the infection event rather than a removal, such as leaky RingVaccination, still has no representation here and is reported with a warning.

Several routes at once

Everything above gives every edge the same status: one kernel, one infectious window. So anything that closes the window cuts all transmission at once — which is the wrong shape for the most ordinary control measure there is. Someone who self-isolates stops mixing in the community and goes on infecting the people they live with.

RoutedNetwork separates the edges into routes. Each is a RouteWindow with its own adjacency, kernel and set of states that end it, so isolation can cut one and leave another running. Households are cliques, community contact a sparser graph over the same people:

function households_and_community(n_households, household_size, rng)
    n = n_households * household_size
    hh = [Int[] for _ in 1:n]
    for h in 0:(n_households - 1)
        members = (h * household_size + 1):(h * household_size + household_size)
        for i in members, j in members
            i != j && push!(hh[i], j)
        end
    end
    comm = [Int[] for _ in 1:n]
    for _ in 1:(3 * n)
        a, b = rand(rng, 1:n), rand(rng, 1:n)
        if a != b && !(b in comm[a])
            push!(comm[a], b); push!(comm[b], a)
        end
    end
    return hh, comm
end

hh_adj, comm_adj = households_and_community(150, 4, StableRNG(99))

clinical2 = clinical_presentation(incubation_period = LogNormal(0.5, 0.3),
    prob_asymptomatic = 0.0)
iso2 = Isolation(onset_to_isolation_delay = Exponential(1.0), test_sensitivity = 1.0)
REM = EpiBranch.INTERVENTION_REMOVAL

# The household route differs from the community route in one tuple: whether
# isolation is allowed to end it.
routes(household_until) = [
    RouteWindow(:household; until = household_until,
        kernel = Weibull(1.5, 4.0), reach = hh_adj),
    RouteWindow(:community; until = (:recovered, REM),
        kernel = Exponential(30.0), reach = comm_adj)]

function mean_size(ws, ivs)
    m = ModelSpec(RoutedNetwork(ws);
        progression = [Transition(:recovered; from = :infection, delay = 10.0,
            terminal = true)],
        interventions = ivs, attributes = clinical2)
    sum(simulate(m; n_initial = 3, rng = StableRNG(s)).cumulative_cases
        for s in 1:80) / 80
end

println("no control:                    ",
    round(mean_size(routes((:recovered,)), AbstractIntervention[]), digits = 1))
println("isolation removes the case:     ",
    round(mean_size(routes((:recovered, REM)), [iso2]), digits = 1))
println("self-isolation at home:         ",
    round(mean_size(routes((:recovered,)), [iso2]), digits = 1))

no control:                    599.2
isolation removes the case:     214.4
self-isolation at home:         484.2

Self-isolation is markedly worse than removing the case outright, because the household route keeps running. A single-window model can only produce the second number, which is why reading window closure as "isolation" overstates what self-isolation achieves.

R is unchanged by any of this. It stays what a case would achieve if never removed; the realised figure falls out of which routes were cut.

Which contacts can be traced

Routes also differ in who a case can name. Everyone in the household can be named, but most community contacts are strangers. traceable on a route is the probability that a case can name a contact made on it. Contact tracing then reaches a named contact with its own probability, so the two multiply. A neighbour on both routes is named with the higher probability, since someone you live with can be named whether or not you also meet them elsewhere.

# Slower isolation and more community contact than above, so that tracing has
# transmission left to prevent.
iso3 = Isolation(onset_to_isolation_delay = Exponential(4.0), test_sensitivity = 1.0)
ct3 = ContactTracing(probability = 0.9, isolation_to_trace_delay = Exponential(0.5))
traced_routes(community_traceable) = [
    RouteWindow(:household; until = (:recovered, REM),
        kernel = Weibull(1.5, 4.0), reach = hh_adj),
    RouteWindow(:community; until = (:recovered, REM),
        kernel = Exponential(15.0), reach = comm_adj,
        traceable = community_traceable)]

println("isolation only:                          ",
    round(mean_size(traced_routes(1.0), [iso3]), digits = 1))
for p in (0.0, 0.5, 1.0)
    println("isolation + tracing, community traceable $p: ",
        round(mean_size(traced_routes(p), [iso3, ct3]), digits = 1))
end

isolation only:                          593.5
isolation + tracing, community traceable 0.0: 542.5
isolation + tracing, community traceable 0.5: 268.4
isolation + tracing, community traceable 1.0: 114.4

Tracing household contacts alone already helps, and the more community contacts a case can name, the more tracing prevents. Treating every route as fully traceable overstates what tracing achieves whenever much of the transmission happens between people who cannot name each other.

Community introductions

Without an external hazard the outbreak starts from the seeded index nodes and spreads only along the edges. An external_hazard adds a community force of infection, so fresh introductions appear over time; a finite obs_end on the process bounds the window over which they can arrive.

model_ext = ModelSpec(NetworkProcess(adjacency, Exponential(3.0);
        external_hazard = 0.02, obs_end = 60.0);
    progression = [Transition(:recovered; from = :infection, delay = 7.0, terminal = true)])
state = simulate(model_ext; rng = StableRNG(11))
df = linelist(state)

println("Cases: ", size(df, 1),
    "; community introductions: ", count(df.index))

Cases: 80; community introductions: 14