A systematic procedure for incorporating separable static heterogeneity into compartmental epidemic models

In this paper, we show how to modify a compartmental epidemic model, without changing the dimension, such that separable static heterogeneity is taken into account. The derivation is based on the Kermack–McKendrick renewal equation.


Introduction
Up to high age, Fred Brauer has been very active in the field of Mathematical Epidemiology.His books [4,5], lecture notes [3] and many papers constitute a valuable heritage.In his paper [2], Fred recognizes that Kermack and McKendrick (KM) introduced in [19] an age-of-infection model that is mathematically represented by a Renewal Equation (RE).Only for very special kernels does the RE reduce to a finite system of ODEs.Or, in other words, compartmental models form a rather restricted subclass of the general KM model.
The paper [9] provides necessary and sufficient conditions for when a delay equation (i.e., a Delay Differential Equation (DDE) or a RE) allows a finite dimensional reduction.Most of that paper concentrates on linear equations, but Section 9.3 is devoted to the nonlinear KM model.
In [8] and in work in progress [1], we formulate the abstract RE that incorporates static heterogeneity of the host population into the general KM model, see Section 5 below.From a modeling point of view, everything is straightforward.But when it comes to analysing the equation, the infinite dimensional character is a great stumbling block.In Section 8.4 of [8], it is shown that under the assumption of separable mixing, (a slightly less restrictive version of (5.3) below), we are back to scalar quantities.This facilitates both the computational and the analytical aspects tremendously, see for instance [22,23,29].
The interpretation of the separability condition can be informally described as follows: whenever the trait/type at the moment of becoming infected is following an a priori given distribution (in particular independently of the trait/type of the infecting individual), newly infected individuals are identical in a stochastic sense and therefore we know how to take averages.From a mathematical point of view, the key point is that various operators have a one-dimensional range.
After the reduction of the abstract RE to a scalar RE, we can further reduce to an ODE system, provided the time-since-infection kernel has the required form.We claim that this seemingly circuitous derivation of compartmental models, that incorporate separable static host heterogeneity, is, due to its systematic character, far more powerful and efficient than a direct approach that starts from the compartmental model describing a homogeneous host population and then adds heterogeneity.Stated differently: the RE formulation is much more amenable to generalization than the ODE formulation.The aim of the present paper is to establish this systematic procedure and to demonstrate its effectiveness.
We restrict to models of an outbreak in a closed population, i.e., we ignore both demographic turnover and loss of immunity.This is quite essential.Indeed, we shall first reconsider the homogeneous KM model and reduce it somewhat differently from the manner described in Section 9.3 of [9].This new reduction involves an integration step that only works in the outbreak situation.And it is this new reduction that easily generalizes to the heterogeneous KM model.
The organization of the paper is as follows.In Section 2, we introduce the original KM model, characterized by a kernel A(τ ), with τ the time elapsed since infection and A the expected contribution to the force-of-infection.We derive, as in [6], a RE for the cumulative force-of-infection.In Section 3, we focus on the special case of a kernel A which is a matrix exponential sandwiched between two vectors.For this special case we deduce the corresponding ODE compartmental system.In Section 4, we explain how the form of the compartmental model derived in Section 3 relates to the standard form.As concrete examples, we consider the elementary SIR and SEIR models (the treatment of more complicated examples is postponed till Section 8).In Section 5, we turn to heterogeneity.Individuals are characterized by a trait x taking values in a measurable space Ω.The kernel is now a function of three variables, τ , x and ξ, with τ the time elapsed since an individual with trait ξ became infected and A being the expected contribution to the force-of-infection on individuals with trait x.Assuming that A is the product of functions a(x), b(τ ) and c(ξ), we derive a scalar renewal equation for the function w such that the cumulative force-of-infection on individuals with trait x equals a(x)w(t).Next we assume that b is a matrix exponential sandwiched between two vectors and reduce to an ODE system.This heterogeneous ODE system only differs from the corresponding homogeneous ODE system in the definition of a function Ψ : R → R. The definition of Ψ involves the functions a, c and the measure Φ describing the distribution of the trait in the host population.The upshot is that one can incorporate separable static heterogeneity into compartmental models by appropriately choosing Ψ. Section 6 is devoted to taking heterogeneity into account in the standard form of a compartmental model.In Section 7, inspired by [13,[21][22][23][24]29], we elaborate the special case that Ω = (0, ∞), Φ is the Gamma Distribution, a(x) = x and c(ξ) is either identically equal to one or equal to ξ. Section 8 is devoted to examples and in the final Section 9 we collect some concluding remarks.

The general Kermack-McKendrick model
Let S be the size of the subpopulation of susceptible individuals, and let F denote the force-of-infection.In a closed population, the incidence is equal to both F S and to −dS/dt, so we have The essence of the KM model is the constitutive equation that expresses F in terms of contributions by individuals that became infected before the current time: Here the one-and-only (apart from N , the total host population size) model ingredient A describes the expected contribution to the force-of-infection as a function of the time τ elapsed since infection took place.In this top-down approach we postpone a specification of the stochastic processes that underly the word expected (in general, these concern both within host processes, in particular the struggle between the pathogen and the immune system, and the between host contact process).Indeed, KM wanted to know what general conclusions can be drawn without providing such a specification.Now imagine that F was negligible in the infinite past.We introduce the cumulative force-of-infection w defined by (2.3) Suppose that S(−∞) = N .Then integration of (2.1) yields t) .
(2.4)So note, incidentally, that one can equivalently characterize w by the relation w = − log(S/N ).
Integration of (2.2) with respect to time over (−∞, t] leads, upon replacing F S by −S ′ and a change of the order of the integrals, to the RE where Ψ(w) := N (1 − e −w ), (2.6) which corresponds to the subpopulation of no longer susceptible individuals, given the cumulative force of infection.
As a side remark, we mention that (2.5) can be considered as a deterministic version of the Sellke construction, as described in, for instance, Section 3.5.2 of [8]. 3 The special case in which reduction to a compartmental model is possible Suppose there exist an integer n, 1×n, n×1-matrices U , V and an n×n-matrix Σ such that A(τ ) = U e τ Σ V, (3.1) then, in a sense, we have a state representation for the (one-time) input-(distributed-time) output map A. The matrix Σ generates the autonomous Markov chain dynamics of the state, the scalar quantity incidence is an input along the fixed vector V and the i-th component of the vector U measures the output, as a contribution to the force-of-infection, of the i-th state.System (4.1) below and the concrete examples (4.5) and (4.8) should clarify this somewhat vague description.
The claim is that the RE (2.5) reduces to an ODE system when (3.1) holds.To substantiate the claim, we define the n-vector valued function Q of time by On the other hand, it follows from (3.1), (3.2) and (2.5) that Combining (3.3) and (3.4) we obtain the closed system of ODE Note that, conversely, given a solution of (3.5), we can define w by (3.4) and verify that w satisfies (2.5).For reasons explained in the next section, we call (3.5) the integrated form of the compartmental model corresponding to Σ, V and U .
For completeness, we like to mention that the assumption (3.1) immediately allows us to calculate basic indices for the KM model as follows: 1. the basic reproduction number is given by and the intrinsic growth rate r is given by its real root.3. the generation time, here denoted by T , is given by

Two ways of formulating compartmental models
We start from (2.1)-(2.2),make assumption (3.1), and introduce to count the individuals, that were infected before time t, on the basis of their state at time t.Then we can deduce, as detailed in Section 9.3 of [9], the ODE with There are good reasons to call this the standard form of the compartmental model corresponding to Σ, V and U .We claim that Q, as introduced in Section 3, is the integral of Y , i.e., In fact, from (4.1) and (4.4), we obtain where we used (2.6).Note that convergence of the integral in (4.4) requires that Y converges to zero at −∞. Differentiating the above equation (4.5) with respect to t, we recover (3.3), i.e., Since w(t) = U Q(t), we arrive at the integrated compartmental model (3.5).
A minor advantage of the integrated form is that the dimension is n, not n+ 1.In the present context, the key favourable feature is that the integrated form extends seamlessly to the separable heterogeneous setting (while the standard form does not).
To illustrate the integrated formalism, we consider the two most basic examples.If we write the standard form of the SIR model as the integrated form reads where and Ψ is defined in (2.6).So here n = 1, V = 1, U = β and Σ = −α.
If we write the standard form of the SEIR model as the integrated form reads where .
So here n = 2 and we have

Incorporating heterogeneity
Let host individuals be characterized by a trait x, with x ranging in a measurable space Ω (the advantage of this somewhat abstract formulation is that x may be a discrete variable, a continuous variable or a mixture of these two possibilities, in the sense that x has both a discrete and a continuous component).The kernel A now has three arguments, τ , x and ξ.The variable x specifies the trait of the individual that is at risk of becoming infected, while ξ and τ specify, respectively, the trait and the time-since-infection of an infected individual.
We want a formalism that captures both the case where Ω is discrete (for instance a finite set) and the case where Ω is continuous.To achieve this, we describe the population composition by a measure Φ on Ω and, for any t, a bounded measurable function s(t, •) such that of the individuals with trait x, a fraction s(t, x) is still susceptible at time t (in other words: s(t, x) is the probability that an individual with trait x is susceptible at time t and s(−∞, x) = 1).
We replace (2.1)-(2.2) by, respectively and When functions a, b and c exist such that we deduce from (5.2) that F is the product of a(x) and a function of time.This motivates us to introduce a function w such that (5.4) Next we integrate (5.2) with respect to time over (−∞, t], while using (5.1) to replace the product F s by the time partial derivative of s.After dividing out the factor a(x) that occurs at both sides, we obtain with Ψ now defined by with Φ the measure that describes the probability distribution of the trait in the host population.Note that from (5.6) we recover the earlier definition (2.6) in the special (homogeneous) case that both a and c are identically equal to 1. Also note that when a is identically equal to one, (5.6) just states that we can simply work with the average value of c.So, as indeed emphasized in Chapter 2 of [8], heterogeneity of infectiousness is, in the large numbers deterministic limit, simple : just take the expected value.
Essentially, (2.5) and (5.5) are the same.When a is not constant, (5.6) differs from (2.6).When b(τ ) = U e τ Σ V, (5.7) we can define Q as before, see (3.2), and retrieve both (3.5) and (3.4).We conclude that in terms of the integrated formulation of a compartmental model, we can incorporate separable heterogeneity by simply redefining the function Ψ.The fact that (5.6) involves an integral is of little to no importance for the theory.But when one wants to study (3.5) numerically, it is a nuissance.As noted before by various authors (cf.[13,[21][22][23][24]29]), in a special case one can replace the integral by an explicit expression.In Section 7 we shall provide the details.6 The standard form : a recipe In Section 8 we present a model for a heterosexually transmitted disease.In that case we have to distinguish between newly infected males and newly infected females.So there are TWO states-at-infection.Here we focus on models in which there is only one state-at-infection, represented by the vector V (but note that V may be a probability vector with several non-zero components, see Example 8.1).
Our starting point is the compartmental system We want to incorporate separable static heterogeneity as described by the trait space Ω, the trait distribution Φ, the relative trait-specific susceptibility a and the relative trait-specific infectiousness c.Here relative means that we single out a representative x in Ω and normalise a and c by requiring that a(x) = 1, c(x) = 1 (note that (5.3) offers the possibility to thus normalise a and c; if Ω is a continuum, it makes sense to choose the mean trait as x).
Then s(t) := s(t, x) denotes the fraction of the individuals with the trait x that escaped infection up to the current time.For arbitrary trait x, the corresponding probability equals because s(t, x) = exp(−a(x)w(t)).Hence the fraction of the total population that is still susceptible is given by Moreover, since a(x) = 1, we have w(t) = − log s(t).Thus knowledge of s is sufficient to determine both s(t, x) and w(t).
On the other hand, in analogy with (4.1), let us define the trait-specific y-variable by and let Y now denote the weighted y-variable: Note that, since the dynamics of infected individuals is independent of the trait and linear, it does not matter whether we apply the weight factor c when computing the output by applying U or right at the start, i.e., immediately after the infection took place.Now we formulate the standard compartmental system which takes into account the static heterogeneity: Claim The heterogeneous compartmental system consisting of (3.5) with Ψ defined by (5.6), has the standard form representation where s(t) = s(t, x), F (t) = F (t, x) and explicitly we have ξ)a(ξ)e −a(ξ)w Φ(dξ).(6.7) Proof It follows from (5.2), (5.3), (6.4) and (6.5) that Differentiation of (3.5) with respect to time yields We identify U Q with − log s , since this is what we obtain when we integrate 7 The Gamma Distribution Let Ω = [0, ∞) and let a(x) = x, i.e., let the trait correspond directly to relative susceptibility.Let Φ be the Gamma Distribution with mean 1 and variance p −1 .In other words, let Φ have density The key feature is that under these assumptions we can evaluate the integral in (5.6) when c is a (low order) polynomial and thus obtain an explicit expression for Ψ.The underlying reason is that we deal with (a derivative of) the Laplace Transform of Φ, which is itself explicitly given by If the trait has no influence on infectiousness, i.e., c is identically equal to 1, we have while if infectiousness too is proportional to the trait, i.e., c(ξ) = ξ, we obtain In [1], we compare and contrast these special cases in terms of R 0 , the Herd Immunity Threshold and the final size.Here we limit ourselves to the observation that (3.5), with either (7.3) or (7.4), enables modelers to study rather easily the impact of separable heterogeneity on the dynamics of their favourite compartmental model, cf.[22].Now recall that the standard form involves s, the value of s in a representative point x.In the special case of the Gamma Distribution, we can work with s tot instead of s and replace (6.6) by The derivation of (7.5) starts by choosing x = 1 (= mean trait), a(x) = 1 and rewriting (6.3) in this special case as s tot = Φ(− log s). (7.7) Since Φ is given explicitly by ( 7.2), one can express derivatives of Φ in terms of Φ itself.And derivatives correspond to incorporating powers of the variable in the function that is Laplace transformed.We refer to [23] for an early derivation of (7.5) and for references to still earlier work.

Example 1
The diagram depicted in Figure 1 gives a concise representation of the compartmental model that we consider in this section as a slightly more complex example.Our first aim is to illustrate how one obtains the model ingredients

model with asymptomatic infection and quarantine
The index 1 denotes asymptomatic individuals, the index 2 symptomatic individuals.As the symptoms get noticed, a diagnosis is possible for symptomatic individuals and subsequently they may be put into quarantaine, i.e., enter Q.The two types of individuals occur with ratio p : 1−p , with p a parameter.So the state-at-infection/birth is a probability vector V with two non-zero components.Immediately following infection an individual is Exposed but not yet Infectious.The sojourn times of the various compartments are all exponentially distributed with a parameter specified in the diagram (in the form of a name label, so as a parameter, not as a number with dimension 1/time).The compartments R and Q aid the bookkeeping, but their contents is irrelevant for future dynamics, so we do not incorporate them into the population state vector.
As is customarily done, we use the same characters to denote a compartment and the contents of this compartment.Define the 4-vector Y by the vector U by and the matrix Σ by Either by a direct consideration, or by verifying that and applying the formula (3.6), we obtain Similarly we obtain from (3.8) that the generation time is given by Example 1, continued Next let us introduce immune system related heterogeneity in the sense that we distinguish between standard individuals, which we label 1, and partly immune individuals, which we label 2. The relative susceptibility of type 2 individuals is given by the parameter ǫ 1 and the relative infectiousness by the parameter ǫ 2 .In the present context it does not matter whether the immunity results from an earlier outbreak or from vaccination.But we assume it exists before the outbreak that we model is initiated or, in other words, that it does not result from control measures during the outbreak.
Let N = N 1 + N 2 with N 1 and N 2 the size of the subpopulation of individuals of type, respectively, 1 and 2. Then Φ has components N 1 /N and N 2 /N .We may choose 1 as x and let a have components 1 and ǫ 1 and let c have components 1 and ǫ 2 .It follows that One can now consider (3.5) or (6.6) with the above definitions of U , Σ, V and Ψ (resp.Ψ ′ ) and study, for instance, numerically how peak size is influenced by the parameters ǫ 1 , ǫ 2 and N 1 (for given N ).
Note that when we choose ǫ 1 = ǫ = ǫ 2 , this example allows for an alternative interpretation: as a result of a control measure, individuals reduce their social activity with a factor ǫ, but only a fraction N 2 /N complies, the complementary fraction does not reduce its social activity.

Example 2
In the formulation of our results, we have restricted to the situation in which one vector V and one vector U suffice.In [9], the generic result actually concerns systems of REs for which one needs as many vectors V and U as the number of components of the system.Here we illustrate how this works by concentrating on the epidemiologically relevant example of a heterosexually transmitted disease (or a disease transmitted by a vector).
Let S i (i = 1, 2) be the size of the susceptible population, where the index 1 denotes males and the index 2 females (or the host and the vector, respectively).We postulate that where Define the cumulative force of infection as Then we obtain the renewal equation: with where N 1 denotes the total size of the male population and N 2 is the total size of the female population.
For the compartmental case, we have Here U i , V i are 1 × n i , n i × 1 matrices and Σ i is an n i × n i matrix.It seems to make sense to assume that n 1 = n 2 and perhaps Σ 1 = Σ 2 .But if within host processes are different for males and females, we should allow for n 1 = n 2 and Σ 1 = Σ 2 .This is anyhow a reasonable assumption for the host-vector situation.Define vectors Q 1 and Q 2 by This is the integrated version of the following standard form: (8.17)

Example 2, continued
We use the same trait x ∈ Ω to characterize males and females.The trait x may represent promiscuity.We start from (5.2) but now F is a 2-vector and A a 2 × 2-matrix with zero's on the diagonal.We assume then which is of the form (8.13) with A = 0 b 2 b 1 0 and Now we assume that (8.16) holds.Then (8.17) holds, with now Ψ as defined above.
In our presentation of this example we avoided to discuss the consistency requirement that there are, in total, as many contacts of males with females as there are contacts of females with males.Partly we did so because transmission risk may be asymmetric, which obviously impairs the symmetry requirement.But the more important reason is that our aim here is just to illustrate the flexibility of the bookkeeping framework (and not to analyse a model of a heterosexually transmitted disease).

Concluding remarks
On November 13, 2022, the KM paper [19] had, according to Google Scholar, 12.590 citations.On that same date the Royal Society listed 60.773 downloads since the beginning of 1997, when the paper became available online.No doubt there are among the authors that cite [19] some who actually read the paper and who refer to the general age-of-infection model, see e.g.[2,6,8,15,21,23,24,28,29]. In the big majority of cases, however, it is explicitly stated that [19] introduced the SIR compartmental model and implicitly suggested that that is it.This both reflects and reinforces an incessant misconception in the mathepi community at large, viz., that [19] is just about the SIR compartmental model.
In fact, as shown in [9], any compartmental model in which the (probability distribution of the) state-at-infection is described by a given fixed vector V corresponds to a (very) special case of the renewal equation.The compartmental system is coded by the triple (V, Σ, U ), with the vector V describing the initial state, the matrix Σ specifying the rates at which state transitions occur and the vector U describing the contribution of the various states to the force of infection.(When there are several possibilities for the state-at-infection, one needs a system of renewal equations and a corresponding number of vectors V and U , see [9] and Example 8.2.) It has long been recognized that host heterogeneity can have a big impact on epidemic dynamics, see, e.g., [7,8,14,15,18,23,24,31,32].In general, the incorporation of heterogeneity necessitates the introduction of kernels and leads to infinite dimensional dynamical systems.In this context too, the question arises whether or not one can reduce to a finite system of ODE.In [23,24] and in the more recent Covid-triggered papers [13,21,22,25] a restricted form of static heterogeneity, in which the traits of the two individuals involved in a contact are assumed to have independent influence on the likelihood of contact and/or transmission, is considered.The impact of such heterogeneity on the outbreak dynamics, is, of course, most easily studied if the heterogeneity is captured by a modification of the compartmental model one is interested in.
Here we have shown that it is straightforward to derive the desired modification if one first formulates the heterogeneous version of the KM RE and relates it to an integrated version of the compartmental model.Once the modified integrated version is available, it is easy to derive the modification of the standard form by differentiation.The end result only involves (the derivative of) a function Ψ from R to R defined in terms of the distribution of the trait and the functions that describe how susceptibility and infectiousness depend on the trait (so one can use the end result without any reference to renewal equations).
From a more general theoretical point of view, our methodology is in the spirit of [11,12] and the much older references in there.In a similar manner, [23] and [24] build on the ideas of G.P. Karev as described in [16,17] and the much older references in there.
In [6], it is shown how to extend the RE formulation to models incorporating demographic turnover and/or temporary immunity.Such models allow for endemic steady states.They do not have an integrated version, for the simple reason that the integrals are bound to diverge.At present it remains unclear whether or not one can include separable static heterogeneity in such models via a simple modification (the authors are not very optimistic ...).
We refer to [29] for an up to date Covid inspired account of the role of heterogeneity in outbreak dynamics.In there it is emphasized that, on top of persistent static heterogeneity, dynamic heterogeneity may play a major role in damping the overshoot of the herd immunity threshold (in much the same way as a prey refuge dampens prey-predator oscillations).An additional variable h is introduced to capture the dynamic heterogeneity.
Our message in this paper is a bit equivocal.On the one hand, our aim was to show how separable static heterogeneity can be easily incorporated into compartmental models.On the other hand, we want to emphasize that the RE formulation is far more general and flexible and that the predominance of compartmental models is rather detrimental for epidemic modeling.No doubt the fact, that user friendly tools for the numerical study of RE are lacking, contributes to their unpopularity.In order to end at a positive note, we call attention to two recent developments: i) discrete time models, see [10,20,27] and ii) pseudospectral approximation, see [26].

-( 4 . 2 )
is the standard representation of the compartmental model specified by the diagram.So if we define Q by (4.4) then we obtain the integrated representation (3.5).