Real-time growth rate for general stochastic SIR epidemics on unclustered networks

0.91 MB·English

by  Lorenzo Pellis

Real-time growth rate for general stochastic SIR
epidemics on unclustered networks
5 Lorenzo Pellis1,∗, Simon Spencer2, Thomas House3,1
1
0
2

r
p
A
Abstract

0
2 Networks have become an important tool for infectious disease epidemiol-
ogy. Most previous theoretical studies of transmission network models have

]
E either considered simple Markovian dynamics at the individual level, or have
P focused on the invasion threshold and final outcome of the epidemic. Here,
.
o we provide a general theory for early real-time behaviour of epidemics on
i
b large configuration model networks (i.e. static and locally unclustered), in
-
q particular focusing on the computation of the Malthusian parameter that
[ describes the early exponential epidemic growth. Analytical, numerical and

2 Monte-Carlo methods under a wide variety of Markovian and non-Markovian
v
assumptions about the infectivity profile are presented. Numerous examples
4
2 provide explicit quantification of the impact of the network structure on the
8 temporal dynamics of the spread of infection and provide a benchmark for
4
validating results of large scale simulations.
0
.
1 Keywords: Epidemic, Malthusian parameter, Basic reproduction number,
0
5 Configuration model, Branching process
1 2010 MSC: 05C82, 60K20, 92D30
:
v
i
X
1. Introduction
r
a
The field of infectious disease epidemiology has benefitted from the use of
networks both as conceptual tools and as a practical representation of inter-
∗
Corresponding author: [email protected]
1Warwick Infectious Disease Epidemiology Research centre (WIDER) and Warwick
Mathematics Institute, University of Warwick, Coventry, CV4 7AL, UK.
2WarwickInfectiousDiseaseEpidemiologyResearchcentre(WIDER)andDepartment
of Statistics, University of Warwick, Coventry, CV4 7AL, UK.
3School of Mathematics, University of Manchester, Manchester, M13 9PL, UK.
Preprint submitted to Mathematical Biosciences April 21, 2015
action between the agents involved in the spread of infections (Bansal et al.,
2007; Danon et al., 2011). From a theoretical perspective, they have been
successfully used to obtain important insight in the behaviour of epidemics
in idealised populations. Most analytical results, however, have either been
derived in the specific case of a Markovian SIR model involving constant
infection and recovery rates (e.g. Diekmann et al., 1998; Eames and Keeling,
2002;Volz,2008;Ball and Neal,2008;Miller,2011;Decreusefond et al.,2012;
Barbour and Reinert, 2013; Janson et al., 2013) or involve quantities that do
not depend on the temporal details of the disease dynamics (e.g. Newman,
2002; Kenah and Robins, 2007; Ball and Neal, 2008; Ball et al., 2009). In
thispaperweconsidergeneralnon-MarkovianSIRepidemicmodelsandfocus
our attention on arguably the most important piece of information concern-
ing the system’s temporal dynamics: the Malthusian parameter, or real-time
growth rate. This quantity corresponds to the rate of exponential growth in
the number of infectives observed in many models when an epidemic takes
off in a large population and the susceptible population is still large enough
that its reduction can be ignored.
Realistic patterns of contact between people typically involve repeated
interactions with the same individual, and significant heterogeneity in the
number of contacts reported (Danon et al., 2012). From an analytical point
of view, such a population structure is associated with three problems that
need to be addressed before the epidemic dynamics can be fully understood.
These problems are called repeated contacts, infection interval contraction,
and generational ordering. We begin by illustrating these problems on a
model scenario.
Consideraninfectionspreadingbetweenindividualswhoareabletotrans-
mit when they enter in contact with each other. Throughout the paper we
use the term contact to mean an infectious contact, i.e. a contact that is ‘suf-
ficiently intense’ to result in an infection whenever the individual that makes
the contact is infectious and the other is susceptible. A precise mathematical
description of the model is given in Section 2, but in short we assume that,
upon infection, individuals are attributed a (potentially time-varying) rate
at which they make infectious contacts. When normalised, this gives the in-
fectious contact interval distribution, i.e. the distribution of times between an
individual becoming infected and the subsequent infectious contacts made by
that individual (see also Svensson, 2007; Scalia Tomba et al., 2010; Kenah,
2011).
When an infectious contact occurs, the infective chooses an individual to
2
P P P
1 2 3
i j i j i j
k k k
Figure 1: Schematic representation of three populations, P ,P and P . In P mixing is
1 2 3 1
homogeneous and hence all individuals are connected to each other. In P and P , i, j
2 3
and k all have a small number of neighbours they can potentially infect, although in P
2
the network is tree-like, while in P the network presents a triangle between i,j and k.
3
be the ‘destination’ of the contact among a set of potential ones. Consider
now a population P where homogeneous mixing is assumed, so that any pair
1
of individuals can make infectious contacts, and compare it with populations
P and P where each individual can only make infectious contacts with
2 3
individuals in a fixed set of neighbours (see Figure 1). In all cases we suppose
that we are working with a large population of size N ≫ 1, meaning that we
will ignore effects that are, for example O(N−1). We represent populations
P and P with a static (undirected) network where edges connect pairs of
2 3
neighbouringindividuals. InSection2weregardtherateofmakinginfectious
contacts across each specified edge as our basic model ingredient. Note that
population P also admits a network representation, namely the complete
1
network where each node is connected to all others; however, due to our
assumption of a large population, the rate of making infectious contacts
across each specific edge vanishes and the neighbouring relation becomes
somewhat ephemeral.
We now focus our attention on a group of 3 individuals, labelled i, j and
k, where we assume that i is infectious and the other two are susceptible,
and we contrast the epidemic dynamics in P with those in P and P (see
1 2 3
Figure 1).
First of all, whilst in P individual i can potentially infect everybody in
1
thepopulation, inP andP , nomatter howinfectious iis, theycannot infect
2 3
more than four other individuals, thus ‘wasting’ part of their infectivity on
repeated contacts with the same individuals.
3
Secondly, in P each infectious contact is made with an individual chosen
1
at random and, because the population is large, it is unlikely for i to contact
the same individual more than once, so the average time at which trans-
mission (to j, say) occurs is equal to the average time at which i makes an
infectious contact with j. In P and P , where i may try to infect j multiple
2 3
times, the time of transmission corresponds to the time of the first infectious
contacts, which will occur on average earlier than the mean time at which
a randomly selected infectious contact from i to j occurs. We call this the
problem of infection interval contraction. In choosing this terminology we
have tried to avoid the very similar term ‘generation interval contraction’,
which is already adopted in Kenah et al. (2008) to denote a slightly different
phenomenon. Their approach based on survival analysis takes into account
automatically what we here call the infection interval contraction. What
they focus on instead is the fact that the time at which a susceptible is in-
fected is the minimum of the times of all first infectious contacts from all
potential infectors who are competing for infection of the same susceptible.
The assumptions we make below avoid this problem as they imply that there
is never more than one potential infector. In other words, the difference be-
tween infection interval contraction and generation interval contraction can
be explained by the two different comparisons being made. In Kenah et al.
(2008), different levels of infection prevalence are compared within a single
epidemic model; on the other hand, here we are comparing a network-based
epidemic(whererepeatedcontactsoccur)toamass-actionmodel(whereeach
infectious contact leads to an infection). In particular, because we are inter-
ested in characterising the real-time growth rate, we assume we are in the
early epidemic phase when the number of cases grows exponentially, i.e. in
our context the prevalence is assumed to be small.
Thirdly, once i has infected j in P and P , then j cannot infect i and
2 3
‘loses’ a neighbour due to generational ordering. Considering the next gen-
eration of infection, if j infects k (an individual with only one contact) then
k actually plays no further role in transmission events due to this effect. In
P , however, j’s infector is a negligible fraction of its total neighbours and so
1
this effect is insignificant.
What differentiates P from P is the presence of a triangle between i,j
3 2
and k. In the presence of short loops in the network (as it is the case, for
example, if a small number of individuals all live in the same household), in
addition to the previous three, two further effects become important.
Fourth, if i infects j in P , the probability that j infects k before i does is
1
4
negligible. Therefore the event that i infects k is not affected by the epidemic
in the neighbourhood of i or what i did before infecting j. The same occurs
in P , because of the lack of a transitive link between i and k. In P , instead,
2 3
j can infect k before i has the chance to, and so the number of susceptibles
around i can decrease because of infectious events not caused directly by i.
We refer to this as the effect of local saturation of susceptibles.
Fifth, if i infects j, j infects k and i subsequently tries to infect k (an
event which does not result in an infection as k is not susceptible any longer),
we have two different ways of attributing cases to different generations of
infection. Of course, assuming that i is in generation 0, the natural choice
wouldbetomaintainthereal-timedescription ofwhoinfectswhomandplace
j ingeneration1andk ingeneration2. However, thisapproachisanalytically
difficult to handle. The alternative is to consider all those that i would have
infected(bothj andk)andplacethembothingeneration1,inwhichcasethe
transmission event from j to k is then overlooked (Ludwig, 1975; Pellis et al.,
2008). We term this phenomenon overlapping generations. In Pellis et al.
(2012) and Ball et al. (2014) the impact of this problem in defining and
computing the basic reproduction number R in models that involve small
0
mixing groups (e.g. households or workplaces) is carefully examined.
In generic social structures, these five effects often occur simultaneously.
However, in this paper we consider scenarios in which they can be discrimi-
nated and progress can be made analytically, namely we focus on population
with the structure of P , where the problems of local saturation of suscepti-
2
bles and of overlapping generations described above need not be considered.
For this reason, we make the strong assumption that that the proportion of
possible transitive links is negligible. This means that the network is locally
tree-like. While some progress is possible for real-time growth rate calcula-
tions in the presence of many short loops in the network (like in P ), this
3
often involves either restriction to compartmental dynamics or approximate
methods (Fraser, 2007; Ross et al., 2010; Pellis et al., 2010). Analysis of gen-
eral stochastic dynamics of the kind we present here in the presence of short
loops would require a major conceptual advance not currently available, al-
though see Ball et al., 2014, for a deeper exploration of these problems in the
context of households models.
In addition to our first assumption of a locally tree-like network, we make
two further assumptions. The second assumption is that the individuals are
sufficiently ‘close’ to each other that the epidemic does indeed grow expo-
nentially during its early phase rather than more slowly (e.g. quadratically,
5
as would be expected on a two-dimensional lattice). More formally, this
assumption requires that mean path lengths on the giant component are
O(lnN). Thirdly, we assume that the distribution of numbers of contacts
is not too heterogeneous (technically, that the second moment of the degree
distribution does not diverge) so that the Malthusian parameter does not
diverge as the system size becomes very large.
The three assumptions detailed above are fundamental to our approach,
in the sense that it is not clear how to analyse real-time behaviour of an
epidemic on a network mathematically if they do not hold. We also make
other assumptions that significantly simplify the analytical results obtained,
but which can be relaxed with a certain amount of elementary but poten-
tially tedious algebra. In particular, throughout we assume individuals are
all identical to each other, in the sense that there is no variation in suscep-
tibility and in the network model they only differ in terms of their degree.
However, a key feature of the present work is a careful treatment of an indi-
vidual’s infectiousness, which is allowed to vary over time, according to some
stochastic infectivity profile. These infectivity profiles need not be the same
for every individual, but we assume that they are independent realisations
from a specified distribution that is the same for every individual. We also
assume that there are no degree-degree correlations, meaning that (given
our other technical assumptions) the configuration model can be used for the
contact structure.
In setting up our framework, we took inspiration from Diekmann et al.
(1998), where the authors focus on the computation of R in the presence of
0
repeated contacts under the same unclustered network approach discussed
above. However, we here extend their work by adding the computation of
the real-time growth rate and by moving from their deterministic framework
to a stochastic one. Any results involving the real-time growth rate has its
foundations on the so called Lotka-Euler equation, originated in the field of
demography anddiscussed, for example, inWallinga and Lipsitch (2007)and
references therein. In the simple case of a epidemic spreading in a large and
homogeneously mixing population, the Lotka-Euler equation is derived as
follows (see Diekmann et al., 2013, p. 10 and 212). Assume the population
has size N and each infective, t units after having been infected, makes
infectious contacts with other individuals at a rate B(t). Then the expected
6
incidence i(τ) at absolute time τ satisfies the renewal equation
∞
S(τ)
i(τ) = B(t)i(τ −t)dt, (1)
N Z
0
whereS(τ)isthenumberofsusceptibles attimeτ. Earlyon,whenalmostthe
entire population is susceptible the incidence grows exponentially, S(τ)/N ≈
1 and B(t) in this linearised process describes the rate at which new cases are
generated by a single infective at time t after their infection. Substituting
the Ansatz i(τ) = i erτ into (1), we deduce that the real-time growth rate r
0
must satisfy the Lotka-Euler equation:
∞
B(t)e−rt dt = 1. (2)
Z
0
Since the left-hand side of Equation (2) is continuous in r and has limits of
0 and ∞ (assuming B is not zero almost everywhere), a solution must exist
by the intermediate value theorem. It is also strictly monotonic in r, and so
any solution must also be unique.
While we believe the question of real-time behaviour of epidemics on
networks isofinherent theoretical interest, themainaimofthisworkisthree-
fold: (i) to make clear the key differences between the epidemic dynamics in
homogeneous mixing populations and on locally unstructured networks; (ii)
to provide a comprehensive list of results in analytically tractable cases, that
can be used by modellers to validate the outputs of complex simulations, for
exampleindividual-basedones; (iii)toprovidethetoolstoassesstheaccuracy
of approximations for the real time growth rate that involve either ignoring
the network structure altogether or neglecting part of the complexities to
achieve significantly simpler analytical results, as done in Fraser (2007) and
Pellis et al. (2010).
The rest of the paper is structured as follows. In Section 2 we present
our modelling approach, spelling out the assumptions behind the network
constructionandtheepidemicspreadingonit, andwedefinethefundamental
quantities of interest in the network model and their counterparts in the
limitinghomogeneously mixing scenario. Wethenconsider alistspecial cases
in Section 3, for which we provide extensive analytical results. In Sections
4 and 5 we discuss more general analytical considerations and fully general
7
Monte-Carlo methods for numerically computing the real-time growth rate.
Finally, extensive numerical illustrations are presented in Section 6 and final
considerations in Section 7.
2. Model definitions
2.1. Configuration network model
We consider a static network representing a population of individuals
and their interactions. We assume the population size N is large relative to
the number of infected individuals, which is in turn significantly larger than
unity, as we are interested in studying the asymptotic early spread of the
infection.
To ensure that the network is locally unclustered, we assume the network
isconstructedaccordingtotheso-calledconfiguration model (Molloy and Reed,
1995)withtheconventionsmostcommonlyusedinepidemiology(Ball and Neal,
2008; Ball et al., 2010). In this model, each individual i (i = 1,2,...,N)
is given a degree d which is an independent realisation of a non-negative
i
integer-valued random variable D. The distribution of D is called the degree
distribution. Individual i is then allocated d stubs (half-edges). If the total
i
number of stubs in the population is odd, the process of attributing a degree
to each individual is repeated. If the number of stubs is even, stubs are
selected at random without replacement and paired. The resulting network
might have some short loops and self-edges, but if a sequence of such net-
works is constructed and indexed with the population size N and we assume
that D has finite mean and variance, short loops and self-edges appear with
probability of order O(N−1) and can be ignored as N → ∞ (Durrett, 2007).
2.2. The infectivity profile
Thefunctionthatdescribesthewayanindividual’sinfectiousness changes
through time is called their infectivity profile. We make the simplifying as-
sumptions that individuals do not differ in terms of biological susceptibility
and that the infectivity profiles of different individuals are independent and
identically distributed according to some known distribution. Inaddition, we
assume that individuals can only transmit the infection across edges, that
each individual’s infectious behaviour applies equally across all links with
neighbours (in the precise way specified below) and that the infectious state
of the individuals does not modify or affect the network structure.
8
When infected, individual i is attributed a realisation λ (t) from the in-
i
fectivity profile distribution, where t ≥ 0 is the time since the infection of
i. For each neighbour in the network, i makes infectious contacts at the
points of an inhomogeneous Poisson process with rate λ (t). Thus, given
i
λ (t) for individual i, the probability that i makes infectious contact with a
i
given neighbour in the interval (t,t+∆t) is λ (t)∆t+o(∆t). If the individual
i
contacted by i is susceptible, then they become infected. We assume that
the Poisson processes governing the infectious contacts between every pair
of individuals are conditionally independent given the infectivity profile of
individual i.
We denote by Λ(t) (t ≥ 0) the random infectivity profile, which is the
same for every individual andof which the λ (t)s areindependent realisations
i
and we assume that E[ ∞Λ(t)dt] < ∞ and ∞E[Λ(t)]dt < ∞ and so by
0 0
Fubini’s Theorem theseRtwo integrals must beRequal.
2.3. The basic reproduction number R
0
We begin by calculating R , defined to be the expected number of sec-
0
ondary infections caused by a typical infective in the early stages of the epi-
∞
demic. Let therandom variableA = Λ(t)dt represent the total infectivity
0
spread by an arbitrary infective acroRss each of the links they have. Condi-
tional on degree D = d and total infectivity A = a, the number of secondary
cases that an infective will cause during the early stages of the epidemic is
distributed Bin(d−1,1−e−a). The number of trials d−1 takes into account
the fact that the infective must have acquired the infection from one of its
neighbours, who is therefore no longer susceptible (the generational ordering
effect mentioned in the introduction). All other neighbours are almost surely
susceptible because in the early stages almost all of the population is suscep-
tible and short loops in the network appear with negligible probability. The
probability of transmission occurring across an edge between the infective
and a susceptible neighbour (1 − e−a) takes into account that at least one
infectious contact across a link is necessary, but all subsequent contacts are
ineffective. Note that the infection of each of the susceptible neighbours are
not independent events, but are conditionally independent given A = a.
A typical infective during the early stages of the epidemic has more con-
nections than the average individual in the network. The probability that
transmission across an edge reaches an individual of degree d is equal to
d×P(D = d)/E[D], whichwecallthedegree-biased distribution. Conditional
9
on having total infectivity a, this individual of degree d then infects on av-
erage (d−1)(1−e−a) new cases. Since an individual’s degree D and their
infectivity A are independent, removing the conditioning gives
E[D(D−1)]
R = E 1−e−A . (3)
0 E[D]
(cid:2) (cid:3)
Hereafter we denote the mean degree of the node reached by a randomly
selected edge by
E[D(D −1)]
d∗ = . (4)
E[D]
2.4. The real-time growth rate
In this paper, our primary interest is studying the real-time growth rate
– the rate of the exponential growth in the number of infectives during the
early phase of an emerging epidemic – which is denoted by r. Unlike R , the
0
real-time growth rate depends on the full infectivity profile Λ(t) and not just
its integral A. Of all the infectious contacts towards a susceptible neighbour,
only the first contact results in an infection. The rate at which infections
occur is therefore equal to the rate at which contacts are made multiplied by
the probability that no contacts have occurred previously. Let Λ (t) denote
1
the rate at which first contacts occur, and hence
t
Λ1(t) = Λ(t)e−R0Λ(s)ds. (5)
Note that ∞Λ (t)dt = 1 − e−A. We can also describe the distribution
0 1
of times sinRce infection at which infectious contacts are made, which has
probability density function (pdf) w(t) = E[Λ(t)]/E[A] and the distribution
of times since infection at which infections are made, which has pdf w (t) =
1
E[Λ (t)]/E 1−e−A . Werefertotheformerastheinfectious contactinterval
1
distribution(cid:2)and the(cid:3) latter as the generation interval distribution. Let the
random variables W and W denote draws from these two distributions,
1
which are both implicitly conditional on the fact that at least one infectious
contact is made across the edge. The mean generation interval, often referred
to in the literature as generation time, of the epidemic on the network is
therefore
T = E[W ]. (6)
g 1
Using again the degree-biased distribution to compute the degree of a
randomly selected new infective, the average infection intensity of a typical
10