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 is a pure kernel. The population it acts on (attributes), the policy in force (interventions), and how cases are observed (observation) are modelling layers you compose with a ModelSpec wrapper below; simulate and loglikelihood then read them off the resulting specification.

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.93

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: 52 (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