Getting started

The core workflow is covered here: defining a model, running simulations, and extracting results.

Defining a model

At minimum, a branching process is defined by an offspring distribution (how many secondary cases each case produces). If you also want timing (which you need for interventions), add a generation time distribution.

using EpiBranch
using Distributions
using StableRNGs

model = BranchingProcess(
    NegBin(2.5, 0.16),         # R₀ = 2.5, dispersion k = 0.16
    LogNormal(1.6, 0.5)        # generation time
)

BranchingProcess(offspring=Distributions.NegativeBinomial{Float64}, generation_time=Distributions.LogNormal{Float64}, population_size=unlimited)

This is already a complete model. The transmission process describes how infection spreads. The population it acts on (attributes), the policy in force (interventions), and how cases are observed (observation) are separate parts you add with a ModelSpec below, and simulate and loglikelihood then account for all of them.

NegBin is a convenience constructor for the Negative Binomial, parameterised by mean (R) and dispersion (k), matching the epidemiological convention.

Running a simulation

rng = StableRNG(123)
state = simulate(model;
    max_cases = 500,
    rng = rng,
)
println("Cases: $(state.cumulative_cases), Extinct: $(state.extinct)")

Cases: 1, Extinct: true

A SimulationState is returned, containing all individuals and outbreak metadata.

Adding a population

To model symptom onset (needed for isolation-based interventions), provide an incubation period distribution via clinical_presentation as an attributes layer on a ModelSpec:

spec = ModelSpec(BranchingProcess(NegBin(2.5, 0.16), LogNormal(1.6, 0.5));
    attributes = clinical_presentation(incubation_period = LogNormal(1.5, 0.5)),
)

rng = StableRNG(42)
state = simulate(spec; max_cases = 500, rng = rng)

# Check onset times are set
ind = state.individuals[1]
println("Infection: $(round(ind.infection_time, digits=1)), Onset: $(round(onset_time(ind), digits=1))")

Infection: 0.0, Onset: 5.6

Adding interventions

Interventions are another layer you compose with a ModelSpec. simulate reads them from the specification, so there is nothing to pass at the call site:

iso = Isolation(onset_to_isolation_delay = Exponential(2.0))
ct = ContactTracing(probability = 0.5, isolation_to_trace_delay = Exponential(1.5))

spec = ModelSpec(BranchingProcess(NegBin(2.5, 0.16), LogNormal(1.6, 0.5));
    attributes = clinical_presentation(incubation_period = LogNormal(1.5, 0.5)),
    interventions = [iso, ct],
)

rng = StableRNG(42)
state = simulate(spec; max_cases = 500, rng = rng)
println("Cases: $(state.cumulative_cases)")
println("Isolated: $(count(is_isolated, state.individuals))")
println("Traced: $(count(is_traced, state.individuals))")

Cases: 1
Isolated: 1
Traced: 0

Batch simulation

Run many replicates to estimate containment probability:

rng = StableRNG(42)
results = simulate(spec, 500; max_cases = 5000, rng = rng)
println("Containment probability: $(round(containment_probability(results), digits=3))")

Containment probability: 0.96

Contacts vs cases

Every potential transmission is tracked. Contacts that weren't successfully infected are stored alongside cases:

n_total = length(state.individuals)
n_infected = count(is_infected, state.individuals)
println("Total contacts: $n_total")
println("Infected (cases): $n_infected")
println("Not infected: $(n_total - n_infected)")

Total contacts: 1
Infected (cases): 1
Not infected: 0

Intervention effort can be tracked this way. For more information, see Interventions.

Outputs

Simulation state can be converted to DataFrames with several output functions.

using DataFrames, Dates

# Line list (cases only)
ll = linelist(state; reference_date = Date(2024, 1, 1))
println("Line list: $(nrow(ll)) rows, $(ncol(ll)) columns")
first(ll, 3)

1×14 DataFrame

Row	id	parent_id	generation	chain_id	date_infection	date_isolation	date_onset	asymptomatic	incubation_period	isolated	isolated_by_isolation	quarantined	test_positive	traced
	Int64	Int64	Int64	Int64	Date	Date?	Date?	Bool	Float64	Bool	Bool	Bool	Bool	Bool
1	1	0	0	1	2024-01-01	2024-01-08	2024-01-06	false	5.60445	true	true	false	true	false

# Contacts table (all contacts, with infected flag)
ct_df = contacts(state; reference_date = Date(2024, 1, 1))
println("Contacts: $(nrow(ct_df)) ($(count(ct_df.infected)) infected)")

Contacts: 0 (0 infected)

# Chain statistics
cs = chain_statistics(state)
cs

1×3 DataFrame

Row	chain_id	size	length
	Int64	Int64	Int64
1	1	1	0

# Effective R per generation
r_df = generation_R(state)
first(r_df, 5)

0×2 DataFrame

Row	generation	offspring_ratio
	Int64	Float64

Analytical functions

Some quantities have closed-form solutions — no simulation needed:

println("P(extinction): $(round(extinction_probability(spec), digits=3))")
println("P(epidemic):   $(round(epidemic_probability(spec), digits=3))")
println("Top 20% cause $(round(proportion_transmission(spec; prop_cases=0.2) * 100, digits=1))% of transmission")

P(extinction): 0.795
P(epidemic):   0.205
Top 20% cause 88.3% of transmission

Conditioned simulation

Generate outbreaks of a specific size via rejection sampling:

rng = StableRNG(42)
state = simulate(spec;
    condition = 50:100,
    max_cases = 200,
    rng = rng,
)
println("Outbreak size: $(state.cumulative_cases) (target: 50-100)")

Outbreak size: 65 (target: 50-100)

Next steps

• Interventions — combining interventions and the competing risks framework

• Multi-type models — age-structured and heterogeneous transmission

• Line lists and contacts — generating epidemiological data

• Chain statistics and likelihood — inference from chain data

• Analytical functions — closed-form solutions

• Extending EpiBranch — writing a custom intervention or a new transmission model