From individual-based epidemic models to McKendrick-von Foerster PDEs: A guide to modeling and inferring COVID-19 dynamics

The transmission of pathogens between species is a global concern [38, 11]. As such zoonotic episodes are expected to become increasingly common in humans, it is critical to develop analytic tools that can quickly transform epidemiological observations into informed public policy in order to mitigate and control epidemics.

A novel coronavirus, SARS-CoV-2, has recently crossed the species barrier into humans and, within months, has rapidly spread to all corners of our planet [60]. The sheer scale of this pandemic has overburdened our medical infrastructure, caused fatalities estimated well into millions, and shut down entire economies. Remarkably, the rapid spread of COVID-19 and its consequences can be attributed to the unique life cycle of a 30,000 base pair single-stranded virus. SARS-CoV-2 is an airborne pathogen transmitted by both symptomatic and asymptomatic carriers in close proximity to non-infected individuals. Milder COVID-19 symptoms include a dry cough, fever, and/or shortness of breath while more serious cases include respiratory failure and possible death. With millions of infections and hundreds of thousands of documented deaths and recoveries, the COVID-19 pandemic is providing a wealth of independent estimates of important clinical characteristics that can help predict health outcomes specific for a country or region.

It quickly became understood that accurate descriptions of the life cycle of this disease needed to distinguish between several stages of the disease, referred to as compartments, depending on whether an infected individual is infectious or not, symptomatic or not, hospitalized, etc. However it remains unclear to what extent making precise predictions of the dynamics of such a complex disease requires to have a precise knowledge of clinical features such as incubation period, generation time, and duration times between infection, symptom establishment, hospitalization, recovery and death, to know how these durations correlate and what are the exact probabilities of transition between stages.

In this work, we consider a fully stochastic, generic epidemiological model with an arbitrary number of compartments, that encompasses life cycles of most complex diseases and that of COVID-19 in particular. We show how structuring the infected population by its infection age, i.e., time elapsed since infection, allows us to decouple dependencies between stages and to time. More specifically, when the population size is large enough, the joint evolution of all compartment sizes can be described by means of a linear, first-order partial differential equation (PDE) known as the McKendrick-von Foerster equation describing the number $n(t,a)$ of infecteds of (infection) age $a$ at time $t$. The boundary condition at age 0 is driven by the infection rate from infecteds of age $a$, averaged over all possible courses of infection, and the number of individuals of age $a$ in compartment $i$ at time $t$ is obtained by thinning $n(t,a)$ by a factor $p(a,i)$ which is the probability of being in compartment $i$ conditional on having age $a$, averaged over all possible courses of infection.

In the case of COVID-19, we display a simple procedure to infer these parameters, some from the biological literature and most from time series of numbers of severe cases, hospitalized cases, discharged patients and deaths that can be applied easily to any regional or national dataset. We also allow for time inhomogeneity in the infection rate to account for temporary mitigation measures such as lockdowns or social distancing. We apply this procedure to French COVID-19 data from March to May 2020 and estimate various parameters of interest including the reproduction number in different phases of the epidemic (before, during, and after lockdown) and biological parameter values that we compare to empirical estimates.

The large population size limit of our stochastic model is a PDE “decorated” with compartments. This point of view extends the usual sets of ODEs used in epidemiology, and allows us to represent in the same framework a large class of deterministic epidemic models. Before describing the stochastic model underlying such decorated PDEs, let us illustrate this notion by recalling the well-known derivation of the classic SIR set of ODEs from an age-structured model.

Consider the solution $(n(t,a);\,t,a\geq 0)$ to the following partial differential equation:

where $0\leq x_{0}\leq 1$, and $g,\tau\geq 0$ fulfill

Equation (1) was first proposed to describe the dynamics of an epidemic where infected individuals are structured by their age of infection, and is known as the Kermack-McKendrick model [36]. Note that in the original work of [36] the model is formulated as the solution to a convolution equation rather than as a PDE, but that the two formulations are equivalent, see Section 2.2. In this context, the age of an individual refers to the time elapsed since its infection, and not to its actual age. Then, $n(t,a)$ is the density at time $t$ of all individuals with age (of infection) $a$, and $S(t)$ the density of individuals still susceptible to the disease. The differential term describes the aging process: the age of an individual increases linearly with time at rate one. The interpretation of the age boundary condition of (1) is that individuals with age $a$ infect susceptible individuals at a rate $\tau(a)$ that only depends on their age.

In order to recover the SIR model, suppose that $\tau$ is given by

for some $\gamma,\beta>0$. Further define

Then a simple calculation shows that $(S,I,R)$ solves the following well-known system of ODEs:

The previous expressions have an interesting probabilistic interpretation. Consider a Markov process $(X(a);\,a\geq 0)$ with two states $I$ and $R$. Suppose that it starts from $I$ and jumps to $R$ at rate $\gamma$. The process $(X(a);\,a\geq 0)$ can be interpreted as describing the sequence of states (infective then recovered) visited by a typical individual in the microscopic model underlying (2). Then, clearly

so that

Furthermore, suppose that a typical infected individual yields new infections at constant rate $\beta$ while it is in state $I$. Then, the mean number of new infections occurring in the time interval $[a,a+\mathrm{d}a]$ is

The picture that emerges from this simple calculation is that, instead of keeping track of the number of individuals in each compartment, one can consider the age structure of the population, given by Eq. (1). The dynamics of the age structure is uniquely prescribed by the average number $\tau(a)$ of infections that an individual yields at age $a$. The individual counts in each compartment can then be recovered by integrating against the age structure the one-dimensional marginals of a process $(X(a);\,a\geq 0)$ that describes the sequence of compartments visited by a typical individual in the population. We say that the PDE is “decorated” with compartments, as the process $(X(a);\,a\geq 0)$ is used to recover the counts in the compartments, and only influences the dynamics of the infection which is described by the sole Eq. (1) through the average infection rate $\tau(a)$.

This viewpoint has several advantages compared to the usual ODE setting of (2). First, we can make sense of (3) for any process $(X(a);\,a\geq 0)$. If this process is not Markovian, the number of individuals in each compartment no longer solves a system of ODEs similar to (2). Hence, this approach allows to go beyond the usual ODE framework. This generalization is of great modeling interest since a hypothesis underlying sets of ODEs is that the sojourn times in each compartment and the time between successive infections are exponentially distributed. In particular they cannot account for sojourn time distributions that are peaked around a value, which have been reported for instance for COVID-19 [41, 58].

Second, regardless of the number of compartments, the age structure of the population is described by the same one-dimensional PDE. This is particularly valuable when models have a large number of compartments, as in the context of COVID-19 [50, 19, 12, 17], as it avoids the use of high-dimensional systems of ODEs that are cumbersome to study mathematically. However, this requires to work with the PDE (1) rather than with ODEs, resulting in a mild computational cost.

Third, Eq. (1) only involves the mean number of infections $\tau(a)$ induced by individuals at age $a$ averaged over all compartments. In particular, it is unnecessary to assess to which compartments individuals belong when they yield new infections. This is in contrast with the usual ODE framework, where for each compartment an infection rate needs to be prescribed. As we will see, $\tau(a)$ relates to well-known epidemiological quantities that can be assessed directly from the data.

The main contribution of our work is to show that such decorated age-structured models arise naturally as the law of large numbers limit of a wide class of general stochastic epidemic models that we now introduce.

We consider a population model in which individuals are either susceptible, if they have not yet met the disease, or infected. Our definition of infected is broader than usual: an individual is infected if it has been infected in the past. In particular, individuals that have recovered or died from the disease are still infected, even if they are not infective anymore. At any point in time, an infected individual is in one of several states, that will also be referred to as compartments, types, classes, or stages. The set of all such states is denoted by $\mathcal{S}$ and is assumed to be finite. Depending on the disease complexity, the number of stages can vary. In the SARS-CoV-2 example, typical stages are asymptomatic, mild case, severe case, hospitalized, intensive care unit, recovered, and dead (see Figure 5).

We assume that upon infection a susceptible individual immediately changes state and never becomes susceptible anymore (ruling out multiple infections, in particular), and that it will eventually end in one of two states: recovered or dead. The sequence of states visited by an individual $x$ is then encoded by a stochastic process $X_{x}\coloneqq(X_{x}(a);\,a\geq 0)$ valued in $\mathcal{S}$, where the random variable $X_{x}(a)$ is the state of $x$ at age of infection $a$. We call $(X_{x}(a);\,a\geq 0)$ the life-cycle process.

Each individual is further endowed with a random point process $\mathcal{P}_{x}$ on $[0,\infty)$, called the infection point process. Each atom of $\mathcal{P}_{x}$ gives the age at which $x$ makes an infectious contact with another individual in the population, and we assume that all atoms are distinct so that secondary infections cannot occur simultaneously. (By an infectious contact, we mean a contact that would lead to a new infection if the target individual is susceptible to the disease. The pair $(\mathcal{P}_{x},X_{x})$ characterizes the course of the infection of individual $x$. We assume that these pairs, for all individuals $x$, are i.i.d. copies of the same pair $(\mathcal{P},X)$ that describes the infection of a typical individual in the population.

In order to make the mathematical treatment of the model easier, we make the simplifying assumption that the number of susceptible individuals is large compared to the number of infected individuals, so that each infectious contact leads to a new infection. (We neglect the saturation in the number of infections due to the finiteness of the population.) The epidemic is then described by a branching process known as a Crump-Mode-Jagers process [31, 54]: an individual $x$ infected at time $\sigma_{x}$ will produce $N_{x}$ new secondary infections at times $\sigma_{x}+A_{1},\dots,\sigma_{x}+A_{N_{x}}$, where $(A_{1},\dots,A_{N_{x}})$ are the atoms of $\mathcal{P}_{x}$, that is ,

(This branching hypothesis is relaxed in a recent work by some of the authors [18].)

Lastly, we superimpose time heterogeneity to this process by means of a contact rate $(c(t);\,t\geq 0)$ valued in $[0,1]$ thinning the infection process. More precisely, if $t$ is a potential time of infection for individual $x$, we ignore the infection with probability $1-c(t)$. This contact rate can model the effect of vaccination, or density-dependence (i.e., relaxing the branching assumption due to an excess of removed or of deceased individuals), or of governmental mitigation measures (i.e., social distancing, lockdown).

The infection process is more formally constructed in Section 2.1.

The stochastic epidemic models we consider here are fairly general and can exhibit quite complex dependencies (i) between states and time, due to the lack of any Markov-type assumption, (ii) between states, due to possibly hidden structuring variables impacting the life cycle, (iii) between state and infection rate, and (iv) between past and future infection events. The main result of this work is that despite this apparent complexity, most of this complexity vanishes when the size of the population is large. More specifically, we show that in the limit of large populations (obtained by starting from a large initial population or as a consequence of natural exponential growth), the population of infected individuals structured by age (of the infection) can be described by means of a one-dimensional PDE, and that the counts in each compartment are recovered by decorating this PDE with the life-cycle process. The limiting expression only depends on:

1. 1.

   The average infection rate

   $$\tau(\mathrm{d}a)\coloneqq\mathbb{E}(\mathcal{P}(\mathrm{d}a)),$$

   formally defined as the intensity measure of $\mathcal{P}$. We make the simplifying assumption that $\tau$ has a density w.r.t. the Lebesgue measure and, with a slight abuse of notation, we still denote it by $\tau(a)$.

2. 2.

   The one-dimensional marginals of the life-cycle process

   $$p(a,i)\coloneqq\mathbb{P}(X(a)=i).$$

We prove two main theorems that are two laws of large numbers for the age and compartment structure in the population: one started from a large number of individuals, and the other from a single individual.

We have illustrated the practical interest of our approach by carrying out parameter inference on data from the early French COVID-19 epidemic. We focus on two important inference aspects of this epidemic: providing estimates for key epidemiological quantities, such as the reproduction number that allows to assess the impact of control measures; and predicting the number of individuals in ICU and hospital to monitor the pressure on the healthcare system.

In our framework, the first task only requires a simple and parsimonious model that can be adjusted on incidence data, whereas fitting the number of individuals in ICU requires a more complex model that better accounts for the population heterogeneity.

Section 1.2 shows that the SIR model can be seen as a decorated age-structured model, with well-chosen Markov life-cycle process and Poisson infection point process. This section extends this representation to a broader class of ODE models. We will be interested in solutions to the following set of ODEs:

The parameter $\beta_{ij}\geq 0$ gives the rate of new infections from individuals in compartment $i$ such that the newly infected individual starts in compartment $j$. The matrix $T$ with entries $(\beta_{ij})$ is referred to as the transmission matrix. For $i\neq j$, $q_{ij}\geq 0$ corresponds to transition rate from compartment $i$ to compartment $j$. The transition matrix with entries $(q_{ij})$ is denoted by $Q$, and we further impose that

This class of ODE models encompasses many common epidemic models, including the SIR and SEIR models, as well as all models described in [8, Chapter 4] for instance.

The previous result provides a simple criterion for a system of ODEs to be represented as a decorated age-structured PDE. This criterion has already been proposed previously and is referred to as the “separable mixing” assumption [15, 14]. A direct consequence of this result is that not all ODE models can be represented using a single decorated age-structured PDE. However many models fulfill the requirement that $\operatorname{rank}(T)=1$, including models with a single infectious state, those where new infected individuals always start in the same state, and all classical models exposed in [8, Chapter 4]. In that sense, our framework greatly extends the usual systems of ODEs widely used in epidemic modeling.

An important situation where $\operatorname{rank}(T)>1$ is that of heterogeneous contact rates in the population. For instance, the contact rate could depend heavily on the (actual) age group to which individuals belong, and more contacts are made within the same age class than between classes. A second example is that of spatial heterogeneity, where contacts are more likely to occur between spatially close individuals. It remains possible to derive a representation of (7) in the case $\operatorname{rank}(T)>1$ using an age-structured model, but this would require to use a multi-dimensional version of (1), which is a straightforward extension of this model.

Deterministic epidemic models where the infectivity depends on the individuals’ age of infection were first introduced in [36], and their mathematical properties have been further studied thoroughly [48, 13, 55], see [30] for a recent account. However, these models have received surprisingly little attention in applications compared to their ODE counterpart, which have been widely used for instance in the context of COVID-19 [49, 50, 19, 12, 17]. In this direction, let us mention [24] which makes use of an epidemic model with memory and [28] were a set of transport PDEs similar to (1) is used. Note that, in contrast with our approach, the PDE in the latter work has two dimensions and is structured according to the time since the entrance in the compartments rather than by infection age. We hope that our work illustrates well the practical potential of such general models. The relation between age-structured and ODE epidemic models exposed in Section 1.2 is known since their very introduction [36], and has been acknowledged multiple times since then [43, 16, 7, 30].

In the most general formulation of a CMJ process, individuals can carry a trait valued in an abstract measure space that encodes all the information about their infection [31, 54]. We have restricted this information to the sequence of compartments visited by each individual, but we could have included some additional details, such as the evolution of the viral load for instance, which could be modeled as a continuous trait following a diffusion. Our result would carry over to the general setting, with a modified limiting equation in Theorem 1 which has already been proposed in [43, Eq. (2.5)] and [44, Chapter IV, Section 1.3]. (Note that none of these works is concerned with the underlying stochastic model.) However, we believe that our current formalism, where the state of an infected individual can be described by a discrete set of compartments, is flexible enough for applications, while being easier to grasp than the general case.

Non-Markov epidemic models have already been investigated, see e.g. [52, 46, 3, 5]. In particular our work shares similarities with a recent series of work in this direction [46, 23, 47]. Let us briefly show how to translate the model in [23] in our framework. Let $\lambda=(\lambda(a);\,a\geq 0)$ be some random function giving the infectiousness of a typical individual in the population. Conditional on $\lambda$, let $\mathcal{P}$ be a time-inhomogeneous Poisson point process with intensity $\lambda$, define

and set

Then, up to the choice of the initial condition and the fact that we make a branching assumption, the model considered in [23] coincides with the model considered here for this specific choice of the pair $(\mathcal{P},X)$. (It coincides exactly with the extension of our model considered in [18].)

The main result in [23] shows that the fraction of individuals in each state ($S$, $E$, $I$, and $R$) converges to the solution of a set of Volterra equations, see their Theorem 2.1. A straightforward computation shows that the latter equations are equivalent to the non-linear version our decorated McKendrick-von Foerster PDE obtained by replacing (5) by (1), with the specific choice of $(\mathcal{P},X)$ described above.

The idea of representing a general branching population by its age structure has a rich history in probability theory [31, 20, 35, 33, 32, 29, 57, 22] and the connection with the McKendrick-von Foerster PDE has been acknowledged several times [20, 29]. In the latter two works, the authors allow for birth and death rates that may depend not only on abundances of each type, but also on the whole age structure of the population. This impressive level of generalization comes at the cost of assuming that the process describing the evolution of the empirical measure on ages and types is Markovian. In particular, birth and death rates are not allowed to depend on past individual birth events. The Markov property then allows the use of a generator for the empirical measure and with some extra finite second moment assumptions on the intensity measure, this approach allows the authors to obtain a law of large numbers and a central limit theorem.

Even if the current work is not as mathematically challenging as that alluded to above, we believe that our point of view does deserve to be highlighted in the current sanitary crisis since it provides both a general modeling framework and an efficient inference methodology. Furthermore, since we ignore finite population effects, our proofs are quite elementary compared to [20, 29] and should be accessible to a much wider audience interested in such a modeling approach. Finally, as far as we can tell, the duality result exposed in Section 2.4 is new and can presumably be extended to more general branching processes where birth and death rates are allowed to be frequency-dependent. In [18], some of the authors of the present work show that this duality result has a natural counterpart in a model with a finite but large population.

This section provides the existence and uniqueness of weak solutions to Eq. (5). Even if similar results are well-known for the time-homogeneous McKendrick-von Foerster equation [30, Chapter 1], we derive them briefly here for the sake of completeness.

In order to motivate our definition of weak solutions, we start by giving a well-known formal resolution of the PDE using the method of characteristics. Fix $a>0$. Let

Then

so that $s\mapsto n(t-s,a-s)$ is conserved along the characteristics, i.e.,

It follows that

where

is the number of new infections at time $t$. We now determine the function $b$. Injecting the previous expression into the “age” boundary condition of the PDE, we obtain a convolution equation for $b$: for every $t>0$

Recall that $\alpha$ denotes the Malthusian parameter defined in (9).

Consider a CMJ with initial shifted law and define

where $Z(t)$ is interpreted as the number of infections between $0$ and $t$. Recall that $R_{0}=\int_{0}^{\infty}\tau(u)\mathrm{d}u<\infty$ guarantees that $B(t)<\infty$ for all $t\geq 0$. Finally, $B$ is non-decreasing and we denote by $\mathrm{d}B$ the Stieljes measure associated to $B$.