Date: 2026-03-25

Model: laser.measles.compartmental (daily SEIR, discrete integer compartments)

Goal: Find parameters (R0, seasonality, season_start) that produce sustained biennial (~730-day) endemic measles cycles in a single spatial node

Population: 500K -> 1M

CBR: 30-33 /1k/yr, no vaccination

This investigation used systematic parameter identifiability sweeps at each stage to empirically map the parameter space, detect failure modes before committing to expensive calibration runs, and confirm when a proposed approach was fundamentally incapable of reaching the target -- rather than just poorly tuned.

Finding: Biennial endemic cycles cannot be achieved in the discrete compartmental model without a background importation term. The root cause is a hard incompatibility between the two requirements for biennial dynamics:

1. Epidemics must be large enough to deplete susceptibles sufficiently that recovery takes ~2 years.
2. But those same large epidemics drive the infectious count I to zero between cycles, and I = 0 is an absorbing state in the discrete model (no spontaneous infection generation).

This is not a calibration failure -- it is a structural property of the model. The sweeps made this visible quickly and prevented wasted compute on a full-scale optimization that could not have succeeded.

In SEIR models with vital dynamics and seasonal forcing, biennial cycles arise from a 2:1 seasonal resonance: the seasonal peak in transmission triggers one epidemic every other year rather than every year.

The mechanism requires three ingredients:

The prerequisite condition is that each epidemic must deplete susceptibles deeply enough that it takes approximately 2 years of births to refill the susceptible pool above the epidemic threshold S* = N / R0.

For N = 500K, CBR = 30/1k/yr, R0 = 8:

• Births per year: 15,000
• S* = 62,500
• For a 2-year recovery: S must be depleted by ~30,000 below S*, requiring S to start the epidemic at ~92,500 (~1.48 x S*)

Setup: Start at S = N (fully susceptible), seed 50 infections at tick 0. VitalDynamics + InfectionProcess + StateTracker. No InitEquil, no importation.

Identifiability sweep: 72 combos -- R0 in {6,8,10,12,14,16,18,20}, seasonality in {0.05-0.45}. N = 500K, CBR = 30, CDR = 15.

See: biennial_single_node_sweep.png

Result: 0/72 combos produced any cycling.

The naive-population epidemic is catastrophic: with S ~= N and R0 = 8, the effective reproduction number R_eff = 8. The epidemic sweeps through the population in ~90 days, depleting I to zero. With I = 0 and no importation, VitalDynamics replenishes S indefinitely but transmission never restarts.

The sweep output confirmed this immediately -- every cell showed --- (< 3 peaks detected) -- and ruled out the naive approach without needing to run a calibration optimizer.

User insight: Rather than starting at S = N (which causes an explosive burnout), initialize near the endemic fixed point (S = N/R0) and let susceptibles accumulate for one year before seeding. This should produce a moderate first epidemic, and the I trough may remain above the integer floor.

Setup: InitializeEquilibriumStatesProcess sets S = N/R0 at tick 0. ImportationPressureProcess delivers ~100 seed cases over ticks 365-382 only (a narrow "one-shot" seed, not ongoing importation). InfectionProcess runs throughout with seasonality.

Key diagnostic: Before sweeping, a single traced run at R0 = 16, seas = 0.25 confirmed the approach works qualitatively:

S at tick 365: 45,497  (= N/R0 + 1 year of births, R_eff = 1.46)
I at peak:     ~124 (moderate epidemic, not catastrophic)
I trough min:  1  (stays above the integer floor!)

The I trace shows the first epidemic fires cleanly, declines to a minimum of I = 1, then grows again as S rebuilds -- exactly the endemic cycling behaviour hypothesised. This was a qualitative validation of the approach before running the full sweep.

Identifiability sweep: 72 combos -- same grid as Approach 1, but with the corrected initialization. Peak detection threshold lowered to I > 5 (the endemic epidemics are much smaller than the naive-population ones).

See: fig2_corrected_sweep.png

Result: Cycling was found, but in the wrong period band.

All cycling solutions are sub-annual (~190-230 days) -- roughly semi-annual, not biennial. This is the natural oscillation frequency of the SEIR system near its endemic fixed point, not the biennial resonance driven by large susceptible depletion.

The sweep made this visible as a clean heatmap: cycling exists only in the top rows (high R0) and only at sub-annual periods. The 730-day contour is nowhere in the grid. A calibration optimizer would have been sent searching in empty space.

The sweep revealed a structural incompatibility, not a parameter tuning problem.

The biennial cycle mechanism requires S to accumulate to ~1.5-2x S* between epidemics -- which means each epidemic must be large enough to drive S back down to or below S*. But large epidemics are exactly what drive I toward zero in the discrete model.

Left panel: As the epidemic trigger level (S_start / S*) increases, the estimated trough I falls. For R0 = 8, N = 500K, CBR = 30:

• Near-equilibrium (S_start ~= S*): trough I ~= I* ~= 300, but recovery time is only ~0 days -- no meaningful inter-epidemic gap.
• Biennial zone (S_start ~= 1.5x S*): recovery time ~= 730 days, but trough I approaches the integer floor.
• Extinction zone (S_start >= 1.6x S*): trough I < 1 -> epidemic burnout.

Right panel: The recovery time reaches 730 days only when S_start/S* is in a regime where the estimated trough I is <= 1. The biennial timing window and the integer-floor extinction window overlap completely -- there is no gap between them in the discrete model.

This is an intrinsic property of the discrete compartmental implementation at these demographic parameters (N = 500K, CBR = 30). It is not resolvable by parameter tuning alone.

Rationale: Doubling N doubles the absolute trough I (I* ~= 300 -> ~= 660), which may open a viable region between the biennial timing window and the integer floor. Adding season_start as a calibration parameter addresses the phase alignment requirement for 2:1 resonance.

Setup: N = 1,000,000, CBR = 33, CDR = 15. Same InitEquil + delayed seed. Three parameters calibrated: R0 in [12, 18], seasonality in [0.05, 0.50], season_start in {0, 90, 180, 270} days.

Identifiability sweep: 80 combos (4 x 5 x 4 grid).

Result: Still no biennial cycles. Maximum period achieved: ~338 days.

season_start = 270 (peak transmission in October) is the most favourable phase for extending the inter-epidemic period -- it delays the seasonal forcing peak to autumn, misaligning with the natural post-epidemic S-recovery window and extending the cycle. But even the best case reaches only 338 days -- less than half the 730-day target.

The 11 combos with I never hitting zero (sustained oscillation, zero extinction fraction) are all in the sub-annual regime. The biennial gap persists at N = 1M.

A 500-trial Optuna calibration was planned for this stage but was not run: the sweep clearly showed no biennial region in the parameter space, making the calibration search futile. This is the key value of the identifiability sweep -- it prevented a ~25-minute compute run that could not have found what we were looking for.

At each stage of this investigation, a parameter identifiability sweep was run before any calibration optimizer. This workflow is deliberate and important:

An optimizer like Optuna TPE explores the parameter space efficiently but only within the space it is given. If the target objective (730d biennial period) has no attractor in the parameter space, the optimizer will:

• Converge to the nearest local minimum (sub-annual oscillations, ~200d)
• Report a "best" trial with a score of ~0.7 (far from 0)
• Give no indication that the fundamental regime is unreachable

The sweep, by contrast, maps the entire grid in one pass and makes the absence of the biennial region immediately visible. The heatmaps show clearly that the 730-day contour does not exist anywhere in the plotted space.

The objective function used throughout:

combined_score = biennial_score + 0.5 * regularity_cv + 2.0 * zero_penalty

where:
  biennial_score  = |median_period - 730| / 730
  regularity_cv   = std(intervals) / mean(intervals)
  zero_penalty    = fraction of post-seed ticks with I = 0
  no-cycle penalty: 9.0 if < 3 peaks detected

The zero_penalty term (weighted 2x) is diagnostic as well as objective: a high zero_penalty in sweep results is a direct signal that the integer-floor absorption problem is dominating, independent of how well the period is tuned. This term revealed the extinction structure clearly in the sweep heatmaps.

The discrete compartmental model in laser.measles.compartmental cannot sustain biennial endemic cycles without a background importation term, regardless of:

• Population size (tested N = 500K and N = 1M)
• Birth rate (CBR = 30-33)
• R0 in the empirical measles range (6-20)
• Seasonality amplitude (0.05-0.55)
• Season phase (0-270 degrees)
• Initialization strategy (naive vs. InitEquil + delayed seed)

The constraint is structural: the discrete integer compartments make I = 0 an absorbing state, and biennial cycles require epidemic sizes that drive I to or near 0 in the inter-epidemic trough.

Add a minimum background importation rate (0.01-0.05 cases/1k/yr, applied to the single node only). This is:

1. Epidemiologically defensible: no real-world measles endemic country has truly zero importation.
2. Numerically minimal: at 0.01/1k/yr in a 500K city, importation contributes ~0.014 cases/day -- negligible compared to epidemic FOI, and far below what could drive wavefront timing in a spatial model.
3. Demonstrated to work: the earlier biennial wavefront biweekly calibration (calib_bwk.db, best score = 1.55) used ImportationPressure and successfully produced biennial endemic cycling.

With minimal importation, the identifiability sweep can be re-run expecting to find a genuine biennial region, after which Optuna calibration of (R0, seasonality, season_start) is well-posed and likely to converge.