GoldsteinGameTheory

Host–pathogen coevolution simulator based on the Goldstein et al. framework for punctuated evolutionary dynamics.

The model asks how evolving conditional response rules (evolved rules, ER) rather than fixed trait values (evolved traits, ET) shapes long-term host–pathogen coevolution. Host clearance c and pathogen virulence v evolve by a continuous-time Markov jump process under strong selection and weak mutation (SSWM), with Kimura fixation probabilities. Four scenarios cross the two strategy types: ET/ET, ET host / ER path, ER host / ET path and ER/ER.

The code is based on J. Goldstein's original Java implementation (Simulate.java, kept for reference) and the original paper draft.

Repository layout

GoldsteinGameTheory/
├── src/                     # the model
│   ├── simulation.py            engine (importable module)
│   ├── run_experiments.py       CLI: configure and launch runs
│   └── timeshift.py             time-shift (cross-inoculation) assay
├── scripts/                 # launchers that reproduce every simulation
│   ├── run_all.sh               the whole pipeline in dependency order
│   ├── run_main.sh  run_zoom.sh  run_tempo_sweep.sh
│   ├── run_tracking_sweep.sh  run_timeshift.sh
│   └── missing_runs.py          prints commands for runs not yet complete
├── analysis/                # everything that reads runs and writes figures
│   ├── Plots.R                  plotting and analysis helpers (definitions only)
│   ├── manuscript_figures.R     builds every manuscript figure
│   ├── timeshift_plots.R        helpers for Fig 5
│   ├── si_extra_figures.R       Figs S12 (escalation) and S13 (sawtooth)
│   ├── si_tables.R              run inventory and per-replicate tables
│   ├── anchor_comparison.R      equilibrium-rule comparison
│   ├── manuscript_numbers.py    every simulation-derived number in the text
│   ├── diag_equilibrium.py      which equilibrium the solver realises
│   └── excursion_attribution.py what drives the punctuated jumps
├── tests/
│   ├── test_root_finder.py      enumeration vs a dense scan, 1e6 parameter sets
│   └── test_dimensional_weight.py  51^2 units of weight per player, every scenario
├── data/                    # simulation output (git-ignored, regenerated)
│   ├── main/                    minimal, acute and tracking runs
│   ├── zoom/                    every-substitution runs
│   ├── tempo/                   gamma sweep, N_H = N_P
│   ├── random/                  robustness: random Gaussian proposals
│   ├── rule_mutant/             robustness: superseded selection rule
│   └── timeshift/               assay summaries written by src/timeshift.py
├── output/                  # git-ignored, regenerated
│   ├── figures/                 manuscript figures (PDF + PNG)
│   └── tables/                  supplementary tables (CSV)
├── superseded/              # runs kept for provenance (see its README)
└── reference/Simulate.java  # original Java implementation

Run every command from the repository root, e.g. Rscript analysis/manuscript_figures.R fig3 or scripts/run_main.sh. Data and output paths are relative to it.

Simulation output and figures are git-ignored; they are regenerated by the steps below.

Requirements

• Python 3, standard library only for everything under src/, scripts/ and tests/ (tested with 3.14). scipy is optional and changes nothing: without it the inverse error function used to build the mutation grid is computed by bisection. analysis/manuscript_numbers.py additionally needs numpy.
• R with tidyverse, ggplot2 ≥ 4.0, patchwork, pracma, scales, zoo and jsonlite (tested with R 4.6, ggplot2 4.0.3, patchwork 1.3.2).

Run every command from the repository root.

Reproducing the manuscript

1. Simulations

Every launcher runs jobs in parallel (JOBS=…), writes one log per run to logs/, and uses deterministic seeds (seed_base + rep − 1). All runs are diploid with mutation step size σ = 0.01.

Script	Output	Design	Used in
scripts/run_main.sh	results/minimal/, results/acute/	2 models × 4 scenarios × 4 reps; 100K substitutions after 10K burn-in, one row per 10	Fig 2–4, S1, S3–S6, S10–S11
scripts/run_tracking_sweep.sh	results/tracking/	k ∈ {0, 0.5, 1, 1.5, 2, 3, 4} × 4 scenarios × 4 reps; 10K substitutions after 10K burn-in	S8–S9
scripts/run_tempo_sweep.sh	results_gamma/	γ ∈ {10⁻⁴, 10⁻², 0.5, 0.99} × 3 scenarios × 4 reps, N_H = N_P = 10⁴; 10K after 10K burn-in	Fig 4D–E, S7
scripts/run_zoom.sh	results_zoom/	minimal model, 4 scenarios × 8 reps; 10K substitutions, every one recorded	Fig 4C, S2, Fig 5

scripts/run_all.sh runs the whole pipeline in dependency order; scripts/missing_runs.py prints the commands for any run in the design that is not already complete on disk, which is the resume-safe way back after a code change that invalidates part of it:

scripts/run_all.sh                                        # everything, in order
python3 scripts/missing_runs.py --count                   # what is missing
python3 scripts/missing_runs.py | tr '\n' '\0' | xargs -0 -n 1 -P 10 bash -c

Two robustness sets sit outside the main design, each run the same way with one flag changed: --proposals random (random Gaussian proposals instead of the deterministic grid, into results_random/) and --eq-selection mutant (the superseded selection rule, into results_rule_mutant/).

ER/ER runs are by far the slowest: each ER player evaluates 51 × 51 candidate rules per substitution, against 51 candidate traits for an ET player.

With equal population sizes the tempo ratio R = (1 − γ)N_P / (γN_H) is set by γ alone (R ≈ 10⁴, 10², 1). At R = 1 the two mixed scenarios must mirror each other with v and c exchanged, which S7 uses as a code check.

2. Time-shift summaries (Fig 5)

After the main and zoom runs have finished:

scripts/run_timeshift.sh main   # results/timeshift/timeshift_main*.csv (panels D–E)
scripts/run_timeshift.sh zoom   # results/timeshift/timeshift_zoom*.csv (panels A–C)

3. Figures

Rscript manuscript_figures.R              # all 16 figures -> figures/
Rscript manuscript_figures.R fig3 figS8   # selected figures

If several large PDFs in one R session hit an intermittent "write failed" error, build each figure in its own process: for f in fig1 fig2 fig3; do Rscript manuscript_figures.R $f; done.

Figure	Content	Data
fig1	Model: fitness landscapes and response-rule geometry	analytical
fig2	Time series of v, c, W_P, W_H, four scenarios	results/minimal
fig3	Mechanism of destabilisation: neutral fraction, genotype/phenotype ratio, own vs opponent-induced step sizes, slope phase space, Nash violation map, time-weighted state occupancy	results/minimal
fig4	Selection and tempo: ω/2N distributions and per-replicate medians, aligned zoom, tempo sweep	results/minimal, results_zoom, results_gamma
fig5	Time-shift assay	results/timeshift
figS1	Trait and fitness density hexbins	results/minimal
figS2	Time series, first 1,000 substitutions	results_zoom
figS3	Time-series statistics (mean, SD, spectral slope, correlation length)	results/minimal
figS4	Response-rule snapshots, ER/ER	results/minimal
figS5	Response-rule parameter time series, ER/ER	results/minimal
figS6	ω/2N time series, full runs	results/minimal
figS7	Tempo symmetry check at R = 1	results_gamma
figS8	Tracking sweep: virulence SD, corr(W_H, W_P), ER-host fitness relative to ET/ET	results/tracking
figS9	Tracking-model ER/ER time series per k	results/tracking
figS10	Acute-model time series	results/acute
figS11	Acute-model diagnostics (analogues of Fig 3C, 3F, 4B)	results/acute

Code

simulation.py

• Fitness models: acute, chronic, minimal, Taylor et al. (2006) and tracking (see the table below). set_fitness_model() switches the active functions and trait domain.
• Simulation class: evolutionary state, mutation proposals (equal-probability Gaussian quantile bins), selection by Kimura fixation probability (haploid or diploid), and the continuous-time substitution clock.
• Gillespie dynamics: at each step every candidate mutation of both players is scored. The player that substitutes is chosen in proportion to its total rate, weighted by the mutation asymmetry γ. The dwell time is exponential.
• Equilibria: for ER players, the realised traits are a mutual fixed point of the two clamped linear response rules. The solver enumerates interior and boundary fixed points. When several exist (only possible when m_c·m_v > 1), the realised one is the equilibrium the behavioural dynamics reach from the phenotype expressed before the substitution (EQ_SELECTION).

Notation (paper → code):

Paper	Code	Description
ET (evolved trait)	evolved_strategy=False	Fixed trait
ER (evolved response rule)	evolved_strategy=True	Linear reaction norm
c	s	Host clearance
v	v	Pathogen virulence
c₀, m_c	bS, mS	Host rule intercept and slope, c = c₀ + m_c v
v₀, m_v	bV, mV	Pathogen rule intercept and slope, v = v₀ + m_v c
W_H, W_P	hostFit, pathFit	Host and pathogen fitness
ω_H, ω_P	omegaHost, omegaPath	Substitution rate, 2N·P_fix (ceiling 2N)
γ	prob_host_mutate	Probability that a mutation arises in the host

Modelling choices that were corrected. Four changes altered the model's behaviour, in roughly increasing order of consequence. Runs made before each are kept in superseded/, and every run records the settings it used in its config.json, so old and new output cannot be mixed by accident.

1. Fixation probability. For a strongly beneficial mutation kimura_fixation_prob returns the limit 1 − e^(−2s). An earlier version returned the small-s approximation 2s, which is unbounded: with s measured against a near-zero-fitness resident it let P_fix exceed 1 and ω exceed its 2N ceiling.
2. Mutation grid (STEP_GRID, default quantile). The 51 steps are the conditional means of the standard normal within equal-probability bins. Until scipy was removed as a dependency, a missing scipy silently substituted evenly spaced steps on ±3σ, whose SD is 1.77 rather than 1.0 — so σ behaved 1.77× larger than stated, and the imbalance scaled with the grid size. check_step_grid() now asserts the grid's mean and SD at every run.
3. Fixed-point enumeration (ROOT_FINDER: exact-piecewise). The composed clamped map is piecewise affine with at most four interior breakpoints, so every fixed point is solved in closed form. This replaced a 256-point sign-change scan that missed roots whenever two lay within one grid cell (~0.7% of states), and is about 20× faster. Checked by test_root_finder.py against a dense scan over 10⁶ parameter sets per model.
4. Dimensional weighting (DIM_WEIGHT: per-player). A player whose slope is pinned proposes 51 trait steps against a plastic partner's 51 × 51, so it carries a ×51 weight — as in Simulate.java, where sFactor is set from fixedS alone. The port had applied it only when both players were ET, leaving the pinned player in a mixed scenario with 51× less mutational supply than its partner. Checked by test_dimensional_weight.py, which asserts 51² units of weight per player in all four scenarios.

Equilibrium selection (EQ_SELECTION, default anchor). Where the paired rules admit several fixed points — only possible when m_c·m_v > 1, about 5% of ER/ER states — the realised one is what the behavioural dynamics reach from the phenotype expressed before the substitution, so the interaction carries over a substitution and neither player chooses its own phenotype. Setting it to mutant reproduces the original rule, in which the mutating player took the fixed point best for itself; results_rule_mutant/ holds that comparison.

run_experiments.py

Command-line runner. It sets the simulation.py globals, names the output directory and saves the full configuration to config.json.

python3 run_experiments.py --list                                    # models and scenarios
python3 run_experiments.py -f minimal -c ERhost_ERpath --diploid \
    --step-sizes 0.01 --rep 1 --gens 100000 --write-every 10         # one main run
python3 run_experiments.py -f minimal --diploid --gens 1000 --bins 11 \
    --tag test -o scratch                                            # quick test

Flag	Description
-f/--fitness	acute (default), chronic, minimal, taylor, tracking
-c/--condition	One scenario (default: all four)
--gens N, --quick	Recorded substitutions (default 1M; --quick = 10K)
--burn-in N	Burn-in before recording (default 10K, but --gens alone shrinks it to gens/10)
--write-every N	Record every Nth substitution (default 100)
--step-sizes σ[,σ…]	Mutation step size(s)
--diploid	Diploid semi-dominant fixation (4Ns) instead of haploid (2Ns)
--gamma γ, --gamma-sweep γ,…	Mutation asymmetry (default 0.01)
--pop-size N, --path-pop-size N	Host N (pathogen defaults to 100× host); pathogen N independently
--k k	Tracking strength for the tracking model (0 = minimal)
--rep N	Replicate: offsets the seed by N − 1 and tags the directory repN
--tag name, -o dir	Directory tag; output root (default results)
--fix-host x, --fix-path x	Pin one player's trait
--nS, --nV	Clearance and virulence cost parameters
--bins N	Mutation quantile bins (default 51; 11 is ~20× faster for tests)
--seed, --nash-sweep	Seed override; Nash-dwell parameter sweep

Output directories are named {condition}[_sigma{σ}][_diploid][_gamma{γ}][_fixH{x}][_fixP{x}][_N{n}][_NP{n}][_nS{x}][_nV{x}][_k{k}][_rep{n}][_tag]/ (γ is omitted at its default 0.01).

timeshift.py

Challenges one antagonist from time t against the other from time t + Δ, contrasting same-lineage (sympatric) with different-lineage (allopatric) pairings (Gaba & Ebert 2009). No re-simulation is needed. The phenotype protocol crosses stored realised traits. The rule protocol crosses response rules and re-solves the fixed point. See the module docstring and python3 timeshift.py --help.

Plots.R, timeshift_plots.R, manuscript_figures.R

Plots.R and timeshift_plots.R only define functions. manuscript_figures.R sources them and holds one short function per figure. The loaders in Plots.R discover runs from their config.json, so no paths are hard-coded. Point them at another tree with options(ggt.results_root = "…"), or pass results_root = to a figure function.

source("Plots.R")
fig_timeseries("minimal", sigma = 0.01, diploid = TRUE, replicates = TRUE,
               highlight_rep = 1, show_omega = FALSE)

Scenarios and fitness models

Scenario	Host	Pathogen
EThost_ETpath	Fixed trait c	Fixed trait v
EThost_ERpath	Fixed trait c	Rule v = v₀ + m_v c
ERhost_ETpath	Rule c = c₀ + m_c v	Fixed trait v
ERhost_ERpath	Rule c(v)	Rule v(c)

Model	Host fitness	Pathogen fitness	Trait domain
minimal	c(1−c)(1−v)	v(1−v)(1−c)	[0, 1]
acute	c / (c + d)	v^β / (c + d)	[0, 1]
tracking	c(1−c)(1−v) + k c²(1−c)v	v(1−v)(1−c) + k v²(1−v)c	[0, 1]
chronic	1 / d, immunity modulates virulence	(1−c) v^β / d	[0, 1]
taylor	[c/(v+c)] · [b/(m₀+c)]	vⁿ / (v+c)	[0.001, 20]

d is the mortality function of each model. The tracking model reduces exactly to the minimal model at k = 0. The manuscript uses minimal (main text), acute and tracking (SI).

Key parameters

Parameter	Default	CLI flag	Description
std_dev_move (σ)	0.1	--step-sizes	Mutation step size (manuscript: 0.01)
prob_host_mutate (γ)	0.01	--gamma	Share of mutations arising in the host
HOST_POP_N	10⁴	--pop-size	Host effective population size
PATH_POP_N	10⁶	--path-pop-size	Pathogen effective population size
DIPLOID_KIMURA	False	--diploid	Diploid fixation (manuscript: True)
nS_HLP / nV_HLP	0.1 / 1.0	--nS / --nV	Clearance / virulence cost
TRACKING_K (k)	0	--k	Tracking strength (tracking model only)
num_step_bins	51	--bins	Mutation quantile bins (51² candidates per ER player)

Output files

Each run directory contains:

• config.json: the full configuration, including scenario, model, population sizes, step size, γ, k, replicate, seeds and run length.
• simulation.csv: the trajectory. Each recorded substitution writes two rows. The pre row holds the state before the substitution and the post row the state after, so post − pre is exactly one substitution. Burn-in is not written. Main columns:
  • gen: substitution count.
  • time: evolutionary (continuous) time.
  • dwell (post rows): time spent in the pre-substitution state.
  • mutator: which player substituted, host or path.
  • v, s: realised traits.
  • bV, mV, bS, mS, vAngle, sAngle: rule parameters.
  • pathFit, hostFit: fitness.
  • omegaPath, omegaHost: substitution rates.
  • mutSelCoeff, mutClass: the substitution's selection coefficient and class.
  • nash, stableFP, hostLineMax, pathLineMax: equilibrium diagnostics.

Because gen counts substitutions, the same number of substitutions spans very different evolutionary time across scenarios. Over 100K substitutions in the minimal model, ET/ET covers roughly 10⁴ times more time than ER/ER.

Citation

Goldstein et al. "Punctuated instabilities of the immune system" (in preparation)

Contact

GaTechBrownLab

Canan Karakoc (GitHub: @canankarakoc, canankarakoc@gmail.com)