We validated the normalized relative risk (nRR) approach using individual-based
epidemiological simulations implemented with the nosoi package (v1.1.3) in R. The
standard scenario (two_cluster) models pathogen transmission across a structured
metapopulation of six discrete demes arranged as two spatial clusters, connected by a
single inter-cluster route. To test nRR's ability to detect edge-specific changes in
connectivity, we compared a null model with no intervention against two staged
movement-restriction policies: one applied uniformly to every route, and one that held
the inter-cluster route under tighter restriction than the rest of the network.
Six demes were arranged as two fully-connected triangles, A–B–C (cluster C1) and X–Y–Z (cluster C2), joined by the single direct edge A ↔ X. There is no intermediate bridge deme: A and X are each other's only route between clusters, and each is also an ordinary member of its own triangle.
B----C Y----Z
\ / \ /
A ---- 0.30 ---- X
movement_network.pdf in the repository root renders this structure with the baseline
probability on every directed edge; it is regenerated by
Rscript scripts/two_cluster/two_cluster_network_plot.R.
Hosts moved at each time step with a baseline probability of PROB_MOVE = 0.10.
Conditional on movement, a host's destination was drawn according to the following
symmetric weight matrix (rows sum to 1):
| → A | → B | → C | → X | → Y | → Z | |
|---|---|---|---|---|---|---|
| A: | — | 0.35 | 0.35 | 0.30 | 0 | 0 |
| B: | 0.50 | — | 0.50 | 0 | 0 | 0 |
| C: | 0.50 | 0.50 | — | 0 | 0 | 0 |
| X: | 0.30 | 0 | 0 | — | 0.35 | 0.35 |
| Y: | 0 | 0 | 0 | 0.50 | — | 0.50 |
| Z: | 0 | 0 | 0 | 0.50 | 0.50 | — |
The intervention is defined on per-edge movement rates,
rate_ij = PROB_MOVE × P(destination j | move), rather than on the destination weights.
This is what makes an edge-specific policy expressible: scaling a single edge changes
both that deme's total mobility and its destination mix, in the correct proportions.
Suppressed travel is lost, not redirected — restricting the A–X edge lowers A's total
mobility and leaves the absolute A→B and A→C rates untouched — so any change observed on
within-cluster edges is not an artefact of the bridge policy.
Movement followed a five-phase schedule, expressed as multipliers on the baseline per-edge rate:
| Phase | Start (time step) | All edges | A–X bridge (heterogeneous) |
|---|---|---|---|
| burn-in | 0 | 1.00 | 1.00 |
| shutdown | 100 | 0.10 | 0.10 |
| initial | 150 | 0.33 | 0.10 |
| moderate | 200 | 0.66 | 0.33 |
| full | 250 | 1.00 | 0.50 |
Three scenarios were run:
none— the null model; baseline movement rates throughout, no intervention.homogeneous— the schedule above applied identically to every edge.heterogeneous— the same schedule, but restoration of the A–X bridge lags a full phase behind and never recovers past half of baseline.
Because nosoi fixes structure.matrix for the duration of a run, time-varying movement
was implemented by patching nosoi:::moveFunction.discrete at load time so that the
destination matrix is recomputed from the current simulation time on every move; the
patch fails loudly if the upstream function signature changes.
Transmission dynamics included density-dependent suppression via a logistic carrying capacity function. The per-contact transmission probability at each time step was drawn as:
p_trans = Beta(α=1, β=4) × capacity_factor(n)
where capacity_factor(n) = 1 / (1 + exp((n − 150) / 30)) is a logistic decline in
transmissibility as the local active host count n approaches the carrying capacity
midpoint of 150 (out of MAX_CAPACITY = 300). Transmission was set to zero during the
incubation period. Transmission is identical in every deme and every phase — movement is
the only quantity the intervention touches. The full parameter set was:
| Parameter | Value | Description |
|---|---|---|
INF_ALPHA |
1 | Beta distribution shape (α) for transmission probability |
INF_BETA |
4 | Beta distribution shape (β) for transmission probability |
MAX_CAPACITY |
300 | Logistic carrying capacity (midpoint 150, scale 30) |
CONTACT_MEAN |
4 | Mean contacts per time step |
CONTACT_K |
0.2 | Negative-binomial dispersion for contacts (overdispersed) |
PROB_EXIT |
1/7 ≈ 0.143 | Per-step probability of recovery/exit (mean infectious period = 7 steps) |
INCUB_MEAN |
3 | Mean incubation period (time steps) |
INCUB_SD |
0.5 | Standard deviation of incubation period |
PROB_MOVE |
0.10 | Baseline per-step movement probability |
Contacts were drawn from a negative binomial distribution (size = 0.2, mu = 4) rather
than a Poisson, so that transmission is strongly overdispersed and superspreading events
dominate — the regime in which identical-sequence clustering is informative.
The incubation period was drawn from a gamma distribution parameterized by method of
moments: shape = (mean/sd)² = 36, rate = shape/mean = 12.
Each simulation was initialized with a single infected host in deme A (a gateway deme of
cluster C1) and run for up to MAX_TIME = 350 time steps or MAX_OUTBREAK = 120,000
total infections. An attempt was accepted only if it reached MIN_INFECTED = 200 total
hosts and MIN_INFECTED_BURN_IN = 100 infections before the intervention began at
time step 100; otherwise it was re-run (up to 50 attempts). Judging acceptance on the
burn-in period only is deliberate: every scenario is identical before the intervention,
so all three accept or reject exactly the same attempts and remain genuinely paired by
seed. Testing anything measured after the intervention starts would let the more
restrictive scenarios reject more often and silently accept a different attempt with a
different burn-in, breaking the cross-scenario comparison. Each scenario was replicated
across 25 independent random seeds by default.
For each simulation replicate, the resulting transmission tree (Newick format) was passed
to phastSim to simulate nucleotide sequence evolution under a molecular clock, using
the SARS-CoV-2 Wuhan-Hu-1 reference genome (data/reference/wuhan.fasta) as a starting
sequence. A fixed substitution rate of 2 × 10⁻⁶ substitutions per site per day
(calibrated for SARS-CoV-2) was applied via the --scale 0.000002 flag, and the phastSim
seed was held constant at 17 across all replicates to isolate stochasticity to the
epidemiological component.
Because pairwise SNP calculation is O(n²), the simulated FASTA was first thinned to a
target tip count (default 8,000) using scripts/common/subsample_fasta.py. The sampling
budget was divided evenly across infection-time bins (width 10 time steps) rather than
uniformly over the file: uniform sampling tracks incidence and leaves the early epidemic
with too few tips for a stable RR estimate. Sampling remains uniform within each bin, so
deme composition inside a time window is unbiased.
Pairwise genomic distances between the subsampled sequences were then computed using
pairsnp. Only sequence pairs with zero SNP distance (i.e., genetically identical
pairs, -d 0) were retained, as these are the most informative signal for recent shared
transmission history.
From the identical-sequence pairs, we computed a geographic relative risk (RR) for each ordered deme pair (X, Y):
RR(X, Y) = (n_XY × N) / (n_X × n_Y)
where n_XY is the number of identical pairs with one sequence from deme X and one from
deme Y, n_X and n_Y are the total identical pairs involving deme X and Y respectively,
and N is the total count across all deme pairs. RR > 1 indicates more genetic similarity
between demes X and Y than expected under random mixing.
To control for variation in within-deme transmission intensity—which can inflate or deflate raw RR values independently of between-deme connectivity—we computed the normalized RR (nRR):
nRR(X, Y) = RR(X, Y) / mean(RR(X, X), RR(Y, Y))
nRR is thus normalized by the arithmetic mean of the two relevant within-deme RR values.
RR and nRR were computed at two spatial resolutions — the full 6×6 deme matrix, and a collapsed 2×2 cluster matrix (C1 vs. C2) — and under two temporal decompositions:
- Per phase. Because the phase boundaries are fixed and known, all pairs falling within a phase were pooled into a single estimate. This is substantially more stable than taking a median over narrow sliding windows, which collapses to zero whenever more than half the windows happen to be empty.
- Sliding window. A window of width ±20 time steps advancing every 5 time steps, to trace nRR continuously and locate the timing of connectivity changes.
Recovery relative to baseline (the recovery_by_phase figures) divides each pair's phase
nRR by that same pair's sliding-window nRR at tpoint = 90, i.e. the window [70, 110]
bracketing the start of the intervention. A ±20 window is wide enough to keep the sparse
inter-cluster pairs from collapsing to zero; strictly pre-intervention windows lose too
many of them to emptiness.
The burn-in phase is excluded from these figures. The seeding period is intrinsically dense in identical sequences -- one introduction, no standing diversity -- so it sits high for reasons unrelated to the intervention and is not comparable to the later phases.
Recovery ratios are plotted on a pseudo-log axis (sigma = 0.005), which is log-like
above sigma and linear through zero. A phase in which a deme pair has no identical
sequences at all yields a ratio of exactly 0 -- a real result, and the most extreme one --
which a plain log axis silently discards. This matters most for the heterogeneous arm,
whose inter-cluster ratio is 0 for the majority of seeds during shutdown and initial.
All simulations were implemented in R (≥4.0) using the following packages: nosoi
(individual-based transmission simulation), ape and treeio/ggtree (phylogenetic tree
handling and visualization), tidyverse (data manipulation), igraph/ggraph
(transmission network plotting), and argparse (command-line argument parsing). Sequence
evolution was simulated with phastSim and pairwise SNP distances were calculated with
pairsnp. Analysis scripts were executed in a reproducible conda environment
(ncov-sim) using Python 3.12.
conda env create -f ncov-sim-env.yml
conda activate ncov-sim
Rscript scripts/common/install_packages.RRun a single replicate of the standard scenario, with diagnostic plots:
Rscript scripts/two_cluster/two_cluster_sim.R --seed 42 --scenario heterogeneous --draw
bash scripts/two_cluster/two_cluster_seq.sh 8000
Rscript scripts/two_cluster/two_cluster_RR.R --seed 42 --scenario heterogeneousOr run a full sweep — two_cluster_batch.sh [SCENARIO] [N_SEEDS] [SUBSAMPLE_N], which
loops sim → phastSim → subsample → pairsnp → RR over every seed and ends by producing the
summary plots:
bash scripts/two_cluster/two_cluster_batch.sh none 25 8000
bash scripts/two_cluster/two_cluster_batch.sh homogeneous 25 8000
bash scripts/two_cluster/two_cluster_batch.sh heterogeneous 25 8000Run the batch script once per scenario. It clears only the accumulators belonging to the
scenario it was given, and every output file is prefixed with the scenario name, so
sequential sweeps do not overwrite each other. The --scenario value must be passed
consistently to _sim.R, _RR.R, and _plot.R: it selects the intervention and
namespaces the output files.
Once all three sweeps are done, two figures are produced outside the per-scenario loop:
Rscript scripts/two_cluster/two_cluster_combined_plot.R # results/two_cluster/recovery_by_phase_combined.pdf
Rscript scripts/two_cluster/two_cluster_network_plot.R # movement_network.pdftwo_cluster_combined_plot.R stacks the three per-scenario recovery panels into one 6x8
inch PDF on a shared axis, labelled No Shutdown, Uniform Recovery, and Staggered
Recovery. It reads the *_recovery_by_phase.tsv tables, so all three scenarios must have
been plotted first; a missing scenario is skipped with a warning rather than failing.
All scripts are written to be run from the repository root, e.g.
Rscript scripts/two_cluster/two_cluster_sim.R --seed 42. Paths inside the scripts are
relative to the root, so never cd into a scenario folder.
scripts/
common/ # shared helpers: calc_rr.R, plot_pair_time.R,
# transmission_pair_analysis.R, sim_analysis.R,
# subsample_fasta.py, install_packages.R
two_cluster/ # six-deme two-cluster model (standard scenario)
basic/ # single-population baseline scenario
shutdown_threebody/ # three-deme model with intervention / movement shutdown
carrying_capacity/ # density-dependent transmission (+ seasonality, isolation, immunization)
movement_network.pdf # movement-structure diagram (model property, not a run output)
data/
reference/ # wuhan.fasta reference genome (tracked)
results/
two_cluster/ # scenario outputs: trees, metadata, sequences, RR tables, plots
basic/
shutdown_threebody/
carrying_capacity/
Scenario folders follow a common naming convention, <scenario>_* (not every scenario has
all five):
| Suffix | Purpose |
|---|---|
_sim.R |
nosoi simulation; writes tree + metadata to results/<scenario>/ |
_seq.sh |
single-replicate phastSim + (subsample) + pairsnp run |
_RR.R |
relative risk / nRR calculation, appends to results/<scenario>/*rr*.tsv |
_plot.R |
summary plots across replicates |
_batch.sh |
full sweep: sim → phastSim → pairsnp → RR → plots |
Two two_cluster scripts sit outside that convention because they are not per-scenario:
two_cluster_combined_plot.R (all three arms in one figure) and
two_cluster_network_plot.R (the movement diagram, which depends only on the model
definition and not on any run).
Simulation outputs under results/ are gitignored; the directories are recreated by the
scripts as needed. data/ holds only tracked static inputs.