Interventions

Interventions in EpiBranch.jl are modelled through a competing risks framework. Potential contacts are generated by the branching process. Each contact's fate is determined by competing risks: when would transmission occur (generation time) vs when is the parent isolated (intervention time)?

There is a connection to survival analysis. The generation time CDF is the survival function of remaining potential transmission, truncated by isolation.

A policy is a layer on a model: wrap the process in a ModelSpec and pass interventions = [iso] to it, and simulate and loglikelihood read the policy from the ModelSpec. This tutorial compares several policies against the same base process, so it defines a small builder (each scenario is its own ModelSpec) and simulates each.

Without interventions

First, the baseline: a supercritical outbreak with no interventions.

using EpiBranch
using Distributions
using StableRNGs

clinical = clinical_presentation(incubation_period = LogNormal(1.5, 0.5))

# A scenario is a ModelSpec: the base process under a given policy.
scenario(interventions = AbstractIntervention[], attributes = clinical) =
    ModelSpec(BranchingProcess(Poisson(3.0), Exponential(5.0)); interventions, attributes)

rng = StableRNG(42)
results_baseline = simulate(scenario(), 200; max_cases = 500, rng = rng)
println("Containment (no interventions): $(round(containment_probability(results_baseline), digits=3))")

Containment (no interventions): 0.04

With R = 3.0, most outbreaks are not contained. Interventions are needed.

Built-in interventions

Isolation

Symptomatic, test-positive individuals are isolated after a delay from symptom onset using Isolation. Clinical state on individuals is required, set by clinical_presentation:

iso = Isolation(onset_to_isolation_delay = Exponential(2.0))

rng = StableRNG(42)
results = simulate(scenario([iso]), 200; max_cases = 500, rng = rng)
println("Containment (isolation): $(round(containment_probability(results), digits=3))")

Containment (isolation): 0.18

Effectiveness depends on how quickly isolation happens relative to the generation time. Faster isolation truncates more of the infectious period:

for d in [0.5, 2.0, 10.0]
    let iso = Isolation(onset_to_isolation_delay = Exponential(d)),
        rng = StableRNG(42)
        results = simulate(scenario([iso]), 200; max_cases = 500, rng = rng)
        println("Delay ~ Exp($d): containment = $(round(containment_probability(results), digits=3))")
    end
end

Delay ~ Exp(0.5): containment = 0.255
Delay ~ Exp(2.0): containment = 0.18
Delay ~ Exp(10.0): containment = 0.12

Leaky isolation

With post_isolation_transmission > 0, isolated individuals still transmit at a reduced rate (e.g. household contacts):

iso_leaky = Isolation(onset_to_isolation_delay = Exponential(2.0), post_isolation_transmission = 0.3)

rng = StableRNG(42)
results = simulate(scenario([iso_leaky]), 200; max_cases = 500, rng = rng)
println("Leaky isolation: $(round(containment_probability(results), digits=3))")

Leaky isolation: 0.095

Contact tracing

Contacts of isolated cases are identified using ContactTracing. With quarantine, traced contacts are isolated before symptom onset:

iso = Isolation(onset_to_isolation_delay = Exponential(2.0))
ct = ContactTracing(probability = 0.7, isolation_to_trace_delay = Exponential(1.0), quarantine_on_trace = true)

rng = StableRNG(42)
results = simulate(scenario([iso, ct]), 200; max_cases = 500, rng = rng)
println("Isolation + tracing: $(round(containment_probability(results), digits=3))")

Isolation + tracing: 0.215

Who gets traced: eligibility policies

The keyword form above uses the default policy: a contact is traced once its infector is both symptomatic and isolated. Real programmes start tracing on different events, so the infector's eligibility is set separately. Each built-in policy tests one thing about the infector:

Policy	Traces when the infector…
OnSymptomOnset	is symptomatic
OnLabConfirmation	has tested positive
OnIsolation	has been isolated
TraceEveryone / TraceNobody	always / never

Pass a policy with the positional constructor (an eligibility policy, a trace probability, and a delay distribution):

# Begin tracing as soon as the infector shows symptoms, without waiting
# for a positive test.
ct_fast = ContactTracing(OnSymptomOnset(), 0.7, Exponential(1.0))

Combine policies with the ordinary boolean operators &, |, !:

# Trace suspected OR lab-confirmed cases.
elig = OnSymptomOnset() | OnLabConfirmation()

# Trace symptomatic infectors who have not yet been isolated.
elig_gap = OnSymptomOnset() & !OnIsolation()

ct_combined = ContactTracing(elig, 0.7, Exponential(1.0))

For logic beyond the built-ins, define a policy type and one is_eligible method. It then composes with the operators like any built-in. Per-individual attributes live in infector.state:

struct SymptomaticOver65 <: EpiBranch.TraceEligibility end

function EpiBranch.is_eligible(::SymptomaticOver65, infector, contact, state)
    !EpiBranch.is_asymptomatic(infector) && get(infector.state, :age, 0) >= 65
end

elig_age = SymptomaticOver65() | OnLabConfirmation()

Asymptomatic cases and test sensitivity

Asymptomatic cases escape symptom-based surveillance. The asymptomatic fraction is set via clinical_presentation. Imperfect testing is a property of isolation — symptomatic cases are missed with probability 1 - test_sensitivity:

disease_hard = clinical_presentation(
    incubation_period = LogNormal(1.5, 0.5),
    prob_asymptomatic = 0.3,
)
iso_imperfect = Isolation(onset_to_isolation_delay = Exponential(2.0), test_sensitivity = 0.8)

rng = StableRNG(42)
results = simulate(scenario([iso_imperfect, ct], disease_hard), 200; max_cases = 500, rng = rng)
println("30% asymptomatic, 80% test sensitivity: $(round(containment_probability(results), digits=3))")

30% asymptomatic, 80% test sensitivity: 0.14

Ring vaccination

Traced contacts are vaccinated using RingVaccination. The vaccination acts as a competing risk against the contact's transmission: if the vaccination has had time to confer immunity by the contact's transmission time, it blocks transmission with probability efficacy. Requires ContactTracing in the intervention stack so contacts are identified.

iso = Isolation(onset_to_isolation_delay = Exponential(2.0))
ct = ContactTracing(probability = 0.7, isolation_to_trace_delay = Exponential(1.0))
rv = RingVaccination(efficacy = 0.8)

rng = StableRNG(42)
results = simulate(scenario([iso, ct, rv]), 200; max_cases = 500, rng = rng)
println("Iso + tracing + ring vaccination: $(round(containment_probability(results), digits=3))")

Iso + tracing + ring vaccination: 0.215

A delay between vaccination and protective immunity can be specified. If transmission occurs before immunity develops, there is no protection:

rv_delayed = RingVaccination(efficacy = 0.9, delay_to_immunity = 7.0)

rng = StableRNG(42)
results = simulate(scenario([iso, ct, rv_delayed]), 200; max_cases = 500, rng = rng)
println("With 7-day delay to immunity: $(round(containment_probability(results), digits=3))")

With 7-day delay to immunity: 0.215

We can also count the number of vaccine doses administered:

rng = StableRNG(42)
state = simulate(scenario([iso, ct, rv]); condition = 50:200, max_cases = 200, rng = rng)
n_vaccinated = count(is_vaccinated, state.individuals)
n_infected = count(is_infected, state.individuals)
println("Vaccinated: $n_vaccinated, Infected: $n_infected")

Vaccinated: 201, Infected: 200

Rings beyond direct contacts

By default tracing reaches a case's direct contacts. ContactTracing takes a depth to widen the ring: depth = 2 traces the contacts of those contacts as well, the level-2 ring that an Ebola-style protocol vaccinates around a confirmed case. Each infected, eligible case seeds a ring of radius depth, and uninfected ring members keep generating contacts for one more hop so the ring can grow past them. The same RingVaccination then vaccinates the whole ring:

ct2 = ContactTracing(probability = 0.7, isolation_to_trace_delay = Exponential(1.0), depth = 2)

rng = StableRNG(42)
state = simulate(scenario([iso, ct2, rv]); condition = 50:200, max_cases = 200, rng = rng)
doses_depth2 = count(is_vaccinated, state.individuals)
println("Doses with a level-2 ring: $doses_depth2")

Doses with a level-2 ring: 245

A wider ring reaches more people, so it costs more doses. Whether the extra reach buys extra control depends on how much transmission the direct-contact ring already caught. For a tightly traced outbreak the second ring is often mostly doses with little added containment. Compare dose counts and containment_probability across depths to see the trade-off for a given setting.

Post-exposure prophylaxis

For PEP (antivirals or antibiotics given to traced contacts), use RingVaccination with delay_to_immunity = 0.0 (the default). PEP only blocks transmission for contacts whose trace time falls before their would-be transmission time — the engine's competing-risks resolution handles this automatically:

pep = RingVaccination(efficacy = 0.9)  # delay_to_immunity defaults to 0

rng = StableRNG(42)
results = simulate(scenario([iso, ct, pep]), 200; max_cases = 500, rng = rng)
println("Iso + tracing + PEP: $(round(containment_probability(results), digits=3))")

Iso + tracing + PEP: 0.215

Mass vaccination

MassVaccination vaccinates contacts on a rolling schedule independent of tracing. Each contact becomes eligible at a time set by the eligibility_time argument; whether vaccination actually blocks their infection then depends on whether eligibility plus delay_to_immunity falls before their transmission time.

Whole population eligible on day 30 with a 14-day delay to immunity:

mv = MassVaccination(efficacy = 0.85, eligibility_time = 30.0,
    delay_to_immunity = 14.0)

rng = StableRNG(42)
results = simulate(scenario([mv]), 200; max_cases = 500, rng = rng)
println("Mass vaccination from day 30: $(round(containment_probability(results), digits=3))")

Mass vaccination from day 30: 0.065

For a rollout where individuals draw their eligibility independently (slow random rollout over weeks/months), pass a distribution:

mv_random = MassVaccination(efficacy = 0.85,
    eligibility_time = Exponential(60.0),
    delay_to_immunity = 14.0)

MassVaccination{Float64, Distributions.Exponential{Float64}, LeakyMode}(0.85, Distributions.Exponential{Float64}(θ=60.0), 14.0, LeakyMode(), :default)

For an age-stratified rollout (older first), pass a function:

mv_age = MassVaccination(
    efficacy = 0.85,
    eligibility_time = (rng, ind) -> ind.state[:age] >= 65 ? 30.0 : 90.0,
    delay_to_immunity = 14.0,
)

MassVaccination{Float64, Main.var"#2#3", LeakyMode}(0.85, Main.var"#2#3"(), 14.0, LeakyMode(), :default)

This requires the :age attribute, so add demographics to the attributes list:

attrs = [clinical, demographics(age_distribution = Uniform(0, 90))]
rng = StableRNG(42)
results = simulate(scenario([mv_age], attrs), 200; max_cases = 500, rng = rng)
println("Age-stratified rollout: $(round(containment_probability(results), digits=3))")

Age-stratified rollout: 0.07

Per-individual efficacy

efficacy accepts the same Real | Distribution | Function set as eligibility_time. A distribution gives independent per-individual draws (heterogeneous immune response); a function lets efficacy depend on individual state. Age-conditional efficacy alongside the age-stratified rollout:

mv_heterogeneous = MassVaccination(
    efficacy = (rng, ind) -> ind.state[:age] >= 65 ? 0.7 : 0.9,
    eligibility_time = (rng, ind) -> ind.state[:age] >= 65 ? 30.0 : 90.0,
    delay_to_immunity = 14.0,
)

MassVaccination{Main.var"#6#7", Main.var"#8#9", LeakyMode}(Main.var"#6#7"(), Main.var"#8#9"(), 14.0, LeakyMode(), :default)

Multi-dose schedules

For prime-and-boost (or longer) schedules, pass multiple MassVaccination (or RingVaccination) instances with different dose_labels. Each dose namespaces its per-contact state (:vaccinated_prime, :vaccinated_boost, etc.) and contributes its own competing risk; the engine composes them via the standard independent-Bernoulli product. Whichever dose's immunity has developed by the contact's transmission time contributes its blocking probability.

prime = MassVaccination(efficacy = 0.6, eligibility_time = 30.0,
    delay_to_immunity = 14.0, dose_label = :prime)
boost = MassVaccination(efficacy = 0.9, eligibility_time = 60.0,
    delay_to_immunity = 14.0, dose_label = :boost)

rng = StableRNG(42)
results = simulate(scenario([prime, boost]), 200; max_cases = 500, rng = rng)
println("Two-dose rollout: $(round(containment_probability(results), digits=3))")

Two-dose rollout: 0.045

Both single-dose and multi-dose state can be queried via the dose-suffixed keys: ind.state[:vaccinated_prime], ind.state[:vaccination_time_boost], ind.state[:vaccine_efficacy_prime], and so on.

Effort tracking

Because all contacts are stored (infected and non-infected), intervention effort is fully trackable:

rng = StableRNG(42)
state = simulate(scenario([iso, ct]); condition = 50:200, max_cases = 200, rng = rng)

total = length(state.individuals)
infected = count(is_infected, state.individuals)
traced = count(is_traced, state.individuals)
println("Contacts: $total, Infections: $infected, Traced: $traced")
println("Contacts per case: $(round(total / infected, digits=1))")

Contacts: 283, Infections: 200, Traced: 201
Contacts per case: 1.4

Time-dependent policies

In real outbreaks, interventions are not active from the start. Testing may begin on day 14, contact tracing may start once cumulative cases exceed a threshold.

The Scheduled wrapper

Scheduled is the single entry point for time-based scheduling. Wrapping an intervention with Scheduled(...; start_time = ...) filters on action time — an individual is only isolated if their computed isolation time falls after the policy start, regardless of when they were infected. This is a competing risk: the testing infrastructure must be available at the time the individual would be tested.

# Testing starts on day 10
iso_delayed = Scheduled(Isolation(onset_to_isolation_delay = Exponential(2.0)); start_time = 10.0)

rng = StableRNG(42)
results = simulate(scenario([iso_delayed]), 200; max_cases = 500, rng = rng)
println("Isolation from day 10: $(round(containment_probability(results), digits=3))")

Isolation from day 10: 0.055

Someone infected on day 8 with symptom onset on day 9 and delay 2 has isolation time = 11, which is after day 10, so they are isolated. Someone with isolation time = 9 is not isolated — testing was not yet available.

Scheduled also handles conditions that cannot be expressed as a fixed time, such as case-count triggers:

# Start contact tracing after 20 cumulative cases
iso = Isolation(onset_to_isolation_delay = Exponential(2.0))
ct_triggered = Scheduled(
    ContactTracing(probability = 0.7, isolation_to_trace_delay = Exponential(1.0));
    start_after_cases = 20,
)

rng = StableRNG(42)
results = simulate(scenario([iso, ct_triggered]), 200; max_cases = 500, rng = rng)
println("Tracing after 20 cases: $(round(containment_probability(results), digits=3))")

Tracing after 20 cases: 0.195

Conditions can be combined:

# Active only between day 5 and day 30
iso_window = Scheduled(Isolation(onset_to_isolation_delay = Exponential(1.0));
    start_time = 5.0, end_time = 30.0)

Scheduled{Isolation{SymptomaticOnly, Distributions.Exponential{Float64}, Float64}, EpiBranch.var"#38#39"{Vector{Function}}}(Isolation{SymptomaticOnly, Distributions.Exponential{Float64}, Float64}(SymptomaticOnly(), Distributions.Exponential{Float64}(θ=1.0), 1.0, 0.0), EpiBranch.var"#38#39"{Vector{Function}}(Function[EpiBranch.var"#32#33"{Float64}(5.0), EpiBranch.var"#34#35"{Float64}(30.0)]), 5.0)

For full flexibility, pass a predicate on SimulationState:

# Start isolation from generation 3 onwards
iso_gen3 = Scheduled(
    Isolation(onset_to_isolation_delay = Exponential(2.0)),
    state -> state.current_generation >= 3,
)

Scheduled{Isolation{SymptomaticOnly, Distributions.Exponential{Float64}, Float64}, Main.var"#11#12"}(Isolation{SymptomaticOnly, Distributions.Exponential{Float64}, Float64}(SymptomaticOnly(), Distributions.Exponential{Float64}(θ=2.0), 1.0, 0.0), Main.var"#11#12"(), 0.0)

How Scheduled enforces start_time

Scheduled enforces start times at two levels. The population gate (is_active) skips the inner resolve_individual! and apply_post_transmission! until the condition turns true. Once active, the wrapper also performs individual-level reset: after each per-individual hook runs, if intervention_time for that individual falls before start_time, Scheduled calls reset! to undo the effect. Individual interventions therefore do not need to know their own scheduling — they just declare their action time and how to undo it.

Writing a custom intervention

Everything above this point is configuration — combining the interventions that ship with EpiBranch, which is the common case and needs no Julia beyond keyword arguments. This section is extension: writing a new intervention in Julia. That is developer work, alongside adding a transmission model — see Extending EpiBranch for the full picture of the extension points.

Custom interventions are defined as structs subtyping AbstractIntervention. One or more of the following methods should be implemented:

• initialise_individual! — set up fields on new contacts

• resolve_individual! — determine state before transmission

• apply_post_transmission! — act on contacts after creation

To make the intervention schedulable via Scheduled(...; start_time = ...), also implement:

• intervention_time — return the time at which the effect occurs for an individual

• reset! — undo the effect when it falls before start_time

Reference: the built-in Isolation intervention

Before writing your own, it helps to see a complete built-in example. The source code of Isolation is shown below via CodeTracking.jl, so it always reflects the current implementation:

using CodeTracking
print(@code_string EpiBranch.resolve_individual!(iso, state.individuals[1], state))

function resolve_individual!(iso::Isolation, individual, state)
    is_isolated(individual) && return nothing

    # Three isolation pathways, each independent:
    #   - test_isolation_time:  onset + delay, fires iff test_positive
    #   - traced_isolation_time: set by ContactTracing's FlagOnly action
    #     (for symptomatic traced contacts), fires iff contact was traced
    # Isolation fires at the earlier of any active pathway. A
    # test-negative-but-traced contact is still isolated via tracing.
    traced_time = get(individual.state, :traced_isolation_time, Inf)
    test_time = if is_test_positive(individual)
        onset_time(individual) + rand(state.rng, iso.onset_to_isolation_delay)
    else
        Inf
    end
    final = min(test_time, traced_time)
    isfinite(final) || return nothing
    set_isolated!(individual, final)
    # Mark provenance so a Scheduled reset undoes only Isolation's own effect.
    individual.state[:isolated_by_isolation] = true
    return nothing
end

# A gathering limit that caps the number of contacts per individual
struct GatheringLimit <: AbstractIntervention
    max_contacts::Int
end

function EpiBranch.apply_post_transmission!(gl::GatheringLimit, state, new_contacts)
    # Count contacts per parent, mark excess as not infected
    parent_counts = Dict{Int, Int}()
    for c in new_contacts
        count = get(parent_counts, c.parent_id, 0) + 1
        parent_counts[c.parent_id] = count
        if count > gl.max_contacts
            c.state[:infected] = false
        end
    end
end

# Test it
gl = GatheringLimit(5)
model = ModelSpec(BranchingProcess(NegBin(2.5, 0.16), Exponential(5.0)); interventions = [gl])
rng = StableRNG(42)
results_gl = simulate(model, 200; max_cases = 500, rng = rng)
println("With gathering limit (max 5): $(round(containment_probability(results_gl), digits=3))")

With gathering limit (max 5): 0.825