The previous chapter built up three layers of uncertainty — parameter, process, and observation — and showed how to interrogate each through simulation. But understanding the baseline epidemic is not the end goal. The questions that matter for policy are counterfactual: what would change if we intervened? And the answers must be robust: do they hold across plausible parameter values?
This chapter introduces camdl’s tools for asking and answering those questions systematically: interventions that modify the model mid-simulation, scenarios that package counterfactual comparisons, parameter sweeps that test robustness, and a batch system that manages hundreds of runs reproducibly.
camdl’s interventions block declares actions that change the population structure at specific times. To model a vaccination campaign, we first add a Vaccinated compartment and then declare the intervention:
compartments { S, I, R, V }
interventions {
vaccination : transfer(fraction = vacc_eff, from = S, to = V) at [15]
}This reads: “on day 15, move a fraction vacc_eff of the susceptible population into the vaccinated compartment.” The transfer action takes individuals from one compartment and puts them in another — here, vaccinating 80% of remaining susceptibles in a single pulse.
Interventions are off by default. The model file declares what’s possible, not what happens in every run. To activate the vaccination, pass --enable:
$ camdl simulate sirv.camdl --params params.toml --enable vaccination --seed 42# camdl 0.1.0+24c1977 (2026-04-27)
t S I R V flow_infection flow_infection_v flow_recovery
0 9990 10 0 0 0 0 0
1 9986 9 5 0 4 0 5
2 9982 10 8 0 4 0 3
3 9977 13 10 0 5 0 2
4 9968 20 12 0 9 0 2
5 9963 21 16 0 5 0 4no_vacc = run_simulate(MODEL, PARAMS, seed=42, cwd=EXAMPLE_DIR)
with_vacc = run_simulate(MODEL, PARAMS, seed=42, cwd=EXAMPLE_DIR,
extra_args=["--enable", "vaccination"])
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 3.5), sharey=True)
for ax, df, title in [(ax1, no_vacc, "Intervention off (default)"),
(ax2, with_vacc, "--enable vaccination")]:
colors = {"S": "#3366cc", "I": "#cc3333", "R": "#2e8b2e", "V": "#9b59b6"}
for col, color in colors.items():
if col in df.columns:
ax.plot(df["t"], df[col], color=color, label=col, linewidth=1.5)
ax.set_xlabel("Time (days)")
ax.set_title(title, fontsize=11, color=GREY)
ax.legend(frameon=False, fontsize=9)
ax.set_xlim(0, df["t"].max())
ax.spines[["top", "right"]].set_visible(False)
ax1.set_ylabel("Individuals")
plt.tight_layout()
plt.show()
--enable vaccination. Left: baseline run with the intervention declared but inactive — the epidemic runs unchecked. Right: enabling the vaccination intervention moves 80% of susceptibles to V on day 15, sharply reducing transmission.
camdl supports several scheduling patterns:
interventions {
# Single pulse
sia_round_1 : transfer(fraction = 0.8, from = S, to = V) at [180]
# Multiple pulses
sia_rounds : transfer(fraction = 0.8, from = S, to = V) at [180, 545, 910]
# Recurring schedule
routine_vacc : transfer(fraction = vacc_rate, from = S, to = V) {
every = 30 'days
from = 0 'days
until = 2 'years
}
}Interventions are off by default because they represent hypothetical policy levers — “what if we vaccinated?” But some scheduled actions aren’t hypothetical. If you’re fitting a model to data from a region that had three historical vaccination campaigns, those campaigns happened — they shaped the data you’re fitting to, and omitting them is wrong.
camdl has a separate events {} block for this. Events use the same syntax as interventions but are always active — they fire in every simulation and during inference without any scenario or --enable flag:
tables {
sia_dates : round = read("data/campaign_dates.tsv")
}
events {
historical_sia[r in round] :
transfer(fraction = sia_coverage, from = S, to = V)
at [sia_dates[r]]
}The schedule comes from a data file — one row per campaign with its date. No hardcoded times, no scenario boilerplate. When you run camdl fit run fit.toml, the events fire automatically at the right times.
The decision rule is simple:
events {}. Always on, zero config in fit.toml.interventions {}. Off by default, enable via scenario.If you later want a counterfactual — “what if the third campaign hadn’t happened?” — you can disable individual events with --disable:
camdl simulate model.camdl --scenario baseline --disable historical_sia_r3We’ll see events {} used for annual cohort entry in the He et al. measles vignette, where each year’s birth cohort enters the susceptible class at school start — a structural demographic process, not a policy choice.
A scenario packages a set of parameter values and intervention switches into a named configuration. This is where the counterfactual becomes explicit.
The --enable flag we used above is fine for quick one-off runs, but for systematic comparisons we want the configurations named and version-controlled. The SIRV model declares three scenarios:
scenarios {
baseline {
label = "No intervention"
set = {
beta = 0.4
gamma = 0.15
N0 = 10_000
I0 = 10
rho = 0.6
k = 10.0
vacc_eff = 0.8
leaky_susc = 0.0
}
}
with_vaccination {
extends = baseline
label = "All-or-nothing vaccination at day 15"
enable = [vaccination]
}
with_leaky_vaccination {
extends = baseline
label = "Leaky vaccination at day 15"
enable = [vaccination]
set = { leaky_susc = 0.3 }
}
}The extends = baseline keyword says “inherit everything from baseline, then apply my overrides.” The with_vaccination scenario only adds enable = [vaccination] — same parameters, same initial conditions, so the comparison isolates the causal effect.
The with_leaky_vaccination scenario goes further: it inherits from baseline, enables the same intervention, and overrides one parameter. With leaky_susc = 0.0 (the default), vaccinated individuals in V are completely protected — the V --> I transition has rate zero. With leaky_susc = 0.3, every vaccinee retains 30% residual susceptibility. This is the distinction between an all-or-nothing vaccine (a fraction are fully protected, the rest get nothing) and a leaky vaccine (everyone gets partial protection). Same coverage, same VE on average, very different epidemic dynamics — a structural question answered by a single parameter.
Run all three with:
$ camdl simulate sirv.camdl --scenario baseline --seed 42# camdl 0.1.0+24c1977 (2026-04-27)
t S I R V flow_infection flow_infection_v flow_recovery
0 9990 10 0 0 0 0 0
1 9986 9 5 0 4 0 5
2 9982 10 8 0 4 0 3
3 9977 13 10 0 5 0 2
4 9968 20 12 0 9 0 2
5 9963 21 16 0 5 0 4$ camdl simulate sirv.camdl --scenario with_vaccination --seed 42# camdl 0.1.0+24c1977 (2026-04-27)
t S I R V flow_infection flow_infection_v flow_recovery
0 9990 10 0 0 0 0 0
1 9986 9 5 0 4 0 5
2 9982 10 8 0 4 0 3
3 9977 13 10 0 5 0 2
4 9968 20 12 0 9 0 2
5 9963 21 16 0 5 0 4$ camdl simulate sirv.camdl --scenario with_leaky_vaccination --seed 42# camdl 0.1.0+24c1977 (2026-04-27)
t S I R V flow_infection flow_infection_v flow_recovery
0 9990 10 0 0 0 0 0
1 9986 9 5 0 4 0 5
2 9982 10 8 0 4 0 3
3 9977 13 10 0 5 0 2
4 9968 20 12 0 9 0 2
5 9963 21 16 0 5 0 4Our usage of baseline as preset parameters might seem to contract camdls “no default parameter values” rule introduced in Chapter 2. However, the distinction is deliberate: the parameters block defines the model’s degrees of freedom — the knobs that may be the focus of estimation or inference from observed data. When you run a simulation without specifying all paramters or a particular scenario like baseline, camdl errors out and asks you to supply parameters explicitly.
Scenarios are a separate mechanism, walled off from inference. They exist so researchers can package named configurations for reproducible simulation, without those presets leaking into the fitting pipeline as implicit defaults. This separation is the same principle as the “no default parameters” rule — every assumption must be visible and explicit.
Mechanically: --scenario baseline supplies all parameter values from the scenario’s set block. If you also pass --params, the scenario values take precedence — the scenario is the final word.
$ camdl simulate sirv.camdl --params params.toml --seed 42 --backend chain_binomial --dt 1.0 --scenario baseline(no output)$ camdl simulate sirv.camdl --params params.toml --seed 42 --backend chain_binomial --dt 1.0 --scenario with_vaccination(no output)$ camdl simulate sirv.camdl --params params.toml --seed 42 --backend chain_binomial --dt 1.0 --scenario with_leaky_vaccination(no output)fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12, 3.5), sharey=True)
for ax, df, title in [(ax1, baseline, "Baseline"),
(ax2, vaccinated, "All-or-nothing"),
(ax3, leaky, "Leaky (30% residual)")]:
colors = {"S": "#3366cc", "I": "#cc3333", "R": "#2e8b2e", "V": "#9b59b6"}
for col, color in colors.items():
if col in df.columns:
ax.plot(df["t"], df[col], color=color, label=col, linewidth=1.5)
ax.set_xlabel("Time (days)")
ax.set_title(title, fontsize=11, color=GREY)
ax.legend(frameon=False, fontsize=8)
ax.set_xlim(0, df["t"].max())
ax.spines[["top", "right"]].set_visible(False)
ax1.set_ylabel("Individuals")
plt.tight_layout()
plt.show()
leaky_susc = 0.3). The leaky vaccine produces a larger, delayed epidemic because every vaccinee can still be infected.
But a single trajectory per scenario is exactly the mistake we warned about. To compare scenarios honestly, we need the distribution of outcomes under each:
$ camdl simulate sirv.camdl --params params.toml --seed 42 --backend chain_binomial --dt 1.0 --replicates 200 --scenario baseline200 replicates = 200 runs$ camdl simulate sirv.camdl --params params.toml --seed 42 --backend chain_binomial --dt 1.0 --replicates 200 --scenario with_vaccination200 replicates = 200 runs$ camdl simulate sirv.camdl --params params.toml --seed 42 --backend chain_binomial --dt 1.0 --replicates 200 --scenario with_leaky_vaccination200 replicates = 200 runsfig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12, 3.5), sharey=True)
ribbon(ax1, baseline_multi, "I", "#cc3333", "Infected")
ax1.set_title("Baseline", fontsize=11, color=GREY)
ax1.set_xlabel("Time (days)")
ax1.set_ylabel("Infected individuals")
ribbon(ax2, vacc_multi, "I", "#cc3333", "Infected")
ax2.axvline(15, color="#9b59b6", linestyle="--", alpha=0.5, label="Vaccination")
ax2.set_title("All-or-nothing", fontsize=11, color=GREY)
ax2.set_xlabel("Time (days)")
ribbon(ax3, leaky_multi, "I", "#cc3333", "Infected")
ax3.axvline(15, color="#9b59b6", linestyle="--", alpha=0.5, label="Vaccination")
ax3.set_title("Leaky (30% residual)", fontsize=11, color=GREY)
ax3.set_xlabel("Time (days)")
for ax in (ax1, ax2, ax3):
ax.legend(frameon=False, fontsize=8)
ax.set_xlim(0, 120)
ax.spines[["top", "right"]].set_visible(False)
plt.tight_layout()
plt.show()
def get_peaks(multi_df):
return (
multi_df.group_by("replicate")
.agg(pl.col("I").max().alias("peak_I"))
.sort("replicate")
)
peaks_base = get_peaks(baseline_multi)["peak_I"].to_numpy()
peaks_vacc = get_peaks(vacc_multi)["peak_I"].to_numpy()
peaks_leaky = get_peaks(leaky_multi)["peak_I"].to_numpy()
fig, ax = plt.subplots(figsize=(7, 3.5))
all_peaks = np.concatenate([peaks_base, peaks_vacc, peaks_leaky])
bins = np.linspace(0, all_peaks.max() * 1.05, 40)
ax.hist(peaks_base, bins=bins, alpha=0.6, color="#cc3333",
label="Baseline", edgecolor="none")
ax.hist(peaks_leaky, bins=bins, alpha=0.6, color="#e67e22",
label="Leaky vaccine", edgecolor="none")
ax.hist(peaks_vacc, bins=bins, alpha=0.6, color="#9b59b6",
label="All-or-nothing", edgecolor="none")
ax.set_xlabel("Peak infected count")
ax.set_ylabel("Frequency")
ax.legend(frameon=False, fontsize=9)
ax.spines[["top", "right"]].set_visible(False)
plt.tight_layout()
plt.show()
The ribbons and histograms tell the story that single trajectories can’t. The all-or-nothing vaccine doesn’t just reduce the peak — it narrows the range of possible outcomes. The leaky vaccine reduces the median peak but the distribution stays wide: because every vaccinee retains some susceptibility, stochastic variation in when those breakthrough infections happen produces a broader range of epidemic trajectories. The mechanism of protection matters, not just the coverage — and that’s a structural question about the model, not a parameter you can sweep over.
The results above assumed \(\beta = 0.4\). But what if transmission is higher or lower? Does vaccination still help at \(\beta = 0.6\)? Does it matter at \(\beta = 0.2\)?
A parameter sweep runs the model across a grid of values, crossed with scenarios and replicates. This is where individual CLI calls become unwieldy — you’d need a loop over \(\beta\) values, each with --param beta=X. camdl’s batch system handles this natively with a TOML manifest:
# sweep.toml
[config]
model = "sirv.camdl"
params = "params.toml"
seeds = { n = 50 }
[[scenario]]
name = "baseline"
[[scenario]]
name = "with_vaccination"
enable = ["vaccination"]
[sweep]
beta = { linspace = { min = 0.2, max = 0.6, n = 5 } }$ camdl batch run sweep.toml --dry-runcamdl batch run (dry run)
model: sirv.camdl
backend: chain_binomial
dt: 1
output_dir: results
parallel: 1
Scenarios (2):
baseline (baseline)
with_vaccination enable=vaccination
Sweep grid (5 points × 2 scenarios × 50 seeds = 500 runs):
point 0:
I0 = 10 params file: params.toml
N0 = 10000 params file: params.toml
beta = 0.2 sweep override
gamma = 0.15 params file: params.toml
k = 10 params file: params.toml
leaky_susc = 0 params file: params.toml
rho = 0.6 params file: params.toml
vacc_eff = 0.8 params file: params.toml
point 1:
beta = 0.3 sweep override
point 2:
beta = 0.4 sweep override
point 3:
beta = 0.5 sweep override
point 4:
[...1 lines omitted...]
Cache status:
0 cache hits → skipped
500 cache misses → would simulate
(dry run — no simulation, no files written.)The batch system manages the Cartesian product: 5 sweep points × 2 scenarios × 50 seeds = 500 runs, parallelized across cores and written to content-addressed output directories.
# Run the sweep via batch subcommand (writes to content-addressed output)
batch_dir = tempfile.mkdtemp()
for f in ["sirv.camdl", "params.toml", "sweep.toml"]:
shutil.copy(os.path.join(EXAMPLE_DIR, f), batch_dir)
display(run_cli(
"camdl batch run sweep.toml",
echo="camdl batch run sweep.toml",
cwd=batch_dir, check=True,
))
# Read manifest and all trajectories
with open(os.path.join(batch_dir, "results", "sims", "manifest.json")) as f:
manifest = json.load(f)
sweep_results = []
for run in manifest["runs"]:
traj_path = os.path.join(batch_dir, "results", "sims", run["run_path"], "traj.tsv")
traj = pl.read_csv(traj_path, separator="\t")
peak = traj["I"].max()
sweep_results.append({
"beta": run["sweep_point"]["beta"],
"scenario": run["scenario"],
"peak_I": peak,
})
sweep_df = pl.DataFrame(sweep_results)
betas = sorted(sweep_df["beta"].unique().to_list())
fig, ax = plt.subplots(figsize=(7, 4))
colors = {"baseline": "#cc3333", "with_vaccination": "#9b59b6"}
labels = {"baseline": "Baseline", "with_vaccination": "With vaccination"}
for scenario, color in colors.items():
sub = sweep_df.filter(pl.col("scenario") == scenario)
for beta in betas:
pts = sub.filter(pl.col("beta") == beta)["peak_I"].to_numpy()
ax.scatter([beta] * len(pts), pts, color=color, alpha=0.15, s=10,
edgecolors="none")
medians = (
sub.group_by("beta")
.agg(pl.col("peak_I").median())
.sort("beta")
)
ax.plot(medians["beta"], medians["peak_I"], color=color,
linewidth=2, marker="o", markersize=5, label=labels[scenario])
ax.set_xlabel("Transmission rate (β)")
ax.set_ylabel("Peak infected count")
ax.legend(frameon=False, fontsize=9)
ax.spines[["top", "right"]].set_visible(False)
plt.tight_layout()
plt.show()$ camdl batch run sweep.tomlSweep: 5 parameter points
point 0: beta=0.2
point 1: beta=0.3
point 2: beta=0.4
... (2 more)
[2026-04-28T05:04:18Z INFO camdl::util] --param beta=0.2 overrides previous value 0.4
[1/500] scenario=baseline seed=1
[2026-04-28T05:04:18Z INFO camdl::util] --param beta=0.2 overrides previous value 0.4
[2/500] scenario=baseline seed=2
[2026-04-28T05:04:18Z INFO camdl::util] --param beta=0.2 overrides previous value 0.4
[3/500] scenario=baseline seed=3
[2026-04-28T05:04:18Z INFO camdl::util] --param beta=0.2 overrides previous value 0.4
[4/500] scenario=baseline seed=4
[2026-04-28T05:04:18Z INFO camdl::util] --param beta=0.2 overrides previous value 0.4
[5/500] scenario=baseline seed=5
[2026-04-28T05:04:18Z INFO camdl::util] --param beta=0.2 overrides previous value 0.4
[6/500] scenario=baseline seed=6
[2026-04-28T05:04:18Z INFO camdl::util] --param beta=0.2 overrides previous value 0.4
[7/500] scenario=baseline seed=7
[2026-04-28T05:04:18Z INFO camdl::util] --param beta=0.2 overrides previous value 0.4
[8/500] scenario=baseline seed=8
[2026-04-28T05:04:18Z INFO camdl::util] --param beta=0.2 overrides previous value 0.4
[9/500] scenario=baseline seed=9
[2026-04-28T05:04:18Z INFO camdl::util] --param beta=0.2 overrides previous value 0.4
[10/500] scenario=baseline seed=10
[2026-04-28T05:04:18Z INFO camdl::util] --param beta=0.2 overrides previous value 0.4
[11/500] scenario=baseline seed=11
[2026-04-28T05:04:18Z INFO camdl::util] --param beta=0.2 overrides previous value 0.4
[12/500] scenario=baseline seed=12
[2026-04-28T05:04:18Z INFO camdl::util] --param beta=0.2 overrides previous value 0.4
[13/500] scenario=baseline seed=13
[2026-04-28T05:04:18Z INFO camdl::util] --param beta=0.2 overrides previous value 0.4
[14/500] scenario=baseline seed=14
[2026-04-28T05:04:18Z INFO camdl::util] --param beta=0.2 overrides previous value 0.4
[...866 lines omitted...]
[498/500] scenario=with_vaccination seed=48
[2026-04-28T05:04:19Z INFO camdl::util] --param beta=0.6 overrides previous value 0.4
[499/500] scenario=with_vaccination seed=49
[2026-04-28T05:04:19Z INFO camdl::util] --param beta=0.6 overrides previous value 0.4
[500/500] scenario=with_vaccination seed=50
Done: 500/500 runs completed. Manifest: results/sims/manifest.json
After a batch completes (or during a long run), camdl batch status shows what’s done:
$ camdl batch status sweep.tomlExperiment status for: sweep.toml
Model: sirv.camdl
Output dir: results
Scenarios: baseline, with_vaccination
Seeds: 50 total (1..=50)
Completed: 500/500
thread 'main' (21552181) panicked at crates/cli/src/hashing.rs:11:10:
model_hash: invalid JSON: Error("expected ident", line: 1, column: 2)
note: run with `RUST_BACKTRACE=1` environment variable to display a backtraceThe sweep reveals that vaccination consistently reduces the peak across all plausible transmission rates. The conclusion is robust — it doesn’t depend on getting \(\beta\) exactly right. If the conclusion had flipped at some plausible \(\beta\), that would mean we need better data to estimate \(\beta\) before the policy comparison is meaningful.
This is sensitivity analysis applied to simulation: not “what is the answer?” but “does the answer change when I vary my assumptions?”
The batch run above generated 500 simulations. Where did they go? camdl list browses the output store:
$ camdl list --scenario baseline | head -12 CREATED MODEL SCENARIO SEED PARAMS SIZE PATH
just now camdl_32677 baseline 33 beta=0.5 2K ./results/sims/camdl_32677-2debeb21/baseline-b9896440/seed_33
just now camdl_32677 baseline 34 beta=0.5 2K ./results/sims/camdl_32677-2debeb21/baseline-b9896440/seed_34
just now camdl_32677 baseline 35 beta=0.5 2K ./results/sims/camdl_32677-2debeb21/baseline-b9896440/seed_35
just now camdl_32677 baseline 32 beta=0.5 2K ./results/sims/camdl_32677-2debeb21/baseline-b9896440/seed_32
just now camdl_32677 baseline 50 beta=0.5 2K ./results/sims/camdl_32677-2debeb21/baseline-b9896440/seed_50
just now camdl_32677 baseline 43 beta=0.5 2K ./results/sims/camdl_32677-2debeb21/baseline-b9896440/seed_43
just now camdl_32677 baseline 44 beta=0.5 2K ./results/sims/camdl_32677-2debeb21/baseline-b9896440/seed_44
just now camdl_32677 baseline 21 beta=0.5 2K ./results/sims/camdl_32677-2debeb21/baseline-b9896440/seed_21
just now camdl_32677 baseline 19 beta=0.5 2K ./results/sims/camdl_32677-2debeb21/baseline-b9896440/seed_19
just now camdl_32677 baseline 26 beta=0.5 2K ./results/sims/camdl_32677-2debeb21/baseline-b9896440/seed_26
just now camdl_32677 baseline 10 beta=0.5 2K ./results/sims/camdl_32677-2debeb21/baseline-b9896440/seed_10Each row is one run. The PARAMS column shows the sweep point (beta=0.2, beta=0.3, …) — every run knows what produced it. The PATH column reveals the layout: runs live under content-addressed hashes derived from their exact inputs (model source, parameter values, backend configuration). This is the same principle behind Git’s object store and Nix’s /nix/store. The hash is the identity.
camdl show displays the full provenance for a single run:
$ camdl show results/sims/camdl_32677-2debeb21/baseline-e69e1aa7/seed_1path
results/sims/camdl_32677-2debeb21/baseline-e69e1aa7/seed_1
model
/var/folders/mp/j_ybkq8j72ndvhfm_1msdj1c0000gn/T/camdl_32677.ir.json
scenario
baseline
seed
1
backend
chain_binomial (dt = 1)
hashes
sim 2debeb21aae747214247b666727dc2217b9db1a19824b92521a8ea476e8de4bb
scen e69e1aa71701d5b6efea94543b2e0bc56bbc4fc135c576b1e102f5f4ee422076
model f50ef45413d00a64da2a2e9c38dd4f1cd532f2413c8a816d0a72b41183a2fa47
created
2026-04-28T05:04:18Z (just now)
version
0.1.0+24c1977
argv
camdl batch run sweep.toml
trajectory
2114 bytesAnd camdl cat emits the cached trajectory — useful for piping into analysis tools without hunting for file paths:
$ camdl cat results/sims/camdl_32677-2debeb21/baseline-e69e1aa7/seed_1 | head -6t S I R V
0 9990 10 0 0
1 9989 9 2 0
2 9986 12 2 0
3 9986 10 4 0
4 9984 12 4 0The hash-based layout has three important properties:
Every output is traceable to its exact inputs. When a collaborator asks “which parameters produced Figure 3?” the answer is in camdl show. When a reviewer asks “can you rerun with updated data?” you change the data, and the hash changes — old results are preserved, new results land in a new directory. No results_v3_final_FINAL2/.
The scenarios above used a toy SIRV to introduce the syntax. Now let’s apply the same machinery to the boarding school outbreak with realistically parameterized interventions. The premise: what if we could apply modern public health tools to the 1978 epidemic?
In January 1978, the available interventions were limited. Amantadine existed but was not routinely deployed in UK schools. Influenza vaccines existed but were not matched to the A/USSR/77 H1N1 strain that had re-emerged after a 20-year absence. The realistic response was non-pharmaceutical: isolating symptomatic boys in the infirmary and restricting dormitory mixing.
We’ll model two interventions — one reactive, one proactive — and then combine them:
Reactive isolation: starting on day 4 (when the school recognizes the outbreak), 50% of currently infectious boys are moved daily to an isolation ward (Q compartment) where they recover but no longer transmit. This reflects realistic detection delays and imperfect compliance (Wu et al., 2006).
Hypothetical vaccination: a pre-epidemic pulse that moves a fraction of susceptibles to a vaccinated compartment V. We use a vaccine effectiveness of 70% (Jefferson et al., 2005) at two coverage levels (40% and 80%), giving effective immune fractions of 28% and 56%.
Combined: 60% coverage vaccination (42% effectively immune) plus reactive isolation.
compartments { S, I, R, Q, V }
interventions {
vaccination : transfer(fraction = vacc_frac, from = S, to = V) at [1]
isolation : transfer(fraction = isol_frac, from = I, to = Q) {
every = 1 'days
from = 4 'days
until = 14 'days
}
}The Q compartment removes individuals from the force of infection (\(\beta S I / N\) depends on \(I\), not \(I + Q\)). Isolated boys still recover at rate \(\gamma\), but they can’t infect anyone.
$ camdl simulate boarding_school_intv.camdl --params boarding_school_params.toml --seed 7 --backend chain_binomial --dt 1.0 --replicates 200 --scenario baseline200 replicates = 200 runsfig, axes = plt.subplots(2, 2, figsize=(10, 7), sharey=True)
configs = [
("isolation_only", "Isolation only (day 4)", "#e67e22", axes[0, 0]),
("vaccination_low", "Vaccination (40% cov × 70% VE)", "#27ae60", axes[0, 1]),
("vaccination_high", "Vaccination (80% cov × 70% VE)", "#2980b9", axes[1, 0]),
("combined", "Vaccination (60%) + isolation", "#8e44ad", axes[1, 1]),
]
for sc, title, color, ax in configs:
ribbon(ax, intv_scenarios["baseline"], "I", "#cc3333", "Baseline")
ribbon(ax, intv_scenarios[sc], "I", color, sc.replace("_", " ").title())
ax.scatter(flu_data["day"], flu_data["in_bed"], color="black",
s=15, zorder=5, alpha=0.6)
ax.set_title(title, fontsize=10, color=GREY)
ax.set_xlabel("Day")
ax.set_xlim(-0.5, 14.5)
ax.spines[["top", "right"]].set_visible(False)
ax.legend(frameon=False, fontsize=7, loc="upper right")
axes[0, 0].set_ylabel("Infected individuals")
axes[1, 0].set_ylabel("Infected individuals")
plt.tight_layout()
plt.show()
def get_peaks_intv(multi_df, active_col="I"):
return (
multi_df.group_by("replicate")
.agg(pl.col(active_col).max().alias("peak"))
.sort("replicate")
)
fig, ax = plt.subplots(figsize=(8, 4))
colors = {
"baseline": ("#cc3333", "Baseline"),
"isolation_only": ("#e67e22", "Isolation only"),
"vaccination_high": ("#2980b9", "Vacc (80% cov)"),
"combined": ("#8e44ad", "Combined"),
}
all_peaks = []
for sc, (color, label) in colors.items():
peaks = get_peaks_intv(intv_scenarios[sc])["peak"].to_numpy()
all_peaks.extend(peaks)
bins = np.linspace(0, max(all_peaks) * 1.05, 40)
for sc, (color, label) in colors.items():
peaks = get_peaks_intv(intv_scenarios[sc])["peak"].to_numpy()
ax.hist(peaks, bins=bins, alpha=0.5, color=color,
edgecolor="none", label=f"{label} (med {np.median(peaks):.0f})")
ax.set_xlabel("Peak infected count")
ax.set_ylabel("Frequency")
ax.legend(frameon=False, fontsize=9)
ax.spines[["top", "right"]].set_visible(False)
plt.tight_layout()
plt.show()
Several lessons emerge from this comparison:
Timing is everything for reactive interventions. Isolation starting on day 4 still allows three days of exponential growth. By then, dozens of boys are already infectious and have seeded transmission chains that isolation can’t fully contain. This is a fundamental constraint: with influenza’s 1–2 day incubation period, a substantial fraction of transmission occurs before symptoms appear (Carrat et al., 2008), setting a ceiling on how much symptom-based isolation can achieve.
The herd immunity threshold is demanding. With \(R_0 \approx 3.8\) for this parameterization (at the high end of estimates for this outbreak; Murray 2002 and Vynnycky & White 2010 give \(R_0 \approx 3.0\)–\(3.5\)), the herd immunity threshold is \(1 - 1/R_0 = 0.74\). To reach this with a 70%-effective vaccine, you’d need to vaccinate \(0.74 / 0.70 = 106\%\) of the population — impossible. Even 80% coverage (56% effectively immune) doesn’t cross the threshold, though it substantially blunts the epidemic.
Layered interventions compound. Neither vaccination at 60% coverage nor isolation alone brings \(R_0\) below 1. But the combination reduces the effective reproduction number from 3.8 to roughly 1.3–1.5. The epidemic is nearly eliminated — not because either intervention is sufficient on its own, but because they multiply. This is the Swiss-cheese model of public health: each layer is imperfect, but stacking them drives transmission toward extinction.
The observed data already reflects partial isolation — sick boys were confined to bed. The “baseline” scenario uses parameters fitted to this partially-mitigated epidemic. A truly unmitigated outbreak would have an even higher effective \(\beta\). The counterfactual comparisons above ask “what if we had done more” — the more interesting and policy-relevant question.
The SIR here is deliberately simple — one pathogen, one population, fourteen days. The same tools scale to substantially more complex models. camdl supports indexed parameters, spatial stratification with coupling matrices, time-varying covariates via forcing functions, and observation models with multiple data streams. In the He et al. measles vignette, we apply the same workflow — --draws prior, --replicates, --scenario, batch — to a model with seasonal forcing, environmental overdispersion, and 60 years of weekly case data. The commands are the same; the model is not.