Unifying incidence and prevalence under a time-varying general branching process

Renewal 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 also 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 serological prevalence of SARS-CoV-2 in the UK and historical incidence data on Influenza, Measles, SARS and Smallpox.


Introduction
Mathematical descriptions of infectious disease outbreaks are fundamental to forecasting and simulating the dynamics of epidemics, as well as to understanding the mechanics of how transmission occurs.Epidemiological quantities of interest include incidence (the number of new infections at a given time point), cumulative incidence (the total number of infections up to a given time point) and prevalence (the number of infected individuals at a given time point).Taking a somewhat reductive perspective, it can be said that two main popular frameworks co-exist when modelling an infectious disease outbreak, namely, individual-based models juxtaposed with governing equations.Individual-based models are not only simple to understand in terms of their fundamental assumptions but have also proven extremely impactful [24].However, mathematical tractability is limited, reliable estimates of expectations may require millions of simulations given the fat-tailed, multiplicative nature of epidemics, and inference can be challenging, with parameter inter-dependence making sensitivity analysis unreliable.In contrast, governing equations tend to have a stronger physical interpretation, are easier to perform inference over, and can be embedded in complex models easily [26].
The most widely known set of governing equations was presented in the seminal work of Kermack and McKendrick [39], where they studied the number and distribution of infections of a transmissible disease as it progresses through a population over time.They constructed classes, called compartments, and modelled the propagation of infectious disease via interactions among these compartments.The result is the popular susceptible-infected-recovered (SIR) model, variants of which are widely used in epidemiology.Stochastic versions of SIR models, formulated either as stochastic differential equations or continuous-time Markov chains, are popular when modelling small populations or stochastic environments [2].Deterministic and stochastic SIR models provide an intuitive mechanism for understanding disease transmission, and in the original derivation of [39], they were noted to be similar to the Volterra equation [47].The Volterra equation (of the second kind), or more commonly, the renewal equation, is another popular governing equation [12,17,27,45].A large body of work in infectious disease epidemiology is based around the renewal equation and many modifications exist [1,14,28,50].There is a connection between specific compartmental models and renewal equations [13,49] but this link has not been established in full generality.The vast majority of renewal frameworks model only incidence, and the explicit link between prevalence and incidence often requires the use of a latent process for incidence [10].
Between individual-based and governing equation models are stochastic branching processes.Branching processes are applied in the modelling of epidemics by first constructing a stochastic process where infected individuals transmit disease according to simple rules, and then deriving a governing equation for the average behaviour.For example Galton-Watson processes, where individuals infect other individuals at generations specified by a fixed time, provide a tractable and intuitive way of modelling the spread of an infectious disease [3,31].In 1948, Bellman and Harris [5] elegantly captured a more complex underlying infection mechanism by formu-Figure 1: Simulation of an age-dependent Bellman-Harris branching process in terms of prevalence.Left plot shows the Monte Carlo mean (red) alongside the theoretical mean (green).Right plot shows both the Monte Carlo and theoretical mean, overlaid on the underlying 1,000 simulated trajectories (translucent black lines).In this example, the time-varying reproduction number is given by R(t) = 1.15 + sin(0.15t), while the generation interval follows the Gamma(3, 1) distribution.Algorithm 2, given below, is used to compute the theoretical mean.lating an age-dependent branching process, where the age-dependence alludes to individuals who infect other individuals after a random interval of time.Interestingly, the expectation of the Bellman-Harris process [5] follows a renewal equation, whereby their framework links the two worlds of individual-based modelling and governing equations.The age-dependence assumption of Bellman and Harris allows, in particular, for the variable time between exposure to a pathogen and subsequent transmission to be modelled more realistically, and provides a framework encoding useful biological characteristics of the infecting pathogen, such as incubation periods and non-monotonic infectiousness.Crump, Mode [18,19] and (independently) Jagers [38] further extended the Bellman-Harris process to a general branching process where individuals not only can infect at random times, but can do so randomly over the duration of their infection (as opposed to the Bellman-Harris process where all subsequent infections generated by each infected individual happen at a single random time).
The original formulation by Bellman and Harris [4], along with subsequent work by Harris [35], the work of Crump, Mode and Jagers [18,19,38] as well as the perspective of Bharucha-Reid [7], with specific application to epidemics, all focused on the simple case of a constant/basic reproduction number R 0 .The form of this renewal equation when only considering R 0 is exactly what is commonly used in epidemic modelling where the incidence of infections I(t) follows a renewal equation given by where g(•) is the probability density function (PDF) of the generation interval.Introducing a time-varying reproduction number R(t) within the Bellman-Harris process in general does not simply entail replacing R 0 with R(t) in the renewal equation.This is not possible because a history of how many secondary infections are created is needed.While justifications based on heuristic arguments such as Lotka's [43] (used in tracking the numbers of females in an agestructured population) or the one given by Fraser [27] are valid within their respective contexts, these arguments lose their validity when considering a stochastic age-dependent branching process with a time-varying reproduction process [6,40,41].Indeed, we will demonstrate that these arguments are only valid for the specific case of incidence, not for prevalence or cumulative incidence.Furthermore, to our knowledge, no one has previously investigated a fully timevarying reproduction process under the more general Crump-Mode-Jagers framework.Besides the work by Kimmel [40] on the time-varying Bellman-Harris process, a generation-dependent life length distribution within the Bellman-Harris process has been studied by Fildes [25] and a generation-dependent offspring distribution by Fearn [23].Moreover, Edler [22] and later Biggins and Götz [8] have analysed a generation-dependent reproduction process in the Crump-Mode-Jagers setting.
In this paper, we introduce an outbreak model based on a time-varying version of the Crump-Mode-Jagers process, which we formulate using random characteristics [41].Notably, Bellman-Harris, Galton-Watson and Markov branching processes are all special cases of this process.In our novel time-varying Crump-Mode-Jagers process, we specifically allow the statistical properties of infections, i.e., "offspring", generated by each individual to vary over time.Building on this model, we lay down a general, stochastic process foundation for incidence, cumulative incidence and prevalence, and characterise the renewal-like integral equations they follow.We show that the equations for prevalence and incidence are consistent with the well-known back-calculation relationship [10,20] used in infectious disease epidemiology.We also show that the common renewal equation used ubiquitously for modelling incidence [17,27] is in fact, under specific conditions, equivalent to the integral equation for incidence in our framework.Additionally, we formulate a novel reproduction process where infections occur randomly over the duration of each individual's infection according to an inhomogeneous Poisson process.The model thus eschews the common assumption that infections happen instantaneously at a random time, as in the Bellman-Harris process, but still admits analytically tractable integral equations for prevalence and incidence.Finally, we introduce an efficient discretisation algorithm for our newly derived integral equations and use this scheme to estimate rates of transmission from serological prevalence of SARS-CoV-2 in the UK and historical incidence data on Influenza, Measles, SARS and Smallpox.

Time-varying Crump-Mode-Jagers outbreak model
Throughout the paper, we shall work with an infectious disease outbreak model based on the Crump-Mode-Jagers (CMJ) branching process, which we extend to allow transmission dynam-ics to vary over time.Our formulation is inspired by Vatutin and Zubkov [56,57], who give an exposition of the corresponding time-invariant CMJ process using random characteristics.In our time-varying CMJ outbreak model, the initial infection occurs at non-random time τ ≥ 0. All subsequent infections are "progeny" of this index case, and we shall denote the set of these infected individuals by I * .We denote the set of all infected individuals (i.e., including the index case) by I.
The index case corresponds to an individual endowed with a collection of random elements indexed by the infection time, where, for any τ ≥ 0, • L τ is a (strictly) positive random variable representing the amount of time the individual remains infected, • χ τ (•) is a stochastic process on [0, ∞) which we shall call the random characteristic of the individual, and keeping track of the new infections, i.e., "offspring", generated by the individual.
For completeness, we set N τ (u) := 0 =: χ τ (u) for u < 0. (We will explain the precise roles of χ τ (•) and N τ (•) shortly.)The objects L τ , χ τ (•), and N τ (•) are typically interdependent, as we shall see below, whilst the interdependence of (L τ , χ τ (•), N τ (•)) and (L τ , χ τ (•), N τ (•)) for different τ and τ is in fact immaterial and will be glossed over.We shall moreover endow each individual i ∈ I * with {L τ i , χ τ i (•), N τ i (•)} τ ≥0 , which is an independent copy of {L τ , χ τ (•), N τ (•)} τ ≥0 .(By an independent copy we mean a new random element which is equal in distribution to the original one and independent of it.)Suppose now that individual i ∈ I is infected at (possibly random) time τ i ≥ τ .Intuitively, the infection time τ i then "selects" L τ i i , χ τ i i (•), and , which the subsequent infection dynamics of this individual will "follow."(Note that the collection ) More concretely, N τ i i (u) now stands for the number of new infections generated by the individual i up to time u + τ i .
Example 1 (Bellman-Harris process).The Bellman-Harris branching model can informally be characterised, in the context of epidemics, by the principle that each individual generates a random number of new infections which occur simultaneously at a random time.Once these new infections have occurred, the individual immediately ceases to be infectious.Let ξ(•) be a stochastic process on [0, ∞) with values in N := {0, 1, . ..}, independent of {L τ } τ ≥0 , and then define This specification gives rise to the time-varying Bellman-Harris branching process studied by Kimmel [40].When the distributions of L τ and ξ(t) do not depend on the time parameters τ ≥ 0 and t ≥ 0, we recover the classical Bellman-Harris process [5].
Example 2 (Inhomogeneous Poisson process model).In contrast to the Bellman-Harris process, we can consider a more realistic epidemiological model where each infected individual generates new infections randomly and one by one according to an inhomogeneous Poisson process until they cease to be infectious.This process, with a constant rate of transmission has been previously studied in the context of the generation time [55].The infinitesimal rate at time t of new infections generated by an individual originally infected at time τ ≤ t is specified as where ρ(•) is a non-negative function that models population-level variation in transmissibility while k(•) is another non-negative function describing how individual-level infectiousness varies over time [55].For example, specifying k(t) to be low or zero for small t can be used to incorporate an incubation period in the model.Let Φ(•) be a unit-rate, homogeneous Poisson process on [0, ∞), independent of {L τ } τ ≥0 .Then we can define this model explicitly by (If ρ(t) ≡ ρ and k(t) ≡ k, both constant, then new infections follow a homogeneous Poisson process with rate ρk until the individual is no longer infected.).
Example 3 (Lévy and Cox process models).In the inhomogeneous Poisson process model of Example 2, tractability does not hinge on the assumption that Φ(•) is a Poisson process.We could in fact replace it with a more general, integer-valued Lévy process (i.e., a process with independent and identically distributed increments), where jumps need not be of unit size (e.g., a compound Poisson process).Similarly, replacing the deterministic function ρ(•) with a stochastic process, as long as it is independent of Φ(•), would be straightforward.In the Poisson case, this would turn N τ (•) into to a doubly-stochastic Cox process.However, for simplicity and concreteness, we shall stick to the simpler setting of Example 2.
Remark 4 (Epidemiological interpretation of L τ and k).In the Bellman-Harris process of Example 1, L τ is directly interpreted as the generation interval [55], that is, the time taken for the secondary cases to be infected by a primary case.In the Bellman-Harris process all infections happen at the same time -for example in Figure 2 we have ξ(τ + L τ ) = 3, after L τ time units has elapsed since the index case was infected at time τ .In contrast, in the inhomogenous Poisson process model of Example 2 (and also the Lévy and Cox process models of Example 3), L τ corresponds to how long an individual remains infected (the duration of infection).During this period, an individual can infect others with rate that depends on ρ(•), which describes the calendar-time variation of overall infectiousness in the population, and on k(•), which in turn describes how the infectiousness of each infected individual varies over the course of their infection.The individual's infectiousness profile k(•) can be set as constant, i.e., variation in the individual's infectiousness is only due to calendar-time variation in overall infectiousness.If k(•) is specified to vary significantly, by contrast, then it is advisable to ensure that infections are most likely to end when infectiousness is low.Concretely, this means that k(•) should then be paired with L τ such that the bulk of its distribution coincides with low values of k(•), as in the empirical application in Section 3.2 below.
The random characteristic is used merely as a book-keeping device, to keep track of an individual's infection status in two ways -whether they have been infected in the past or, alternatively, whether they are infected at the moment.It is fundamental to obtaining a unified derivation of both cumulative incidence and prevalence in what follows.
Example 5 (Cumulative incidence and prevalence).The random characteristic (in fact non-random!) determines whether the individual has been infected by time u + τ and is therefore used to derive cumulative incidence.The random characteristic determines whether the individual remains infected at time u + τ and is used to derive prevalence.

Cumulative incidence and prevalence
We will now derive integral equations for cumulative incidence and prevalence under this model.To this end, we study the stochastic process recalling that τ is the infection time of the index case.For the given random characteristic (6), Z(t, τ ) counts the number of infections occurred by time t and for (7) the number of infected individuals at time t, respectively.Our goal is to derive an equation for the expectation of Z(t, τ ), covering both cases.
Before embarking on the derivation of the equation governing E[Z(t, τ )], we shall first introduce technical assumptions ensuring E[Z(t, τ )] < ∞, which a fortiori guarantees that Z(t, τ ) is finite with probability one, a property known as regularity in the branching process literature [53].
Regarding N τ (•), we write and henceforth assume that there is a non-decreasing, right-continuous function (We will give sufficient conditions that imply this assumption in the context of Examples 1 and 2 below in Examples 13 and 16, respectively.)Moreover, we assume that the random characteristic χ τ (•) satisfies 0 ≤ χ τ (u) ≤ 1 for any τ ≥ 0 and u ≥ 0, which evidently accommodates both ( 6) and ( 7) from Example 5.Under these assumptions, straightforward adaptation of the proof of Lemma 4.2 in [18] (cf. the proof of Theorem 2.1 in [22]) yields E[Z(t, τ )] < ∞ for any t ≥ τ ≥ 0. Now, singling out the index case, we can write The key insight in the analysis of ( 9) is to stratify the infected individuals in I * according to their (unique) "ancestor" among the individuals infected by the index case.More concretely, let i 1 , i 2 , . . .∈ I * label the "offspring" of the index case in chronological order, i.e., so that τ ≤ τ i 1 ≤ τ i 2 ≤ • • • , and let I k ⊂ I * for each k = 1, 2, . . .denote the set consisting of i k and its "progeny."We can then write .
This is an analogue of the principle of first generation for the Bellman-Harris process [35, Theorem 6.1] (see also [40, p. 5]).
Conditional on the random times τ i 1 , τ i 2 , . .., the random variables Z 1 (t), Z 2 (t), . . .can be shown to be mutually independent, with Z k (t) equal in distribution to Z(t, τ i k ), where Z(•, τ ) τ ≥0 is an independent copy of {Z(•, τ )} τ ≥0 (independent of τ i 1 , τ i 2 , . .., in particular).Thus, where, using the law of total expectation, For the random characteristic (6), f (t, τ ) is the cumulative incidence at time t, and we shall denote it by CI(t, τ ).Since E[χ τ (t − τ )] = 1 in this case for t ≥ τ , the equation (10 In the case (7), f (t, τ ) is the prevalence at time t, which we henceforth denote by Pr(t, τ ).In this case, for the survival function associated with G τ (•), we have then We will henceforth assume that the maximal reproduction number R := sup t≥0 R(t) is finite and that the function Intuitively, this condition ensures that the distribution of L τ does not become too concentrated near zero over time.We can then define a non-decreasing, right-continuous function Λ(u) := R G(u), u ≥ 0 which, in view of ( 14) and ( 15), satisfies the assumption (8).Hence, we deduce E[Z(t, τ )] < ∞, and the regularity of the branching process follows.Inserting ( 14) into (11) and (12), respectively, we obtain which agree with [40,Theorem 5.1].When G τ (•) admits a PDF g τ (•), the most relevant case in practice, we can simplify the equations further to Example 16 (Inhomogeneous Poisson process model, cont'd).To analyse the Poisson process model of Example 2, we note first that for any u ≥ 0, If we assume, say, that ρ(•) and k(•) are bounded, then Λ(•) is non-decreasing, continuous and satisfies (8), implying E[Z(t, τ )] < ∞ and the regularity of the branching process.(The CDF G τ (•) does not play a role in regularity for this model since, unlike in the Bellman-Harris process, the random variable L τ cannot precipitate secondary infections.)To work out an expression for Λ τ (u), we shall further assume that G τ (•) admits a PDF g τ (•), as above.Invoking the independence between Φ(•) and {L τ } τ ≥0 and the law of total expectation, we obtain and Integrating (17) by parts, whereby the equations for cumulative incidence and prevalence read as respectively, in this case.
Example 20 (Alternative Poisson process model).An alternative version of the inhomogeneous Poisson process model, suggested by an anonymous referee, can be formulated by assuming that an infected individual's infectiousness evolves at a time scale determined by where k(•) is a continuous function, (strictly) positive on the interval [0, 1] and zero elsewhere.Concretely, we then define where ρ(•) and Φ(•) are as in Examples 2 and 16.Note that, by the properties of k(•), we have Since k(•) is necessarily bounded, the regularity of the resulting branching process is guaranteed.Moreover, we note that cumulative incidence and prevalence for the model can be analysed along the lines of the original Poisson process model simply by substituting k(u)G τ (u) with gτ (u) in the integral equations ( 18) and (19).
Remark 21 (Probability generating functions).In the Bellman-Harris case of Examples 1 and 13, we can also analyse the distribution of Z(t, τ ) via its generating function φ(s; t, τ , letting us study, e.g., higher moments.Concretely, one can show that φ( • ; t, τ ) satisfies the integral equations for random characteristics ( 6) and (7), respectively, where ψ(s; t) . These are special cases of [40,Equations (3.3) and (3.4)], whilst self-contained re-derivations in the case where G τ (•) does not depend on the infection time τ are given in [5].
Remark 22 (Relationship between ρ(t) and R(t)).The quantity R(t) in the context of the Bellman-Harris process (Examples 1 and 13) is more precisely the instantaneous reproduction number, i.e., the expected number of secondary cases arising from a primary case when those infections occur at time t.In the context of a real-time epidemic, R(t) is generally interpreted as the average number of secondary cases that would arise from a primary case infected at time t if conditions remained the same after time t [27].The quantity ρ(t) in the Poisson process model (Examples 2 and 16), in contrast, is a time varying transmission rate, i.e., scaled by time, and therefore exists on a different scale.An alternative way of analysing R(•) is to use the case reproduction number R(t) [33,58], which represents the average number of secondary cases arising from a primary case infected at time t, i.e., transmissibility after time t.It is similarly possible to also analyse ρ(•) through the case reproduction number and therefore compare the rates of transmission in both models commensurably.Namely, given ρ(•) and R(•), they can be transformed into R(•) and be comparable on the same scale via Bellman-Harris process, the precise interpretation of g τ (•) is the PDF of the time between an individual becoming infected and occurrence of all subsequent infections generated by the individual, i.e., the generation time or interval [55].In contrast, the Poisson process model is parameterised by the product of the infectiousness profile k(•), which broadly corresponds to the generation time [15], and the survival function G τ (•) of the duration of the infection.
Generally, these two models differ in terms of their behaviour.That said, they give rise to equivalent cumulative incidence and prevalence provided Hence, cumulative incidence and prevalence roughly agree between the two models when the infectiousness profile k(•) approximates the hazard function of L τ , i.e., the right-hand side of (24).Even in this case, the higher moments of the models typically do not agree, however.

Incidence
Incidence is defined as the time-derivative of cumulative incidence.To derive an integral equation for incidence à la (11) and ( 12), we shall assume that the function Λ τ (•) is continuously differentiable, that is, for some continuous function λ τ (•).The function λ τ (•) is necessarily non-negative since N τ (•) is a counting process.The assumption (25) rules out infections occurring in a discrete time grid.It is satisfied with Cumulative incidence, by definition, equals zero before the index case is infected at time τ , whilst it then jumps to one.Hence, cumulative incidence, when understood as a function on the entire real line, satisfies (When t < τ we will interpret the integral, and similar integrals in what follows, as zero.)Incidence is then defined as the time-derivative Before deriving incidence in full generality, let us however study the time-derivative of a related quantity CI(t, τ which omits the initial jump and, in view of ( 26), satisfies Applying the Leibniz integral rule to the second integral on the right-hand side of ( 27) formally (see Remark 33 below), we obtain We have taken (0, t − τ ) as the integration domain since ∂ ∂t CI(t, u + τ ) and I(t, u + τ ) do not agree at u = t − τ for reasons that will become clear in the next paragraph.
Whilst (29) already describes incidence for t > τ , for further developments in Sections 2.4 and 2.5 it is essential that we have an equation characterising incidence for any t ≥ τ .Thus, we need to also deal with the case t = τ where the time-derivative ∂ ∂t CI(t, τ ) cannot be defined in the classical sense due to the jump in cumulative incidence.To this end, it is helpful to note that the derivative of t → 1 [0,∞) (t − τ ) may be understood as a Dirac delta function δ( • − τ ) in a distributional sense.We recall that the Dirac delta function is a generalised function with the characteristic property R f (x)δ(y − x)dx = f (y).Now, In particular, formally I(τ, τ ) = δ(0) + λ τ (0).(30) Note that u)du = 0. Thus, we can write the right-hand side of ( 29) as a single integral over (0, t − τ ], i.e., I(t, τ ) = (0,t−τ ] I(t, u + τ )λ τ (u)du, t > τ.
Consequently, we find that incidence is generally governed by the equation Remark 32.In (31), we have adjusted the integration domain from (0, t − τ ] to [0, t − τ ] to ensure that the equation agrees with (30) for t = τ .(This adjustment is immaterial for t > τ .)To see why this is the case, note that the right-hand side of (31) consists of the generalised function δ(t − τ ) and the integral [0,t−τ ] I(t, u + τ )λ τ (u)du, the latter of which is an ordinary function in t regardless of what the nature of I(t, u + τ ) is.Once we integrate I(t, u + τ )λ τ (u) with respect to u over the singleton {0} in the case t = τ , integration will only pick up the generalised function part of I(τ, u + τ ), i.e., δ τ − (u + τ ) = δ(u), producing the term λ τ (0), as intended.
Remark 33.When applying the Leibniz integral rule in (28), we have not attempted to verify its assumptions.In fact, doing so would be difficult since we do not know a priori that cumulative incidence is differentiable with respect to time.Proving its differentiability from first principles using Lebesgue's dominated convergence theorem would similarly be difficult since it is not straightforward to derive sufficiently sharp a priori estimates for the increments of t → CI(t, τ ).However, there is an alternative way of proving ( 29) and (31) rigorously, which can be outlined as follows.We first treat these equations as an educated guess and show they have a (unique) solution.We can then show that the time-integral of the solution satisfies the equation ( 11) for cumulative incidence.Finally, it is straightforward to prove uniqueness of solutions for (11) using Grönwall's lemma (cf.Appendix B), which then lets us conclude that the timederivative of cumulative incidence indeed follows (31).We will elaborate on the remaining mathematical details of this argument, including rigorous treatment of the Dirac delta function as a generalised function, in a separate paper.

Example 34 (Incidence for the Bellman-Harris process and Poisson process model).
Under the aforementioned assumptions, equations ( 29) and ( 31) read as respectively, for the Bellman-Harris process of Examples 1 and 13, and as respectively, for the Poisson process model of Examples 2 and 16.

Consistency with back-calculation
Back-calculation is a standard method to recover prevalence from incidence by convolving the survival function of the generation interval with incidence [10,20].We will now show that the equations we have obtained for prevalence and incidence are consistent with the back-calculation relationship under the assumption (25) and the additional assumption that the CDF G τ (•) does not depend on the infection time τ , in which case we write G(•) and G(•) in lieu of G τ (•) and G τ (•), respectively.
Let f and f be two functions, one of which may be a generalised function, such that f (t) = 0 for any t < 0 and f (t) = 0 for any t < τ .Their convolution can be expressed as for any t ≥ τ and equals zero otherwise.We proceed now to show that the back-calculation relationship holds, with the convention G(t) := 0 for any t < 0. Starting from (31), we have where the first term on the right-hand side can be computed as The second term on the right-hand side of (37) vanishes for any argument t ≤ τ , so it suffices to consider t > τ .In this case, switching the order of integration, we obtain Combining ( 38) and ( 39), we have altogether Matching this with equation ( 12) under the assumption (25), we deduce By an application of Grönwall's inequality, as outlined in Appendix B, we can finally conclude that the back-calculation relationship (36) holds.
Remark 40 (Modelling HIV incidence from prevalence).HIV is an example of a disease where, due to long incubation times, routine surveillance generally returns prevalence -not incidence [21].However, what is of interest to policy makers is incidence, not prevalence [11].
Common approaches all make use of the back-calculation relationship through convolving a latent function for incidence with the survival function G(•) [11,44,52].Our argument above shows that there is no need to model incidence as a latent function, rather one can fit ρ(•) or R(•) directly to prevalence data using the prevalence integral equation for Pr(t, τ ), after which I(t, τ ) can be computed directly without need for a latent incidence function.This relationship therefore can help facilitate simpler or more pragmatic modelling choices.

Consistency with a common renewal equation model for incidence
The key difference between our newly derived integral equations and the common renewal equation used [17,27,45] is the inclusion of the parameter τ that initially arises due to the timing of the index case.The inclusion of τ means that we need to work with I(t, τ ), not simply I(t), and also gives rise to terms outside of the integral depending on whether one is interested in incidence, cumulative incidence or prevalence.
As in Section 2.4, we assume that G τ (•) does not depend on τ , i.e., we work with G(•), and we moreover assume that G(•) has a PDF g(•).In this context, when extended to accommodate the general initial infection time τ , the common renewal equation for incidence is tantamount to the integral equation We show that the renewal equation ( 41) in fact agrees with the integral equation ( 31) in the Bellman-Harris case, that is, While we focus on the Bellman-Harris process (Examples 1 and 13) here for notational simplicity, the argument also applies to the Poisson process model (Examples 2 and 16) simply by replacing R(•) with ρ(•) and g(•) with k(•)G(•), respectively, throughout.
To this end, we first introduce so that, given (41), Applying (43) to the integrand in (42) yields Additionally, we introduce Subsequently, by applying (41) to the integrand in (45) and switching the order of integration, we obtain The integral equations ( 44) and ( 46) then imply the bound and applying Grönwall's inequality as outlined in Appendix B we deduce that Finally, by ( 43) and ( 47), Given (31) in the Bellman-Harris case, we then have the bound for any t ≥ τ ≥ 0. Applying the result in Appendix B again we conclude that, indeed, I(t, τ ) = I Ren (t, τ ) holds for any t ≥ τ ≥ 0.
Remark 48 (Equivalence does not extend beyond incidence).In the case of prevalence or cumulative incidence, the equivalence between the common renewal equation and our newly derived integral equations is broken.This is easy to see by examining the derivations leading to ( 44) and ( 46).If we considered cumulative incidence, for example, a constant one instead of a Dirac delta function would appear and the leading terms in ( 44) and ( 46) would no longer agree, rendering the rest of the argument impossible to carry through.This illustrates why the common renewal equation is a special case of our integral equations only when the index case is infected at time τ = 0 and when considering incidence.Simpler renewal equations that do not involve varying dependence on the parameter τ for prevalence or cumulative incidence are not possible.
3 Numerical implementation and empirical application

Discretisation of integral equations
The integral equations for cumulative incidence, prevalence and incidence under the assumption ( 25) are all special cases of a generic equation with the choices Recall that for the Bellman-Harris process of Examples 1 and 13, we substitute λ τ (u)du := R(u + τ )g τ (u)du and for the Poisson process model of Examples 2 and 16, λ τ (u)du := ρ(u + τ )k(u)G τ (u)du.
A key hurdle in solving the equation ( 49) is that on the right-hand side, we get f (t, u) for τ ≤ u ≤ t and not f (u, τ ) for τ ≤ u ≤ t.What this means is that in order to solve f (t, 0) for t ≥ 0, say, we need to actually solve f (t, τ ) for any pair (t, τ ) such that t ≥ τ ≥ 0. This is in fact why we left the initial infection time τ as a free parameter.(Alternatively, we could view (49) as a system of coupled integral equations, indexed by τ , that need to be solved simultaneously.) Solving the equation ( 49) numerically is greatly facilitated if we introduce the auxiliary quantity f c (t) := f (c, c − t) for any c ≥ t ≥ 0. From ( 49) we can deduce that, for fixed c ≥ 0, the singleargument function f c (•) is governed by the renewal-like integral equation We then recover f (t, 0) for t ≥ 0 via f (t, 0) = f t (t).In practice, we are interested in solving f (t, 0) discretely for t = 0, ∆, . . ., N ∆ for some N ∈ N and ∆ > 0. To this end, we approximate f n∆ (•) recursively by for i = 0, . . ., n do 3: else 6: end if end for 9: end for 10: return f n∆ (n∆) (≈ f (n∆, 0)) for any n = 0, . . ., N Algorithm 2 Discretisation of integral equations, vectorised for any n = 0, . . ., N .For clarity, we present the entire procedure in pseudo-code in Algorithm 1.
Example 51.Concretely, in the Bellman-Harris case, we set while in the case of the Poisson process model, A simplified version of the algorithm for cumulative incidence in the Bellman-Harris case is found in Appendix A.
In practice, the double for-loop in Algorithm 1 may lead to computational inefficiency when N is large and an interpreted language is used, so it is useful to refine it by vectorisation.To this end, for an m × n matrix A and 1 ≤ s ≤ m and 1 ≤ t ≤ n, we denote by A[s, t] the s-th row, t-th column element of A. Moreover, for 1 ≤ i ≤ j ≤ m and 1 ≤ k ≤ l ≤ n, we write A[i : j, k : l] for the sub-matrix consisting of each element A[s, t] where i ≤ s ≤ j and k ≤ t ≤ l.
(If i = j, we simply write i in lieu of i : j.)We also denote by element-wise (Hadamard) multiplication of matrices.The vectorised version of Algorithm 1 is given as Algorithm 2. This matrix computation is possible by observing that all relevant values of the functions h(•, •) and λ • (•) can be stored in the matrices H and L, respectively.Algorithm 2 can be further vectorised with respect to parameters to produce simultaneously discretisations for multiple parameter values.Additional computational savings could be attained in Algorithm 2 by observing that the topleft corner of the matrix L typically contains very small values since g τ (u) and G τ (u) are small with large u.Therefore the matrices L and F could in practice be truncated with a small error in the computation of diag(F ).
We illustrate the use of these algorithms in Figure 1, where we compute prevalence using Algorithm 2 and compare the results with statistical estimates of prevalence from a Monte Carlo simulation.Python implementations of Algorithms 1 and 2, including a version of the latter vectorised over parameters, are provided as fully documented Jupyter notebooks in: https://github.com/mspakkanen/integral-equations

Bayesian inference on empirical data
We perform Bayesian inference to estimate the time-varying case reproduction number R(t), as defined in Remark 22, for historical incidence data for Influenza [29], Measles [34], SARS [42] and Smallpox [30] and for recent SARS-CoV-2 serological prevalence data in the United Kingdom [48].
∼ Normal(0, σ)   [29], Measles [34], SARS [42] and Smallpox [30].Plots show the case reproduction number R(t), the distribution g(•) in discretised form and incidence for the process.Solid black lines in all plots are means, and the two red envelopes are the interquartile and 95% credible intervals.The horizontal blue line indicates R = 1.The x-axis in all plots is time measured in days.

Historical incidence data
Historical incidence data for Influenza [29], Measles [34], SARS [42] and Smallpox [30] have been extensively used in validating renewal equation frameworks [17].We fit an integral equation for the Bellman-Harris process.We work with G τ (•) = G(•) that does not depend on τ .As demonstrated in Section 2.5, the corresponding integral equation agrees with the common renewal equation ubiquitously used in the modelling of incidence [17].
We first introduce a probabilistic model for the function R(•) through a stochastic random walk process.To aid comparability to alternative methods [58], we transform R(•) to the case reproduction number R(t), which represents the average number of secondary cases arising from a primary case infected at time t, i.e., transmissibility after time t.In Table 1, the Negative Binomial likelihood is re-parameterised to the mean-variance formulation, y is the observed count data (number of infections), φ is the overdispersion parameter and σ is the random walk variance parameter.Therein, we write Normal + (0, a) for a normal distribution Normal(0, a) constrained to the positive real axis.The observed count data and generation intervals were obtained from [16,17].The priors were selected to be weakly informative and were generally robust to change.
Algorithm 2 was used to discretise and solve t → I(t, 0) -recall that τ is a parameter that is intrinsically involved in the solution of the integral equation, although we can ultimately restrict our attention to t → I(t, 0) only, having assumed that the first infection occurs at time τ = 0.For all data sets, an arbitrary seeding period of 10 days was used to correct for poor surveillance in the early epidemic.The seeding period was not included in the likelihood and we found our fits to be robust to different choices of seeding duration.Posterior sampling was using Hamiltonian Monte Carlo (1000 warmup/1000 sampling with multiple chains) in the Bayesian probabilistic programming language Numpyro [9,46].Posterior predictive checks were performed by examining R-hat and K-hat distributions. Figure 3 shows the estimated case reproduction numbers R(t), which, as expected, match those previously estimated [17].

Serological prevalence data
The ONS infection survey, is a weekly, household cross-sectional survey of blood samples which are used to test for the presence of COVID-19 antibodies, led by the Office for National Statistics (ONS) and the Department of Health and Social Care of the United Kingdom.At any point in time the ONS infection survey provides an estimate for the number of individuals currently infected with SARS-CoV-2, i.e., the prevalence of infection/positivity rates.Estimation of incidence from the ONS infection survey is done using a bespoke deconvolution approach, and estimating R(t) or incidence directly from prevalence, to our knowledge, has not been attempted.
We study estimates of prevalence from the ONS infection survey over the period 5th April 2021 to 15th November 2021.Our choice for this period arose from the requirement of wide spread, easily accessible SARS-CoV-2 PCR testing in the general population, which is required to ensure comparability between the ONS infection survey and reported case data (which we compare our estimates to).Prevalence estimates are reported weekly and we therefore use smoothing splines to interpolate these weekly estimates to daily estimates through a log linear generalised additive spline model [36].ONS infection survey results are generally reported on the Friday of any given week, with the sampling period covering Wednesday to Wednesday -a period of 10 days.We therefore incorporate this observation lag by convolving daily prevalence with Normal(10, 0.3) distribution to adjust for these reporting lags and incorporate some uncertainty in this lag.We fit a Poisson process model, detailed in Table 2, to these lagged prevalence data assuming the infectiousness profile to be analogous to the generation time such that k(•) is given by the PDF of Gamma(4.84,1.73) distribution [54].The CDF G(•) of the infection duration was assumed to follow the CDF of Normal(10, 1.5) distribution [59].
We did fit an aggregated likelihood where a daily Poisson process was aggregated to weekly averages, but found little difference in results.
The top left panel of Figure 4 shows the estimate case reproduction number R(t), which fluctuates around 1 over the period of study.The top right panel exhibits an excellent posterior fit to the daily smoothed ONS infection survey prevalence.Moreover, the bottom right panel shows infection incidence using the estimated R(t) from fitting prevalence, and bars indicate the reported case data.Note we do not fit directly to the case data, but only to the prevalence as estimated by the ONS survey -however including a second likelihood would be trivial to add.Lagging the time series and estimating the maximum cross correlation suggest a lag of approximately 7 days between infections and reported cases -a lag that is in line with previous studies [54].Finally, correcting for this lag between infections and cases, we see a reasonably stable (aside from weekly reporting cycles) infection ascertainment ratio (bottom right panel of Figure 4) with a mean of approximately 2.5, implying that for most of the study period there were 2.5 times more infections than reported and that this is relatively stable given testing policies over the period.This example demonstrates how our framework can fit prevalence directly without the need of deconvolution type approaches.

Discussion
Our primary goal in this paper is to bridge the worlds of individual-based models and mechanistic models to gain from the best of both.To this end, we began by choosing the most general branching process available -the Crump-Mode-Jagers process [18,19,38].In the Crump-Mode-Jagers process, an epidemic is created at an individual level where, from a single infected individual, subsequent infections occur at random times according to their level of infectiousness.To our knowledge, for the first time, we generalise the Crump-Mode-Jagers process to allow for fully time-varying reproduction process for new infections.Indeed, rather than assuming the distribution of new infections to be constant (corresponding to a basic reproduction number) we allow it to change over time, which is essential in the modelling of real outbreaks [33] beyond their early phase.We find that under this generalisation, a general integral equation arises from the Crump-Mode-Jagers process.Our framework also allows us to specify the dynamics of how new infections arise (in addition to them changing over time).Studying first the case where each infected individual produces all of their secondary cases, or "offspring", at the same random time, we recover the well known Bellman-Harris process [5].Studying a more complex assumption where each infection can give rise to its offspring over the duration of its infection (an inhomogenous Poisson process) we derive a new integral equation, which to our knowledge, has not been previously presented.Remarkably, we find that despite the Poisson process model being much more complex than the simple Bellman-Harris assumption, the resultant integral equation has exactly the same form as the Bellman-Harris integral equation, only instead of the generation interval CDF, the survival probability is used.
Through starting from a stochastic process, we are able to define prevalence, incidence and cumulative incidence as summary statistics (via moments) of an individual-based infection process.The benefit of defining these well known epidemiological quantities from a single stochastic process is that they are, by design, consistent with one another -i.e., they are parameterised with the same generation interval and transmission rate (either ρ(t) or R(t)).This allows practitioners to fit to prevalence for example, and easily recover incidence with no additional fitting.We mathematically show that this is the case and prove our equations for prevalence and incidence are consistent under the commonly used back-calculation technique in epidemiology [10].Given ever increasing amount of infectious disease surveillance, being able to model prevalence and incidence simultaneously under the same process can greatly improve estimates of the rates of the reproduction number.A recent example is the COVID-19 pandemic where several countries collected high quality data on both cases (incidence) and serology (prevalence) [26].
We also show that the incidence integral equations we recover from the Bellman-Harris process and from the Poisson process model are in fact in agreement with the renewal equation commonly used in the modelling of incidence [17].Specifically, the common renewal equation is a special case of our incidence equations under the scenario where the first infection occurs at a specific, non-random, time.We also show that our equations are more general, and accommodate the modelling of prevalence, cumulative incidence, complex importation functions, and time-varying generation times [40].The common renewal equation is computationally simpler as it does not involve the time τ of the first infection and simplifies the problem from two-dimensional to one-dimensional.We have however introduced an efficient algorithm which relies on straightforward matrix algebra to compute our more general integral equations.Given the ability of modern computers to perform matrix operations efficiently, we do not believe the computational overhead of our integral equations is meaningfully greater than that of the simple renewal equation.However, our integral equations allow for a far greater range of modelling choices with explicitly stated assumptions.
In this work, we have attempted to put the modelling of infectious diseases using renewal equations on firm mathematical ground.These mathematical foundations are broad enough to cover a variety of model specifications for transmission dynamics, and from them we can extract information about a wide range of relevant epidemiological quantities.In doing so, we have once again made explicit the connection between branching processes [5,18] and renewal equations.Explicit links between renewal equations and SEIR models [14] and Hawkes processes [49] have been previously noted.It is likely other such relationships exist, and this is an interesting area of further study.Of additional interest is to use our framework to study the more complex Lévy and Cox process models, which may produce renewal equations with even more realistic dynamics.Equally, recent frameworks [32,51] have extended the seminal work of [58] to estimate case reproduction number on graphs -connecting these two approaches is an interesting area of future research.Finally, our framework, and the vast majority of previous frameworks, only consider the mean integral equation and ignore the dynamics of higher-order moments.Using our framework, we can recover these moments from our stochastic process and formulate more accurate likelihoods for model fitting.

Figure 2 :
Figure 2: Schematic of infections generated under a Bellman-Harris process and an inhomogeneous Poisson process model.In a Bellman-Harris process, after a generation interval has elapsed, new infections happen at the same time (instantaneously).In the inhomogeneous Poisson process model, an individual is infectious for a period, over which their infectiousness varies, and they produce infections one by one.

Remark 23 (
When do the Bellman-Harris process and the Poisson process model agree?).The fundamental difference between the Bellman-Harris (Example 13) and the Poisson process model (Example 16) integral equations is that the Bellman-Harris integral equations are parameterised by g τ (•), and the Poisson process model equations by k(•)G τ (•).Within the u) in Example 16 provided ρ(•) and k(•) are continuous, and with λ τ (u) = R(u + τ )g τ (u) in Example 13 provided R(•) is continuous and G τ (•) has a continuous PDF g τ (•).

Figure 3 :
Figure3: Bayesian modelling of incidence for Influenza[29], Measles[34], SARS[42] and Smallpox[30].Plots show the case reproduction number R(t), the distribution g(•) in discretised form and incidence for the process.Solid black lines in all plots are means, and the two red envelopes are the interquartile and 95% credible intervals.The horizontal blue line indicates R = 1.The x-axis in all plots is time measured in days.
Example 13 (Bellman-Harris process, cont'd).Let us consider the Bellman-Harris case of Example 1 and write R(t) := E[ξ(t)] for the (time-varying) reproduction number at time t ≥ 0. Let us also denote the indicator function of a set A by 1 A .Using the law of total expectation and the independence between ξ(•) and {L τ } τ ≥0 , we get then

Table 1 :
Hierarchical Bayesian model for estimating incidence for a Bellman-Harris process

Table 2 :
Hierarchical Bayesian model for estimating prevalence for a Poisson process model Bayesian modelling of the ONS COVID-19 infection survey for prevalence.Top left show the case reproduction number R(t), top right prevalence, bottom left incidence and bottom right the ascertainment ratio (incidence/reported cases).Solid black lines in all plots are means, and the two red envelopes are the interquartile and 95% credible intervals.The horizontal blue line indicates R = 1.The x-axis in all plots is decimal calendar time.The ascertainment ratio in the bottom right is adjusted for the reporting delay between infections and cases, and this delay is estimated as the maximal lagged cross correlation.