epidemic-test.js

152 lines · in natura-src-1.0.0.zip as natura/epidemic-test.js · a test

// epidemic-test.js -- the EPIDEMIC as real network epidemiology, every emergence verified by RUNNING (never by
// reasoning): mass conserved, an outbreak above R0=1 and a fizzle below, a curve that turns over BEFORE the
// susceptibles run out (and overshoots), the herd-immunity threshold, targeted-beats-random immunisation,
// flatten-the-curve, waning immunity going endemic, and mortality tracking the attack rate.
const NAT = require('./natura.js');
const { makeEpidemic } = require('./epidemic.js');
const { World, UNITS } = NAT;

let pass = 0, fail = 0;
const ok = (c, m) => { if (c) pass++; else { fail++; console.log('  FAIL: ' + m); } };

function run(seed, cfg, days, mutate) {
  const eco = makeEpidemic(cfg || {});
  const w = World({ seed, systems: [eco] });
  const st = w.get('epidemic');
  if (mutate) mutate(st, eco);
  const series = [];
  for (let i = 0; i < days; i++) {
    w.step(UNITS.DAY); const f = w.fields().epidemic;
    series.push({ nI: f.nI, nS: f.nS, nE: f.nE, nR: f.nR, nD: f.nD, attack: f.attackRate, r0: f.r0, reff: f.reff, peakI: f.peakI, peakDay: f.peakDay, sAtPeak: f.sAtPeak, dead: f.dead, susceptible: f.susceptible });
  }
  return { st, f: w.fields().epidemic, series, eco };
}
const peakOf = (s) => s.reduce((m, r) => Math.max(m, r.nI), 0);

console.log('epidemic assertions:');

// 1. DETERMINISM -- same seed, identical trajectory (network + transmission all off one seeded rng).
{
  const a = run(5, { N: 300 }, 90).f, b = run(5, { N: 300 }, 90).f;
  ok(a.nI === b.nI && a.nR === b.nR && a.nD === b.nD && a.attackRate === b.attackRate && a.r0 === b.r0, 'same seed -> identical epidemic (reproducible)');
}

// 2. MASS CONSERVED -- nobody vanishes: S+E+I+R+D == N at every single step. The correctness oracle.
{
  const eco = makeEpidemic({ N: 380 }); const w = World({ seed: 7, systems: [eco] }); const st = w.get('epidemic');
  let held = true; for (let i = 0; i < 120; i++) { w.step(UNITS.DAY); if (st.nS + st.nE + st.nI + st.nR + st.nD !== st.N) held = false; }
  ok(held, 'S+E+I+R+D == N at every step (no one is minted or lost)');
}

// 3. OUTBREAK above R0=1, FIZZLE below -- the epidemic threshold, measured by RUNNING both.
{
  const hi = run(11, { N: 400, beta: 0.12, infDur: 6 }, 160).f;   // high transmissibility
  const lo = run(11, { N: 400, beta: 0.012, infDur: 6 }, 160).f;  // low transmissibility
  ok(hi.r0 > 1.4 && hi.attackRate > 0.4, 'above threshold (R0=' + hi.r0.toFixed(2) + ') a real OUTBREAK sweeps the network (attack ' + (hi.attackRate * 100).toFixed(0) + '%)');
  ok(lo.r0 < 0.95 && lo.attackRate < 0.12, 'below threshold (R0=' + lo.r0.toFixed(2) + ') the pathogen FIZZLES (attack ' + (lo.attackRate * 100).toFixed(0) + '%)');
}

// 4. THE CURVE TURNS OVER -- an epidemic is a WAVE, not a monotone climb: I is tiny, rises to an interior peak,
//    then breaks back down. (No waning, so it ends near zero.)
{
  const { series, f } = run(13, { N: 420, beta: 0.11 }, 200);
  const peak = peakOf(series), last = series[series.length - 1].nI;
  ok(peak > series[2].nI * 5 && last < peak * 0.4 && f.peakDay > 4 && f.peakDay < 190, 'the epidemic CURVE rises to an interior peak (' + peak + ' at day ' + f.peakDay.toFixed(0) + ') then breaks (ends ' + last + ')');
}

// 5. THE TURNING POINT IS HERD IMMUNITY, AND IT OVERSHOOTS -- the peak arrives when Reff falls through 1, i.e.
//    when the susceptible fraction reaches ~1/R0 (still well above zero). And the epidemic does NOT stop there:
//    infectious people keep circulating, so it OVERSHOOTS, ending with fewer susceptible than at the peak.
{
  const { series, f } = run(17, { N: 440, beta: 0.11, infDur: 6, incub: 3 }, 220);
  const finalS = series[series.length - 1].susceptible, inv = 1 / f.r0;
  ok(Math.abs(f.sAtPeak - inv) < 0.16, 'the curve PEAKS when susceptibles reach ~1/R0 (sAtPeak ' + f.sAtPeak.toFixed(2) + ' vs 1/R0 ' + inv.toFixed(2) + ')');
  ok(finalS < f.sAtPeak - 0.04, 'the epidemic OVERSHOOTS herd immunity -- it does not stop at the turning point (final susceptible ' + finalS.toFixed(2) + ' < at-peak ' + f.sAtPeak.toFixed(2) + ')');
}

// 6. HERD IMMUNITY THRESHOLD -- vaccinate ABOVE 1-1/R0 up front and the outbreak never ignites; vaccinate a bit
//    BELOW it and it still sweeps. The threshold, proven by crossing it.
{
  const probe = run(19, { N: 440, beta: 0.11 }, 5).f; const thr = 1 - 1 / probe.r0;
  const above = run(19, { N: 440, beta: 0.11 }, 150, (st, eco) => eco.vaccinate(st, Math.min(0.95, thr + 0.12), false)).f;
  const below = run(19, { N: 440, beta: 0.11 }, 150, (st, eco) => eco.vaccinate(st, Math.max(0, thr - 0.25), false)).f;
  ok(above.attackRate < 0.10, 'vaccinating ABOVE the herd threshold (' + (thr * 100).toFixed(0) + '%) stops ignition (attack ' + (above.attackRate * 100).toFixed(0) + '%)');
  ok(below.attackRate > above.attackRate + 0.15, 'vaccinating BELOW it still lets the outbreak spread (attack ' + (below.attackRate * 100).toFixed(0) + '% >> ' + (above.attackRate * 100).toFixed(0) + '%)');
}

// 7. TARGET THE SUPERSPREADERS -- the network's payoff: immunising the SAME fraction of people does far more when
//    you pick the highest-degree hubs than when you pick at random. (Fraction chosen below the random-herd level.)
{
  const frac = 0.18;
  const rand = run(23, { N: 440, beta: 0.11 }, 150, (st, eco) => eco.vaccinate(st, frac, false)).f;
  const targ = run(23, { N: 440, beta: 0.11 }, 150, (st, eco) => eco.vaccinate(st, frac, true)).f;
  ok(targ.peakI < rand.peakI * 0.6, 'TARGETING superspreaders crushes the peak far more than random immunisation at equal cost (peak ' + targ.peakI + ' vs ' + rand.peakI + ')');
  ok(targ.attackRate < rand.attackRate - 0.08, 'and it lowers the final size too (attack ' + (targ.attackRate * 100).toFixed(0) + '% vs ' + (rand.attackRate * 100).toFixed(0) + '%)');
}

// 8. FINAL SIZE grows with transmissibility -- a more infectious pathogen infects a larger share before it burns out.
{
  const a = run(29, { N: 420, beta: 0.05 }, 220).f.attackRate;
  const b = run(29, { N: 420, beta: 0.16 }, 220).f.attackRate;
  ok(b > a + 0.12, 'a more infectious pathogen reaches a larger FINAL SIZE (attack ' + (a * 100).toFixed(0) + '% -> ' + (b * 100).toFixed(0) + '%)');
}

// 9. FLATTEN THE CURVE -- distancing (a lower effective beta) lowers AND delays the peak, even at a similar final
//    size: the hospital-capacity argument, shown by running with and without it.
{
  const normal = run(31, { N: 440, beta: 0.11 }, 240);
  const distanced = run(31, { N: 440, beta: 0.11 }, 240, (st) => { st.distancing = 0.55; });
  ok(distanced.f.peakI < normal.f.peakI * 0.8, 'DISTANCING flattens the peak (peak ' + normal.f.peakI + ' -> ' + distanced.f.peakI + ')');
  ok(distanced.f.peakDay > normal.f.peakDay, 'and DELAYS it (peak day ' + normal.f.peakDay.toFixed(0) + ' -> ' + distanced.f.peakDay.toFixed(0) + ')');
}

// 10. WANING IMMUNITY -> ENDEMIC -- if immunity fades, recovered people become susceptible again and the disease
//     does not burn out but settles into a persistent, recurring presence.
{
  const permanent = run(37, { N: 440, beta: 0.12 }, 300).f;
  const waning = run(37, { N: 440, beta: 0.12, waning: 0.02 }, 300).f;
  ok(permanent.nI === 0, 'with permanent immunity the epidemic BURNS OUT (final infectious ' + permanent.nI + ')');
  ok(waning.nI > 0, 'with WANING immunity it goes ENDEMIC -- never fully clearing (final infectious ' + waning.nI + ')');
}

// 11. MORTALITY tracks the attack rate -- with an infection-fatality ratio, the dead fraction is about the fatality
//     ratio times the fraction ever infected; with none, nobody dies.
{
  const none = run(41, { N: 420, beta: 0.12, ifr: 0 }, 220).f;
  const deadly = run(41, { N: 420, beta: 0.12, ifr: 0.05 }, 220).f;
  ok(none.nD === 0, 'with no fatality ratio, no one DIES of it (dead ' + none.nD + ')');
  ok(deadly.dead > 0 && Math.abs(deadly.dead - deadly.attackRate * 0.05) < 0.02, 'DEATHS ~ fatality-ratio x attack-rate (dead ' + (deadly.dead * 100).toFixed(1) + '% ~ ' + (deadly.attackRate * 5).toFixed(1) + '%)');
}

// 12. SUPERSPREADERS EXIST -- the contact network is heavy-tailed: the most-connected person has far more contacts
//     than the average, which is the whole reason the network R0 exceeds the naive contact-rate estimate.
{
  const f = run(43, { N: 440, meanDegree: 8 }, 2).f;
  ok(f.maxDegree > f.meanDegree * 3, 'the network has SUPERSPREADERS -- max degree ' + f.maxDegree + ' >> mean ' + f.meanDegree.toFixed(1));
}

// 13. QUARANTINE -- isolating a fraction of cases so they cannot transmit reduces the outbreak, and enough of it
//     crushes it entirely. Like vaccination but acting on the INFECTIOUS instead of the susceptible.
{
  const none = run(11, { N: 440, beta: 0.11 }, 200, (st) => { st.quarantine = 0; }).f.attackRate;
  const some = run(11, { N: 440, beta: 0.11 }, 200, (st) => { st.quarantine = 0.5; }).f.attackRate;
  const lots = run(11, { N: 440, beta: 0.11 }, 200, (st) => { st.quarantine = 0.8; }).f.attackRate;
  ok(some < none - 0.1, 'QUARANTINE (isolating cases) reduces the outbreak (attack ' + (none * 100).toFixed(0) + '% -> ' + (some * 100).toFixed(0) + '%)');
  ok(lots < 0.1, 'enough quarantine CRUSHES the outbreak before it spreads (attack ' + (lots * 100).toFixed(0) + '%)');
}

// 14. SUPERSPREADING / OVERDISPERSION -- transmission is not evenly shared: most infected people infect nobody
//     or one, while a few cause many. Measured over the causal tree the sim now records (secondary cases per
//     person): the top fifth of spreaders account for well over half of all infections. The 80/20 of outbreaks.
{
  const { st } = run(11, { N: 440, beta: 0.11 }, 220);
  const sec = []; for (let i = 0; i < st.N; i++) if (st.state[i] !== 0) sec.push(st.secondary[i]);   // ever-infected (permanent immunity: not still susceptible)
  sec.sort((a, b) => b - a);
  const total = sec.reduce((s, v) => s + v, 0), zero = sec.filter(v => v === 0).length;
  let topSum = 0, top = Math.ceil(sec.length * 0.2); for (let i = 0; i < top; i++) topSum += sec[i];
  ok(total > 0 && topSum / total > 0.5, 'SUPERSPREADING: the top 20% of spreaders cause > half of all infections (they caused ' + (100 * topSum / total).toFixed(0) + '%)');
  ok(zero / sec.length > 0.4, 'most infected people infect NOBODY (overdispersion: ' + (100 * zero / sec.length).toFixed(0) + '% caused zero further cases)');
}

console.log(`\n${pass} passed, ${fail} failed`);
process.exit(fail ? 1 : 0);