在瀏覽器中求解流感與免疫模型：常微分方程、RK4 與三個數值陷阱

When a virus or bacterium gets into the body, how does the immune system fight it? I made a page that computes the course of an infection with a model made of ordinary differential equations and shows it as graphs and a microscopic cross-section. Canvas 2D, no libraries.

https://sishiya.github.io/manabi-lab/immune-battle/?lang=en

You can pick influenza, COVID-19, norovirus, Staphylococcus aureus, pneumococcus or E. coli, and try medicines and vaccines. This post covers how the influenza model is built and where the numerics tripped me up.

The numbers in this app are rough values, tuned to literature model shapes so that well-known courses (days of fever and so on) come out. They cannot be used for diagnosis or treatment decisions.

Base: a target-cell-limited model

Influenza builds on the "target-cell-limited model with an eclipse phase" of Baccam et al. (2006). Cells of the nose and throat lining are split into five fractions (summing to 1):

The total is N0 = 4×10⁸ cells (Baccam 2006, upper airway). Virus V is kept as a count.

Main values:

• eclipse kE = 4/day (about 6 hours)
• infected-cell lifespan del = 1.5/day
• virus released per cell Y = 2000 (particles; far fewer are infectious)

On top of that I added interferon, NK cells, macrophages, killer T cells, B cells and antibodies.

A small dose does not infect: saturating clearance

In the original model, a single virus will always grow under some conditions. In reality, a small amount of virus does not cause infection.

So clearance by macrophages and mucus was made saturating:

clearance rate = Cm · M · V / (1 + V/Vs)     (Cm = 150, Vs = 2000)

While V is small, each virus is cleared quickly and clearance beats growth. When V gets large, clearance levels off and cannot keep up.

As a result, in a healthy adult illness starts above roughly 100,000 viruses (an estimate). Human nasal challenge studies use around 10⁵–10⁷ TCID50, so I think the order of magnitude is not off.

For a healthy adult with 100,000 viruses:

• fever on day 2, up to 38.8 °C
• peak on days 3–4, fever gone on day 6
• virus gone on day 8

Gotcha 1: adding antibodies blew up to 10¹³

I was using ordinary 4th-order Runge–Kutta (RK4) with a fixed step. When antibody neutralization was fast (around 1000 per day), the integration diverged and counts jumped to 10¹³.

When something decays quickly, a too-large step overshoots (a stiff equation). So inside each step, the solver now subdivides automatically so that the fastest decay rate × step < 0.8:

let left = h, n = 0;
while (left > 1e-12 && n++ < 400) {
  path.f(y, P, k1);
  let rate = 0;
  for (const v of vs) if (k1[v] < 0 && y[v] > 1e-30) rate = Math.max(rate, -k1[v] / y[v]);
  const s = Math.min(left, 0.8 / Math.max(rate, 1e-9));
  // ...an ordinary RK4 step of size s (intermediate values clamped at 0)
  left -= s;
}

The decay rate -k1[v] / y[v] is "what fraction is lost per day", so its inverse is a rough timescale for that quantity.

Gotcha 2: 10⁻⁵⁰ of a cell growing back

While testing norovirus, the virus started growing again months after it should have been gone.

It turned out infected cells were left at absurdly small amounts like 10⁻⁵⁰. In the equations, even that can slowly grow.

So when almost nothing is left, it is set to 0:

// less than one particle / cell left = none (otherwise 10⁻⁵⁰ of a cell can grow back in the equations)
if (path.kind === 'virus') { if (y.V < 1e-3 && (y.I + y.E) * path.N0 < 1e-3) { y.V = 0; y.I = 0; y.E = 0; } }

A model that uses continuous quantities for counts needs "counts are whole numbers" put in somewhere.

Gotcha 3: inflammation without bacteria

In the S. aureus model, inflammation never ended even after the bacteria were gone.

The inflammation signal also came from "tissue damage", and damage also came from neutrophils' collateral harm. Signal → neutrophils → damage → signal made a loop. The fix: damage only produces a signal while bacteria are present.

Medicine: make "time" part of the state

How well a drug works depends on time since the first dose. So time t itself is one of the state variables (dt/dt = 1). Then any equation can compute the drug effect from y.t.

Oseltamivir blocks 85% of virus release from cells (5 days, half-life about 8 hours). Starting about 5 hours after fever begins shortens the fever by about a day, roughly the shortening reported in clinical trials.

How I check it

Every time, all six pathogens × four bodies × every memory option run to the end, checking for errors:

for (const pk of PATHOGEN_ORDER) for (const body of BODY_ORDER)
  for (const mem of Object.keys(PATHOGENS[pk].memOptions)) {
    const B = BODIES[body];
    __ib.set({ pk, body, mem, inn: B.inn, adp: B.adp, dose: PATHOGENS[pk].doseDefault, speed: 2, playing: true });
    __ib.frame(1200);
  }
__ib.err

Main sources

• P. Baccam, C. Beauchemin, C. A. Macken, F. G. Hayden, A. S. Perelson (2006) "Kinetics of influenza A virus infection in humans", Journal of Virology 80, 7590
• P. F. M. Teunis et al. (2008) Journal of Medical Virology (how many noroviruses can infect)
• S. D. Elek (1957) Annals of the New York Academy of Sciences (how many S. aureus cause skin infection)
• Immunology textbooks (such as Janeway's Immunobiology)

This page is part of Manabi Lab, a series of small learning apps that run in the browser (English and Japanese). The code and this post were written together with AI (Claude).

https://sishiya.github.io/manabi-lab/?lang=en