key: cord-0923488-cmtmp0xv
authors: Diekmann, O.; Othmer, H. G.; Planque, R.; Bootsma, M. C.
title: On discrete time epidemic models in Kermack-McKendrick form
date: 2021-03-26
journal: nan
DOI: 10.1101/2021.03.26.21254385
sha: 7c4432514783d768d2d240c60ba47de57821a587
doc_id: 923488
cord_uid: cmtmp0xv

Surprisingly, the discrete-time version of the general 1927 Kermack-McKendrick epidemic model has, to our knowledge, not been formulated in the literature, and we rectify this omission here. The discrete time version is as general and flexible as its continuous-time counterpart, and contains numerous compartmental models as special cases. In contrast to the continuous time version, the discrete time version of the model is very easy to implement computationally, and thus promises to become a powerful tool for exploring control scenarios for specific infectious diseases. To demonstrate the potential, we investigate numerically how the incidence-peak size depends on model ingredients. We find that, with the same reproduction number and initial speed of epidemic spread, compartmental models systematically predict lower peak sizes than models that use a fixed duration for the latent and infectious periods.

the mistake will have had its day. 48 To put some flesh on the bones, we show that the qualitative behaviour of 49 the discrete time model is the spitting image of the qualitative behaviour of the 50 continuous time model. In order not to ignore the popularity of compartmental 51 variants (which is unwarranted, in our opinion, as there is neither evidence that 52 the length of, for instance, the infectious period is exponentially/geometrically 53 distributed nor that infectiousness is constant during this period), we put them 54 in the spotlight. To conclude, we illustrate the relevance (in particular for public 55 health policy) of the lesser known members of the Kermack-McKendrick family, 56 by investigating numerically how the peak of the incidence curve varies among 57 members that are identical with respect to both the initial growth rate ρ and 58 the basic reproduction number R 0 , but differ in assumptions about the duration 59 of the exposed and the infectious period. 60 To avoid misunderstanding, we now clarify what is, and what is not, stochas-61 tic in the models formulated and analyzed below. All models are deterministic 62 at the population level. (One can think of them, in Kurtz spirit, see [16] , as 63 the large initial population size limit of a stochastic model for finitely many 64 individuals.) infectiousness A(τ ) at time τ after becoming infected equals βe −γτ . This re-74 flects that after time τ has elapsed, the probability to be still infectious equals 75 e −γτ . At the population level, we translate this into: a fraction e −γτ of those 76 infected at time t is still infectious at time t + τ . In the discrete time setting, 77 we should replace the exponential distribution by the geometric distribution. 78 As a final elucidation we mention that the parameter β and the function A(τ ) 79 implicitly incorporate information about the contact (between individuals) pro-80 cess that underlies transmission. A key feature is that contacts are assumed 81 to be uniformly at random (so neither spatial nor age nor social structure are 82 incorporated). 

So a susceptible is infected with probability 1 − (1 − p) N rather than with 125 "probability" pN .

From a numerical point of view, the exponential has the disadvantage of 127 being expensive in terms of calculation costs. It may therefore be tempting 128 to reduce the step size in order to work safely with the linear approximation. 129 We actually wonder whether solving the relevant ODEs with for instance a is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint Often this is done after first restricting to a parameterized family of functions 147 A (but see [9] for an alternative methodology.) Note that, with N denoting the 148 population size, we have

and that, in the initial phase of the epidemic, the distribution of the generation-150 interval (see [17, 18] , and the references given in there) has as density the renor-151 malized (to have integral 1) function A. 152 The incidence at time t is given by Λ(t)S(t). Under the "permanent immu-153 nity and no demographic turnover" assumption, the incidence equals −S ′ (t), 154 cf. (2.2). Substituting this into (3.1), we deduce by integration (and upon 155 changing the order of integration) the identity

The discrete time counterpart reads

where now A j is the expected contribution to the cumulative force of infection 158 in (t, t + 1] of an individual who itself became infected in the time window is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint 

that follows by combining (2.2) with (3.2) and incidence equal to −S ′ , see [3, 6] . expressing that in the infinite past all host individuals were susceptible.

To derive (3.3), first note that iteration of (2.1) yields, if (3.7) holds, the 

If we use this last identity in (3.8), divide both sides of (3.8) by N and adopt is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted March 26, 2021. ; https://doi.org/10.1101/2021.03.26.21254385 doi: medRxiv preprint into (3.3) and assume that x is so small that it makes sense to replace the 202 exponential by the zero'th and first order terms of its Taylor expansion. This 203 yields the linearized equation

We define

and interpret, based on the last identity, R 0 as the expected number of secondary 206 cases caused by a primary case in a totally susceptible host population.

In order to show that positive solutions of (4.2) grow when R 0 > 1, but

By substitution of (4.4) into (4.2) we find that x defined by (4.4) is indeed a 210 solution if and only if λ is a real root of the discrete time characteristic equation

known as the Euler-Lotka equation. The non-negativity ofà k , k = 1, 2, . . ., 212 guarantees that (4.5) has at most one real root ρ and that it does indeed have a 213 real root when the right hand side assumes a value bigger than one for some real 214 λ, so in particular when R 0 > 1 (when R 0 < 1 andà k has power-like behaviour 215 for k → ∞, the value of the right hand side may jump from a value less than one 216 to infinity when λ is decreased; whenà k = 0 for large k, this cannot happen and 217 ρ exists). Readers who wonder (or even worry) about the potential importance 218 of complex roots can consult [6, Section 8.2] and the references given there, to 219 be eased.

So we see that a key point is that the right hand side of (4.5) is a monotone . 249 We conclude that at the level of theory, there is an exact parallel. We shall use the standard notational convention (or should one say "ambigu-253 ity"?) that a compartment and its contents are denoted by the same symbol. 254 We start with SIR and after that generalize to SEIR, hoping that these two is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint These assumptions lead to the system of recurrence relations

We show that the system (5.1) may be reduced to the scalar recurrence (3.3) 269 by choosing theà k appropriately:

The first step corresponds to the derivation of (3.8): by iteration of the first 271 equation of (5.1) we obtain

Rewriting the second equation of (5.1) as

we obtain by summation the identity

and by substitution of this identity repeatedly at the right hand side,

Finally, substitution of this last identity in the formula for S(t + 1) above yields 276 (3.8).

Conversely, starting from (3.3) withà k given by (5.2), we easily recover

. CC-BY 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint 

We emphasize that if one replaces e −βI(t) by 1 − βI(t) , the reduction to a 283 higher order scalar recursion relation fails (we invite readers to convince them-284 selves of this fact)!

In order to capture a latent period, we next change the assumptions. Upon 286 infection, an individual now enters the compartment E of exposed (i.e., infected 287 but not yet infectious) individuals. When the length of the latent period is 288 geometrically distributed with parameter γ, we have to replace (5.1) by 289 S(t + 1) = e −βI(t) S(t), 

(with the convention that the sum equals zero when the upper index does not 301 exceed or equal the lower index). The parametersà k are again defined by (3.5).

And whenà k has the form defined by (3.5) and (5.5), then (5.4) follows from is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint compartmental epidemic models to a scalar higher order recursion. See Section 307 9.3 of [7] for a detailed elaboration of the continuous time case. 308 We conclude that there is a multitude of compartmental models that corre- 

Note that (6.1) is based on the debatable assumption that at the day of its de- is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint For feasibility, we want a finite dimensional system. To achieve this, we 350 make the very reasonable assumption that the indices j for which A j is strictly Alternatively we might start from (3.9) and choose as before

which corresponds to the incidence in time window t + 1 − j. This leads to the update rules

X j (t + 1) = X j−1 (t) for j > 2. (7.8) In this formulation too, we can allowà k to depend on time t.

This seems a good moment to point out that the use of labels like 'exposed' is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The ongoing outbreak of Covid-19 generates much interest in peaks, largely 391 because of concern that hospitals may be overwhelmed with patients, leading 392 to healthcare breakdown. As far as we know, there is no analytical method to 393 determine the height and timing of the peak from model parameters (except, 394 perhaps, in the oversimplified SIR system of differential equations). So one has 395 to rely on numerical calculations.

The key question addressed in this section is: how much is peak height 397 influenced by model details? Here, we systematically compare the discrete-time 

As initial condition we took a short history of decreasing fractions of sus-408 ceptibles, reflecting an exponential increase in new cases at the rate ρ. We is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted March 26, 2021. ; Figure 2 : Comparing incidences between two types of models. In the block model, the lengths of the latent and infectious periods are deterministic (fixed for all individuals); in the SEIR model, the lengths of these periods are stochastic (independently exponentially distributed with identical parameters, respectively, 1/T E and 1/T I ). See Appendix for details. Top row: maximum incidence of the block model (with deterministic periods; left), SEIR model (with stochastic periods; middle) and the relative ratio between the two (i.e., block−SEIR SEIR ; right), as a function of T E , the (actual resp. expected) time individuals are exposed and T I , the (actual resp. expected) time individuals are infectious. Models were compared after ensuring that they have the same R 0 and initial speed ρ. Note that the incidence of the deterministic model always reaches a higher peak within the ranges of T E and T I considered, by about 8-15%, than the corresponding SEIR model. Middle and bottom rows: example simulations with the deterministic model (blue) and corresponding SEIR model (red) with the same R 0 and ρ. The middle row corresponds to the parameters at which the ratio of peak heights is minimal, (T E , T I ) = (3, 4); the bottom row to when this ratio is maximal, (T E , T I ) = (6, 4). One can clearly see that the incidence grows initially at the same rate.

14 . CC-BY 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted March 26, 2021. ; 6, 4) . The corresponding incidences can be found in Figure 2 . Note that although the peak incidences are not so very different (see Fig. 2 ), there is a large difference in the fraction of infected-andnot-yet-removed individuals, due to the comparatively fatter tail in the expected future contribution to the force of infection in the SEIR model.

• the difference is, for reasonable parameter values, in the order of 10%.

(Note that, since compartmental models have fatter tails, they need, for 418 given R 0 , to have an earlier peak of infectiousness in order to have the same ρ. 419 This is clearly visible in Figures 2 and 4. ) 420 After the first conclusion emerged, we aspired to find a somewhat mecha-421 nistic explanation. This led to the following observation. Roughly speaking, 422 an outbreak reaches its peak when S is reduced to the level corresponding to 423 R 0 = 1. How many more cases there will be after the peak depends largely on time. This is illustrated in Figure 3 . And as reservoir size correlates with peak 435 size, we should expect lower peaks for compartmental models, exactly as found 436 in our numerical results. 437 Just to elucidate that the higher-peak-phenomenon matters in Covid-19 con-438 text, we chose, on the one hand, the parameters A k as integrals over one day SEIR models that has both R 0 and ρ equal to these quantities for the Weibull.

The results of a comparison are presented in Figure 4 . The peak heights differ 443 5 to 10%. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted March 26, 2021. ; Figure 4 : Comparing a model with parameters corresponding to the Weibull distribution derived from SARS-CoV-2 generation-interval data in [8] to the SEIR model with the same reproduction number R 0 and the same initial growth rate ρ. Top row: left: relative ratio between the maximal incidences (−1, as in Figure 2 ) as a function of T E , the expected length of the exposed period (note that the SARS-CoV-2 model does not depend on T E , see the Appendix). right: expected contribution to force of infection for the SARS-CoV-2 model. Middle row: left: incidence of the SARS-CoV-2 model (blue) and the SEIR model (red) with T E = 1. At this T E , the relative difference in the peak incidences are maximal. right: the expected contribution of the SEIR-model for T E = 1. Bottom row: same, but now for T E = 3, at which the difference in peak incidence is minimal. For T E > 4.12, the SEIR-model cannot be parameterized to have the same R 0 and ρ as the SARS-CoV-2 model. For further details, see the Appendix.

. CC-BY 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint Here we introduced a discrete time version that has many advantages:

• the generality and flexibility is retained;

• computing the epidemic time course is super easy;

• the time step can be adjusted to the time interval between data points 455 (e.g., one day or one week).

In Section 8 we showed that precise assumptions about the latent and infec-457 tious period matter for predicting the peak of the incidence curve, a quantity of 458 interest from a public health perspective. So, we claim, the generality matters 459 for practical issues and is not just an academic fancy. Of course in a practical 460 context all kinds of heterogeneity (e.g., reflecting age) matter as well. These 461 have been neglected here, but in the text book [6] is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint In Figure 2 we numerically compare two types of models that have the same initial rate of increase ρ and basic reproduction number R 0 but different expected contributions to the force of infection A k . Both models are defined by prescribing the values of the rescaled expected contribution to the force of infection,Ã k , where k is the number of days that have elapsed since becoming infected. Let us choose T E , the number of days an individual is exposed but not yet infectious, and T I , the number of days an individual is infectious. In the 'block model', in which these periods are assumed to be deterministic, so the same for all individuals, we set

In this way, ∑ ∞ k=1Ã k = R 0 . The initial exponential rate ρ with which the 531 epidemic spreads is the solution to (4.5).

With R 0 , ρ, and T E , we can set up now the corresponding SEIR model.

Using the explicit expression of the A k for the SEIR model (5.5), we can de- Lastly, the link between γ and T E is simply γ = 1 T E . So in all, T E , ρ and R 0 539 define the following parameters for the SEIR model,

so that theà k are given by

) k−l−1 .

. CC-BY 4.0 International license It is made available under a perpetuity.

is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint The numerical simulations were carried out using (3.9), using an initial con-544 dition in which the epidemic has started increasing at rate ρ for six days: Then we setà k = R 0 g k , k = 1, 2, . . .. The initial rate of increase is again 551 estimated using 4.5, and gives ρ = 1.1919.

The corresponding SEIR-model is now given by using α, β and γ as in 553 (10.3) and defining theà k as before.

The initial condition is the same as in the previous illustration. is the author/funder, who has granted medRxiv a license to display the preprint in (which was not certified by peer review) preprint

The copyright holder for this this version posted March 26, 2021. ; https://doi.org/10.1101/2021.03.26.21254385 doi: medRxiv preprint

On the formulation of epidemic models (an appraisal of Kermack and McK-478 endrick)

Limiting behaviour in an epidemic model

Mathematical Tools for

Finite dimensional state rep-488 resentation of linear and nonlinear delay systems

Quantifying SARS-CoV-2 transmis-492 sion suggests epidemic control with digital contact tracing

Nonparametric Estimation under Shape 495

Constraints. Cambridge Series in Statistical and Probabilistic Mathematics

A brief history of R 0 and a recipe for its calculation

The law of mass-action in epidemiology: a his-500 torical perspective

Discrete epidemic 504 models with arbitrary stage distributions and applications to disease con-505 trol

A contribution to the mathematical 507 theory of epidemics

Proposal of a recursive compartment model of 509 epidemics and applications to the covid-19 pandemic

A discrete Kermack-McKendrick model and why 512 it is better adapted to Covid-19 than standard SIR. submitted

Approximation of Population Processes

A practical 517 generation-interval-based approach to inferring the strength of epidemics 518 from their speed

Forward-looking serial intervals