Click here for the IPython Notebook used to generate this post
Revised on 2023-07-05
Since late 2020, the dynamics of SARS-CoV-2 incidence is driven by the emergence of novel variants that are either more transmissible than competing variants, or have changed in such a way that they can (partially) evade population immunity. Immune evasive variants can spread faster due to a larger pool of susceptible hosts. The spread of such contagious diseases and the competition between variants is typically modeled by Susceptible-Infected-Recovered (SIR) models. These SIR models come in many flavors, but even more complex versions are rather crude approximations of the spread of a diverse virus in a population that is heterogenous in many aspects, including immunological, demographic, mobility, contact patters, etc. As a result these models are of limited use for quantitative predictions and inference, but they can inform intuition and illuminate how observed patterns depend on phenomenological parameters that capture general aspects of spread. The latter is the objective of this note and it is mostly reiterating well known material from textbooks. But I find myself having to look them up now and then and it is useful for me to have them all in the same place. In particular, I will discuss
- how the final epidemic size in absence of waning depends on model parameters
- how the fixed point, frequency of waves, and dampening of oscillations depend on the rate of waning, and
- the size of variant waves and how they depend on population heterogeneity.
The last point arose in frequently when we observed a novel variant growing in proportion at a certain rate: How many people will get infected in such a variant driven wave. It turns out that the size of a wave depends strongly on population heterogeneity.
I will close with some thoughts on how to incorporate gradual waning of immunity in the population, for example via gradually decreasing neutralization titers, into SIR models. Waves generated by immune evasive in these models fall between the extreme cases expected from simple SIR models without population heterogeneity.
Basic SIR model with waning
The basic SIR model is given by $$ \dot{S} = -\beta I S + \gamma (1-S) $$ $$ \dot{I} = \beta I S - \nu I $$
where \(\nu\) is the recovery rate of an infected individual. We will use this rate to set the unit of time, that is we measure time in units of pathogen generations and set to \(\nu=1\). For diseases of interest here, \(\nu \approx 1/\mathrm{week}\) or \(2/\mathrm{week}\) such that there are 50 to 100 generations per year. The parameters \(\beta\) and \(\gamma\) are the rate at which the disease spreads from infected to susceptible individuals and the rate at which recovered individuals loose their immunity and become susceptible again, respectively. The basic reproductive number is given by \(R_0 = \beta/\nu\), while the initial growth rate of an outbreak in a completely susceptible population is given by \(\beta - \nu\). We ignore population turn-over in this note.
The behavior of this model with and without waning (\(\gamma=0\) or \(\gamma>0\)) is shown below. The model is solved numerically, see the notebook for the code.
As expected, the numerical solution above shows that the initial dynamics do not depend on \(\gamma\). This holds as long as waning is slow compared to the duration of the first wave. After the initial peak, the solution with waning exhibits a damped oscillation around a fixed point. Without waning, the pool of susceptibles is irreversibly reduced and so small that it does not support pathogen spread. The system hence settles to a disease free state \(I(t)=0\) with 94% of the population immune (in the above example). Both the final epidemic size in absence of waning and the equilibrium properties with waning can be explored analytically.
Final epidemic size
Without waning, the residual susceptibility is simply the fraction of people that didn't get infected during the wave. To solve for the final size (See Miller (2012) for a detailed discussion) for \(\gamma=0\), one divides \(\dot{I}\) by \(\dot{S}\) to obtain $$ \frac{dI}{dS} = -1 + \frac{1}{\beta S} $$ which can be integrated $$ I(t) - I(0) = S(0) - S(t) + \frac{\log(S(0)/S(t))}{\beta} $$ For \(t\to \infty\), \(I(t) \to 0\) and \(I(0)\) can be chosen infinitesimally small such that the final susceptibility is the solution of this transcendental equation: $$ S(t_\infty) = S(0)e^{-\beta (S(0) - S(t_\infty))} $$
For large \(\beta\) and an initially completely susceptible population \(S(0)=1\), \(S(\infty)\approx e^{-\beta}\). For \(\beta\) close to 1 we can expand the exponential and solve for \(R(t_\infty) = 1-S(t_\infty)\) $$ \begin{split} R(t_\infty) & = 1 - e^{-\beta R(t_\infty)} \approx \beta R(t_\infty) - \frac{\beta^2 R(t_\infty)^2}{2} + \cdots \end{split} $$ This linear approximation yields \(R(t_\infty) \approx 2(\beta-1)\) confirming the notion that a similar of people are infected before and after the peak of the wave. The full equation can be solved iteratively.
beta = 1.4
def final_epidemic_size(beta, max_iter=None, S0=1):
if max_iter is None: max_iter = 10 + int(10/(beta-1))
Sf = 0.5*S0
for i in range(max_iter):
Sf = S0*np.exp(-beta*(S0-Sf))
return 1-Sf
beta = 3
print(f"final size for {beta=} is: {final_epidemic_size(beta, S0=1.0):1.3f}")
final size for beta=3 is: 0.940
Equilibrium with waning
With positive \(\gamma\), there is a non-trivial fix point at $$ \bar{S} = \beta^{-1} \quad \mathrm{and} \quad \bar{I} = \gamma(1-\beta^{-1}) = \gamma(1-\bar{S}) $$ This fixed point has the intuitive property that average number of secondary cases \(R_e = \beta\bar{S} = 1\) and that the average incidence is rate at which immune individuals become susceptible (the product of the rate of waning and the size of the immune compartment \(1-\bar{S}\)). Deviations \(\delta S\) and \(\delta I\) from this steady state (\(\bar{S}, \bar{I}\)) can be analyzed by linear perturbation theory. In matrix form, the linearized system is given by $$ \begin{pmatrix} \dot{\delta S}\\ \dot{\delta I} \end{pmatrix} = \begin{pmatrix} -\gamma \beta & -1\\ \gamma (\beta-1) & 0 \end{pmatrix} \begin{pmatrix} \delta S\\ \delta I \end{pmatrix} $$ with characteristic equation \(-\lambda(\gamma\beta - \lambda) + \gamma (\beta-1)=0\) and eigenvalues $$ \lambda_{1/2} = -\frac{\gamma \beta}{2} \pm \sqrt{\gamma^2 \beta^2/4 - \gamma(\beta -1)} $$ The corresponding right-eigenvectors are $$ v_{1/2} = \begin{pmatrix} -1 \\ \lambda_{1/2} + \gamma \beta \end{pmatrix} $$ When \(\gamma < \frac{4}{\beta(1- \beta)}\), as is true in cases of interest, the system exhibits dampened oscillations with a frequency $$ \omega = \sqrt{\gamma(\beta -1) - \gamma^2 \beta^2/4} \approx \sqrt{\gamma(\beta-1)}\left(1 - \frac{\gamma\beta^2}{8(\beta-1)}+\ldots \right) $$ and a decay time \(\tau = \frac{2}{\gamma\beta}\).
A period \(T\) is thus related to \(\gamma\) and \(\beta\) via \(2\pi = \omega T \approx \sqrt{\gamma(\beta-1)}T\) such that the rate of waning is \(\gamma = \frac{4\pi^2}{T^2(\beta-1)}\). With \(\beta\) somewhere around 5 and waves every 4 months (20-25 serial intervals), we obtain \(\gamma = \frac{4\pi^2}{4\cdot{} 25^2}\approx \frac{1}{60}\) per serial interval. Note that this rate of waning includes contributions from viral evolution, which will be discontinuous jumps associated with the emergence of novel variants, and individual loss of immunity which will be gradual at the population level.
Linear solution around the fixed point
To calculate the coefficients of the linear solution with the above eigenvalues and eigenvectors, we need the left eigenvectors $$ u_{1/2} = \frac{1}{\gamma (\beta-1)}\left(\begin{matrix} \gamma (\beta-1) \ \lambda_{1/2} + \gamma \beta \end{matrix}\right) $$ which scaled such that \(v_i \cdot u_i = 1\). Coefficients of the solution $$ \begin{pmatrix} S\\ I \end{pmatrix} = \begin{pmatrix} \bar{S}\\ \bar{I} \end{pmatrix} + e^{-t/\tau} \left(a_1 e^{i\omega t} + a_2 e^{-i \omega t}\right) $$ where \(a_i = u_i \left( \delta S(0), \delta I(0) \right)^t \)
Size of waves caused by invading variants
After having explored the size of an initial outbreak and the equilibrium in a population where immunity wanes, let's consider a new, initially rare, variant that grows at a rate \(\alpha\) in a population where the resident variant is at equilibrium. This growth rate is related to \(\beta^{*} \) and the susceptibility \(S^{*} \) of the new variant via \(\alpha = \beta^{*} S^{*} - 1\).
Within this simple model, a novel variant could grow with rate \(\alpha\) for different reasons
- Transmissibility: The variant is more transmissible with increased \(\beta^{*}\) but unchanged \(S^{*}=S\), no immune evasion. In this case \(\beta^* = \frac{1+\alpha}{S} = (1 + \alpha)\beta\).
- Zero-One immune evasion: The variant grows because a fraction of the population is completely susceptible again but \(\beta^*=\beta\) is unchanged. In this case, the new susceptibility is \(S^* = S (1+\alpha)\).
- Uniform immune evasion: The variant grows because immune protection of the entire population is reduced. This would correspond to a model with \(S^* = 1\) and \(\beta^* = 1+\alpha\) (This ignores the previously susceptible fraction, but this will be small for a disease with large \(R_0\)).
These three scenarios (mixtures of them are of course possible) lead to waves of rather different size. Since we are concerned with the wave caused by a single invading variant, we can ignore waning and set \(\gamma=0\). Despite identical initial growth rates, the size of the wave is much higher when transmission is supported by slow spread in the entire population compared to fast spread in a small completely susceptible population. When significant waning occurs over the course of the wave caused by the new variant, the epidemic sizes of the different scenarios become more similar.
Over the past 18 months, we have frequently seen variants with doubling times between one and two weeks, corresponding to growth rates between 0.25 and 0.5 per serial interval. The above example is for \(\alpha=0.25\) per serial interval with attack rates varying from 10% to above 35%. For faster growing variants, e.g. \(\alpha=0.5\), the attack rates range from 17% for the zero-one immune evasion scenario to 58% for the uniform evasion scenario. In both cases, the size of resulting wave differs about 3-fold depending on population heterogeneity.
Incorporating heterogeneous population immunity in SIR models
The models above classify individuals as either completely susceptible or completely protected with transitions between these "compartments". More generally, different individuals differ in their response to the virus and have been exposed to different variants at different times in the past with gradual waning of immunity. We therefore expect there to be a distribution of protective immunity in the population. Related extensions of SIR models have been proposed by Mohamed El Khalifi and Tom Britton.
A simple model of the probability of infection upon exposure could be $$ p = e^{-x} $$ where \(x\) could be interpreted as a neutralization titer or some other correlate of protection. In the population, there is a distribution \(\phi(x,t)\) that is governed by waning with rate \(v\), boosting by infection with rate \(\nu I(t)\): $$ \frac{\partial \phi}{\partial t} = D \frac{\partial^2 \phi(x,t)}{\partial^2 x} + v\frac{\partial \phi(x,t)}{\partial x} + I(t) \left(\delta(x-x_0) - \beta \phi(x,t) e^{-x}\right) $$ The first diffusive term accounts for variation in the waning rate \(v\) and boosting through cross-immunity with other pathogens.
The infection dynamics is given by $$ \frac{dI}{dt} = I(t) \left(\beta \int_0^\infty dx \phi(x,t) e^{-x} - 1\right) $$
As before, we use the recovery rate to set the unit of time. Other scales and numbers that play a role are:
- \(\beta\) is the infectivity and related to \(R_0\) via \(R_0 = \frac{\beta}{\nu}\)
- \(v\) is the velocity at which titers drop. We generally expect this to be on the order of 0.01 to 0.1/per week corresponding to 0.5 to 5 antigenic units per year.
- \(x_0\) is the typical post-infection titer.
- \(D\) the titer diffusivity. This diffusivity will capture heterogeneity in the rate of titer drops and other variation. I don't have an expectation for its value. But I imagine something like \(0.02\) per serial interval is sensible, resulting in an 2 unit^2 variance after a year.
The average susceptibility in this model is then given by $$ \langle S \rangle = \int_0^{x_0} dx\, \phi(x) e^{-x} $$ And the resulting growth rate of the outbreak is $$ \beta \langle S \rangle - 1 = \beta \int_0^{x_0} dx\, \phi(x) e^{-x} - 1 $$
The simple discretized implementation of this model is given by:
def SusceptI(Y, t, beta, gamma, D, v, xgrid, dx):
'''
time derivative of the prevalence and the susceptibility distribution.
'''
# convert D and v to right/left jump rates. stuff all waning into the left jumps
r = D/2/dx**2 #+ v/2/dx
l = D/2/dx**2 + v/dx
I, phi = Y[0], Y[1:]
dY = np.zeros_like(Y)
dY[0] = I*(beta*dx*np.sum(phi*np.exp(-xgrid))-1)
dY[1] = -(r + gamma + beta*I*np.exp(-xgrid[0]))*phi[0] + gamma/dx + l*phi[1]
dY[-1] = -(l + gamma + beta*I*np.exp(-xgrid[-1]))*phi[-1] + I/dx + r*phi[-2]
dY[2:-1] = - (r + l + gamma + beta*I*np.exp(-xgrid[1:-1]))*phi[1:-1] + r*phi[0:-2] + l*phi[2:]
return dY
This discretized model can be readily solved numerically, see the figure below. At first, this model exhibits narrow waves with increasing frequency before settling into a pattern of damped oscillations similar to the compartmental model above (top left panel).
Initially, the is no population heterogeneity and the time between waves is given by the time it takes for population to loose immunity from \(x=x_0\) to \(x<1\), which is \(t=x_0/v\). These intervals are indicated by the gray vertical bars. Since the gradual transmission from an immune state to a susceptible state takes time, this model effectively implements a delay which could give rise to sustained oscillations. For parameter values investigate here, oscillations are damped.
The population distribution of immunity \(x\) is shown in the lower two panels. The lower left panel show the distribution at three early time points, the right panel the approximate equilibrium distribution at long times. Once the population distribution of \(x\) has equilibrated, the system has a tendency to oscillate around a fixed point similar to the SIR model. But these oscillations decay less rapidly than in the equivalent SIR model with \(\gamma=v/x_0\), which would correspond to the inverse of the time over which people loose their immunity. The period of the oscillations is also longer than the \(t = 2\pi/\sqrt{(\beta-1)\gamma}\) expected in the compartmental model.
This pattern of accelerating waves that gradually dampen is also observed in COVID-19, for example in the hospitalization data of the UK:
Figure: Hospitalization data from England visualized by OurWorldInData.
However, the first two peaks here correspond to initial wave spring 2020 and the Alpha wave and their timing is not determined by immune waning dynamics but the emergence of variants. Since late 2021, the more or less continuous circulation of Omicron variants with a 3-4 months period is more suggestive of a dynamics determined by immune waning and reinfection.
Immune evasive variants
Now consider an immune evasive variant invading. The susceptibility of the population to this variant can differ from that of the original variant in different ways.
- the titers in the population could be uniformly reduced by a factor, that is \(x \rightarrow \max(0, x-\Delta x)\) and \(\phi'(x) = \phi(x+\Delta x) + \phi_0 \delta(x)\), where \(\delta(x) \) is a Dirac-Delta function and \(\phi_0 = \int_0^{\Delta x} \phi(x)dx \). We'll refer to this case as shifted susceptibility.
- a fraction \(\epsilon\) of the population becomes completely susceptible, i.e. \(\phi'(x) = \epsilon \delta(x) + (1-\epsilon)\phi(x)\). We'll refer to this case as partial loss.
With increasing number of variants, susceptibilities and cross-immunities, model become a lot more complex. But we will consider a simple case where with only one immune escape variant with the modified susceptibility \(\phi'(x)\) and ignore the resident variant. In this case, we still have a single variant model with the same parameters. The variant initially grows with a rate \(\beta \langle S' \rangle - 1\).
To compare how the size of a wave depends on the model assumptions, we numerically solve for a wave with an initial susceptibility \(\phi\) according the shift and partial loss models and compare this to the waves in standard SIR models. For the latter, we assume that (i) every body in the population is partially susceptible (uniform evasion), that that a fraction of the population is fully susceptible and the remainder fully immune (zero/one). In all cases, parameters are adjusted such that the initial growth rates are equal:
# copy the equilibrium susceptibility distribution
phi_base_line = np.copy(sol[-1,1:])
# calculate the shifted probability distribution
delta_x = 1
phi_shift = np.zeros_like(phi_base_line)
phi_shift[:(xgrid>delta_x).sum()] = phi_base_line[xgrid>delta_x]
phi_shift[0] += phi_shift[xgrid<=delta_x].sum()
S_avg_shift = np.sum(np.exp(-xgrid)*phi_shift)*dx
# calculate the probability for the partial_loss case
eps = 0.21
phi_loss = phi_base_line*(1-eps)
phi_loss[0] += phi_base_line.sum()*eps
S_avg_loss = np.sum(np.exp(-xgrid)*phi_loss)*dx
growth_rate = beta*S_avg_shift-1
print(f"average susceptibility -- shift:{S_avg_shift:1.3f}, loss:{S_avg_loss:1.3f}")
print(f"growth rate: {growth_rate:1.3f}")
I = 0.001
T = np.linspace(0,30,31)
Y = np.array([I] + list(phi_shift))
sol_shift = odeint(SusceptI, Y, T, (beta, gamma, D, v, xgrid, dx))
Y = np.array([I] + list(phi_loss))
sol_loss = odeint(SusceptI, Y, T, (beta, gamma, D, v, xgrid, dx))
average susceptibility -- shift:0.365, loss:0.366
growth rate: 0.823
As expected, models with a continuous susceptibility distribution generate wave sizes in between the extreme cases where the entire population is uniformly susceptible, and where a fraction is completely susceptible. The shifted titer-case is closer to the former, the partial loss closer to the latter.