Model specification

fleet is not yet validated for real use; see the status note.

This article specifies the daily update that fleet applies. §1 states the design constraint the rest follows from, §2 what the state space does and does not carry, sections A to H the P. falciparum equations as implemented, and §V the P. vivax model. The implementation is inst/odin/malaria_daily.R, with its inputs built in R/build_inputs.R and R/interventions.R. All references to malariasimulation are to v3.0.0.

Getting started introduces the package with runnable code, vignette("parameters") covers each malariasimulation argument, and vignette("using") covers what to do differently because the model is a mean field.

1. Scope and design

fleet is a deterministic mean-field counterpart to the malariasimulation individual-based model (IBM) of Plasmodium falciparum and P. vivax, advanced one day at a time on the IBM’s own clock.

Mean field means that the model tracks the distribution of the population over states, not individual people, so the state cannot carry per-individual attributes. A quantity that depends on two attributes of the same person, such as owning a net and being vaccinated, is unavailable by construction. Most of §2.2 and §2.3 follows from this.

fleet follows the IBM’s code rather than its published equations. malariasimulation resolves a day as a sequence of processes that all read the state at the start of the day and queue their changes to its end: immunity decays, the day’s bites are drawn and deduplicated, infection and each state’s progression are resolved as one competing draw per person, the mosquito model is stepped across the day, and the dead are replaced by newborns. fleet applies the same steps in the same order to the distribution of the population over states. The IBM’s per-day arithmetic is therefore reproduced by construction rather than approximated: whole-day sojourns, integer refractory windows, the delay lines behind the EIR, human infectivity and mosquito incubation, and a boost overwriting the day’s immunity decay.

The IBM records for each person their net, their spray status, their vaccination history, their drug and the date they took it, and the date each of their four immunity functions last boosted. A mean-field model cannot afford a state dimension for each. fleet keeps two human structuring dimensions, age and biting heterogeneity (the vivax model adds a third, §V.1), and represents everything else as a time-varying coefficient in the daily update, a closure that carries its average effect, or a discrete state jump between days. §2.1 to §2.3 set out which is which.

Abbreviations and terms

Term	Meaning
EIR	Entomological inoculation rate, infectious bites per person per year. fleet reports it per adult per year
FOI, FOIM	Force of infection on humans; force of infection on mosquitoes
EIP	Extrinsic incubation period, the delay in the mosquito between an infectious blood meal and becoming infectious
LM, PCR	Light microscopy; polymerase chain reaction. The two prevalence diagnostics reported
ETF, SPC	Early treatment failure; slow parasite clearance. The two antimalarial-resistance arms malariasimulation implements
PEV, TBV	Pre-erythrocytic vaccine (RTS,S, R21); transmission-blocking vaccine
TRA, TBA	Transmission-reducing activity (an antibody-level quantity); transmission-blocking activity (the resulting per-state reduction in infectivity)
MDA, SMC, PMC	Mass drug administration; seasonal malaria chemoprevention; perennial malaria chemoprevention
EPI	The Expanded Programme on Immunization routine schedule. Unrelated to set_epi_outputs(), where “epi” means epidemiological
IRS	Indoor residual spraying
Stage chain	A linear chain of $k$ sequential stages, each left with probability $p$ a day, so each lasts a geometric number of days. Its total has mean $k/p$ and variance $k(1-p)/p^2$
Delay line	The last $n$ days of a quantity, one value per day, so a lag of $n$ days is read off directly, as the IBM’s LaggedValue and incubation queue hold it
Gauss–Hermite quadrature	A deterministic replacement for Monte Carlo sampling: a small fixed set of representative values and weights that reproduce the distribution’s moments. A tensor rule is the product grid over several such axes
Knots	The times at which an interpolated series may change value
Seed	The initial state vector, never a random-number seed. fleet has no random component

2. State space

2.1 What has a state dimension

Every human compartment is a two-dimensional array over age group $i$ and biting-heterogeneity node $j$; every mosquito compartment is indexed by species $s$. Human compartments are per-capita densities (fractions of the total human population), so run time does not depend on human_population.

Dimension	Symbol	Default size	Set by	Indexes
Age group	$i = 1,\dots,n_a$	118	age_lower = / default_age_lower()	The graded age grid described below. Individuals age by a linear-chain (McKendrick) flow $r_i$ between adjacent groups
Biting heterogeneity	$j = 1,\dots,n_z$	5	parameters$n_heterogeneity_groups, enable_heterogeneity	Gauss–Hermite nodes $\zeta_j$ and weights $w_j$ of the log-normal relative-biting-rate distribution with variance $\sigma^2$ (sigma_squared); $\zeta_j = \exp(z_j\sigma - \sigma^2/2)$. Fixed at birth and never mixed: a person stays in their stratum for life, as in the IBM
Vector species	$s = 1,\dots,n_v$	1	set_species()	Independent aquatic and adult mosquito populations, each with its own biting rate, mortality, carrying capacity and vector-control response, coupled only through the shared human population
EIR delay line	$1,\dots,\lfloor\tau_E\rfloor+1$ days	13	de	The last days of the per-species EIR $a_s I_{M,s}$, from which the lag $\tau_E$ is read, interpolating between two days when it is fractional (§B.2)
Infectivity delay line	$1,\dots,\lfloor\tau_l\rfloor+1$ days	13	delay_gam	The last days of the population’s onward infectivity, from which the gametocyte lag $\tau_l$ is read (§C). The default 12.5 d is read halfway between the values saved 12 and 13 days ago
Incubation delay line	$1,\dots,\lceil\tau_M\rceil-1$ days	9	dem	The last days of the per-species flux into incubation $S_M\Lambda^M$: the IBM’s queue of $\lceil\tau_M\rceil$ days releases the flux that entered $\lceil\tau_M\rceil - 1$ days ago (§D)

The default age grid is fine in childhood and coarse in adulthood. Group edges are pinned at the conventional reporting ages (0, 1, 2, 3, 5, 7, 10, 15, 20, 30, 40 and 60 y), so output bands fall on group boundaries, and at 21 y, so that whole groups make up the year from which maternal immunity’s mothers are drawn (§B.3). Between two anchors the groups are of equal width. The budget of groups is shared between the anchor intervals in proportion to their width in log(1 + age), because immunity rises with log age rather than with age, weighted 1.25 below 5 y and 0.5 from 21 y to move groups to where the grid’s error arises. An absorbing top group runs from 80 y.

Human compartments are all $[n_a \times n_z]$. S_count counts the uninfected, including the drug-protected, as the IBM counts them; Ph_count gives how many of them a drug is protecting.

Compartment	odin	Output column	Meaning
$S$	S	S_count (with the $P$ chains)	Susceptible
$D$	D	D_count	Clinical disease, untreated
$T$	Tr	Tr_count (with $T_s$)	Successfully treated, standard clearance
$T_s$	Tr_slow	in Tr_count	Successfully treated, slow parasite clearance. A separate compartment, so the two clearance speeds form a mixture of two sojourns rather than one sojourn at the blended mean (§B.4)
$A$	A	A_count	Asymptomatic patent infection
$U$	U	U_count	Sub-patent infection
$T_c$, $T_{c,s}$	Tr_c, Tr_cs	in Tr_count	The chemoprevention treated phase: the clinical and LM-detectable asymptomatic cases a chemoprevention round treats, which stay detectable and infectious for the Tr stay (fast and slow clearance). Each has a stock, J_c or J_cs, of the infectivity it carries (§E)
$P_{1..k_P}$	Ph	Ph_count, and in S_count	Post-treatment prophylaxis: a stage chain of $k_P$ stages, so that the whole $T + P$ sojourn matches the mean and variance of the IBM’s Weibull protection rather than lasting a single geometric time (§B.4). $k_P$ is 9 for SP-AQ, 10 for DHA-PQP, 1 for AL (whose 10-day protection is less variable than $T$ itself), and 1 with no clinical drugs
$P_{c,1..k_{P_c}}$	Ph_c	as $P$	Chemoprevention prophylaxis (MDA, SMC, PMC). A separate chain, because its drug, and so its mean duration and shape, can differ from the clinical first line: 10 stages for SP-AQ, 9 for DHA-PQP
$P_{ct,1..k_{P_{ct}}}$	Ph_ct	as $P$	The protection left after the chemoprevention treated phase, sized as the post-treatment chain is (§E)
$I_B$	IB	ib_mean	Pre-erythrocytic (anti-infection) immunity
$I_{CA}$	ICA	ica_mean	Acquired clinical immunity
$I_D$	ID	id_mean	Anti-parasite (detection) immunity
$I_{VA}$	IVA	iva_mean	Acquired severe-disease immunity

The four immunity arrays hold the mean immunity of the individuals in each (age, heterogeneity) cell, not an extensive density, so their ageing term is a difference $\alpha_{ij}(I_{i-1,j} - I_{ij})$, with $\alpha_{ij}$ the day’s inflow per head (§B.5), rather than a flux.

Mosquito compartments are all $[n_v]$:

Compartment	Meaning
$E$, $L$, $P_L$	Aquatic early larval, late larval and pupal (odin names ME, ML, MP). Pupae emerge at $P_L/d_{P_L}$, and the half that are female enter $S_M$. Males are not modelled, so the factor of $\tfrac12$ on emergence is a sex ratio, not a loss
$S_M$	Susceptible adult females
$E_M$	Incubating (exposed) adults (odin Em)
$I_M$	Infectious adults

With the defaults ($n_a = 118$, $n_z = 5$, $n_v = 1$) and no drugs, so that every chain has one stage, there are 10,030 falciparum human states, $(14 + k_P + k_{P_c} + k_{P_{ct}})\,n_a n_z$, and 6 mosquito states per species. To these are added the three delay lines (35 days of history) and the day’s recorded outputs, which are carried as state because a daily update has nowhere else to keep them. Each extra prophylaxis stage adds $n_a n_z = 590$ states. The vivax block (§V) multiplies its compartments by the eight levels of its hypnozoite ladder, which hold the eleven batch counts (§V.1).

2.2 What is captured without a dimension

These mechanisms are per-individual attributes in the IBM and time-varying (in some cases age-varying) coefficients in fleet. Each row is a dimension that was not added.

Mechanism	IBM representation	fleet representation	What is lost
Clinical treatment coverage	Per-person Bernoulli draw at a clinical episode	Scalar $f_t(t)$, step-interpolated over the union of all drugs’ timesteps, splitting the clinical inflow between $T/T_s$ and $D$	Nothing at the flow level
Which drug a person received	Drug index stored per treated individual, determining their efficacy, infectivity and prophylaxis duration	Three scalar series $\text{eff}(t)$, $c_T(t)$, $r_P(t)$, each a coverage-share-weighted blend across the active clinical drugs, step-interpolated over the coverage change times	A first-line switch is modelled (the blend moves in time); a mixed first line is represented by its mean, not its mixture
Post-treatment prophylaxis clock	Weibull survival $W(t - t_{\text{drug}})$ multiplying each treated person’s infection probability, running from the dose	A stage chain entered after $T$, its mean and stage count moment-matched to $W$ (§B.4)	The chain reproduces the integrated protection exactly and, for SP-AQ and DHA-PQP, the survival curve to within a few points (§B.4). It protects fully while it lasts, where the IBM’s protection is partial and fades, and AL’s single stage spreads its sharp 10-day protection
Early treatment failure	Per-person draw returning a treated case to clinical	Folded into the effective treated fraction $f_t^{\text{eff}}(t)$ (§B.4)	Nothing
Bed net ownership, net age, retention, repeat rounds	Per-person net-receipt timestep; efficacy decays from that date; nets lost stochastically	Two scalars per species, $a_s(t)$ and $\mu_s(t)$, from the mean-field mixture over all past distributions, each weighted by most-recent-receipt probability $\times$ retention survival and decaying with its own $d_{n0}$ and $r_n$	Protection is population-averaged rather than correlated within individuals across bites, so nets and IRS tend to under-suppress transmission relative to the IBM
IRS spray status and spray age	Per-person house-spray timestep	Folded into the same $a_s(t)$, $\mu_s(t)$: a mixture over all past rounds weighted by most-recent-spray probability, with no retention factor (sprayed protection never expires in the IBM)	As for bed nets
Vaccination status, dose count, booster stratum, antibody titre	Per-person dose history and last-vaccination date. Antibody parameters are not stored: the IBM re-draws all four from their profile distributions every timestep, for every vaccinated person	One $[n_a \times n_t]$ multiplier $v_i(t)$ on the infection hazard. A campaign’s protection stays with the cohort it vaccinated, its age band moving up with the time since; repeated rounds combine as the most recent dose; the vaccinated are partitioned analytically by the most recent booster each has reached; each stratum’s efficacy is the average of the Hill efficacy over the IBM’s own 4-D antibody distribution (create_pev_profile())	The antibody average is exact, since the IBM’s draws are independent. The dose history is lost: vaccination status is uncorrelated with anything else (nets, treatment) and between rounds, min_wait re-vaccination exclusion is not applied (set_mass_pev()), and seasonal boosters become a fixed days-since-primary schedule (with a warning)
TBV antibody titre	Per-person, mapped to transmission-reducing activity	Four $[n_a \times n_t]$ infectivity multipliers, one per infectious state ($U$, $A$, $D$, $T$), because the IBM’s TRA-to-TBA map is state-specific (§C)	Only the mean-field averaging
Immunity boosting refractory clock	Per-person last_boosted timestep; a boost fires only if $\ge u$ days have passed	For bites ($I_B$), the renewal probability on the day’s bite probability, with the wait rounded as the IBM rounds it; for infections ($I_{CA}$, $I_D$, $I_{VA}$), a quasi-steady closure for where the recently boosted have got to (§B.5)	The renewal form is exact for the mean inter-boost interval. The closure omits the refractory who are back in $S$ within the window, of order $10^{-3}$ of the refractory days in $A$
Age-specific mortality, demographic transition	Per-person death hazard by age band and year	$\mu_i(t)$, a constant-interpolated $[n_a \times n_t]$ series; its $t = 0$ row also sets the seed’s age structure (§F)	Ages above the top deathrate_agegroups take the top rate, where the IBM removes them; a warning is given
Maternal immunity	Sampled from a mother aged 20 to 21 drawn at birth	Algebraic, from the population-weighted mean acquired immunity of those aged 20 to 21 (§B.3), with no state variable	Nothing beyond the mean-field average over mothers
Seasonality	Daily rainfall drives larval carrying capacity	$K_s(t)$ on a daily grid from the same truncated Fourier series, held at the day’s start value across the day as the IBM holds it (§D)	Nothing (identical functional form)
Repeated bites in one day	Bitten individuals are collected in a set, so a person bitten several times counts once and can be infected at most once a day	The day’s infection probability saturates, $1 - e^{-\varepsilon}$ bitten times $b$, rather than growing with every bite (§B.2)	Nothing: it is the IBM’s own daily probability

Mass drug administration (MDA, SMC, PMC) is also carried without a dimension, but as an event rather than a coefficient: fleet stops at the end of every day that carries a scheduled round, applies a state jump (§E), and resumes.

2.3 What is not captured

Feature	Reason	Behaviour
MDA, SMC or PMC under P. vivax	malariasimulation itself fails on the combination (update_mass_drug_admin() reads the falciparum-only id)	Error at input
Individual mosquitoes (individual_mosquitoes = TRUE)	fleet is always compartmental	Ignored
Intervention correlation (get_correlation_parameters(), run_simulation(correlations =))	Individual-level; the mean field assumes independence between interventions	Warned and ignored
Per-individual stochastic variation	Deterministic model	By construction
Partner-drug resistance, late clinical or parasitological failure, reinfection during prophylaxis	malariasimulation itself rejects non-zero values for these in set_antimalarial_resistance()	Cannot reach the model

A. Notation

Abbreviations and terms are listed in §1.

Symbol	odin name	Meaning
$r_i$	r_age	Fraction of age group $i$ ageing into the next each day, exponentially fitted (§B.1); $r_{n_a} = 0$ (absorbing)
$\mu_i(t)$	mu_age	Age- and time-specific death rate
$\rho_i(t) = r_i + \mu_i(t)$	re	Fraction of age group $i$ leaving it each day
$\psi_i$	psi	Relative biting rate by age, $1 - \rho\,e^{-a_i/a_0}$
$\zeta_j$, $w_j$	zeta, het_wt	Heterogeneity node and quadrature weight
$N_{ij}$	Npop	Total population density in cell $(i,j)$
$\varepsilon_{ij}$	EPS	Expected infectious bites per person per day
$p_{ij}$	p_inf	Probability a person in $S$, $A$ or $U$ is infected today
$h^b_{ij} = -\log(1 - b_{ij}v_i)$	hb	The infection hazard of someone bitten today, for the competing draw
$b_{ij}$	b	Probability an infectious bite establishes infection
$\phi_{ij}$	phi	Probability an infection is clinical
$\theta_{ij}$	theta	Probability an infection is severe
$q_{ij}$	q	Probability an asymptomatic infection is LM-detectable
$c_A, c_D, c_U, c_T$	cA, cD, cU, cT	Onward infectivity to mosquitoes, by state
$f_t(t)$	ft	Clinical treatment coverage
$v_i(t)$	pev_factor	PEV hazard multiplier
$a_s(t)$, $\mu_s(t)$	a_spp, mum	Human blood-meal rate and adult mortality, per species
$K_s(t)$	Kcap	Larval carrying capacity
$\tau_E, \tau_l, \tau_M$	de, delay_gam, dem	EIR lag, gametocyte lag, extrinsic incubation period

B. Human model

Each day takes the state $Y(t)$ at the start of day $t+1$ (the IBM’s timestep $t+1$) to $Y(t+1)$ at its end. Every quantity on the right-hand side of an update is evaluated at the start of the day, as every IBM process reads the state before any of the day’s changes land.

B.1 Demography

Ageing is a linear chain over the age grid. Deaths are recycled as births into age group 1, distributed across heterogeneity nodes by the quadrature weights, so the total population is conserved exactly:

$$N_{ij} = S_{ij}+D_{ij}+A_{ij}+U_{ij}+T_{ij}+T_{s,ij}+T_{c,ij}+T_{cs,ij}+P_{ij}+P_{c,ij}+P_{ct,ij}, \qquad B = \sum_{i,j}\mu_i(t)\,N_{ij}.$$

Every human compartment $Y$ carries the same structural terms

$$Y_{ij}(t+1) \;\supset\; Y_{ij} + r_{i-1}Y_{i-1,j} \;-\; \rho_i(t)\,Y_{ij},$$

where at $i = 1$ the ageing inflow is replaced by $B\,w_j$ for $Y = S$ and by zero for every other compartment.

The ageing rate is exponentially fitted: $r_i = \mu_i(0) / (e^{\mu_i(0) h_i} - 1)$ for band width $h_i$ in days, computed as mu/expm1(mu*h) to avoid cancellation when $\mu h$ is small. This makes the stationary ratio between adjacent groups exactly $e^{-\mu_i h_i}$, where the obvious $r_i = 1/h_i$ gives $1/(1 + \mu_i h_i)$ and so always decays too slowly. $r_i$ tends to $1/h_i$ as $\mu_i h_i \to 0$. It is fitted at the baseline mortality $\mu_i(0)$ and held constant, so it does not follow set_demography() in time. A group must be wider than a day, so that $\rho_i < 1$. The stationary age structure of the daily update is the same as the continuous chain’s, so the fitting carries over to the daily clock.

B.2 Exposure and the day’s infection probability

The IBM saves each species’ $a_s I_{M,s}$ every day and draws today’s bites from the value saved $\tau_E$ days ago (lagged_eir$get(t - de)). fleet keeps the same delay line and, for a fractional lag, interpolates linearly between the two saved days as LaggedValue does. (For a fractional de the IBM’s buffer is one day too short and returns its initial value indefinitely; fleet does not reproduce this. The default de = 12 is whole.)

$$\bar\varepsilon(t) = \sum_s \left[(1-\varphi_E)\,a_s I_{M,s}\big(t-\lfloor\tau_E\rfloor\big) + \varphi_E\,a_s I_{M,s}\big(t-\lfloor\tau_E\rfloor-1\big)\right], \qquad \varphi_E = \tau_E - \lfloor\tau_E\rfloor.$$

The IBM draws the day’s bites as a Poisson total of $\bar\varepsilon$ times the population mean of $\psi$ and shares them in proportion to $\zeta\psi$ over the live population (human_pi()), so a person’s expected bites are

$$\varepsilon_{ij} = \bar\varepsilon\,\zeta_j\,\psi_i\,\frac{\overline{\psi}}{\overline{\zeta\psi}},$$

with both means taken over the current population. ($\overline{\zeta}$ is the quadrature’s $\sum_j w_j\zeta_j$, 0.99972 at five nodes, not 1.)

The IBM collects the day’s bitten in a set, so a person bitten several times in one timestep can be infected at most once. They are bitten at least once with probability $q^b_{ij} = 1 - e^{-\varepsilon_{ij}}$, and a bitten person in $S$, $A$ or $U$ is infected with probability $b_{ij}v_i(t)$, so

$$p_{ij} = q^b_{ij}\, b_{ij}\, v_i(t).$$

Setting parameters$bite_dedup = 0 lets every bite infect independently instead, with hazard $b_{ij}\varepsilon_{ij}v_i$. This agrees with the IBM to within 2% at low exposure and exceeds it at seasonal peaks and in high-$\zeta$ strata.

One competing draw. In the IBM, infection and disease progression are resolved together, once a day per person. Someone facing infection at hazard $h^b$ and progression at hazard $h$ has an event with probability $1 - e^{-(h^b + h)}$, and it is the infection with probability $h^b/(h^b + h)$, so an infection pre-empts that day’s progression. Only the day’s bitten carry an infection hazard, $h^b_{ij} = -\log(1 - b_{ij}v_i)$; the unbitten face progression alone. $S$ has no progression and $D$ and $T$ cannot be infected, so the draw matters for $A$ ($A \to U$ at $h_A = 1/d_A$) and $U$ ($U \to S$ at $h_U = 1/d_U$):

$$i^A_{ij} = q^b_{ij}\left(1 - (1 - b_{ij}v_i)\,e^{-h_A}\right)\frac{h^b_{ij}}{h^b_{ij} + h_A}, \qquad g^{AU}_{ij} = q^b_{ij}\left(1 - (1 - b_{ij}v_i)\,e^{-h_A}\right)\frac{h_A}{h^b_{ij} + h_A} + \left(1 - q^b_{ij}\right)\left(1 - e^{-h_A}\right),$$

and $i^U$, $g^{US}$ likewise with $h_U$. The day’s infections are $n_{ij} = p\,S + i^A A + i^U U$; a fraction $\phi$ of them is clinical, so the daily clinical count is exactly the IBM’s.

B.3 Immunity-dependent probabilities

All four probabilities are Hill functions of the immunity states. $\delta$ is acquired_immunity_offset: 0 by default for falciparum, and 0.5 reproduces the IBM’s per-individual +0.5. The IBM applies that offset only where acquired immunity is strictly positive, so never-boosted individuals get none. This is one reason a flat +0.5 on a stratum mean is the worse match, and why 0 is the default.

$$b_{ij} = b_0\left[b_1 + \frac{1-b_1}{1+\left(\frac{I_{B,ij}+\delta}{I_{B0}}\right)^{k_b}}\right]$$

$$\phi_{ij} = \phi_0\left[\phi_1 + \frac{1-\phi_1}{1+\left(\frac{I_{CA,ij}+\delta+I_{CM,ij}}{I_{C0}}\right)^{k_c}}\right]$$

$$\theta_{ij} = \theta_0\left[\theta_1 + \frac{1-\theta_1}{1+f_{v,i}\left(\frac{I_{VA,ij}+\delta+I_{VM,ij}}{I_{V0}}\right)^{k_v}}\right]$$

$$q_{ij} = d_1 + \frac{1-d_1}{1+\left(\frac{I_{D,ij}}{I_{D0}}\right)^{k_d} f_{d,i}}, \qquad c_{A,ij} = c_U + (c_D-c_U)\,q_{ij}^{\gamma_{\text{inf}}}$$

with the age-dependent detection and severity modifiers

$$f_{d,i} = 1 - \frac{1-f_{d0}}{1+(a_i/a_{d0})^{\gamma_d}}, \qquad f_{v,i} = 1 - \frac{1-f_{v0}}{1+(a_i/a_v)^{\gamma_v}} .$$

$q$ takes no offset, matching the IBM.

Maternal immunity is algebraic, not a state. The IBM draws each newborn’s mother uniformly from those aged $\lfloor \text{age}/365 \rfloor = 20$ in its heterogeneity group, so a newborn inherits, in expectation, their population-weighted mean immunity. fleet weights each age group by its share $m_i$ of that year; on the default grid, with edges at 20 and 21, $m_i$ is 1 for the groups making up $[20, 21)$ and 0 elsewhere. Maternal immunity decays with age from birth:

$$\overline{I_{CA}}_{20,j} = \sum_i m_i I_{CA,ij}, \qquad I_{CM,ij} = P_M\,\overline{I_{CA}}_{20,j}\,\kappa^{C}_i, \qquad I_{VM,ij} = P_{VM}\,\overline{I_{VA}}_{20,j}\,\kappa^{V}_i,$$

where $\kappa^{C}_i$ is the exact mean of $e^{-a/d_m}$ over age band $i$,

$$\kappa^{C}_i = \frac{d_m}{h_i}\left(e^{-a^-_i/d_m} - e^{-a^+_i/d_m}\right),$$

and $\kappa^V_i$ is the same with $d_{vm}$.

Mean-field bias in severe incidence

Of the three acquired-immunity Hill functions, $\theta$ has the shallowest exponent, not the steepest:

Function	Exponent	Scale
$b$	$k_b = 2.160$	$I_{B0} = 43.9$
$\phi$	$k_c = 2.369$	$I_{C0} = 18.02$
$\theta$	$k_v = 2.000$	$I_{V0} = 1.096$

What distinguishes $\theta$ is its scale. A Hill function is flat well below its scale parameter and approaches its full log-log slope only well above it, so its local slope depends on position. Writing $X = f_{v,i}\left(\left(I_{VA}+\delta+I_{VM}\right)/I_{V0}\right)^{k_v}$ for the denominator term,

$$\frac{\partial\log\theta}{\partial\log I_{VA}} = -\,\frac{k_v\left(1-\theta_1\right)X}{\left(1+X\right)\left(1+\theta_1X\right)},$$

which is near 0 for $X \ll 1$, reaches $-1.96$ at $X = \theta_1^{-1/2} \approx 92$, and returns towards 0 once $\theta$ reaches its $\theta_0\theta_1$ floor. The same expression with $k_c$ and $\phi_1$ governs $\phi$, whose steepest attainable slope, $-2.24$, is the steeper of the two.

Position separates them. $I_{V0}$ is sixteen times smaller than $I_{C0}$, while the two acquired immunities take nearly the same values at every age, so they differ almost entirely in what they are divided by. At fleet’s own seed at the reference EIR of 20:

Age	$I_{VA}/I_{V0}$	$\partial\log\theta/\partial\log I_{VA}$	$I_{CA}/I_{C0}$	$\partial\log\phi/\partial\log I_{CA}$
1	1.80	$-0.65$	0.12	$-0.02$
3	6.62	$-1.80$	0.44	$-0.30$
5	12.20	$-1.95$	0.82	$-0.91$
10	26.31	$-1.87$	1.79	$-1.88$

$\theta$ is within 5% of its steepest possible slope from age 3.4 to age 10.5. $\phi$ at age 3 is less than a fifth as steep, and does not overtake $\theta$ until age 10.

On that tail, $\theta$ is convex in $I_{VA}$ across almost the whole population, whereas $\phi$ is concave in $I_{CA}$ through early childhood at this EIR. This fixes the sign of the mean-field error within a stratum. fleet evaluates $\theta$ once per (age group $\times$ heterogeneity group), at that stratum’s mean $I_{VA}$; the IBM averages $\theta$ over individuals whose bite histories differ within the stratum, and by Jensen’s inequality $\mathbb{E}[\theta(I_{VA})] \geq \theta(\mathbb{E}[I_{VA}])$, so the mean field returns the smaller value. fleet resolves the log-normal biting strata explicitly (§2.1), so only the spread within an (age, $\zeta$) cell is lost, and that is small: refined to 417 age groups, fleet is within 1% of the IBM on all-age severe incidence at EIR 20, 50 and 120.

On a coarser grid the remaining error is the grid’s, with the same sign for the same reason one level up: a group’s severe incidence is evaluated at its mean immunity, and because $\theta$ is convex this under-counts wherever immunity changes across the group. It changes fastest in the young groups, where severe disease concentrates. The error is first order in the group width: fleet’s own age profile halves its distance from the converged one with each doubling of the number of groups. Against the IBM median, severe incidence at EIR 120 is 7% low on 53 groups and 2% low on 118, the default. vignette("using") covers what a coarser grid costs.

B.4 Disease-state equations

Writing $n_{ij} = p_{ij}S_{ij} + i^A_{ij}A_{ij} + i^U_{ij}U_{ij}$ for the day’s infections, $\text{SPC}(t)$ for the slow-parasite-clearance fraction of the treated (set_antimalarial_resistance()), and $f_t^{\text{eff}}(t) = f_t(t)\,\text{eff}(t)\,(1-\text{ETF}(t))$ for the effective treated fraction, and omitting the ageing terms of §B.1 that every compartment carries:

$$S'_{ij} = S_{ij} + g^{US}_{ij}U_{ij} + k_P r_P P_{k_P,ij} + k_{P_c} r_{P_c} P_{c,k_{P_c},ij} + k_{P_{ct}} r_{P_{ct}} P_{ct,k_{P_{ct}},ij} - p_{ij}S_{ij}$$

$$D'_{ij} = D_{ij} + \left(1 - f_t^{\text{eff}}\right)\phi_{ij}n_{ij} - r_D D_{ij}$$

$$T'_{ij} = T_{ij} + (1-\text{SPC})\,f_t^{\text{eff}}\,\phi_{ij}n_{ij} - r_T T_{ij}, \qquad T'_{s,ij} = T_{s,ij} + \text{SPC}\,f_t^{\text{eff}}\,\phi_{ij}n_{ij} - r_T^{\text{slow}} T_{s,ij}$$

$$A'_{ij} = A_{ij} + (1-\phi_{ij})\left(p_{ij}S_{ij} + i^U_{ij}U_{ij}\right) - \phi_{ij}\,i^A_{ij}A_{ij} + r_D D_{ij} - g^{AU}_{ij}A_{ij}$$

$$U'_{ij} = U_{ij} + g^{AU}_{ij}A_{ij} - i^U_{ij}U_{ij} - g^{US}_{ij}U_{ij}$$

$$P'_{1,ij} = P_{1,ij} + r_T T_{ij} + r_T^{\text{slow}} T_{s,ij} - k_P r_P P_{1,ij}, \qquad P'_{m,ij} = P_{m,ij} + k_P r_P\left(P_{m-1,ij} - P_{m,ij}\right), \quad m = 2,\dots,k_P$$

$$P'_{c,1,ij} = P_{c,1,ij} - k_{P_c} r_{P_c} P_{c,1,ij}, \qquad P'_{c,m,ij} = P_{c,m,ij} + k_{P_c} r_{P_c}\left(P_{c,m-1,ij} - P_{c,m,ij}\right), \quad m = 2,\dots,k_{P_c}$$

$$T'_{c,ij} = T_{c,ij} - r_T T_{c,ij}, \qquad T'_{cs,ij} = T_{cs,ij} - r_{T,c}^{\text{slow}} T_{cs,ij}, \qquad P'_{ct,1,ij} = P_{ct,1,ij} + r_T T_{c,ij} + r_{T,c}^{\text{slow}} T_{cs,ij} - k_{P_{ct}} r_{P_{ct}} P_{ct,1,ij},$$

with the later stages of $P_{ct}$ as those of $P$, and the infectivity stocks $J_c$ and $J_{cs}$ leaving with $T_c$ and $T_{cs}$.

A clinical case is treated successfully with probability $f_t^{\text{eff}}$ and goes to $T$ or $T_s$; otherwise it goes to $D$. A non-clinical infection of someone in $S$ or $U$ goes to $A$, and an $A$ infected sub-clinically stays in $A$.

Prophylaxis chains. The IBM has no prophylaxis state. It multiplies each treated person’s infection probability by the Weibull survival $W(t - t_{\text{drug}})$, and a dose acts from the day after it is given, so a treated cohort’s mean protection on day $n$ is $W(n)$. Read as the survival of a whole number of protected days $T_W$, $P(T_W \ge n) = W(n)$, its mean is $M = \sum_n W(n)$ and its variance $V = \sum_n (2n-1)W(n) - (\sum_n W(n))^2$. A chain of $k$ stages, each left with probability $k\,r$ a day, lasts a sum of $k$ geometric times, with mean $1/r$ and variance $1/(k r^2) - 1/r$, so it matches that mean and variance at $k = M^2/(V + M)$. For SP-AQ as chemoprevention this gives 10 stages, which put the chain’s survival within a few points of $W$ (0.67 at day 30, against the Weibull’s 0.70 and a single stage’s 0.42). A stage lasts at least a day, so $k \le M$, and $k$ is capped at 20.

$P$ is entered only after the $T$ sojourn, whereas the IBM’s clock runs from the dose in parallel with $T$: someone treated on day 0 is still in $T$ on day $n$ with probability $(1-r_T)^{n-1}$. The chain therefore carries the protected days left once out of $T$, $m = \sum_n \left(1 - (1-r_T)^{n-1}\right)W(n)$, and its stage count matches the variance of the whole $T + P$ sojourn, $k_P = m^2/(V - (1-r_T)/r_T^2 + m)$. This gives 9 stages for SP-AQ and 1 for AL, whose 10-day protection is already less variable than $T$ itself, so no chain length can sharpen it.

People in a chain are fully protected and are not infectious. This is the chain’s one departure from the IBM: its protection is all or nothing, where the IBM’s is partial and fades, so at high transmission the chain protects dosed children somewhat more than the IBM’s curve does (see the intervention-impact claim in fleetcheck).

Two treated compartments. The IBM assigns each successfully treated individual to standard or slow parasite clearance by a Bernoulli draw, giving a mixture of two sojourns, not one at the blended mean. $T$ and $T_s$ reproduce that split structurally.

Whole-day sojourns. The IBM advances states once a day with exit probability $1 - e^{-1/d}$, so its realised mean dwell is $1/(1-e^{-1/d})$ (5.52 d for $d_D = 5$, not 5). fleet uses the same probabilities, $r_X = 1 - e^{-1/d_X}$ for $X \in \{A, D, U, T, T_{\text{slow}}\}$, so the dwell distributions, not only their means, are the IBM’s. In the IBM, slow clearance is resolved through the same per-day path as ordinary treatment.

Chemoprevention compartments. $P_c$ and the treated phase $T_c$, $T_{cs}$ have no inflow in the daily update; they are filled only by the chemoprevention pulses of §E. The chains carry neither infection risk nor infectivity; the treated phase is detectable and infectious, as the IBM’s Tr is.

B.5 Immunity equations

The IBM boosts immunity on a per-day event with an integer refractory window: boost_immunity() fires only if timestep - last_boosted >= u, and the timestep is an integer, so a person boosted on day $t_0$ is refractory for the next $u_{\text{eff}} = \lceil u \rceil - 1$ days ($u = 0$ means no window).

The event differs between $I_B$ and the other three. For $I_B$ it is being bitten at least once, $q^b_{ij}$, whatever the person’s state, so its boosts are a renewal process with mean gap $u_{\text{eff}} + 1/q^b$ and daily boosting probability

$$\pi^B_{ij} = \frac{q^b_{ij}}{q^b_{ij}\,u^B_{\text{eff}} + 1}.$$

For $I_{CA}$, $I_D$ and $I_{VA}$ the event is infection, which also moves the person: someone boosted today is in $D$, $A$ or $T$ tomorrow, and whether they are still refractory when next infected depends on where they have got to. The renewal form, applied in every state, treats someone in $S$ or $U$ as though reinfected as often as someone in $A$, and so under-boosts. A window lasts days and the sojourns in $A$ and $U$ last months, so fleet uses a quasi-steady closure. From today’s rates, let $a^c_{ij,d}$ be the probability that someone boosted $d$ days ago is in state $c \in \{A, D, U\}$ today:

$$a^A_{1} = 1 - \phi, \quad a^D_{1} = \phi\,(1 - f_t^{\text{eff}}), \quad a^U_1 = 0,$$

$$a^A_{d+1} = a^A_d\left(1 - g^{AU} - \phi\,i^A - \mu\right) + a^D_d\,r_D + (1 - \phi)\,a^U_d\,i^U,$$

$$a^D_{d+1} = a^D_d\left(1 - r_D - \mu\right) + \phi\,(1 - f_t^{\text{eff}})\left(a^A_d\,i^A + a^U_d\,i^U\right), \qquad a^U_{d+1} = a^U_d\left(1 - i^U - g^{US} - \mu\right) + a^A_d\,g^{AU}$$

(the treated leave for $T$ and the protection chain, and cannot be infected before the window ends). $K^c_{ij} = \sum_{d=1}^{u_{\text{eff}}} a^c_{ij,d}$ is the number of refractory days spent in $c$ per boost. With $n^X_{ij}$ boosting infections a day, the refractory stock is $K^A n^X$ in $A$ and $K^U n^X$ in $U$. Those people are infected at $A$’s and $U$’s own rates without boosting, so $n^X = n - n^X(K^A i^A + K^U i^U)$ and

$$\pi^X_{ij} = \frac{1}{N_{ij}}\,\frac{n_{ij}}{1 + K^{A,X}_{ij}\,i^A_{ij} + K^{U,X}_{ij}\,i^U_{ij}}, \qquad X \in \{CA, D, VA\},$$

each with its own window. The closure omits the refractory who are back in $S$ within the window: from $U$ they are of order $10^{-3}$ of the refractory days in $A$, and a treated case reaches $S$ only through a prophylaxis chain.

A boost replaces the day’s decay. The IBM queues each immunity’s decay, $I \to I\,e^{-1/d}$, first in the day, and a boost, $I \to I + 1$ from the start-of-day value, later. The queue applies both in order and the later update overwrites the earlier for the people it covers (boost_immunity() is called with no decay of its own), so a boosted person gains 1 and does not decay that day. Averaged over a cell,

$$I'_{X,ij} = e^{-1/d_X} I_{X,ij} + \pi^X_{ij}\left(1 + I_{X,ij}\left(1 - e^{-1/d_X}\right)\right) + \alpha_{ij}\left(I_{X,i-1,j} - I_{X,ij}\right), \qquad X \in \{B, CA, D, VA\},$$

where $\alpha_{ij}$ is the day’s inflow into the cell per head (births at $i = 1$, ageing from the group below otherwise), with $I_{X,0,j} \equiv 0$ (immunity is zero at birth). The immunities hold the mean immunity of the people in a cell, so ageing moves them as a difference, and the dead carry the cell mean, so death leaves it unchanged. At the boosting levels of high transmission the skipped decay raises equilibrium immunity by a factor of up to $1/(1 - \pi)$, at most about a sixth; the continuous-time formula $\dot I = \beta - I/d$ misses this.

C. Human-to-mosquito infectivity

Onward infectivity is the state-weighted, biting-weighted human infectiousness. The IBM saves it every day and infects mosquitoes with the value from $\tau_l$ days ago (lagged_infectivity$get(t - delay_gam)). fleet keeps the same delay line, interpolating linearly between the two days around a fractional lag as LaggedValue does ($\tau_l = 12.5$ d is read halfway between the values saved 12 and 13 days ago). TBV enters here as per-state, per-age multipliers $\tau^U_i, \tau^A_i, \tau^D_i, \tau^T_i$:

$$\text{inf}_{ij} = c_D D_{ij}\tau^D_i + c_{A,ij} A_{ij}\tau^A_i + c_U U_{ij}\tau^U_i + \left(c_T(t)\left(T_{ij}+T_{s,ij}\right) + J_{c,ij} + J_{cs,ij}\right)\tau^T_i, \qquad \overline{\text{inf}}(t) = \frac{1}{\overline{\zeta\psi}}\sum_{i,j}\zeta_j\psi_i\,\text{inf}_{ij},$$

$$\Lambda^M_s = a_s(t)\left[(1-\varphi_l)\,\overline{\text{inf}}\big(t-\lfloor\tau_l\rfloor\big) + \varphi_l\,\overline{\text{inf}}\big(t-\lfloor\tau_l\rfloor-1\big)\right], \qquad \varphi_l = \tau_l - \lfloor\tau_l\rfloor,$$

where $\overline{\zeta\psi}$ is the mean of $\zeta\psi$ over the live population, as the IBM normalises it.

D. Mosquito model

The IBM’s mosquito model is a compartmental ODE, integrated across each day with that day’s biting rate, adult mortality and force of infection held fixed; fleet does the same. Per species $s$, with adult total $M_s = S_{M,s} + E_{M,s} + I_{M,s}$ and larval total $n_{L,s} = E_s + L_s$, the equations over the day are

$$\dot{E}_s = \beta_s M_s - \frac{E_s}{d_E} - \mu_E E_s\left(1 + \frac{n_{L,s}}{K_s(t)}\right)$$

$$\dot{L}_s = \frac{E_s}{d_E} - \frac{L_s}{d_L} - \mu_L L_s\left(1 + \gamma\frac{n_{L,s}}{K_s(t)}\right)$$

$$\dot{P}_{L,s} = \frac{L_s}{d_L} - \frac{P_{L,s}}{d_P} - \mu_P P_{L,s}$$

$$\dot{S}_{M,s} = \tfrac{1}{2}\frac{P_{L,s}}{d_P} - S_{M,s}\Lambda^M_s - \mu_s(t)S_{M,s}$$

Extrinsic incubation. malariasimulation keeps a queue of $\lceil\tau_M\rceil$ days. Each day it pushes that day’s flux into incubation, $S_M\Lambda^M$ at the start of the day, and pops the oldest. The day’s integration releases the value at the front, which is the flux that entered $\lceil\tau_M\rceil - 1$ days ago (9 days at the default $\tau_M = 10$), with the incubation survival $e^{-\mu\tau_M}$ at today’s $\mu$ (adult_mosquito_eqs.cpp). fleet keeps the same delay line:

$$m_s = S_{M,s}\Lambda^M_s\big(t - \lceil\tau_M\rceil + 1\big)\,e^{-\mu_s(t)\tau_M},$$

$$\dot{E}_{M,s} = S_{M,s}\Lambda^M_s - m_s - \mu_s(t)E_{M,s}, \qquad \dot{I}_{M,s} = m_s - \mu_s(t)I_{M,s},$$

with $m_s$ fixed across the day.

Carrying capacity. The IBM’s aquatic model evaluates $K_s$ at the solver’s time truncated to a whole day (carrying_capacity() takes an integer timestep), so $K_s$ is constant across each day at its start-of-day value. fleet holds it there too.

Integrating the day. The larval equations are stiff: late-larval density-dependent mortality runs at about 5 a day at equilibrium, and far faster when a run starts with larvae far above a seasonal carrying capacity, as a seasonal run seeded at the annual mean does. fleet crosses the day in $n_s$ sub-steps (ode_tuning(n_sub =), default 32) of an exponential midpoint rule. Every compartment has the form $\dot X = I - R X$, and a step of length $h$ takes it to $X e^{-Rh} + (I/R)(1 - e^{-Rh})$, which is exact for frozen $I$ and $R$ and stays stable and positive however large $R$ is. Each sub-step predicts the half-step state and then takes the full step with $I$ and $R$ from it, which makes the rule second order: once a run is on its cycle the error falls fourfold with each doubling of the sub-steps. At 32 sub-steps a seasonal run’s daily EIR is within $2.4\times10^{-4}$ of the converged solution at the seasonal trough and $5.5\times10^{-6}$ at its peak, comparable to the per-step tolerance to which the IBM integrates its own mosquito model (r_tol = a_tol = 1e-4).

Vector control enters only through $a_s(t)$ and $\mu_s(t)$. For each species and time, fleet uses malariasimulation’s own net and IRS decay and combination algebra, averaged over the net-using and sprayed fractions, to obtain a mean repellency $\bar{r}$ and a mean survival $\bar{p}$; then

$$W = (1-Q_0) + Q_0\,\bar{p}, \qquad Z = Q_0\,\bar{r}, \qquad f = \frac{1}{\dfrac{t_f}{1-Z} + \left(\dfrac{1}{f_0} - t_f\right)},$$

$$a_s = \left(1 - \frac{1-Q_0}{W}\right) f, \qquad \mu_s = -f\log\!\left(\frac{p_1^{(0)}W}{1 - Zp_1^{(0)}}\;e^{-\mu^{(0)}(1/f_0 - t_f)}\right), \qquad p_1^{(0)} = e^{-\mu^{(0)}t_f}.$$

Oviposition is unaffected: the IBM’s eggs_laid(beta, mu, f) is algebraically identical to $\beta$ for all $\mu, f$, so $\beta_s$ is constant and vector control acts only through $\mu$ and $a$.

$a_s(t)$ and $\mu_s(t)$ are read once a day, as the IBM recomputes them every timestep, from series that follow each spray round’s logistic decay and each net distribution’s retention and decay.

E. Chemoprevention

MDA, SMC and PMC are applied as state jumps between days, as the IBM applies them in update_mass_drug_admin(). The run stops at the end of every day that carries a scheduled round. At each round the target age band $[\ell, u)$ is mapped onto the model age grid by fractional overlap

$$o_i = \frac{\max\!\left(0,\ \min(u, a^+_i) - \max(\ell, a^-_i)\right)}{a^+_i - a^-_i},$$

so narrow bands (for example PMC dose ages) are captured rather than dropped, and coarse groups are not over-treated. A fraction

$$\text{frac}_i = o_i \times \text{coverage} \times \text{drug efficacy} \times (1 - \text{ETF})$$

of every state in each affected cell is treated successfully and routed by state, as the IBM routes it:

• the clinical, $D$, and the LM-detectable asymptomatic, $q\,A$, go to the treated phase $T_c$ (or $T_{cs}$, the drug’s slow-clearing share under resistance). They stay detectable and infectious there, at their infectivity times the drug’s drug_rel_c, for the Tr stay, $1/r_T$ (or its slow counterpart), and then join the chain $P_{ct}$;
• everyone else treated moves to the first stage of $P_c$: $S$, $U$, the undetected $(1-q)A$, and the existing $T$, $T_s$, $T_c$, $T_{cs}$ and prophylaxis stages. The IBM resets drug_time for everyone it treats, so the protection of those already covered is renewed.

$P_c$ runs down its chain to $S$, each stage left with probability $k_{P_c} r_{P_c}$ a day: $1/r_{P_c}$ is the protected days of the coverage-weighted mixture of the chemoprevention drugs, and $k_{P_c}$ is matched to its variance as in §B.4. $P_{ct}$ carries the protection left once out of the Tr stay, sized as the post-treatment chain is (§B.4), since the dose’s clock runs from the round in parallel with that stay. Both are distinct from the clinical chain.

A round takes effect the day after its scheduled timestep, as an IBM event does: it fires after the day’s processes, and its updates land at the end of the day. PMC, which the IBM delivers on an age trigger, is approximated as pulses every 30 days over each dose-age band.

F. Initial conditions

1. The malariasimulation list is back-translated to malariaEquilibrium parameter names (R/translate_params.R), and parameters$eq_params is merged over the result, taking precedence. set_equilibrium() always stores its own back-translation there, even when called with eq_params = NULL, so after any set_equilibrium() call it is that stored set, not R/translate_params.R, that supplies every shared constant.
2. For each heterogeneity node $j$, malariaEquilibrium::human_equilibrium_no_het() is solved on the model age grid at $\text{EIR} = \text{init_EIR}\times\zeta_j$ and treated fraction $f_t$, with $f_t$ the coverage in force at timestep 0, as the IBM draws its initial humans (calculate_eq()). This is raw coverage, as the IBM passes it to the solver, which sends that share of clinical cases to treatment; efficacy acts during the run, so a treated seed relaxes onto the model’s own treated state, as the IBM’s does. The solution’s $\Lambda$, $\phi$, prop, $r$, IB, ICA, ID, IVA and cA columns are used directly. The $S/T/D/A/U/P$ block is re-solved with a corrected prophylaxis ageing recursion, $b_P = (r_T b_T + r_{i-1}P_{i-1})/\beta_P$ (upstream divides only the $r_{i-1}P_{i-1}$ term by $\beta_P$, leaving $r_T b_T$ undivided), generalised to the $k_P$-stage chain: stage 1 is fed by $r_T T$, stage $m$ by $k_P r_P P_{m-1}$, and every stage ages at $r_{i-1}$. The re-solve makes the seed the fixed point of the prophylaxis chain implemented here, whose balance is the same on the daily clock as in continuous time.
3. The equilibrium age structure is recomputed under $\mu_i(0)$ (which may be the custom-demography baseline row), and the constant-hazard seed is rescaled onto it.
4. $T_0$ is split between $T$ and $T_s$ by the $t = 0$ slow-clearance fraction.
5. Mosquito compartments are seeded from initial_mosquito_counts() at the FOIM implied by the human equilibrium under the coverage in force at timestep 1, as set_equilibrium() sizes the IBM’s (its init_foim), with total_M chosen so that $\sum_s a_s I_{M,s} = \text{EIR}_{\text{seed}}/365$ exactly. Under the default demography $\text{EIR}_{\text{seed}} = \text{init_EIR}$. Under a custom demography total_M is fixed to the IBM’s (ibm_total_M(); see set_equilibrium()) and $\text{EIR}_{\text{seed}}$ is the root of total_M(EIR) = total_M_IBM, found by uniroot() on $\log \text{EIR}$, with steps 2 to 4 re-solved at each trial EIR. $K_s(0)$ is set at the corresponding carrying capacity, floored at $10^{-9}$ of the maximum for species with proportion 0, which would otherwise give $0/0$ in the larval term.
6. The three delay lines start full of the seed’s own values (the EIR $a_s I_{M,s}$, the incubation flux $S_{M,s}\Lambda^M_s$ and the humans’ infectivity), as the IBM’s lagged values default to their initial values and its incubation queue starts full.

The clocks in steps 2 and 5 coincide unless treatment starts after timestep 0. With the usual set_clinical_treatment(timesteps = 1), the IBM, and so fleet, starts untreated humans among mosquitoes sized for a treated population, and both relax from there.

malariaEquilibrium encodes the continuous-time forms: a linear force of infection, refractory boosting on the raw rate with the unrounded $u$, no at-risk restriction on boosting, continuous state exits with no competing draw, and immunity that decays through a boost. fleet reproduces the IBM’s per-day semantics, so the seed is a close approximation to fleet’s fixed point, not the fixed point itself. An undisturbed run’s PfPR(2-10) relaxes by under 0.5% and then holds (0.36% at EIR 20, flat to 0.006% over the last five of fifteen years). Incidence relaxes further, mostly within five years: under-5 clinical incidence by -0.3% at EIR 1, +1.5% at EIR 20 and +6% at EIR 120, and severe incidence by at most 2%. Burn a run in before reading its incidence. malariasimulation relaxes off its own malariaEquilibrium seed for the same reason, and by more, although the two seeds are not identical: the IBM solves on a fixed 0–99.9 y grid with heterogeneity handled inside human_equilibrium(), and does not correct the prophylaxis recursion.

G. Outputs

There is one row per day, 1..timesteps, as malariasimulation returns. Row $t$ holds the state at the start of day $t$ and the incidence during it, which is what the IBM renders at timestep $t$, so row 1 is the seed. The columns are the IBM’s: a parameter list gives the columns malariasimulation would give it, with the same names and meanings (tests/testthat/test-output-parity.R checks this against the IBM). Outputs are computed per age group, aggregated to output bands in R, as fractions of the total human population, and multiplied by human_population on output:

Output	Definition
n_age_*	$\sum_j N_{ij}$
n_detect_lm_*, p_detect_lm_*	$\sum_j \left(D + T + T_s + T_c + T_{cs} + q_{ij}A\right)$. The IBM samples each asymptomatic person’s detection into n_detect_lm_* and sums the probabilities into p_detect_lm_*; the expected count is both. Neither is a prevalence, which is n_detect_lm_* / n_age_*
n_detect_pcr_*	$\sum_j \left(D + T + T_s + T_c + T_{cs} + A + U\right)$: the IBM convention, not malariaEquilibrium’s sub-patent-weighted pos_PCR
n_inc_clinical_*, p_inc_clinical_*	$\sum_j \phi_{ij}n_{ij}$: the day’s clinical episodes, the IBM’s count and its expectation
n_inc_severe_*, p_inc_severe_*	$\sum_j \theta_{ij}n_{ij}$: $\theta$ times the day’s infections, because the IBM draws severe cases from the same infected set as clinical ones
n_inc_*	$\sum_j n_{ij}$, the day’s new infections, $pS + i^A A + i^U U$
p_inc_*	0: malariasimulation 3.0.0 renders it from incidence_min_ages, a field no parameter list sets
n_infections	The day’s new infections at all ages
ft	$f_t(t)$, rendered only when clinical treatment is deployed, as in the IBM
S_count, A_count, D_count, U_count, Tr_count	Population totals by state as the IBM counts them: S_count is $S$ plus every prophylaxis chain, Tr_count all four treated compartments
ica_mean, icm_mean, ib_mean, iva_mean, ivm_mean, id_mean	Each immunity’s mean over the population, $\sum_{ij} I_{ij} N_{ij} / \sum_{ij} N_{ij}$, and by age as <name>_mean_* with the band’s population as n_*, over <name>_rendering_*. The IBM seeds its infants’ maternal immunity from the equilibrium bin 0.1 years older than they are, so its icm_mean and ivm_mean run below these until the seeded infants age out
EIR, FOIM, Ph_count	Not in the IBM’s table. $365\,\bar\varepsilon$, the EIR per adult per year biting humans today (the IBM’s EIR_<species> is the same lagged quantity in total bites a day); $\Lambda^M_1$; and the prophylaxis chains, the drug-protected part of S_count

The vivax table differs where the IBM’s does (§V.6).

Mapping age groups onto bands. A band [lower, upper] holds the whole-day ages lower to upper, both ends included, as the IBM’s in_age_range() does. Nobody is under a day old when the IBM renders a day, since births land at the end of one, so on fleet’s age axis, which starts at 0 at birth, the band spans $[\max(1, \lceil\text{lower}\rceil) - 1, \lfloor\text{upper}\rfloor)$ days: 0_1824 and 1825_5474 tile the axis, and 0_1825 and 1825_5475 share one day’s cohort, as in the IBM. A model age group contributes the fraction of its own width inside that span, so over bands that tile a stretch of the axis the weights sum to exactly 1 per group inside it: no family double-counts, and none loses anyone it renders. The absorbing top group is the exception, since it is not a finite interval: a band reaching above its lower edge takes all of it, and a band ending exactly on that edge takes none of it, with a warning.

Per-family bands. Each output family is emitted only over its own *_rendering_ages band list, as in malariasimulation, with n_age_* over the union; an empty list renders nothing. postie treats every output column as an independent age stratum, so if a family were emitted over the union of all bands, the same ages would appear in two columns and any person-day-weighted aggregate would count them twice, inflating all-age clinical incidence by a factor of about 1.21.

H. Daily clock

The model is compiled from odin2 to C++ (src/malaria_daily.cpp) and advanced by dust2 as a discrete-time system, one update a day; only the mosquito model is sub-stepped within the day (§D). There is no solver tolerance, and the sub-step count is the only time-stepping setting. The update from time $t$ to $t + 1$ is the IBM’s timestep $t + 1$ and records that day’s counts, which is why the output rows start at day 1.

Every time-varying input is read at the start of the day. Every intervention series starts at its baseline value at $t = 0$, so the equilibrium seed is preserved. Two kinds of schedule meet the daily clock differently:

• Deployments (nets, spraying, vaccination, a carrying-capacity change) act in the IBM from the day after their timestep: nets, spraying and EPI vaccination are processes whose updates land at the end of the day, mass campaigns are events that fire after the day’s processes, and the carrying capacity is read at the start of the day. Their series carry a knot at $t_k$ that includes the deployment and one at $t_k - 1$ that does not, so the update that reads time $t_k$ (day $t_k + 1$) is the first to see it. The value at a knot is day $t_k + 1$’s, so a net’s, a spray’s or a vaccine’s efficacy is read one day after deployment, as in the IBM. The vector-control, PEV and TBV series are interpolated linearly between knots, which suits the within-round decay each carries; the $t_k - 1$ knot stops a deployment ramping in from the previous knot of a coarse background grid.
• Settings read on the day (treatment coverage, the drug mix, resistance, death rates) take effect in the IBM on their own timestep (match_timestep(ts, timestep)), so their knots are moved one day earlier and day $t_k$ reads the value scheduled for timestep $t_k$.

The PEV, TBV, carrying-capacity and (constant-demography) mortality grids extend a year past timesteps; the vector-control grid ends exactly at timesteps.

V. P. vivax model

A list from get_parameters(parasite = "vivax") runs a second human block in the same compiled model. The falciparum block is held empty, and under falciparum the vivax block collapses to one empty cell. The biting model (§B.2), the delay lines and the mosquito model (§D) are shared, fed by whichever block holds the population. A vivax list uses the falciparum names wherever the two parasites share a quantity (da, dd, dt, cd, cu, ct, rc, uc, rm, pcm and the clinical curve phi0, phi1, ic0, kc), so these reach the vivax block through the shared parameters.

V.1 The hypnozoite dimension

malariasimulation gives each person an integer count of hypnozoite batches, $k \in \{0, \dots, k_{\max}\}$. Relapse is a hazard $k f$, linear in $k$, but immunity, detectability and the U → S rate are not, and a mean count would not carry them: people with many batches relapse often, so they are more immune, and more often inside a refractory window (§V.3), than the average of their age and heterogeneity group. Giving every count its group’s pooled immunity puts clinical incidence 28% high at EIR 3. Every vivax compartment and immunity stock is therefore indexed $[i, j, \ell]$ (age, heterogeneity, level). The first levels hold the batch counts. When a radical-cure drug’s liver-stage protection outlasts the blood-stage protection given with it, further levels hold people whose liver stage is still drug-protected (§V.5).

Counts 0 to 3 and $k_{\max}$ have a level each, and the counts between share levels in pairs: $\{4, 5\}$, $\{6, 7\}$, $\{8, 9\}$ for the default $k_{\max} = 10$, giving eight levels for eleven counts. Within a pair the people follow a Poisson distribution. This is the shape a cell’s count takes when batches arrive with bites and are each lost at $\gamma_L$; the cap at $k_{\max}$ truncates it without changing the ratio of one count to the next. Its rate $\theta$ is fitted to the occupancies of the levels either side of the pair: for neighbours whose lowest counts are $a$ below and $c$ above,

$$\frac{N_{\text{above}}}{N_{\text{below}}} = \frac{\sum_{k \in \text{above}} \theta^k / k!}{\sum_{k \in \text{below}} \theta^k / k!},$$

solved by three fixed-point steps in $\log\theta$, each reducing the error at least $(c - a)$-fold. The pair’s top count then holds $\theta / (\theta + t)$ of its people, $t$ being the top count. These are the people a bite moves up a level; the rest, at the bottom count, are the only ones a lost batch moves down one. The level’s relapse hazard is $\bar k f$ at its people’s mean count $\bar k$. Fitting the rate to the neighbours, rather than reading it from counts 0 and 1, keeps the shape right under radical cure, which empties the high counts into 0. Against a ladder with a level per count (ode_tuning(batch_bins = FALSE)), the pairs move no output of the fleetcheck vivax suite by more than 0.4%, for a quarter less run time. Pairing $\{2, 3\}$ as well would move under-five clinical incidence at EIR 30 by 0.9%, against a replicate band 2.3% either side of the IBM’s median.

Each day, everyone:

• moves up one batch on an infection that came from a bite, capped at $k_{\max}$, and stays at the same count on a relapse. The bite share of a cell’s infections is $\lambda / (\lambda + k f)$, the IBM’s relative_rates, with $\lambda = b\,\varepsilon$ (§V.2);
• loses a batch with probability $1 - e^{-k \gamma_L}$, whatever their disease state (hypnozoite_batch_decay_process), unless a bite formed a batch that day. The IBM queues the decay first and the bite’s $\min(k + 1, k_{\max})$, from the start-of-day $k$, later, and the later update overwrites the earlier, at the top batch too, where the bite writes $k_{\max}$ back.

The decay is a separate draw, independent of the day’s infection and progression, so it is applied after them, to the day’s outcome less its batch-formers. Someone who relapses and loses a batch on the same day ends a batch down and infected, as in the IBM, and the dead, reset as newborns by a later update, do not decay at all.

V.2 Infection and progression

The day’s expected infectious bites per person, $\varepsilon$, are those of §B.2, shared over both blocks’ population. Vivax counts every bite rather than deduplicating them: a person bitten $n$ times is infected with probability $1 - (1 - b)^n$, which averaged over Poisson bites is $1 - e^{-b\varepsilon}$. Relapse adds $k f$, and the total is reduced by PEV, as calculate_vivax_infections() does, so the vaccine also blocks relapses:

$$\iota_S = \left(1 - e^{-(b\varepsilon + k f)}\right)(1 - \text{PEV}), \qquad h = -\log(1 - \iota_S).$$

The IBM resolves infection and a person’s own progression ($r$: A → U at $1/d_a$, D → A at $1/d_d$, U → S at $1/d_{\text{PCR}}(\text{IAA})$) in one competing draw a day: an event happens with probability $1 - e^{-(h + r)}$, and it is the infection with share $h / (h + r)$. So, per state,

$$\iota_X = \left(1 - e^{-(h + r_X)}\right) \frac{h}{h + r_X}, \qquad \rho_X = \left(1 - e^{-(h + r_X)}\right) \frac{r_X}{h + r_X},$$

with $r_S = 0$, and $d_{\text{PCR}}$ read at IAA + IAM. The bites are averaged before the draw. This is exact for S, whose infection has no competitor, but for D, A and U the IBM’s per-person draw is nonlinear in the day’s bite count (§V.7).

Everyone in S, U, A and D can be infected; an infected D stays in D. Infected S and U pass first through LM detectability $\phi_{LM}(\text{IAA} + \text{IAM})$: the rest of S go to U, and the rest of U stay in U. LM-detectable S and U, and all infected A, then pass through the clinical draw $\phi_D(\text{ICA} + \text{ICM})$. Clinical cases are treated with coverage × efficacy as for falciparum; the rest go to D. A is LM-detectable by definition, so LM prevalence counts D, Tr and all of A, and PCR prevalence adds U.

V.3 Immunity, its refractory windows and its spread

IAA (anti-parasite: LM detectability and the U → S duration) and ICA (clinical) are boosted by infection once the refractory window since the last boost has passed (44 days for IAA, 4 for ICA), and decay by $e^{-1/d}$ a day, a boosted person skipping the day’s decay as in §B.5. Immunity is held as stocks $J = \sum I$ per cell, not means, because the ladder has cells that are exactly empty, where a mean’s inflows would divide by zero.

Who is inside their window is also carried, as a refractory stock $R$ per cell. It moves with the people like any other stock (ageing, dying, losing batches, moved up a batch by a bite, taken by radical cure) and is released when the window ends; a day’s infections then boost in the share $1 - R/N$. A renewal rate $p/(pu + 1)$, the same in every cell, is correct for someone whose infection risk has been steady for the whole window, but not for vivax, because a bite both boosts a person and moves them up a batch into a cell whose relapse hazard is higher by $f$. A child moving to one batch relapses within IAA’s window two times in three, inside the window of the bite that brought the batch. The renewal rate cannot see that a cell’s newcomers are freshly boosted, and put young children’s IAA 10–11% above the IBM’s; the stock captures it. Each window is a chain of stages passed on with probability $n/u$ a day, on a clock applied after the day’s other moves. ICA’s window runs exactly, one stage a day, and IAA’s as four 11-day stages: an exact 44-day chain changes infants’ IAA by 2.7% of itself and no other output by more than 0.7%, at 2.4 times the cost.

People who share a cell do not share an immunity: their batch history differs, and years with few batches mean far fewer infections. In the IBM at EIR 10 the coefficient of variation of IAA within a cell is 0.5 in infancy and 0.13 by age 20, and the vivax curves are steep ($\kappa_C = 5.4$), so the mean clinical probability can be several times the probability at the mean immunity. Each cell therefore also carries $K = \sum I^2$, transported exactly as $J$ is, with a boost adding $2I + 1$ and decay multiplying $I^2$ by $e^{-2/d}$. The transport is exact because nothing that moves a person between cells depends on their own immunity once their cell is known. The variance a boost injects is that of the refractory renewal process, whose Fano factor $(1 - p)\,t^2$ ($t$ the share of infections that boost) is far below one for IAA’s 44-day window. The curves are averaged over a gamma distribution with each cell’s mean and variance, at five Gauss–Hermite nodes, through the Wilson–Hilferty map. Against the IBM’s individuals at EIR 10, the gamma average of the LM-detection curve is within 1% of the true cell average at every age, where the curve at the mean is up to 10% off, and the within-cell variance of IAA is 0.96–1.01 times the IBM’s from age 3 to 20.

A boost also moves its person’s mean, by the day’s boost probability $pt$, which adds $(pt)^2$ to $I^2$ over and above the variance. ICA carries that term too: with it the injection is exactly one boost’s worth per person-day with no window, and ICA’s spread stays within 6% of the IBM’s rule up to $p = 0.2$ a day. IAA’s long window makes its boosts so regular that the term overshoots, so IAA keeps the renewal term alone. That holds at the daily infection probabilities of EIR 10 and below; above about $p = 0.05$ it drifts from the IBM’s spread, in either direction depending on how steady a cell’s exposure has been (§V.7).

IAA and ICA share nodes: they are boosted by the same infections and correlate at 0.92–0.97 within a cell below age 40, so “LM-detectable and clinical” is averaged jointly. The nodes stand for people, so the clinical curve takes the IBM’s per-person $+0.5$ (acquired_immunity_offset = 0.5, the vivax default). parameters$immunity_spread = FALSE evaluates the curves at the cell mean, with no offset.

Maternal immunity is algebraic, as in §B.3: IAM and ICM are $p_{cm}$ times the population-weighted mean IAA and ICA of those aged 20 to 21 in the same heterogeneity group (the mothers the IBM draws from), waning at $1/r_m$.

V.4 Treatment

The treated leave Tr at $1 - e^{-1/d_t}$ a day into a post-treatment chain that carries the level dimension, since people keep losing batches while protected. Its stage count follows the daily rule of §B.4: eight stages for chloroquine.

V.5 Radical cure

Of everyone treated, a share drug_hypnozoite_efficacy lose every batch, whether or not the blood stage cleared, since the IBM draws the two independently. They start liver-stage protection, which blocks new batches (not infection) for a Weibull time from the dose. Protection is an attribute of the person, not of their disease state, so it sits on the level dimension as a chain of levels matched to the protected days, $M^2 / (V + M)$ stages: eleven for tafenoquine’s 27.

Liver-stage protection matters only while a bite can infect, and nobody in Tr or the post-treatment chain is infected, so protection that ends inside the blood-stage protection given with it guards only the few whose blood stage the drug did not clear. Primaquine’s 4.7 days end inside chloroquine’s 16 after Tr, so primaquine has no protected levels: its four would add half again to every vivax compartment, for an output change below 0.2%. A schedule whose longest liver-stage protection outlasts its shortest blood-stage protection keeps them, sized for the shortest-protecting radical-cure drug on it. A switch from primaquine to tafenoquine therefore keeps four stages: tafenoquine’s mean protection is kept and its spread widened, a standard deviation of 12.5 days against the Weibull’s 6.3. The last level returns to batch 0.

A primaquine stage on such a schedule is left with probability 0.85 a day, faster than Tr is left (0.63), so the chain’s clock acts on what the day leaves: someone who leaves Tr and moves a stage on the same day does both, and no cell can lose more than it holds. The day’s radically cured join the first level after the clock has acted, so a stage is left no earlier than the day after the dose. They take with them their immunity, their refractory windows, and the spread their boosts injected. The immunity they take is that of the nodes that fell ill, not their cell’s mean, since the clinical draw is taken node by node and the less immune fall ill more.

V.6 Seeding and outputs

The seed is the IBM’s own: malariaEquilibriumVivax on malariasimulation’s durations, solved on fleet’s age grid under the treatment coverage at timestep 0, with the mosquitoes sized under the coverage at timestep 1, as set_equilibrium() does. Both use raw coverage, as the IBM passes it to the solver. With heterogeneity off, the IBM still solves the heterogeneous equilibrium, sizes the mosquitoes from all of it, and draws everyone from its first group, the least bitten; fleet does the same. The IBM draws each person’s state and batch count from the equilibrium’s S, D, A, U and T within their age and heterogeneity group, omitting its prophylaxis state, so nobody starts drug-protected; fleet does the same. Every cell starts at its mean immunity (no spread), as initial_immunity() starts the IBM, and nobody starts inside a refractory window, as the IBM’s last_boosted starts at -1.

Neither model holds that state: both build up a within-cell spread the equilibrium does not have, so both relax from it together, and further than a falciparum run does. From EIR 10 down to 0.3 a vivax run’s LM prevalence settles 2 to 14% above the seed, its clinical incidence 15 to 21% above it and its realised EIR 4 to 17% above init_EIR, within about eight years at EIR 3 and over two decades at 0.3, where it overshoots first. The IBM moves as far: its realised EIR ends 18%, 13% and 8% above init_EIR at EIR 0.3, 1 and 3. With immunity_spread = FALSE, fleet holds its seed to within 1.3%, so the spread is what moves both. There is no infectivity lag: vivax has delay_gam = 0.

The output table is the IBM’s vivax table: no severe disease and no p_detect_lm_*, and in their place n_relapses, n_with_hypnozoites, the immunity means ica_mean, icm_mean, iaa_mean, iam_mean and hypnozoites_mean, and their age bands when incidence_relapse_rendering_*, n_with_hypnozoites_rendering_* and <immunity>_rendering_* are set (§G).

V.7 Limitations

• Immunity is correlated with disease state within a cell. The IBM’s immunity-dependent transitions sort people by immunity: clearance from U is faster with more IAA, so S holds more of the immune, and the clinical draw takes the least immune into D (their IAA is 0.6–0.8 times their cell’s mean from age 5). A mean field with one immunity distribution per cell cannot carry that sorting, and carrying immunity per disease state would double or triple the cost of a vivax run. At EIR 10 this leaves LM prevalence and clinical incidence in school-age children 1–3% above the IBM’s.
• Relapses on days without a bite. malariasimulation 3.0.0 evaluates relapses only on days when someone is bitten, so on a day with no infectious bite anywhere the IBM loses every relapse; fleet keeps them. This is a defect in the IBM and is not replicated. Such a day occurs with probability about $e^{-\bar\varepsilon N}$, for $\bar\varepsilon$ the mean infectious bites per person per day and $N$ the population: in a 10,000-person run it exceeds 5% only below an EIR of about 0.14, and in 50,000 below about 0.03, unless vector control reduces the day’s bites further.
• The day’s bite count is averaged before the competing draw. The IBM’s infection-or-progression draw and its bite-or-relapse split are per person, nonlinear in the day’s bite count $n$; fleet gives everyone in a cell the draw at the averaged hazard. This is exact for S, puts D’s infection probability 0.1–0.5% low and A’s and U’s by less, and puts the bite share of infections 0.2–1% low.
• Refractory people are infected at their cell’s average rate, although those in Tr and the drug-protected cannot be infected, so under treatment a cell’s infections fall inside a window slightly too often: an estimated up to 1.6% too few IAA boosts under age 2 at 50% chloroquine coverage and EIR 10.
• IAA’s within-cell spread at high transmission. The renewal closure (§V.3) is checked against the IBM at EIR 10 and below. Above a daily infection probability of about 0.05 it drifts from the IBM’s spread, and at the highest EIRs a cell’s variance can compute below zero and is read as none.
• Tr is fully protected, whereas the IBM’s prophylaxis on the day after treatment is $W(1)$, just below one for the vivax drugs.