J. Math. Biol. (1990) 28:365-382
dournalof
Mathematical
Biology
© Springer-Verlag 1990
On the definition and the computation of the basic
reproduction ratio R o in models for infectious diseases
in heterogeneous populations
O. Diekmann 1'2, J. A. P. Heesterbeek 1, and J. A. J. Metz 2,1
Centre for Mathematics and Computer Science, P.O. Box 4079, 1009 AB Amsterdam,
The Netherlands
2 Institute of Theoretical Biology, Leiden University, Kaiserstraat 63, 2311 GP Leiden,
The Netherlands
Abstract.
The expected number of secondary cases produced by a
typical
infected individual during its entire period of infectiousness in a completely
susceptible population is mathematically defined as the dominant eigenvalue
of a positive linear operator. It is shown that in certain special cases one
can easily compute or estimate this eigenvalue. Several examples involving
various structuring variables like age, sexual disposition and activity are
presented.
Key words: Epidemic models -- Heterogeneous populations -- Basic repro-
ductive number -- Invasion
1. Introduction
Suppose we want to know whether or not a contagious disease can "invade" into
a population which is in a steady (at the time scale of disease transmission)
demographic state with all individuals susceptible. To decide about this question
we first of all
linearize,
i.e. we ignore the fact that the density of susceptibles
decreases due to the infection process. It has become common practice in the
analysis of the simplest models to consider next the associated
generation process
and to define the
basic reproduction ratio I Ro
as the
expected number of secondary
cases
produced, in a completely susceptible population, by a
typical
infected
1 Quite often this is wrongly called the reproduction 'rate'. Many authors use 'reproductive' instead
of 'reproduction', but, as was pointed out to us by M. Gyllenberg, the latter is grammatically more
correct. We have followed the advice of I. Nfisell to use 'ratio' rather than 'number' in order to
emphasize that R o does not even have a quasi-dimension (R o ~ cases/case!)
366 o. Diekmann et al.
individual during its entire period of infectiousness. The famous
threshold
criterion
then states:
the disease can invade if R0 > 1, whereas it cannot if Ro < 1.
It is the aim of this note to demonstrate how these ideas extend to less simple
(though probably still highly oversimplified) models involving
heterogeneity
in
the population and to explain the meaning of "typical" in the "definition" of Ro
above. Subsequently we shall deal with the actual computation of Ro in certain
special cases, in particular the so-called "separable" or "weighted homogeneous
mixing" case.
2. The definition
Let the individuals be characterized by a variable 4, which we shall call the
h-state variable (h for heterogeneity). Let S = S(~) denote the density function of
susceptibles describing the
steady demographic state in the absence of the disease
(in order to avoid confusion we emphasize that S is not a probability density
function; that is, its integral equals the total population size in the steady
demographic state, and not one). Define A(r, ~, r/) to be the expected infectivity
of an individual which was infected r units of time ago, while having h-state t/,
towards a susceptible which has h-state 4. The expected number of infections
produced during its entire infective life by an individual which was itself infected
while having h-state t/is then given by
fo;o
s(~) A(r, 4, ~) dr d~,
where O denotes the h-state space, i.e. the domain of definition of 4. We may call
this quantity the next generation factor of ~.
Remark.
In order to have a unified notation for various cases we write integrals
to denote sums whenever f2 is discrete (completely or just with respect to some
component of 4). A precise mathematical justification involves a dominant
measure and Radon-Nikodym derivatives.
Since the new cases arise, in general, with h-states different from t/, these
numbers do not tell us exactly what happens under iteration, i.e. in subsequent
generations (although it is clear that the supremum with respect to ~ yields an
upper estimate for Ro).
So we abandon the idea of introducing an infected individual with a
particular well-defined h-state and start instead with a "distributed" individual
described by a density 4~. The next-generation operator
K(S)
defined by
(K(S)~b)(~) = S(¢) A(r, ~, t/) dr q~(~/) d~/ (2.1)
tells us both how many secondary cases arise from ~b and how they are
distributed over the h-state space. Ignoring the task of writing down conditions
The basic reproduction ratio R o 367
on S and A which guarantee that
K(S)
is a bounded operator on LI (f2), we note
that the next generation factor of q~ is simply the Ll(f2)-norm of
K(S)c~,
i.e.,
f S(~)f~fo°°A(~,~,~)d~q~(q)dtld~
(note that we do not have to write absolute value signs since the biological
interpretation requires all functions to be positive). If we take the supremum of
the next generation factor over all ~b with I] 4~ [1 = 1 we obtain, by definition, the
operator norm of
K(S).
This yields an upper estimate for Ro for the same reason
as above: the distribution with respect to ~ is changed in the next generation and
consequently the factor of q~ does not predict accurately what happens under
iteration.
As a concrete example consider a host-vector model. (For the purpose of
exposition we adopt here the strict version of the law of mass action even though
this does not necessarily yield a good model in this case, see, e.g., Bailey (1982),
Chap. 7.) Taking O = {1, 2} and
S~ A(z, i,j) dT = a o
with ag > 0 if and only if
i #j we find that
K(S)
is represented by the matrix
a,;xl)
and the operator norm is
max{alzS~, a21S2}.
These two numbers correspond to
vector ~ host and host ~ vector transmission, respectively. No matter which of
the two is the larger one, in the next generation it is necessarily the other of the
two numbers which is the relevant
factor. 2
Therefore the operator norm of
K(S)
is not a good definition of Ro. Since
a~2Sla2~ Sz
is the two-generation factor, the
average
next generation factor is
X/al2 SI a21 $2.
How can we define such a quantity in general?
After m generations the magnitude of the infected population is (in the linear
approximation)
K(S)"c~
and consequently the per-generation growth factor is
[[g(s)m[[ l/m.
We
want to know what happens to the population in the long run,
so we let m ~ ~. The so-called
spectral radius
(Schaefer, 1974)
r(K(S))
is defined
by
r(K(S))
=
inf
I]g(s)m[[ 1/m=
lim
[[g(a)m[[ 1/m.
(2.2)
rn~>l m~co
Starting from the zeroth generation ~b, the mth generation
K(S)mc~
converges to
zero for m ~ oo if
r(K(S))
< 1 whereas it can be made arbitrarily large by a
suitable choice of ~b and m when
r(K(S))
> 1. Moreover, the
positivity
guarantees
that in the latter case there is not really a restriction on ~b. Indeed,
K(S)
is a
2 One could argue, as MacDonald did (see Bailey 1982, p. 100 and the references given there), that
one should consider the average number of cases in the host population arising from one case in
the host population via vector cases. From our point of view this amounts to looking
two
generations ahead. Indeed one obtains exactly MacDonald's result if one writes out
al2Sla21 S 2
in
terms of biting rates etc
368 O. Diekmann et al.
positive operator (i.e. nonnegative functions are mapped onto nonnegative
functions) and one can specify conditions on K (for example compactness) that
guarantee that
r(K(S))
is an eigenvalue (Schaefer 1960, 1974), which we shall call
the
dominant eigenvalue
(since ]hi ~<
r(K(S))
for all h in the spectrum of
K(S))
and
denote by Qa. Under minor technical conditions on A and S (see Remark 4
below) one has in addition that
K(S)m~ ) "~ C(~))Ond (Pd
for m --* or, (2.3)
where q~a is the corresponding eigenvector (which is positive) and c(~b) a scalar
which is positive whenever ~b is nonnegative and not identically zero. So after a
certain period of transient behaviour each generation is (in an approximation
which improves as time proceeds) 0a times as big as the preceding one and
distributed over h-state space as described by ~ba.
If we rephrase this as: "the
typical
number of secondary cases is Qa", we are
ready for the
Definition. Ro =
r(K(S))
= 0a = the dominant eigenvalue of
K(S).
With this definition the threshold criterion remains valid, as can be verified
as follows. The threshold criterion relates the generation process to the develop-
ment of the epidemic in real time, both in the linearized version. The linearized
real time equation is
for0
i(t, 3) = S(¢) A(z, ~, r/)i(t - z, r/) dz dr~,
(2.4)
where
i(t, 3)
is the rate at which susceptibles with h-state ~ are infected at time
t. This equation has a solution of the form
i(t,
3) = eatS(C) if and only if ~b is an
eigenvector of the operator Kz defined by
;of0
(Kz$)(¢) = S({) A(% 3, r/) e-~ dz q~(r/) dr/ (2.5)
with eigenvalue one. Positivity arguments can be used to show that among the
set of such h with
largest real part
there is a real one, which we shall denote by
ha (and the corresponding eigenvector by $a). Monotonicity arguments then
imply that
ha>0 ¢~ R0>l and ha<0 ~ Ro<l.
(Heijmans 1986, Sects. 4-6, works this out in detail for a different but similar
example and gives appropriate references. Hethcote and van Ark (1987) contains
a proof for the finite dimensional case.)
Remarks.
(1) Whereas Ro is a number, ha is a rate.
(2) Note that ha and ~b a describe the growth and the h-state distribution in the
exponential phase of an epidemic, when the influence of the precise manner in
which the epidemic started has died off and the influence of the nonlinearity is
not yet perceptible.
The basic reproduction ratio R 0 369
(3) In order to guarantee that
any
introduction of infectivity in the population
leads to an epidemic when R0 > 1 we need to make an
irreducibility
hypothesis
(Schaefer 1974).
(4) Similar parameters determine the asymptotic behaviour in branching
processes with general state space. See Jagers and Nerman (1984); Mode
(1971).
(5) To obtain a complete model one has to specify the demographic processes,
and in particular how per capita birth- and death rates are affected by the
disease. If one makes the obvious assumption that the disease leads to a lower
(or equal) birth rate and to a higher (or equal) death rate one can use the
linearized problem to obtain upper estimates for the nonlinear problem. Thus
one can prove, in general,
global
rather than local stability for R0 < 1. Or, in
other words, endemic states are impossible when Ro < 1.
(6) Let S denote the susceptible population in a steady endemic state. Then
necessarily
r(K(S))
= 1. See Example 4.3.
(7) We have restricted our attention to the bilinear case. However, replacing
S(¢) in the definition (2.1) of
K(S)
by
h(S(~))
or
S(~)/(I+SaS(~I)dq)
or
something similar does not make any essential difference. See Examples 4.1 and
4.2 below. Note that for the invasion problem one will always have an expres-
sion involving the (known) function S only. Of course things are different if
one wants to characterize endemic states, like in Remark 6 above.
3. Computational aspects: easy special cases
3.1. Separable mixing rate
To compute the dominant eigenvalue of a positive operator is, in general, not
an easy task. However, there is one special case in which the task is trivial:
when the operator has one-dimensional range. Biologically this corresponds to
the situation in which the distribution (over the h-state space f2) of the
"offspring" (i.e. the ones who become infected) is
independent
of the state of
the "parent" (i.e. the one who transmits the infection). In this case we speak of
a separable mixing rate
or separable infectivity and susceptibility, or (separably)
weighted homogeneous mixing.
Assume that
then
fo ~ A(z, ~, tl) cl~ = a(~)b(q)
(3.1)
(K(S)~b)(~) =
S(~)a(~) fa b(tl)c~(q) dtl.
(3.2)