MetroFluSim Mathematical Formulation

The mathematical framework is inspired by the immunoSEIRS model of the Meyers Lab (see Bi and Bandekar et al. 2023, Bi et al. 2022 and Bouchnita et al. 2021 for some related recent publications).

Please note that the MetroFluSim model requires extensive notation -- please review the entirety of each section for provided definitions.

Flu model: diagram

flu_model_diagram

Flu model: deterministic differential equations

We start off with common indices and arguments: - $t \in \mathbb N$: current simulation day - $\ell$: location, i.e. subpopulation, $\mathcal L$: set of all locations/subpopulations - $a$: age group, $\agegroups$: set of all age groups - $r$: risk group, $\riskgroups$: set of all risk groups - $i$: type of immunity-inducing event, $\mathcal{I} := \left\{\text{infection}, \text{vaccine}\right\}$: the set of all types of immunity-inducing events: infection and vaccination, respectively.

Population-level immunity

For each $\ell \in \mathcal L$, $a \in \agegroups$, $r \in \riskgroups$:

\[\begin{align*} \frac{dM\locationell\agerisktime}{dt} &= \frac{\rateRtoS(t) R\locagerisktime}{N\locationell_{a, r}} \cdot \left(1 - oM\locationell\agerisktime \right) - wM\locationell\agerisktime \tag{M1} \\ \frac{dMV\locationell\agerisktime}{dt} &= V\locationell\agerisk(t - \delta) - w_v MV\locationell\agerisktime \tag{M2} \end{align*}\]

where

$\rateRtoS$: rate at which recovered individuals become susceptible, so that $1/\rateRtoS$ is the average number of days a person is totally immune from reinfection until being susceptible again.
$V\locationell\agerisktime$: proportion of individuals residing in location $\ell \in \mathcal L$ in age-risk group $a$, $r$ who receive vaccines at time $t$.
$o$: positive constant modeling the saturation of antibody production in individuals who have infection-induced immunity.
$\delta$: number of days after dose for vaccine to become effective.
$w$: rate at which infection-induced immunity wanes.
$w_v$: rate at which vaccine-induced immunity wanes.

Compartment equations

Note that the following are all $\numagegroups \times \numriskgroups$ matrices:

$\boldsymbol{K}^I_{I}$: infection-induced reduction in infection risk.
$\boldsymbol{K}^H_{I}$: infection-induced reduction in hospitalization risk.
$\boldsymbol{K}^D_{I}$: infection-induced reduction in death risk.
$\boldsymbol{M}\locationell(t)$: infection-induced population-level immunity in location $\ell \in \mathcal L$.
$\boldsymbol{K}^I_{V}$: vaccine-induced reduction in infection risk.
$\boldsymbol{K}^H_{V}$: vaccine-induced reduction in hospitalization risk.
$\boldsymbol{K}^D_{V}$: vaccine-induced reduction in death risk.
$\boldsymbol{MV}\locationell(t)$: vaccine-induced population-level immunity in location $\ell \in \mathcal L$.

To simplify notation, we have the following terms that characterize the effect of population-level infection-induced immunities for a given subpopulation $\ell$, age $a$, and risk $r$:

\[\begin{align*} \LambdaIlocagerisktime &= \left[\frac{\boldsymbol{K_{I,a,r}^{I}}}{\boldsymbol{1}_{\lvert \mathcal L \rvert \times 1} - \boldsymbol{K_{I,a,r}^{I}}}\right]^T \boldsymbol{M_{a,r}\locationell(t)} \tag{I1} \\ \LambdaHlocagerisktime &= \left[\frac{\boldsymbol{K_{I,a,r}^{H}}}{\boldsymbol{1}_{\lvert \mathcal L \rvert \times 1} - \boldsymbol{K_{I,a,r}^{H}}}\right]^T \boldsymbol{M_{a,r}\locationell(t)} \tag{I2} \\ \LambdaDlocagerisktime &= \left[\frac{\boldsymbol{K_{I,a,r}^{D}}}{\boldsymbol{1}_{\lvert \mathcal L \rvert \times 1} - \boldsymbol{K_{I,a,r}^{D}}}\right]^T \boldsymbol{M_{a,r}\locationell(t)} \tag{I3} \end{align*}\]

where $\boldsymbol{1}_{\lvert \mathcal L \rvert \times 1}$ is an $\lvert \mathcal L \rvert \times 1$ vector of $1$'s and the fraction notation indicates element-wise division.

The vaccine-induced reductions in risk $\boldsymbol{K}^I_{V}$, $\boldsymbol{K}^H_{V}$, and $\boldsymbol{K}^D_{V}$ are calculated based on reported vaccine effectiveness over an entire season to take into account waning:

\[\begin{align*} K^{X}_{V} &= \frac{w_v \cdot T \cdot VE^{X}_{\text{season}}}{(1 - e^{-w_v}) \sum_{\tau=t_0}^{T-1} \frac{\sum_{u = t_0}^{\tau} p_{\text{prot}}(u) e^{-w_v (\tau - u)}}{\sum_{u = t_0}^{\tau} p_{\text{prot}}(u)}} \tag{K1} \end{align*}\]

where

$X \in \{ I, H, D\}$.
$VE^{X}_{\text{season}}$: reported seasonal average vaccine effectiveness.
$T$: length of season in days.
$t_0$: first day an individual is protected through vaccination (first vaccination date plus vaccine effectiveness delay $\delta$).
$p_{\text{prot}}(u)$: proportion of all vaccinated individuals who become protected on day $u$.

If vaccine effectiveness does not wane ($w_v = 0$) we simply have $K^{X}_{V} = VE^{X}_{\text{season}}$.

For each $\ell \in \mathcal L$, $a \in \agegroups$, $r \in \riskgroups$, we have the following equations that characterize transitions between compartments:

\[\begin{align} \frac{dS\locagerisktime}{dt} &= \underbrace{\rateRtoS R\locagerisktime}_{\text{$R$ to $S$}} -\underbrace{S\locagerisktime \frac{\beta\locationell(t)}{\left(1 + \LambdaIlocagerisktime\right)} (1 - K^I_{V,a,r} \cdot MV_{a,r}\locationell(t)) \cdot \totalforceofinfection}_{\text{$S$ to $E$}} \tag{C1} \\[1.5em] \frac{dE\locagerisktime}{dt} &= \underbrace{S\locagerisktime \frac{\beta\locationell(t)}{\left(1 + \LambdaIlocagerisktime\right)} (1 - K^I_{V,a,r} \cdot MV_{a,r}\locationell(t)) \cdot \totalforceofinfection }_{\text{$S$ to $E$}} \\ & \qquad \qquad - \underbrace{\rateEtoI (1-\propIA) E\locagerisktime}_{\text{$E$ to $IP$}} - \underbrace{\rateEtoI \propIA E\locagerisktime}_{\text{$E$ to $IA$}} \tag{C2} \\[1.5em] \frac{dIP\locagerisktime}{dt} &= \underbrace{\rateEtoI (1-\propIA) E\locagerisktime}_{\text{$E$ to $IP$}} - \underbrace{\rateIPtoIS IP\locagerisktime}_{\text{$IP$ to $ISR$ and $ISH$}} \tag{C3} \\[1.5em] \frac{dISR\locagerisktime}{dt} &= \underbrace{\left(1-\frac{\pi^H }{1 + \LambdaHlocagerisktime}(1 - K^H_{V,a,r} \cdot MV_{a,r}\locationell(t)) \right) \rateIPtoIS IP\locagerisktime}_{\text{$IP$ to $ISR$}} \\ & \qquad \qquad - \underbrace{\rateISRtoR ISR\locagerisktime}_{\text{$ISR$ to $R$}} \tag{C4} \\[1.5em] \frac{dISH\locagerisktime}{dt} &= \underbrace{\frac{\pi^H }{1 + \LambdaHlocagerisktime} (1 - K^H_{V,a,r} \cdot MV_{a,r}\locationell(t)) \cdot \rateIPtoIS IP\locagerisktime}_{\text{$IP$ to $ISH$}} \\ & \qquad \qquad - \underbrace{\rateISHtoH ISH\locagerisktime}_{\text{$ISH$ to $HR$ and $HD$}} \tag{C5} \\[1.5em] \frac{dIA\locagerisktime}{dt} &= \underbrace{\rateEtoI \propIA E\locagerisktime}_{\text{$E$ to $IA$}} - \underbrace{\rateIAtoR IA\locagerisktime}_{\text{$IA$ to $R$}} \tag{C6} \\[1.5em] \frac{dHR\locagerisktime}{dt} &= \underbrace{\left(1-\frac{\pi^D_{a, r}}{1 + \LambdaDlocagerisktime} (1 - K^D_{V,a,r} \cdot MV_{a,r}\locationell(t)) \right) \rateISHtoH ISH\locagerisktime}_{\text{$ISH$ to $HR$}} \\ & \qquad \qquad - \underbrace{\rateHRtoR HR\locagerisktime}_{\text{$HR$ to $R$}} \tag{C7} \\[1.5em] \frac{dHD\locagerisktime}{dt} &= \underbrace{\frac{\pi^D_{a, r} }{1 + \LambdaDlocagerisktime} (1 - K^D_{V,a,r} \cdot MV_{a,r}\locationell(t)) \cdot \rateISHtoH ISH\locagerisktime}_{\text{$ISH$ to $HD$}} \\ & \qquad \qquad - \underbrace{\rateHDtoD HD\locagerisktime}_{\text{$HD$ to $D$}} \tag{C8} \\[1.5em] \frac{dR\locagerisktime}{dt} &= \underbrace{\rateISRtoR ISR\locagerisktime}_{\text{$ISR$ to $R$}} + \underbrace{\rateIAtoR IA\locagerisktime}_{\text{$IA$ to $R$}} + \underbrace{\rateHRtoR HR\locagerisktime}_{\text{$HR$ to $R$}} - \underbrace{\rateRtoS R\locagerisktime}_{\text{$R$ to $S$}} \tag{C9} \\[1.5em] \frac{dD\locagerisktime}{dt} &= \underbrace{\rateHDtoD HD\locagerisktime}_{\text{$HD$ to $D$}}. \tag{C10} \end{align}\]

where

The $\lambda$-terms are location/subpopulation mixing terms that we define in the next section on the travel model.
$\boldsymbol{N}\locationell$: $\numagegroups \times \numriskgroups$ matrix corresponding to total population in location $\ell \in \mathcal L$, where element $N\locagerisk$ is the total population of age group $a$ and risk group $\ell$ in location $\ell$.
$\beta\locationell(t) = \beta\locationell_0 (1 + q(t))$: time-dependent transmission rate per day for individuals residing in location $\ell \in \mathcal L$.
$q(t) = \xi \cdot \exp^{-180 * h(t)}$: humidity adjustment, where $\xi$ is the humidity impact factor and $h(t)$ is absolute humidity. This formula is taken from this paper.
$\propIA$: proportion exposed who are completely asymptomatic when infectious.
$r_{IP}$, $r_{IA}$: relative infectiousness (compared to infected symptomatic people) of infected presymptomatic and infected asymptomatic people respectively.
$\rateISRtoR, \rateHRtoR, \rateIAtoR$: recovery rates for infected symptomatic ($ISR$), hospital ($HR$), and infected asymptomatic ($IA$) compartments respectively, so that $1/\gamma$ is the average number of days it takes for an infected person not in the hospital to recover, and $1/\rateHtoR$ is analogous, but for an infected person in the hospital.
$\rateEtoI$: infection rate (both exposed to infected presymptomatic transition rate and exposed to infected asymptomatic transition rate), so that $1/\rateEtoI$ is the average number of days after exposure before a person becomes infectious.
$\rateIPtoIS$: infected presymptomatic to infected symptomatic (both $ISR$ and $ISH$) transition rate, so that $1/\rateIPtoIS$ is the average number of days that an infected person is presymptomatic before becoming symptomatic.
$\rateISHtoH$: hospitalization rate (infected to hospital transition rate), so that $1/\rateISHtoH$ is the average number of days a person is infected before going to the hospital.
$\rateHDtoD$: death rate from hospital, so that $1/\rateHDtoD$ is the average number of days a person spends in the hospital before dying.

$\boldsymbol{\propH}$: $\numagegroups \times \numriskgroups$ proportion hospitalized based on age-risk group $a, r$.

$\boldsymbol{\propD}$: $\numagegroups \times \numriskgroups$ in-hospital mortality rate (proportion who die based on age-risk group $a, r$).

Flu model: travel model

Exposure intensity

For each $\ell \in \mathcal L$, $k \in \mathcal L \setminus \{\ell\}$, $a \in \agegroups$, we have

\[ \totalforceofinfection = \lambda^{(\ell), \text{local}}\agetime + \lambda^{(\ell), \text{visitors}}\agetime + \lambda^{(\ell), \text{residents traveling}}\agetime. \tag{T1} \]

This can loosely can be interpreted as exposure intensity: the (weighted) proportion of the population that interacts with $\ell,a$ individuals that are infectious.

The decompositions model the following phenomenon:

$\lambda^{(\ell), \text{local}}\agetime$: rate at which individuals in location $\ell$ get exposed to infected people who live in location $\ell$ (these contacts occur in location $\ell$).
$\lambda^{(\ell), \text{visitors from } k}\agetime$: rate at which individuals in location $\ell$ get exposed to infected people who live in location $k$ but travel to location $\ell$ (these contacts occur in location $\ell$).
$\lambda^{(\ell), \text{residents traveling to } k}\agetime$: rate at which individuals in location $\ell$ get exposed to infected people while visting location $k$, which includes both people who live in location $k$ and people from all other locations visting location $k$ (these contacts occur in location $k$).

Note that we assume this exposure intensity is the same for a given age group regardless of risk group, so we do not have the $r$-subscript here.

Specifically, we have

\[\begin{align*} \lambda^{(\ell), \text{local}}\agetime &= \psi_a \cdot \propdaytravelelltoell_a \cdot \sum \limits_{a^\prime \in \agegroups} \phi\locationell_{a,a^\prime}(t) \frac{\tilde{I}\locationell_{a^\prime}(t)}{\tilde{N}\locationell_{a^\prime}(t)} \propdaytravelelltoell_{a^\prime} \tag{T2} \\[1.5em] \lambda^{(\ell), \text{visitors}}\agetime &= \psi_a \cdot \propdaytravelelltoell_a \cdot \sum_{k \in \mathcal{L} \setminus \{\ell\}} \underbrace{\left( \sum\limits_{a^\prime \in \agegroups} \phi\locationell_{a, a^\prime}(t) \frac{\tilde{I}\locationk_{a^\prime}(t)}{\tilde{N}\locationell_{a^\prime}(t)} m_{a^\prime}(t) \cdot \propdaytravelktoell \right)}_{\text{Each summand: } \lambda^{(\ell), \text{visitors from } k}\agetime} \tag{T3} \\[1.5em] \lambda_{a,r}^{(\ell), \text{residents traveling}}(t) &= \psi_a \cdot \sum_{k \in \mathcal{L} \setminus \{\ell\}} \left( \proptravelelltok \cdot m_a \cdot \sum\limits_{a^\prime \in \agegroups} \phi\locationk_{a, a^\prime}(t) \right. \notag \\ & \qquad \qquad \qquad \qquad \underbrace{ \left. \frac{ \tilde{I}\locationk_{a^\prime}(t) \cdot\propdaytravelktok_{a^\prime} + \sum_{k^\prime \in \mathcal{L} \setminus \{k\}} \tilde{I}\locationkprime_{a^\prime}(t) \cdot m_{a^\prime} \cdot \proptravelkprimetok }{\tilde{N}\locationk_{a^\prime}(t)} \right)}_{\text{Each summand: } \lambda\agerisk^{(\ell), \text{residents traveling to } k}(t)} \tag{T4} \end{align*}\]

where

\[ \tilde{I}\locationell_{a}(t) := \sum_{r \in \riskgroups} \left[ISR\locationell\agerisk(t) + ISH\locationell\agerisk(t) + r_{IP} IP\locationell\agerisk(t) + r_{IA} IA\locationell\agerisk(t)\right] \tag{T5} \]

and where $\effectiveNlocagetime$ defined below is the effective population in location $\ell \in \mathcal L$ and age group $a \in \agegroups$, at time $t$.

\[\begin{align*} \effectiveNlocagerisktime &= \Nlocagerisk + m_a \cdot \sum_{k \in \mathcal L \setminus \{\ell\}} \propdaytravelktoell \cdot (N^{(k)}_{a, r} - HR^{(k)}_{a,r}(t) - HD^{(k)}_{a,r}(t)) \\ &\quad\quad\quad - m_a \cdot \sum_{k \in \mathcal L \setminus \{\ell\}} \proptravelelltok \cdot (N\locagerisk - HR\locationell_{a,r}(t) - HD\locationell_{a,r}(t)) \tag{T6} \\[1.5em] \effectiveNlocagetime &= \sum_{r \in \riskgroups} \tilde{N}^{(l)}_{a, r} (t) \tag{T7} \end{align*}\]

Contact matrix

The contact matrix is defined as

$\phi_{a, a^\prime}\locationell(t)$: the number of contacts that individuals in age group $a \in \agegroups$ residing in location $\ell \in \mathcal L$ have with other individuals (regardless of location) in age group $a^\prime \in \agegroups$ on day $t$.

Let $\phi^{(\ell), \text{total}}$, $\phi^{(\ell), \text{work}}$, and $\phi^{(\ell), \text{school}}$ represent the total contact matrix, school contact matrix, and work contact matrix, respectively for subpopulation $\ell$.

Then the contact matrix for the subpopulation at time $t$ is $$ \phi^{(\ell)}(t) := \phi^{(\ell), \text{total}} - (1 - d_{\text{work}}(t)) \phi^{(\ell), \text{work}} - (1 - d_{\text{school}}(t)) \phi^{(\ell), \text{school}} $$ where $d_{\text{work}}(t)$ is $1$ if the real-world date corresponding to simulation time $t$ is a work day and $0$ otherwise, and $d_{\text{school}}(t) \in [0, 1]$ is defined analogously for school days and represents the proportion of schools open on day $t$.

Other travel parameters

We have

$\psi_a \in [0, 1]$: relative susceptibility of individuals in age group $a \in \agegroups$.
$m_a$: a positive scalar modifying travel intensity depending on age $a$.
$\propdaytravelktoell$: (on average) proportion of people time-weighted by proportion of the day that residents of $k$ spend traveling to location $\ell \in \mathcal L$.

Note that the arrows in $\propdaytravelktoell$ correspond to direction of travel (e.g. $k \rightarrow \ell$ represents residents of location $k$ traveling to location $\ell$). The $\propdaytravelktoell$ values are calculated from mobility data, corresponding to

\[\begin{align*} p^{k \rightarrow \ell} &= \sum_{\poiell} c^{\poiell} \cdot v^{k \rightarrow \poiell} \tag{T8} \\[1.5em] \propdaytravelelltoell_a &= 1 - m_a \sum_{k \in \mathcal{L} \setminus \{\ell\}} \proptravelelltok \tag{T9} \end{align*}\]

where

$c^{\poiell}$: average proportion of a day spent at $\poiell$ on a given visit.
$v^{k \rightarrow \poiell}$: average number of visits per day per resident of $k$ to $\poiell$.

Flu model: discretized stochastic implementation

To actually implement/simulate this compartmental model, we discretize the deterministic differential equations and treat transitions between compartments as stochastic to model uncertainty. We extend the notation from the deterministic differential equations to capture the stochastic elements.

Let $\boldsymbol{\mathcal X}(t) = \left\{\boldsymbol{S}(t), \boldsymbol{E}(t), \boldsymbol{IA}(t), \boldsymbol{IP}(t), \boldsymbol{ISR}(t), \boldsymbol{ISH}(t), \boldsymbol{HR}(t), \boldsymbol{HD}(t), \boldsymbol{R}(t), \boldsymbol{D}(t), \boldsymbol{M}(t), \boldsymbol{MV}(t), q(t), \boldsymbol{\phi}(t), V\locationell(t)\right\}$ be the "simulation state" at time $t$. $\boldsymbol{\mathcal X}(t)$ is a set of matrices.

Let $\boldsymbol{\Theta}$ be the set of fixed parameters.

Then given initial state $\boldsymbol{\mathcal X}_0 = \boldsymbol{\mathcal X}(0)$, we can formulate our discretized stochastic implementation as

\[ \boldsymbol{\mathcal X}(t + \Delta t) = \boldsymbol{\mathcal X}(t) + f\left(\boldsymbol{\mathcal X}(t), \Delta t, \omega; \boldsymbol{\Theta}\right) \quad \text{for} \quad t \ge 0, \]

where $f$ is parametrized by $\boldsymbol{\Theta}$, and depends on the step size of discretization $\Delta t$ and a sample path $\omega$. We assume that each sample path $\omega$ is realized from a random process that does not depend on $\boldsymbol{\mathcal X}(t)$ or $\Delta t$ for each $t$. When we are discussing a single model with a fixed set of parameters $\boldsymbol{\Theta}$, we drop the $\boldsymbol{\Theta}$ notation for simplicity.

Now we formulate how we implement discretized stochastic transitions. We assume that $q(t)$, $\boldsymbol{\phi}(t)$, and $V\locationell(t)$ are updated deterministically according to some "schedule."

We model stochastic transitions between compartments using "transition variables." Transition variables correspond to incoming and outgoing flows of epidemiological compartments (see the compartment equations above).

Below we formulate the discretized stochastic transitions. Note that the population-level immunity variables behave as aggregate epidemiological metrics. They are deterministic functions of the simulation state and transitions between compartments.

IMPORTANT NOTE_: all $y$ and $y^*$-variables depend on $\left(\boldsymbol{\mathcal X}(t), \Delta t, \omega\right)$. For notation simplicity, we define $\simstate := \left(\boldsymbol{\mathcal X}(t), \Delta t, \omega\right)$ and write $y$ and $y^*$-variables as functions of $\simstate$.

For each $\ell \in \mathcal L$, $a \in \agegroups$, and $r \in \riskgroups$:

\[\begin{align} M_{a, r}\locationell(t + \Delta t) &= M\locationell\agerisktime + \frac{dM\locationell\agerisktime}{dt} \cdot \Delta t \\ MV_{a,r}\locationell(t + \Delta t) &= MV\agerisktime\locationell + \frac{dMV\agerisktime\locationell}{dt} \cdot \Delta t \\ S\locagerisk(t + \Delta t) &= S_{a, r}\locationell(t) + \underbrace{\tvarloc_{R\rightarrow S, a, r}(\Xi_t)}_{\text{$R$ to $S$}} - \underbrace{\tvarloc_{S\rightarrow E, a, r}(\Xi_t)}_{\text{$S$ to $E$}} \\ E\locagerisk(t + \Delta t) &= E_{a, r}\locationell(t) + \underbrace{\tvarloc_{S\rightarrow E, a, r}(\Xi_t)}_{\text{$S$ to $E$}} - \underbrace{\jointtvarloc_{E\rightarrow IP, a, r}(\Xi_t)}_{\text{$E$ to $IP$}} - \underbrace{\jointtvarloc_{E\rightarrow IA, a, r}(\Xi_t)}_{\text{$E$ to $IA$}} \\ IP\locagerisk(t + \Delta t) &= IP\locagerisktime + \underbrace{\jointtvarloc_{E\rightarrow IP, a, r}(\Xi_t)}_{\text{$E$ to $IP$}} - \underbrace{\tvarloc_{IP \rightarrow ISR, a, r}(\Xi_t)}_{\text{$IP$ to $ISR$}} - \underbrace{\tvarloc_{IP \rightarrow ISH, a, r}(\Xi_t)}_{\text{$IP$ to $ISH$}} \\ IA_{a, r}\locationell(t + \Delta t) &= IA_{a, r}\locationell(t) + \underbrace{\jointtvarloc_{E\rightarrow IA, a, r}(\Xi_t)}_{\text{$E$ to $IA$}} - \underbrace{\tvarloc_{IA \rightarrow R, a, r}(\Xi_t)}_{\text{$IA$ to $R$}} \\ ISR\locagerisk(t + \Delta t) &= ISR_{a, r}\locationell(t) + \underbrace{\tvarloc_{IP \rightarrow ISR, a, r}(\Xi_t)}_{\text{$IP$ to $ISR$}} - \underbrace{\jointtvarloc_{ISR \rightarrow R, a, r}(\Xi_t)}_{\text{$ISR$ to $R$}} \\ ISH\locagerisk(t + \Delta t) &= ISH_{a, r}\locationell(t) + \underbrace{\tvarloc_{IP \rightarrow ISH, a, r}(\Xi_t)}_{\text{$IP$ to $ISH$}} - \underbrace{\jointtvarloc_{ISH \rightarrow HR, a, r}(\Xi_t)}_{\text{$ISH$ to $HR$}} - \underbrace{\jointtvarloc_{ISH \rightarrow HD, a, r}(\Xi_t)}_{\text{$ISH$ to $HD$}} \\ HR\locagerisk(t + \Delta t) &= HR_{a, r}\locationell(t) + \underbrace{\jointtvarloc_{ISH \rightarrow HR, a, r}(\Xi_t)}_{\text{$ISH$ to $HR$}} - \underbrace{\jointtvarloc_{HR\rightarrow R, a, r}(\Xi_t)}_{\text{$HR$ to $R$}} \\ HD\locagerisk(t + \Delta t) &= HD_{a, r}\locationell(t) + \underbrace{\jointtvarloc_{ISH \rightarrow HD, a, r}(\Xi_t)}_{\text{$ISH$ to $HD$}} - \underbrace{\jointtvarloc_{HD\rightarrow D, a, r}(\Xi_t)}_{\text{$HD$ to $D$}} \\ R\locagerisk(t + \Delta t) &= R_{a, r}\locationell(t) + \underbrace{\jointtvarloc_{IA \rightarrow R, a, r}(\Xi_t)}_{\text{$IA$ to $R$}} + \underbrace{\jointtvarloc_{ISR \rightarrow R, a, r}(\Xi_t)}_{\text{$ISR$ to $R$}} + \underbrace{\jointtvarloc_{HR\rightarrow R, a, r}(\Xi_t)}_{\text{$HR$ to $R$}} - \underbrace{\tvarloc_{R\rightarrow S, a, r}(\Xi_t)}_{\text{$R$ to $S$}} \\ D\locagerisk(t + \Delta t) &= D_{a, r}\locationell(t) + \underbrace{\jointtvarloc_{HD\rightarrow D, a, r}(\Xi_t)}_{\text{$HD$ to $D$}}. \end{align}\]

IMPORTANT NOTE: the "$*$" superscript indicates that the transition variable has a joint distribution with another transition variable. In general, if a compartment has more than one outgoing transition variable, these transition variables must be modeled jointly.

Consider the two transitions out of the infected symptomatic compartment, for example. Given that a patient is infected and symptomatic, exactly one outcome occurs: they recover (from home) or they go to the hospital. Since one and only one of these outcomes must occur, we must model these two transition variables jointly. Joint distribution derivation details are provided in the next section on transition types.

Each transition variable depends on a "base count" and a "rate" (which both depend on the current state of the system). This decomposition is displayed in the table below. Note that these transitions are for a given location $\ell \in \mathcal L$.

Transition	Transition variable	Base count	Rate
$R$ to $S$	$\tvarloc_{R \rightarrow S, a, r}(\simstate)$	$R_{a, r}\locationell(t)$	$\rateRtoS$
$S$ to $E$	$\tvarloc_{S \rightarrow E, a, r}(\simstate)$	$S_{a, r}\locationell(t)$	$\frac{\beta(t) \totalforceofinfection}{1 + \LambdaIlocagerisktime}$
$E$ to $IP$	$\jointtvarloc_{E \rightarrow IP, a, r}(\simstate)$	$E_{a, r}\locationell(t)$	$\rateEtoI (1 - \propIA)$
$E$ to $IA$	$\jointtvarloc_{E \rightarrow IA, a, r}(\simstate)$	$E_{a, r}\locationell(t)$	$\rateEtoI \propIA$
$IP$ to $ISR$	$\tvarloc_{IP \rightarrow ISR, a, r}(\simstate)$	$IP_{a, r}\locationell(t)$	$\rateIPtoIS \left(1-\frac{\pi^H_{a,r}}{1 + \LambdaHlocagerisktime}\right)$
$IP$ to $ISH$	$\tvarloc_{IP \rightarrow ISH, a, r}(\simstate)$	$IP_{a, r}\locationell(t)$	$\rateIPtoIS \frac{\pi^H_{a,r}}{1 + \LambdaHlocagerisktime}$
$ISR$ to $R$	$\jointtvarloc_{ISR \rightarrow R, a, r}(\simstate)$	$ISR_{a, r}\locationell(t)$	$\gamma^{ISR\rightarrow R}$
$ISH$ to $HR$	$\jointtvarloc_{ISH \rightarrow HR, a, r}(\simstate)$	$ISH_{a, r}\locationell(t)$	$\rateISHtoH \left(1-\frac{\pi^D_{a,r}}{1 + \LambdaDlocagerisktime}\right)$
$ISH$ to $HD$	$\jointtvarloc_{ISH \rightarrow HD, a, r}(\simstate)$	$ISH_{a, r}\locationell(t)$	$\rateISHtoH \frac{\pi^D_{a,r}}{1 + \LambdaDlocagerisktime}$
$IA$ to $R$	$\tvarloc_{IA \rightarrow R, a, r}(\simstate)$	$IA_{a, r}\locationell(t)$	$\rateIAtoR$
$HR$ to $R$	$\jointtvarloc_{HR \rightarrow R, a, r}(\simstate)$	$HR_{a, r}\locationell(t)$	$\rateHRtoR$
$HD$ to $D$	$\jointtvarloc_{HD \rightarrow D, a, r}(\simstate)$	$HD_{a, r}\locationell(t)$	$\rateHDtoD$

The base count and rate of a transition variable parameterize the distribution that defines its realization.

See this page for mathematical formulations of marginal and joint stochastic transitions between compartments.

General model: discretized stochastic implementation

We make the important note that the flu model's discretized stochastic implementation can be generalized to models with different structures. More broadly, we let $\boldsymbol{\mathcal C}(t)$ be a model's set of epidemiological compartments, $\boldsymbol{\mathcal M}(t)$ its set of aggregate epidemiological metrics, and $\boldsymbol{S (t)}$ its set of schedule-dependent (time-dependent) deterministic values. Then the above formulation still holds.

In fact, in our code, we model $\boldsymbol{\mathcal C(t)}$ using an Compartment class, $\boldsymbol{\mathcal M}(t)$ using an EpiMetric class, and $\boldsymbol{\mathcal S(t)}$ using a Schedule class. We handle stochastic transitions using TransitionVariable and TransitionVariableGroup classes. These classes form some of the building blocks of the base model code.

Updated 8/31/2026. Documentation written by LP and updated by Rémy Pasco, mathematical notation by LP (advised by Lauren Meyers and Dave Morton, edited by Susan Ptak, Meyers Lab, and Shiyuan Liang), travel model conceptualized by Rémy and Susan Ptak, immunity formulation by Anass Bouchnita.

Transition	Transition variable	Base count	Rate
\(R\) to \(S\)	\(\tvarloc_{R \rightarrow S, a, r}(\simstate)\)	\(R_{a, r}\locationell(t)\)	\(\rateRtoS\)
\(S\) to \(E\)	\(\tvarloc_{S \rightarrow E, a, r}(\simstate)\)	\(S_{a, r}\locationell(t)\)	\(\frac{\beta(t) \totalforceofinfection}{1 + \LambdaIlocagerisktime}\)
\(E\) to \(IP\)	\(\jointtvarloc_{E \rightarrow IP, a, r}(\simstate)\)	\(E_{a, r}\locationell(t)\)	\(\rateEtoI (1 - \propIA)\)
\(E\) to \(IA\)	\(\jointtvarloc_{E \rightarrow IA, a, r}(\simstate)\)	\(E_{a, r}\locationell(t)\)	\(\rateEtoI \propIA\)
\(IP\) to \(ISR\)	\(\tvarloc_{IP \rightarrow ISR, a, r}(\simstate)\)	\(IP_{a, r}\locationell(t)\)	\(\rateIPtoIS \left(1-\frac{\pi^H_{a,r}}{1 + \LambdaHlocagerisktime}\right)\)
\(IP\) to \(ISH\)	\(\tvarloc_{IP \rightarrow ISH, a, r}(\simstate)\)	\(IP_{a, r}\locationell(t)\)	\(\rateIPtoIS \frac{\pi^H_{a,r}}{1 + \LambdaHlocagerisktime}\)
\(ISR\) to \(R\)	\(\jointtvarloc_{ISR \rightarrow R, a, r}(\simstate)\)	\(ISR_{a, r}\locationell(t)\)	\(\gamma^{ISR\rightarrow R}\)
\(ISH\) to \(HR\)	\(\jointtvarloc_{ISH \rightarrow HR, a, r}(\simstate)\)	\(ISH_{a, r}\locationell(t)\)	\(\rateISHtoH \left(1-\frac{\pi^D_{a,r}}{1 + \LambdaDlocagerisktime}\right)\)
\(ISH\) to \(HD\)	\(\jointtvarloc_{ISH \rightarrow HD, a, r}(\simstate)\)	\(ISH_{a, r}\locationell(t)\)	\(\rateISHtoH \frac{\pi^D_{a,r}}{1 + \LambdaDlocagerisktime}\)
\(IA\) to \(R\)	\(\tvarloc_{IA \rightarrow R, a, r}(\simstate)\)	\(IA_{a, r}\locationell(t)\)	\(\rateIAtoR\)
\(HR\) to \(R\)	\(\jointtvarloc_{HR \rightarrow R, a, r}(\simstate)\)	\(HR_{a, r}\locationell(t)\)	\(\rateHRtoR\)
\(HD\) to \(D\)	\(\jointtvarloc_{HD \rightarrow D, a, r}(\simstate)\)	\(HD_{a, r}\locationell(t)\)	\(\rateHDtoD\)