CICADAS

CICADAS (Causal Inference for Critical-Care Anti-Seizure Treatment with Disease Dynamics, Automated PKPD/Trial Simulation, and Survival Optimization) is a simulation-based framework for emulating randomized trials of anti-seizure treatment in critically ill ICU patients, combining the parametric g-formula with mechanistic pharmacokinetic–pharmacodynamic models and disease-dynamics simulations.

This repository accompanies:

McCauley, Westover T., Wynn, Sartipi, et al. CICADAS: A Simulation-Based Framework for Randomized Trial Emulation in Critical-Care Seizure Treatment. iScience, in revision, 2026.

Repository layout

cicadas/
├── matlab/                           # MATLAB analysis pipeline
│   ├── a0_GenerateTrialData.m        # Generate RCT + observational cohorts
│   ├── a0_GenerateDoseSwitchingData.m# Dose-switching cohort (for PKPD estimation)
│   ├── a1_EstimatePKPD.m             # PKPD parameter estimation (EKF/RTS/EM, 2-stage w/ bias correction)
│   ├── a2_CausalSurvivalAnalysis.m   # g-formula causal survival + bootstrap CIs
│   ├── a3_ThreeTreatmentTargets.m    # Compare three treatment targets
│   ├── a4_HeatMap_Agressive.m        # Heterogeneous-effects heatmap (aggressive)
│   ├── a4_OptimalTreatmentTarget.m   # Personalized optimal threshold (parallel bootstrap)
│   ├── a4_Optimize_Heatmap.m         # Personalized optimization heatmap
│   ├── fcn*.m                        # Helpers (PKPD, disease dynamics, plotting, bootstrap)
│   └── run_all.m                     # Orchestrator
├── CICADA_FIGURES/                   # MATLAB figure-generation scripts (one per manuscript figure)
│   └── (a1_SingleTraces, a2_EvaluatePKPD_estimates_figures, a3a/b/c swimmers,
│        a4_HeatMaps_Combined, a5_OptimizationCurve, run_all_figures, helpers)
├── python/                           # Python port (parity-verified; see PARITY_REPORT.md)
│   ├── a0/a1/a2/a3/a4 analogues of matlab/*.m
│   ├── fcn*.py                       # Helpers
│   ├── run_all.py                    # Pipeline orchestrator
│   ├── test_smoke.py                 # Quick import + tiny-N sanity checks
│   ├── README.md                     # Python-specific usage
│   └── requirements.txt
├── CICADA_FIGURES_PY/                # Python figure-generation scripts (ports of CICADA_FIGURES/*.m)
├── sensitivity/                      # Sensitivity analyses (unmeasured confounding, measurement error, alt censoring)
│   ├── alt_censoring/run_alt_censoring.{m,py}
│   ├── measurement_error/run_measurement_error.{m,py}
│   └── nuc_injection/run_nuc_sensitivity.{m,py} + fcnSimulate_N_Patients_withU
├── benchmarks/                       # IPTW / MSM benchmark implementations (Python)
├── matlab_outputs/                   # MATLAB reference pipeline outputs (for parity verification)
├── python_outputs/                   # Python pipeline outputs (for parity verification)
├── verify_parity.py                  # Automated MATLAB-vs-Python comparison harness
├── run_all_verify.sh                 # Bash driver to run MATLAB pipeline + snapshot outputs
├── run_all_verify.m                  # (alternate MATLAB wrapper, now superseded by .sh)
├── PARITY_REPORT.md                  # Full per-artifact parity table (generated by verify_parity.py)
├── docs/
│   └── figure_map.md                 # Script → manuscript-figure mapping
├── parmsTrue.mat                     # Ground-truth simulation parameters (tracked)
├── LICENSE                           # MIT
├── CITATION.cff                      # Citation metadata
└── README.md

Requirements

• MATLAB R2024b or R2025b. The paper was generated under R2024b; R2025b produces qualitatively identical but byte-divergent results due to minor floating-point drift in fmincon and the EKF. For byte-level reproducibility, pin R2024b.
  • Optimization Toolbox (fmincon)
  • Statistics and Machine Learning Toolbox (pooled logistic regression, bootstrap)
  • Parallel Computing Toolbox (parfor is used by a4_OptimalTreatmentTarget.m)
• Python ≥3.10 for the Python port and for the sensitivity/benchmarking analyses added in the revision. See python/requirements.txt.
  • numpy, scipy, pandas, matplotlib, statsmodels, lifelines, h5py (optional, for MATLAB v7.3 .mat files)

Quickstart (MATLAB)

From the repository root, add all subdirectories to the MATLAB path and then run the orchestrator:

addpath(genpath(pwd));

% Full end-to-end pipeline (includes a0 -> a4 -> figure generation):
run('matlab/run_all.m')

Or step through manually (same path setup):

addpath(genpath(pwd));

% 1. Generate simulated cohorts
a0_GenerateTrialData                 % ~25 s  -> trialData0.csv, trialData1.csv, parmsTrue.mat
a0_GenerateDoseSwitchingData         % ~90 s  -> trialDataDoseChanging.csv

% 2. Estimate PKPD parameters
a1_EstimatePKPD                      % ~85 s  -> PKPD_estimation_results.mat

% 3. g-formula causal survival analysis with bootstrap CIs (1000 replicates)
a2_CausalSurvivalAnalysis            % ~40 min -> bootstrap_confidence_bands.mat

% 4. Treatment-target analyses
a3_ThreeTreatmentTargets             % ~15 s  -> ThreeCurves.mat
a4_HeatMap_Agressive                 % ~3 min -> HeatMapAggressive.mat
a4_OptimalTreatmentTarget            % ~1 h   -> A01Data.mat       (parallel)
a4_Optimize_Heatmap                  % ~45 min -> HeatMapData.mat

% 5. Regenerate figures
cd CICADA_FIGURES
run_all_figures

Quickstart (Python)

From the repository root:

pip install -r python/requirements.txt

# Full end-to-end Python pipeline (a0 -> a4 -> figures -> sensitivity):
python3 python/run_all.py

# Or step through stages individually:
python3 python/run_all.py --only a2_CausalSurvivalAnalysis
python3 python/run_all.py --skip-figures --skip-sensitivity

# Quick sanity check (imports + tiny-N sims, ~5 s):
python3 python/test_smoke.py

Individual Python analysis scripts can also be run directly — they read/write from the repo root the same way the MATLAB scripts do (trialData*.csv, parmsTrue*.mat, etc.).

Expected paper results (from a0 and a2)

Quantity	Paper	MATLAB (R2025b)	Python
RCT ATE at 168h	+14.8%	+14.8%	+14.8%
Naive Kaplan–Meier ATE	−4.8%	−4.8%	−4.9%
g-formula ATE	+14.2%	+14.4%	+14.1%
Treated survival at 168h	46.3%	46.6%	44.6%
Untreated survival at 168h	32.1%	32.2%	30.5%

Python and MATLAB agree on every paper-headline quantity within the ok stochastic-parity tier (mean-z < 0.5, sign match, quantile-relative < 20%). See PARITY_REPORT.md for the full per-artifact comparison.

Parity verification

# 1. MATLAB reference run (writes to matlab_outputs/)
bash run_all_verify.sh

# 2. Python pipeline (writes to python_outputs/)
python3 python/run_all.py

# 3. Compare (regenerates PARITY_REPORT.md)
python3 verify_parity.py

Tolerance tiers: deterministic (abs < 1e-8); optimizer (rel|Δmean| < 1%); stochastic (mean-z < 0.5, sign match, quantile-relative < 20%).

Reproducibility notes

• RNG seeding. a0_GenerateTrialData calls rng(0) / np.random.seed(0) and produces byte-identical CSVs within each language. Cross-language bit-equal parity is NOT achievable because MATLAB uses Mersenne Twister and NumPy uses PCG64; the parity harness compares statistical summaries, not bit equality.
• Bootstrap replicates. Both a2_CausalSurvivalAnalysis and a4_OptimalTreatmentTarget use Nboot = 1000 in both languages.
• Paper figures. The authoritative figures in the manuscript were generated under MATLAB R2024b. Regeneration under R2025b or via Python produces visually equivalent figures with minor per-pixel drift; Fig_Combined_PKPD_Analysis shows visible drift in Joint-method estimation errors under R2025b due to floating-point changes in fmincon.
• Known residual drift. a4_OptimalTreatmentTarget threshold-sweep ATE magnitudes at extreme θ (θ<0.1) differ ~2× between MATLAB and Python (both find the same optimal θ* = 0.02 and same qualitative shape). Not affecting the paper's primary endpoint; documented in PARITY_REPORT.md.
• nuc_injection sensitivity. Defaults to "quick" mode (N=400, 1 seed/cell, ~6 min). Set CICADAS_NUC_FULL=1 to run the full 5×5×3-seed N=1500 grid (~3–5 hours).

License

MIT (see LICENSE).

Citation

See CITATION.cff.

Contact

M. Brandon Westover, MD, PhD — mwestove@bidmc.harvard.edu Brain Data Science Platform, Beth Israel Deaconess Medical Center / Harvard Medical School