Epidemic growth in a fast-switching environment
The early phase of SIR and SEIR epidemics whose transmission and recovery rates follow a fast hidden regime: the growth of the expected epidemic, the growth of a typical one, and the gap between them.
The model
Early in an epidemic almost everyone is susceptible, and the infected compartments grow linearly. With a latent stage $E$ and an infectious stage $I$, the linearized SEIR model is
$$E' = \beta\,I - \sigma E, \qquad I' = \sigma E - \gamma I .$$Here $\beta$ is the transmission rate, $1/\sigma$ the mean latent period and $1/\gamma$ the mean infectious period. The linearized SIR model has no latent stage:
$$I' = (\beta - \gamma)\,I .$$Now let the rates depend on a hidden regime $y_t$: seasons, swings in behaviour, or a policy that tightens and relaxes. The regime is a Markov chain with generator $Q$ and stationary law $\pi$.
It starts from $\pi$, since an observer cannot see it. In matrix form, with $x = (E, I)$,
$$x' = A(y_t)\,x, \qquad A(y) = \begin{pmatrix} -\sigma & \beta(y) \\ \sigma & -\gamma(y) \end{pmatrix}.$$There are two time scales. The disease moves through its stages over $1/\sigma$ and $1/\gamma$. The regime changes every $1/\lambda$, where $\lambda$ is a switching rate, and it changes much faster.
Parameters: rates per week; latent period two weeks, $\sigma = 0.5$; infectious period four weeks, $\gamma = 0.25$; transmission $\beta = (0.455,\ 0.055)$ on a symmetric two-state chain switching at rate $\lambda$, so $\bar\beta = 0.255$ and $\bar\beta/\gamma = 1.02$. The three-regime example below has its own parameters.
Only ratios of rates matter. Multiplying every rate by seven turns weeks into days and describes a faster infection under faster switching.
Two growth rates
Two growth rates describe the linear phase. The mean growth rate is the growth of the expected number infected:
$$\Lambda_{\mathrm{mean}} = \lim_{t\to\infty} \frac1t \log \mathbb E\big[\,I_t\,\big].$$The almost-sure growth rate is the growth of a typical path, the top Lyapunov exponent of the switched system:
$$\Lambda_{\mathrm{as}} = \lim_{t\to\infty} \frac1t \log I_t .$$By Jensen's inequality $\Lambda_{\mathrm{mean}} \ge \Lambda_{\mathrm{as}}$. The mean is carried by rare regime paths with long runs of high transmission.
The almost-sure rate decides whether an introduction takes off. For a branching process in a random environment, extinction is certain exactly when the typical log growth is not positive (Smith and Wilkinson 1969; Athreya and Karlin 1971). For epidemics with Markov switching the same exponent decides extinction (Gray, Greenhalgh, Mao and Pan 2012; Prodhomme and Strickler 2024).
The mean has an exact linear equation. Let $m_i(t) = \mathbb E[x_t;\ y_t = i]$. Then
$$m_i' = A_i\,m_i + \sum_j Q_{ji}\,m_j ,$$so $\Lambda_{\mathrm{mean}}$ is the principal eigenvalue of the combined matrix $Q^\top \otimes I + \operatorname{blockdiag}(A_i)$. The almost-sure rate has no such equation.
Averaging
When the regime switches infinitely fast, the rates are replaced by their averages under $\pi$:
$$\bar A = \begin{pmatrix} -\sigma & \bar\beta \\ \sigma & -\bar\gamma \end{pmatrix}.$$Both growth rates tend to its principal eigenvalue $r$, which has the sign of $\bar\beta/\bar\gamma - 1$. For SEIR,
$$r = \frac{-(\sigma + \bar\gamma) + \sqrt D}{2}, \qquad D = (\sigma - \bar\gamma)^2 + 4\sigma\bar\beta .$$With the parameters above, $\bar\beta/\gamma = 1.02$ and the averaged SEIR epidemic grows at $r = 0.00332$ per week. The averaged SIR epidemic grows at $\bar\beta - \gamma = 0.005$. The question is how fast switching moves each of the two growth rates away from these values.
One infected compartment: SIR
In the SIR linearization the infectious count solves a scalar equation, so given the regime path
$$I_t = I_0\,\exp\Big(\int_0^t \big(\beta(y_s) - \gamma\big)\,ds\Big).$$By the ergodic theorem the typical growth rate is the average rate, at every switching speed:
$$\Lambda_{\mathrm{as}} = \bar\beta - \gamma .$$The mean given the starting regime is $\mathbb E[I_t \mid y_0 = i] = I_0\,a_i(t)$, where
$$a' = \big(Q + \operatorname{diag}(\beta - \gamma)\big)\,a, \qquad a(0) = \mathbf 1 .$$This is the engine's equation $a' = (Q + \operatorname{diag} g)\,a$ with a constant forcing $g_i = \beta_i - \gamma$. Its growth rate is the principal eigenvalue of $Q + \operatorname{diag}(\beta - \gamma)$. To first order it is the averaged rate plus the Green–Kubo integral of the transmission rate, as for a hazard rate on the any-equation page, with the sign reversed.
For a symmetric two-state chain the eigenvalue is a quadratic root.
For any chain, $\Lambda_{\mathrm{mean}} = \bar\beta - \gamma + K_{\beta\beta} + O(\lambda^{-2})$, with $K_{\beta\beta} = -\pi\cdot(\tilde\beta\,Q^{\#}\tilde\beta)$ and $Q^{\#}$ the group inverse.
The gap $K_{\beta\beta}$ is the gap of Lewontin and Cohen (1969): the log of an average of exponentials exceeds the average exponent by half the variance of the integrated rate, and that variance grows like $2tK_{\beta\beta}$.
With the parameters above and $\lambda = 1$, so that each regime lasts a week on average, the expected number infected grows 26.9% over ten weeks, while the median path grows 5.0%:
| ten weeks, $\lambda = 1$ | $I_{10}/I_0$ |
|---|---|
| mean, Monte Carlo over 200,000 regime paths | 1.2686 ± 0.0019 |
| mean, combined system | 1.2691 |
| median, Monte Carlo | 1.0497 |
| typical path, $e^{10(\bar\beta - \gamma)}$ | 1.0513 |
Several infected compartments: the first-order rule
With a latent stage the infected state is a vector, and the direction of that vector matters. Write the regime matrix as its average plus fluctuations:
$$A(y) = \bar A + \sum_k \phi_k(y)\,B_k .$$Here $\phi_k(y)$ is the deviation of the $k$-th switched rate from its $\pi$-average, and $B_k$ is the matrix it multiplies. For SEIR, transmission multiplies $B_\beta$ and recovery multiplies $B_\gamma$:
$$B_\beta = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}, \qquad B_\gamma = \begin{pmatrix} 0 & 0 \\ 0 & -1 \end{pmatrix}.$$Let $\bar x$ and $\bar y$ be the right and left principal eigenvectors of $\bar A$, scaled so that $\bar y^\top \bar x = 1$. The right vector is the stable mix of latent and infectious cases. The left vector is the reproductive value of each stage: $\bar y^\top x$ grows at exactly $r$ under the averaged model.
The Green–Kubo matrix is
$$K_{lk} = \int_0^\infty \operatorname{Cov}\big(\phi_l(y_0),\,\phi_k(y_t)\big)\,dt = -\pi\cdot\big(\phi_l\,Q^{\#}\phi_k\big).$$The first index is the earlier parameter. The matrix need not be symmetric.
The mean
The first-order rule adds $\sum K_{lk}$ times a product of operators to the averaged generator. The population runs forward in time, so the operator of the later parameter stands on the left, as in a product of propagators:
$$\mathbb E[x_t] \approx e^{t(\bar A + M)}\,x_0, \qquad M = \sum_{k,l} K_{lk}\,B_k B_l .$$The principal eigenvalue of $\bar A + M$ moves, to first order, by $\bar y^\top M \bar x$:
$$\Lambda_{\mathrm{mean}} \approx r + \sum_{k,l} K_{lk}\;\bar y^\top B_k B_l\,\bar x .$$The typical path
Follow the size $z = \bar y^\top x$ and the shape $v = x/z$, so that $\bar y^\top v = 1$. Differentiating,
$$(\log z)' = \bar y^\top A(y_t)\,v, \qquad v' = A(y_t)\,v - v\;\bar y^\top A(y_t)\,v .$$Write $v = \bar x + w$ with $\bar y^\top w = 0$. Since $\bar y^\top \bar A = r\,\bar y^\top$, the size grows at
$$(\log z)' = \rho(y_t) + \sum_k \phi_k(y_t)\;\bar y^\top B_k\,w_t, \qquad \rho(y) = \bar y^\top A(y)\,\bar x = r + \sum_k \phi_k(y)\;\bar y^\top B_k\,\bar x .$$The first term, the projected rate, averages to $r$ because each $\phi_k$ averages to zero. The second term averages to something non-zero only if the shape $w$ remembers recent fluctuations.
It does. To leading order the shape obeys
$$w' = (\bar A - r)\,w + \sum_l \phi_l(y_t)\,P B_l\,\bar x, \qquad P = I - \bar x\,\bar y^\top .$$Each fluctuation pushes the shape along $PB_l\bar x$, the part of $B_l\bar x$ that is not a change of size. Over one correlation time of the chain, $\bar A - r$ does little, so $w_t$ is the accumulated push of recent fluctuations. Its correlation with the current $\phi_k$ is
$$\mathbb E\big[\phi_k(y_t)\,w_t\big] \approx \sum_l \Big(\int_0^\infty \operatorname{Cov}\big(\phi_l(y_0),\,\phi_k(y_s)\big)\,ds\Big)\,P B_l\,\bar x = \sum_l K_{lk}\,P B_l\,\bar x .$$Averaging the growth of the size gives the almost-sure rate:
$$\Lambda_{\mathrm{as}} \approx r + \sum_{k,l} K_{lk}\;\bar y^\top B_k\,P\,B_l\,\bar x .$$This is the fast-switching formula of Monmarché, Schreiber and Strickler (2025) for Markov switching, written with the Green–Kubo matrix. Their formula uses the group inverse of the jump chain, from Monmarché and Strickler (2025). The derivation above is a heuristic for it.
The gap
The two formulas differ only by the projection $P$. Since $B_l\bar x = PB_l\bar x + \bar x\,(\bar y^\top B_l \bar x)$, the difference is the part of each fluctuation that changes the size of the epidemic without changing its shape.
The typical path takes the log of the size, and those parts average out. The mean takes the average of the size, and they compound.
The gap is a structured Lewontin–Cohen term. The number $\bar y^\top B_k \bar x$ is the derivative of $r$ with respect to the $k$-th rate, so $\rho(y)$ is the growth rate linearized regime by regime.
The gap is the Green–Kubo integral of that rate, and it sees only the symmetric part of $K$. The antisymmetric part moves both growth rates by the same amount.
With one infected compartment $P = 0$. The almost-sure rate then has no correction and the mean gains $K_{\beta\beta}$, as in the SIR section.
SEIR with switched transmission
Only $\beta$ switches. Its matrix is nilpotent:
$$B_\beta^2 = 0 .$$Transmission moves infection from $I$ into $E$, and $E$ is not transmitting. A burst of transmission cannot compound within one correlation time, because the new cases are still latent. So the mean has no first-order correction.
The almost-sure rate still has one. Since $\bar y^\top B_\beta P B_\beta \bar x = -(\bar y^\top B_\beta \bar x)^2$, it falls by $K_{\beta\beta}(\bar y_E\,\bar x_I)^2$. For SEIR the principal vectors are explicit:
$$\bar x = \begin{pmatrix} r + \gamma \\ \sigma \end{pmatrix}, \qquad \bar y = \frac{1}{\sigma\sqrt D}\begin{pmatrix} \sigma \\ r + \sigma \end{pmatrix}, \qquad \bar y_E\,\bar x_I = \frac{\sigma}{\sqrt D} = \frac{\partial r}{\partial\bar\beta}.$$For a symmetric two-state chain the mean has an exact equation too. Write $s = m_1 + m_2$ and $d = m_1 - m_2$. Then
$$s' = \bar A\,s + \tilde\beta\,B_\beta\,d, \qquad d' = (\bar A - 2\lambda)\,d + \tilde\beta\,B_\beta\,s .$$Eliminating $d$ from the eigenvalue problem leaves the averaged matrix with $\bar\beta$ raised by $\tilde\beta^2$ times the $(I, E)$ entry of $(\Lambda_{\mathrm{mean}} + 2\lambda - \bar A)^{-1}$. That gives a scalar equation for the mean growth rate, solved by iteration, and its expansion.
An exact typical-growth formula for two regimes
The first-order formula is not needed when only two transmission rates switch. Put $z=E/I$. In regime $i$,
$$z' = F_i(z) = \beta_i + (\gamma-\sigma)z-\sigma z^2 = -\sigma(z-a_i)(z-b_i), \qquad b_i<0<a_i.$$Suppose $\beta_h>\beta_\ell>0$, with switching rates $q_{h\ell}$ and $q_{\ell h}$. Then $a_\ell<z<a_h$ is invariant and attracting: the high-transmission flow points right and the low-transmission flow points left. Write
$$d_i=\sqrt{(\gamma-\sigma)^2+4\sigma\beta_i}, \qquad \kappa_h=\frac{q_{h\ell}}{d_h}, \qquad \kappa_\ell=\frac{q_{\ell h}}{d_\ell},$$ $$H(z)=\left(\frac{a_h-z}{z-b_h}\right)^{\kappa_h} \left(\frac{z-a_\ell}{z-b_\ell}\right)^{\kappa_\ell}.$$The invariant densities of the joint regime–ratio process are exactly
where $C>0$ makes the two densities integrate to one in total.
Indeed, the stationary forward equations and zero boundary flux give $F_h p_h+F_\ell p_\ell=0$. If $u=F_h p_h$, the remaining equation is $u'=-u(q_{h\ell}/F_h+q_{\ell h}/F_\ell)$; integration gives $H$. The endpoint powers $\kappa_h,\kappa_\ell>0$ make the densities integrable, and their regime marginals are necessarily the stationary probabilities of the two-state chain.
Finally $(\log I)'=\sigma z-\gamma$, so the ergodic average gives the displayed Lyapunov exponent. This proves the formula for every positive initial infected vector.
The slow and fast switching endpoints
The same formula also resolves the opposite, slow-switching limit. Scale the transition rates as $q_{h\ell}=c\alpha$ and $q_{\ell h}=c\delta$, with fixed $\alpha,\delta>0$, and let $c\downarrow0$.
Put $r_i=\sigma a_i-\gamma$ and $\pi_h=\delta/(\alpha+\delta)$, $\pi_\ell=\alpha/(\alpha+\delta)$. Then
Here is a direct renewal proof. A high-to-low switch occurs at stationary flux $c\alpha\delta/(\alpha+\delta)$ and contributes the deterministic relaxation area
$$\int_0^\infty\{z_\ell(t;a_h)-a_\ell\}\,dt =\frac1\sigma\log\frac{a_h-b_\ell}{a_\ell-b_\ell}.$$A low-to-high switch contributes
$$\int_0^\infty\{z_h(t;a_\ell)-a_h\}\,dt =-\frac1\sigma\log\frac{a_h-b_h}{a_\ell-b_h}.$$Multiplying their sum by the switch flux and by the outer factor $\sigma$ gives $cJ$. The projective flows approach their attracting roots exponentially.
Consequently the expected displacement of the entry point from the previous root is $O(c)$; multiplying its bounded relaxation-area effect by the $O(c)$ switch flux gives the stated $O(c^2)$ remainder. Positivity is exact: the numerator minus denominator inside the logarithm is $(a_h-a_\ell)(b_\ell-b_h)>0$.
At the other endpoint, with $\bar\beta=\pi_h\beta_h+\pi_\ell\beta_\ell$ and $\bar D=(\gamma-\sigma)^2+4\sigma\bar\beta$, the Green–Kubo expansion gives
Since $r(\beta)$ is strictly concave, the fast limit $r(\bar\beta)$ strictly exceeds the slow limit $\pi_h r_h+\pi_\ell r_\ell$ when the transmission rates differ. Here the endpoint signs are part of a stronger global result.
Strict global monotonicity by convex order
For this two-state, transmission-only model the almost-sure exponent is strictly increasing in the common switching-speed scale $c$. To see why, remove the constant damping from the infectious component and use the Riccati coordinate
$$y=\sigma z+\frac{\sigma-\gamma}{2},\qquad \dot y=A_i-y^2,\qquad A_i=\frac{(\sigma-\gamma)^2}{4}+\sigma\beta_i.$$Put $\rho_i=\sqrt{A_i}$ and $\bar A=\pi_hA_h+\pi_\ell A_\ell$. The total invariant density of $y$ on $(\rho_\ell,\rho_h)$ can be written
Let $X_c=y^2$. Stationarity applied to $y$ gives the speed-independent mean
$$\mathbb E X_c=\mathbb E A_{Y}=\bar A.$$If $\psi(x)=\phi(\sqrt x)$, direct differentiation gives
$$\psi'(x)=\frac{(\alpha+\delta)(\bar A-x)} {2\sqrt x\,(A_h-x)(x-A_\ell)}.$$Thus $\psi$ rises strictly to its unique maximum at the common mean $\bar A$ and then falls strictly. For $c_2>c_1$, the likelihood ratio of $X_{c_2}$ to $X_{c_1}$ is proportional to $e^{(c_2-c_1)\psi(x)}$: it is larger on one central interval and smaller in both tails. Because the two laws have the same mean, their distribution functions cross once, so $X_{c_2}\preceq_{\rm cx}X_{c_1}$.
This is the classical convex-order crossing argument of Ohlin (1969). Since $\sqrt{x}$ is strictly concave,
$$\mathbb E\sqrt{X_{c_2}}>\mathbb E\sqrt{X_{c_1}},\qquad \Lambda_{\rm as}(c)=\mathbb E\sqrt{X_c}-\frac{\sigma+\gamma}{2}.$$This proves strict global monotonicity whenever $\beta_h\ne\beta_\ell$. It is special to the ordered two-environment model in which only transmission changes; it does not assert monotonicity when recovery also switches or for a general family of cooperative matrices.
The same speed-invariant mean gives a quantitative bound at every finite switching speed, rather than only an endpoint expansion. Let $\Lambda_\infty=\sqrt{\bar A}-(\sigma+\gamma)/2$. The tangent remainder of the square root is the exact identity
The equality follows from $\sqrt{\bar A}+(x-\bar A)/(2\sqrt{\bar A})-\sqrt{x} =(x-\bar A)^2/[2\sqrt{\bar A}(\sqrt{x}+\sqrt{\bar A})^2]$ after taking expectations. The inequalities use only $A_\ell<X_c<A_h$, so they are non-asymptotic and uniform over every $c>0$.
Convex-order concentration also makes $\operatorname{Var}(X_c)$ strictly decrease with speed. These bounds quantify convergence to the fast limit; they do not supply a rate in $c$ until the variance itself is estimated.
After $z=a_\ell+(a_h-a_\ell)t$, the two integrals have Jacobi weights $t^{\kappa_\ell}(1-t)^{\kappa_h-1}$ and $t^{\kappa_\ell-1}(1-t)^{\kappa_h}$. The certificate therefore evaluates them by Gauss–Jacobi quadrature, without time discretization or endpoint clipping.
For the parameters here and symmetric switching, the exact extinction–invasion threshold is
$$\lambda_* = 2.3126879896\ \text{per week}.$$The slow limit is negative and the fast limit positive, so strict monotonicity proves that this extinction–invasion crossing is unique. The first-order formula predicts $2.6316808690$, about $13.8\%$ too high. This exact reduction is special to two regimes with only $\beta$ switched; with three or more regimes, or with switched recovery, the projective invariant law is no longer this one-dimensional two-flow problem.
For any chain the almost-sure formula holds with $K_{\beta\beta} = -\pi\cdot(\tilde\beta\,Q^{\#}\tilde\beta)$, and the mean still has no first-order term.
The latent stage reverses the SIR picture. In SIR the fluctuation raises the mean and leaves the typical path alone.
In SEIR it leaves the mean nearly alone and lowers the typical path. With $\bar\beta/\gamma = 1.02$ and regimes lasting half a week, the averaged model and the mean both grow, yet almost every introduction dies out.
SEIR with switched transmission and recovery: the order of regimes
Now both $\beta$ and $\gamma$ switch, for instance a testing campaign that shortens the infectious period. The products of the two matrices are
$$B_\beta B_\gamma = -B_\beta, \qquad B_\gamma B_\beta = 0, \qquad B_\gamma^2 = \begin{pmatrix} 0 & 0 \\ 0 & 1 \end{pmatrix}.$$So the mean picks up $K_{\gamma\gamma}$ and $K_{\gamma\beta}$, recovery followed by transmission, but not $K_{\beta\gamma}$. Fast recovery removes the infectious who would transmit next. Transmission followed by fast recovery does nothing at first order, because the new cases are latent while recovery is fast.
That asymmetry makes the direction of a cycle matter. Take three regimes visited in the order $1 \to 2 \to 3 \to 1$ at rate $a$, with reverse steps at rate $b$.
The stationary law is uniform. The Green–Kubo integral of two rates $f$ and $g$ on this cycle is a short calculation with the discrete Fourier transform:
$$K_{fg} = \frac{3(a + b)\,\tilde f\cdot\tilde g + (a - b)\,W(f, g)}{18\,(a^2 + ab + b^2)}, \qquad W(f, g) = \sum_{j} \big(f_j\,g_{j+1} - f_{j+1}\,g_j\big),$$with indices taken mod 3. The first term is symmetric.
The second is antisymmetric: $W(\beta, \gamma)$ is twice the signed area of the triangle traced by the points $(\beta_j, \gamma_j)$ in the order the cycle visits them. Reversing the cycle swaps $a$ and $b$ and flips it.
Putting these into the general rule, with $\bar y^\top B_\beta \bar x = \sigma/\sqrt D$ and $\bar y^\top B_\gamma \bar x = -(r + \sigma)/\sqrt D$:
The averages, and the symmetric part of $K$, are the same in both directions. Only the order differs, and it moves the mean and the typical growth rate by the same amount.
Parameters: regimes open, testing and cautious, with $\beta = (0.605,\ 0.105,\ 0.055)$ and $\gamma = (0.1,\ 0.55,\ 0.1)$ per week; $\sigma = 0.5$; so $\bar\beta = 0.255$, $\bar\gamma = 0.25$ and $\bar\beta/\bar\gamma = 1.02$ as before. The forward cycle runs open, testing, cautious at rate $a = 2c$ with reverse steps at $b = 0.1c$; the reversed cycle swaps $a$ and $b$.
At $c = 2$ each regime lasts a little under two days and a lap of the cycle takes about five and a half days. In the forward order the testing push follows the open phase and finds its new infections still latent. In the reversed order it comes just before the open phase and removes the people who would have transmitted in it.
A worked example
Transmission switched on two states at $\lambda = 2$, so that each regime lasts half a week:
| quantity | value |
|---|---|
| $\sqrt D$ | 0.756637 |
| $\bar y_E\,\bar x_I = \sigma/\sqrt D$ | 0.660819 |
| $K_{\beta\beta} = \tilde\beta^2/(2\lambda)$ | 0.010000 |
| averaged SEIR growth $r$ | 0.003319 |
| SEIR mean, second order | 0.004145 |
| SEIR mean, numerical | 0.004012 |
| SEIR almost-sure, first order $r - K_{\beta\beta}\sigma^2/D$ | −0.001048 |
| SEIR almost-sure, exact projective density | −0.000448 |
| SEIR almost-sure, Monte Carlo | −0.000450 ± 3e−6 |
| SIR mean, first order $\bar\beta - \gamma + K_{\beta\beta}$ | 0.015000 |
| SIR almost-sure $\bar\beta - \gamma$ | 0.005000 |
The averaged model says the epidemic grows. So does the mean.
The typical path declines, and the exact value is −0.0004478 per week, while Monte Carlo gives −0.000450. The first-order error falls by four for each doubling of $\lambda$, as the checks show.
The three-regime cycle at $c = 2$:
| quantity | forward | reversed |
|---|---|---|
| $a$, $b$ | 4, 0.2 | 0.2, 4 |
| $W(\beta, \gamma)$ | 0.2475 | −0.2475 |
| $K_{\beta\beta}$, $K_{\gamma\gamma}$ | 0.007690, 0.005612 | 0.007690, 0.005612 |
| $K_{\beta\gamma}$, $K_{\gamma\beta}$ | 0.000297, −0.005909 | −0.005909, 0.000297 |
| averaged growth $r$ | 0.003319 | 0.003319 |
| mean, first order | 0.010956 | 0.006855 |
| mean, numerical | 0.010650 | 0.007269 |
| gap $\Lambda_{\mathrm{mean}} - \Lambda_{\mathrm{as}}$, first order | 0.008308 | 0.008308 |
| almost-sure, first order | 0.002648 | −0.001453 |
| almost-sure, Monte Carlo | 0.002592 ± 3e−6 | −0.000903 ± 3e−6 |
The averaged rates are the same and $\bar\beta/\bar\gamma = 1.02$ in both cases. In the forward cycle a typical introduction takes off.
In the reversed cycle it dies out. The mean grows in both.
In the forward cycle the mean grows by 1.07% a week and doubles in about 65 weeks. The typical path grows by 0.26% a week and doubles in about five years.
The mean is pulled up by the few regime sequences in which open phases keep landing at the right moments. Growth compounds, so those paths end up far larger than the rest.
Checks
The mean growth rate is checked against the principal eigenvalue of the combined matrix $Q^\top \otimes I + \operatorname{blockdiag}(A_i)$, computed numerically. The almost-sure rate is checked by Monte Carlo of the switched system with no time discretization: each holding interval of the chain is propagated with the matrix exponential of its regime.
Each estimate uses 2,000 paths of 4,000 to 40,000 regime changes. The control variate $\int_0^t \rho(y_s)\,ds$, whose time average is exactly $r$, is subtracted from $\log \bar y^\top x_t$ to reduce the variance.
That the combined matrix gives the mean is checked directly for SIR at $\lambda = 1$ over ten weeks: 200,000 regime paths give $\mathbb E[I_{10}]/I_0 = 1.2686 \pm 0.0019$, against 1.2691 from the matrix exponential.
SIR with switched transmission. The quadratic root agrees with the numerical eigenvalue to $10^{-12}$. For the symmetric chain the first-order error is third order, falling by eight per doubling:
| $\lambda$ | $K_{\beta\beta}$ | mean, numerical | mean, first order | error | almost-sure |
|---|---|---|---|---|---|
| 1 | 0.02000 | 0.024804 | 0.025000 | 2.0e−4 | 0.005000 |
| 2 | 0.01000 | 0.014975 | 0.015000 | 2.5e−5 | 0.005000 |
| 4 | 0.00500 | 0.009997 | 0.010000 | 3.1e−6 | 0.005000 |
| 8 | 0.00250 | 0.007500 | 0.007500 | 3.9e−7 | 0.005000 |
| 16 | 0.00125 | 0.006250 | 0.006250 | 4.9e−8 | 0.005000 |
SEIR with switched transmission. The scalar mean equation agrees with the numerical eigenvalue to $10^{-11}$. The exact projective-density exponent lies within four Monte Carlo standard errors in every row; doubling the Gauss–Jacobi rule from 96 to 192 nodes changes the value at $\lambda=2$ by less than $2\times10^{-14}$.
The mean moves at second order, and $(\Lambda_{\mathrm{mean}} - r)\lambda^2$ tends to $\tilde\beta^2\sigma^2/(4\sqrt D) = 0.00330$. The almost-sure error falls by four per doubling:
| $\lambda$ | mean, numerical | $(\Lambda_{\mathrm{mean}} - r)\lambda^2$ | almost-sure, exact | Monte Carlo | first order | first-order error |
|---|---|---|---|---|---|---|
| 1 | 0.005703 | 0.00238 | −0.003310 | −0.003307 | −0.005415 | −2.1e−3 |
| 2 | 0.004012 | 0.00278 | −0.000448 | −0.000450 | −0.001048 | −6.0e−4 |
| 4 | 0.003507 | 0.00302 | 0.001297 | 0.001301 | 0.001135 | −1.6e−4 |
| 8 | 0.003368 | 0.00315 | 0.002269 | 0.002268 | 0.002227 | −4.2e−5 |
| 16 | 0.003331 | 0.00323 | 0.002784 | 0.002783 | 0.002773 | −1.1e−5 |
Across 29 switching speeds from $0.01$ to $64$, for both symmetric rates and the asymmetric ratio $q_{h\ell}:q_{\ell h}=1.7:0.4$, the certificate checks the exact weighted-variance identity and its bounds, and verifies that $\operatorname{Var}(X_c)$ strictly decreases. The lower bound is at least $53.9\%$ of the exact gap, while the upper bound is at most $1.933$ times the gap over this grid.
The three-regime cycle. The closed-form $K$ agrees with $-\pi\cdot(f\,Q^{\#}g)$ to $10^{-13}$. Both first-order formulas have errors falling at second order in each direction:
| $c$ | forward mean, numerical | first order | reversed mean, numerical | first order |
|---|---|---|---|---|
| 1 | 0.017530 | 0.018593 | 0.011917 | 0.010392 |
| 2 | 0.010650 | 0.010956 | 0.007269 | 0.006855 |
| 4 | 0.007055 | 0.007137 | 0.005195 | 0.005087 |
| 8 | 0.005207 | 0.005228 | 0.004230 | 0.004203 |
| 16 | 0.004268 | 0.004273 | 0.003768 | 0.003761 |
| $c$ | forward almost-sure, Monte Carlo | first order | reversed almost-sure, Monte Carlo | first order |
|---|---|---|---|---|
| 1 | 0.001816 | 0.001977 | −0.004200 | −0.006224 |
| 2 | 0.002592 | 0.002648 | −0.000903 | −0.001453 |
| 4 | 0.002968 | 0.002983 | 0.001078 | 0.000933 |
| 8 | 0.003147 | 0.003151 | 0.002162 | 0.002126 |
| 16 | 0.003234 | 0.003235 | 0.002731 | 0.002722 |
Monte Carlo standard errors are below $6\times10^{-6}$ in every row.
The forward-minus-reversed difference against the prediction $(a - b)\,\sigma W/(9(a^2 + ab + b^2)\sqrt D)$, with $W = 0.2475$:
| $c$ | mean, numerical | almost-sure, Monte Carlo | predicted | numerical / predicted |
|---|---|---|---|---|
| 1 | 0.005613 | 0.006017 ± 8e−6 | 0.008201 | 0.684 |
| 2 | 0.003382 | 0.003495 ± 5e−6 | 0.004101 | 0.825 |
| 4 | 0.001860 | 0.001890 ± 3e−6 | 0.002050 | 0.907 |
| 8 | 0.000976 | 0.000985 ± 1e−6 | 0.001025 | 0.952 |
| 16 | 0.000500 | 0.000503 ± 9e−7 | 0.000513 | 0.976 |
The ratio tends to one, and the error of the prediction falls by four per doubling of $c$. The almost-sure difference follows the mean difference, since the gap between the two rates sees only the symmetric part of $K$.
The general formulas on a three-compartment model with two switched rates and the irreversible three-state generator $Q_A$ of the any-equation page sped up by a factor $c$. The gap equals the Green–Kubo integral of $\rho$ to machine precision at every speed:
| $c$ | mean, numerical | first order | almost-sure, Monte Carlo | first order | gap, first order |
|---|---|---|---|---|---|
| 2 | 0.053513 | 0.053229 | 0.051206 ± 3e−6 | 0.050775 | 0.002455 |
| 4 | 0.053906 | 0.053832 | 0.052717 ± 1e−6 | 0.052604 | 0.001227 |
| 8 | 0.054152 | 0.054133 | 0.053550 ± 8e−7 | 0.053519 | 0.000614 |
| 16 | 0.054288 | 0.054283 | 0.053984 ± 5e−7 | 0.053976 | 0.000307 |
Where epidemics turn over, with the immunity factor measured from contact surveys, US county waves and serology in six places, is developed at epidemiology.microprediction.org.
Prior art
The gap between the growth seen early in an epidemic and the threshold at which it turns over was the subject of Cotton (2020), first set out in May 2020. In a mixture of independent sub-populations the mean early growth is raised by a growth convexity, computed from the moment generating function of the transmission rate, or from Vasicek's bond formula when that rate follows an Ornstein–Uhlenbeck process.
The susceptible fraction at the peak is raised by a harmonic-mean convexity. Together they explain an epidemic that turns over well below the herd-immunity threshold its early growth implies.
The growth convexity is the mean growth rate above. For a fast chain it is the averaged rate plus $K_{\beta\beta}$.
Cotton (2020) attenuates transmission by the probability that a contact is new, so that how fast acquaintances go stale changes the final size by an order of magnitude.
Britton, Ball and Trapman (2020) and Gomes et al. (2022) lower the threshold with persistent heterogeneity in activity or susceptibility. Tkachenko et al. (2021) let individual activity change over time and find a transient collective immunity below the long-run threshold, with an immunity factor near 4 in New York City and Chicago.
The distinction between the two growth rates is old. Lewontin and Cohen (1969) showed that a population can grow in mean while almost every path goes extinct. Tuljapurkar (1982) gave small-noise approximations of the stochastic growth rate for structured populations in correlated environments, and Arnold (1984) connected moment and sample stability of linear stochastic systems.
For epidemics, Gray, Greenhalgh, Mao and Pan (2012) treat SIS with telegraph-noise switching of the rates. They give a reproduction number for the mean and prove that almost-sure extinction is decided by a different threshold, the $\pi$-average of the transmission rate against the $\pi$-average of the removal rate. That is the scalar statement $\Lambda_{\mathrm{as}} = \bar\beta - \gamma$ above.
Li, Liu and Cui (2017) extend it to SIRS. Bacaër and Khaladi (2013) define a reproduction number in a random environment that decides growth or decay of typical populations, and contrast it with the parameter of Hernandez-Suarez, Rabinovich and Hernandez (2012), which decides growth of the expectation. Prodhomme and Strickler (2024) show that the sign of the top Lyapunov exponent of the linearized switched system decides the fate of a multitype SIS epidemic.
With white noise in place of switching, Gray, Greenhalgh, Hu, Mao and Pan (2011) find that white noise of intensity $\sigma_W$ in the transmission term, in Itô form, lowers the almost-sure threshold. In the linear phase it lowers the typical growth rate by $\sigma_W^2/2$ and leaves the mean unchanged. A fast telegraph rate with $2K_{\beta\beta} = \sigma_W^2$ converges to Stratonovich noise (Wong and Zakai 1965).
In both models the mean exceeds the typical growth by $K_{\beta\beta}$. They differ only in which of the two carries the average rate.
For fast switching of a general linear system $x' = A(y_t)x$ with Metzler matrices, Monmarché, Schreiber and Strickler (2025) give the first-order correction to the top Lyapunov exponent, for periodic and for Markov switching, building on Monmarché and Strickler (2025). The almost-sure formula above is their Markov result, rewritten with the Green–Kubo matrix and the projection $P$. For periodic switching they find that the first-order term is a commutator, which vanishes for two environments.
Their Markov formula contains the full non-symmetric group inverse, and so the antisymmetric part of $K$, without separating it. Benaïm, Le Borgne, Malrieu and Zitt (2014) show that switching between stable planar systems can be unstable.
For periodic environments, Wang and Zhao (2008) define the reproduction number of compartmental models and prove it is the invasion threshold. That threshold is the periodic, deterministic counterpart of the invasion question here.
Not found in the literature reviewed: the first-order mean growth rate in the Green–Kubo form, the gap between the two rates as the Green–Kubo integral of the projected rate, the absence of a first-order mean correction in SEIR when only transmission switches, the explicit two-regime SEIR projective density and exact threshold above, and the closed-form effect of the direction of a regime cycle through $W(\beta, \gamma)$.
All of this is the linear phase. Depletion of susceptibles, final size and endemic levels are outside it.
Certificate: verify_epidemics.py.