Uncertainty can be classified as either aleatoric (intrinsic randomness) or epistemic (imperfect knowledge of parameters). The majority of frameworks assessing infectious disease risk consider only epistemic uncertainty. We only ever observe a single epidemic, and therefore cannot empirically determine aleatoric uncertainty. Here, we characterise both epistemic and aleatoric uncertainty using a time-varying general branching process. Our framework explicitly decomposes aleatoric variance into mechanistic components, quantifying the contribution to uncertainty produced by each factor in the epidemic process, and how these contributions vary over time. The aleatoric variance of an outbreak is itself a renewal equation where past variance affects future variance. We find that, superspreading is not necessary for substantial uncertainty, and profound variation in outbreak size can occur even without overdispersion in the offspring distribution (i.e. the distribution of the number of secondary infections an infected person produces). Aleatoric forecasting uncertainty grows dynamically and rapidly, and so forecasting using only epistemic uncertainty is a significant underestimate. Therefore, failure to account for aleatoric uncertainty will ensure that policymakers are misled about the substantially higher true extent of potential risk. We demonstrate our method, and the extent to which potential risk is underestimated, using two historical examples.

§.§ Introduction

Infectious diseases remain a major cause of human mortality. Understanding their dynamics is essential for forecasting cases, hospitalisations, and deaths, and to estimate the impact of interventions. The sequence of infection events defines a particular epidemic trajectory – the outbreak – from which we infer aggregate, population-level quantities. The mathematical link between individual events and aggregate population behaviour is key to inference and forecasting. The two most common analytical frameworks for modelling aggregate data are susceptible-infected-recovered (SIR)  models <cit.> or renewal equation models <cit.>. Under certain specific assumptions, these frameworks are deterministic and equivalent to each other <cit.>. Several general stochastic analytical frameworks exist <cit.>,  and to ensure analytical tractability make strong simplifying assumptions (e.g. Markov or Gaussian) regarding the probabilities of individual events that lead to emergent aggregate behaviour.

We can classify uncertainty as either aleatoric (due to randomness) or epistemic (imprecise knowledge of parameters) <cit.>. The study of uncertainty in infectious disease modelling has a rich history in a range of disciplines, with many different facets <cit.>. These frameworks commonly propose two general mechanisms to drive the infectious process. The first is the infectiousness, which is a probability distribution for how likely an infected individual is to infect someone else. The second is the infectious period, i.e. how long a person remains infectious. The infectious period can also be used to represent isolation, where a person might still be infectious but no longer infects others and therefore is considered to have shortened their infectious period. Consider fitting a renewal equation to observed incidence data <cit.>, where infectiousness is known but the rate of infection events $\rho(\cdot)$ must be fitted. The secondary infections produced by an infected individual will occur randomly over their infectious period $g$, depending on their infectiousness $\nu$. The population mean rate of infection events is given by $\rho(t)$, and we assume that this mean does not differ between individuals (although each individual has a different random draw of their number of secondary infections). In Bayesian settings, inference yields multiple posterior estimates for $\rho$, and therefore multiple incidence values. This is epistemic uncertainty: any given value of $\rho$ corresponds to a single realisation of incidence. However, each posterior estimate of $\rho$ is in fact only the mean of an underlying offspring distribution (i.e. the distribution of the number of secondary infections an infected person produces). If an epidemic governed by identical parameters were to happen again, but with different random draws of infection events, each realisation would be different, thus giving aleatoric uncertainty.

When performing inference, infectious disease models tend to consider epistemic uncertainty only due to the difficulties in performing inference with aleatoric uncertainty (e.g. individual-based models) or analytical tractability. There are many exceptions such as the susceptible-infected-recovered model, which has stochastic variants that are capable of determining aleatoric uncertainty <cit.> and have been used in extensive applications (e.g. <cit.>). However, we will show that this model can underestimate uncertainty under certain conditions. An empirical alternative is to characterise aleatoric uncertainty by the final epidemic size from multiple historical outbreaks <cit.> but these are confounded by temporal, cultural, epidemiological, and biological context, and therefore parameters vary between each outbreak. Here, following previous approaches <cit.>, we analyse aleatoric uncertainty by studying an epidemiologically-motivated stochastic process, serving as a proxy for repeated realisations of an epidemic. Within our framework, we find that using epistemic uncertainty alone is a vast underestimate, and accounting for aleatoric uncertainty shows potential risk to be much higher.  We demonstrate our method using two historical examples: firstly the 2003 severe acute respiratory syndrome (SARS) outbreak in Hong Kong, and secondly the early 2020 UK COVID-19 epidemic.

§.§ Results

§.§.§ An analytical framework for aleatoric uncertainty

A time-varying general branching processes proceeds as follows: first, an individual is infected, and their infectious period is distributed with probability density function $g$ (with corresponding cumulative distribution function $G$). Second, while infectious, individuals randomly infect others (via a counting process with independent increments), driven by their infectiousness $\nu$ and a rate of infection events $\rho$. That is, an individual infected at time $l$, will, at some later time while still infectious $t$, generate secondary infections at a rate $\rho(t)\nu(t-l)$. $\rho(t)$ is a population-level parameter closely related to the time-varying reproduction number $R(t)$ (see Methods and  <cit.> for further details), while $\nu(t-l)$ captures the individual's current infectiousness (note that $t-l$ is the time since infection).  We allow multiple infection events to occur simultaneously, and assume individuals behave independently once infected, thus allowing mathematical tractability <cit.>. Briefly, we model an individual's secondary infections using a stochastic counting process, which gives rise to secondary infections (i.e. offspring) that are either Poisson or Negative Binomial distributed in their number, and Poisson distributed in their timing (see Supplementary Notes 3.3 and 3.4). We study the aggregate of these events (prevalence or incidence) through closed-form probability generating functions and probability mass functions. Our approach models epidemic evolution through intuitive individual-level characteristics while retaining analytical tractability. Importantly, the mean of our process follows a renewal equation <cit.>. Our formulation unifies mechanistic and individual-based modelling within a single analytical framework based on branching processes. Figure <ref> shows a schematic of this process. Formal derivation is in Supplementary Note 3.

Schematic of a time-varying general branching process. (a) shows schematics for the infectious period, an individual's time-varying infectiousness (both functions of time post infection $t^*$), and the population-level mean rate of infection events. The infectious period is given by probability density function $g$. For each individual their (time-varying) infectiousness and rate of infection events are given by $\nu$ and $\rho$ respectively. In an example (b), an individual is infected at time $l$, and infects three people (random variables $K$, purple dashed lines) at times $l+K_1$, $l+K_2$ and $l+K_3$. The times of these infections are given by a random variable with probability density function $\sim \frac{\rho(t)\nu(t-l)}{\int_l^{t}\rho(u)\nu(u-l) du}$. Each new infection then has its own infectious period and secondary infections (thinner coloured lines).

Randomness occurs at individual level, and there is a distribution of possible realisations of the epidemic given identical parameters. Simulating our general branching process would be cumbersome using the standard approach of Poisson thinning <cit.>, and inference from simulation is more challenging still. Using probability generating functions, we analytically derive important quantities from the distribution of the number of infections, including the (central) moments and marginal probabilities given $\rho,g$ and $\nu$ (with or without epistemic uncertainty). We additionally use the probability generating function to prove general, closed-form, analytical results such as the decomposition of variance into mechanistic components, and the conditions under which overdispersion exists (i.e. where variance is greater than the mean). Finally, we derive a general probability mass function (likelihood function) for incidence.

If infection event $k = 0, \ldots, n$ occurred at time $\tau_k$ and produced $y_k$ infections, let $x_{kj}$ denote the end time of the infectious period of the $j^{th}$ infection at event $k$. Note that $\tau_0 = l$ is the time of the first infection event and $y_0 = 1$. Then the likelihood $L_{\text{InfPeriod}}$ of each infected person's infectious period is a product over all infections given by
\begin{equation}
L_{\text{InfPeriod}} = \prod_{k=0}^n \prod_{j=1}^{y_k} g(x_{kj} - \tau_k,\tau_k).
\end{equation}
The likelihood of there being $y_k$ infections at time $\tau_k$ is given by
\begin{equation}
L_{\text{InfTime}} = \prod_{k=1}^n \bigg( \sum_{i = 0}^{k-1}\sum_{j=1}^{y_i} \mathds{1}_{\{x_{ij} < \tau_k\}} p_{y_k}(\tau_k,\tau_i) \bigg),
\end{equation}
where $p_{y_k}(\tau_k,\tau_i)$ is the (infinitesimal) rate at which an individual infected at $\tau_i$ causes $y_k$ infections at time $\tau_k$, provided it is still infectious. Finally, the probability that no other infections occurred between the infection events at times $(\tau_k)_{k=0}^n$ is given by
\begin{equation}
L_{\text{Only}} = \exp\bigg(-\sum_{i=0}^n\sum_{j=1}^{y_i}\int_{\tau_i}^{\min(t,x_{ij})}r(u,\tau_i)du\bigg),
\end{equation}
where $r$ is the infection event rate and $t$ is the current time. Note the term $\text{exp}(-x)$ comes from a Poisson assumption.  Our full likelihood $L_{\text{Full}}$ is then
\begin{equation}
L_{\text{Full}} = L_{\text{InfPeriod}} \times L_{\text{InfTime}} \times L_{\text{Only}}.
\label{eq:likelihood}
\end{equation}

Full derivations of these quantities are provided in Supplementary Note 3. If discrete time is assumed, equation <ref> simplifies to a likelihood commonly used for inference <cit.>. Markov Chain Monte Carlo can be used on equation <ref> to sample aleatoric incidence realisations, but it is often simpler to solve the probability generating function with complex integration. The probability generating function, equations for the variance, and derivations of the probability mass function are found in Supplementary Notes 3,4,5 and 6, and a summary of the main analytical results is found in the Methods.

§.§.§ The dynamics of Uncertainty

We derive the mean and variance of our branching process. The general variance Equation <ref> (see Methods) captures uncertainty in  prevalence over time, where individual-level parameters govern each infection event. This equation comprises three terms: the timing of secondary infections from the infectious period (Equation <ref>a); the offspring distribution (Equation <ref>b); and propagation of uncertainty through the descendants of the initial individual (Equation <ref>c). Importantly, this last term depends on past variance, showing that the infection process itself contributes to aleatoric variance, and does not arise only from uncertainty in individual-level events. In short, unlike common Gaussian stochastic processes, the general variance in disease prevalence is described through a renewal equation. Therefore, future uncertainty depends on past uncertainty, and so the uncertainty around subsequent epidemic waves has memory. Additionally, uncertainty is driven by a complex interplay of time-varying factors, and not simply proportional to the mean. For example, a large first wave of infection can increase the variance of the second wave. As such, the general variance equation <ref> disentangles and quantifies the causes of uncertainty, which remain obscured in brute-force simulation experiments <cit.>.

Consider a toy simulated epidemic with $\rho(t) = 1.4 + \text{sin}(0.15t)$, where the offspring distribution is Poisson in both timing and number of secondary infections, and where infectiousness $\nu$ is given by the probability density function $\nu\sim \text{Gamma}(3,1)$, and, similarly, the infectious period $g \sim \text{Gamma(5,1)}$. Here the parameters of the Gamma distribution are the shape and scale respectively. The resulting variance is counterintuitive. We prove analytically that overdispersion emerges despite a non-overdispersed Poisson offspring distribution. The second wave has a lower mean but a higher variance than the first wave (Figure <ref>), because uncertainty is propagated. If the variance were Poisson, i.e. equal to the mean, the second wave would instead have a smaller variance due to fewer infections. Initially, uncertainty from individuals is largest, but as the epidemic progresses, compounding uncertainty propagated from the past dominates [Figure <ref>, bottom right]. Note that in this example with zero epistemic uncertainty (we know the parameters perfectly), aleatoric uncertainty is large.

Aleatoric uncertainty without overdispersed offspring distribution. Plots show simulated epidemic where $\rho(t) = 1.4 + \text{sin}(0.15t)$, with a Poisson offspring distribution. We use infectiousness $\nu \sim \text{Gamma}(3,1)$, and infectious period $g \sim \text{Gamma(5,1)}$. (a) Overlap between $g$ and the infectiousness $\nu$, where $g$ controls when the infection ends e.g. by isolation. (b) Predicted mean and 95% aleatoric uncertainty intervals for prevalence. Note there is no epistemic uncertainty as the parameters are known exactly (c) Phase plane plot showing the mean plotting against the variance. (d) Proportional contribution to the variance from the individual terms in Equation <ref>. Compounding uncertainty from past events is the dominant contributor to overall uncertainty.

In Equation <ref>, the first two terms account for uncertainty in the infectious periods of all infected individuals. The third term denotes the uncertainty from the offspring distribution. By construction, the timing of infections is an inhomogenous Poisson process, where at each infection time the number of infections is random. The third term (Equation <ref>b) contains the second moment of the offspring distribution, which is the variability around its mean (i.e. $\rho(t)$). The second moment quantifies the extent of possible superspreading. In contrast to other studies <cit.>, we find that individual-level overdispersion in the offspring distribution is less important than explosive epidemics. Under a null Poisson model, with no overdispersion (see Poisson case in Figure <ref>), substantial aleatoric uncertainty arises from a Poisson offspring distribution combined with variance propagation. We rigorously prove via the Cauchy-Schwarz inequality that, under a mild condition on the possible spread of the epidemic, the variance of number of infections at a given time is always greater than the mean, and hence is overdispersed. Overdispersion in the offspring infection distribution is therefore not necessary  for high aleatoric uncertainty, although it still increases variance at both individual-level and population-level.

We derive the conditional variance, with known past events but unknown future events. Conditional variance grows proportionally to the square of the mean, with  additional terms containing the previous variance. Therefore aleatoric uncertainty grows and forecasting exercises based only on epistemic uncertainty greatly underestimates the risk of very large epidemics, and this underestimation becomes more severe as the forecast horizon expands or as the epidemic grows.

§.§.§ Aleatoric uncertainty over the SARS 2003 epidemic

To demonstrate the importance of aleatoric uncertainty, we analyse daily incidence of symptom onset in Hong Kong during the 2003 severe acute respiratory syndrome (SARS) outbreak <cit.>. The epidemic struck Hong Kong in March-May 2003, with a case fatality ratio of 15%. We fit a Bayesian renewal equation assuming a random walk prior distribution for the rate of infection events $\rho$ <cit.>, using Equation <ref> for inference. We ignore $g$ and assume that the distribution of generation times mirrors the distribution of infectiousness, i.e. that the infectiousness $\nu$ equals the generation time <cit.>. Note these parameter choices are illustrative and do not affect our main conclusions. The fitted $\rho(t)$ in Figure <ref> (top left) shows two major peaks, consistent with the major transmission events in the epidemic <cit.>. Figure <ref> (top right) shows the mean epistemic fit, with epistemic (posterior) uncertainty tightly distributed around the data. Figure <ref> (bottom left) shows the aleatoric uncertainty under optimistic and pessimistic scenarios (i.e. the upper and lower bounds of $\rho(t)$ in Figure <ref> (top right)). The pessimistic scenario includes the possibility of extinction, but also an epidemic that could have been more than six times larger than that observed. The optimistic scenario suggests we would observe an epidemic of at worst comparable size to that observed. Finally, Figure <ref> (bottom right) shows epistemic and aleatoric forecasts at day 60 of the epidemic, fixing $\rho(t)$ using the 95% epistemic uncertainty interval to be constant at either $\rho(t\geq 60)=0.38$ or $\rho(t\geq60) = 0.83$ and simulating forwards. While the epistemic forecast does contain the true unobserved outcome of the epidemic, it underestimates true forecast uncertainty, which is 1.3 times larger. The range of the constant $\rho$ for forecast is below 1, and yet we still see substantial aleatoric uncertainty. If $\rho$ were above 1 for a sustained period, aleatoric uncertainty would play a smaller role <cit.>, but this is rare with real epidemics, where susceptible depletion, behavioural changes or interventions keep $\rho$ around 1. Our results therefore highlight that epistemic uncertainty drastically underestimates potential epidemic risk.

The 2003 SARS epidemic in Hong Kong <cit.>. (a) $\rho(t)$ with 95% epistemic uncertainty. (b) Fitted incidence mean, 95% epistemic uncertainty with observational noise from using Equation <ref>. Data is daily incidence of symptom onset. (c) Aleatoric uncertainty from the start of the epidemic under an optimistic and pessimistic $\rho(t)$. (d) Epistemic (blue) and epistemic and aleatoric uncertainty (red) while keeping $\rho$ constant at the forecast data (dotted line). Forecasting is from day 60.

§.§.§ Aleatoric risk assessment in the early 2020 COVID-19 pandemic in the UK

To demonstrate the practical application of our model, we retrospectively examine the early stage of the COVID-19 pandemic in the UK, using only information available at the time. While the date of the first locally transmitted case in the UK remains unknown (likely mid-January 2020 <cit.>), COVID-19 community transmission was confirmed in the UK by late January 2020, and we therefore start our simulated epidemic on January 31st 2020. We consider uncertainty in the predicted number of deaths on March 16th 2020 <cit.>, during which time decisions regarding non-pharmaceutical interventions were made. Testing was extremely limited during this period, and COVID-19 death data were unreliable. For this illustration, we assume that we did not know the true number of COVID-19 deaths, as was the case for many countries in early 2020. Policymakers then needed estimates of the potential death toll, given limited knowledge of COVID-19 epidemiology and unreliable national surveillance.

We simulated an epidemic from a time-varying general branching process with a  Negative Binomial offspring distribution, using parameters that were largely known by March 16th 2020 (Table <ref>). The infection fatality ratio, infection-to-onset distribution and onset-to-death distribution were convoluted with incidence <cit.> to estimate numbers of deaths. Estimated COVID-19 deaths and uncertainty estimates between January 31st and March 16th 2020 are shown in Figure <ref> (Top). While the epistemic uncertainty contains the true number of deaths, it is still an underestimate, and including aleatoric uncertainty, we find that the epidemic could have had more than four times as many deaths. Consider a hypothetical intervention on March 17th 2020 (Figure <ref> (bottom)) that completely stops transmission. Deaths would still occur from those already infected but no new infections would arise. In this hypothetical case, the aleatoric uncertainty would still be 2.5 times the actual deaths that occurred (when in fact transmission was never zero or close to it). This hypothetical scenario highlights the scale of aleatoric uncertainty, and demonstrates that our method can be useful in assessing risk in the absence of data by giving a reasonable worst case. Further, we observe that using only epistemic uncertainty provides a reasonably good fit in a relatively short time-horizon (Figure <ref>, Top), but soon afterwards greatly underestimates uncertainty (Figure <ref>, Bottom). The fits using aleatoric uncertainty provide a more reasonable assessment of uncertainty. While we concentrate on the upper bound, the lower bound on the worst-case scenario still exceeds zero, and therefore the epidemic going extinct by March 16th in the worst-case with no external seeding would have been very unlikely. Aleatoric uncertainty highlights a more informative reasonable worst-case estimate than epistemic uncertainty alone, and could be a useful metric for a policymaker in real time, with low-quality data, without requiring simulations from costly, individual-based models.

Epidemiological Parameter     Value or Distribution               Citation
Infection Fatality Ratio                0.9%                                        <cit.>
Basic Reproduction Number              $2-4$                                        <cit.>
Serial Interval Distribution           $\sim \text{Gamma}(7.82,\frac{1}{0.62})$     <cit.>
Onset-to-Death Distribution            $\sim \text{Gamma}(1.45, 10.43)$             <cit.>
Infection-to-Onset Distribution        $\sim \text{Gamma}(35.16, 6.9)$              <cit.>
Overdispersion Coefficient             0.53                                         <cit.>

Epidemiological parameters available on March 16th 2020 used in branching process simulation

Early 2020 COVID-19 pandemic in the UK. (a) shows a simulated epidemic using parameters available on March 16th 2020 (Table <ref>), for a plausible range of $\rho=R_0$ between 2 and 4. Blue bars indicate actual COVID-19 deaths, which we assume no knowledge of. The purple line is March 17th 2020, we set transmission to zero i.e. $\rho=0$, to simulate an intervention that stops transmission completely. The grey envelope is the epistemic uncertainty and the red envelope the aleatoric uncertainty. (b) is the same as the top plot, except time is extended past March 17th with transmission being zero.  Note aleatoric uncertainty is presented but is very close to zero.

§.§ Discussion

Stochastic models more realistically model natural phenomena than deterministic equations <cit.>, and particularly so with infection processes <cit.>. Accordingly, individual-based models have found much success <cit.> in capturing the complex dynamics that emerge from infectious disease outbreaks, and have been highly influential in policy <cit.>. However, despite a plethora of alternatives, many analytical frameworks still tend to be deterministic <cit.>, and only consider statistical, epistemic parameter uncertainty. Frameworks that expand deterministic, mechanistic equations to include stochasticity use a Gaussian noise process <cit.>, or restrict the process to be Markovian. Markovian branching processes require the infection period or generation time to be exponentially distributed - a fundamentally unrealistic choice for most infectious diseases. Further, a Gaussian noise process is unlikely to be realistic <cit.>.

Our results show that individual-level uncertainty is overshadowed by uncertainty in the infection process itself. Profound overdispersion in infectious disease epidemics is not simply a result of overdispersion in the offspring distribution, but is fundamental and inherent to the branching process. We rigorously prove that even with a Poisson offspring distribution (not characterized by overdispersion), overdispersion in resulting prevalence or incidence is still virtually always guaranteed. We show that forecast uncertainty increases rapidly, and therefore common forecasting methods almost certainly underestimate true uncertainty. Similar to other existing frameworks, our approach provides a different methodological tool to evaluate uncertainty in the presence of little to no data, assess uncertainty in forecasting, and  retrospectively assess an epidemic. Other approaches, such as agent based models, could also be readily used. However, the framework we present permits the unpicking of dynamics analytically and from first principles without a black box simulator. Equally, this is also a limitation, since new and flexible mechanisms cannot be easily integrated or considered.

We have considered only a small number of mechanisms that generate uncertainty. Cultural, behavioural and socioeconomic factors could introduce even greater randomness. Therefore our framework may underestimate true uncertainty in infectious disease epidemics. The converse is also likely, contact network patterns and spatial heterogeneity also limit the routes of transmission, such that the variability in anything but a fully connected network will be lower. Furthermore, our assumption of homogeneous mixing and spatial independence overestimates uncertainty. A sensible next step for future research to to study the dynamics of these branching processes over complex networks. Finally at the core of all branching frameworks in an assumption of independence, which is unlikely to be completely valid (people mimic other people in their behaviour) but is necessary for analytical tractability. Studying the effect of this assumption compared to agent based models would also be a useful area of future research.

We provide one approach to determining aleatoric uncertainty. Other approaches based on stochastic differential equations, Markov processes, reaction kinetics, or Hawkes processes all have their respective advantages and disadvantages. The differences in model specific aleatoric uncertainty and how close the models come to capturing the true, unknown, aleatoric uncertainty is a fundamental question moving forwards. In this paper we have provided yet another approach to characterise aleatoric uncertainty, where this approach is most useful and how it can be reconciled with existing approaches will be an interesting area of study.

§.§ Methods

Detailed derivations of the methods can be found in the Supplementary Notes, with a high level description of the content found in Supplementary Note 1.

A time-varying general branching process proceeds as follows: first, a single individual is infected at some time $l$, and their infectious period $L$ is distributed with probability density function $g$ (and cumulative distribution function $G$). Second, during their infectious period, they randomly infect other individuals, affected by their infectiousness $\nu(t-l)$, and their mean number of secondary infections, which is assumed to be equal to the population-level rate of infection events $\rho(t)$. $\rho(t)$ is closely related to the time-varying reproduction number $R(t)$ (see <cit.> for details). The infectious period $g$ accounts for variation in individual behaviour. If people take preventative action to reduce onward infections, their reduced infection period can stop transmission despite remaining infectious. Where infectious individuals do not change their behaviour, $g$ can be ignored and individual-level transmission is controlled by infectiousness $\nu$ only. Each newly infected individual then proceeds independently by the same mechanism as above.  Specifics can be found in Supplementary Notes 2.1-2.5.

Formally, if an individual is infected at time $s$, their number of secondary infections is given by a stochastic counting process $\{N(t,s)\}_{t\geq s}$, which is independent of other individuals and has independent increments. We assumehere that the epidemic occurs in continuous time, and hence that $N(t,s)$ is continuous in probability, although we consider discrete-time epidemics in Supplementary Note 7. To aid calculation, we suppose $N(t,s)$ can be defined from a Lévy Process $\mathcal{N}(t)$ - that is, a process with both independent and identically distributed increments - via $N(t,s) = \mathcal{N}\bigg(\int_s^tr(k,s)dk\bigg)$ for some non-negative rate function $r$. It is assumed that each counting process $\{N(t,s)\}_{t\geq s}$ is defined from an independent copy of $M(t)$. This formulation has two advantages: first, the dependence of $N(t,s)$ on $s$ is restricted to the rate function $r$; and second, if $J_{\mathcal{N}}(t)$ counts the number of infection events in $\mathcal{N}(t)$ (where here infection events refer to an increase, of any size, in $N(t,s)$), then $J_{\mathcal{N}}(t)$ is a Poisson process with some rate $\kappa$ <cit.>. We can then define $J(t,s)$ to be the counting process of infection events in $N(t,s)$, and $Y(v)$ to be size of the infection event (i.e. the number of secondary infections that occur) t time $v$. We assume that $Y$ is independent of $s$, although such a dependence would curtail superspreading to depend on infectiousness, and could be incorporated into the framework. Therefore $J(t,s)$ is an inhomogeneous Poisson Process (and so $N(t,s)$ has been characterised as an inhomogeneous compound Poisson Process). We consider the cases where $N(t,s)$ is itself an inhomogeneous Poisson process, and where $N(t,s)$ is a Negative Binomial process. This allows us to examine effects of overdispersion in the number of secondary infections, although our framework allows for more complicated distributions.

Here, $r(t,l) = \rho(t)\nu(t-l)$ where $\rho(t)$ models the population-level rate of infection events, and $\nu(t-l)$ models the infectiousness of an individual infected at time $l$. If $\nu(t-l)$ is sufficiently well characterised by the generation time (i.e. where the timing of secondary infections mirrors tracks their infectiousness) , and the infectious period can be ignored, then the integral $\int_l^t r(s,l) ds$ has the same scale as the commonly used reproduction number $R(t)$ <cit.>. The branching process yields a series of birth and death times for each individual (i.e. the time of infection and the end of the infectious period respectively), from which prevalence (the number of infections at any given time) or cumulative incidence (the total number of infections up to any time) can be defined.

§.§.§ Probability generating function

We derive the probability generating function for a time-varying age-dependent branching process, allowing derivation of the mean and higher-order moments (full derivations can be found in Supplementary Notes 3.1-3.7). We consider two special cases for the number of new infections $Y(v)$ at each infection event: a Poisson distribution and a logarithmic (log series) distribution. In both cases, we assume that the distribution of $Y(v)$ is equal for all values of $v$. In the Poisson case, the number of new infections at each infection time is, by definition, one. Therefore the number of infections an individual creates is Poisson distributed, and closely clustered around the mean rate of infection events. The logarithmic case, which causes $N(t,l)$ to be a Negative Binomial process, more realistically allows multiple infections to occur at each infection time, and so the number of infections an individual causes is overdispersed. The pgf (probability generating function), $F(t,l;s) = E(s^{Z(t,l)})$, can be derived by conditioning on the lifetime, $L$, of the first individual. That is,

\begin{equation}
E\bigg(s^{Z(t,l)}\bigg)  = \bigg(1-G(t-l,l)\bigg) E\bigg(s^{Z(t,l)}\bigg|L \geq t-l\bigg) + \int_0^{t-l}E\bigg(s^{Z(t,l)}\bigg|L = u\bigg) g(u,l)du
\end{equation}

Note that if the individuals directly infected by the initial individual are infected at times $l+t_1,...,l+t_n$, then

\begin{equation}
Z(t,l) = 1 + \sum_{i=1}^nZ(t,l+t_i)
\end{equation}

This observation allows us to write the generating function $F(t,l)$ as a function of $F(t,u)$ for $u \in (t,l)$. As $F(t,t) = s$, this allows us to iteratively find the value of $F(t,l)$. Explicitly, we have

\begin{align}
\label{eq:pgf}
\underbrace{F(t,l;s)}_{\text{pgf}} &= \underbrace{\left(1-G(t-l,l)\right)}_{P(L > t-l)}\underbrace{q_1\left(\int_0^{t-l}\overbrace{f(\underbrace{ F(t,l+k;s)}_{\text{pgf of process started at $l+k$}})}^{\text{pgf of $Y$}}\overbrace{\rho(l+k)\nu(k)dk}^\text{infection rate at time $l+k$}  \right)}_{\text{pgf of $J(t,l)$}} \nonumber\\
\nonumber \\
&+\int_0^{t-l}\underbrace{q_2\left( \int_0^u \overbrace{f(\underbrace{ F(t,l+k;s)}_{\text{pgf of process started at $l+k$}})}^{\text{pgf of $Y$}} \overbrace{\rho(l+k)\nu(k)}^{\text{infection rate at time $l+k$}}dk \right)}_{\text{pgf of $J(l+u,l)$}}\underbrace{g(u,l) du}_{P(L = l+u)},
\end{align}

where $q_1(z;s) = se^z$, and where $q_2(z)=e^z$ in the case where $Z(t,l)$ refers to prevalence, whereas $q_2(z;s)=se^z$ in the case where $Z(t,l)$ refers to cumulative incidence. Note also that $f(z) = z-1$ in the Poisson case and $f(z) = -\phi\left[\log\left(1-(1  -\frac{\phi}{1+\phi}z)\right) - \log\left(\frac{\phi}{1+\phi}\right)\right] $ in the log-series case and that the constant $\kappa$ is absorbed into $\rho$.

The key intuition in understanding Equation <ref> is that for an integer random variable $X$ and iid (independent and identically distributed) random variables $Y_i$, $E(s^{\sum_{i=1}^X Y_i}) = G_X(G_Y(s))$, where $G_X$ and $G_Y$ are the generating functions of $X$ and $Y_i$ respectively. Thus, we expect the pgfs of the various parts of our model to combine via composition, as occurs in the equation above.

Mean incidence can recovered from both prevalence (via back calculation <cit.>) and cumulative incidence. In Equation <ref> for the Negative Binomial case, $\phi$ is the degree of overdispersion. Equation <ref> is solvable using via quadrature and the fast Fourier transform via a result from complex analysis <cit.> and scales easily to populations with millions of infected individuals, and the probability mass function can be computed to machine precision (a full derivation is available in Supplementary Note 3.7).

§.§.§ Variance decomposition

For simplicity, we only summarise the decomposition for prevalence, but an analogous and highly similar derivation for cumulative incidence can be found in Supplementary Note 3.5.  We can derive an analytical equation for the mean and variance of the entire branching process (full derivations can be found in Supplementary Notes 4.1-4.7 and the mathematical properties of the variance equations can be found in Supplementary Notes 6.1-6.3) .
The mean prevalence $M(t,l)$ is given by
\begin{equation}
M(t,l) = (1-G(t-l,l)) + \int_0^{t-l} M(t,l+u)\rho(l+u)\nu(u)\mathbb{E}(Y) (1-G(u,l))du.
\label{eq:mean_renewal}
\end{equation}

Note, $\rho$ can be scaled to absorb the $E(Y)$ and $\kappa$ constants. Equation <ref> is consistent with that previously derived in <cit.>. The second moment, $W(t,l) := \mathbb{E}(Z(t,l)(Z(t,l) - 1))$ allows us to determine the variance, $V(t,l)$ as $V(t,l) = W(t,l) + M(t,l) - M(t,l)^2$.
The variance can be decomposed into three mechanistic components.

\begin{align}
V(t,l) &= \underbrace{\int_0^{t-l} \left[\int_0^u M(t,l+k)\rho(l+k)\nu(k) dk\right]^2 g(u,l)du - M(t,l)^2 }_{\text{(9a): uncertainty from the infectious period
&+\underbrace{(1-G(t-l,l))\bigg[1 + 2\int_0^{t-l} M(t,l+u)\rho(l+u)\nu(u)du +  \bigg(\int_0^{t-l}M(t,l+u)\rho(l+u)\nu(u)du \bigg)^2\bigg]}_{\text{(9a continued): uncertainty from the infectious period}}\nonumber\\
&+ \underbrace{\int_0^{t-l}M(t,l+u)^2\mathbb{E}(Y^2)\rho(l+u)\nu(u)(1-G(u,l))du}_{\text{(9b): uncertainty from the offspring distribution}}\nonumber\\
&+\underbrace{\int_0^{t-l}V(t,l+u)\rho(l+u)\nu(u)(1-G(u,l))du}_{\text{(9c): uncertainty propagated from the past}}.
\label{eq:variance_renewal}
\end{align}

The general variance equation <ref> captures the evolution of uncertainty in population-level disease prevalence over time, where fixed individual-level disease transmission parameters govern each infection event. Unlike the simple Galton-Watson process, we find that previously unknown factors also determine aleatoric variation in disease prevalence. Specifically, the general variance equation <ref> comprises three terms, one for the infectious period (Equation <ref>a), one for the number and timing of secondary infections (Equation <ref>b), and a term that propagates uncertainty through descendants of the initial individual (Equation <ref>c). Importantly, the last term (Equation <ref>c) depends on past variance, showing that the infection process itself contributes to aleatoric variance, and this is distinct from the uncertainty in individual infection events. In short, and unlike Gaussian stochastic processes, the general variance in disease prevalence is described through a renewal equation. Intuitively then, uncertainty in an epidemic's future trajectory is contingent on past infections, and that the uncertainty around consecutive epidemic waves are connected.  As such, the general variance equation <ref> allows us to disentangle important aspects of infection dynamics that remain obscured in brute-force simulations <cit.>.

§.§.§ Overdispersion

We define an epidemic to be expanded if at time $t$ there is a non-zero probability that the prevalence, not counting the initial individual or its secondary infections, is non-zero.

Note that this is a very mild condition on an epidemic - in a realistic setting, the only way for an epidemic to not be expanded is if it is definitely extinct by time $t$, or if $t$ is small enough that tertiary infections have not yet occurred.

Large aleatoric variance intrinsic to our branching process implies that the prevalence of new infections (that is, prevalence excluding the deterministic initial case) is always strictly overdispersed at time $t$, providing the epidemic is expanded at time $t$. A full proof is given in Supplementary Note 4.4, but we provide here a simpler justification in the special case that $G(t-l,l) = 1$.

In this case,  prevalence of new infections is equal to  standard prevalence, and the equations for $M(t,l)$ and $V(t,l)$ simplify significantly. Switching the order of integration in the equation for $M(t,l)$ gives

\begin{equation}
M(t,l) =  \int_0^{t-l} M(t,l+u)\rho(l+u)\nu(u) \mathbb{E}(Y)(1-G(u,l))du =  \int_0^{t-l}\bigg[\int_0^u M(t,l+k)\rho(l+k)\mathbb{E}(Y)\nu(k)dk\bigg]g(u,l)du
\end{equation}

and hence, the Cauchy-Schwarz Inequality shows that

\begin{equation}
\label{eq:CS}
M(t,l)^2 \leq \bigg(\int_0^{t-l}\bigg[\int_0^u M(t,l+k)\mathbb{E}(Y)\rho(l+k)\nu(k)dk\bigg]^2g(u,l)du\bigg)
\end{equation}

as $\int_0^{t-l}g(u,l)du = 1$. Thus, the first term, (<ref>a), in the variance equation is non-negative.

The remaining terms can be dealt with as follows. (<ref>a) is equal to zero, and the sum of (<ref>c) is (using $Y(l+u,l)^2 \geq Y(l+u,l)$) bounded below by $\int_0^{t-l} \mathbb{E}(Z(t,l+u)^2)\mathbb{E}(Y)\rho(l+u)\nu(u)(1-G(u,l)) du$. Finally, noting that $Z(t,l+u)^2 \geq Z(t,l+u)$, this is bounded below by $\int_0^{t-l} M(t,l+u)\mathbb{E}(Y)\rho(l+u)\nu(u)(1-G(u,l))du = M(t,l)$. Hence, $V(t,l) \geq M(t,l)$ holds.

To show strict overdispersion, note that for $V(t,l) = M(t,l)$ to hold, it is necessary that

\begin{equation}
\int_0^{t-l} \mathbb{E}(Z(t,l+u)^2)\mathbb{E}(Y)\rho(l+u)\nu(u) (1-G(u,l))du = \int_0^{t-l} M(t,l+u)\mathbb{E}(Y)\rho(l+u)\nu(u) (1-G(u,l))du
\end{equation}

and hence, for each $u$ (as $\mathbb{E}(Y) > 0$)

\begin{equation}
\label{eq:condition}
\mathbb{E}\bigg[Z(t,l+u)(Z(t,l+u)-1)\bigg] = 0 \quad \text{or} \quad \rho(l+u)\nu(u) (1-G(u,l))=0
\end{equation}

If new infections can be caused, then more than one new infection can be caused. Thus, if an individual infected at $l+u$ has $\mathbb{E}\bigg[Z(t,l+u)(Z(t,l+u)-1)\bigg] = 0$, this individual cannot cause new infections whose infection trees have non-zero prevalence at time $l+u$. Hence, the condition <ref> is equivalent to the epidemic being non-expanded at time $t$, as at each time $l+u$, either no infections are possible from the initial individual, or any individuals that are infected at time $l+u$ contribute zero prevalence at time $t$ from the new infections they cause.

Hence, $Z(t,l)$ is strictly overdispersed for expanded epidemics. This means that Gaussian approximations are unlikely to be useful.

§.§.§ Variance midway through an epidemic

It is important to calculate uncertainty starting midway through an epidemic, conditional on previous events. This derivation is significantly more algebraically involved than the other work in this paper. For simplicity, we assume that $N(t,l)$ is an inhomogeneous Poisson Process, and that $L = \infty$ for each individual.

Suppose that prevalence (here equivalent to cumulative incidence) $Z(t,l) = n+1$. We create a strictly increasing sequence $l = B_0 < B_1<\cdots<B_n$ of $n+1$ infection times, which has probability density function

\begin{equation}
\label{eq:probability_mass_function}
\underbrace{f_{\boldsymbol{B}}(\boldsymbol{b})}_{\text{Joint pdf}} =\underbrace{\frac{1}{P(Z(t,l) = n+1)}}_{\text{normalising constant}} \underbrace{\prod_{i=1}^n\bigg(\rho(b_i) \sum_{j=0}^{i-1}\nu(b_i-b_j)\bigg)}_{\text{infection rates at each $b_i$}}\underbrace{\exp\bigg[-\sum_{i=0}^n \int_0^{t-b_i} \rho(s+l)\nu(s)ds\bigg]}_{\text{probability of no other infections}},
\end{equation}

where pdf is short for probability mass function. Then, the variance at time $t+s$ is given by

\begin{align}
&\underbrace{\text{var}(Z(t+s,l))}_{\text{variance}} =\underbrace{\int_{b=0}^t\sum_{i=0}^nV^*(t+s,b)f_{B_i}(b)db}_{\text{variance from subsequent cases}}...\nonumber\\
&...+ \underbrace{\int_{b=0}^t\int_{c=0}^t\sum_{i=0}^n\sum_{j=0}^nM^*(t+s,b)M^*(t+s,c)(f_{B_i,B_j}(b,c) - f_{B_i}(b)f_{B_j}(c))dbdc}_{\text{variance from unknown infection times}},
\end{align}

where $M^*(t+s,b)$ and $V^*(t+s,b)$ are the mean and variance of the size of the infection tree (i.e. prevalence or cumulative incidence) at time $t+s$, caused by an individual infected at time $b$, ignoring all individuals they infected before time $t$. These quantities are calculated from $M$ and $V$. Note also that $f_{B_i}$ and $f_{B_i,B_j}$ are the one-and-two-dimensional marginal distributions from $f_{\boldsymbol{B}}$.

§.§.§ Bayesian inference and for SARS epidemic in Hong Kong

The data for the SARS epidemic in Hong Kong consist of 114 daily measurements of incidence (positive integers), and an estimate of the generation time <cit.> obtained via the R package EpiEstim <cit.>. We ignore the infectious period $g$ and set the infectiousness $\nu$ to the generation interval. The inferential task is then to estimate a time varying function $\rho$ from these data using Equation <ref>. As we note in Equation <ref> and in Supplementary Note 5 and 7.1-7,4, discretisation simplifies this task considerably. Our prior distributions are as follows

\begin{align*}
\phi & \sim \text{Normal}^+(0,1)\\
\sigma &\sim \text{Exponential(100)}\\
\epsilon & \sim \text{Normal}(0,\sigma)\\
\rho(t) &= \rho(t-1) + \epsilon_t
\end{align*}

where $\rho$ is modelled as a discrete random walk process. The renewal likelihood in Equation <ref> is vectorised using the approach described in <cit.>. Fitting was performed in the probabilistic programming language Numpyro, using Hamiltonian Monte Carlo<cit.> with 1000 warmup steps and 6000 sampling steps across two chains. The target acceptance probability was set at 0.99 with a tree depth of 15. Convergence was evaluated using the RHat statistic<cit.>.

Forecasts were implemented through sampling using MCMC from Equation <ref>. In order to use Hamiltonian Markov Chain Monte Carlo, we relax the discrete constraint on incidence and allow it to be continuous with a diffuse prior. We ran a basic sensitivity analysis using a Random Walk Metropolis with a discrete prior to ensure this relaxation was suitable. In a forecast setting, incidence up to a time point ($T=60$) is known exactly and given as $y^{t\leq T}$. and we have access to an estimate for $\rho(t>T)$ in the future. In our case we fix $\rho(t>T)=\rho(T)$.

Our code is available at available at <https://github.com/MLGlobalHealth/uncertainity_infectious_diseases.git>.

§.§.§ Numerically calculating the probability mass function via the probability generating function

Following <cit.> and <cit.> (originally from <cit.>), the probability mass function $p$ can be recovered through a pgf $F$'s derivatives at $s=0$. i.e. $\mathbb{P}(n) = \frac{1}{n!}\left(\frac{d}{ds}\right)^n F(s;t,\tau) |_{s=0}$
This is generally computationally intractable.
A well-known result from complex analysis <cit.> holds that $ f^{(n)}(a) = \frac{n!}{2\pi i} \oint \frac{f(z)}{\left(z-a\right)^{n+1}}\, dz$ and therefore
$ \mathbb{P}(n) = \frac{1}{2\pi i}\oint\frac{F(z;t,\tau)}{z^{n+1}} dz$
This integral can be very well approximated via trapezoidal sums as
$\mathbb{P}(n) = \frac{1}{Mr^n}\sum_{m=0}^{M-1}F(re^{2\pi i m/M};t,\tau) e^{ - 2\pi i n m/M}
where $r=1$<cit.>. The probability mass function for any time and $n$ can be determined numerically. One needs $M\geq n$, which requires solving $n$ renewal equations for the generating function and performing a fast Fourier transform. This is computationally fast, but may become slightly burdensome for epidemics with very large numbers of infected individuals (millions). A derivation of this approximation is provided in the Supplementary Note 3.7.

§.§.§ Competing interests

All authors declare no competing interests

§.§.§ Data availability

Data from Figure 3 is available via the R-Package EpiEstim<cit.>, and data from Figure 4 is available at <https://imperialcollegelondon.github.io/covid19local> and via official UK Government reporting (<https://www.ons.gov.uk/>).

§.§.§ Code availability

All model code to reproduce Figures 2, 3 and 4 is available at <https://github.com/MLGlobalHealth/uncertainity_infectious_diseases.git>.

§.§.§ Acknowledgements

S.B.  C.A.D  and D.J.L. acknowledge support from the MRC Centre for Global Infectious Disease Analysis (MR/R015600/1), jointly funded by the UK Medical Research Council (MRC) and the UK Foreign, Commonwealth & Development Office (FCDO), under the MRC/FCDO Concordat agreement, and also part of the EDCTP2 programme supported by the European Union. S.B. acknowledges support from the Novo Nordisk Foundation via The Novo Nordisk Young Investigator Award (NNF20OC0059309), which also supports S.M.Ṡ.B. acknowledges support from the Danish National Research Foundation via a chair position. S.B. and C.M. acknowledges support from The Eric and Wendy Schmidt Fund For Strategic Innovation via the Schmidt Polymath Award (G-22-63345). S.B. acknowledges support from the  National Institute for Health Research (NIHR) via the Health Protection Research Unit in Modelling and Health Economics. D.J.L. acknowledges funding from Vaccine Efficacy Evaluation for Priority Emerging Diseases (VEEPED) grant, (ref. NIHR:PR-OD-1017-20002) from the National Institute for Health Research. M.J.P. acknowledges funding from a EPSRC DTP Studentship. C.W. acknowledges support from the Wellcome Trust.

§.§.§ Author Contributions

S.B. and M.J.P conceived and designed the study. S.B. performed analysis with assistance from M.J.P. M.J.P, D.J.L and S.B drafted the original manuscript. M.J.P drafted the supplementary information with assistance from J.P. M.J.P, D.J.L, J.P, C.W, C.M, O.R, S.M, M.P, C.A.D, and S.B revised the manuscript and contributed to its scientific interpretation. S.B, M.P and C.A.D supervised the work.

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

titleEpiNow2: Estimate Real-Time Case Counts and Time-Varying Epidemiological Parameters

abstractSome mathematical methods for formulation and numerical simulation of stochastic epidemic models are presented. Specifically, models are formulated for continuous-time Markov chains and stochastic differential equations. Some well-known examples are used for illustration such as an SIR epidemic model and a host-vector malaria model. Analytical methods for approximating the probability of a disease outbreak are also discussed.
journaltitleInfect Dis Model
titleA primer on stochastic epidemic models: Formulation, numerical simulation, and analysis
60H10; 60J28; 92D30; Branching process; Continuous-time Markov chain; Minor outbreak; Stochastic differential equation

American Association for the Advancement of Science (AAAS)

titleMore is different

Cambridge University Press

abstractLévy processes form a wide and rich class of random process, and have many applications ranging from physics to finance. Stochastic calculus is the mathematics of systems interacting with random noise. Here, the author ties these two subjects together, beginning with an introduction to the general theory of Lévy processes, then leading on to develop the stochastic calculus for Lévy processes in a direct and accessible way. This fully revised edition now features a number of new topics. These include: regular variation and subexponential distributions; necessary and sufficient conditions for Lévy processes to have finite moments; characterization of Lévy processes with finite variation; Kunita's estimates for moments of Lévy type stochastic integrals; new proofs of Ito representation and martingale representation theorems for general Lévy processes; multiple Wiener-Lévy integrals and chaos decomposition; an introduction to Malliavin calculus; an introduction to stability theory for Lévy-driven SDEs.
titleLévy Processes and Stochastic Calculus

Institute of Mathematical Statistics
Bernoulli Society

abstractMany models of epidemic spread have a common qualitative structure. The numbers of infected individuals during the initial stages of an epidemic can be well approximated by a branching process, after which the proportion of individuals that are susceptible follows a more or less deterministic course. In this paper, we show that both of these features are consequences of assuming a locally branching structure in the models, and that the deterministic course can itself be determined from the distribution of the limiting random variable associated with the backward, susceptibility branching process. Examples considered includea stochastic version of the Kermack & McKendrick model, the Reed-Frost model, and the Volz configuration model.
titleApproximating the epidemic curve
60J85; 60K35; 92H30; branching processes; configuration model; deterministic approximation; epidemics; Reed–Frost;

Cambridge University Press

abstractThis paper is concerned with negative binomial processes which are essentially mixed Poisson processes whose intensity parameter is given by the sum of squares of a finite number of independently and identically distributed Gaussian processes. A study is made of the distribution of the number of points of a k-dimensional negative binomial process in a compact subset of Rk, and in particular in the case where the underlying Gaussian processes are independent Ornstein-Uhlenbeck processes when more detailed results may be obtained.
journaltitleJ. Appl. Probab.
titleNegative binomial processes

Springer Science
Business Media LLC

journaltitleFound. Comut. Math.
titleAccuracy and stability of computing high-order derivatives of analytic functions by Cauchy integrals

abstractBackground: The infection fatality ratio (IFR) is a key statistic for estimating the burden of coronavirus disease 2019 (COVID-19) and has been continuously debated throughout the COVID-19 pandemic. The age-specific IFR can be quantified using antibody surveys to estimate total infections, but requires consideration of delay-distributions from time from infection to seroconversion, time to death, and time to seroreversion (i.e. antibody waning) alongside serologic test sensitivity and specificity. Previous IFR estimates have not fully propagated uncertainty or accounted for these potential biases, particularly seroreversion. Methods: We built a Bayesian statistical model that incorporates these factors and applied this model to simulated data and 10 serologic studies from different countries. Results: We demonstrate that seroreversion becomes a crucial factor as time accrues but is less important during first-wave, short-term dynamics. We additionally show that disaggregating surveys by regions with higher versus lower disease burden can inform serologic test specificity estimates. The overall IFR in each setting was estimated at 0.49-2.53%. Conclusion: We developed a robust statistical framework to account for full uncertainties in the parameters determining IFR. We provide code for others to apply these methods to further datasets and future epidemics.
journaltitleCommun. Med.
titleEstimating the COVID-19 infection fatality ratio accounting for seroreversion using statistical modelling
Computational biology and bioinformatics; Respiratory tract diseases

journaltitleProceedings of the National Academy of Sciences
titleThe turning point and end of an expanding epidemic cannot be precisely forecast

journaltitleSIAM Journal on Applied Mathematics
titleEquivalence of the Erlang-distributed SEIR epidemic model and the renewal equation
Differential equations,Epidemic models,Erlang distribution,Generation-interval distribution,Renewal equation,SEIR

abstractWe use two modelling approaches to forecast synthetic Ebola epidemics in the context of the RAPIDD Ebola Forecasting Challenge. The first approach is a standard stochastic compartmental model that aims to forecast incidence, hospitalization and deaths among both the general population and health care workers. The second is a model based on the renewal equation with latent variables that forecasts incidence in the whole population only. We describe fitting and forecasting procedures for each model and discuss their advantages and drawbacks. We did not find that one model was consistently better in forecasting than the other.
titleTwo approaches to forecast Ebola synthetic epidemics
Compartmental model; Epidemics; Forecasting; Renewal equation

Nature Publishing Group

abstractThe COVID-19 pandemic has been a sobering reminder of the extensive damage brought about by epidemics, phenomena that play a vivid role in our collective memory, and that have long been identified as significant sources of risk for humanity. The use of increasingly sophisticated mathematical and computational models for the spreading and the implications of epidemics should, in principle, provide policy- and decision-makers with a greater situational awareness regarding their potential risk. Yet most of those models ignore the tail risk of contagious diseases, use point forecasts, and the reliability of their parameters is rarely questioned and incorporated in the projections. We argue that a natural and empirically correct framework for assessing (and managing) the real risk of pandemics is provided by extreme value theory (EVT), an approach that has historically been developed to treat phenomena in which extremes (maxima or minima) and not averages play the role of the protagonist, being the fundamental source of risk. By analysing data for pandemic outbreaks spanning over the past 2500 years, we show that the related distribution of fatalities is strongly fat-tailed, suggesting a tail risk that is unfortunately largely ignored in common epidemiological models. We use a dual distribution method, combined with EVT, to extract information from the data that is not immediately available to inspection. To check the robustness of our conclusions, we stress our data to account for the imprecision in historical reporting. We argue that our findings have significant implications, including on the extent to which compartmental epidemiological models and similar approaches can be relied upon for making policy decisions. This Perspective argues that an approach called extreme value theory is appropriate for understanding the so-called tail risk of epidemic outbreaks, in particular by demonstrating that the distribution of fatalities due to epidemic outbreaks over the past 2500 years is fat-tailed and dominated by extreme events.
journaltitleNat. Phys.
titleTail risk of contagious diseases

journaltitleAmerican Journal of Epidemiology
titleA new framework and software to estimate time-varying reproduction numbers during epidemics
Incidence,Influenza,Measles,Reproduction number,SARS,Smallpox,Software

abstractThe quantification of transmissibility during epidemics is essential to designing and adjusting public health responses. Transmissibility can be measured by the reproduction number R, the average number of secondary cases caused by an infected individual. Several methods have been proposed to estimate R over the course of an epidemic; however, they are usually difficult to implement for people without a strong background in statistical modeling. Here, we present a ready-to-use tool for estimating R from incidence time series, which is implemented in popular software including Microsoft Excel (Microsoft Corporation, Redmond, Washington). This tool produces novel, statistically robust analytical estimates of R and incorporates uncertainty in the distribution of the serial interval (the time between the onset of symptoms in a primary case and the onset of symptoms in secondary cases). We applied the method to 5 historical outbreaks; the resulting estimates of R are consistent with those presented in the literature. This tool should help epidemiologists quantify temporal changes in the transmission intensity of future epidemics by using surveillance data.
journaltitleAm. J. Epidemiol.
titleA new framework and software to estimate time-varying reproduction numbers during epidemics
Incidence; Influenza; Measles; Reproduction number; SARS; Smallpox; Software

Elsevier BV

journaltitleJ. Math. Anal. Appl.
titleA general age-dependent branching process. II

journaltitleJournal of Mathematical Analysis and Applications
titleA general age-dependent branching process. I

abstractCases of SARS-CoV-2 infection in Manaus, Brazil, resurged in late 2020, despite previously high levels of infection. Genome sequencing of viruses sampled in Manaus between November 2020 and January 2021 revealed the emergence and circulation of a novel SARS-CoV-2 variant of concern. Lineage P.1, acquired 17 mutations, including a trio in the spike protein (K417T, E484K and N501Y) associated with increased binding to the human ACE2 receptor. Molecular clock analysis shows that P.1 emergence occurred around mid-November 2020 and was preceded by a period of faster molecular evolution. Using a two-category dynamical model that integrates genomic and mortality data, we estimate that P.1 may be 1.7-2.4-fold more transmissible, and that previous (non-P.1) infection provides 54-79% of the protection against infection with P.1 that it provides against non-P.1 lineages. Enhanced global genomic surveillance of variants of concern, which may exhibit increased transmissibility and/or immune evasion, is critical to accelerate pandemic responsiveness.
titleGenomics and epidemiology of the P.1 SARS-CoV-2 lineage in Manaus, Brazil

journaltitleThe Annals of Mathematical Statistics
titleOn the Integral Equation of Renewal Theory

Imperial College London

abstractThe global impact of COVID-19 has been profound, and the public health threat it represents is the most serious seen in a respiratory virus since the 1918 H1N1 influenza pandemic. Here we present the results of epidemiological modelling which has informed policymaking in the UK and other countries in recent weeks. In the absence of a COVID-19 vaccine, we assess the potential role of a number of public health measures – so-called non-pharmaceutical interventions (NPIs) – aimed at reducing contact rates in the population and thereby reducing transmission of the virus. In the results presented here, we apply a previously published microsimulation model to two countries: the UK (Great Britain specifically) and the US. We conclude that the effectiveness of any one intervention in isolation is likely to be limited, requiring multiple interventions to be combined to have a substantial impact on transmission. Two fundamental strategies are possible: (a) mitigation, which focuses on slowing but not necessarily stopping epidemic spread – reducing peak healthcare demand while protecting those most at risk of severe disease from infection, and (b) suppression, which aims to reverse epidemic growth, reducing case numbers to low levels and maintaining that situation indefinitely. Each policy has major challenges. We find that that optimal mitigation policies (combining home isolation of suspect cases, home quarantine of those living in the same household as suspect cases, and social distancing of the elderly and others at most risk of severe disease) might reduce peak healthcare demand by 2/3 and deaths by half. However, the resulting mitigated epidemic would still likely result in hundreds of thousands of deaths and health systems (most notably intensive care units) being overwhelmed many times over. For countries able to achieve it, this leaves suppression as the preferred policy option. We show that in the UK and US context, suppression will minimally require a combination of social distancing of the entire population, home isolation of cases and household quarantine of their family members. This may need to be supplemented by school and university closures, though it should be recognised that such closures may have negative impacts on health systems due to increased absenteeism. The major challenge of suppression is that this type of intensive intervention package – or something equivalently effective at reducing transmission – will need to be maintained until a vaccine becomes available (potentially 18 months or more) – given that we predict that transmission will quickly rebound if interventions are relaxed. We show that intermittent social distancing – triggered by trends in disease surveillance – may allow interventions to be relaxed temporarily in relative short time windows, but measures will need to be reintroduced if or when case numbers rebound. Last, while experience in China and now South Korea show that suppression is possible in the short term, it remains to be seen whether it is possible long-term, and whether the social and economic costs of the interventions adopted thus far can be reduced.
titleReport 9: Impact of non-pharmaceutical interventions (NPIs) to reduce COVID19 mortality and healthcare demand
COVID19; Coronavirus; Non-pharmaceutical interventions; healthcare demand; Mortality; Report

abstractHighly pathogenic H5N1 influenza A viruses are now endemic in avian populations in Southeast Asia, and human cases continue to accumulate. Although currently incapable of sustained human-to-human transmission, H5N1 represents a serious pandemic threat owing to the risk of a mutation or reassortment generating a virus with increased transmissibility. Identifying public health interventions that might be able to halt a pandemic in its earliest stages is therefore a priority. Here we use a simulation model of influenza transmission in Southeast Asia to evaluate the potential effectiveness of targeted mass prophylactic use of antiviral drugs as a containment strategy. Other interventions aimed at reducing population contact rates are also examined as reinforcements to an antiviral-based containment policy. We show that elimination of a nascent pandemic may be feasible using a combination of geographically targeted prophylaxis and social distancing measures, if the basic reproduction number of the new virus is below 1.8. We predict that a stockpile of 3 million courses of antiviral drugs should be sufficient for elimination. Policy effectiveness depends critically on how quickly clinical cases are diagnosed and the speed with which antiviral drugs can be distributed.
titleStrategies for containing an emerging influenza pandemic in Southeast Asia
family=Imperial College COVID-19 Response Team,

abstractFollowing the detection of the new coronavirus1 severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and its spread outside of China, Europe has experienced large epidemics of coronavirus disease 2019 (COVID-19). In response, many European countries have implemented non-pharmaceutical interventions, such as the closure of schools and national lockdowns. Here we study the effect of major interventions across 11 European countries for the period from the start of the COVID-19 epidemics in February 2020 until 4 May 2020, when lockdowns started to be lifted. Our model calculates backwards from observed deaths to estimate transmission that occurred several weeks previously, allowing for the time lag between infection and death. We use partial pooling of information between countries, with both individual and shared effects on the time-varying reproduction number (Rt). Pooling allows for more information to be used, helps to overcome idiosyncrasies in the data and enables more-timely estimates. Our model relies on fixed estimates of some epidemiological parameters (such as the infection fatality rate), does not include importation or subnational variation and assumes that changes in Rt are an immediate response to interventions rather than gradual changes in behaviour. Amidst the ongoing pandemic, we rely on death data that are incomplete, show systematic biases in reporting and are subject to future consolidation. We estimate that-for all of the countries we consider here-current interventions have been sufficient to drive Rt below 1 (probability Rt < 1.0 is greater than 99%) and achieve control of the epidemic. We estimate that across all 11 countries combined, between 12 and 15 million individuals were infected with SARS-CoV-2 up to 4 May 2020, representing between 3.2% and 4.0% of the population. Our results show that major non-pharmaceutical interventions-and lockdowns in particular-have had a large effect on reducing transmission. Continued intervention should be considered to keep transmission of SARS-CoV-2 under control.
titleEstimating the effects of non-pharmaceutical interventions on COVID-19 in Europe

journaltitlePLoS ONE
titleEstimating individual and household reproduction numbers in an emerging epidemic

Chapman & Hall/CRC

abstractIncorporating new and updated information, this second edition of THE bestselling text in Bayesian data analysis continues to emphasize practice over theory, describing how to conceptualize, perform, and critique statistical analyses from a Bayesian perspective. Its world-class authors provide guidance on all aspects of Bayesian data analysis and include examples of real statistical analyses, based on their own research, that demonstrate how to solve complicated problems. Changes in the new edition include: ·Stronger focus on MCMC·Revision of the computational advice in Part III·New chapters on nonlinear models and decision analysis·Several additional applied examples from the authors' recent research·Additional chapters on current models for Bayesian data analysis such as nonlinear models, generalized linear mixed models, and more·Reorganization of chapters 6 and 7 on model checking and data collectionBayesian computation is currently at a stage where there are many reasonable ways to compute any given posterior distribution. However, the best approach is not always clear ahead of time. Reflecting this, the new edition offers a more pluralistic presentation, giving advice on performing computations from many perspectives while making clear the importance of being aware that there are different ways to implement any given iterative simulation computation. The new approach, additional examples, and updated information make Bayesian Data Analysis an excellent introductory text and a reference that working scientists will use throughout their professional life.
titleBayesian Data Analysis

Springer Berlin

titleThe theory of branching processes

journaltitleJournal of Machine Learning Research
titleStochastic variational inference
Bayesian inference,Bayesian nonparametrics,Stochastic optimization,Topicmodels,Variational inference

journaltitleJ. R. Soc. Med.
titleThe SARS epidemic in Hong Kong: what lessons have we learned?

abstract1) A mathematical investigation has been made of the progress of an epidemic in a homogeneous population. It has been assumed that complete immunity is conferred by a single attack, and that an individual is not infective at the moment at which he receives infection. With these reservations the problem has been investigated in its most general aspects, and the following conclusions have been arrived at. (2) In general a threshold density of population is found to exist, which depends upon the infectivity, recovery and death rates peculiar to the epidemic. No epidemic can occur if the population density is below this threshold value. (3) Small increases of the infectivity rate may lead to large epidemics; also, if the population density slightly exceeds its threshold value the effect of an epidemic will be to reduce the density as far below the threshold value as initially it was above it. (4) An epidemic, in general, comes to an end, before the susceptible population has been exhausted. (5) Similar results are indicated for the case in which transmission is through an intermediate host.
journaltitleProc. R. Soc. Lond. A Math. Phys. Sci.
titleA contribution to the mathematical theory of epidemics

journaltitleAdvances in Applied Probability
titleThe point-process approach to age- and time-dependent branching processes

abstractThe sources and characters of uncertainties in engineering modeling for risk and reliability analyses are discussed. While many sources of uncertainty may exist, they are generally categorized as either aleatory or epistemic. Uncertainties are characterized as epistemic, if the modeler sees a possibility to reduce them by gathering more data or by refining models. Uncertainties are categorized as aleatory if the modeler does not foresee the possibility of reducing them. From a pragmatic standpoint, it is useful to thus categorize the uncertainties within a model, since it then becomes clear as to which uncertainties have the potential of being reduced. More importantly, epistemic uncertainties may introduce dependence among random events, which may not be properly noted if the character of uncertainties is not correctly modeled. Influences of the two types of uncertainties in reliability assessment, codified design, performance-based engineering and risk-based decision-making are discussed. Two simple examples demonstrate the influence of statistical dependence arising from epistemic uncertainties on systems and time-variant reliability problems.
journaltitleStruct. Saf.
titleAleatory or epistemic? Does it matter?
Aleatory; Epistemic; Ergodicity; Parameter uncertainty; Predictive models; Probability distribution choice; Statistical dependence; Systems; Time-variant reliability; Uncertainty

journaltitleLancet Infect Dis
titleEarly dynamics of transmission and control of {COVID}-19: a mathematical modelling study

abstractSevere acute respiratory syndrome (SARS) is a recently described illness of humans that has spread widely over the past 6 months. With the use of detailed epidemiologic data from Singapore and epidemic curves from other settings, we estimated the reproductive number for SARS in the absence of interventions and in the presence of control efforts. We estimate that a single infectious case of SARS will infect about three secondary cases in a population that has not yet instituted control measures. Public-health efforts to reduce transmission are expected to have a substantial impact on reducing the size of the epidemic.
titleTransmission dynamics and control of severe acute respiratory syndrome

journaltitleJ. Travel Med.
titleThe reproductive number of COVID-19 is higher compared to SARS coronavirus
COVID-19; China; Coronavirus; Epidemic potential; Public health emergency of international concern; R0; SARS; Wuhan; R 0

abstractPopulation-level analyses often use average quantities to describe heterogeneous systems, particularly when variation does not arise from identifiable groups. A prominent example, central to our current understanding of epidemic spread, is the basic reproductive number, R(0), which is defined as the mean number of infections caused by an infected individual in a susceptible population. Population estimates of R(0) can obscure considerable individual variation in infectiousness, as highlighted during the global emergence of severe acute respiratory syndrome (SARS) by numerous 'superspreading events' in which certain individuals infected unusually large numbers of secondary cases. For diseases transmitted by non-sexual direct contacts, such as SARS or smallpox, individual variation is difficult to measure empirically, and thus its importance for outbreak dynamics has been unclear. Here we present an integrated theoretical and statistical analysis of the influence of individual variation in infectiousness on disease emergence. Using contact tracing data from eight directly transmitted diseases, we show that the distribution of individual infectiousness around R(0) is often highly skewed. Model predictions accounting for this variation differ sharply from average-based approaches, with disease extinction more likely and outbreaks rarer but more explosive. Using these models, we explore implications for outbreak control, showing that individual-specific control measures outperform population-wide measures. Moreover, the dramatic improvements achieved through targeted control policies emphasize the need to identify predictive correlates of higher infectiousness. Our findings indicate that superspreading is a normal feature of disease spread, and to frame ongoing discussion we propose a rigorous definition for superspreading events and a method to predict their frequency.
titleSuperspreading and the effect of individual variation on disease emergence

Association for Computing Machinery

abstractSince its introduction in the early part of the nineteenth century, the theory of complex variables has played a steadily increasing role in mathematics, and in scientific research. In some fields complex algebra is used to simplify the description of a physical system. The use of a complex impedance Z in network theory is an example of this. In other fields complex algebra seems to be a basic ingredient of the physical laws. In Wave Mechanics for example a probability density P(x,t) is related to the square modulus of a wave function $\psi$(x,t) which is itself complex, being obtained from a wave equation whose coefficients may be complex.In mathematical research itself, it is rare to find a topic which is naturally restricted to real variables, and in many topics the extension to complex variables results in a simpler theory. For example a polynomial of degree n has exactly n zeros in the field of complex numbers.
booktitleProceedings of the 1967 22nd national conference
seriesACM '67
titleNumerical algorithms based on the theory of complex variable

journaltitleInfectious Disease Modelling
titleA primer on the use of probability generating functions in infectious disease modeling

abstractRenewal processes are a popular approach used in modelling infectious disease outbreaks. In a renewal process, previous infections give rise to future infections. However, while this formulation seems sensible, its application to infectious disease can be difficult to justify from first principles. It has been shown from the seminal work of Bellman and Harris that the renewal equation arises as the expectation of an age-dependent branching process. In this paper we provide a detailed derivation of the original Bellman Harris process. We introduce generalisations, that allow for time-varying reproduction numbers and the accounting of exogenous events, such as importations. We show how inference on the renewal equation is easy to accomplish within a Bayesian hierarchical framework. Using off the shelf MCMC packages, we fit to South Korea COVID-19 case data to estimate reproduction numbers and importations. Our derivation provides the mathematical fundamentals and assumptions underpinning the use of the renewal equation for modelling outbreaks.
titleOn the derivation of the renewal equation from an age-dependent branching process: an epidemic modelling perspective

journaltitleMathematics: frontiers and perspectives
titleThe dawning of the age of stochasticity

abstractMost popular statistical models in epidemic evolution focus on the dynamics of average relevant quantities and overlooks the role of small fluctuations on the model parameters. Models for Covid-19 are no exception. In this paper we show that the role of time-correlated fluctuations, far from being negligible, can in fact determine the spreading of an epidemic and, most importantly, the resurgence of the exponential diffusion in the presence of time-limited episodes in promiscuity behaviours. The results found in this work are not only relevant and specific for the Covid-19 epidemic but are more general and can be applied to other epidemics.
journaltitleSci. Rep.
titleRole of fluctuations in epidemic resurgence after a lockdown

abstractA simple and efficient method of simulation is discussed for point processes that are specified by their conditional intensities. The method is based on the thinning algorithm which was introduced recently by Lewis and Shedler for the simulation of nonhomogeneous Poisson processes. Algorithms are given for past dependent point processes containing multivariate processes. The simulations are performed for some parametric conditional intensity functions, and the accuracy of the simulated data is demonstrated by the likelihood ratio test and the minimum Akaike information criterion (AIC) procedure.
journaltitleIEEE Trans. Inf. Theory
titleOn Lewis' simulation method for point processes

abstractRenewal equations are a popular approach used in modelling the number of new infections, i.e., incidence, in an outbreak. We develop a stochastic model of an outbreak based on a time-varying variant of the Crump-Mode-Jagers branching process. This model accommodates a time-varying reproduction number and a time-varying distribution for the generation interval. We then derive renewal-like integral equations for incidence, cumulative incidence and prevalence under this model. We show that the equations for incidence and prevalence are consistent with the so-called back-calculation relationship. We analyse two particular cases of these integral equations, one that arises from a Bellman-Harris process and one that arises from an inhomogeneous Poisson process model of transmission. We outline an argument to show that the incidence integral equations that arise from both of these specific models agree with the renewal equation used ubiquitously in infectious disease modelling. We present a numerical discretisation scheme to solve these equations, and use this scheme to estimate rates of transmission from both a Bellman-Harris process and Poisson process for Influenza, Measles, Smallpox and SARS incidence.
titleUnifying incidence and prevalence under a time-varying general branching process

abstractThe effective reproduction number, Rt, is a key time-varying prognostic for the growth rate of any infectious disease epidemic. Significant changes in Rt can forewarn about new transmissions within a population or predict the efficacy of interventions. Inferring Rt reliably and in real-time from observed time-series of infected (demographic) data is an important problem in population dynamics. The renewal or branching process model is a popular solution that has been applied to Ebola and Zika virus disease outbreaks, among others, and is currently being used to investigate the ongoing COVID-19 pandemic. This model estimates Rt using a heuristically chosen piecewise function. While this facilitates real-time detection of statistically significant Rt changes, inference is highly sensitive to the function choice. Improperly chosen piecewise models might ignore meaningful changes or over-interpret noise-induced ones, yet produce visually reasonable estimates. No principled piecewise selection scheme exists. We develop a practical yet rigorous scheme using the accumulated prediction error (APE) metric from information theory, which deems the model capable of describing the observed data using the fewest bits as most justified. We derive exact posterior prediction distributions for infected population size and integrate these within an APE framework to obtain an exact and reliable method for identifying the piecewise function best supported by available epidemic data. We find that this choice optimises short-term prediction accuracy and can rapidly detect salient fluctuations in Rt, and hence the infected population growth rate, in real-time over the course of an unfolding epidemic. Moreover, we emphasise the need for formal selection by exposing how common heuristic choices, which seem sensible, can be misleading. Our APE-based method is easily computed and broadly applicable to statistically similar models found in phylogenetics and macroevolution, for example. Our results explore the relationships among estimate precision, forecast reliability and model complexity.
journaltitlePLoS Comput. Biol.
titleUsing information theory to optimise epidemic models for real-time prediction and estimation

abstractAs countries in Europe gradually relaxed lockdown restrictions after the first wave, test-trace-isolate strategies became critical to maintain the incidence of coronavirus disease 2019 (COVID-19) at low levels1,2. Reviewing their shortcomings can provide elements to consider in light of the second wave that is currently underway in Europe. Here we estimate the rate of detection of symptomatic cases of COVID-19 in France after lockdown through the use of virological3 and participatory syndromic4 surveillance data coupled with mathematical transmission models calibrated to regional hospitalizations2. Our findings indicate that around 90,000 symptomatic infections, corresponding to 9 out 10 cases, were not ascertained by the surveillance system in the first 7 weeks after lockdown from 11 May to 28 June 2020, although the test positivity rate did not exceed the 5% recommendation of the World Health Organization (WHO)5. The median detection rate increased from 7% (95% confidence interval, 6-8%) to 38% (35-44%) over time, with large regional variations, owing to a strengthening of the system as well as a decrease in epidemic activity. According to participatory surveillance data, only 31% of individuals with COVID-19-like symptoms consulted a doctor in the study period. This suggests that large numbers of symptomatic cases of COVID-19 did not seek medical advice despite recommendations, as confirmed by serological studies6,7. Encouraging awareness and same-day healthcare-seeking behaviour of suspected cases of COVID-19 is critical to improve detection. However, the capacity of the system remained insufficient even at the low epidemic activity achieved after lockdown, and was predicted to deteriorate rapidly with increasing incidence of COVID-19 cases. Substantially more aggressive, targeted and efficient testing with easier access is required to act as a tool to control the COVID-19 pandemic. The testing strategy will be critical to enable partial lifting of the current restrictive measures in Europe and to avoid a third wave.
titleUnderdetection of cases of COVID-19 in France threatens epidemic control

journaltitleVirological. org
titlePreliminary analysis of SARS-CoV-2 importation & establishment of UK transmission lineages

journaltitleArXiv e-prints
titleOn the predictability of infectious disease outbreaks
Physics - Physics and Society,Quantitative Biology - Populations and Evolution

Springer Science
Business Media LLC

abstractEuropean governments use non-pharmaceutical interventions (NPIs) to control resurging waves of COVID-19. However, they only have outdated estimates for how effective individual NPIs were in the first wave. We estimate the effectiveness of 17 NPIs in Europe's second wave from subnational case and death data by introducing a flexible hierarchical Bayesian transmission model and collecting the largest dataset of NPI implementation dates across Europe. Business closures, educational institution closures, and gathering bans reduced transmission, but reduced it less than they did in the first wave. This difference is likely due to organisational safety measures and individual protective behaviours-such as distancing-which made various areas of public life safer and thereby reduced the effect of closing them. Specifically, we find smaller effects for closing educational institutions, suggesting that stringent safety measures made schools safer compared to the first wave. Second-wave estimates outperform previous estimates at predicting transmission in Europe's third wave.
journaltitleNat. Commun.
titleUnderstanding the effectiveness of government interventions against the resurgence of COVID-19 in Europe

abstractThe time between the infection of a primary case and one of its secondary cases is called a generation time. The distribution (and mean) of the generation times is derived for a rather general class of epidemic models. The relation to assumptions on distributions of latency times and infectious times or more generally on random time varying infectiousness, is investigated. Serial times, defined as the times between occurrence of observable events in the progress of an infectious disease (e.g., the onset of clinical symptoms), are also considered.
journaltitleMath. Biosci.
titleA note on generation times in epidemic models

journaltitleLancet Infect Dis
titleEstimates of the severity of {COVID}-19 disease
volumein press

journaltitleBMC Infect. Dis.
titleLessons from a decade of individual-based models for infectious disease transmission: a systematic review (2006-2015)
Agent-based; Dynamics; Emerging diseases; Endemic diseases; Individual-based; Mathematical epidemiology; Modeling; Networks; ODD protocol; Transmission

journaltitleProceedings of the National Academy of Sciences
titleEvidence that coronavirus superspreading is fat-tailed

abstractFrom an analysis of the distributions of measures of transmission rates among hosts, we identify an empirical relationship suggesting that, typically, 20% of the host population contributes at least 80% of the net transmission potential, as measured by the basic reproduction number, R0. This is an example of a statistical pattern known as the 20/80 rule. The rule applies to a variety of disease systems, including vector-borne parasites and sexually transmitted pathogens. The rule implies that control programs targeted at the “core” 20% group are potentially highly effective and, conversely, that programs that fail to reach all of this group will be much less effective than expected in reducing levels of infection in the population as a whole.
journaltitleProc. Natl. Acad. Sci. U. S. A.
titleHeterogeneities in the transmission of infectious agents: implications for the design of control programs

Supplementary Information
Intrinsic Randomness in Epidemic Modelling Beyond Statistical Uncertainty

§ SUMMARY

This supplement provides full derivations of the results from the main text. The results are, as in the main text, presented for an epidemic occurring in continuous time, although some additional results on discrete epidemics are given in the final note of this supplement. The supplement is structured as follows.

* The first note, “Modelling”, provides a precise definition of the branching process model used throughout the paper.

* The second note, “Probability generating functions” derives probability generating functions (pgfs) for prevalence and cumulative incidence. It also discusses their efficient solution, including some special cases in which one can speed up the solution process

* The third note, “Properties of the prevalence variance”, derives the equation for the variance (via the previously derived equations for the pgf) and explores its properties, providing explanations for the various terms and proving that the prevalence of new infections is (under a mild condition on the possible spread of the epidemic) overdispersed.

* The fourth note, “Likelihood functions” contains the derivations of the pgf of the infection event times and the likelihood function presented in the main text.

* The fifth note, “Assessing future variance during an epidemic” derives the equation for variance of future cases when the cumulative incidence is known at some point in time.

* Finally, the sixth note, “Discrete epidemics” provides a range of similar results in the discrete setting, and shows the convergence of the pgf to its continuous equivalent as the step-size tends to zero.

§ BACKGROUND LITERATURE ON RENEWAL EQUATIONS

A common approach to modelling infectious diseases is to use the renewal equation. The early theory on the properties of the renewal equation can be found here <cit.>. Epidemiologically derived descriptions can be found here <cit.> where the renewal equation is framed in an epidemiological framework with reference to infection processes. The link between the renewal equation and the popular susceptible-infected-recovered models can be found here <cit.>. The basics of branching processes can be found here <cit.>. In what follows, we will arrive at a renewal equation from first principles by first starting with the probability generating function of a general branching process.

§ MODELLING

§.§ Branching process framework

We present a general time-varying age-dependent branching process that is most similar to the general branching process initially proposed by Crump, Mode and Jagers <cit.>. Following <cit.>, in our process, we begin with a single individual infected at some time $l$ whose infectious period is a random variable distributed by cumulative distribution function $G(\cdot, l)$, admitting a probability density $g(\cdot,l)$. During this individual's life length, the individual gives rise to an integer-valued random number of secondary infections according to a counting processes $\{N(t,l)\}_{t\geq l}$ ($\{N(t,l)\}$ is the number of secondary infections) where $t$ is a global “calendar” time. The amount of time for which the individual has been infected before time $t$ is therefore $t-l$.

For each infection event time - that is, for each $v$ such that

\begin{equation}
v \in\bigg\{ u \leq t : \lim_{s \to u_{-}}(N(s,l)) \neq \lim_{s \to u_{+}}(N(s,l)) \bigg\}
\end{equation}

we then define a random variable

\begin{equation}
Y(v,l) :=  \lim_{s \to v_{+}}(N(s,l)) - \lim_{s \to v_{-}}(N(s,l))
\end{equation}

to be the size of the infection event at time $v$; that is, this is the number of individuals that are infected (by the initial individual) at time $v$. Throughout this paper, it will be assumed that $Y = Y(v)$, so that $Y$ does not depend on the length of time for which an individual has been infected. However, this assumption could be removed from the model if desired.

Each newly infected individual then proceeds, independently, in the same way as the initial individual. The only change is that the time at which they are infected will be different (but, for example, the infection tree rooted at an individual infected at time $s > l$ is equal in distribution to the full infection tree if one started an epidemic with $l = s$). This self-similarity property underpins the derivations in the subsequent notes, as it allows an epidemic to be characterised purely by the “first generation” of infected individuals (and hence, the equations are derived using the “first generation principle”).

§.§ The counting process,

Our framework relies on the assumption that the counting processes $N(t,l)$ has independent increments and is continuous in probability:

\begin{equation}
\lim_{\delta \to 0}\bigg[\mathbb{P}\bigg(N(t+\delta,l) - N(t,l)\bigg)\bigg] = 0 \quad \forall t\geq l \geq 0
\end{equation}

This condition excludes any discrete formulations of the epidemic process. It will be shown later in the supplement that discrete epidemics (which are not continuous in probability), are structurally different as extra terms appear in the equations for the pgf. However, the equations in the continuous case are recovered as the step-size of the discrete process tends to zero.

A further assumption on $N(t,l)$ is that it can be constructed from a Lévy Process - that is, there is some non-negative rate function $r(t,l)$ and some Lévy Process $\mathcal{N}(t)$ such that

\begin{equation}
N(t,l) = \mathcal{N}\bigg(\int_l^tr(s,l)ds\bigg)
\end{equation}

Note that the counting processes relating to different individuals are independent, and hence will come from different independent copies of the base process $\mathcal{N}$.

This assumption is important because it means that the counting process of “infection events“ (that is, points in time such that the value of $N(t,l)$ changes) is an inhomogeneous Poisson Process, which can be shown as follows. Consider a counting process, $J_\mathcal{N}(t,l)$ that counts the increases in $\mathcal{N}$. That is,

\begin{equation}
J_{\mathcal{N}}(t) := \bigg|\bigg\{ u \leq t : \lim_{s \to u_{-}}(\mathcal{N}(s)) \neq \lim_{s \to u_{+}}(\mathcal{N}(s)) \bigg\}\bigg|
\end{equation}

where here $|\cdot|$ denotes the number of elements in a set. Then, as $\mathcal{N}$ is a Lévy Process, $J_{\mathcal{N}}(t)$ has iid (independent and identically distributed) increments and is non-decreasing in $t$ with jumps of size 1 and thus follows a Poisson Process with some rate $\kappa$ <cit.>. Thus, if $J(t,l)$ is the counting process of infection events in $\mathcal{N}(t,l)$, then

\begin{equation}
J(t,l) = J_{\mathcal{N}}\bigg(\int_l^tr(s,l)ds\bigg)
\end{equation}

and hence, $J(t,l)$ is an inhomogeneous Poisson Process with rate $\kappa r(t,l)$ as required. In particular, defining

\begin{equation}
\lambda(t,l) := \int_l^{t} r(s,l)ds,
\end{equation}

$J(t,l)$ has a generating function of

\begin{equation}
\mathcal{J}_{(t,l)}(s) = e^{\kappa\lambda(t,l)(s-1)}
\end{equation}

§.§ The rate function,

Throughout the examples in this paper, the rate function $r(t,l)$ will be given as

\begin{equation}
r(t,l) = \rho(t)\nu(t-l)
\end{equation}

Here, $\rho(t)$ is a population-level infection event rate. Note that, because the number of infections caused at each infection rate may be greater than 1 (that is one may have $J(t,l) < N(t,l)$), $\rho(t)$ cannot necessarily be interpreted in direct analogue to the reproduction number. $\nu(t-l)$ gives the infectiousness of an individual after it has been infected for time $(t-l)$. It will be assumed that $\int_0^{\infty}\nu(s)ds = 1$ so that it $\rho$ can be interpreted as the infection event rate.

§.§ Smoothness assumptions

Note that, throughout the derivations of this paper, the smoothness of $\rho$, $\nu$ and $g$ will not be explicitly considered when taking limits - it will be assumed that they are sufficiently smooth for “natural” results to hold. The authors believe that the results of this paper will hold for any piecewise continuous choices for these functions, although more detailed analysis would be needed to provide a rigorous proof of this. It is possible that they hold for much wider classes of functions, but this seems to the authors to be outside the realm of epidemiological interest, as it appears implausible that any of these functions would not be piecewise continuous in a realistic setting.

Moreover, it will be assumed that unique solutions to the equations for the pgf, mean and variance exist. Again, a proof of this property is beyond the scope of this work, although the classes of equations presented in this paper are common across the literature, and it is likely that interested readers with a pure mathematical background could find applicable results to address this issue.

§.§ Special cases for

Throughout this paper, two special cases for $N(t,l)$ are considered - the case where $N(t,l)$ is itself an inhomogeneous Poisson Process, and the case where $N(t,l)$ is a Negative Binomial process. These were used to construct the figures in the paper and explanations as to how they can be used will be presented throughout this supplement.

§ PROBABILITY GENERATING FUNCTIONS

§.§ General case

Define $F(t,l;s):=E\bigg(s^{Z(t,l)}\bigg)$ to be the generating function of $Z(t,l)$. For simplicity of notation the dependence of $F$ on $s$ will be suppressed.

To derive the generating function $F(t,l)$, we condition on  the infection period (lifetime) of the initial case, $L$.

\begin{align}
E\bigg(s^{Z(t,l)}\bigg) &= \int_0^{\infty}E\bigg(s^{Z(t,l)}\bigg|L = u\bigg) g(u,l)du\\
\label{eq:pgf_splitint} &= \int_{t-l}^{\infty}E\bigg(s^{Z(t,l)}\bigg|L = u\bigg) g(u,l)du + \int_0^{t-l}E\bigg(s^{Z(t,l)}\bigg|L = u\bigg) g(u,l)du
\end{align}

The counting process of the first individual, $N(t,l)$ is independent of this first individual's infection period $L$. If $L>t-l$ then this individual is still infectious and able to infect others at time $t$. Therefore, conditional on $L>t-l$, the number of people they have infected before time $t$ is independent of $L$ (as all infections from $N(s,l)_{l \leq s \leq t}$ are counted, irrespectively of the value of $L$). That is (the first term in Equation 11)

\begin{equation}
\int_{t-l}^{\infty} E\bigg(s^{Z(t,l)}\bigg|L = u\bigg) g(u,l)du =  \int_{t-l}^{\infty} E\bigg(s^{Z(t,l)}\bigg|L\geq t-l\bigg) g(u,l)du
\end{equation}

and hence, the first integral in Supplementary Equation <ref> can be simplified to give

\begin{equation}
E\bigg(s^{Z(t,l)}\bigg)  = \bigg(1-G(t-l,l)\bigg) E\bigg(s^{Z(t,l)}\bigg|L \geq t-l\bigg) + \int_0^{t-l}E\bigg(s^{Z(t,l)}\bigg|L = u\bigg) g(u,l)du
\label{eq:pgf_splitint_full}
\end{equation}

Let us consider the second part of Supplementary Equation <ref>. Suppose first that $L = u$ for some $u < t-l$ so that the index case is no longer alive at time $t$. Thus, the number of infection events caused by the index case is given by $J(l+u,l)$.

Define the set of times at which these infected events occurred to be $\{K_1,...,K_{J(l+u,l)}\}$ where here, importantly, the $K_i$ are labelled in a random order (so it is not necessarily the case that $K_1 < ... < K_{J(l+u,l)}$). As $J$ is an homogeneous Poisson Process and $N(t,l)$ is continuous in probability, the $K_i$ are therefore iid with pdf (probability density function)

\begin{equation}
f_{K}(k) = \frac{r(l+k,l)}{\int_0^{u}r(l+s,l)ds}
\end{equation}

It is perhaps helpful to note that this is the step which relies on $N$ being continuous in probability. If this were not the case and $N(t,l)$ had non-zero probability of increasing at some time $s$, then the knowledge that $K_1 = s$ would give some information about $K_2$, as the fact that $K_2 \neq s$ would change its probability distribution, meaning $K_1$ and $K_2$ would not be independent. Conversely, in the continuous case, $K_1 = s$ removes an event of zero measure from the probability space of $K_2$, and hence $K_1$ and $K_2$ are still independent.

Now, by the self-similarity property (<cit.>) we have

\begin{equation}
Z(t,l) = \sum_{i=1}^{J(l+u,l)}\sum_{j=1}^{Y(l+K_i(l+u,l))}Z_{ij}(t,l+K_i(l+u,l))
\label{eq:self_similarity}
\end{equation}

where each $Z_{ij}$ is an independent copy of $Z$ that is equal in distribution.  $Z_{ij}$ denotes the $j$th individual corresponding to infection event time $i$. The two summations, from all previous infections, sum over all the infection events and their sizes. This summation is valid as each individual behaves independently once it has been infected.

Recall that if $X_i$ are iid random variables (with a generating function, $G_X(s)$) and if $Y$ is a non-negative integer-valued random variable (again with a generating function, $G_Y(s)$), then,

\begin{equation}
E\bigg(s^{\sum_{i=1}^{Y}X_i}\bigg) = G_Y(G_X(s))
\label{eq:pgf_lemma}
\end{equation}

By defining $\mathcal{J}_{(t,l)}$ to be the generating function of $J(t,l)$, this relationship allows us to write $\mathbb{E}(s^{Z(t,l)}|L=u)$ as

\begin{align}
\mathbb{E}(s^{Z(t,l)}|L=u) &= \mathcal{J}_{(l+u,l)}\bigg(E\bigg[s^{\sum_{j=1}^{Y(l+K(l+u,l))}Z_j(t,l+K(l+u,l))}\bigg]\bigg)
\end{align}

where here, $K$ is equal in distribution to the $K_i$. Conditioning on the value of $K$,

\begin{equation}
E\bigg[s^{\sum_{j=1}^{Y(l+K)}Z_j(t,l+K)}\bigg] = \int_0^{u} E\bigg[s^{\sum_{j=1}^{Y(l+k)}Z_j(t,l+k)}\bigg]\frac{r(l+k,l)}{\lambda(l+u,l)}dk
\end{equation}

Thus, defining $\mathcal{Y}_{(l+k)}$ to be the generating function of $Y(l+k)$

\begin{equation}
E\bigg[s^{\sum_{j=1}^{Y(l+K)}Z_j(t,l+K)}\bigg] =   \int_0^{u} \mathcal{Y}_{(l+k)}(F(t,l+k))\frac{r(l+k,l)}{\lambda(l+u,l)}dk
\end{equation}

We can equivalently write this as an exponential, using the fact that $J(t,l)$ is Poisson distributed:

\begin{align}
\mathbb{E}(s^{Z(t,l)}|L=u)&= \mathcal{J}_{(l+u,l)}\bigg(\int_0^{u}  \mathcal{Y}_{(l+k)}(F(t,l+k))\frac{r(l+k,l)}{\lambda(l+u,l)}dk\bigg)\\
&=\exp\bigg[\kappa\lambda(l+u,l)\bigg(\int_0^{u}  \mathcal{Y}_{(l+k)}(F(t,l+k))\frac{r(l+k,l)}{\lambda(l+u,l)}dk - 1\bigg)\bigg]
\end{align}

An identical derivation can be performed on the first integral in Supplementary Equation <ref> (swapping $t-l$ for $u$ and multiplying by $s$ to account for the initial case, which is counted in the prevalence at $t$ when $L > t-l$), resulting in

\begin{align}
\mathbb{E}(s^{Z(t,l)}|L\geq t-l)&= s\mathcal{J}_{(t,l)}\bigg(\int_0^{t-l}  \mathcal{Y}_{(l+k)}(F(t,l+k))\frac{r(l+k,l)}{\lambda(t,l)}dk\bigg)\\
&= s\exp\bigg[\kappa\lambda(t,l)\bigg(\int_0^{t-l}  \mathcal{Y}_{(l+k)}(F(t,l+k))\frac{r(l+k,l)}{\lambda(t,l)}dk -1\bigg)\bigg]
\end{align}

and therefore, this yields an overall pgf

\begin{align}
\nonumber F(t,l) &= s\bigg(1-G(t-l,l)\bigg)\mathcal{J}_{(t,l)}\bigg(\int_0^{t-l}  \mathcal{Y}_{(l+u)}(F(t,l+u))\frac{r(l+u,l)}{\lambda(t,l)}du\bigg)... \\
&...+\int_0^{t-l}\mathcal{J}_{(l+u,l)}\bigg(\int_0^{u}  \mathcal{Y}_{(l+k)}(F(t,l+k))\frac{r(l+k,l)}{\lambda(l+u,l)}dk\bigg)g(u,l)du
\end{align}

or, equivalently

\begin{align}
\nonumber F(t,l) &= s\bigg(1-G(t-l,l)\bigg)\exp\bigg[\kappa\lambda(t,l)\bigg(\int_0^{t-l}  \mathcal{Y}_{(l+k)}(F(t,l+k))\frac{r(l+k,l)}{\lambda(t,l)}dk -1\bigg)\bigg]... \\
&...+\int_0^{t-l}\exp\bigg[\kappa\lambda(l+u,l)\bigg(\int_0^{u}  \mathcal{Y}_{(l+k)}(F(t,l+k))\frac{r(l+k,l)}{\lambda(l+u,l)}dk - 1\bigg)\bigg] g(u,l) du
\label{eq:pgf_general}
\end{align}

Note that by absorbing $\kappa$ into the rate function $r(l+k,l)$, it can be assumed that $\kappa = 1$. Intuitively this is simply scaling the probability density by the number of points.

§.§ Solving the pgf equation

Practically, one will always set $l=0$ for an epidemic, and so only the values $F(t,0)$ are directly relevant. However, it is still necessary to solve for $F(t,l)$ for $0\leq l \leq t$. In the language of PDEs (partial differential equations) and, specifically, the Cauchy problem, this can be explained by the fact that the “data curve” is the line $t=l$ (as the values of $F(t,t)$ are known to be equal to $s$) and the “characteristics” of the system are the lines $t = \text{constant}$. Thus, to calculate the value of $F(t,0)$, it is necessary to follow the characteristic from $(t,t)$ to $(t,0)$ and hence calculate $F(t,l)$ for $0\leq l \leq t$.

Hence, following <cit.>, solving Supplementary Equation <ref> can be greatly facilitated by defining an auxiliary equation $ F_c(t) = F(c,c-t)$ and allows us to write Supplementary Equation <ref> an equation in one variable. This is

\begin{align}
\nonumber F_c(t) &= s\bigg(1-G(t,l)\bigg)\mathcal{J}_{(c,c-t)}\bigg(\int_0^{t}  \mathcal{Y}_{(c-t+u)}(F_c(t-u))\frac{r(c-t+u,c-t)}{\lambda(c,c-t)}du\bigg)... \\
&...+\int_0^{t}\mathcal{J}_{(c-t+u,c-t)}\bigg(\int_0^{u}  \mathcal{Y}_{(c-t+k)}(F_c(t-k))\frac{r(c-t+k,c-t)}{\lambda(u,c-t)}dk\bigg) g(u,l) du
\end{align}

or, equivalently

\begin{align}
\nonumber F_c(t) &= s\bigg(1-G(t,l)\bigg)\exp\bigg[\lambda(c,c-t)\kappa\bigg(\int_0^{t}  \mathcal{Y}_{(c-t+u)}(F_c(t-u))\frac{r(c-t+u,c-t)}{\lambda(c,c-t)}du-1\bigg)\bigg]... \\
&...+\int_0^{t}\exp\bigg[\lambda(u,c-t)\kappa\bigg(\int_0^{u}  \mathcal{Y}_{(c-t+k)}(F_c(t-k))\frac{r(c-t+k,c-t)}{\lambda(u,c-t)}dk-1\bigg)\bigg]g(u,l) du
\end{align}

§.§ Poisson case

If $N(t,l)$ is an inhomogeneous Poisson Process, then, as the infection event size for a Poisson Process is always 1 <cit.>, one has $\mathcal{Y}_{(t)}(s) = s$. To aid understanding below in the Negative Binomial case, it is helpful to note that the Lévy Process, $\mathcal{N}$, can hence be characterised by

\begin{align}
\mathbb{P}( \mathcal{N}(t+dt) - \mathcal{N}(t) = 0) &= 1-\kappa dt \nonumber\\
\mathbb{P}(\mathcal{N}(t+dt) - \mathcal{N}(t) = 1) &= \kappa dt\nonumber\\
\mathbb{P}(\mathcal{N}(t+dt) - \mathcal{N}(t) > 1) &= o(dt)\nonumber
\end{align}

Setting $\kappa = 1$ as discussed above, the generating function equation becomes

\begin{align}
F(t,l) &= s\bigg(1-G(t-l,l)\bigg)\exp\bigg[\bigg(\int_0^{t-l}  F(t,l+k)\rho(l+k)\nu(k)dk -\lambda(t,l)\bigg)\bigg]... \\
&...+\int_0^{t-l} \exp\bigg[\bigg(\int_0^{u}  F(t,l+k)\rho(l+k)\nu(k)dk - \lambda(l+u,l)\bigg)\bigg] g(u,l) du
\end{align}

This equation can be further simplified by recalling that

\begin{align}
\lambda(t,l) &:= \int_l^{t} r(u,l) du = \int_0^{t-l}r(u+l,l) du = \int_0^{t-l} \rho(u+l)\nu(u) du
\end{align}

\begin{align}
F(t,l) &= s\bigg(1-G(t-l,l)\bigg)\exp\bigg[\bigg(\int_0^{t-l}  F(t,l+k)\rho(l+k)\nu(k)dk -\int_0^{t-l} \rho(l+k)\nu(k) dk\bigg)\bigg]... \nonumber\\
&...+\int_0^{t-l} \exp\bigg[\bigg(\int_0^{u}  F(t,l+k)\rho(l+k)\nu(k)dk - \int_0^{u} \rho(l+k)\nu(k) du\bigg)\bigg] g(u,l) du \nonumber\\
\nonumber\\
&= s\bigg(1-G(t-l,l)\bigg)\exp\bigg[\bigg(\int_0^{t-l}  \rho(l+k)\nu(k) \left(F(t,l+k)-1\right)dk \bigg)\bigg]... \nonumber\\
&...+\int_0^{t-l} \exp\bigg[\bigg(\int_0^{u}  \rho(l+k)\nu(k)dk \left(F(t,l+k)-1\right) \bigg)\bigg] g(u,l) du
\end{align}

For computational ease the auxiliary function equation is then

\begin{align}
F_c(t) &=s \bigg(1-G(t,l)\bigg)  \text{exp}\bigg[\bigg(\int_0^{t} \left(F_c(t-u)-1\right)\rho(c-t+u)\nu(u)du \bigg)\bigg]...    \nonumber
\\
&...+\int_0^{t}\text{exp}\bigg[\bigg(\int_0^u \left(F_c(t-k)-1\right)\rho(c-t+k)\nu(k)dk \bigg)\bigg] g(u,l)du
\end{align}

§.§ Inhomogeneous Negative Binomial case

Our derivation follows from the well-known relationship that the Negative Binomial distribution arises from a compound Poisson distribution. For $p\in(0,1)$ and $\phi\in\mathbb{R}^+$, if

\begin{equation}
X=\sum_{i=1}^N Y_i
\end{equation}

\begin{equation}
\end{equation}

and each $Y_i$ is independent of $N$, iid, and follows a logarithmic series distribution

\begin{equation}
\end{equation}

then the random variable $X$ is Negative Binomial distributed. This can easily be proven using pgfs. Therefore we have  $\kappa=-\ln(p)\phi$ and can calculate the pgf for $Y$ as $\mathcal{Y}(s) =  \frac{\ln(1-(1-p)s)}{\ln(p)}$. These can then be substituted into our general Supplementary Equation <ref>.

For clarity we re-derive this relationship explicitly. We have

\begin{equation}
\mathcal{N}(t) \sim \text{NB}(\phi t, p)
\end{equation}

As $M(t)$ has iid increments,

\begin{equation}
\mathbb{P}\bigg(\mathcal{N}(t+dt) - \mathcal{N}(t) = k\bigg) =  \mathbb{P}\bigg(\mathcal{N}(dt)= k\bigg) = \frac{(k+\phi dt-1)(k+\phi dt-2)...\phi dt}{k!}(1-p)^kp^{\phi dt}
\end{equation}

Thus, to leading order, for $k > 0$, one has

\begin{equation}
\mathbb{P}\bigg(\mathcal{N}(t+dt) - \mathcal{N}(t) = k\bigg) = \frac{(1-p)^k\phi dt}{k} + o(dt)
\end{equation}

while if $k=0$,

\begin{equation}
\mathbb{P}\bigg(\mathcal{N}(t+dt) - \mathcal{N}(t) = 0\bigg) = p^{\phi dt} = 1 + \ln(p)\phi dt + o(dt)
\end{equation}

(noting that $\ln(p) < 0$). This means that the infection event process $J_{\mathcal{N}}$ satisfies

\begin{equation}
\mathbb{P}\bigg(J_{\mathcal{N}}(t+dt) - J_{\mathcal{N}}(t) = 0\bigg) =  1 + \ln(p)\phi dt + o(dt)
\end{equation}

\begin{align}
\mathbb{P}\bigg(J_{\mathcal{N}}(t+dt) - J_{\mathcal{N}}(t) = 1\bigg) &=  \sum_{k=1}^{\infty}\frac{(1-p)^k\phi dt}{k} + o(dt)\\
&= - \ln(p)\phi dt + o(dt)
\end{align}

and hence, $J_{\mathcal{N}}$ is a Poisson Process of rate $-\ln(p)\phi$ <cit.> . Thus, one has

\begin{equation}
\kappa = -\ln(p)\phi
\end{equation}

as expected. Moreover, the pmf (probability mass function) of a infection event size, $Y$ is given by

\begin{equation}
\label{eq:logpdf}
\mathbb{P}(Y = k) = \frac{(1-p)^k }{- k \ln(p)}
\end{equation}

One can hence find the generating function as

\begin{equation}
\mathcal{Y}(s) = \sum_{k=1}^{\infty} \frac{((1-p)s)^k \phi}{-k \ln(p)}
\end{equation}

Noting that

\begin{equation}
\sum_{k=1}^{\infty} \frac{(1-p)^k}{-k \ln(p)} = 1
\end{equation}

one has

\begin{equation}
\mathcal{Y}(s) =  \frac{\ln(1-(1-p)s)}{\ln(p)}\sum_{k=1}^{\infty} \frac{(1 - (1 - (1-p)s))^k }{-k \ln(1-(1-p)s)} =  \frac{\ln(1-(1-p)s)}{\ln(p)}
\end{equation}

These results can be substituted into the general formula to give

\begin{align}
F(t,l) = & s\bigg(1-G(t-l,l)\bigg)\exp\bigg[ -\phi \bigg(\int_0^{t-l}\ln(1-(1-p)F(t,l+u))  \rho(u+l)\nu(u)du  + \ln(p)\lambda(t,l)\bigg)\bigg]...\nonumber\\
&...+   \int_0^{t-l}\exp\bigg[ -\phi\bigg(\int_0^{u}\ln(1-(1-p)F(t,l+k))  \rho(k+l)\nu(k)dk  + \ln(p)\lambda(u,l)\bigg)\bigg]g(u,l)du
\end{align}

As in the Poisson case, this equation can be simplified by factoring $\lambda$
\begin{align}
F(t,l) = & s\bigg(1-G(t-l,l)\bigg)\exp\bigg[ -\phi \bigg(\int_0^{t-l}\left(\ln(1-(1-p)F(t,l+u)) - \ln(p) \right)  \rho(u+l)\nu(u)du \bigg)\bigg]...\nonumber\\
&...+   \int_0^{t-l}\exp\bigg[ -\phi\bigg(\int_0^{u} \left(\ln(1-(1-p)F(t,l+k)) -\ln(p)\right) \rho(k+l)\nu(k)dk  \bigg)\bigg]g(u,l)du
\end{align}

The easier-to-solve auxiliary function is given by

\begin{align}
F_c(t) = &s\bigg(1-G(t-l,l)\bigg)\exp\bigg[ -\phi \bigg(\int_0^{t}\left(\ln(1-(1-p)F_c(t-u)) -\ln(p)\right) \rho(c-t+u)\nu(u)du  \bigg)\bigg]...\nonumber\\
&...+    \int_0^{t}\exp\bigg[ -\phi  \bigg(\int_0^{u}\left(\ln(1-(1-p)F_c(t-k)) -\ln(p)\right)  \rho(c-t+k)\nu(k)dk  \bigg)\bigg]g(u,l)du
\label{eq:pgf_negativebin_auxiliary}
\end{align}

If $p = \frac{\phi}{1 + \phi}$, then the Poisson case (with $\kappa = 1$) is recovered in the $\phi \to \infty$ limit.

Note that $\mathbb{E}[N(t,l)] = \frac{\phi \lambda(t,l) (1-p)}{p}$ while in our case, we impose that $\mathbb{E}[N(t,l)] = \lambda(t,l)$. Solving for $p$ we can see $p = \frac{\phi}{1+\phi}$ and this relation can be substituted into Supplementary Equation <ref>. Note that this agrees with the definition of $p$ in the Poisson limit.

§.§ Cumulative incidence

Similar to prevalence, cumulative incidence can be calculated by counting all previous infections as well as current ones. Following an identical derivation to prevalence the pgf for cumulative incidence simply requires multiplying the second integral by $s$ as the initial infection is counted in the cumulative incidence regardless of the value of $L$.

\begin{align}
\nonumber F(t,l) &= s\bigg(1-G(t-l,l)\bigg)\mathcal{J}_{(t,l)}\bigg(\int_0^{t-l}  \mathcal{Y}_{(l+u)}(F(t,l+u))\frac{r(l+u,l)}{\lambda(t,l)}du\bigg)... \\
&...+s\int_0^{t-l}\mathcal{J}_{(l+u,l)}\bigg(\int_0^{u}  \mathcal{Y}_{(l+k)}(F(t,l+k))\frac{r(l+k,l)}{\lambda(t,l)}dk\bigg)g(u,l)du
\end{align}

§.§ A simplified pgf ignoring

By assuming $g(u,l)=0 ~~\forall~~ u$ and therefore $G(u,l)=0~~ \forall ~~u$, the pgf for prevalence (or, in this case, equivalently, cumulative incidence) simplifies to

\begin{align}
\nonumber F(t,l) &= s\mathcal{J}_{(t,l)}\bigg(\int_0^{t-l}  \mathcal{Y}_{(l+u,l)}(F(t,l+u))\frac{r(l+u,l)}{\lambda(t,l)}du\bigg)
\label{eq:pdf_incidence_approx}
\end{align}

Additional computational savings can be gained in our case $r(t,l)=\rho(t)\nu(t-l)$ if the infectiousness $\nu$ decays to zero quickly. This means that the auxiliary equation used for computation can be truncated to some time $\text{min}(t,T)$. For example, in the Poisson case this becomes,

\begin{align}
F_c(t) &= \exp\bigg[\bigg(\int_0^{\text{min}(t,T)}  \left(F_c(t-u)-1\right)\rho(c-t+u)\nu(u)du \bigg)\bigg]s
\end{align}

and in the Negative Binomial case this becomes,

\begin{align}
F_c(t) = & \exp\bigg[ -\phi \bigg(\int_0^{\text{min}(t,T)} \left(\ln(1-(1-p)F_c(t-u)) -\ln(p)\right) \rho(c-t+u)\nu(u)du  \bigg)\bigg]s
\end{align}

These computational savings allow computation of the pgf for millions of iterations in minutes.

§.§ Calculating the probability mass function via the pgf

Following <cit.> and <cit.> (originally from <cit.>), by the properties of pgfs, the probability mass function $p$ can be recovered through a pgf $F$'s derivatives at $s=0$

\begin{align*}
\mathbb{P}(n) = \frac{1}{n!}\left(\frac{d}{ds}\right)^n F(s;t,\tau) |_{s=0}
\end{align*}

This is generally computationally intractable.
A well-known result from complex analysis <cit.> holds that

\begin{align}
f^{(n)}(a) = \frac{n!}{2\pi i} \oint \frac{f(z)}{\left(z-a\right)^{n+1}}\, dz.
\end{align}

\begin{align}
\mathbb{P}(n) = \frac{1}{2\pi i}\oint\frac{F(z;t,\tau)}{z^{n+1}} dz
\end{align}

This integral can be done on a closed circle around the origin such that $z=re^{i\theta}$ and $dz=izd\theta$ - i.e.

\begin{align}
\mathbb{P}(n) = \frac{1}{2\pi}\int_0^{2\pi}\frac{F(re^{i\theta};t,\tau)}{(re^{i\theta})^n} d\theta
\end{align}

Finally through substitution $\theta = 2\pi u$ such that $d\theta=2\pi du$, where $u\in{[0,1]}$ we find

\begin{align}
\mathbb{P}(n) = \int_0^{1}\frac{F(re^{2\pi i u};t,\tau)}{r^n e^{2\pi i u n}} du
\end{align}

Since trapezoidal sums are known to converge geometrically for periodic analytic
functions (Davis 1959) a simple approximation becomes

\begin{align}
\mathbb{P}(n) = \frac{1}{Mr^n}\sum_{m=0}^{M-1}F(re^{2\pi i m/M};t,\tau) e^{ - 2\pi i n m/M}
\end{align}

Bornemann<cit.> suggest using $r=1$.

The probability mass function for any time and $n$ can be determined numerically. One needs $M\geq n$, which requires solving $n$ renewal equations for the generating function and performing a fast Fourier transform. This is generally computationally fast, but may become slightly burdensome for epidemics with very large numbers of infected individuals.

§ PROPERTIES OF THE PREVALENCE VARIANCE

§.§ Derivation of equation for mean prevalence

Before deriving the equation for the prevalence variance, it is important to derive the equation governing the mean prevalence. This has been previously derived in <cit.>, although here, we re-derive it from our new pgfs. First note that

\begin{align}
&\frac{\partial}{\partial s}\bigg(\mathcal{J}_{(t,l)}\bigg(\int_0^{t-l}  \mathcal{Y}_{(l+u,l)}(F(t,l+u))\frac{r(l+u,l)}{\lambda(t,l)}du\bigg)\bigg)...\nonumber\\
&= \bigg[\int_0^{t-l} F_s(t,l+u)\frac{r(l+u,l)}{\lambda(t,l)} \mathcal{Y}_{(l+u,l)}'(F(t,l+u))du\bigg]\bigg[\mathcal{J}'_{(t,l)}\bigg(\int_0^{t-l}  \mathcal{Y}_{(l+u,l)}(F(t,l+u))\frac{r(l+u,l)}{\lambda(t,l)}du\bigg)\bigg)\bigg]
\end{align}

Now, setting $s = 1$ so that $F(\cdot,\cdot) = 1$ and $F_s(\cdot,\cdot) = M(\cdot,\cdot)$, one has

\begin{equation}
\bigg[\int_0^{t-l} M(t,l+u)\frac{r(l+u,l)}{\lambda(t,l)} \mathcal{Y}_{(l+u,l)}'(1)du\bigg]\bigg[\mathcal{J}'_{(t,l)}\bigg(\int_0^{t-l}  \mathcal{Y}_{(l+u,l)}(1)\frac{r(l+u,l)}{\lambda(t,l)}du\bigg)\bigg)\bigg]
\end{equation}

Now, define $B(t) = \mathbb{E}(Y(t))$ so that $\mathcal{Y}_{(l+u)}'(1) = B(l+u)$. Moreover, $\mathcal{Y}_{(l+u)}(1) = 1$ so the equation becomes

\begin{equation}
\bigg[\int_0^{t-l} M(t,l+u)\frac{r(l+u,l)}{\lambda(t,l)} B(l+u)du\bigg]\bigg[\mathcal{J}'_{(t,l)}\bigg(\int_0^{t-l} \frac{r(l+u,l)}{\lambda(t,l)}du\bigg)\bigg)\bigg]
\end{equation}

Now, necessarily

\begin{equation}
\int_0^{t-l} \frac{r(l+u,l)}{\lambda(t,l)}du = 1 \Rightarrow \mathcal{J}'_{(t,l)}\bigg(\int_0^{t-l} \frac{r(l+u,l)}{\lambda(t,l)}du\bigg) = \mathbb{E}(J(t,l)) = \kappa \lambda(t,l)
\end{equation}

and so, this results in

\begin{equation}
\int_0^{t-l} M(t,l+u)\frac{r(l+u,l)}{\lambda(t,l)} B(l+u)\kappa \lambda(t,l)du
\end{equation}

Moreover, evaluating

\begin{equation}
\mathcal{J}_{(t,l)}\bigg(\int_0^{t-l}  \mathcal{Y}_{(l+u)}(F(t,l+u))\frac{r(l+u,l)}{\lambda(t,l)}du\bigg)
\end{equation}

at $s=1$ gives

\begin{align}
\mathcal{J}_{(t,l)}\bigg(\int_0^{t-l}  \mathcal{Y}_{(l+u)}(1)\frac{r(l+u,l)}{\lambda(t,l)}du\bigg) &=   \mathcal{J}_{(t,l)}\bigg(\int_0^{t-l} 1\times \frac{r(l+u,l)}{\lambda(t,l)}du\bigg)\\
&= \mathcal{J}_{(t,l)}(1)\\
\end{align}

Thus, the derivative of the full generating function equation gives

\begin{align}
\label{eq:meanold}
M(t,l) &= (1-G(t-l,l))\bigg[1 + \int_0^{t-l} M(t,l+u)\frac{r(l+u,l)}{\lambda(t,l)} B(l+u)\kappa \lambda(t,l)du\bigg]...\\
&...+ \int_0^{t-l}\int_0^uM(t,l+k)\frac{r(l+k,l)}{\lambda(l+u,l)}B(l+k)\kappa \lambda(l+u,l)g(u,l) dk du
\end{align}

This can be simplified significantly. Note that,

\begin{align}
&\nonumber \int_0^{t-l}\int_0^uM(t,l+k)\frac{r(l+k,l)}{\lambda(l+u,l)}B(l+k)\kappa \lambda(l+u,l)g(u,l) dk du =\\ &\int_0^{t-l}\int_0^uM(t,l+k)r(l+k,k)B(l+k)\kappa g(u,l) dk du
\end{align}

Moreover, one can change the order of integration to get

\begin{equation}
\int_0^{t-l}\int_k^{t-l}M(t,l+k)r(l+k,k)B(l+k)\kappa g(u,l) du dk = \int_0^{t-l}M(t,l+k)r(l+k,k)B(l+k)\kappa (G(t-l,l) - G(k,l))
\end{equation}

and hence, one can write the equation for $M(t,l)$ as

\begin{equation}
\label{eq:Meaneq}
M(t,l) = (1-G(t-l,l)) + \int_0^{t-l} M(t,l+u)r(l+u,l) \left(B(l+u)\kappa\right) (1-G(u,l))du
\end{equation}

Note that, for the Poisson special case, $B(l+u,l) = 1$ and for the Negative Binomial special case, $B(l+u,l) = \frac{p-1}{p\ln(p)} = -\frac{1}{\ln(p)\phi}$. In both cases, it may improve the epidemiological interpretation of $\rho$ to absorb the $B(l+u,l)\kappa$ term into $\rho$ (so that $\rho$ becomes a measure of the rate of new infections). This gives the simpler equation

\begin{equation}
M(t,l) = (1-G(t-l,l)) + \int_0^{t-l} M(t,l+u)\rho(l+u)\nu(u) (1-G(u,l))du
\end{equation}

which agrees with <cit.>.

§.§ Derivation of equation for prevalence variance

The equation for variance can now be found by taking the second derivative of the pgf. Define $W(t,l) := \mathbb{E}(Z(t,l)(Z(t,l) - 1))$. Note that this then gives the variance, $V(t,l)$ as $V(t,l) = W(t,l)+M(t,l) - M(t,l)^2$.

Consider first the term

\begin{equation}
s\bigg(1-G(t-l,l)\bigg)\mathcal{J}_{(t,l)}\bigg(\int_0^{t-l}  \mathcal{Y}_{(l+u)}(F(t,l+u))\frac{r(l+u,l)}{\lambda(t,l)}du\bigg)
\end{equation}

The first derivative of this term is equal to

\begin{align}
&\bar{G}(t-l,l)\mathcal{J}_{(t,l)}\bigg(\int_0^{t-l}  \mathcal{Y}_{(l+u)}(F(t,l+u))\frac{r(l+u,l)}{\lambda(t,l)}du\bigg) + ...\nonumber \\
&s\bar{G}(t-l,l)\bigg[\int_0^{t-l}F_s(t,l+u)\mathcal{Y}'_{(l+u)}(F(t,l+u))\frac{r(l+u,l)}{\lambda(t,l)}du\bigg]\mathcal{J}'_{(t,l)}\bigg(\int_0^{t-l}  \mathcal{Y}_{(l+u)}(F(t,l+u))\frac{r(l+u,l)}{\lambda(t,l)}du\bigg)
\end{align}

Then, the second derivative is equal to

\begin{align}
&2\bar{G}(t-l,l)\bigg[\int_0^{t-l}F_s(t,l+u)\mathcal{Y}'_{(l+u)}(F(t,l+u))\frac{r(l+u,l)}{\lambda(t,l)}du\bigg]\mathcal{J}'_{(t,l)}\bigg(\int_0^{t-l}  \mathcal{Y}_{(l+u)}(F(t,l+u))\frac{r(l+u,l)}{\lambda(t,l)}du\bigg) +\nonumber\\
&+s\bar{G}(t-l,l)\bigg[\int_0^{t-l}F_s(t,l+u)\mathcal{Y}'_{(l+u)}(F(t,l+u))\frac{r(l+u,l)}{\lambda(t,l)}du\bigg]^2\mathcal{J}''_{(t,l)}\bigg(\int_0^{t-l}  \mathcal{Y}_{(l+u)}(F(t,l+u))\frac{r(l+u,l)}{\lambda(t,l)}du\bigg)\nonumber\\
&+s\bar{G}(t-l,l)\bigg[\int_0^{t-l}F_{ss}(t,l+u)\mathcal{Y}'_{(l+u)}(F(t,l+u))\frac{r(l+u,l)}{\lambda(t,l)}du\bigg]\mathcal{J}'_{(t,l)}\bigg(\int_0^{t-l}  \mathcal{Y}_{(l+u)}(F(t,l+u))\frac{r(l+u,l)}{\lambda(t,l)}du\bigg)\nonumber\\
&+s\bar{G}(t-l,l)\bigg[\int_0^{t-l}F^2_s(t,l+u)\mathcal{Y}''_{(l+u)}(F(t,l+u))\frac{r(l+u,l)}{\lambda(t,l)}du\bigg]\mathcal{J}'_{(t,l)}\bigg(\int_0^{t-l}  \mathcal{Y}_{(l+u)}(F(t,l+u))\frac{r(l+u,l)}{\lambda(t,l)}du\bigg)
\end{align}

Now, one can evaluate this as $s=1$. Note that

\begin{align}
\int_0^{t-l}  \mathcal{Y}_{(l+u)}(F(t,l+u))\frac{r(l+u,l)}{\lambda(t,l)}du &= \int_0^{t-l}  \mathcal{Y}_{(l+u)}(1)\frac{r(l+u,l)}{\lambda(t,l)}du\nonumber\\
&= \int_0^{t-l} 1 \times \frac{r(l+u,l)}{\lambda(t,l)}du\nonumber\\
&= 1
\end{align}

Moreover, define $B^W(t) := \mathbb{E}(Y(t)(Y(t) - 1))$ and $C^W(t,l) := \mathbb{E}(J(t,l)(J(t,l)-1))$. Note also $\mathbb{E}(J(t,l))=\lambda(t,l)$. Thus, the second derivative evaluated at $s=1$ is

\begin{align}
\label{eq:cumincdifference} &2\bar{G}(t-l,l)\bigg[\int_0^{t-l}M(t,l+u)B(l+u)\frac{r(l+u,l)}{\lambda(t,l)}du\bigg]\kappa \lambda(t,l)\nonumber\\
&+ \bar{G}(t-l,l)\bigg[\int_0^{t-l}W(t,l+u)B(l+u)\frac{r(l+u,l)}{\lambda(t,l)}du\bigg]\kappa \lambda(t,l)\nonumber\\
&+  \bar{G}(t-l,l)\bigg[\int_0^{t-l}M(t,l+u)^2B^W(l+u)\frac{r(l+u,l)}{\lambda(t,l)}du\bigg]\kappa \lambda(t,l)
\end{align}

Noting that $J(t,l)$ is Poisson, one has

\begin{equation}
C^W(t,l) + \mathbb{E}(J(t,l)) - \mathbb{E}(J(t,l)^2) = \text{var}(J(t,l)) = \mathbb{E}(J(t,l))
\end{equation}

and hence

\begin{equation}
C^W(t,l) = \mathbb{E}(J(t,l))^2
\end{equation}

\begin{align}
\chi(t,l,k) :=\kappa \bigg[W(t,l+k)B(l+k)r(l+k,l)+ M(t,l+k)^2B^W(l+k)r(l+k,l)\bigg]
\end{align}

Then, the same process can be carried out for the second part of the equation to give

\begin{align}
W(t,l) = &2\bar{G}(t-l,l)\bigg[\int_0^{t-l}M(t,l+u)B(l+u)r(l+u,l)du\bigg] + \bar{G}(t-l,l)\int_0^{t-l}\chi(t,l,k)dk...\nonumber\\
...&+\bar{G}(t-l,l)\bigg[\int_0^{t-l}M(t,l+u)B(l+u)\kappa r(l+u,l)du\bigg]^2...\nonumber\\
...&+ \int_0^{t-l}\bigg[\int_0^{u}M(t,l+u)B(l+u)\kappa r(l+u,l)du\bigg]^2g(u,l)du...\nonumber\\
...& +  \int_0^{t-l}\int_0^u\chi(t,l,k)dkg(u,l)du
\end{align}

For ease of notation, define

\begin{equation}
S(t,l,u) := \bigg[\int_0^{u}M(t,l+k)B(l+k)\kappa r(l+k,l)dk\bigg]^2
\end{equation}

so that

\begin{align}
W(t,l) = &2\bar{G}(t-l,l)\bigg[\int_0^{t-l}M(t,l+u)B(l+u)r(l+u,l)du\bigg] + \bar{G}(t-l,l)\int_0^{t-l}\chi(t,l,k)dk...\nonumber\\
...&+\bar{G}(t-l,l)S(t,l,t-l) + \int_0^{t-l}S(t,l,u)g(u,l)du+  \int_0^{t-l}\int_0^u\chi(t,l,k)dkg(u,l)du
\end{align}

Now, changing the order of integration in the final term (as was done in the derivation of the mean prevalence), this can be rewritten as

\begin{align}
W(t,l) &= 2\bar{G}(t-l,l)\bigg[\int_0^{t-l}\kappa M(t,l+u)B(l+u)r(l+u,l)du\bigg] +\bar{G}(t-l,l)S(t,l,t-l)...\nonumber\\
&...+ \int_0^{t-l}S(t,l,u)g(u,l)du+ \int_0^{t-l}\chi(t,l,k)\bar{G}(k,l)dk
\end{align}

From this, we can create an equation for $\mathbb{E}(Z(t,l)^2) := X(t,l) = W(t,l) + M(t,l)$ by defining

\begin{equation}
\chi^X(t,l,k) = \kappa \bigg[X(t,l+k)B(l+k,l)r(l+k,l)+ M(t,l+k)^2B^W(l+k,l)r(l+k,l)\bigg]
\end{equation}

and then simply adding the equation for $M$ to give

\begin{align}
X(t,l) &=\bar{G}(t-l,l)+ 2\bar{G}(t-l,l)\bigg[\int_0^{t-l}\kappa M(t,l+u)B(l+u)r(l+u,l)du\bigg] +\bar{G}(t-l,l)S(t,l,t-l)...\nonumber\\
\label{eq:meanadd} &...+ \int_0^{t-l}S(t,l,u)g(u,l)du+ \int_0^{t-l}\chi^X(t,l,k)\bar{G}(k,l)dk
\end{align}

Finally, to form the equation for the variance $V(t,l) = X(t,l) - M(t,l)^2$, note that

\begin{align}
\chi^X(t,l,k) &= \kappa \bigg[X(t,l+k)B(l+k)r(l+k,l)+ M(t,l+k)^2(\mathbb{E}(Y(l+k)^2) -  B(l+k))r(l+k,l)\bigg]\\
&= \kappa\bigg[V(t,l+k)B(l+k)r(l+k,l)+ M(t,l+k)^2\mathbb{E}(Y(l+k)^2)r(l+k,l)\bigg]\\
&:= \chi^V(t,l,k)
\end{align}

and hence, subtracting $M(t,l)^2$ from both sides of the equation for $X(t,l)$ gives

\begin{align}
\label{eq:Vareq}
V(t,l) &=\bar{G}(t-l,l)+ 2\bar{G}(t-l,l)\bigg[\int_0^{t-l}\kappa M(t,l+u)B(l+u)r(l+u,l)du\bigg] +\bar{G}(t-l,l)S(t,l,t-l)...\nonumber\\
&...+ \int_0^{t-l}S(t,l,u)g(u,l)du+ \int_0^{t-l}\chi^V(t,l,k)\bar{G}(k,l)dk - M(t,l)^2
\end{align}

§.§ An explanation of the variance equation

There are two main sources of uncertainty in the infection process - the infectious period of an individual, and the number and timing of infections that occur during this infectious period. One can show that the variance splits into three terms - one for each of these two sources of uncertainty from the initial individual, and one which propagates the uncertainty through the descendants of the initial individual.

Each term will be derived by assuming that all other parts of the model are deterministic. To begin, suppose that the infectious period of the initial individual is random but all other parts of the model are deterministic, so that, given that the initial individual is infectious at time $l+u$, it will infect $B(l+u)r(l+u,l)dt$ people in the interval $[u,u+dt]$ (note that this is an abstraction to illustrate the source of this variance, as it is impossible for non-integer numbers of infections to occur). Moreover, it is assumed that each of these individuals have given rise to exactly $M(t,l+u)$ infections at time $t$. Then, note that

\begin{align}
\text{var}(Z(t,l)) &= \mathbb{E}(Z(t,l)^2) - \mathbb{E}(Z(t,l))^2\\
&= \int_0^{\infty}\mathbb{E}(Z(t,l)^2 | L = u)g(u,l)du -M(t,l)^2\\
&= \int_0^{t-l}\bigg[\int_0^{u}M(t,l+k)B(l+k)r(l+k,l)dk\bigg]^2g(u,l)du ...\\
\nonumber&...+\bar{G}(t-l,l)\bigg(1 + \int_0^{t-l}M(t,l+k)B(l+k)r(l+k,l)dk\bigg)^2 - M(t,l)^2\\
&=  \bar{G}(t-l,l) + 2\bar{G}(t-l,l)\int_0^{t-l}M(t,l+k)B(l+k)r(l+k,l)dk +...\\
\nonumber&...+\bar{G}(t-l,l)S(t,l,t-l) + \int_0^{t-l}S(t,l,u)g(u,l)du - M(t,l)^2
\end{align}

which recovers all the terms of the variance equation except for $\int_0^{t-l}\chi^V(t,l,k)\bar{G}(k,l)dk$.