Transmission model

We considered a discrete-time generalised renewal equation model that accounts for heterogeneity in infectiousness between different infected individuals. In the transmission model, each infected individual is assigned an individual infectiousness factor, \(\alpha\), from a gamma distribution with shape parameter \(k\) and mean one (although, in principle, other distributions could be used). If the reproduction number, \(R\), remains constant during that individual’s infection, the realised number of secondary infections generated by the individual is then assumed to follow a Poisson distribution with mean \(\alpha R\). This Poisson-gamma formulation is equivalent to assuming that the realised number of secondary infections from a randomly chosen infected individual follows a negative binomial offspring distribution with mean \(R\) and dispersion parameter \(k\)¹¹\(.\) Smaller values of \(k\) correspond to a greater degree of superspreading (Supplementary Fig. 1). The offspring distribution is a geometric distribution if \(k=1\), while a Poisson distribution is obtained in the limit \(k\to {\infty }\)¹¹\(.\)

More generally, accounting for temporal changes in the reproduction number, \({R}_{t}\), between different times, \(t\), the model-simulated incidence of newly infected individuals, \({I}_{t}\), is generated as follows:

Here, \({w}_{\tau }\) is the probability mass function of the discrete generation time distribution (the probability that the interval between the dates of infection of an infector-infectee pair is \(\tau\) days), and \({J}_{t}={\sum }_{i=1}^{{I}_{t}}{\alpha }_{t,i}\) represents the sum of the individual infectivity factors (\({\alpha }_{t,i}\)) of the individuals infected on day \(t\) (labelled \(i=1,\ldots,{I}_{t}\)). We note that this formulation is valid because the sum of independent gamma-distributed random variables with identical scale parameters also follows a gamma distribution (with the same scale parameter).

Outbreak risk

Analytical approach

The outbreak risk, \({p}_{t}\), is defined to be the probability (under the above transmission model, and assuming that the reproduction number is known at all times) that an outbreak, initiated with a single incident infection at time \(t\) (i.e., taking \({I}_{t}=1\), with \({I}_{s}=0\) for \(s ), does not become extinct. As shown in the Supplementary Methods, the outbreak risk satisfies the following equation:

In our numerical analyses, we assumed that \({R}_{t}\) is a periodic function of time (so that \({p}_{t}\) is also periodic with the same period). In that scenario, a closed system of equations for the values of \({p}_{t}\) (across all values of \(t\)), which can be solved numerically, can be obtained from the above expression (see the Supplementary Methods for details).

If \({R}_{t}=R\) is constant in time, then the implicit equation for the outbreak risk, \({p}_{t}=p\), derived in ref. ¹¹ is recovered:

The outbreak risk is then given by the largest solution less than one of Eq. (4)¹¹.

Simulation-based approach

An alternative method for calculating the outbreak risk at time \(t\) involves simulating the renewal model a large number of times, initiating each simulation with a single incident infection at time \(t\). The outbreak risk is then calculated as the proportion of model simulations in which the daily incidence of new infections ever exceeds a specified threshold. Each simulation is run until either the outbreak ends (i.e., there are no new infections for a period exceeding the maximal generation time) or the incidence threshold is exceeded. When we used the simulation-based approach in Fig. 1E, we conducted 20,000 model simulations for each time of pathogen introduction and used an incidence threshold of 30 infections.

Preliminary analysis of seasonal transmission

We initially considered a simple scenario in which the reproduction number varies periodically with period 365 days:

Over the course of each year, \({R}_{t}\) therefore varies between one and three. We compared the outbreak risk, calculated using either the analytical or simulation-based approach, between different values of the dispersion parameter, \(k\) (Fig. 1E), under the generation time distribution for COVID-19 described below and shown in Fig. 1C.

COVID-19 booster vaccination case study

Antibody dynamics model

We considered a model of individual-level antibody dynamics following a dose of an mRNA COVID-19 vaccine^29,30,31. In this model, the IgG(S) antibody titre, \(A\left(\tau \right)\), of an individual at time since (most recent) vaccination \(\tau\) evolves according to the ordinary differential equation,

where

is the quantity of mRNA inoculated by the vaccine dose remaining. Here, \({\tau }_{d}\) represents a delay before the vaccine elicits an antibody response; for \(0\le \tau , the antibody titre of a previously vaccinated individual is assumed to continue to vary based on the remaining inoculated mRNA and induced antibody titre from their previous vaccination dose (for a previously unvaccinated individual, the antibody titre is assumed to be zero until time \({\tau }_{d}\)). Interpretations of the model parameters \(H\), \(m\), \(K\), \(\mu\), \(D\) and \(\delta\) are given in Supplementary Table 1. We note that we here slightly simplified the model compared to previous studies^29,30,31 by assuming that residual mRNA from previous vaccination doses is negligible (after the delay, \({\tau }_{d}\)).

In a previous study³¹, the antibody dynamics model was fitted to longitudinal data collected following a booster (third overall) dose of the BNT162b2 or mRNA-1273 vaccine. In that study³¹, densely sampled blood data from 12 healthcare workers were considered first, indicating that only the parameters \(H\) and \(m\) exhibit substantial variability between individuals. Then, sparsely sampled data from a larger cohort of 1618 individuals were used to obtain individual estimates of \(H\) and \(m\), with remaining parameters fixed based on model fits to the healthcare worker data. Here, we used the mean and standard deviation of logarithms of individual parameter estimates for the 1618 individuals (Supplementary Table 1) to generate synthetic antibody profiles (see below and the Supplementary Methods), assuming each parameter value varies between individuals according to a lognormal distribution.

Individual susceptibility model

We assumed that individual relative susceptibility (i.e. one minus the level of vaccine protection against infection), \(S\left(\tau \right)\), is given in terms of the antibody titre by

This functional relationship has previously been applied to relate susceptibility to neutralising antibody titres^38,39,40,41, and the application of the same formula here is justified by the high correlation between neutralising antibodies and the IgG(S) antibody titre³¹. We assumed a value of \(\kappa=3.1\) as estimated in ref. ⁴⁰. The parameter \({A}_{1/2}\), representing the antibody titre at which 50% protection is conferred, is likely to depend on both vaccine type and the SARS-CoV-2 variant under consideration^25,38,39,41. Therefore, we took a default value of \({A}_{1/2}={{\mathrm{1,000}}}\) AU/mL—this gave a maximum average vaccine protection against infection of 79% (Fig. 2B), which is comparable to previous model estimates of the effectiveness of a fourth vaccine dose against symptomatic disease due to the omicron SARS-CoV-2 variant²⁵. We also considered alternative values of \({A}_{1/2}\) in Supplementary Fig. 5.

Annual booster vaccination

We suppose that booster vaccines are received by a proportion, \(\theta\), of the population, each year between calendar times \({t}_{s}\) and \({t}_{e}\) (assuming that vaccinations are spread equally between each day in this time window).

To calculate the mean population susceptibility, \({\eta }_{t}\), on each calendar day \(t\), we considered a synthetic cohort of 10,000 individuals. Considering individual variation in antibody dynamics model parameters and assuming an interval of one year between successive vaccine doses, we calculated periodic solutions of the antibody dynamics model for each individual. This allowed us to calculate the expected susceptibility, \(\bar{S}\left(\tau \right)\), of individuals included in booster vaccination campaigns, as a function of time since most recent vaccine dose, \(\tau\). Details of this calculation are given in the Supplementary Methods. Assuming individuals who do not receive booster vaccinations have relative susceptibility one, and for simplicity neglecting variation in susceptibility within each day, we then calculated

In our main analyses, we took \(\theta=0.6\) to represent a vulnerable population with relatively high booster uptake. Alternative values of \(\theta\) are considered in Fig. 4D–F. Initially, we considered booster vaccination rollout between \({t}_{s}=\,\)1 October and \({t}_{e}=\,\)15 December, which is similar to the original planned timetable for 2023 in England⁴² (the 2023 booster vaccination start date was later moved forwards in response to the emergence of the BA.2.86 omicron subvariant⁴³). We then explored how the rollout timing could be optimised to minimise the annual peak outbreak risk.

Seasonal transmissibility

In the absence of booster vaccination (i.e., assuming all individuals are entirely susceptible), we assumed the instantaneous basic reproduction number to be given by^27,44

We assumed that the basic reproduction number is highest on 1 January each year, which is consistent with the respiratory infectious disease burden in the UK typically being highest in December and January⁴².

We took default values of \(\varepsilon=0.2\), as previously assumed for COVID-19 based on estimates for other coronaviruses in ref. ²⁷, and \(\overline{{R}_{0}}=2.5\), giving a maximum annual basic reproduction number (\({R}_{0,t}=3\)) consistent with estimates of the reproduction number of the omicron SARS-CoV-2 variant (BA.1 subvariant) shortly after it emerged in late 2021⁴⁵. Different values of \(\overline{{R}_{0}}\) and \(\varepsilon\) are considered in Supplementary Figs. 3, 4, respectively.

The instantaneous reproduction number accounting for booster vaccination was then calculated as \({R}_{t}={{\eta }_{t}R}_{0,t}\).

Heterogeneity in infectiousness

In our initial analysis, we assumed a dispersion parameter of \(k=0.41\), as estimated for COVID-19 in ref. ¹⁵. Alternative values of \(k\) are considered in Fig. 4A–C.

Generation time distribution

Throughout our analyses, we assumed a lognormal continuous generation time distribution (i.e., distribution of the interval between the exact times of infection of an infector-infectee pair) with mean 3.0 days and standard deviation 1.5 days (this corresponds to lognormal distribution parameters, describing the mean and standard deviation of the natural logarithm of the continuous generation time, of \(\mu=0.98\) and \(\sigma=0.47\), respectively), as estimated for household SARS-CoV-2 omicron variant transmissions in ref. ³². We then discretised this distribution using the method described in refs. ^{46,47, truncating the discretised distribution at 14 days (after which the remaining probability mass is negligible)}. The discretised generation time distribution is shown in Fig. 1C.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.