epidemic.js
188 lines · in natura-src-1.0.0.zip as natura/epidemic.js · a module · its entry in the module reference · its test: epidemic-test.js
// epidemic.js -- Natura/EXISTA EPIDEMIC (the SOCIETY domain's involuntary verb: CATCHING). A contagion moving
// through a population, run as real network epidemiology -- not a well-mixed scalar SIR but a POPULATION on a
// CONTACT NETWORK, where a pathogen passes only along the edges between people, and a few highly-connected
// SUPERSPREADERS carry the outbreak far past what the average contact-rate would predict. Each person is
// SUSCEPTIBLE, then (if a contact infects them) EXPOSED and incubating, then INFECTIOUS, then RECOVERED-and-
// immune (or DEAD). Out of nothing but that per-contact chance comes the whole shape of a season: the epidemic
// CURVE that rises, peaks, and breaks; the peak that arrives BEFORE the susceptibles run out (the turning point
// is herd immunity, and the epidemic OVERSHOOTS it); the HERD-IMMUNITY THRESHOLD past which a vaccinated
// population never ignites; the way immunising the few superspreaders stops far more than immunising the many
// at random; and -- with waning immunity -- the drop into recurrent, ENDEMIC waves. This is not mood-contagion
// (see society.js, where fear spreads); this is a real disease with immunity and mortality, run as physics.
//
// A Natura system (init/step/fields/events), deterministic + node-testable. Randomness only from init's seeded
// rng (the network + initial seed) and ctx.rng (per-contact transmission, recovery, mortality). UMD.
(function (root, factory) {
if (typeof module !== 'undefined' && module.exports) module.exports = factory();
else root.Epidemic = factory();
})(typeof self !== 'undefined' ? self : this, function () {
'use strict';
const clamp = (v, a, b) => v < a ? a : v > b ? b : v;
const S = 0, E = 1, I = 2, R = 3, D = 4; // susceptible / exposed / infectious / recovered / dead
function makeEpidemic(config) {
config = config || {};
const N = config.N || 420; // people
const meanDegree = config.meanDegree != null ? config.meanDegree : 8; // average contacts per person
const hubSigma = config.hubSigma != null ? config.hubSigma : 0.6; // degree heterogeneity (lognormal spread of sociability): real superspreaders, not one mega-hub
const beta = config.beta != null ? config.beta : 0.10; // per-contact chance of transmission per infectious day
const incub = config.incub != null ? config.incub : 3; // incubation (E->I) in days; 0 = SIR (no latent stage)
const infDur = config.infDur != null ? config.infDur : 6; // infectious period (I->R/D) in days
const ifr = config.ifr != null ? config.ifr : 0.012; // infection fatality ratio (I -> D instead of R)
const waning = config.waning != null ? config.waning : 0; // immunity loss rate (R->S per day); 0 = permanent immunity, >0 = endemic
const seed0 = config.seed0 != null ? config.seed0 : 3; // initial infectious people (patient zero + a few)
// ---- network: a heavy-tailed contact graph so a few SUPERSPREADERS emerge (weighted edge sampling) ----
function buildNetwork(rng, x, y) {
const weight = new Float64Array(N);
let wsum = 0; const cum = new Float64Array(N);
for (let i = 0; i < N; i++) { // lognormal sociability: moderate heterogeneity -> superspreaders WITHOUT a single pathological mega-hub
const g = Math.sqrt(-2 * Math.log(1e-12 + rng())) * Math.cos(6.2831853 * rng());
const w = Math.exp(hubSigma * g); weight[i] = w; wsum += w; cum[i] = wsum;
}
const pick = () => { const r = rng() * wsum; let lo = 0, hi = N - 1; while (lo < hi) { const m = (lo + hi) >> 1; if (cum[m] < r) lo = m + 1; else hi = m; } return lo; };
const set = Array.from({ length: N }, () => new Set());
const Etot = Math.round(meanDegree * N / 2);
let made = 0, tries = 0, cap = Etot * 40;
while (made < Etot && tries < cap) {
tries++; const a = pick(), b = pick();
if (a === b || set[a].has(b)) continue;
set[a].add(b); set[b].add(a); made++;
}
const adj = set.map(s => Int32Array.from(s));
const deg = Int32Array.from(adj.map(a => a.length));
return { adj, deg };
}
return {
name: config.name || 'epidemic',
domain: 'society',
scale: 'population',
params: config,
init(rng) {
// fixed scatter positions (a disc), used only by the viz
const x = new Float32Array(N), y = new Float32Array(N);
for (let i = 0; i < N; i++) { const ang = rng() * 6.2831853, rad = Math.sqrt(rng()); x[i] = Math.cos(ang) * rad; y[i] = Math.sin(ang) * rad; }
const { adj, deg } = buildNetwork(rng, x, y);
// network-corrected basic reproduction number depends on the EXCESS degree <k^2-k>/<k>, not the mean --
// this is exactly why superspreaders make an outbreak bigger than the average contact-rate suggests.
let sk = 0, sk2 = 0; for (let i = 0; i < N; i++) { sk += deg[i]; sk2 += deg[i] * deg[i]; }
const meanK = sk / N, meanK2 = sk2 / N, excessK = meanK > 0 ? (meanK2 - meanK) / meanK : 0;
const state = new Uint8Array(N); // all S
const eT = new Float32Array(N), iT = new Float32Array(N), iso = new Uint8Array(N);
// seed the outbreak: infect the seed0 people of MEDIAN-ish connectivity (not the top hub, so ignition is
// a fair test) -- pick by a fixed stride so it is deterministic and spread across the population.
let planted = 0; for (let i = 0; i < N && planted < seed0; i += Math.max(1, Math.floor(N / seed0))) { state[i] = I; iT[i] = infDur; planted++; }
return {
N, adj, deg, x, y, state, eT, iT, iso,
meanK, meanK2, excessK,
distancing: 1, // global transmissibility multiplier (a viz/scenario lever; 1 = normal, <1 = distancing)
quarantine: 0, // fraction of newly-infectious who self-isolate (do not transmit)
nS: N - planted, nE: 0, nI: planted, nR: 0, nD: 0,
everInfected: planted, peakI: planted, peakDay: 0, sAtPeak: 1,
reff: 0, r0: 0, day: 0, prevI: planted,
parent: (() => { const p = new Int32Array(N); p.fill(-1); return p; })(), // who infected each person (-1 = a seed / never infected)
gen: new Int16Array(N), // transmission generation: seeds 0, their cases 1, ...
secondary: new Int32Array(N), // how many people each person went on to infect (the superspreading distribution)
_newExp: new Uint8Array(N), _newExpBy: new Int32Array(N), justInfected: [], evs: []
};
},
step(st, ctx, dt) {
const rng = ctx && ctx.rng ? ctx.rng : Math.random; // seeded world rng
const dd = clamp((dt || 86400) / 86400, 1e-3, 30); // step length in days
const adj = st.adj, state = st.state, eT = st.eT, iT = st.iT, iso = st.iso, newExp = st._newExp;
st.evs.length = 0;
newExp.fill(0);
const pInf = clamp(beta * st.distancing * dd, 0, 1); // per-contact transmission chance this step
// --- 1. TRANSMISSION along the network: every infectious (non-isolated) person rolls against each of its
// susceptible contacts. Freshly-exposed people are staged in newExp and do NOT transmit this step
// (clean generation separation), so the outbreak advances one contact-hop at a time.
const newExpBy = st._newExpBy;
for (let i = 0; i < st.N; i++) {
if (state[i] !== I || iso[i]) continue;
const nb = adj[i];
for (let e = 0; e < nb.length; e++) { const j = nb[e]; if (state[j] === S && !newExp[j] && rng() < pInf) { newExp[j] = 1; newExpBy[j] = i; } } // remember WHO infected j (first contact wins)
}
// --- 2. CLINICAL PROGRESSION (using the state BEFORE this step's new exposures): E->I after incubation,
// I->R (or ->D by the fatality ratio) after the infectious period, R->S if immunity wanes. Periods
// carry per-person spread so the population does not transition in lockstep (real waves are smooth).
for (let i = 0; i < st.N; i++) {
const s = state[i];
if (s === E) { eT[i] -= dd; if (eT[i] <= 0) { state[i] = I; iT[i] = infDur * (0.6 + 0.8 * rng()); iso[i] = (st.quarantine > 0 && rng() < st.quarantine) ? 1 : 0; } }
else if (s === I) { iT[i] -= dd; if (iT[i] <= 0) { state[i] = (rng() < ifr) ? D : R; iso[i] = 0; } }
else if (s === R) { if (waning > 0 && rng() < waning * dd) state[i] = S; }
}
// --- 3. apply the new exposures: S -> E (or straight to I if there is no latent stage) ---
st.justInfected.length = 0;
const parent = st.parent, gen = st.gen, secondary = st.secondary;
for (let i = 0; i < st.N; i++) {
if (newExp[i] && state[i] === S) {
st.everInfected++; st.justInfected.push(i);
const by = newExpBy[i]; parent[i] = by; gen[i] = (gen[by] + 1) | 0; secondary[by]++; // record the transmission link, the generation, and credit the infector a secondary case
if (incub > 0) { state[i] = E; eT[i] = incub * (0.6 + 0.8 * rng()); }
else { state[i] = I; iT[i] = infDur * (0.6 + 0.8 * rng()); iso[i] = (st.quarantine > 0 && rng() < st.quarantine) ? 1 : 0; }
}
}
// --- 4. tally, and the LIVE reproduction number. R0 is the network-corrected transmissibility x excess
// degree; Reff = R0 x (susceptible fraction) -- it starts above 1 (outbreak grows) and falls
// through exactly 1 AT THE PEAK, which is why the curve turns over before everyone is infected.
let nS = 0, nE = 0, nI = 0, nR = 0, nD = 0;
for (let i = 0; i < st.N; i++) { const s = state[i]; if (s === S) nS++; else if (s === E) nE++; else if (s === I) nI++; else if (s === R) nR++; else nD++; }
st.nS = nS; st.nE = nE; st.nI = nI; st.nR = nR; st.nD = nD;
const T = 1 - Math.pow(1 - clamp(beta * st.distancing, 0, 1), infDur); // chance of transmission across a contact over the whole infectious period
st.r0 = T * st.excessK;
st.reff = st.r0 * (nS / st.N);
st.day += dd;
if (nI > st.peakI) { st.peakI = nI; st.peakDay = st.day; st.sAtPeak = nS / st.N; }
// --- events ---
if (st.prevI < st.N * 0.02 && nI >= st.N * 0.02) st.evs.push({ kind: 'outbreak', nI });
if (st.prevI > st.peakI * 0.98 && nI < st.prevI && st.day > 2) st.evs.push({ kind: 'peak', nI, day: st.day });
if (st.prevI > 0 && nI === 0) st.evs.push({ kind: 'faded', attack: st.everInfected / st.N });
st.prevI = nI;
return st;
},
fields(st) {
return {
N: st.N, x: st.x, y: st.y, deg: st.deg, adj: st.adj, state: st.state, justInfected: st.justInfected,
parent: st.parent, gen: st.gen, secondary: st.secondary,
nS: st.nS, nE: st.nE, nI: st.nI, nR: st.nR, nD: st.nD,
susceptible: st.nS / st.N, infectious: st.nI / st.N,
everInfected: st.everInfected, attackRate: st.everInfected / st.N, dead: st.nD / st.N,
r0: st.r0, reff: st.reff, herdThreshold: st.r0 > 1 ? 1 - 1 / st.r0 : 0,
peakI: st.peakI, peakDay: st.peakDay, sAtPeak: st.sAtPeak,
meanDegree: st.meanK, maxDegree: (() => { let m = 0; for (let i = 0; i < st.N; i++) if (st.deg[i] > m) m = st.deg[i]; return m; })(),
day: st.day
};
},
events(st) { return st.evs; },
// ---- scenario helpers (a reader/viz/test may reach in for a controlled experiment) ----
// vaccinate a fraction of the CURRENTLY-susceptible population, either at RANDOM or TARGETED at the
// highest-degree superspreaders. Deterministic (no rng needed) so experiments reproduce. Returns the count.
vaccinate(st, frac, targeted) {
const k = Math.floor(clamp(frac, 0, 1) * st.N);
let order;
if (targeted) order = Array.from({ length: st.N }, (_, i) => i).sort((a, b) => st.deg[b] - st.deg[a]);
else order = Array.from({ length: st.N }, (_, i) => i).sort((a, b) => (((a * 2654435761) >>> 0) % 100003) - (((b * 2654435761) >>> 0) % 100003));
let done = 0;
for (const i of order) { if (done >= k) break; if (st.state[i] === S) { st.state[i] = R; done++; } }
st.nS = 0; for (let i = 0; i < st.N; i++) if (st.state[i] === S) st.nS++; // keep the tally honest for an immediate read
st.nR += done;
return done;
},
_S: S, _E: E, _I: I, _R: R, _D: D
};
}
return { makeEpidemic };
});