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.04With 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.18Effectiveness 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
endDelay ~ Exp(0.5): containment = 0.255
Delay ~ Exp(2.0): containment = 0.18
Delay ~ Exp(10.0): containment = 0.12Leaky 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.095Contact 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.27Who 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.155Ring 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.27efficacy adds nothing under this tracing
That number is the same as for isolation and tracing alone, as expected. By default ContactTracing traces a contact once its infector has been isolated, and that isolation already blocks every later transmission to the contact. efficacy protects a contact only against exposure after its immunity arrives, so it has nothing left to prevent, with or without quarantine (quarantine_on_trace = false). It acts only where a contact can still be infected after being traced: under leaky isolation (post_isolation_transmission > 0), when tracing starts before the infector is isolated (for example eligibility = OnSymptomOnset()), or in a depth > 1 ring passing through members who keep transmitting after they are traced. onward_efficacy acts on the traced contact's own later transmission. A quarantine already blocks that transmission, so onward_efficacy acts when tracing does not quarantine. So does post_exposure_efficacy, which acts on the infection the contact already has; see Protecting a contact who has already been exposed. Until that section the examples keep the default tracing and set only efficacy, so they show how the machinery works but no vaccine effect.
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.27We 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: 249, Infected: 200Rings 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: 410A wider ring reaches more people, so it costs more doses, and the extra doses protect only contacts still exposed after they are traced. Here every infected ring member with symptoms is isolated and seeds a ring of its own, so its contacts are traced after that isolation and the second ring adds doses without adding protection. The second ring's doses can prevent exposures when the ring passes through members who keep transmitting after they are traced, such as asymptomatic members under tracing without quarantine. 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 that stops an infection a contact already has is set by post_exposure_efficacy, covered below. Set through efficacy, 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, and under the tracing used here PEP set this way leaves the result unchanged (see the warning above):
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.27Mass 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.065For 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}, Float64, Float64, Nothing, LeakyMode}(0.85, Distributions.Exponential{Float64}(θ=60.0), 14.0, 0.0, nothing, 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", Float64, Float64, Nothing, LeakyMode}(0.85, Main.var"#2#3"(), 14.0, 0.0, nothing, 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.07Per-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", Float64, Float64, Nothing, LeakyMode}(Main.var"#6#7"(), Main.var"#8#9"(), 14.0, 0.0, nothing, 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.045Both 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.
Two doses in a ring
Ring doses are given at the trace, so a second dose sets dose_delay (days from the trace to that dose) and names the dose it follows with requires_dose. Only contacts who have had the earlier dose by the time the later one falls due receive the later one, so the boost's coverage is the retention between doses:
prime_ring = RingVaccination(efficacy = 0.6, delay_to_immunity = 21.0,
coverage = 0.8, dose_label = :prime)
boost_ring = RingVaccination(efficacy = 0.5, dose_delay = 28.0,
delay_to_immunity = 14.0, coverage = 0.9,
requires_dose = :prime, dose_label = :boost)
rng = StableRNG(42)
results = simulate(scenario([iso, ct, prime_ring, boost_ring]), 200;
max_cases = 500, rng = rng)
println("Two-dose ring: $(round(containment_probability(results), digits=3))")Two-dose ring: 0.255A protocol that puts the boost four to six weeks after the trace says so with a distribution: dose_delay = Uniform(28.0, 42.0), drawn once per contact when the dose is scheduled. delay_to_immunity, post_exposure_efficacy and onward_efficacy take a distribution or a function (rng, contact) -> Real in the same way, each drawn once when the dose is given, so a contact meets every exposure with the same immunity time and the same efficacy.
The boost's efficacy is the protection it adds among those the prime left unprotected, because doses compose as competing risks. A schedule described as 60% after one dose and 80% after two therefore needs efficacy = 0.5 on the boost: (0.8 - 0.6) / (1 - 0.6).
Counting doses is cheap, and the count shows what a schedule's second dose costs:
doses(key) = mean(count(i -> i.state[key], s.individuals) for s in results)
println("Primed: $(round(doses(:vaccinated_prime), digits = 1)), ",
"boosted: $(round(doses(:vaccinated_boost), digits = 1))")Primed: 488.3, boosted: 439.4Whether those doses prevent anything depends on timing. Through efficacy, a dose protects a ring member only against infection after its immunity arrives. Under the tracing used here no ring member is infected after being traced (see the warning above), so neither dose changes the outbreak. Where infections do follow the trace, as when tracing starts at the infector's symptom onset and does not quarantine, they follow within days:
ct_onset = ContactTracing(probability = 0.7,
isolation_to_trace_delay = Exponential(1.0),
eligibility = OnSymptomOnset(), quarantine_on_trace = false)
for (label, tracing) in (("after isolation", ct), ("at symptom onset", ct_onset))
let rng = StableRNG(42)
runs = simulate(scenario([iso, tracing]), 200; max_cases = 500, rng = rng)
# Days from the trace to infection, for contacts infected after their trace
days_after = [ind.infection_time - ind.state[:trace_time]
for s in runs for ind in s.individuals
if is_traced(ind) && is_infected(ind) &&
ind.infection_time > get(ind.state, :trace_time, Inf)]
println("Tracing $label: $(length(days_after)) infected after their trace, ",
"$(count(>(7.0), days_after)) of them more than a week after")
end
endTracing after isolation: 0 infected after their trace, 0 of them more than a week after
Tracing at symptom onset: 4320 infected after their trace, 14 of them more than a week afterThose infections come from infectors who have not yet isolated, which they do within days of onset. A dose on the day of the trace can still prevent them. Protection arriving a week later comes too late for almost all of them, and a boost given 28 days after the trace comes later still. Across many simulated outbreaks under this tracing, a trace-day dose raises containment and the same efficacy a week later does not.
Adding an intervention is not a controlled comparison
Re-running with the boost removed will not show this, and may suggest the opposite. boost_ring has coverage = 0.9, so it draws from the rng for every primed contact and shifts the whole stream. The two runs are then different samples, not the same outbreak with and without a dose. Here containment is 0.255 with the boost and 0.295 without, and that difference is Monte Carlo noise, not harm from the boost. A dose with coverage = 1.0, scalar efficacy and no eligibility window draws from the rng only when it meets an exposure it could block. Such a dose with nothing to block leaves the stream untouched, which is why the RingVaccination(efficacy = 0.8) case above reproduces its baseline exactly. For any other dose, compare across many seeds, or reason from the timing as here.
Protecting a contact who has already been exposed
Under default tracing a contact is traced only after its infector has been isolated, so by the time of the trace it has already been exposed and can no longer be infected. efficacy, which protects only against exposure after immunity arrives, then has nothing to block. The dose can still act on the infection the contact already has, in two ways:
post_exposure_efficacyaborts the infection, with that probability, if immunity arrives before the contact's symptom onset. The contact transmits as usual until immunity arrives and not at all afterwards. It has no symptom onset, and its clinical course ends when immunity arrives: nothing in theprogression(hospitalisation, death, recovery) happens from then on. It still counts as a case.onward_efficacyleaves the infection and its disease alone and blocks each of the contact's transmissions after immunity with that probability, whenever its onset falls.
The two compose: a dose setting both aborts some infections and makes the rest less infectious. Both act only on transmission after immunity arrives, so neither does much for a contact that has infected most of its own contacts by then. A quarantine from the trace already blocks that later transmission, so the examples here trace without quarantine. Infectiousness here starts at infection, and the clearest measure of the effect is the number of people each traced case goes on to infect:
ct_flag = ContactTracing(probability = 0.7,
isolation_to_trace_delay = Exponential(1.0), quarantine_on_trace = false)
# Infections per traced case, over cases whose own contacts were simulated
function onward_per_traced(runs)
mean(count(id -> is_infected(s.individuals[id]), ind.secondary_case_ids)
for s in runs for ind in s.individuals
if is_infected(ind) && is_traced(ind) && ind.generation < s.current_generation)
end
for (label, rv) in [
("no vaccine", nothing),
("efficacy = 0.9", RingVaccination(efficacy = 0.9)),
("onward_efficacy = 0.9",
RingVaccination(efficacy = 0.0, onward_efficacy = 0.9)),
("post_exposure_efficacy = 0.9",
RingVaccination(efficacy = 0.0, post_exposure_efficacy = 0.9)),
]
stack = rv === nothing ? [iso, ct_flag] : [iso, ct_flag, rv]
rng = StableRNG(42)
results = simulate(scenario(stack), 400; max_cases = 500, rng = rng)
println(rpad(label, 30), "infections per traced case ",
round(onward_per_traced(results), digits = 2),
", containment ", round(containment_probability(results), digits = 3))
end┌ Warning: Assignment to `rng` in soft scope is ambiguous because a global variable by the same name exists: `rng` will be treated as a new local. Disambiguate by using `local rng` to suppress this warning or `global rng` to assign to the existing global variable.
└ @ interventions.md:504
┌ Warning: Assignment to `results` in soft scope is ambiguous because a global variable by the same name exists: `results` will be treated as a new local. Disambiguate by using `local results` to suppress this warning or `global results` to assign to the existing global variable.
└ @ interventions.md:505
no vaccine infections per traced case 1.93, containment 0.215
efficacy = 0.9 infections per traced case 1.93, containment 0.215
onward_efficacy = 0.9 infections per traced case 1.66, containment 0.242
post_exposure_efficacy = 0.9 infections per traced case 1.65, containment 0.225Either parameter cuts the infections each traced case causes by about a seventh, while efficacy changes nothing. Containment moves much less, because each traced contact has usually made a good share of its transmissions before the trace. Across many seeds both raise it from about 0.21 to about 0.24; at 400 replicates that difference is barely larger than the Monte Carlo noise, so a single run of this block can place the two in either order.
How much post_exposure_efficacy achieves depends on speed. Immunity has to arrive before onset, and incubation periods here average about five days:
for days in [0.0, 2.0, 5.0, 21.0]
let rv = RingVaccination(efficacy = 0.0, post_exposure_efficacy = 0.9,
delay_to_immunity = days),
rng = StableRNG(42)
results = simulate(scenario([iso, ct_flag, rv]), 400;
max_cases = 500, rng = rng)
cases = [ind for s in results for ind in s.individuals
if is_infected(ind) && is_traced(ind)]
aborted = count(ind -> haskey(ind.state, :infection_aborted_time), cases)
println("Immunity after $(lpad(Int(days), 2)) days: ",
round(Int, 100 * aborted / length(cases)), "% of traced cases aborted, ",
"infections per traced case ",
round(onward_per_traced(results), digits = 2))
end
endImmunity after 0 days: 43% of traced cases aborted, infections per traced case 1.65
Immunity after 2 days: 26% of traced cases aborted, infections per traced case 1.83
Immunity after 5 days: 8% of traced cases aborted, infections per traced case 1.92
Immunity after 21 days: 0% of traced cases aborted, infections per traced case 1.95With same-day immunity the dose aborts over two in five traced cases. A two-day delay brings that to about a quarter and removes most of the reduction in onward infections. With a five-day delay little is left, and with immunity three weeks after the trace nothing is. An aborted case appears in the line list with a date_infection_aborted and no onset date.
Where immunity is in place before the exposure, post_exposure_efficacy blocks the infection itself. That is all it can do for a contact with no onset to race (an asymptomatic one, whose incubation period is NaN). It also means post_exposure_efficacy covers every contact efficacy would protect, so setting both double-counts. post_exposure_efficacy requires :incubation_period, set by clinical_presentation.
Under quarantine neither parameter has transmission left to block, and post_exposure_efficacy can lower containment: an aborted case never has symptoms, so the contacts it infected before its dose are no longer traced from it.
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: 360, Infections: 200, Traced: 249
Contacts per case: 1.8Time-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.055Someone 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.125Conditions 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"#44#45"{Vector{Function}}}(Isolation{SymptomaticOnly, Distributions.Exponential{Float64}, Float64}(SymptomaticOnly(), Distributions.Exponential{Float64}(θ=1.0), 1.0, 0.0), EpiBranch.var"#44#45"{Vector{Function}}(Function[EpiBranch.var"#38#39"{Float64}(5.0), EpiBranch.var"#40#41"{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"#35#36"}(Isolation{SymptomaticOnly, Distributions.Exponential{Float64}, Float64}(SymptomaticOnly(), Distributions.Exponential{Float64}(θ=2.0), 1.0, 0.0), Main.var"#35#36"(), 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 contactsresolve_individual!— determine state before transmissionapply_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 individualreset!— undo the effect when it falls beforestart_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)
# An isolation already standing on the individual is a quarantine written by
# `ContactTracing`. On the generation-based path that cannot happen here,
# because isolation resolves before tracing runs; on the continuous-time
# path it routinely does, because a case's contacts are traced when the
# *infector* is finalised, which is before the contact itself resolves. The
# quarantine is then a competing pathway rather than a reason to stop:
# without this the contact would keep a trace time later than the onset it
# would have self-reported on, and tracing would delay isolation instead of
# advancing it.
if is_isolated(individual)
is_test_positive(individual) || return nothing
self_t = onset_time(individual) + rand(state.rng, iso.onset_to_isolation_delay)
self_t < isolation_time(individual) || return nothing
# Remember what we are overwriting. Claiming provenance below tells a
# `Scheduled` reset that this isolation is Isolation's to undo, but the
# standing quarantine underneath it belongs to ContactTracing and must
# survive that reset, so stash it for `reset!` to restore.
individual.state[:isolation_time_before_isolation] = isolation_time(individual)
set_isolated!(individual, self_t)
individual.state[:isolated_by_isolation] = true
return nothing
end
# Three isolation pathways, each independent:
# - test_isolation_time: onset + delay, fires iff test_positive
# - traced_isolation_time: set by ContactTracing's FlagOnly action
# for traced contacts, fires iff contact was traced and has an onset
# Isolation fires at the earlier of any active pathway. A
# test-negative-but-traced contact is still isolated via tracing.
#
# The traced pathway isolates a flagged contact once it has symptoms, so a
# contact with no onset never isolates through it. That includes a contact
# whose infection was aborted before onset, after the trace had already
# recorded its expected onset. A continuous-time model traces a contact
# before its own infection is settled, when the onset is not yet known, so
# the recorded time is held back to the onset.
onset = onset_time(individual)
traced_time = isnan(onset) ? Inf :
max(get(individual.state, :traced_isolation_time, Inf), onset)
test_time = if is_test_positive(individual)
onset + 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