11footnotetext: School of Mathematical Sciences, The University of Nottingham,
University Park, Nottingham NG7 2RD, UK

# Central limit theorems for SIR epidemics and percolation on configuration
model random graphs

Frank Ball1

###### Abstract

We consider a stochastic SIR (susceptible $\to$ infective $\to$ recovered)
epidemic defined on a configuration model random graph, in which infective
individuals can infect only their neighbours in the graph during an infectious
period which has an arbitrary but specified distribution. Central limit
theorems for the final size (number of initial susceptibles that become
infected) of such an epidemic as the population size $n$ tends to infinity,
with explicit, easy to compute expressions for the asymptotic variance, are
proved assuming that the degrees are bounded. The results are obtained for
both the Molloy-Reed random graph, in which the degrees of individuals are
deterministic, and the Newman-Strogatz-Watts random graph, in which the
degrees are independent and identically distributed. The central limit
theorems cover the cases when the number of initial infectives either (a)
tends to infinity or (b) is held fixed as $n\to\infty$. In (a) it is assumed
that the fraction of the population that is initially infected converges to a
limit (which may be $0$) as $n\to\infty$, while in (b) the central limit
theorems are conditional upon the occurrence of a large outbreak (more
precisely one of size at least $\log n$). Central limit theorems for the size
of the largest cluster in bond percolation on Molloy-Reed and Newman-Strogatz-
Watts random graphs follow immediately from our results, as do central limit
theorems for the size of the giant component of those graphs. Corresponding
central limit theorems for site percolation on those graphs are also proved.

##  1 Introduction

There has been considerable work in the past two decades on models for the
spread of epidemics on random networks; see, for example, the recent book
[Kiss et al.(2017)]. The usual paradigm is that individuals in a population
are represented by nodes in a random graph and infected individuals are able
to transmit infection only to their neighbours in the graph. The graph is
often constructed using the configuration model (see, for example, [van der
Hofstad(2016)], Chapter 7), which allows for an arbitrary but specified degree
distribution. The most-studied type of epidemic model is the SIR (susceptible
$\to$ infective $\to$ recovered) model. In this model individuals are
classified into three types: susceptibles, infectives and recovered. If a
susceptible individual is contacted by an infective then it too becomes an
infective and remains so for a time, called its infectious period, that is
distributed according to a non-negative random variable $I$ having an
arbitrary but specified distribution. An infective individual recovers at the
end of its infectious period and is then immune to further infection. During
its infectious period, an infective contacts its susceptible neighbours in the
graph independently at the points of Poisson processes having rate $\lambda$.
The graph is assumed to be static and the population closed (i.e. there are no
births or deaths), so eventually the epidemic process terminates. The final
size of the epidemic is the number of initial susceptibles that are infected
during its course. The final size is a key epidemic statistic, not only as a
measure of the impact of an epidemic but also in an inferential setting, since
often it can be observed more reliably than the precise temporal spread. The
main aim of this paper is to develop central limit theorems for the final size
of an SIR epidemic on configuration model graphs as the population size
$n\to\infty$.

The configuration model, which was introduced by [Bollobás(1980)], is a model
for random graphs with a given degree sequence. There are two distinct
approaches for constructing configuration model graphs with a given degree
distribution as $n\to\infty$. In both approaches, individuals are assigned a
number of half-edges, corresponding to their degree, and then these half-edges
are paired uniformly at random. In [Molloy and Reed(1995)], the degrees of
indiviudals are prescribed deterministically whilst in [Newman et al.(2001)]
they are i.i.d. (independent and identically distributed) copies of a random
variable $D$, that describes the limiting degree distribution. We refer to the
former as the MR random graph and to the latter as the NSW random graph.
Subject to suitable conditions on the degree sequences in the MR model and $D$
in the NSW model, law of large number limits for SIR epidemics on the two
graphs are the same. That is not the case for central limit theorems, as for
finite $n$, there is greater variability in the degrees of individuals in the
NSW model than in the MR model; indeed, in the NSW model, such variability is
of the same order of magnitude as the variability in the epidemic. Thus,
though the asymptotic means are the same, the asymptotic variances in the
central limit theorems for final size are greater for the epidemic on the NSW
random graph.

There have been numerous studies, some fully rigorous and some heuristic, of
SIR epidemics on configuration model networks making various assumptions
concerning the infectious period random variable $I$. For example, assuming
$I$ is constant, [Andersson(1998)] derives a law of large numbers for the
final size of an epidemic on an NSW random graph when a strictly positive
fraction of the population is initially infected in the limit as $n\to\infty$
and [Britton et al.(2007)] obtain a similar result for epidemics on an MR
random graph initiated by a single infective. In the latter case, a large
outbreak is possible only if the basic reproduction number $R_{0}>1$; see
(2.8) in Section 2.4. In a highly influential paper, [Newman(2002)] uses
heuristic percolation arguments to obtain a number of results, including the
fraction of the population infected by a large outbreak, for SIR epidemics on
NSW random graphs with $I$ having an arbitrary but specified distribution.
Several authors have studied the case when $I$ has an exponential
distribution, so the model becomes Markovian. [Decreusefond et al.(2012)]
obtain a law of large numbers type result for the epidemic process on an NSW
random graph with a strictly positive fraction initially infected, which
yields a rigorous proof of the deterministic approximation of [Volz(2008)]
(see also [Miller(2011)] and [Miller et al.(2012)]). [Bohman and
Picollelli(2012)] obtain law of large numbers results for both the process and
final size of an epidemic with one intial infective on an MR random graph with
bounded degrees. [Janson et al.(2012)] obtain similar results under weaker
conditions on the degree sequences considering the cases when the fraction
initially infected, in the limit as $n\to\infty$, is either strictly positive
or zero (assuming of course there is at least one initial infective). In the
latter case, the limiting “deterministic” process involves a random time
translation reflecting the time taken for the number of infectives to reach
order $n$; a similar result is obtained by [Barbour and Reinert(2013)]
assuming a bounded degree sequence and an arbitrary but specified distribution
for $I$.

There has been very little work to date on central limit theorems for SIR
epidemics on configuration model networks. A functional central limit theorem
for the SI epidemic (in which ${\rm P}(I=\infty)=1$, so infectives remain
infectious forever) on an MR random graph with unbounded degrees is obtained
by [KhudaBukhsh et al.(2017)], who note that their method is not
straightforward to extend to an SIR model. Assuming that $I$ follows an
exponential distribution and bounded degrees, [Ball et al.(2018)] use an
effective degree approach ([Ball and Neal(2008)]) and density dependent
population processes ([Ethier and Kurtz(1986)], Chapter 11) to obtain
functional central limit theorems for SIR epidemics on MR and NSW random
graphs, in which susceptible individuals can also drop their edges to
infective neighbours. They also conjecture central limit theorems for the
final size of such epidemics (and hence as a special case for the final size
of standard SIR epidemics, without dropping of edges), assuming either a
strictly positive fraction or a constant number of initial infectives (in
which case the central limit theorem is conditional on the occurrence of a
large outbreak). However, the arguments are not fully rigorous and the result
for a constant number of initial infectives is based purely on the existence
of equivalent results for other (non-network) SIR epidemic models. Another
limitation of [Ball et al.(2018)] is the assumption that $I$ is exponentially
distributed, which is unrealistic for most real-life diseases.

In the present paper, we address these shortcomings and derive fully rigorous
central limit theorems for the final size of SIR epidemics on MR and NSW
random graphs having bounded degrees, when the infectious period $I$ follows
an arbitrary but specified distribution. We consider the cases when the
limiting fraction of the population is (i) strictly positive and (ii) zero.
For the latter we treat the situtations where the number of initial infectives
either (i) is held fixed independent of $n$ or (ii) tends to $\infty$ as
$n\to\infty$. The mean parameter $\rho$ in the central limit theorems, which
coincides with the corresponding law of large numbers limit, depends on the
solution $z$ of a non-linear equation (see (2.2) and (2.9) in Section 2.4).
Given $z$, the variance parameter in the central limit theorems is fully
explicit and hence easy to compute.

If $I$ is constant, say ${\rm P}(I=1)=1$ then the above SIR model is
essentially bond percolation with probability $\pi=1-{\rm e}^{-\lambda}$ and
if ${\rm P}(I=\infty)=\pi=1-{\rm P}(I=\infty)=0$ then it is closely related to
site percolation with probability $\pi$; see, for example, [Durrett(2007)],
page 15, and [Janson(2009a)]. Central limit theorems for the size of the giant
component (largest cluster) of bond percolation on the MR and NSW random
graphs follow immediately from our results. (Corresponding theorems for site
percolation are also obtained using our methodology.) Further, setting $\pi=1$
yields central limit theorems for the giant component of those graphs (Remark
2.6); cf. [Barbour and Röllin(2017)] who obtain a central limit theorem for
the giant component of the MR random graph and [Ball and Neal(2017)] who
derive the asymptotic variance of the giant component of MR and NSW random
graphs, all allowing for unbounded degrees.

The proofs involve constructing the random graph and epidemic on it
simultanaeously, modifying the infection mechanism so that when a susceptible
is infected it decides which of its half-edges it will try to infect along
with (its remaining half-edges becoming recovered half-edges), with the times
of those infection attempts (relative to the time of infection of the
susceptible) being realisations of i.i.d. exponential random variables. The
distribution of the final outcome of the epidemic, and hence also its final
size, is invariant to this modification. The process describing the evolution
of the numbers of susceptibles of different degrees, infective half-edges and
recovered half-edges is an asymptotically density dependent population process
([Ethier and Kurtz(1986)], Chapter 11, and [Pollett(1990)]). The asymptotic
distribution of the final outcome of the epidemic is studied by considering a
boundary crossing problem for a random time-scale transformation of that
process. The proofs extend, at least in principle, to SIR epidemics and
percolation on extensions of the configuration model that include fully-
connected cliques ([Trapman(2007)], [Gleeson(2009)], [Ball et al.(2010)] and
[Coupechoux and Lelarge(2014)]), though explicit calculation of the asymptotic
variances may be difficult.

The remainder of the paper is organised as follows. The MR and NSW random
graphs are defined in Section 2.1 and the SIR epidemic model is described in
Section 2.2. The main central limit theorems (Theorems 2.1-2.3 for SIR
epidemics and Theorem 2.7 for percolation) are stated in Section 2.4, together
with some remarks giving comparisons of their variance parameters and
applications to giant components indicated above. Some numerical
illustrations, which show that the central limit theorems can yield good
approximations even for relatively small graphs, are given in Section 3. The
proofs are given in Section 5. They make extensive use of asymptotically
density dependent population processes and in particular require a version of
the functional central limit theorem for such processes to include
asymptotically random initial conditions. For ease of reference, the required
results for such processes are collected together in Section 4. Some brief
concluding comments are given in Section 6. Calculation of the asymptotic
variances for the central limit theorems is lengthy, though straightforward,
so this and a few other details are deferred to an appendix.

###  1.1 Notation

All vectors are row vectors and ⊤ denotes transpose. With the dimension being
obvious from the context, $I$ denotes an identity matrix and $\boldsymbol{0}$
and $\boldsymbol{1}$ denotes vectors all of whose elements are $0$ and $1$,
respectively. For $x\in\mathbb{R}$, the usual floor and ceiling functions are
denoted by $\left\lfloor{x}\right\rfloor$ and $\left\lceil{x}\right\rceil$,
respectively. Thus $\left\lfloor{x}\right\rfloor$ is the greatest integer
$\leq x$ and $\left\lceil{x}\right\rceil$ is the smallest integer $\geq x$.
For a positive integer $k$, the $k$th derivative of a real-valued function $f$
is denoted by $f^{(k)}$. The cardinality of a set $A$ is denoted by $|A|$.
Sums are zero if vaccuous. We use $\stackrel{{\scriptstyle
p}}{{\longrightarrow}}$, $\stackrel{{\scriptstyle a.s.}}{{\longrightarrow}}$
and $\stackrel{{\scriptstyle D}}{{\longrightarrow}}$ to denote convergence in
probability, convergence almost sure and convergence in distribution,
respectively.

Further, ${\rm U}(0,1)$ denotes a uniform random variable on $(0,1)$; ${\rm
Exp}(1)$ denotes an exponential random variable with mean $1$; ${\rm
N}(0,\sigma^{2})$ denotes a univariate normal random variable with mean $0$
and variance $\sigma^{2}$; and ${\rm N}(\boldsymbol{0},\Sigma)$ denotes a
multivariate normal random variable with mean vector $\boldsymbol{0}$ and
variance matrix $\Sigma$, whose dimension again is obvious from the context.
For a positive integer $n$ and $p\in[0,1]$, ${\rm Bin}(n,p)$ denotes a
binomial random variable with $n$ trials and success probability $p$. Also, if
$I$ is a non-negative random variable and $\lambda\in(0,\infty)$ then ${\rm
Bin}(n,1-{\rm e}^{-\lambda I})$ denotes a mixed-Binomial random variable
obtained by first sampling $I_{1}$ from the distribution of $I$ and then,
given $I_{1}$, sampling independently from ${\rm Bin}(n,1-{\rm e}^{-\lambda
I_{1}})$. Similarly, if $I$ is a non-negative random variable and $D$ is a
non-negative integer-valued random variable, then ${\rm Bin}(D,1-{\rm
e}^{-\lambda I})$ denotes a mixed-Binomial random variable, where the
realisations of $D$ and $I$ are independent. Thus, if $X\sim{\rm Bin}(D,1-{\rm
e}^{-\lambda I})$, then

${\rm P}(X=k)=\sum_{d=k}^{\infty}{\rm P}(D=d){\rm E}\left[\binom{d}{k}(1-{\rm
e}^{-\lambda I})^{k}{\rm e}^{-(d-k)\lambda I}\right]\qquad(k=0,1,\dots).$

Note that we allow the possibility $D=0$.

##  2 Model and main results

###  2.1 Random graph

Consider a population of $n$ indivdiuals labelled $1,2,\dots,n$. For
$i=1,2,\dots,n$, let $D^{(n)}_{i}$ denote the degree of individual $i$. We
assume that $0\leq D^{(n)}_{i}\leq d_{\rm max}$ for all $i$, i.e. that there
is a maximum degree $d_{\rm max}$. In the MR random graph the degrees are
prescribed, while in the NSW random graph
$D^{(n)}_{1},D^{(n)}_{2},\dots,D^{(n)}_{n}$ are i.i.d. copies of a random
variable $D$ having probability mass function given by ${\rm P}(D=i)=p_{i}$
$(i=0,1,\dots,d_{\rm max})$. Under both models, the network (random graph) is
formed by attaching $D^{(n)}_{i}$ half-edges to individual $i$, for
$i=1,2,\dots,n$, and then pairing up the
$D^{(n)}_{1}+D^{(n)}_{2}+\dots+D^{(n)}_{n}$ half-edges uniformly at random to
give the edges in the random graph, which we denote by $\mathcal{G}^{(n)}$. In
the NSW model, if $D^{(n)}_{1}+D^{(n)}_{2}+\dots+D^{(n)}_{n}$ is odd there is
a left-over stub, which is ignored. (Of course in the MR model the prescribed
degrees can be chosen so that $D^{(n)}_{1}+D^{(n)}_{2}+\dots+D^{(n)}_{n}$ is
even.)

We are interested in asymptotic results as the number of individuals
$n\to\infty$. In the MR random graph, for $i=0,1,\dots,d_{\rm max}$, let
$v^{(n)}_{i}=\sum_{k=1}^{n}1_{\\{D^{(n)}_{k}=i\\}}$ be the number of
individuals having degree $i$. We assume that

$\lim_{n\to\infty}\sqrt{n}\left(n^{-1}v^{(n)}_{i}-p_{i}\right)=0\qquad(i=0,1,\dots,d_{\rm
max}).$ (2.1)

Note that (2.1) implies $\lim_{n\to\infty}n^{-1}v^{(n)}_{i}=p_{i}$
$(i=0,1,\dots,d_{\rm max})$.

In both models the random graph may have some imperfections, specifically
self-loops and multiple edges, but they are sparse in the network as
$n\to\infty$; more precisely, the number of such imperfections converges in
distrubution to a Poisson random variable as $n\to\infty$ ([Durrett(2007)],
Theorem 3.1.2). Moreover, [Janson(2009b)] implies that the probability that
the random graph is simple is bounded away from $0$, as $n\to\infty$, and law
of large numbers results continue to hold if the graph is conditioned on being
simple. However, that is not necessarily the case for convergence in
distribution; see [Janson(2010)], Remark 1.4, and [Barbour and Röllin(2017)],
Remark 2.5. Thus whether or not our central limit theorems continue to hold
when the random graph is conditioned on being simple is an open question.

###  2.2 SIR epidemic

An SIR epidemic, denoted by $\mathcal{E}^{(n)}$, is constructed on the above
network as follows. Initially, at time $t=0$, a number of individuals are
infective and the remaining individuals are susceptible. (Precise statements
concerning the initial infectives are made later.) Distinct infectives behave
independently of each other and of the construction of the network. Each
infective remains infectious for a period of time that is distributed
according to a random variable $I$, having an arbitrary but specified
distribution, after which it becomes recovered. During its infectious period,
an infective contacts its neighbours in the network independently at the
points of Poisson processes each having rate $\lambda$, so the probability
that a given neighbour is contacted is $p_{I}=1-\phi(\lambda)$, where
$\phi(\theta)=\rm E[\exp(-\theta I)]$ $(\theta\geq 0)$ is the Laplace
transform of $I$. If a contacted individual is susceptible then it becomes an
infective, otherwise nothing happens. The epidemic ends when there is no
infective individual in the population.

For $t\geq 0$, let $X^{(n)}_{i}(t)$ be the number of degree-$i$ susceptible
individuals at time $t$ ($i=0,1,\dots,d_{\rm max}$) and let $Y^{(n)}(t)$ be
the total number of infectives at time $t$. Let $\tau^{(n)}=\inf\\{t\geq
0:Y^{(n)}(t)=0\\}$ be the time of the end of the epidemic. Then
$T^{(n)}_{i}=X^{(n)}_{i}(0)-X^{(n)}_{i}(\tau^{(n)})$ is the total number of
degree-$i$ susceptibles that are infected by the epidemic. Let
$T^{(n)}=\sum_{i=0}^{d_{\rm max}}T^{(n)}_{i}$ be the total number of
susceptibles infected by the epidemic, i.e. the final size of the epdiemic. We
are primarily interested in the asymptotic distribution of $T^{(n)}$ as
$n\to\infty$.

The proofs allow for the possibility that $p_{I}=1$, i.e. ${\rm
P}(I=\infty)=1$. In that case the set of individuals that are infected during
the epidemic $\mathcal{E}^{(n)}$ comprises all individuals in the components
of $\mathcal{G}^{(n)}$ that contain at least one initial infective. Thus
central limit theorems for the size of the giant component in MR and NSW
random graphs follow immediately from our results; see Remark 2.6 below.

###  2.3 Bond and site percolation

In bond percolation on $\mathcal{G}^{(n)}$, each edge in $\mathcal{G}^{(n)}$
is deleted indpendently with probability $1-\pi$, while in site percolation
$\mathcal{G}^{(n)}$, each vertex (together with all incident edges) is deleted
independently with probability $1-\pi$ ([Janson(2009a)]). Interest is often
focused on the size, $C^{(n)}$ say, of the largest connected component in the
resulting graph.

###  2.4 Main results

For $i=0,1,\dots,d_{\rm max}$, let $a^{(n)}_{i}$ be the number of degree-$i$
initial infectives in the epidemic $\mathcal{E}^{(n)}$ and let
$a^{(n)}=\sum_{i=0}^{d_{\rm max}}a^{(n)}_{i}$ denote the total number of
initial infectives. In the epidemic on the MR random graph, we assume that
$a^{(n)}_{0},a^{(n)}_{1},\dots,a^{(n)}_{d_{\rm max}}$ are prescribed. In the
epidemic on the NSW random graph, we assume that $a^{(n)}$ is prescribed and
that the $a^{(n)}$ initial infectives are chosen by sampling uniformly at
random without replacement from the $n$ individuals in the population.

Let $\epsilon^{(n)}=n^{-1}a^{(n)}$ and $\epsilon^{(n)}_{i}=n^{-1}a^{(n)}_{i}$
$(i=0,1,\dots,d_{\rm max})$. Suppose that $\epsilon^{(n)}\to\epsilon$ as
$n\to\infty$ and that $\lim_{n\to\infty}\sqrt{n}(\epsilon^{(n)}-\epsilon)=0$.
For the epidemic on the MR random graph, suppose further that, for
$i=0,1,\dots,d_{\rm max}$, there exists $\epsilon_{i}$ such that
$\lim_{n\to\infty}\sqrt{n}(\epsilon^{(n)}_{i}-\epsilon_{i})=0$. For the
epidemic on the NSW random graph, let $\epsilon_{i}=\epsilon p_{i}$
($i=0,1,\dots,d_{\rm max}$). Let $\mu_{D}={\rm E}[D]$, $\sigma_{D}^{2}=\rm
var(D)$ and, for $s\in[0,1]$, $f_{D}(s)=\sum_{i=0}^{d_{\rm max}}p_{i}s^{i}$
and $f_{D_{\epsilon}}(s)=\sum_{i=0}^{d_{\rm max}}(p_{i}-\epsilon_{i})s^{i}$.
Let $q_{I}=1-p_{I}=\phi(\lambda)$ be the probability that an infective fails
to contact a given neighbour and $q_{I}^{(2)}=\phi(2\lambda)$ be the
probability that a given infective fails to contact two given neighbours.

The first theorem concerns the case $\epsilon>0$, so in the limit as
$n\to\infty$ a strictly positive fraction of the population is initially
infective. Let $T^{(n)}_{\rm MR}$ and $T^{(n)}_{\rm NSW}$ denote the final
size of the epidemic $\mathcal{E}^{(n)}$ on the MR and NSW random graphs,
respectively. Let $z$ be the unique solution in $[0,1)$ of

$z-q_{I}=\mu_{D}^{-1}p_{I}f_{D_{\epsilon}}^{(1)}(z)$ (2.2)

and

$\rho=1-\epsilon-f_{D_{\epsilon}}(z).$ (2.3)

###### Theorem 2.1

Suppose that $p_{i}\epsilon_{i}>0$ for at least one $i>0$ and, if $p_{I}=1$
then $p_{1}-\epsilon_{1}>0$. Then, as $n\to\infty$,

$\sqrt{n}\left(n^{-1}T^{(n)}_{\rm MR}-\rho\right)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}{\rm N}(0,\sigma^{2}_{\rm MR})$

and

$\sqrt{n}\left(n^{-1}T^{(n)}_{\rm NSW}-\rho\right)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}{\rm N}(0,\sigma^{2}_{\rm NSW}),$

where

$\displaystyle\sigma^{2}_{\rm MR}$
$\displaystyle=h(z)^{2}\left\\{\left[(p_{I}q_{I}+2(z-q_{I})^{2}\right]\mu_{D}-p_{I}^{2}\left[f_{D_{\epsilon}}^{(1)}(z^{2})+z^{2}f_{D_{\epsilon}}^{(2)}(z^{2})\right]\right\\}$
(2.4)
$\displaystyle\,+h(z)\left[2p_{I}zf_{D_{\epsilon}}^{(1)}(z^{2})-(z-q_{I})\mu_{D}\right]+1-\epsilon-\rho-
f_{D_{\epsilon}}(z^{2})$
$\displaystyle\,+\left(q_{I}^{(2)}-q_{I}^{2}\right)h(z)^{2}\left[f_{D}^{(2)}(1)-f_{D_{\epsilon}}^{(2)}(z)\right],$

with

$h(z)=\frac{p_{I}^{-1}(q_{I}-z)}{1-p_{I}\mu_{D}^{-1}f_{D_{\epsilon}}^{(2)}(z)},$
(2.5)

and

$\displaystyle\sigma^{2}_{\rm NSW}$
$\displaystyle=\frac{\rho(1-\epsilon-\rho)}{1-\epsilon}-h(z)(z-q_{I})\left(1-2\epsilon+\frac{2\epsilon\rho}{1-\epsilon}\right)\mu_{D}$
(2.6)
$\displaystyle\,+h(z)^{2}\left\\{\left[p_{I}q_{I}-2z(z-q_{I})\right]\mu_{D}+(z-q_{I})^{2}\left(\sigma_{D}^{2}+\frac{1-2\epsilon}{1-\epsilon}\mu_{D}^{2}\right)\right\\}$
$\displaystyle\,+\left(q_{I}^{(2)}-q_{I}^{2}\right)h(z)^{2}\left[f_{D}^{(2)}(1)-(1-\epsilon)f_{D}^{(2)}(z)\right],$

with

$h(z)=\frac{p_{I}^{-1}(q_{I}-z)}{1-p_{I}\mu_{D}^{-1}(1-\epsilon)f_{D}^{(2)}(z)}.$
(2.7)

The second theorem concerns the case when the number of initial infectives is
held fixed as $n\to\infty$, so $\epsilon=0$. More specifically, in the
epidemic on the MR random graph, we assume that $a^{(n)}_{i}=a_{i}$
$(i=0,1,\dots,d_{\rm max})$ for all $n\geq a=\sum_{i=1}^{d_{\rm max}}a_{i}$
and in the epidemic on the NSW random graph, we assume that $a^{(n)}=a$ for
all $n\geq a$. It is well known that, for large $n$, the process of infectives
in the early stages of such an epidemic can be approximated by a Galton-Watson
branching process, $\mathcal{B}$ say, in which, except for the initial
generation, the offspring distribution is ${\rm Bin}(\tilde{D}-1,1-{\rm
e}^{-\lambda I})$, where $\tilde{D}$ and $I$ are independent and $\tilde{D}$
has the size-biased degree distribution ${\rm
P}(\tilde{D}=k)=\mu_{D}^{-1}kp_{k}$ $(k=1,2,\dots,d_{\rm max})$; see, for
example, [Ball and Sirl(2013)]. This offspring distribution has mean

$R_{0}=\mu_{\tilde{D}-1}p_{I}=\left(\mu_{D}+\mu_{D}^{-1}\sigma_{D}^{2}-1\right)p_{I},$
(2.8)

where $\mu_{\tilde{D}-1}={\rm E}[\tilde{D}-1]$. The quantity $R_{0}$ is called
the basic reproduction number of the epidemic.

For $\mathcal{E}^{(n)}$, we say that a major outbreak occurs if and only if
the event $G^{(n)}=\\{T^{(n)}\geq\log n\\}$ occurs. Now $\lim_{n\to\infty}{\rm
P}(G^{(n)})={\rm P}(B=\infty)$, where $B$ is the total progeny (not including
the initial generation) of the branching process $\mathcal{B}$ (cf. Theorem
5.1). Thus, in the limit as $n\to\infty$, a major outbreak occurs with non-
zero probability if and only if $R_{0}>1$.

Suppose that that $R_{0}>1$. Now let $z$ be the unique solution in $[0,1)$ of

$z-q_{I}=\mu_{D}^{-1}p_{I}f_{D}^{(1)}(z)$ (2.9)

and $\rho=1-f_{D}(z)$. Note that if $p_{I}\in(0,1)$, or $p_{I}=1$ and
$p_{1}>0$, then $z>0$.

###### Theorem 2.2

Suppose that $R_{0}>1$ and, if $p_{I}=1$ then $p_{1}>0$. Then, as
$n\to\infty$,

$\sqrt{n}\left(n^{-1}T^{(n)}_{\rm
MR}-\rho\right)|G^{(n)}\stackrel{{\scriptstyle D}}{{\longrightarrow}}{\rm
N}(0,\tilde{\sigma}^{2}_{\rm MR})$

and

$\sqrt{n}\left(n^{-1}T^{(n)}_{\rm
NSW}-\rho\right)|G^{(n)}\stackrel{{\scriptstyle D}}{{\longrightarrow}}{\rm
N}(0,\tilde{\sigma}^{2}_{\rm NSW}),$

where $\tilde{\sigma}^{2}_{\rm MR}$ is given by (2.4) and (2.5), with
$f_{D_{\epsilon}}$ replaced by $f_{D}$, and $\tilde{\sigma}^{2}_{\rm NSW}$ is
obtained by setting $\epsilon=0$ in (2.6) and (2.7).

The next theorem concerns the case when $\epsilon=0$ but the number of initial
infectives $a^{(n)}\to\infty$ as $n\to\infty$.

###### Theorem 2.3

Suppose that $R_{0}>1$, $\epsilon=0$, $\sum_{i=1}^{d_{\rm
max}}a^{(n)}_{i}\to\infty$ as $n\to\infty$ and, if $p_{I}=1$ then $p_{1}>0$.
Then, as $n\to\infty$,

$\sqrt{n}\left(n^{-1}T^{(n)}_{\rm MR}-\rho\right)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}{\rm N}(0,\tilde{\sigma}^{2}_{\rm MR})$

and

$\sqrt{n}\left(n^{-1}T^{(n)}_{\rm NSW}-\rho\right)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}{\rm N}(0,\tilde{\sigma}^{2}_{\rm NSW}),$

where $\tilde{\sigma}^{2}_{\rm MR}$ and $\tilde{\sigma}^{2}_{\rm NSW}$ are as
in Theorem 2.2.

###### Remark 2.4

Note that $q_{I}^{(2)}=q_{I}^{2}$ if $I$ is almost surely constant, otherwise
$q_{I}^{(2)}>q_{I}^{2}$ by Jensen’s inequality. Also
$f_{D}^{(2)}(1)-f_{D_{\epsilon}}^{(2)}(z)>0$, so as one would expect on
intuitive grounds, if $p_{I}$ is held fixed, the asymptotic variance
$\sigma^{2}_{\rm MR}$ is smallest when the infectious period is constant. A
similar comment holds for $\sigma^{2}_{\rm NSW},\tilde{\sigma}^{2}_{\rm MR}$
and $\tilde{\sigma}^{2}_{\rm NSW}$.

###### Remark 2.5

Although it is not transparent from (2.4) and (2.6), it is seen easily from
the proof that, again as one would expect on intuitive grounds,
$\sigma^{2}_{\rm NSW}\geq\sigma^{2}_{\rm MR}$ and $\tilde{\sigma}^{2}_{\rm
NSW}\geq\tilde{\sigma}^{2}_{\rm MR}$, with strict inequalities unless the
support of the degree random variable $D$ is concentrated on a single point
(when the two models are identical); see Appendix A.3.4.

###### Remark 2.6

In the setting of Theorem 2.2, if $a^{(n)}=1$ and $p_{I}=1$ then with
probability tending to $1$ as $n\to\infty$, the event $G^{(n)}$ occurs if and
only if the initial infective belongs to the giant component of
$\mathcal{G}^{(n)}$. Thus setting $p_{I}=1$ in Theorem 2.2 yields central
limit theorems for giant components of MR and NSW random graphs.

The final theorem is concerned with percolation. Let $R_{0}$ be given by (2.8)
with $p_{I}=\pi$. Let $C^{(n)}_{\rm MR,B}$ and $C^{(n)}_{\rm MR,S}$ denote the
size of the largest connected component after bond and site percolation,
respectively, on an MR graph; define $C^{(n)}_{\rm NSW,B}$ and $C^{(n)}_{\rm
NSW,S}$ analogously for the NSW graph. Then, if $R_{0}>1$, for each of these
four choices for $C^{(n)}$, there exists $\epsilon>0$ such that
$\lim_{n\to\infty}{\rm P}(C^{(n)}\geq\epsilon n)=1$; see [Janson(2009a)],
Theorems 3.5 and 3.9, which also give law of large number limits for the
$C^{(n)}$ under weaker conditions than here. The following theorem gives
associated central limit theorems.

###### Theorem 2.7

Suppose that $\pi\in(0,1)$ and $R_{0}>1$. Let $z$ be the unique solutiuon in
$(0,1)$ of

$(z-1+\pi)=\mu_{D}^{-1}\pi f_{D}^{(1)}(z)$ (2.10)

and $\rho=1-f_{D}(z)$. Then, as $n\to\infty$,

$\sqrt{n}\left(n^{-1}C^{(n)}_{\rm MR,B}-\rho\right)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}{\rm N}(0,\sigma^{2}_{\rm MR,B}),$ (2.11)
$\sqrt{n}\left(n^{-1}C^{(n)}_{\rm NSW,B}-\rho\right)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}{\rm N}(0,\sigma^{2}_{\rm NSW,B}),$ (2.12)
$\sqrt{n}\left(n^{-1}C^{(n)}_{\rm MR,S}-\pi\rho\right)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}{\rm N}(0,\sigma^{2}_{\rm MR,S}),$ (2.13)

and

$\sqrt{n}\left(n^{-1}C^{(n)}_{\rm NSW,S}-\pi\rho\right)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}{\rm N}(0,\sigma^{2}_{\rm NSW,S}),$ (2.14)

where

$\displaystyle\sigma^{2}_{\rm MR,B}$
$\displaystyle=h(z)^{2}\left\\{\left[\pi(1-\pi)+2(z-1+\pi)^{2}\right]\mu_{D}-\pi^{2}\left[f_{D}^{(1)}(z^{2})+z^{2}f_{D}^{(2)}(z^{2})\right]\right\\}$
(2.15) $\displaystyle\,+h(z)\left[2\pi
zf_{D}^{(1)}(z^{2})-(z-1+\pi)\mu_{D}\right]+1-\rho-f_{D}(z^{2})$

and

$\displaystyle\sigma^{2}_{\rm NSW,B}$
$\displaystyle=\rho(1-\rho)-h(z)(z-1+\pi)\mu_{D}$ (2.16)
$\displaystyle\,+h(z)^{2}\left\\{\left[\pi(1-\pi)-2z(z-1+\pi)\right]\mu_{D}+(z-1+\pi)^{2}\left(\sigma_{D}^{2}+\mu_{D}^{2}\right)\right\\},$

with

$h(z)=\frac{\pi^{-1}(1-\pi-z)}{1-\pi\mu_{D}^{-1}f_{D}^{(2)}(z)};$
$\displaystyle\sigma^{2}_{\rm MR,S}$
$\displaystyle=\pi^{2}\left\\{\sigma^{2}_{\rm
MR,B}+\pi(1-\pi)h(z)^{2}\left[f_{D}^{(2)}(1)-f_{D}^{(2)}(z)\right]\right\\}$
(2.17) $\displaystyle\qquad+\pi(1-\pi)[\rho-2h(z)(1-z)\mu_{D}]$ (2.18)

and $\sigma^{2}_{\rm NSW,S}$ is given by (2.17), with $\sigma^{2}_{\rm MR,B}$
replaced by $\sigma^{2}_{\rm NSW,B}$.

##  3 Illustrations

In this section, we illustrate the central limit theorems in Section 2.4 by
using simulations to explore their applicability to graphs with finite $n$. We
also investigate briefly the dependence of the limiting variances in the
central limit theorems on the degree distribution, graph type and infectious
period distribution. We consider four degree distributions:

  1. (i)

$D\equiv d$, i.e. $D$ is constant with $p_{d}=1$;

  2. (ii)

$D\sim{\rm Po}(\mu_{D})$, i.e. $D$ is Poisson with mean $\mu_{D}$;

  3. (iii)

$D\sim{\rm Geom}(p)$, i.e. $p_{k}=(1-p)^{k}p$ $(k=0,1,\dots)$;

  4. (iv)

$D\sim{\rm Power}(\alpha,\kappa)$, i.e. $p_{k}=ck^{-\alpha}{\rm
e}^{-\frac{k}{\kappa}}$ $(k=1,2,\dots)$, where $\alpha,\kappa\in(0,\infty)$
and the normalising constant $c={\rm Li}_{\alpha}({\rm
e}^{-\frac{1}{\kappa}})$, with ${\rm Li}_{\alpha}$ being the polylogarithm
function.

The fourth distribution is a power law with exponential cut-off (see, for
example, [Newman(2002)]) that has been used extensively in the physics
literature. Note that, with $\theta={\rm e}^{-\frac{1}{\kappa}}$ and
$\beta={\rm Li}_{\alpha}(\theta)$, $f_{D}(s)=\beta^{-1}{\rm
Li}_{\alpha}(\theta s)$, $f_{D}^{(1)}(s)=\frac{1}{\beta s}{\rm
Li}_{\alpha-1}(\theta s)$ and $f_{D}^{(2)}(s)=\frac{1}{\beta s^{2}}[{\rm
Li}_{\alpha-2}(\theta s)-{\rm Li}_{\alpha-1}(\theta s)]$, which enables
$R_{0}$ and the asymptotic means and variances in the central limit theorems
to be calculated. Distributions (ii)-(iv) have unbounded support, so do not
satisfy the conditions of our theorems. The asymptotic ditributions in this
section are calculated under the assumption that the theorems still apply.

Each simulation consists of first simulating one graph, by simulating
$D^{(n)}_{1},D^{(n)}_{2},\dots,D^{(n)}_{n}$ and pairing up the half-edges, and
then simulating one epidemic, or percolation process, on it. For an MR graph
with (limiting) degree random variable $D$ on $n$ nodes, the degrees are given
by $D^{(n)}_{i}=\inf\\{d\in\mathbb{Z}_{+}:F_{D}(d)>i/(n+1)\\}$
$(i=1,2,\dots,n)$, where $F_{D}$ is the distribution function of $D$. For
heavy-tailed $D$, in particular, this choice of MR degrees converges faster
with $n$ to the intended $D$ than one based on rounding $np_{i}$
$(i=0,1,\dots)$ to nearest integers. Two choices of infectious period
distribution are used in the simulations: (i) $I\equiv 1$ (i.e. $\rm
P(I=1)=1$) and (ii) $I$ is $0$ or $\infty$, matched to have a common $p_{I}$.
Note that these yield the minimum and maximum asymptotic variances for fixed
$p_{I}$. The parameters of the degree distributions are chosen so that
$\mu_{D}=5$.

We first consider epidemics on an NSW network in which a fraction
$\epsilon=0.05$ of the population is initially infective. Table 1 shows
estimates of $\rho_{n}=n^{-1}\rm E[T^{(n)}_{\rm NSW}]$ and
$\sigma_{n}^{2}=n^{-1}\rm var[T^{(n)}_{\rm NSW}]$ for epidemics with
$p_{I}=0.3$ and various population size $n$, based on 100,000 simulations for
each set of parameters, together with the corresponding asymptotic
($n=\infty$) values given by Theorem 2.1. The table indicates that the
asymptotic approximations are useful for even moderate $n$. The approximations
are better for the model with $I\equiv 1$ than that with $I=0$ or $\infty$, as
one might expect since there is less variability in the process with $I\equiv
1$, and improve with increasing $\rm var(D)$. The approximations are
noticeably worse when $n=200$ than for the other values of $n$. Histograms of
the final size of 100,000 simulated epidemics on an NSW network with $n=200$
and $D\sim{\rm Geom}(1/6)$, together with the corresponding density of ${\rm
N}(n\rho,n\sigma^{2}_{\rm NSW})$ with $\rho$ and $\sigma^{2}_{\rm NSW}$ given
by Theorem 2.1, are shown in Figure 1. For the epidemic with $I\equiv 1$, the
asymptotic normal distrbution gives an excellent approximation, even though
$n$ is rather small. The approximation is markedly worse for model with $I=0$
or $\infty$, owing to the increased likelihood of small outbreaks. If
$\epsilon$ is held fixed, the probability of a small outbreak decreases
approximately exponentially with $n$, as the number of initial infectives is
proportional to $n$, and the approximation improves significantly,
particularly in the $I=0$ or $\infty$ model.

$\begin{array}[]{|c|c|c|c|c|c|}\hline\cr&&\vrule\lx@intercol\hfil I\equiv
1\hfil\lx@intercol\vrule\lx@intercol&\vrule\lx@intercol\hfil I=0\mbox{ or
}\infty\hfil\lx@intercol\vrule\lx@intercol\\\
\cline{3-6}\cr\mbox{D}&n&\rho_{n}&\sigma^{2}_{n}&\rho_{n}&\sigma^{2}_{n}\\\
\hline\cr&200&0.5097&2.6355&0.4259&9.9771\\\
&500&0.5284&2.3239&0.4957&9.6846\\\ &1000&0.5334&2.2180&0.5200&8.0686\\\
D\equiv 5&2000&0.5360&2.1794&0.5295&7.2082\\\
&5000&0.5374&2.1377&0.5351&6.7611\\\ &10000&0.5379&2.1177&0.5366&6.6402\\\
&\infty&0.5384&2.1187&0.5384&6.5200\\\
\hline\cr&200&0.5659&1.3149&0.4624&10.2610\\\
&500&0.5761&1.0808&0.5504&6.6926\\\ &1000&0.5789&1.0416&0.5708&4.0793\\\
\mbox{Po}(5)&2000&0.5803&1.0204&0.5766&3.5567\\\
&5000&0.5811&1.0095&0.5798&3.3708\\\ &10000&0.5814&1.0058&0.5807&3.3139\\\
&\infty&0.5817&1.0044&0.5817&3.2505\\\
\hline\cr&200&0.5160&0.4029&0.4132&8.5412\\\
&500&0.5184&0.3698&0.5055&3.0416\\\ &1000&0.5189&0.3690&0.5163&1.1780\\\
\mbox{Geom}(1/6)&2000&0.5194&0.3635&0.5180&1.0824\\\
&5000&0.5196&0.3635&0.5190&1.0607\\\ &10000&0.5196&0.3653&0.5194&1.0438\\\
&\infty&0.5197&0.3650&0.5197&1.0381\\\
\hline\cr&200&0.4957&0.4826&0.3900&8.8942\\\
&500&0.4985&0.4260&0.4828&4.1147\\\ &1000&0.4992&0.4237&0.4960&1.8536\\\
\mbox{Power}(1,13.796)&2000&0.4996&0.4177&0.4981&1.6862\\\
&5000&0.4999&0.4164&0.4992&1.6534\\\ &10000&0.4999&0.4182&0.4996&1.6483\\\
&\infty&0.5000&0.4180&0.5000&1.6372\\\ \hline\cr\end{array}$ Table 1:
Simulation results against theoretical (asymptotic) calculations for final
size of epidemics with $\epsilon=0.05$ and $p_{I}=0.3$ on NSW networks.

Turning to Theorem 2.2, Figure 2 shows simulations of the final size of
epidemics on an NSW network with one initial infective and $I$ constant. Note
that with a single initial infective there is always a non-negligible
probability of a minor outbreak, even if $n$ is large. As $n\to\infty$, the
probability of a major epidemic converges to the survival probability, $p_{\rm
maj}$ say, of a branching process; see, for example, [Ball and Sirl(2013)],
Section 2.1.1. Theorem 5.1. Moreover, as $I$ is constant, $p_{\rm maj}=\rho$;
see, for example, [Kenah and Robins(2007)]. Superimposed on each histogram in
Figure 2 is the density of ${\rm N}(n\rho,n\tilde{\sigma}^{2}_{\rm NSW})$
multiplied by $p_{\rm maj}$, which approximates the component of the
distribution of $T^{(n)}_{\rm NSW}$ corresponding to a major outbreak. Even
with the very small $n$, the approximation is very good when $D\sim{\rm
Geom}(1/6)$ or $D\sim{\rm Power}(1,13.796)$. For both of these degree
distributions there is a clear distiction between major and minor outbreaks.
That is not the case for the other two degree distributions, though the
approximation is still useful.

Figure 1: Histograms of 100,000 simulations of final size for epidemics with
$n=200,\epsilon=0.05$ and $p_{I}=0.3$ on NSW networks with $D\sim{\rm
Geom}(1/6)$, with asymptotic normal approximation superimposed.

Figure 2: Histograms of 100,000 simulations of final size for epidemics on NSW
networks with $I\equiv 1,p_{I}=1,n=100$ and one initial infective, with a
normal approximation superimposed; see text for details.

The upper panels of Figure 3 shows the dependence of $\rho$ (left panel) and
$\tilde{\sigma}^{2}_{\rm MR}$ and $\tilde{\sigma}^{2}_{\rm NSW}$ (right panel)
on $p_{I}$. The latter two are for the model with $I$ constant. (Recall that,
given $p_{I}$, the scaled asymptotic mean $\rho$ is independent of the
distribution of $I$.) The asymptotic scaled variances all decrease with
$p_{I}$ and converge to the asymptotic scaled variance of the giant component
on the relevant graph as $p_{I}\uparrow 1$ (cf. Remark 2.6). The asymptotic
scaled variances tend to $\infty$ as $p_{I}\downarrow p_{C}$, where
$p_{C}=\mu_{\tilde{D}-1}^{-1}$ is the critical value of $p_{I}$ so that
$R_{0}=1$. Note that $\tilde{\sigma}^{2}_{\rm NSW}\geq\tilde{\sigma}^{2}_{\rm
MR}$ for all choices of $D$, cf. Remark 2.5. The lower panels of Figure 3 show
plots of $(\tilde{\sigma}^{2}_{\rm NSW}-\tilde{\sigma}^{2}_{\rm
MR})/\tilde{\sigma}^{2}_{\rm MR}$ against $p_{I}$. Note the plots for the
Poisson and geometric degree distributions are both increasing with $p_{I}$,
while that for the ${\rm Power}(1,13.796)$ distribution is unimodal.

Figure 3: Dependence of asymptotic means and variances on $p_{I}$; see text
for details.

The final illustrations are concerned with percolation. Figure 4 shows plots
of estimates $\hat{\sigma}_{n}^{2}$ of the scaled variance $n^{-1}\rm
var(C^{(n)})$ of the size of the largest component, based on $n_{\rm
sim}=10,000$ simulations for each choice of parameters, together with $95$%
equal-tailed confidence intervals given by $\left[(n_{\rm
sim}-1)/q_{2},(n_{\rm sim}-1)/q_{1}\right]$, where $q_{1}$ and $q_{2}$ are
respectively the 2.5% and 97.5% quantiles of the $\chi^{2}_{n_{{\rm sim}-1}}$
distribution. In all cases $\pi=0.3$. The filled and unfilled markers
correspond to percolation on NSW and MR networks, respectively. The
$n\to\infty$ scaled variances, given by Theorem 2.7, are shown by horizontal
dashed and solid lines for NSW and MR networks, respectively. The figure
suggests that the asymptotic approximations are again generally good, even for
moderately-sized networks. For fixed $n$, the approximation is better for
${\rm Power}(1,13.796)$ distribution than for the ${\rm Po}(5)$ distribution.
The plot for site percolation when $D\sim{\rm Po}(5)$ is a bit odd, as
$\hat{\sigma}_{n}^{2}$ is not monotonic in $n$. This is explored further in
Figure 5, which is for percolation on NSW random graphs with $D\sim{\rm
Po}(5)$. Note that the distribution of $C^{(n)}$ is clearly bimodal for site
percolation with $n=200$ and $\pi=0.3$. The lower plots in Figure 4
demonstrate that increasing $n$ or $\pi$ alleviates the issue of small largest
components.

Figure 4: Plots of scaled variances $n^{-1}\rm var(C^{(n)})$ of the size of
the largest component in bond and site percolation on configuration model
random graphs; see text for details.

Figure 5: Histograms 100,000 simulations of size of largest component in
percolation on NSW random graphs having $D\sim{\rm Po}(5)$, with asymptotic
normal approximation superimposed.

Finally, Figure 6 shows histograms of the size of the largest component in
100,000 simulated bond and site percolations on an MR random graph with
$n=200,\pi=0.3$ and $D\sim{\rm Power}(1,13.796)$. Two asymptotic normal
approximations are superimposed. The solid lines are the densities of ${\rm
N}(n\rho,n\sigma^{2}_{\rm MR,B})$ (bond percolation) and ${\rm
N}(n\pi\rho,n\sigma^{2}_{\rm MR,S})$ (site percolation), with
$\rho,\sigma^{2}_{\rm MR,B}$ and $\sigma^{2}_{\rm MR,S}$ obtained by setting
$D\sim{\rm Power}(1,13.796)$ in Theorem 2.7. An improved approximation (dashed
lines) is obtained by instead setting $D$ to be the empirical ditribution of
$D^{(n)}_{1},D^{(n)}_{2},\dots,D^{(n)}_{n}$. The difference between the
approximations is more noticeable for bond percolation. (The support of the
histogram has been truncated slightly to make the difference clearer.) The
difference decreases with $n$ and is appreciably greater with heavy-tailed
degree distributions. Observe from Figures 5 and 6 that the asymptotic normal
approximation underestimates the left tail and overestimated the right tail of
the distribution of $C^{(n)}$. This phenomenon is present also in the
asymptotic normal approximation of the epidemic final size $T^{(n)}$.

Figure 6: Histograms 100,000 simulations of size of largest component in
percolation on MR random graphs having $n=200,\pi=0.3$ and $D\sim{\rm
Power}(1,13.796)$, with two asymptotic normal approximation superimposed; see
text for details.

##  4 Density dependent population processes

This section collects together some results for density dependent population
processes that are required in the paper. It is based on [Ethier and
Kurtz(1986)], Chapter 11, and [Pollett(1990)], though the statement of the
functional central limit theorem is slightly more general than that in those
references. The notation is local to this section.

For $n=1,2,\dots$, let $\\{\boldsymbol{X}^{(n)}(t):t\geq
0\\}=\\{(X^{(n)}_{1}(t),X^{(n)}_{2}(t),\dots,X^{(n)}_{p}(t)):t\geq 0\\}$ be a
continuous-time Markov chain with state space $S^{(n)}\subseteq\mathbb{Z}^{p}$
and transition intensities of the form

$q^{(n)}(\boldsymbol{i},\boldsymbol{i}+\boldsymbol{l})=n\beta^{(n)}_{\boldsymbol{l}}(n^{-1}\boldsymbol{i})\qquad(\boldsymbol{i}\in
S^{(n)},\boldsymbol{l}\in\Delta),$ (4.1)

where $\Delta$ is the set of possible jumps from a typical state
$\boldsymbol{i}=(i_{1},i_{2},\dots,i_{p})$ and the
$\beta^{(n)}_{\boldsymbol{l}}:E\to\mathbb{R}$ are nonnegative functions
defined on an open set $E\subseteq\mathbb{R}^{p}$. We assume that $\Delta$ is
finite. The theory in [Ethier and Kurtz(1986)], Chapter 11, and
[Pollett(1990)] allows $\Delta$ to be infinite but only the finite case is
required in our application and the results are easier to state in that
setting. Suppose that
$\beta_{\boldsymbol{l}}(\boldsymbol{x})=\lim_{n\to\infty}\beta^{(n)}_{\boldsymbol{l}}(\boldsymbol{x})$
exists for all $\boldsymbol{l}\in\Delta$ and all $\boldsymbol{x}\in E$; the
corresponding family of processes is then called asymptotically density
dependent by [Pollett(1990)]. In [Ethier and Kurtz(1986)], Chapter 11, a
family of processes which satisfies (4.1) with $\beta^{(n)}_{\boldsymbol{l}}$
replaced by $\beta_{\boldsymbol{l}}$ is called a density dependent family, and
it is noted that the results usually carry over with little additional effort
to the more general form where

$q^{(n)}(\boldsymbol{i},\boldsymbol{i}+\boldsymbol{l})=n\left[\beta_{\boldsymbol{l}}(n^{-1}\boldsymbol{i})+O(n^{-1})\right]\qquad(\boldsymbol{i}\in
S^{(n)},\boldsymbol{l}\in\Delta).$

Let

$F(\boldsymbol{x})=\sum_{\boldsymbol{l}\in\Delta}\boldsymbol{l}\beta_{\boldsymbol{l}}(\boldsymbol{x})\quad\mbox{and}\quad
F^{(n)}(\boldsymbol{x})=\sum_{\boldsymbol{l}\in\Delta}\boldsymbol{l}\beta_{\boldsymbol{l}}^{(n)}(\boldsymbol{x})\qquad(\boldsymbol{x}\in
E).$ (4.2)

The following weak law of large numbers follows from [Ethier and Kurtz(1986)],
Theorem 11.2.1, allowing for random initial conditions and asymptotic density
dependence; see also [Kurtz(1970)], Theorem 3.1 and [Pollett(1990)], Theorem
3.1.

###### Theorem 4.1

Suppose that for each compact $K\in H$,

$\sum_{\boldsymbol{l}\in\Delta}\sup_{\boldsymbol{x}\in
K}|\boldsymbol{l}|\beta_{\boldsymbol{l}}(\boldsymbol{x})<\infty,\qquad\lim_{n\to\infty}\sup_{\boldsymbol{x}\in
K}|F^{(n)}(\boldsymbol{x})-F(\boldsymbol{x})|=0$

and there exists $M_{K}>0$ such that

$|F(\boldsymbol{x})-F(\boldsymbol{y})|<M_{K}|\boldsymbol{x}-\boldsymbol{y}|\qquad\mbox{for
all }\boldsymbol{x},\boldsymbol{y}\in E.$

Suppose also that $n^{-1}\boldsymbol{X}^{(n)}(0)\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\boldsymbol{x}_{0}$ as $n\to\infty$. Let
$\boldsymbol{x}(t)$ $(t\geq 0)$ be given by

$\boldsymbol{x}(t)=\boldsymbol{x}_{0}+\int_{0}^{t}F(\boldsymbol{x}(u))\,{\rm
d}u.$

Then, for all $t\geq 0$,

$\sup_{0\leq u\leq
t}|n^{-1}\boldsymbol{X}^{(n)}(u)-\boldsymbol{x}(u)|\stackrel{{\scriptstyle
p}}{{\longrightarrow}}0\qquad\mbox{as }n\to\infty.$

Write
$F(\boldsymbol{x})=\left(F_{1}(\boldsymbol{x}),F_{2}(\boldsymbol{x}),\dots,F_{p}(\boldsymbol{x})\right)$
and $\boldsymbol{x}=(x_{1},x_{2},\dots,x_{p})$. Let $\partial
F(\boldsymbol{x})$ and $G(\boldsymbol{x})$ be the $p\times p$ matrix functions
defined by

$\partial F(\boldsymbol{x})=\left[\dfrac{\partial F_{i}}{\partial
x_{j}}(\boldsymbol{x})\right]\qquad\mbox{and}\qquad
G(\boldsymbol{x})=\sum_{\boldsymbol{l}\in\Delta}\boldsymbol{l}^{\top}\boldsymbol{l}\beta_{\boldsymbol{l}}(\boldsymbol{x}).$

Further, let $\Phi(t,u)=[\phi_{ij}(t,u)]$ $(0\leq u\leq t<\infty)$ be the
solution of the matrix differential equation

$\dfrac{\partial}{\partial t}\Phi(t,u)=\partial
F(\boldsymbol{x}(t))\Phi(t,u),\quad\Phi(u,u)=I.$

Note that

$\phi_{ij}(t,u)=\dfrac{\partial x_{i}(t-u)}{\partial
x_{j}(0)}\qquad(i,j=1,2,\dots,p).$

The following functional central limit theorem follows from [Ethier and
Kurtz(1986)], Theorem 11.2.3, (see also [Kurtz(1970)], Theorem 3.5, and
[Pollett(1990)], Theorem 3.2) allowing for asymptotically random initial
conditions. Let $\Rightarrow$ denote weak convergence in the space of right-
continuous functions $f:[0,\infty)\to\mathbb{R}^{p}$ having limits from the
left (i.e. càdlàg functions), endowed with the Skorohod metric.

###### Theorem 4.2

Suppose that for each compact $K\in H$,

$\sum_{\boldsymbol{l}\in\Delta}|\boldsymbol{l}|^{2}\sup_{\boldsymbol{x}\in
K}\beta_{\boldsymbol{l}}(\boldsymbol{x})<\infty\qquad\mbox{and}\qquad\lim_{n\to\infty}\sqrt{n}\sup_{\boldsymbol{x}\in
K}|F^{(n)}(\boldsymbol{x})-F(\boldsymbol{x})|=0,$

and that $\beta_{\boldsymbol{l}}$ $(\boldsymbol{l}\in\Delta)$ and $\partial F$
are continuous. Suppose also that

$\sqrt{n}\left(n^{-1}\boldsymbol{X}^{(n)}(0)-\boldsymbol{x}(0)\right)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}\boldsymbol{V}(0)\qquad\mbox{as }n\to\infty,$

where $\boldsymbol{V}(0)\sim N(\boldsymbol{0},\Sigma_{0})$. Then, as
$n\to\infty$,

$\left\\{\sqrt{n}\left(n^{-1}\boldsymbol{X}^{(n)}(t)-\boldsymbol{x}(t)\right):t\geq
0\right\\}\Rightarrow\\{\boldsymbol{V}(t):t\geq 0\\},$ (4.3)

where $\\{\boldsymbol{V}(t):t\geq 0\\}$ is a zero-mean Gaussian process with
covariance function given, for $t_{1},t_{2}\geq 0$, by

$\displaystyle{\rm
cov}\left(\boldsymbol{V}(t_{1}),\boldsymbol{V}(t_{2})\right)=$
$\displaystyle\Phi(t_{1},0)\Sigma_{0}\Phi(t_{2},0)^{\top}$ (4.4)
$\displaystyle\qquad+\int_{0}^{\min(t_{1},t_{2})}\Phi(t_{1},u)G(\boldsymbol{x}(u))\Phi(t_{2},u)^{\top}\,{\rm
d}u.$

We also require the following theorem concerned with the first exit of
$\\{\boldsymbol{X}^{(n)}(t):t\geq 0\\}$ from a region of its state space. The
theorem is a subset of [Ethier and Kurtz(1986)], Theorem 11.4.1, expressed now
for asymptotically density dependent processes having random initial
conditions.

###### Theorem 4.3

Suppose that the conditions of Theorem 4.2 are satisfied. Let
$\varphi:\mathbb{R}^{p}\to\mathbb{R}$ be continuously differentiable. Suppose
that $\varphi(\boldsymbol{x}(0))>0$. Let $\tau^{(n)}=\inf\\{t\geq
0:\varphi(n^{-1}\boldsymbol{X}^{(n)}(t))\leq 0\\}$ and $\tau=\inf\\{t\geq
0:\varphi(\boldsymbol{x}(t))\leq 0\\}$. Suppose that $\tau<\infty$ and
$\nabla\varphi(\boldsymbol{x}(\tau))\cdot F(\boldsymbol{x}(\tau))<0$, where
$\cdot$ denotes inner vector product. Then, as $n\to\infty$,

$\sqrt{n}\left(n^{-1}\boldsymbol{X}^{(n)}(\tau^{(n)})-\boldsymbol{x}(\tau)\right)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}\boldsymbol{V}(\tau)-\frac{\nabla\varphi(\boldsymbol{x}(\tau)\cdot\boldsymbol{V}(\tau)}{\nabla\varphi(\boldsymbol{x}(\tau)\cdot
F(\boldsymbol{x}(\tau))}F(\boldsymbol{x}(\tau)).$

###### Corollary 4.4

Suppose that the conditions of Theorem 4.3 are satisfied. Let

$B=I-\frac{F(\boldsymbol{x}(\tau))\bigotimes\nabla\varphi(\boldsymbol{x}(\tau))}{\nabla\varphi(\boldsymbol{x}(\tau))\cdot
F(\boldsymbol{x}(\tau))},$

where $\bigotimes$ denotes outer vector product. Then,

$\sqrt{n}\left(n^{-1}\boldsymbol{X}^{(n)}(\tau^{(n)})-\boldsymbol{x}(\tau)\right)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}{\rm
N}(\boldsymbol{0},B\Sigma(\tau)B^{\top})\qquad\mbox{as }n\to\infty,$

where

$\Sigma(\tau)=\Phi(\tau,0)\Sigma_{0}\Phi(\tau,0)^{\top}+\int_{0}^{\tau}\Phi(\tau,u)G(\boldsymbol{x}(u))\Phi(\tau,u)^{\top}\,{\rm
d}u.$

  * Proof.

Corollary 4.4 follows immediately from Theorems 4.2 and 4.3 on noting that

$\boldsymbol{V}(\tau)-\frac{\nabla\varphi(\boldsymbol{x}(\tau)\cdot\boldsymbol{V}(\tau)}{\nabla\varphi(\boldsymbol{x}(\tau)\cdot
F(\boldsymbol{x}(\tau))}F(\boldsymbol{x}(\tau))=\boldsymbol{V}(\tau)B^{\top}.$

$\Box$

##  5 Proofs

###  5.1 Alternative construction of final size $T^{(n)}$

We describe first the well-known construction of the final outcome of the
epidemic $\mathcal{E}^{(n)}$ using a random directed graph. We then use that
construction to show that $T^{(n)}$ can be realised using the location of the
first exit of an asymptotically density dependent population process from a
given region.

Given a realisation of $\mathcal{G}^{(n)}$, construct a directed random graph
$\tilde{\mathcal{G}}^{(n)}$, having vertex set
$\mathcal{N}^{(n)}=\\{1,2,\dots,n\\}$, as follows. For each $i=1,2,\dots,n$,
by sampling from its infectious period distribution $I$ and then the relevant
Poisson processes, draw up a list of individuals $i$ would make contact with
if $i$ were to become infected. Then, for each ordered pair of individuals,
$(i,j)$ say, a directed edge from $i$ to $j$ is present in
$\tilde{\mathcal{G}}^{(n)}$ if and only if $j$ is in $i$’s list. Let
$\mathcal{I}^{(n)}$ denote the set of initial infectives in
$\mathcal{E}^{(n)}$. For distinct $i,j\in\mathcal{N}^{(n)}$, write $i\leadsto
j$ if and only if there is a chain of directed edges from $i$ to $j$ in
$\tilde{\mathcal{G}}^{(n)}$. Let $\mathcal{T}^{(n)}$ be the set of initial
susceptibles that are infected in $\mathcal{E}^{(n)}$. Then
$\mathcal{T}^{(n)}=\\{j\in\mathcal{N}^{(n)}\setminus\mathcal{I}^{(n)}:i\leadsto
j\mbox{ for some i }\in\mathcal{I}^{(n)}\\}$ and the final size $T^{(n)}$ of
$\mathcal{E}^{(n)}$ is given by the cardinality of $\mathcal{T}^{(n)}$.

Note that $\mathcal{T}^{(n)}$ is determined purely by the presence/absence of
directed edges in $\tilde{\mathcal{G}}^{(n)}$ and does not depend on the times
of the corresponding infections in $\mathcal{E}^{(n)}$. (This implies that the
distribution of $\mathcal{T}^{(n)}$, and hence also $T^{(n)}$ is invariant to
the introduction of a latent/exposed period into the model, i.e. the time
elapsing after infection of an individual before it is able to infect other
individuals.) It follows that the final outcome $\mathcal{T}^{(n)}$ has the
same distribution as that of a related epidemic $\tilde{\mathcal{E}}^{(n)}$,
with set of initial infectives $\mathcal{I}^{(n)}$, in which for any
infective, $i$ say, it is determined upon infection which, if any, of its
neighbours $i$ will contact and those contacts take place at the first points
of independent Poisson processes, each having rate $1$. More precisely,
suppose $i$ is infected at time $t_{0}$ and $i$ has $d$ neighbours,
$i_{1},i_{2},\dots,i_{d}$ say. Let $I_{i}$ be a realisation of $I$ and, given
$I_{i}$, let $\chi_{i1},\chi_{i2},\dots,\chi_{id}$ be i.i.d. Bernoulli random
variables with success probability $1-\exp(-\lambda I_{i})$. Let
$W_{i1},W_{i2},\dots,W_{id}$ be an independent set of i.i.d. unit-mean
exponential random random variables. Then $i$ contacts $i_{j}$ if and only if
$\chi_{ij}=1$, and in that case the contact occurs at time $t_{0}+W_{ij}$. Of
course, the $I$ and $W$ random variables for any set of distinct infectives
are mutually independent.

Given the degrees $D^{(n)}_{1},D^{(n)}_{2},\dots,D^{(n)}_{n}$ and the set of
initial infectives $\mathcal{I}^{(n)}$, the random graph $\mathcal{G}^{(n)}$
and the epidemic on it $\tilde{\mathcal{E}}^{(n)}$ can be constructed
simultaneously as follows (cf. [Ball and Neal(2008)]). The process starts with
no half-edge paired. The individuals in $\mathcal{I}^{(n)}$ become infected at
time $t=0$ and all other individuals are susceptible. When an individual is
infected, it determines immediately which, if any, of its half-edges it will
transmit infection along and when it will make those contacts, according to
the probabilistic law described above; the half-edges that the indivdiual will
infect along and not infect along are then called infective and recovered
half-edges, respectively. When infection is transmitted along a half-edge that
half-edge is paired with a half-edge chosen uniformly at random from all
unpaired half-edges, forming an edge in the network. If the chosen half-edge
is attached to a susceptible individual then that individual is infected and
determines immediately along which, if any, of its remaining half-edges it
will transmit infection. If the chosen half-edge is infective or recovered
then nothing happens, apart from the two half-edges being replaced by an edge.
The process continues until there is no infective half-edge remaining. (In the
epidemic on the NSW random graph, if
$D^{(n)}_{1}+D^{(n)}_{2}+\dots+D^{(n)}_{n}$ is odd then it is possible for the
process to end with one unpaired half-edge, which is infective, but we can
ignore that possibility because under the conditions of the theorems its
probability tends to $0$ as $n\to\infty$.)

Note that as we are interested only in the final outcome of the epidemic, it
is not necessary to keep track of the degrees of individuals to which
infective and recovered half-edges are attached; we need to know just the
total numbers of such half-edges. For $t\geq 0$, let $X^{(n)}_{i}(t)$ be the
number of degree-$i$ susceptible individuals at time $t$ $(i=0,1,\dots,d_{\rm
max})$, and let $Y^{(n)}_{E}(t)$ and $Z^{(n)}_{E}(t)$ be the total numbers of
infective and recovered half-edges, respectively, at time $t$. Let
$\\{\boldsymbol{W}^{(n)}(t)\\}=\\{\boldsymbol{W}^{(n)}(t):t\geq 0\\}$, where
$\boldsymbol{W}^{(n)}(t)=(\boldsymbol{X}^{(n)}(t),Y^{(n)}_{E}(t),Z^{(n)}_{E}(t))$
and
$\boldsymbol{X}^{(n)}(t)=(X^{(n)}_{0}(t),X^{(n)}_{1}(t),\dots,X^{(n)}_{d_{\rm
max}}(t))$.

The process $\\{\boldsymbol{W}^{(n)}(t)\\}$ is a continuous-time Markov chain,
whose initial state $\boldsymbol{W}^{(n)}(0)$ is random, even for the epidemic
on the MR random graph as the numbers of infective and recovered half-edges
created by the initial infectives are random. In the epidemic on the NSW
random graph, $X^{(n)}_{i}(t)$ $(i=0,1,\dots,d_{\rm max})$ are also random.
Before giving the possible transition intensities of
$\\{\boldsymbol{W}^{(n)}(t)\\}$ some more notation is required.

Let $H^{(n)}$ denote the state space of $\\{\boldsymbol{W}^{(n)}(t)\\}$ and
$\boldsymbol{n}=(n^{X}_{0},n^{X}_{1},\dots,n^{X}_{d_{\rm
max}},n^{Y}_{E},n^{Z}_{E})$ denote a typical element of $H^{(n)}$. Let
$n^{X}_{E}=\sum_{i=1}^{d_{\rm max}}in^{X}_{i}$. Define the unit vectors
$\boldsymbol{e}^{\rm S}_{i}$ $(i=0,1,\dots,d_{\rm max})$, $\boldsymbol{e}^{\rm
I}$ and $\boldsymbol{e}^{\rm R}$ on $H^{(n)}$, where, for example,
$\boldsymbol{e}^{\rm S}_{i}$ has a one in the $i$th suscpetible component and
zeros elsewhere. For $i=1,2,\dots,d_{\rm max}$ and $k=0,1,\dots,i-1$, let
$p_{i,k}$ be the probability that if a degree-$i$ susceptible is infected it
subsequently transmits infection along $k$ of its remaining $i-1$ half-edges.
(Note that a degree-$0$ susceptible cannot be infected.) Thus $p_{i,k}={\rm
P}(X=k)$, where $X\sim{\rm Bin}(i-1,1-\exp(-\lambda I))$. The transition
intensities of $\\{\boldsymbol{W}^{(n)}(t)\\}$ are as follows.

  1. (i)

For $i=1,2,\dots,d_{\rm max}$ and $k=0,1,\dots,i-1$, an infective half-edge is
paired with a degree-$i$ susceptible yielding $k$ infective half-edges and
$i-1-k$ recovered half-edges

$q^{(n)}(\boldsymbol{n},\boldsymbol{n}-\boldsymbol{e}^{\rm
S}_{i}+(k-1)\boldsymbol{e}^{\rm I}+(i-k-1)\boldsymbol{e}^{\rm
R})=\frac{n^{Y}_{E}in^{X}_{i}p_{i,k}}{n^{X}_{E}+n^{Y}_{E}+n^{Z}_{E}-1};$

  2. (ii)

an infective half-edge is paired with an infective half-edge

$q^{(n)}(\boldsymbol{n},\boldsymbol{n}-2\boldsymbol{e}^{\rm
I})=\frac{n^{Y}_{E}(n^{Y}_{E}-1)}{n^{X}_{E}+n^{Y}_{E}+n^{Z}_{E}-1};$

  3. (iii)

an infective half-edge is paired with a recovered half-edge

$q^{(n)}(\boldsymbol{n},\boldsymbol{n}-\boldsymbol{e}^{\rm
I}-\boldsymbol{e}^{\rm
R})=\frac{n^{Y}_{E}n^{Z}_{E}}{n^{X}_{E}+n^{Y}_{E}+n^{Z}_{E}-1}.$

Note that these transition intensities are independent of the population size
$n$. We index then by $n$ to connect with theory of density dependent
population processes.

Let $\tau^{(n)}=\inf\\{t\geq 0:Y^{(n)}_{E}(t)=0\\}$. Then
$X^{(n)}_{i}(\tau^{(n)})$ $(i=0,1,\dots,d_{\rm max})$ give the numbers of
susceptibles of the different degrees at the end of the epidemic and
$T^{(n)}=\sum_{i=1}^{d_{\rm max}}(X^{(n)}_{i}(0)-X^{(n)}_{i}(\tau^{(n)})$.

###  5.2 Random time-scale transformation

We wish to apply Corollary 4.4 to obtain a central limit theorem for $T^{(n)}$
but that corollary cannot be applied directly to
$\\{\boldsymbol{W}^{(n)}(t)\\}$ as $\tau^{(n)}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\infty$ as $n\to\infty$. Thus we consider the following
random time-scale transformation of $\\{\boldsymbol{W}^{(n)}(t)\\}$; cf.
[Ethier and Kurtz(1986)], page 467, [Watson(1980)] and [Janson et al.(2012)].

For $t\in[0,\tau^{(n)}]$, let

$A^{(n)}(t)=\int_{0}^{t}\frac{Y^{(n)}_{E}(u)}{X^{(n)}_{E}(u)+Y^{(n)}_{E}(u)+Z^{(n)}_{E}(u)-1}\,{\rm
d}u,$

where $X^{(n)}_{E}(u)=\sum_{i=1}^{d_{\rm max}}iX^{(n)}_{i}(u)$. Let
$\tilde{\tau}^{(n)}=A^{(n)}(\tau^{(n)})$. For $t\in[0,\tilde{\tau}^{(n)}]$,
let $U^{(n)}(t)=\inf\\{u\geq 0:A^{(n)}(u)=t\\}$ and
$\tilde{\boldsymbol{W}}^{(n)}(t)=\boldsymbol{W}^{(n)}\left(U^{(n)}(t)\right)$.
Then
$\\{\tilde{\boldsymbol{W}}^{(n)}(t)\\}=\\{\tilde{\boldsymbol{W}}^{(n)}(t):0\leq
t\leq\tilde{\tau}^{(n)}\\}$ is a continuous-time Markov chain with transition
intensities:

  1. (i)

for $i=1,2,\dots,d_{\rm max}$ and $k=0,1,\dots,i-1$,

$\tilde{q}^{(n)}(\boldsymbol{n},\boldsymbol{n}-\boldsymbol{e}^{\rm
S}_{i}+(k-1)\boldsymbol{e}^{\rm I}+(i-k-1)\boldsymbol{e}^{\rm
R})=in^{X}_{i}p_{i,k};$

  2. (ii)

$\tilde{q}^{(n)}(\boldsymbol{n},\boldsymbol{n}-2\boldsymbol{e}^{\rm
I})=n^{Y}_{E}-1;$

  3. (iii)

$\tilde{q}^{(n)}(\boldsymbol{n},\boldsymbol{n}-\boldsymbol{e}^{\rm
I}-\boldsymbol{e}^{\rm R})=n^{Z}_{E}.$

For $t\geq 0$, let
$\tilde{\boldsymbol{X}}^{(n)}(t)=(\tilde{X}^{(n)}_{0}(t),\tilde{X}^{(n)}_{1}(t),\dots,\tilde{X}^{(n)}_{d_{\rm
max}}(t))$, so
$\tilde{\boldsymbol{W}}^{(n)}(t)=(\tilde{\boldsymbol{X}}^{(n)}(t),\tilde{Y}^{(n)}_{E}(t),\tilde{Z}^{(n)}_{E}(t))$.

For reasons that will become clear later, we need to extend
$\\{\tilde{\boldsymbol{W}}^{(n)}(t)\\}$ so that it is defined beyond time
$\tilde{\tau}^{(n)}$ and allow $\tilde{Y}^{(n)}_{E}(t)$ and
$\tilde{Z}^{(n)}_{E}(t)$ to take negative values. Thus we add two new
transition intensities:

  1. (ii′)

for $n^{Y}_{E}<0$,

$\tilde{q}^{(n)}(\boldsymbol{n},\boldsymbol{n}+2\boldsymbol{e}^{\rm
I})=-(n^{Y}_{E}+1);$

  2. (iii′)

for $n^{Z}_{E}<0$,

$\tilde{q}^{(n)}(\boldsymbol{n},\boldsymbol{n}+\boldsymbol{e}^{\rm
I}+\boldsymbol{e}^{\rm R})=-n^{Z}_{E}.$

Note that the transition intensities (ii) and (iii) are defined only for
$n^{Y}_{E}>0$ and $n^{Z}_{E}>0$, respectively.

The state space of $\\{\tilde{\boldsymbol{W}}^{(n)}(t)\\}$ is a subset of

$\tilde{H}^{(n)}=\left([0,n]^{d_{\rm max}+1}\times(-\infty,nd_{\rm
max}]^{2}\right)\cap\mathbb{Z}^{d_{\rm max}+3}.$

Let $\boldsymbol{l}_{ik}^{(1)}=-\boldsymbol{e}^{\rm
S}_{i}+(k-1)\boldsymbol{e}^{\rm I}+(i-k-1)\boldsymbol{e}^{\rm R}$
$(i=1,2,\dots,d_{\rm max};k=0,1,\dots,i-1)$,
$\boldsymbol{l}^{(2)}_{+}=-2\boldsymbol{e}^{\rm I}$,
$\boldsymbol{l}^{(2)}_{-}=2\boldsymbol{e}^{\rm I}$,
$\boldsymbol{l}^{(3)}_{+}=-\boldsymbol{e}^{\rm I}-\boldsymbol{e}^{\rm R}$ and
$\boldsymbol{l}^{(3)}_{+}=\boldsymbol{e}^{\rm I}+\boldsymbol{e}^{\rm R}$. The
set of possible jumps of $\\{\tilde{\boldsymbol{W}}^{(n)}(t)\\}$ from a
typical state $\boldsymbol{n}$ is
$\Delta=\Delta_{1}\cup\Delta_{2}\cup\Delta_{3}$, where
$\Delta_{1}=\\{\boldsymbol{l}^{(1)}_{ik}:i=1,2,\dots,d_{\rm
max};k=0,1,\dots,i-1\\}$,
$\Delta_{2}=\\{\boldsymbol{l}^{(2)}_{+},\boldsymbol{l}^{(2)}_{-}\\}$ and
$\Delta_{3}=\\{\boldsymbol{l}^{(3)}_{+},\boldsymbol{l}^{(3)}_{-}\\}$. Let
$\boldsymbol{w}=(\boldsymbol{x},y_{E},z_{E})$, where
$\boldsymbol{x}=(x_{0},x_{1},\dots,x_{d_{\rm max}})$, and $H=[0,1]^{d_{\rm
max}+1}\times(-\infty,d_{\rm max}]^{2}$. The intensities of the jumps of
$\\{\tilde{\boldsymbol{W}}^{(n)}(t)\\}$ admit the form

$\tilde{q}^{(n)}(\boldsymbol{n},\boldsymbol{n}+\boldsymbol{l})=n\tilde{\beta}_{\boldsymbol{l}}^{(n)}(n^{-1}\boldsymbol{n})\qquad(\boldsymbol{n}\in\tilde{H}^{(n)},\boldsymbol{l}\in\Delta),$
(5.1)

where the functions $\tilde{\beta}_{\boldsymbol{l}}^{(n)}:H\to\mathbb{R}_{+}$
$(\boldsymbol{l}\in\Delta)$ are given by

$\tilde{\beta}_{\boldsymbol{l}}^{(n)}(\boldsymbol{w})=\begin{cases}\tilde{\beta}_{\boldsymbol{l}}(\boldsymbol{w})&\text{
if
}\boldsymbol{l}\in\Delta\setminus\\{\boldsymbol{l}^{(2)}_{+},\boldsymbol{l}^{(2)}_{-}\\},\\\
\tilde{\beta}_{\boldsymbol{l}}(\boldsymbol{w})-n^{-1}&\text{ if
}\boldsymbol{l}=\boldsymbol{l}^{(2)}_{+},\\\
\tilde{\beta}_{\boldsymbol{l}}(\boldsymbol{w})+n^{-1}&\text{ if
}\boldsymbol{l}=\boldsymbol{l}^{(2)}_{-},\end{cases}$ (5.2)

with

$\tilde{\beta}_{\boldsymbol{l}}(\boldsymbol{w})=\begin{cases}\tilde{\beta}_{ik}^{(1)}(\boldsymbol{x},y_{E},z_{E})=ix_{i}p_{i,k}&\text{
for }\boldsymbol{l}=\boldsymbol{l}_{ij}^{(1)}\in\Delta_{1},\\\
\tilde{\beta}^{(2)}_{+}(\boldsymbol{x},y_{E},z_{E})=y_{E}1_{\\{y_{E}>0\\}}&\text{
for }\boldsymbol{l}=\boldsymbol{l}^{(2)}_{+},\\\
\tilde{\beta}^{(2)}_{-}(\boldsymbol{x},y_{E},z_{E})=-y_{E}1_{\\{y_{E}<0\\}}&\text{
for }\boldsymbol{l}=\boldsymbol{l}^{(2)}_{-},\\\
\tilde{\beta}^{(3)}_{+}(\boldsymbol{x},y_{E},z_{E})=z_{E}1_{\\{z_{E}>0\\}}&\text{
for }\boldsymbol{l}=\boldsymbol{l}^{(3)}_{+},\\\
\tilde{\beta}^{(3)}_{-}(\boldsymbol{x},y_{E},z_{E})=-z_{E}1_{\\{z_{E}<0\\}}&\text{
for }\boldsymbol{l}=\boldsymbol{l}^{(3)}_{-}.\end{cases}$ (5.3)

The family of processes $\\{\tilde{\boldsymbol{W}}^{(n)}(t)\\}$ is
asymptotically density dependent (see Section 4).

###  5.3 Proof of Theorem 2.1

Note that
$\boldsymbol{W}^{(n)}(\tau^{(n)})=\tilde{\boldsymbol{W}}^{(n)}(\tilde{\tau}^{(n)})$,
so the final size of the epidemic is given by

$T^{(n)}=\sum_{i=0}^{d_{\rm
max}}\left(\tilde{X}^{(n)}_{i}(0)-\tilde{X}^{(n)}_{i}(\tilde{\tau}^{(n)})\right).$
(5.4)

Note also that $\tilde{\tau}^{(n)}=\inf\\{t\geq
0:\tilde{Y}^{(n)}_{E}(t)=0\\}$. We use Corollary 4.4 to obtain a central limit
theorem for $\tilde{\boldsymbol{W}}^{(n)}(\tilde{\tau}^{(n)})$, and hence for
$T^{(n)}$. The asymptotic variance matrix in the central limit theorem for
$\tilde{\boldsymbol{W}}^{(n)}(\tilde{\tau}^{(n)})$ is not in closed form.
However, we derive a closed-form expression for the asymptotic variance of
$T^{(n)}$. The main concepts of the proof are presented here, with some
detailed but elementary calculations deferred to Appendix A.

####  5.3.1 Deterministic model

As at (4.2), define the drift function
$\tilde{F}(\boldsymbol{w})=\sum_{\boldsymbol{l}\in\Delta}\boldsymbol{l}\tilde{\beta}_{\boldsymbol{l}}(\boldsymbol{w})$.
Using (5.3),

$\displaystyle\tilde{F}(\boldsymbol{w})$ $\displaystyle=-\sum_{i=1}^{d_{\rm
max}}ix_{i}\boldsymbol{e}^{\rm S}_{i}+\left\\{\sum_{i=1}^{d_{\rm
max}}ix_{i}[(i-1)p_{I}-1]-2y_{E}-z_{E}\right\\}\boldsymbol{e}^{\rm I}$ (5.5)
$\displaystyle\qquad\qquad+\left[\sum_{i=1}^{d_{\rm
max}}i(i-1)x_{i}q_{I}-z_{E}\right]\boldsymbol{e}^{\rm R}.$

For $t\geq 0$, let
$\tilde{\boldsymbol{w}}(t)=(\tilde{x}_{0}(t),\tilde{x}_{1}(t),\dots,\tilde{x}_{d_{\rm
max}}(t),\tilde{y}_{E}(t),\tilde{z}_{E}(t))$ be defined by

$\tilde{\boldsymbol{w}}(t)=\tilde{\boldsymbol{w}}(0)+\int_{0}^{t}\tilde{F}(\boldsymbol{w}(u))\,{\rm
d}u.$ (5.6)

Thus $\tilde{\boldsymbol{w}}(t)$ satisfies the differential equations

$\displaystyle\dfrac{d\tilde{x}_{i}}{dt}$
$\displaystyle=-i\tilde{x}_{i}\quad(i=0,1,\dots,d_{\rm max}),$ (5.7)
$\displaystyle\dfrac{d\tilde{y}_{E}}{dt}$ $\displaystyle=\sum_{i=2}^{d_{\rm
max}}i(i-1)p_{I}\tilde{x}_{i}-\tilde{x}_{E}-2\tilde{y}_{E}-\tilde{z}_{E},$
(5.8) $\displaystyle\dfrac{d\tilde{z}_{E}}{dt}$
$\displaystyle=\sum_{i=2}^{d_{\rm
max}}i(i-1)q_{I}\tilde{x}_{i}-\tilde{z}_{E},$ (5.9)

where $\tilde{x}_{E}=\sum_{i=1}^{d_{\rm max}}i\tilde{x}_{i}$, having solution
(see Appendix A.1)

$\displaystyle\tilde{x}_{i}(t)$ $\displaystyle=\tilde{x}_{i}(0){\rm
e}^{-it}\quad(i=0,1,\dots,d_{\rm max}),$ (5.10)
$\displaystyle\tilde{y}_{E}(t)$
$\displaystyle=(\tilde{x}_{E}(0)+\tilde{y}_{E}(0)+\tilde{z}_{E}(0)){\rm
e}^{-2t}-[\tilde{z}_{E}(0)+q_{I}\tilde{x}_{E}(0)]{\rm e}^{-t}$ (5.11)
$\displaystyle\qquad-p_{I}\sum_{i=1}^{d_{\rm max}}i\tilde{x}_{i}(0){\rm
e}^{-it},$ $\displaystyle\tilde{z}_{E}(t)$
$\displaystyle=[\tilde{z}_{E}(0)+q_{I}\tilde{x}_{E}(0)]{\rm
e}^{-t}-q_{I}\sum_{i=1}^{d_{\rm max}}i\tilde{x}_{i}(0){\rm e}^{-it}.$ (5.12)

Let $\tilde{\eta}_{E}(t)=\tilde{x}_{E}(t)+\tilde{y}_{E}(t)+\tilde{z}_{E}(t)$.
Note that, for $t\geq 0$,

$\tilde{\eta}_{E}(t)=\tilde{\eta}_{E}(0){\rm e}^{-2t}.$ (5.13)

####  5.3.2 Initial conditions

Consider first the epidemic on the MR random graph. Recall that in the
epidemic $\mathcal{E}^{(n)}$, for $i=0,1,\dots,d_{\rm max}$, there are
$v^{(n)}_{i}$ individuals of degree $i$, of whom $a^{(n)}_{i}$ are initially
infected. Also, $\left(\tilde{Y}^{(n)}_{E}(0),\tilde{Z}^{(n)}_{E}(0)\right)$
is given by the total numbers of infective and recovered half-edges created by
the initial infectives. Thus $\tilde{X}^{(n)}_{i}(0)=v^{(n)}_{i}-a^{(n)}_{i}$
($i=0,1,\dots,d_{\rm max}$) and

$\left(\tilde{Y}^{(n)}_{E}(0),\tilde{Z}^{(n)}_{E}(0)\right)=\sum_{i=1}^{d_{\rm
max}}\sum_{j=1}^{a^{(n)}_{i}}\left(Y_{ij},i-Y_{ij}\right),$

where $Y_{ij}$ $(i=1,2,\dots,d_{\rm max};j=1,2,\dots)$ are independent and
$Y_{ij}\sim{\rm Bin}(i,1-\exp(-\lambda I))$. Let
$\epsilon_{E}=\sum_{i=1}^{d_{\rm max}}i\epsilon_{i}$. Then, recalling that
$\epsilon^{(n)}_{i}=n^{-1}a^{(n)}_{i}$,

$\displaystyle\sqrt{n}$
$\displaystyle\left[n^{-1}\left(\tilde{Y}^{(n)}_{E}(0),\tilde{Z}^{(n)}_{E}(0)\right)-\epsilon_{E}(p_{I},q_{I})\right]$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{i=1}^{d_{\rm
max}}\left[\sum_{j=1}^{a^{(n)}_{i}}\left(Y_{ij},i-Y_{ij}\right)-a^{(n)}_{i}i(p_{I},q_{I})\right]+\sqrt{n}\sum_{i=1}^{d_{\rm
max}}i(\epsilon^{(n)}_{i}-\epsilon_{i}).$

For $i=1,2,\dots,d_{\rm max}$, let $\sigma_{Y,i}^{2}={\rm var}(Y_{i1})$. Now
$\lim_{n\to\infty}\sqrt{n}(\epsilon^{(n)}_{i}-\epsilon_{i})=0$
$(i=0,1,\dots,d_{\rm max})$, by assumption, so the central limit theorem and
Slutsky’s theorem imply that

$\sqrt{n}\left[n^{-1}\left(\tilde{Y}^{(n)}_{E}(0),\tilde{Z}^{(n)}_{E}(0)\right)-\epsilon_{E}(p_{I},q_{I})\right]\stackrel{{\scriptstyle
D}}{{\longrightarrow}}{\rm
N}\left(\boldsymbol{0},\begin{bmatrix}\sigma_{Y}^{2}&-\sigma_{Y}^{2}\\\
-\sigma_{Y}^{2}&\sigma_{Y}^{2}\end{bmatrix}\right)$

as $n\to\infty$, where

$\sigma_{Y}^{2}=\sum_{i=1}^{d_{\rm max}}\epsilon_{i}\sigma_{Y,i}^{2}.$ (5.14)

(A closed-form expression for $\sigma_{Y}^{2}$ is given by (A.27) in Appendix
A.3.2.) Since $\tilde{\boldsymbol{X}}^{(n)}(0)$ is non-random, it follows
using (2.1) that

$\sqrt{n}\left(n^{-1}\tilde{\boldsymbol{W}}^{(n)}(0)-\tilde{\boldsymbol{w}}(0)\right)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}{\rm N}\left(\boldsymbol{0},\Sigma^{\rm
MR}_{0}\right)\qquad\mbox{as }n\to\infty,$ (5.15)

where

$\tilde{\boldsymbol{w}}(0)=(p_{1}-\epsilon_{1},p_{2}-\epsilon_{2},\cdots,p_{d_{\rm
max}}-\epsilon_{d_{\rm max}},\epsilon_{E}p_{I},\epsilon_{E}q_{I})$ (5.16)

and

$\Sigma^{\rm MR}_{0}=\begin{bmatrix}0&\boldsymbol{0}&\boldsymbol{0}\\\
\boldsymbol{0}&\sigma_{Y}^{2}&-\sigma_{Y}^{2}\\\
\boldsymbol{0}&-\sigma_{Y}^{2}&\sigma_{Y}^{2}\end{bmatrix}.$ (5.17)

Turning to the epidemic on the NSW random graph, recall that in
$\mathcal{E}^{(n)}$ the number of initial infectives $a^{(n)}$ is prescribed
and the $a^{(n)}$ initial infectives are chosen by sampling uniformly at
random without replacement from the population. Thus the network can be
constructed using two independent sets of i.i.d. copies of $D$, viz.
$D_{1}^{\prime},D_{2}^{\prime},\dots,D_{n-a^{(n)}}^{\prime}$ for the initial
susceptibles and $D_{1},D_{2},\dots,D_{a^{(n)}}$ for the initial infectives.
Let $(Y_{E},Z_{E})$ be the bivariate random variable obtained by first
sampling $D$ and then letting $(Y_{E},Z_{E})=(Y_{E},D-Y_{E})$, where
$Y_{E}|D~{}\sim{\rm Bin}(D,1-\exp(-\lambda I))$. Let $\sigma_{Y_{E}}^{2}={\rm
var}(Y_{E})$, $\sigma_{Y_{E},Z_{E}}={\rm cov}(Y_{E},Z_{E})$ and
$\sigma_{Y_{E}}^{2}={\rm var}(Y_{E})$. (Closed-form expressions for
$\sigma_{Y_{E}}^{2}$, $\sigma_{Z_{E}}^{2}$ and $\sigma_{Y_{E},Z_{E}}$ are
given in (A.3.3)-(A.40) in Appendix A.3.3.) Let
$\boldsymbol{p}=(p_{0},p_{1},\dots,p_{d_{\rm max}})$ and $\Sigma_{XX}$ be the
$(d_{\rm max}+1)\times(d_{\rm max}+1)$ matrix with elements

$(\Sigma_{XX})_{ij}=\begin{cases}-p_{i}p_{j}&\text{ if }i\neq j,\\\
p_{i}(1-p_{i})&\text{ if }i=j.\end{cases}$ (5.18)

Recalling that $\lim_{n\to\infty}\sqrt{n}(\epsilon^{(n)}-\epsilon)=0$, where
$\epsilon^{(n)}=n^{-1}a^{(n)}$, a similar argument to the above shows that

$\sqrt{n}\left(n^{-1}\tilde{\boldsymbol{W}}^{(n)}(0)-\tilde{\boldsymbol{w}}(0)\right)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}{\rm N}\left(\boldsymbol{0},\Sigma^{\rm
NSW}_{0}\right)\qquad\mbox{as }n\to\infty,$ (5.19)

where

$\tilde{\boldsymbol{w}}(0)=\left((1-\epsilon)\boldsymbol{p},\epsilon\mu_{D}p_{I},\epsilon\mu_{D}q_{I}\right)$
(5.20)

and

$\Sigma^{\rm
NSW}_{0}=\begin{bmatrix}(1-\epsilon)\Sigma_{XX}&\boldsymbol{0}&\boldsymbol{0}\\\
\boldsymbol{0}&\epsilon\sigma_{Y_{E}}^{2}&\epsilon\sigma_{Y_{E},Z_{E}}\\\
\boldsymbol{0}&\epsilon\sigma_{Y_{E},Z_{E}}&\epsilon\sigma_{Z_{E}}^{2}\end{bmatrix}.$
(5.21)

####  5.3.3 Central limit theorem

Noting from (5.2) that

$\max_{\boldsymbol{l}\in\Delta}\sup_{\boldsymbol{w}\in
H}|\tilde{\beta}_{\boldsymbol{l}}^{(n)}(\boldsymbol{w})-\tilde{\beta}_{\boldsymbol{l}}(\boldsymbol{w})|=n^{-1},$

it is easily checked that $\\{\tilde{\boldsymbol{W}}^{(n)}(t)\\}$ satisfies
the conditions of Theorem 4.3. Now $\tilde{\tau}^{(n)}=\inf\\{t\geq
0:\varphi(n^{-1}\tilde{\boldsymbol{W}}^{(n)}(t))\leq 0\\}$, where
$\varphi(\boldsymbol{w})=\varphi(\boldsymbol{x},y_{E},z_{E})=y_{E}$.

Suppose that $\tilde{x}_{i}(0)=p_{i}-\epsilon_{i}$ $(i=0,1,\dots,d_{\rm
max})$, $\tilde{y}_{E}(0)=p_{I}\sum_{i=1}^{d_{\rm max}}i\epsilon_{i}$ and
$\tilde{z}_{E}(0)=q_{I}\sum_{i=1}^{d_{\rm max}}i\epsilon_{i}$, so
$\tilde{\eta}_{E}(0)=\mu_{D}$. Let $\tilde{\tau}=\inf\\{t\geq
0:\varphi(\tilde{\boldsymbol{w}}(t))\leq 0\\}=\inf\\{t\geq
0:\tilde{y}_{E}(t)=0\\}$. Then it follows from (5.11) that $\tilde{\tau}$
satisfies the equation

${\rm e}^{-\tilde{\tau}}-q_{I}-\mu_{D}^{-1}p_{I}f_{D_{\epsilon}}^{(1)}({\rm
e}^{-\tilde{\tau}})=0.$ (5.22)

We show in Appendix A.2 that, under the conditions of Theorem 2.1, the
equation (5.22) has a unique solution in $(0,\infty)$. Note that $z={\rm
e}^{-\tilde{\tau}}$, where $z$ is defined at (2.2). Also, using (5.10) the
deterministic final fraction of the population that is susceptible is given by

$\sum_{i=0}^{d_{\rm max}}\tilde{x}_{i}(\tilde{\tau})=\sum_{i=0}^{d_{\rm
max}}(p_{i}-\epsilon_{i}){\rm e}^{-i\tilde{\tau}}=f_{D_{\epsilon}}({\rm
e}^{-\tilde{\tau}}).$

The corresponding deterministic final size is $\rho=1-\epsilon-
f_{D_{\epsilon}}({\rm e}^{-\tilde{\tau}})$, agreeing with (2.3).

Let
$a(\tilde{\tau})=\nabla\varphi(\boldsymbol{w}(\tilde{\tau}))\cdot\tilde{F}(\boldsymbol{w}(\tilde{\tau}))$.
Then, using (5.6), and noting that $\tilde{y}_{E}(\tilde{\tau})=0$,

$\displaystyle a(\tilde{\tau})$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{d_{\rm
max}}i\tilde{x}_{i}(\tilde{\tau})[(i-1)p_{I}-1]-2\tilde{y}_{E}(\tilde{\tau})-\tilde{z}_{E}(\tilde{\tau})$
$\displaystyle=$ $\displaystyle p_{I}\sum_{i=1}^{d_{\rm
max}}i(i-1)\tilde{x}_{i}(\tilde{\tau})-\tilde{\eta}_{E}(\tilde{\tau}).$

Thus, using (5.10) and (5.13),

$a(\tilde{\tau})={\rm
e}^{-2\tilde{\tau}}\left(p_{I}f_{D_{\epsilon}}^{(2)}({\rm
e}^{-\tilde{\tau}})-\mu_{D}\right),$ (5.23)

since $\tilde{\eta}_{E}(0)=\mu_{D}$. We show in Appendix A.2 that
$a(\tilde{\tau})<0$, so we may apply Corollary 4.4.

Writing
$\tilde{\boldsymbol{w}}=(\tilde{w}_{0},\tilde{w}_{1},\dots,\tilde{w}_{p})$,
where $p=d_{\rm max}+2$, let $\tilde{\Phi}(t,u)=[\tilde{\phi}_{ij}(t,u)]$
$(0\leq u\leq t<\infty)$, where

$\tilde{\phi}_{ij}(t,u)=\dfrac{\partial\tilde{w}_{i}(t-u)}{\partial\tilde{w}_{j}(0)}\qquad(i,j=0,1,\dots,p).$
(5.24)

Also, let

$\tilde{\Sigma}(\tilde{\tau})=\tilde{\Phi}(\tilde{\tau},0)\Sigma_{0}\tilde{\Phi}(\tilde{\tau},0)^{\top}+\int_{0}^{\tilde{\tau}}\tilde{\Phi}(\tilde{\tau},u)\tilde{G}(\tilde{\boldsymbol{w}}(u))\tilde{\Phi}(\tilde{\tau},u)^{\top}\,{\rm
d}u,$ (5.25)

where

$\tilde{G}(\tilde{\boldsymbol{w}}(u))=\sum_{\boldsymbol{l}\in\Delta}\boldsymbol{l}^{\top}\boldsymbol{l}\tilde{\beta}_{\boldsymbol{l}}(\tilde{\boldsymbol{w}}(u)),$
(5.26)

and $\Sigma_{0}=\Sigma^{\rm MR}_{0}$ or $\Sigma^{\rm NSW}_{0}$ dependening on
whether the epidemic is on an MR or an NSW random graph. Then application of
Corollary 4.4 yields

$\sqrt{n}\left(n^{-1}\tilde{\boldsymbol{W}}^{(n)}(\tilde{\tau}^{(n)})-\tilde{\boldsymbol{w}}(\tilde{\tau})\right)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}{\rm
N}(\boldsymbol{0},B\tilde{\Sigma}(\tilde{\tau})B^{\top})\qquad\mbox{as
}n\to\infty,$ (5.27)

where

$B=I-\frac{\tilde{F}(\tilde{\boldsymbol{w}}(\tilde{\tau}))\bigotimes\nabla\varphi(\tilde{\boldsymbol{w}}(\tilde{\tau}))}{\nabla\varphi(\tilde{\boldsymbol{w}}(\tilde{\tau}))\cdot\tilde{F}(\tilde{\boldsymbol{w}}(\tilde{\tau}))}.$

Thus, recalling (5.4),

$\sqrt{n}\left(n^{-1}T^{(n)}-\rho\right)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}{\rm N}(0,\sigma^{2})\qquad\mbox{as }n\to\infty,$

where

$\sigma^{2}=(\boldsymbol{1},0,0)B\tilde{\Sigma}(\tilde{\tau})B^{\top}(\boldsymbol{1},0,0)^{\top}.$
(5.28)

Now $\nabla\varphi(\tilde{\boldsymbol{w}}(\tilde{\tau}))=(\boldsymbol{0},1,0)$
and, using (5.6),
$(\boldsymbol{1},0,0)[\tilde{F}(\tilde{\boldsymbol{w}}(\tilde{\tau}))]^{\top}=-\tilde{x}_{E}(\tilde{\tau})$,
so $(\boldsymbol{1},0,0)B=(\boldsymbol{1},b(\tilde{\tau}),0)$, where

$b(\tilde{\tau})=a(\tilde{\tau})^{-1}\tilde{x}_{E}(\tilde{\tau}).$ (5.29)

Further, it follows from (5.10)-(5.12) and (5.24) that

$\left[(\boldsymbol{1},0,0)\tilde{\Phi}(\tilde{\tau},u)\right]_{i}={\rm
e}^{-i(\tilde{\tau}-u)}\qquad(i=0,1\dots,d_{\rm max})$

and

$\left[(\boldsymbol{0},1,0)\tilde{\Phi}(\tilde{\tau},u)\right]_{i}=\begin{cases}i\left[{\rm
e}^{-2(\tilde{\tau}-u)}-q_{I}{\rm e}^{-(\tilde{\tau}-u)}-p_{I}{\rm
e}^{-i(\tilde{\tau}-u)}\right]&\text{ if }i=0,1,\dots,d_{\rm max},\\\ {\rm
e}^{-2(\tilde{\tau}-u)}&\text{ if }i=d_{\rm max}+1,\\\ -{\rm
e}^{-(\tilde{\tau}-u)}(1-{\rm e}^{-(\tilde{\tau}-u)})&\text{ if }i=d_{\rm
max}+2.\end{cases}$ (5.30)

Further, let

$\displaystyle\boldsymbol{c}(\tilde{\tau},u)$ $\displaystyle=$
$\displaystyle(\boldsymbol{1},0,0)B\tilde{\Phi}(\tilde{\tau},u)$ (5.31)
$\displaystyle=$
$\displaystyle\left(\boldsymbol{c}_{S}(\tilde{\tau},u),c_{I}(\tilde{\tau},u),c_{R}(\tilde{\tau},u)\right),$

where

$\displaystyle\boldsymbol{c}_{S}(\tilde{\tau},u)$
$\displaystyle=\left(c_{0}(\tilde{\tau},u),c_{1}(\tilde{\tau},u),\dots,c_{d_{\rm
max}}(\tilde{\tau},u)\right),$ (5.32) $\displaystyle c_{I}(\tilde{\tau},u)$
$\displaystyle=b(\tilde{\tau}){\rm e}^{-2(\tilde{\tau}-u)},$ (5.33)
$\displaystyle c_{R}(\tilde{\tau},u)$ $\displaystyle=-b(\tilde{\tau}){\rm
e}^{-(\tilde{\tau}-u)}(1-{\rm e}^{-(\tilde{\tau}-u)})$ (5.34)

and, for $i=0,1,\dots,d_{\rm max}$,

$c_{i}(\tilde{\tau},u)={\rm e}^{-i(\tilde{\tau}-u)}+b(\tilde{\tau})i\left[{\rm
e}^{-2(\tilde{\tau}-u)}-p_{I}{\rm e}^{-i(\tilde{\tau}-u)}-q_{I}{\rm
e}^{-(\tilde{\tau}-u)}\right].$ (5.35)

Noting that $\boldsymbol{c}(\tilde{\tau},u)\boldsymbol{l}^{\top}$ is a scalar,
it follows from (5.25) and (5.28) that

$\sigma^{2}=\boldsymbol{c}(\tilde{\tau},0)\Sigma_{0}\boldsymbol{c}(\tilde{\tau},0)^{\top}+\sum_{\boldsymbol{l}\in\Delta}\int_{0}^{\tilde{\tau}}\left(\boldsymbol{c}(\tilde{\tau},u)\boldsymbol{l}^{\top}\right)^{2}\tilde{\beta}_{\boldsymbol{l}}(\tilde{\boldsymbol{w}}(u))\,{\rm
d}u.$ (5.36)

The asymptotic variances, $\sigma^{2}_{\rm MR}$ and $\sigma^{2}_{\rm NSW}$ in
Theorem 2.1 can be obtained by substituting $\Sigma_{0}=\Sigma^{\rm MR}_{0}$
and $\Sigma_{0}=\Sigma^{\rm NSW}_{0}$, respectively, in (5.36) and using
(5.31)-(5.35), together with (5.3) and (5.10)-(5.12), to evaluate the second
term on the right-hand side of (5.36). The details are lengthy and may be
found in Appendix A.3.

###  5.4 Proof of Theorem 2.2

We prove Theorem 2.2 for the epidemic on an NSW random graph. The proof for
the epidemic on an MR random graph is similar but simpler, as there is no
randomness in the degrees of individuals, and is thus omitted. The proof
proceeds in two stages. First, in Section 5.4.1, we couple the early stages of
the epidemic $\tilde{\mathcal{E}}^{(n)}$, defined in Section 5.1, to a two-
type branching process $\tilde{\mathcal{B}}^{(n)}$ which assumes that all
infective half-edges in $\tilde{\mathcal{E}}^{(n)}$ are paired with
susceptible half-edges. The branching processes $\tilde{\mathcal{B}}^{(n)}$
$(n=1,2,\dots)$ are coupled to a limiting branching process
$\tilde{\mathcal{B}}$. The couplings and standard properties of the limiting
branching process $\tilde{\mathcal{B}}$ show that, with probability tending to
$1$ as $n\to\infty$, a major outbreak occurs if and only if the branching
process $\tilde{\mathcal{B}}$ does not go extinct, and yield weak convergence
results concerning the composition of the population in
$\tilde{\mathcal{E}}^{(n)}$ when the number of infective half-edges first
reaches $\log n$ in the event of a major outbreak (see Theorem 5.1). Then, in
Section 5.4.2, we use the random time-scale transformation introduced in
Section 5.2 to determine the asymptotic distribution of the final size of a
major outbreak. The argument proceeds as in the proof of Theorem 2.1 but the
equation defining $\tilde{\tau}$ now has a solution at $0$ and one at
$\tilde{\tau}>0$ (see the discussion following (5.49)) and a lower bounding
branching process for the epidemic $\tilde{\mathcal{E}}^{(n)}$ is used to show
that $\tilde{\tau}>0$ is the relevant asymptotic hitting time.

For ease of presentation we assume that, for $n=1,2,\dots$, there is one
initial infective in $\tilde{\mathcal{E}}^{(n)}$ (i.e. that $a=1$), who is
chosen by sampling a half-edge uniformly at random from all
$D^{(n)}_{1}+D^{(n)}_{2}+\dots+D^{(n)}_{n}$ half-edges and infecting the
individual who owns that half-edge. The proofs are easily extended to $a>1$
and other ways of choosing the initial infective(s) but the details are more
complicated.

####  5.4.1 Coupling of epidemic and branching processes

Let $(\Omega,\mathcal{F},{\rm P})$ be a probability space on which are defined
the following independent sets of random varaibles:

  1. (i)

$D_{1},D_{2},\dots$ i.i.d. $\sim D$;

  2. (ii)

$U_{0},U_{1},\dots$ i.i.d. $\sim{\rm U}(0,1)$;

  3. (iii)

$L_{1},L_{2},\dots$ i.i.d. $\sim{\rm Exp}(1)$;

  4. (iv)

$Y_{01},Y_{02},\dots Y_{0d_{\rm max}}$, where $Y_{0i}\sim{\rm Bin}(i,1-{\rm
e}^{-\lambda I})$;

  5. (v)

for $i=1,2,\dots,d_{\rm max}$, $Y_{i1},Y_{i2},\dots$ i.i.d. $\sim{\rm
Bin}(i-1,1-{\rm e}^{-\lambda I})$.

For $n=1,2,\dots$, let $p^{(n)}_{i}=n^{-1}\sum_{k=1}^{n}1_{\\{D_{k}=n\\}}$
$(i=0,1,\dots,d_{\rm max})$ and
$\tilde{p}^{(n)}_{i}=ip^{(n)}_{i}/\mu_{D}^{(n)}$ $(i=1,2,\dots,d_{\rm max})$,
where $\mu_{D}^{(n)}=n^{-1}\sum_{k=1}^{n}D_{k}$. Note that by the strong law
of large numbers $\mu_{D}^{(n)}>0$ (and $\tilde{p}^{(n)}_{i}$ is well defined)
for all sufficiently large $n$ almost surely. Let
$\tilde{c}^{(n)}_{i}=\sum_{j=1}^{i}\tilde{p}^{(n)}_{j}$ $(i=1,2,\dots,d_{\rm
max})$. For $x\in(0,1)$, let
$\tilde{d}^{(n)}(x)=\min\\{i:x\leq\tilde{c}^{(n)}_{i}\\}$. Similarly, let
$\tilde{p}_{i}=ip_{i}/\mu_{D}$ and $\tilde{c}_{i}=\sum_{j=1}^{i}\tilde{p}_{j}$
$(i=1,2,\dots,d_{\rm max})$, and let
$\tilde{d}(x)=\min\\{i:x\leq\tilde{c}_{i}\\}$ $(0<x<1)$.

For $n=1,2,\dots$, construct on $(\Omega,\mathcal{F},{\rm P})$ a realisation
of a two-type continuous-time Markov branching process
$\tilde{\mathcal{B}}^{(n)}$, which approximates the process of infected and
recovered half-edges in the epidemic $\tilde{\mathcal{E}}^{(n)}$, as follows.
The types are denoted $I$ and $R$ depending on whether the individual
corresponds to an infective or recovered half-edge. Only type-$I$ individuals
have offspring and they do so at their moment of death. Type-$R$ individuals
live forever. For $t\geq 0$, let $\hat{Y}^{(n)}_{E}(t)$ and
$\hat{Z}^{(n)}_{E}(t)$ denote respectively the number of type-$I$ and type-$R$
individuals alive in $\tilde{\mathcal{B}}^{(n)}$ at time $t$. The initial
state $(\hat{Y}^{(n)}_{E}(0),\hat{Z}^{(n)}_{E}(0))$ is determined as follows.
Let $d^{(n)}_{0}=\tilde{d}^{(n)}(U_{0})$, which corresponds to the initial
infective in $\tilde{\mathcal{E}}^{(n)}$ having degree $d^{(n)}_{0}$. If
$d^{(n)}_{0}=0$ then $(\hat{Y}^{(n)}_{E}(0),\hat{Z}^{(n)}_{E}(0))=(0,0)$ and
$\tilde{\mathcal{B}}^{(n)}$ goes extinct immediately. Alternatively, if
$d^{(n)}_{0}>0$ then
$(\hat{Y}^{(n)}_{E}(0),\hat{Z}^{(n)}_{E}(0))=(Y_{0d^{(n)}_{0}},d_{0}-Y_{0d^{(n)}_{0}})$.
In that case, for $k=1,2,\dots$, the $k$th type-$I$ individual born in
$\tilde{\mathcal{B}}^{(n)}$ (including the initial individuals) has degree
$d^{(n)}_{k}=\tilde{d}^{(n)}(U_{k})$ and lives until age $L_{k}$, when it
dies. Denote this individual by $i^{*}$. Suppose that $d^{(n)}_{k}=i$ and
$i^{*}$ is the $l$th degree-$i$ individual (excluding the initial individuals)
born in $\tilde{\mathcal{B}}^{(n)}$. Then, when $i^{*}$ dies, it leaves
$Y_{il}$ type-$I$ and $i-1-Y_{il}$ type-$R$ offspring. Of course, reproduction
stops in $\tilde{\mathcal{B}}^{(n)}$ if $\hat{Y}^{(n)}_{E}(t)=0$. Construct
also on $(\Omega,\mathcal{F},{\rm P})$ a realisation of a two-type continuous-
time branching process $\tilde{\mathcal{B}}$, defined analogously to
$\tilde{\mathcal{B}}^{(n)}$ but using the function $\tilde{d}$ instead of
$\tilde{d}^{(n)}$. For $t\geq 0$, let $Y_{E}(t)$ and $Z_{E}(t)$ denote
respectively the numbers of type-$I$ and type-$R$ individuals alive in
$\tilde{\mathcal{B}}$ at time $t$.

For $n=1,2,\dots$, construct on $(\Omega,\mathcal{F},{\rm P})$ a realisation
of the epidemic $\tilde{\mathcal{E}}^{(n)}$, defined in Section 5.1, as
follows. Give the $n$ individuals in $\tilde{\mathcal{E}}^{(n)}$ the labels
$1,2,\dots,n$ in increasing order of degree. Now label the $n\mu_{D}^{(n)}$
half-edges $1,2,\dots,n\mu_{D}^{(n)}$, starting with the half-edges (if any)
attached to individual $1$, then the half-edges (if any) attached to
individual $2$ and so on. Thus half-edges attached to the same individual have
consecutive labels. As in Section 5.1, for $t\geq 0$, let $Y^{(n)}_{E}(t)$ and
$Z^{(n)}_{E}(t)$ denote respectively the number of infective and recovered
half-edges in $\tilde{\mathcal{E}}^{(n)}$ at time $t$. The initial infective
in $\tilde{\mathcal{E}}^{(n)}$ is the individual who owns the half-edge
labelled $\left\lfloor{n\mu_{D}^{(n)}U_{0}}\right\rfloor+1$. Note that this
individual has degree $d^{(n)}_{0}=\tilde{d}^{(n)}(U_{0})$. If $d^{(n)}_{0}=0$
then the epidemic stops immediately. Alternatively, if $d^{(n)}_{0}>0$ then
the initial infective infects along $k^{(n)}_{0}=Y_{0d^{(n)}_{0}}$ of its
half-edges with its remaining $d^{(n)}_{0}-k^{(n)}_{0}$ half-edges becoming
recovered half-edges. The epidemic stops if $k^{(n)}_{0}=0$, otherwise the
initial infective transmits infection along its infective half-edges at times
$L_{1},L_{2},\dots,L_{k_{0}}$.

For $k=1,2,\dots,$ let
$l_{k}^{(n)}=\left\lfloor{n\mu_{D}^{(n)}U_{k}}\right\rfloor+1$ and note that
the half-edge having label $l_{k}^{(n)}$ is attached to an individual having
degree $d^{(n)}_{k}=\tilde{d}^{(n)}(U_{k})$. When infection is transmitted
along a half-edge that half-edge, $l_{*}$ say, is attempted to be paired with
the half-edge having label $l_{k}^{(n)}$, where $k$ is the number of the
$U_{0},U_{1},\dots$ that have been used already in the construction of
$\tilde{\mathcal{E}}^{(n)}$. (Thus, for example, the first half-edge emanating
from the initial infective is attempted to be paired with the half-edge having
label $l_{1}^{(n)}$.) If the half-edge $l_{k}^{(n)}$ has already been paired
or $l_{k}^{(n)}=l_{*}$ then the attempt fails and $l_{*}$ is attempted to be
paired with the half-edge $l_{k+1}^{(n)}$, and so on until a valid pairing is
obtained and an edge is formed. Suppose that a valid pairing is made with the
half-edge having label $l_{V}$. Let $i^{*}$ be the individual that owns the
half-edge $l_{V}$ and $i$ be the degree of $i^{*}$. If $i^{*}$ is susceptible
then it becomes an infective, otherwise nothing happens apart from the
formation of the edge. Suppose that $i^{*}$ is susceptible and is the $l$th
degree-$i$ susceptible to be infected in $\tilde{\mathcal{E}}^{(n)}$,
excluding the initial infective. Then $i^{*}$ infects along $Y_{il}$ of its
half-edges and its remaining $i-1-Y_{il}$ half-edges become recovered half-
edges. (When $i^{*}$ was infected one of its $i$ half-edges was paired to its
infector.) Let $k_{1}=Y_{il}$. The times of these $k_{1}$, infections,
relative to the infection time of $i^{*}$, are given by
$L_{k^{*}+1},L_{k^{*}+2},\dots,L_{k^{*}+k_{1}}$, where $k^{*}$ is the number
of infective half-edges created in $\tilde{\mathcal{E}}^{(n)}$ prior to the
infection of $i^{*}$. The epidemic terminates when $Y^{(n)}_{E}(t)=0$, or when
$Y^{(n)}_{E}(t)=1$ and all other half-edges have been paired.

As in Section 5.1, for $t\geq 0$, let
$\boldsymbol{X}^{(n)}(t)=(X^{(n)}_{0}(t),X^{(n)}_{1}(t),\dots,X^{(n)}_{d_{\rm
max}}(t))$, where $X^{(n)}_{i}(t)$ is the number of degree-$i$ susceptible
individuals at time $t$ in $\tilde{\mathcal{E}}^{(n)}$. For $n=1,2,\dots$, let
$\hat{\tau}^{(n)}=\inf\\{t\geq 0:Y^{(n)}_{E}(t)\geq\log n\\}$, where
$\hat{\tau}^{(n)}=\infty$ if $Y^{(n)}_{E}(t)<\log n$ for all $t\geq 0$. Let
$A_{\rm ext}=\\{\omega\in\Omega:\lim_{t\to\infty}Y_{E}(t)=0\\}$ denote the set
on which type-$I$ individuals become extinct in the branching process
$\tilde{\mathcal{B}}$ and let $\alpha$ denote the Malthusian parameter of
$\tilde{\mathcal{B}}$. Then $\alpha$ is given by the unique real solution of
the equation

$\int_{0}^{\infty}{\rm e}^{-\alpha t}\mu_{\tilde{D}-1}p_{I}{\rm e}^{-t}\,{\rm
d}t=1,$

so $\alpha=R_{0}-1$, where $R_{0}$ is defined at (2.8).

###### Theorem 5.1

  1. (a)

$\lim_{n\to\infty}{\rm P}\left(\hat{\tau}^{(n)}=\infty|A_{\rm ext}\right)=1$.

  2. (b)

Suppose that $R_{0}>1$, so ${\rm P}(A_{\rm ext}^{C})>0$. Then, as
$n\to\infty$,

    1. (i)

${\rm P}\left(\hat{\tau}^{(n)}<\infty|A_{\rm ext}^{C}\right)\to 1$;

    2. (ii)

$\frac{1}{\log n}Y^{(n)}_{E}(\hat{\tau}^{(n)})\left|\right.A_{\rm
ext}^{C}\stackrel{{\scriptstyle p}}{{\longrightarrow}}1$;

    3. (iii)

$\frac{1}{\log n}Z^{(n)}_{E}(\hat{\tau}^{(n)})\left|\right.A_{\rm
ext}^{C}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\alpha^{-1}q_{I}\mu_{\tilde{D}-1}$;

    4. (iv)

$\sqrt{n}\left(n^{-1}\boldsymbol{X}^{(n)}(\hat{\tau}^{(n)})-\boldsymbol{p}\right)|A_{\rm
ext}^{C}\stackrel{{\scriptstyle D}}{{\longrightarrow}}{\rm
N}(\boldsymbol{0},\Sigma_{XX})$, where $\Sigma_{XX}$ is defined at (5.18).

  * Proof.

The key observations underlying the proof are that (i) the processes
$\\{(\hat{Y}^{(n)}_{E}(t),\hat{Z}^{(n)}_{E}(t)):t\geq 0\\}$ and
$\\{(Y^{(n)}_{E}(t),Z^{(n)}_{E}(t)):t\geq 0\\}$ coincide up until at least the
first time that an attempt is made to pair a half-edge with a half-edge
belonging to an individual previously used in the construction of
$\tilde{\mathcal{E}}^{(n)}$ and (ii) the branching processes
$\tilde{\mathcal{B}}^{(n)}$ and $\tilde{\mathcal{B}}$ coincide up until the
first time that $\tilde{d}^{(n)}(U_{k})\neq\tilde{d}(U_{k})$. Thus we show
that the probability that the processes
$\\{(\hat{Y}^{(n)}_{E}(t),\hat{Z}^{(n)}_{E}(t)):0\leq
t\leq\hat{\tau}^{(n)}\\}$ and $\\{(Y_{E}(t),Z_{E}(t)):0\leq
t\leq\hat{\tau}^{(n)}\\}$ coincide converges to $1$ as $n\to\infty$. The
theorem then follows using standard results concerning the asymptotic
behaviour of the branching process $\tilde{\mathcal{B}}$.

For $n=1,2,\dots$ and $k=1,2,\dots,n\mu_{D}^{(n)}$, let $\mathcal{C}^{(n)}(k)$
be the set of half-edges attached to the individual owning half-edge $k$ in
$\tilde{\mathcal{E}}^{(n)}$. Thus $k\in\mathcal{C}^{(n)}(k)$ and
$|\mathcal{C}^{(n)}(k)|\leq d_{\rm max}$. For $n=1,2,\dots$, let

$M^{(n)}=\min\left\\{k\geq
1:l_{k}^{(n)}\in\bigcup_{i=0}^{k-1}\mathcal{C}^{(n)}(l_{i}^{(n)})\right\\}$

be the number of pairings required in $\tilde{\mathcal{E}}^{(n)}$ until an
attempt is made to pair a half-edge with a half-edge belonging to a previoulsy
used individual. Let $\boldsymbol{D}=(D_{1},D_{2},\dots)$. Then, for
$m=1,2,\dots$,

$\displaystyle{\rm P}(M^{(n)}\leq m|\boldsymbol{D})$ $\displaystyle=$
$\displaystyle{\rm
P}\left(\bigcup_{k=1}^{m}\left\\{l_{k}^{(n)}\in\bigcup_{i=0}^{k-1}\mathcal{C}^{(n)}(l_{i}^{(n)})\right\\}|\boldsymbol{D}\right)$
$\displaystyle\leq$ $\displaystyle\frac{m(m+1)}{2}\frac{d_{\rm
max}}{n\mu_{D}^{(n)}}.$

Now $\mu_{D}^{(n)}\stackrel{{\scriptstyle a.s.}}{{\longrightarrow}}\mu_{D}$ as
$n\to\infty$ by the strong law of large numbers. Hence, for any
$\gamma\in(0,1/2)$, ${\rm P}(M^{(n)}\leq
n^{\gamma}|\boldsymbol{D})\stackrel{{\scriptstyle p}}{{\longrightarrow}}0$ as
$n\to\infty$. Thus,

$\lim_{n\to\infty}{\rm P}(M^{(n)}\leq n^{\gamma})=\lim_{n\to\infty}{\rm
E}[{\rm P}(M^{(n)}\leq n^{\gamma}|\boldsymbol{D})]=0,$ (5.37)

as ${\rm P}(M^{(n)}\leq n^{\gamma}|\boldsymbol{D})$ $(n=1,2,\dots)$ is
uniformly bounded.

Turning to the coupling of $\tilde{\mathcal{B}}^{(n)}$ and
$\tilde{\mathcal{B}}$, for $n=\,2,\dots$,

$\displaystyle\max_{i=1,2,\dots,d_{\rm
max}-1}|\tilde{c}^{(n)}_{i}-\tilde{c}_{i}|$ $\displaystyle\leq$
$\displaystyle\frac{1}{2}\sum_{i=1}^{d_{\rm
max}}|\tilde{p}^{(n)}_{i}-\tilde{p}_{i}|$ $\displaystyle\leq$
$\displaystyle\frac{1}{\mu_{D}}\sum_{i=1}^{d_{\rm max}}i|p^{(n)}_{i}-p_{i}|.$

The first inequality in (5.4.1) follows from the total variation distance
between the distributions $\tilde{p}^{(n)}_{i}$ $(i=1,2,\dots,d_{\rm max})$
and $\tilde{p}_{i}$ $(i=1,2,\dots,d_{\rm max})$, and the second inequality
follows using the triangle inequality and elementary manipulation ([Ball and
Neal(2017)], equation (3.2)). For $n=1,2,\dots$, let
$\hat{M}^{(n)}=\min\\{k\geq 0:\tilde{d}^{(n)}(U_{k})\neq\tilde{d}(U_{k})\\}$.
Then, for $k=0,1,\dots$,

${\rm P}\left(\hat{M}_{n}\leq k\right)\leq\frac{(k+1)(d_{\rm
max}-1)}{\mu_{D}}{\rm E}\left[\sum_{i=1}^{d_{\rm
max}}i|p^{(n)}_{i}-p_{i}|\right].$

Now, for $i=1,2,d_{\rm max}$, as $np^{(n)}_{i}\sim{\rm Bin}(n,p_{i})$,

${\rm E}\left[|p^{(n)}_{i}-p_{i}|\right]\leq\sqrt{{\rm
var}(p^{(n)}_{i})}=\frac{1}{\sqrt{n}}\sqrt{p_{k}(1-p_{k})}.$

Thus, for any $\gamma\in(0,1/2)$,

$\lim_{n\to\infty}{\rm P}\left(\hat{M}^{(n)}\leq n^{\gamma}\right)=0.$ (5.39)

For $t\geq 0$, let $T_{E}(t)$ be the total number of individuals (of either
type) alive in the branching process $\tilde{\mathcal{B}}$ during $[0,t]$ and
let $T_{E}(\infty)=\lim_{t\to\infty}T_{E}(t)$. Thus
$T_{E}(\infty,\omega)<\infty$ if and only if $\omega\in A_{\rm ext}$. For
$n=1,2,\dots$ and $t\geq 0$, let

$A^{(n)}_{\rm
couple}(t)=\left\\{\omega\in\Omega:\left(Y^{(n)}_{E}(u),Z^{(n)}_{E}(u)\right)=\left(Y_{E}(u),Z_{E}(u)\right)\mbox{
for all }u\in[0,t]\right\\}.$

Let $\tau=\min\\{t\geq 0:Y_{E}(t)=0\\}$, where $\tau(\omega)=\infty$ if
$\omega\in A_{\rm ext}^{C}$. Then (5.37) and (5.39) imply that, for
$k=0,1,\dots$,

$\lim_{n\to\infty}{\rm P}\left(A^{(n)}_{\rm
couple}(\tau)\cap\left\\{T_{E}(\infty)=k\right\\}\right)=1.$

Part (a) of the theorem follows since $A_{\rm
ext}=\bigcup_{k=0}^{\infty}\\{T_{E}(\infty)=k\\}$.

Note that although $\tilde{\mathcal{B}}$ is a two-type branching process,
since only type-$I$ individuals produce offspring, it can be treated as a
single-type branching process consisting only of type-$I$ individuals, in
which attached to each (type-$I$) individual is a random characteristic (see
[Nerman(1981)]) given by its number of type-$R$ offspring. Then it follows
from [Nerman(1981)], Theorem 5.4, that there exists a random variable $W\geq
0$, where $W(\omega)=0$ if and only if $\omega\in A_{\rm ext}$ such that, as
$t\to\infty$,

${\rm e}^{-\alpha t}Y_{E}(t)\stackrel{{\scriptstyle
a.s.}}{{\longrightarrow}}W,$ (5.40) ${\rm e}^{-\alpha
t}Z_{E}(t)\stackrel{{\scriptstyle
a.s.}}{{\longrightarrow}}\alpha^{-1}q_{I}\mu_{\tilde{D}-1}W$ (5.41)

and

${\rm e}^{-\alpha t}T_{E}(t)\stackrel{{\scriptstyle
a.s.}}{{\longrightarrow}}\alpha^{-1}\mu_{\tilde{D}-1}W.$ (5.42)

(It is easily checked that (5.40)-(5.42) hold by calculating the the
appropriate $m_{\infty}^{\phi}$ in [Nerman(1981)], Theorem 5.4, and using
Corollary 3.2 of that paper. For a heuristic argument note, for example, that
if $Y_{E}(u)\approx W{\rm e}^{\alpha u}$ $(u\geq 0)$ then
$T_{E}(t)\approx\int_{0}^{t}\mu_{\tilde{D}-1}{\rm e}^{\alpha u}W\,{\rm
d}u\approx\alpha^{-1}\mu_{\tilde{D}-1}W{\rm e}^{\alpha t}$, since type-$I$
individuals die at rate $1$ and have on average $\mu_{\tilde{D}-1}$
offspring.)

For $n=1,2,\dots$, let $\bar{\tau}^{(n)}=\inf\\{t\geq 0:Y_{E}(t)\geq\log
n\\}$, where $\bar{\tau}^{(n)}(\omega)=\infty$ if $\omega\in A_{\rm ext}$.
Then it follows from (5.40) that $\bar{\tau}^{(n)}(\omega)<\infty$ and
$\lim_{n\to\infty}\bar{\tau}^{(n)}(\omega)=\infty$, for ${\rm P}$-almost all
$\omega\in A_{\rm ext}^{C}$. Thus, for ${\rm P}$-almost all $\omega\in A_{\rm
ext}^{C}$,

$\lim_{n\to\infty}\frac{1}{\log n}Y_{E}(\bar{\tau}^{(n)})=1,$ (5.43)

since an individual has at most $d_{\rm max}$ offspring, and using
(5.40)-(5.42),

$\lim_{n\to\infty}\frac{1}{\log
n}Z_{E}(\bar{\tau}^{(n)})=\alpha^{-1}q_{I}\mu_{\tilde{D}-1}$ (5.44)

and

$\lim_{n\to\infty}\frac{1}{\log
n}T_{E}(\bar{\tau}^{(n)})=\alpha^{-1}\mu_{\tilde{D}-1}.$ (5.45)

Now (5.37), (5.39) and (5.45) imply that

$\lim_{n\to\infty}{\rm P}\left(A^{(n)}_{\rm couple}(\bar{\tau}^{(n)})|A_{\rm
ext}^{C}\right)=1,$ (5.46)

and (i)-(iii) of part (b) of the theorem follow using (5.41) and (5.43).

To prove part (b)(iv), note by the multivariate central limit theorem that

$\sqrt{n}\left(n^{-1}\boldsymbol{X}^{(n)}(0)-\boldsymbol{p}\right)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}{\rm N}(\boldsymbol{0},\Sigma_{XX})\qquad\mbox{as
}n\to\infty.$

Also, (5.45) and (5.46) imply that, for $i=0,1,\dots,d_{\rm max}$,

$\lim_{n\to\infty}{\rm
P}\left(X^{(n)}_{i}(\hat{\tau}^{(n)})-X^{(n)}_{i}(0)\leq
2\alpha^{-1}\mu_{\tilde{D}-1}\log n|A_{\rm ext}^{C}\right)=1,$

so
$\frac{1}{\sqrt{n}}\left(X^{(n)}_{i}(\hat{\tau}^{(n)})-X^{(n)}_{i}(0)\right)|A_{\rm
ext}^{C}\stackrel{{\scriptstyle p}}{{\longrightarrow}}0$ as $n\to\infty$. Part
(b)(iv) now follows using Slutsky’s theorem. $\Box$

####  5.4.2 Epidemic starting with $\left\lceil{\log n}\right\rceil$
infectives

Turning to the proof of Theorem 2.2, note that by Theorem 5.1,
$\lim_{n\to\infty}{\rm P}\left(G^{(n)}\triangle A_{\rm ext}\right)=0$, where
$\triangle$ denotes symmetric difference. Hence, using the strong Markov
property, we can determine the asymptotic distribution of $T^{(n)}_{\rm
NSW}|G^{(n)}$ by considering the random time-scale transformed process
$\\{\tilde{\boldsymbol{W}}^{(n)}(t)\\}$, defined in Section 5.2, with initial
state
$\tilde{\boldsymbol{W}}^{(n)}(0)=(\boldsymbol{X}^{(n)}(\hat{\tau}^{(n)}),Y^{(n)}_{E}(\hat{\tau}^{(n)}),Z^{(n)}_{E}(\hat{\tau}^{(n)}))$.
Thus, by Theorem 5.1(b),

$\sqrt{n}\left(n^{-1}\tilde{\boldsymbol{W}}^{(n)}(0)-(\boldsymbol{p},0,0)\right)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}{\rm N}(\boldsymbol{0},\Sigma^{\rm
NSW}_{0})\qquad\mbox{as }n\to\infty,$ (5.47)

where $\Sigma^{\rm NSW}_{0}$ is obtained by setting $\epsilon=0$ in (5.21).

By Theorem 4.1, for any $t>0$,

$\sup_{0\leq u\leq
t}|n^{-1}\tilde{Y}^{(n)}_{E}(u)-\tilde{y}_{E}(u)|\stackrel{{\scriptstyle
p}}{{\longrightarrow}}0\qquad\mbox{as }n\to\infty,$ (5.48)

where, using (5.11)

$\tilde{y}_{E}(t)=\mu_{D}\left({\rm e}^{-2t}-q_{I}{\rm
e}^{-t}\right)-p_{I}f_{D}^{(1)}({\rm e}^{-t}).$

Note that $\tilde{y}_{E}(t)=0$ if and only if

${\rm e}^{-t}-q_{I}-\mu_{D}^{-1}p_{I}f_{D}^{(1)}({\rm e}^{-t})=0$ (5.49)

and that $t=0$ is a solution of (5.49). If $R_{0}\leq 1$ then $t=0$ is the
only solution in $[0,\infty)$, but if $R_{0}>1$ then we show in Appendix A.2
that, under the conditions of Theorem 2.2, there is a unique solution in
$(0,\infty)$ which we denote by $\tilde{\tau}$.

As in Section 5.3.3, let $\tilde{\tau}^{(n)}=\inf\\{t\geq
0:\varphi(n^{-1}\tilde{\boldsymbol{W}}^{(n)}(t))\leq 0\\}$. Recall that
$R_{0}>1$ is assumed in Theorem 2.2. It follows from (5.48) that,
$\min\left(\tilde{\tau}^{(n)},|\tilde{\tau}^{(n)}-\tilde{\tau}|\right)\stackrel{{\scriptstyle
p}}{{\longrightarrow}}0$ as $n\to\infty$. We show below that there exists
$\epsilon_{0}>0$ such that

$\lim_{n\to\infty}{\rm P}(\tilde{\tau}^{(n)}<\epsilon_{0})=0.$ (5.50)

Theorem 4.3 and Corollary 4.4, as stated, cannot be applied in the present
setting as $\varphi(\tilde{\boldsymbol{w}}(0))=0$. However, in the terminology
of Theorem 4.3, the proof of [Ethier and Kurtz(1986)], Theorem 11.4.1, extends
easily to the situation when
$\tau^{(n)}=\inf\\{t\geq\epsilon_{0}:\varphi(n^{-1}\boldsymbol{X}^{(n)}(t))\leq
0\\}$ and $\tau=\inf\\{t\geq\epsilon_{0}:\varphi(\boldsymbol{x}(t))\leq 0\\}$,
for fixed $\epsilon_{0}>0$. Hence, so do Theorem 4.3 and Corollary 4.4.

Let

$\tilde{T}^{(n)}=\sum_{i=0}^{d_{\rm
max}}\left(\tilde{X}^{(n)}_{i}(0)-\tilde{X}^{(n)}_{i}(\tilde{\tau}^{(n)})\right).$
(5.51)

Then using (5.50) and the above extension of Corollary 4.4 (we show in
Appendix A.2 that
$a(\tilde{\tau})=\nabla\varphi(\boldsymbol{w}(\tilde{\tau}))\cdot\tilde{F}(\boldsymbol{w}(\tilde{\tau}))<0$),
setting $\epsilon=0$ in the proof of Theorem 2.1 yields

$\sqrt{n}\left(n^{-1}\tilde{T}^{(n)}-\rho\right)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}{\rm N}(0,\tilde{\sigma}^{2}_{\rm NSW})\qquad\mbox{as
}n\to\infty.$

Theorem 2.2 follows using Slutsky’s theorem since, by (5.45),

$\frac{1}{\sqrt{n}}\left(T^{(n)}_{\rm
NSW}-\tilde{T}^{(n)}\right)\stackrel{{\scriptstyle
p}}{{\longrightarrow}}0\qquad\mbox{as }n\to\infty.$

To complete the proof of Theorem 2.2 we use a device, first introduced by
[Whittle(1955)] in the setting of a homegeneously mixing epidemic, to show
that there exists $\epsilon_{0}>0$ such that (5.50) holds. For $n=1,2,\dots$,
construct on $(\Omega,\mathcal{F},{\rm P})$ an epidemic
$\check{\mathcal{E}}^{(n)}$ analogously to $\tilde{\mathcal{E}}^{(n)}$ except
having initial state given in an obvious notation by

$\left(\check{\boldsymbol{X}}^{(n)}(0),\check{Y}^{(n)}_{E}(0),\check{Z}^{(n)}_{E}(0)\right)=\left(\boldsymbol{X}^{(n)}(\hat{\tau}^{(n)}),Y^{(n)}_{E}(\hat{\tau}^{(n)}),Z^{(n)}_{E}(\hat{\tau}^{(n)})\right).$

Let $\check{\mathcal{B}}^{(n)}$ be the corresponding branching process which
assumes all attempted pairing in $\check{\mathcal{E}}^{(n)}$ are valid. Let

$\check{T}^{(n)}=\sum_{i=0}^{d_{\rm
max}}\left(\check{T}^{(n)}_{i}(0)-\check{T}^{(n)}_{i}(\infty)\right)$

be the final size of $\check{\mathcal{E}}^{(n)}$. Fix $\delta\in(0,1)$. Then,
provided $\check{T}^{(n)}\leq n\delta$, for a given pairing, the probability
that an invalid pairing is attempted is at most $p^{(n)}(\delta)=\min(\delta
nd_{\rm max}/(n\mu_{D}^{(n)}),1)=\min(\delta d_{\rm max}/\mu_{D}^{(n)},1)$. It
follows that

${\rm P}\left(\check{T}^{(n)}\geq n\delta\right)\geq{\rm
P}\left(\check{T}^{(n)}_{B}(\delta)\geq n\delta\right),$ (5.52)

where $\check{T}^{(n)}_{B}(\delta)$ is the total number of type-$I$
individuals (excluding the initial individuals) ever alive in the branching
process $\check{\mathcal{B}}^{(n)}(\delta)$ that is obtained from
$\check{\mathcal{B}}^{(n)}$ by aborting each type-$I$ individual independently
with probability $p^{(n)}(\delta)$.

Let $R^{(n)}_{0}(\delta)$ be the mean number of type-$I$ individuals spawned
by a typical type-$I$ individual in $\check{\mathcal{B}}^{(n)}(\delta)$ and
let $\pi^{(n)}(\delta)$ be the extinction probability for type-$I$ individuals
in $\check{\mathcal{B}}^{(n)}(\delta)$ if the process were to start with one
type-$I$ individual. As $n\to\infty$, $p^{(n)}_{i}\stackrel{{\scriptstyle
a.s.}}{{\longrightarrow}}p_{i}$ ($i=0,1,\dots,d_{\rm max}$), by the strong law
of large numbers, so $\tilde{p}^{(n)}_{i}\stackrel{{\scriptstyle
a.s.}}{{\longrightarrow}}\tilde{p}_{i}$ ($i=1,2,\dots,d_{\rm max}$) and
$\mu^{(n)}_{\tilde{D}-1}\stackrel{{\scriptstyle
a.s.}}{{\longrightarrow}}\mu_{\tilde{D}-1}$, where
$\mu^{(n)}_{\tilde{D}-1}=\sum_{i=1}^{d_{\rm max}}(i-1)\tilde{p}^{(n)}_{i}$.
Thus,

$R^{(n)}_{0}(\delta)=(1-p^{(n)}(\delta))\mu^{(n)}_{\tilde{D}-1}p_{I}\stackrel{{\scriptstyle
a.s.}}{{\longrightarrow}}(1-p(\delta))R_{0}\qquad\mbox{as }n\to\infty,$

where $p(\delta)=\min(\delta d_{\rm max}/\mu_{D},1)$. Further, as
$n\to\infty$, the offspring distribution of
$\check{\mathcal{B}}^{(n)}(\delta)$ converges almost surely to that of
$\tilde{\mathcal{B}}(\delta)$, where $\tilde{\mathcal{B}}(\delta)$ is the
branching process obtained from $\tilde{\mathcal{B}}$ by aborting each
type-$I$ individual independently with probability $p(\delta)$. Thus, by a
simple extension of [Britton et al.(2007)], Lemma 4.1,
$\pi^{(n)}(\delta)\stackrel{{\scriptstyle
a.s.}}{{\longrightarrow}}\pi(\delta)$ as $n\to\infty$, where $\pi(\delta)$ is
the extinction probability for type-$I$ individuals in
$\tilde{\mathcal{B}}(\delta)$ starting from one type-$I$ individual.

Recall that $R_{0}>1$. Thus $\delta\in(0,1)$ can be chosen sufficently small
such that $(1-p(\delta))R_{0}>1$, whence $\pi(\delta)<1$.

Now, $\check{Y}^{(n)}_{E}(0)\geq\log n$ by construction, so using (5.52), for
such $\delta$,

$\lim_{n\to\infty}{\rm P}\left(\check{T}^{(n)}\geq
n\delta|\boldsymbol{D}\right)\geq
1-\lim_{n\to\infty}\left(\pi^{(n)}\right)^{\log n}=1,$

for ${\rm P}$-almost all $\omega\in A_{\rm ext}^{C}$. Hence

$\lim_{n\to\infty}{\rm P}\left(\check{T}^{(n)}\geq n\delta|A_{\rm
ext}^{C}\right)=1.$ (5.53)

Finally, for $t\geq 0$, let $\tilde{x}(t)=\sum_{i=0}^{d_{\rm
max}}\tilde{x}_{i}(t)=f_{D}({\rm e}^{-t})$, using (5.10) with
$\tilde{x}_{i}(0)=p_{i}$ $(i=0,1,\dots,d_{\rm max})$. By construction,
$\boldsymbol{W}^{(n)}(\tau^{(n)})=\tilde{\boldsymbol{W}}^{(n)}(\tilde{\tau}^{(n)})$,
so (5.53) and Theorem 4.1 applied to the process
$\\{\tilde{\boldsymbol{W}}^{(n)}(t)\\}$, with initial state satisfying (5.47),
imply that (5.50) holds with $\epsilon_{0}$ being given by the unique solution
in $(0,\infty)$ of $1-f_{D}({\rm e}^{-\epsilon_{0}})=\delta$.

###  5.5 Proof of Theorem 2.3

The proof of Theorem 2.3 for the epidemic on the NSW random graph parallels in
an obvious fashion the argument in Section 5.4.2 above without the need to
condition on $A_{\rm ext}^{C}$, so the details are omitted. Again, the proof
for the epidemic on the MR random graph is similar but simpler.

###  5.6 Proof of Theorem 2.7

Consider first bond percolation and note that when $I\equiv 1$ and
$\lambda=-\log(1-\pi)$ then in the directed random graph
$\tilde{\mathcal{G}}^{(n)}$ defined at the start of Section 5.1, all the
possible directed edges in $\tilde{\mathcal{G}}^{(n)}$ (i.e. between pairs of
neighbours in $\mathcal{G}^{(n)}$) are present independently, each with
probability $\pi$. Further, when constructing the final outcome of the
corresponding epidemic $\mathcal{E}^{(n)}$ using $\tilde{\mathcal{G}}^{(n)}$,
for any pair $(i,j)$ of distinct individuals use is made of at most one of the
directed edges $i\to j$ and $j\to i$ (if $i$ infects $j$, whether or not $j$
tries to infect $i$ is immaterial). Thus, in this situation
$\tilde{\mathcal{G}}^{(n)}$ can be replaced by $\mathcal{G}^{(n)}_{\rm bond}$,
obtained using bond percolation on $\mathcal{G}^{(n)}$, and in
$\mathcal{E}^{(n)}$, an initial susceptible is ultimately infected if and only
if in $\mathcal{G}^{(n)}_{\rm bond}$ there is a chain of edges connecting it
to an initial infective.

Suppose that there is one initial infective. Then the final size of the
epidemic (including the initial infective) is given by the size of the
connected component of $\mathcal{G}^{(n)}_{\rm bond}$ that contains the
initial infective. With probability tending to $1$ as $n\to\infty$, a major
outbreak occurs in $\mathcal{E}^{(n)}$ if and only if the initial infective
belongs to the largest connected component of $\mathcal{G}^{(n)}_{\rm bond}$,
and the final size of a major epidemic is given by the size $C^{(n)}$ of that
connected component. Further $C^{(n)}$ has the same asymptotic distribution as
$T^{(n)}|G^{(n)}$. Thus (2.11) and (2.12) follow immediately on setting
$I\equiv 1$ and $\lambda=-\log(1-\pi)$ in Theorem 2.2. (Note that $p_{I}=\pi$
and, since $I$ is constant, $q_{I}^{(2)}=q_{I}^{2}$.)

Turning to site percolation, consider the epidemic $\mathcal{E}^{(n)}$ with
${\rm P}(I=\infty)=\pi=1-{\rm P}(I=0)$, so each infective infects all of its
neighbours with probability $\pi$ and none of them otherwise. Suppose that
there is one initial infective, $i^{*}$ say, and let
$\mathcal{T}^{(n)}=\\{j\in\mathcal{N}^{(n)}\setminus\\{i^{*}\\}:i^{*}\leadsto
j\\}$ using the directed random graph $\tilde{\mathcal{G}}^{(n)}$. Then
$\mathcal{T}^{(n)}\cup\\{i^{*}\\}$ differs from the connected component
containing $i^{*}$ in $\mathcal{G}^{(n)}_{\rm site}$ (site percolation on
$\mathcal{G}^{(n)}$) in that $\mathcal{T}^{(n)}$ also includes infected
indiviudals having $I=0$, which are deleted in $\mathcal{G}^{(n)}_{\rm site}$.
Thus, to obtain a central limit theorem for $C^{(n)}$, we need one for
$V^{(n)}=|\\{j\in\mathcal{N}^{(n)}\setminus\\{i^{*}\\}:I_{j}=\infty\mbox{ and
}i^{*}\leadsto j\\}|$, the final size of $\mathcal{E}^{(n)}$ counting only
individuals with $I=\infty$. This can be obtained by augmenting the process
$\\{\boldsymbol{W}^{(n)}(t)\\}$ as we now outline.

Let
$\\{\hat{\boldsymbol{W}}^{(n)}(t)\\}=\\{\hat{\boldsymbol{W}}^{(n)}(t):t\geq
0\\}$, where

$\hat{\boldsymbol{W}}^{(n)}(t)=(\boldsymbol{W}^{(n)}(t),V^{(n)}(t))=(\boldsymbol{X}^{(n)}(t),Y^{(n)}_{E}(t),Z^{(n)}_{E}(t),V^{(n)}(t))$

and $V^{(n)}(t)$ is the total number of initial susceptibles that are infected
during $(0,t]$ and have $I=\infty$. A typical element of $\hat{H}^{(n)}$, the
state space of $\\{\hat{\boldsymbol{W}}^{(n)}(t)\\}$, is now
$\hat{\boldsymbol{n}}=(n^{X}_{0},n^{X}_{1},\dots,n^{X}_{d_{\rm
max}},n^{Y}_{E},n^{Z}_{E},n^{V})$. The transition intensities of
$\\{\hat{\boldsymbol{W}}^{(n)}(t)\\}$ are essentially the same of those of
$\\{\boldsymbol{W}^{(n)}(t)\\}$, given at the end of Section 5.1, except now
transitions of type (i) are partitioned according to whether or not the
infected susceptible has $I=\infty$. For $i=1,2,\dots,d_{\rm max}$, we have

$\hat{q}^{(n)}(\hat{\boldsymbol{n}},\hat{\boldsymbol{n}}-\boldsymbol{e}^{\rm
S}_{i}-\boldsymbol{e}^{\rm I}+(i-1)\boldsymbol{e}^{\rm
R})=\frac{n^{Y}_{E}in^{X}_{i}p_{i,k}(1-\pi)}{n^{X}_{E}+n^{Y}_{E}+n^{Z}_{E}-1},$

corresponding to the degree-$i$ infected susceptible having $I=0$, and

$\hat{q}^{(n)}(\hat{\boldsymbol{n}},\hat{\boldsymbol{n}}-\boldsymbol{e}^{\rm
S}_{i}+(i-1)\boldsymbol{e}^{\rm I}+\boldsymbol{e}^{\rm
V})=\frac{n^{Y}_{E}in^{X}_{i}p_{i,k}\pi}{n^{X}_{E}+n^{Y}_{E}+n^{Z}_{E}-1},$

corresponding to the degree-$i$ infected susceptible having $I=\infty$.

We apply the same random time-scale transformation to
$\\{\hat{\boldsymbol{W}}^{(n)}(t)\\}$ as done to
$\\{\boldsymbol{W}^{(n)}(t)\\}$ in Section 5.2. Denote the time-transformed
process by $\\{\bar{\boldsymbol{W}}^{(n)}(t)\\}$. Let
$\hat{\boldsymbol{l}}^{(1)}_{i,0}=\boldsymbol{l}_{i,0}^{(1)}$ and
$\hat{\boldsymbol{l}}^{(1)}_{i,i-1}=\boldsymbol{l}_{i,i-1}^{(1)}+\boldsymbol{e}^{\rm
V}$ $(i=1,2,\dots,d_{\rm max})$, and
$\hat{\boldsymbol{l}}^{(j)}_{+}=\boldsymbol{l}^{(j)}_{+}$ and
$\hat{\boldsymbol{l}}^{(j)}_{-}=\boldsymbol{l}^{(j)}_{-}$ $(j=2,3)$. The set
of possible jumps of $\\{\bar{\boldsymbol{W}}^{(n)}(t)\\}$ from a typical
state $\hat{\boldsymbol{n}}$ is now
$\hat{\Delta}=\hat{\Delta}_{1}\cup\hat{\Delta}_{2}\cup\hat{\Delta}_{3}$, where
$\hat{\Delta}_{1}=\bigcup_{i=1}^{d_{\rm
max}}\\{\hat{\boldsymbol{l}}^{(1)}_{i,0}\\}\cup\\{\hat{\boldsymbol{l}}^{(1)}_{i,i-1}\\}$,
$\hat{\Delta}_{2}=\\{\hat{\boldsymbol{l}}^{(2)}_{+},\hat{\boldsymbol{l}}^{(2)}_{-}\\}$
and
$\hat{\Delta}_{3}=\\{\hat{\boldsymbol{l}}^{(3)}_{+},\hat{\boldsymbol{l}}^{(3)}_{-}\\}$.

Let $\hat{\boldsymbol{w}}=(\boldsymbol{x},y_{E},z_{E},v)$. The family of
processes $\\{\bar{\boldsymbol{W}}^{(n)}(t)\\}$ is again asymptotically
density dependent with corresponding functions
$\bar{\beta}_{\boldsymbol{l}}(\hat{\boldsymbol{w}})$
$(\boldsymbol{l}\in\hat{\Delta})$ given by (cf. (5.3))

$\bar{\beta}_{\boldsymbol{l}}(\hat{\boldsymbol{w}})=\begin{cases}\bar{\beta}_{i,0}^{(1)}(\boldsymbol{x},y_{E},z_{E},v)=ix_{i}(1-\pi)&\text{
for }\boldsymbol{l}=\hat{\boldsymbol{l}}_{i,0}^{(1)}\in\hat{\Delta}_{1},\\\
\bar{\beta}_{i,i-1}^{(1)}(\boldsymbol{x},y_{E},z_{E},v)=ix_{i}\pi&\text{ for
}\boldsymbol{l}=\hat{\boldsymbol{l}}_{i,i-1}^{(1)}\in\hat{\Delta}_{1},\\\
\bar{\beta}^{(2)}_{+}(\boldsymbol{x},y_{E},z_{E},v)=y_{E}1_{\\{y_{E}>0\\}}&\text{
for }\boldsymbol{l}=\hat{\boldsymbol{l}}^{(2)}_{+},\\\
\bar{\beta}^{(2)}_{-}(\boldsymbol{x},y_{E},z_{E},v)=-y_{E}1_{\\{y_{E}<0\\}}&\text{
for }\boldsymbol{l}=\hat{\boldsymbol{l}}^{(2)}_{-},\\\
\bar{\beta}^{(3)}_{+}(\boldsymbol{x},y_{E},z_{E},v)=z_{E}1_{\\{z_{E}>0\\}}&\text{
for }\boldsymbol{l}=\hat{\boldsymbol{l}}^{(3)}_{+},\\\
\bar{\beta}^{(3)}_{-}(\boldsymbol{x},y_{E},z_{E},v)=-z_{E}1_{\\{z_{E}<0\\}}&\text{
for }\boldsymbol{l}=\hat{\boldsymbol{l}}^{(3)}_{-}.\end{cases}$ (5.54)

The associated drift function is (cf. (5.6))

$\displaystyle\bar{F}(\hat{\boldsymbol{w}})$
$\displaystyle=-\sum_{i=1}^{d_{\rm max}}ix_{i}\boldsymbol{e}^{\rm
S}_{i}+\left\\{\sum_{i=1}^{d_{\rm
max}}ix_{i}[(i-1)\pi-1]-2y_{E}-z_{E}\right\\}\boldsymbol{e}^{\rm I}$
$\displaystyle\qquad\qquad+\left[\sum_{i=1}^{d_{\rm
max}}i(i-1)(1-\pi)x_{i}-z_{E}\right]\boldsymbol{e}^{\rm
R}+\left[\sum_{i=1}^{d_{\rm max}}ix_{i}\pi\right]\boldsymbol{e}^{\rm V}.$

For $t\geq 0$, let
$\bar{\boldsymbol{w}}(t)=(\bar{x}_{0}(t),\bar{x}_{1}(t),\dots,\bar{x}_{d_{\rm
max}}(t),\bar{y}_{E}(t),\bar{z}_{E}(t),\bar{v}(t))$ be defined analogously to
$\tilde{\boldsymbol{w}}(t)$ at (5.6). Noting that $p_{I}=\pi$ and
$q_{I}=1-\pi$, the corresponding deterministic model for
$\bar{\boldsymbol{w}}(t)$ is given by (5.7)-(5.9), with $\tilde{x}_{i}$
replaced by $\bar{x}_{i}$ etc., augmented with

$\dfrac{d\bar{v}}{dt}=\pi\sum_{i=1}^{d_{\rm max}}i\bar{x}_{i}.$ (5.55)

Thus, (5.10)-(5.12) still hold, with the above change of notation, and, using
(5.10),

$\bar{v}(t)=\bar{v}(0)+\pi\sum_{i=1}^{d_{\rm max}}\bar{x}_{i}(0)(1-{\rm
e}^{-it}).$ (5.56)

The stopping time $\tilde{\tau}^{(n)}$ is unchanged, except now
$\varphi(\hat{\boldsymbol{w}})=\varphi(\boldsymbol{x},y_{E},z_{E},v)=y_{E}$.
Recall that site percolation corresponds to one initial infective, so under
both the MR and NSW models, $\epsilon_{i}=0$ and $\tilde{x}_{i}(0)=p_{i}$
$(i=0,1,\dots,d_{\rm max})$. Thus $\tilde{\tau}=\inf\\{t\geq
0:\bar{y}_{E}(t)=0\\}$ ($=\inf\\{t\geq 0:\tilde{y}_{E}(t)=0\\}$) is given by
the unique solution of (5.49) in $(0,\infty)$. Hence $z={\rm
e}^{-\tilde{\tau}}$ is given by the unique solution in $(0,1)$ of (2.10). Now
$\tilde{v}(0)=0$, since $V^{(n)}(0)=0$ for all $n=1,2,\dots$, so using (5.56)
and recalling that $\rho=1-f_{D}(z)$,

$\bar{v}(\tilde{\tau})=\pi\left[1-f_{D}({\rm
e}^{-\tilde{\tau}})\right]=\pi\rho.$

Let $\bar{V}^{(n)}=\bar{V}^{(n)}(\tilde{\tau}^{(n)})$, so in the site
percolation setting, $C^{(n)}\stackrel{{\scriptstyle
D}}{{=}}1+\bar{V}^{(n)}|G^{(n)}$. Now application of Corollary 4.4, as at
(5.27), yields

$\sqrt{n}\left(n^{-1}\bar{\boldsymbol{W}}^{(n)}(\tilde{\tau}^{(n)})-\bar{\boldsymbol{w}}(\tilde{\tau})\right)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}{\rm
N}(\boldsymbol{0},\bar{B}\bar{\Sigma}(\tilde{\tau})\bar{B}^{\top})\qquad\mbox{as
}n\to\infty,$

where

$\bar{B}=I-\frac{\bar{F}(\bar{\boldsymbol{w}}(\tilde{\tau}))\bigotimes\nabla\varphi(\bar{\boldsymbol{w}}(\tilde{\tau}))}{\nabla\varphi(\bar{\boldsymbol{w}}(\tilde{\tau}))\cdot\bar{F}(\bar{\boldsymbol{w}}(\tilde{\tau}))}$

and $\bar{\Sigma}(\tilde{\tau})$ is obtained by making obvious changes to
(5.24)-(5.26) to account for the extra dimension. Further,
$\bar{V}^{(n)}=\bar{\boldsymbol{W}}^{(n)}(\tilde{\tau}^{(n)})(\boldsymbol{0},0,0,1)^{\top}$,
so

$\sqrt{n}\left(n^{-1}\bar{V}^{(n)}-\pi\rho\right)\stackrel{{\scriptstyle
D}}{{\longrightarrow}}{\rm N}(0,\bar{\sigma}^{2})\qquad\mbox{as }n\to\infty,$

where

$\hat{\sigma}^{2}=(\boldsymbol{0},0,0,1)\bar{B}\bar{\Sigma}(\tilde{\tau})\bar{B}^{\top}(\boldsymbol{0},0,0,1)^{\top}.$

A simple calculation yields
$(\boldsymbol{0},0,0,1)\bar{B}=(\boldsymbol{0},-\pi b(\tilde{\tau}),0,1)$,
where $b(\tilde{\tau})$ is given by (5.29) with $a(\tilde{\tau})$ obtained by
replacing $f_{D_{\epsilon}}$ by $f_{D}$ in (5.23). Define $\bar{\Phi}(t,u)$
analogously to $\tilde{\Phi}(t,u)$ at (5.24). Then it follows using
(5.10)-(5.12) and (5.56) that

$\left[(\boldsymbol{0},0,0,1)\bar{\Phi}(\tilde{\tau},u)\right]_{i}=\begin{cases}\pi\left(1-re^{-i(\tilde{\tau}-u)}\right)&\text{
if }i=0,1,\dots,d_{\rm max},\\\ 0&\text{ if }i=d_{\rm max}+1,d_{\rm max}+2,\\\
1&\text{ if }i=d_{\rm max}+3,\end{cases}$

and $\left[(\boldsymbol{0},1,0,0)\bar{\Phi}(\tilde{\tau},u)\right]_{i}$ is
$0$, if $i=d_{\rm max}+3$, and given by the right-hand side of (5.30) (with
$p_{I}=\pi$ and $q_{I}=1-\pi$, if $i=0,1,\dots,d_{\rm max}+2$.

It follows that
$\bar{\boldsymbol{c}}(\tilde{\tau},u)=(\boldsymbol{0},0,0,1)\bar{B}\bar{\Phi}(\tilde{\tau},u)$
is given by

$\bar{\boldsymbol{c}}(\tilde{\tau},u)=-\pi(\boldsymbol{c}(\tilde{\tau},u),0)+\pi(\boldsymbol{1},0,0,0)+(\boldsymbol{0},0,0,1),$
(5.57)

where $\boldsymbol{c}(\tilde{\tau},u)$ is defined at (5.31), and (cf. (5.36))

$\bar{\sigma}^{2}=\bar{\boldsymbol{c}}(\tilde{\tau},0)\bar{\Sigma}_{0}\bar{\boldsymbol{c}}(\tilde{\tau},0)^{\top}+\sum_{\hat{\boldsymbol{l}}\in\hat{\Delta}}\int_{0}^{\tilde{\tau}}\left(\bar{\boldsymbol{c}}(\tilde{\tau},u)\hat{\boldsymbol{l}}^{\top}\right)^{2}\bar{\beta}_{\boldsymbol{l}}(\bar{\boldsymbol{w}}(u))\,{\rm
d}u,$

where $\bar{\Sigma}_{0}=\bar{\Sigma}^{\rm MR}_{0}$ or $\bar{\Sigma}^{\rm
NSW}_{0}$, depending on whether the random graph is MR or NSW. Here,
$\bar{\Sigma}^{\rm MR}_{0}$ and $\bar{\Sigma}^{\rm NSW}_{0}$ are the
asymptotic variance matrices of
$n^{-\frac{1}{2}}\bar{\boldsymbol{W}}^{(n)}(0)$ for the MR and NSW random
graphs, respectively; cf. (5.15) and (5.19).

We show in Appendix A.4 that

$\displaystyle\sum_{\hat{\boldsymbol{l}}\in\hat{\Delta}}\int_{0}^{\tilde{\tau}}\left(\bar{\boldsymbol{c}}(\tilde{\tau},u)\hat{\boldsymbol{l}}^{\top}\right)^{2}\bar{\beta}_{\boldsymbol{l}}(\bar{\boldsymbol{w}}(u))\,{\rm
d}u=\pi(1-\pi)[\rho-2h(z)(1-z)\mu_{D}]$ (5.58)
$\displaystyle\qquad\qquad\qquad+\pi^{2}\sum_{\boldsymbol{l}\in\Delta}\int_{0}^{\tilde{\tau}}\left(\boldsymbol{c}(\tilde{\tau},u)\boldsymbol{l}^{\top}\right)^{2}\tilde{\beta}_{\boldsymbol{l}}(\tilde{\boldsymbol{w}}(u))\,{\rm
d}u$

and, for both the MR and NSW random graphs, that

$\bar{\boldsymbol{c}}(\tilde{\tau},0)\bar{\Sigma}_{0}\bar{\boldsymbol{c}}(\tilde{\tau},0)^{\top}=\boldsymbol{c}(\tilde{\tau},0)\Sigma_{0}\boldsymbol{c}(\tilde{\tau},0)^{\top}.$
(5.59)

Equations (2.13) and (2.14) then follow using Theorem 2.2, noting that
$q_{I}^{(2)}=\pi(1-\pi)$.

##  6 Concluding comments

A shortcoming of our results, from a mathematical though not a practical
viewpoint, is the requirement of a maximal degree $d_{\rm max}$. It seems
likely that Theorems 2.1-2.3 continue to hold when that requirement is
relaxed, subject to appropriate conditions on the degree sequence (MR model)
or degree distribution $D$ (NSW model). This conjecture is supported by the
recent work of [Barbour and Röllin(2017)], who use Stein’s method to obtain a
central limit theorem for local graph statistics in the MR model, allowing for
unbounded degree, and as an application obtain a central limit theorem for the
size of the giant component with asymptotic variance given by
$\tilde{\sigma}^{2}_{\rm MR}$ with $p_{I}=1$. To extend the present proof to
models with unbounded degree would require a functional central limit theorem
for density dependent population processes with countable state spaces.
[Barbour and Luczak(2012)] give such a theorem but it is not applicable in our
setting as it would require a finite upper bound on the number of neighbours
an individual can infect.

The central limit theorems can be extended, at least in principle, to allow
for the infection rate $\lambda$ to depend on the degree of an infective, and
also to more general infection processes in which the set of its neighbours
that are contacted by a given infective is a symmetric sampling procedure
([Martin-Löf(1986)]). In both cases it is straightforward to determine the
limiting deterministic model in Section 5.3.1 and the equation corresponding
to (5.22), which governs $\tilde{\tau}$, but calculation of the asymptotic
variances in the central limit theorems is likely to be prohibitive.

The configuration model does not display clustering in the limit as
$n\to\infty$ and several authors have considered modifications of the
configuration model that introduce clustering. In [Trapman(2007)] and
[Coupechoux and Lelarge(2014)], in the configuration model construction, for
$d=1,2,\dots$, some individuals having $d$ half-edges are replaced by fully
connected cliques, each of size $d$, with each member of a clique having
exactly one half-edge. The half-edges are then paired up in the usual fashion.
In [Gleeson(2009)] and [Ball et al.(2010)], the network is formed as in the
configuration model and the population is also partitioned into fully
connected cliques. In both models, the set of edges in the network is the
union of those in cliques and the paired half-edges. The methodology in this
paper can be extended to this general class of models as follows.

As in Section 5.1, the network and epidemic are constructed simultaneously.
The objects counted are now fully susceptible cliques, typed by their size and
degree composition, and infective and recovered half-edges. When infection is
transmitted down a half-edge that half-edge is paired with a uniformly chosen
half-edge as before. If it is paired with a susceptible half-edge, then an
epidemic is triggered within the corresponding clique and associated half-
edges, leading to the creation of further infective and recovered half-edges
and, unless the clique epidemic infects the entire clique, a new susceptible
clique having reduced size. Central limit theorems for the final size of
epidemics on MR and NSW versions of such random graphs should follow using
similar arguments to before but again calculation of the asymptotic variances
may be difficult. If the infectious period is constant, so the epidemic model
is equivalent to bond percolation on the network, the analysis may perhaps be
simplified by first splitting the cliques into components determined by bond
percolation and then using similar methods to the present paper treating the
components as super-individuals.

##  Appendix A Detailed derivations

###  A.1 Deterministic solution (5.10)-(5.12)

First note that (5.10) follows immediately from (5.7). Substituting (5.10)
into (5.9) yields

$\dfrac{d\tilde{z}_{E}}{dt}=q_{I}\left[\sum_{i=2}^{d_{\rm
max}}i(i-1)\tilde{x}_{i}(0){\rm e}^{-it}\right]-\tilde{z}_{E}$

and (5.12) follows using the integrating factor ${\rm e}^{t}$.

Multiplying (5.7) by $i$ and adding over $i=1,2,\dots,d_{\rm max}$ yields

$\dfrac{d\tilde{x}_{E}}{dt}=-\sum_{i=2}^{d_{\rm
max}}i(i-1)\tilde{x}_{i}-\tilde{x}_{E}.$ (A.1)

Recall that
$\tilde{\eta}_{E}(t)=\tilde{x}_{E}(t)+\tilde{y}_{E}(t)+\tilde{z}_{E}(t)$.
Summing (A.1), (5.8) and (5.9) gives

$\dfrac{d\tilde{\eta}_{E}}{dt}=-2\tilde{\eta}_{E},$

and (5.13) follows. Equation (5.11) then follows since
$\tilde{y}_{E}(t)=\tilde{\eta}_{E}(t)-\tilde{x}_{E}(t)-\tilde{z}_{E}(t)$ and,
using (5.10), $\tilde{x}_{E}(t)=\sum_{i=1}^{d_{\rm max}}i\tilde{x}_{i}(0){\rm
e}^{-it}$.

###  A.2 Properties of $\tilde{\tau}$ and $a(\tilde{\tau})$

In this appendix we show that $\tilde{\tau}\in(0,\infty)$ and
$a(\tilde{\tau})<0$, first in the setting of Theorem 2.1 and then in the
setting of Theorems 2.2 and 2.3. Let
$G(s)=p_{I}f_{D_{\epsilon}}^{(1)}(s)-\mu_{D}(s-q_{I})$ $(0\leq s\leq 1)$.
Then, from (5.22), $z={\rm e}^{-\tilde{\tau}}$ satisfies $G(z)=0$. Now
$G(0)=p_{I}f_{D_{\epsilon}}^{(1)}(0)+\mu_{D}q_{I}>0$, unless $q_{I}=0$ and
$p_{1}-\epsilon_{1}=0$. Also $f_{D_{\epsilon}}^{(1)}(1)<\mu_{D}$, since
$p_{i}\epsilon_{i}>0$ for at least one $i>0$, so $G(1)<0$. Thus, under the
conditions of Theorem 2.1, $G(s)$ has at least one zero in $(0,1)$. Moreover
it has precisely one zero, $z$ say, as $G(s)$ is convex on $[0,1]$, since
$G^{(2)}(s)=p_{I}f_{D_{\epsilon}}^{(3)}(s)\geq 0$ for all $s\in[0,1]$. Hence
$\tilde{\tau}=-\log z\in(0,\infty)$, as required, and it follows from (5.23)
that $a(\tilde{\tau})<0$ if and only if $G^{(1)}(z)<0$. Suppose, for
contradiction, that $G^{(1)}(z)\geq 0$. Then, since $G^{(2)}(s)\geq 0$ for all
$s\in[z,1]$,

$G(1)\geq G(z)+\int_{z}^{1}G^{(1)}(s)\,{\rm d}s\geq 0.$

However, $G(1)<0$, so $G^{(1)}(z)<0$, as required.

Turning to the setting of Theorems 2.2 and 2.3, now let
$G(s)=p_{I}f_{D}^{(1)}(s)-\mu_{D}(s-q_{I})$ $(0\leq s\leq 1)$. Then, from
(5.49), $z={\rm e}^{-\tilde{\tau}}$ satisfies $G(z)=0$. Now $G(1)=0$ and,
under the conditions of Theorems 2.2 and 2.3, $G(0)>0$. Also, using (2.8),
$G^{(1)}(1)=p_{I}f_{D}^{(2)}(1)-\mu_{D}>0$, since $R_{0}>1$. Thus, since
$G(s)$ is convex on $[0,1]$, there exists a unique $z\in(0,1)$ such that
$G(z)=0$, so $\tilde{\tau}\in(0,\infty)$. Moreover, $G^{(2)}(s)>0$ for
$s\in[z,1]$ and a similar contradiction argument to before shows that
$G^{(1)}(z)<0$, whence $a(\tilde{\tau})<0$.

###  A.3 Asymptotic variances $\sigma^{2}_{\rm MR}$ and $\sigma^{2}_{\rm NSW}$

In this appendix we derive the expressions for $\sigma^{2}_{\rm MR}$ and
$\sigma^{2}_{\rm NSW}$ given in Theorem 2.1. Recalling the partition
$\Delta=\Delta_{1}\cup\Delta_{2}\cup\Delta_{3}$ defined just before (5.1), it
follows from (5.36) that

$\sigma^{2}_{\rm
MR}=\sigma_{0,MR}^{2}+\sum_{i=1}^{3}\sigma_{i}^{2}\qquad\mbox{and}\qquad\sigma^{2}_{\rm
NSW}=\sigma_{0,NSW}^{2}+\sum_{i=1}^{3}\sigma_{i}^{2},$ (A.2)

where

$\sigma_{0,MR}^{2}=\boldsymbol{c}(\tilde{\tau},0)\Sigma^{\rm
MR}_{0}\boldsymbol{c}(\tilde{\tau},0)^{\top},\qquad\sigma_{0,NSW}^{2}=\boldsymbol{c}(\tilde{\tau},0)\Sigma^{\rm
NSW}_{0}\boldsymbol{c}(\tilde{\tau},0)^{\top}$ (A.3)

and, for $i=1,2,3$,

$\sigma_{i}^{2}=\int_{0}^{\tilde{\tau}}\sum_{\boldsymbol{l}\in\Delta_{i}}\left(\boldsymbol{c}(\tilde{\tau},u)\boldsymbol{l}^{\top}\right)^{2}\tilde{\beta}_{\boldsymbol{l}}(\tilde{\boldsymbol{w}}(u))\,{\rm
d}u.$

We calculate $\sum_{i=1}^{3}\sigma_{i}^{2}$ in Appendix A.3.1,
$\sigma_{0,MR}^{2}$ and then $\sigma^{2}_{\rm MR}$ in Appendix A.3.2, and
$\sigma_{0,NSW}^{2}$ and then $\sigma^{2}_{\rm NSW}$ in Appendix A.3.3.

####  A.3.1 Calculation of $\sum_{i=1}^{3}\sigma_{i}^{2}$

Noting that $\tilde{y}_{E}(t)\geq 0$ and $\tilde{z}_{E}(t)\geq 0$ for $0\leq
t\leq\tilde{\tau}$, it follows from (5.3), (5.33) and (5.34) that

$\sigma_{2}^{2}=\int_{0}^{\tilde{\tau}}4c_{I}(\tilde{\tau},u)^{2}\tilde{y}_{E}(u)\,{\rm
d}u$ (A.4)

and

$\displaystyle\sigma_{3}^{2}$ $\displaystyle=$
$\displaystyle\int_{0}^{\tilde{\tau}}(c_{I}(\tilde{\tau},u)+c_{R}(\tilde{\tau},u))^{2}\tilde{z}_{E}(u)\,{\rm
d}u$ $\displaystyle=$
$\displaystyle\int_{0}^{\tilde{\tau}}(2c_{I}(\tilde{\tau},u)-c_{J}(\tilde{\tau},u))^{2}\tilde{z}_{E}(u)\,{\rm
d}u,$

where

$c_{J}(\tilde{\tau},u)=b(\tilde{\tau}){\rm e}^{-(\tilde{\tau}-u)}.$ (A.6)

Now

$\sigma_{1}^{2}=\int_{0}^{\tilde{\tau}}\tilde{g}(u)\,{\rm d}u,$ (A.7)

where, using (5.3),

$\displaystyle\tilde{g}(u)$ $\displaystyle=$
$\displaystyle\sum_{\boldsymbol{l}\in\Delta_{1}}\left(\boldsymbol{c}(\tilde{\tau},u)\boldsymbol{l}^{\top}\right)^{2}\tilde{\beta}_{\boldsymbol{l}}(\tilde{\boldsymbol{w}}(u))$
(A.8) $\displaystyle=$ $\displaystyle\sum_{i=1}^{d_{\rm
max}}\sum_{k=0}^{i-1}\left[\boldsymbol{c}(\tilde{\tau},u)(\boldsymbol{l}_{ik}^{(1)})^{\top}\right]^{2}p_{i,k}i\tilde{x}_{i}(u).$

Recalling (5.35), it follows using (5.3) that, for $i=1,2,\dots,d_{\rm max}$
and $k=0,1,\dots,i-1$,

$\boldsymbol{c}(\tilde{\tau},u)(\boldsymbol{l}_{ik}^{(1)})^{\top}=(p_{I}b(\tilde{\tau})-1){\rm
e}^{-i(\tilde{\tau}-u)}-2c_{I}(\tilde{\tau},u)+[iq_{I}-(i-k-1)]c_{I}(\tilde{\tau},u).$
(A.9)

Recall that $p_{i,k}={\rm P}(X=k)$, where $X\sim{\rm Bin}(i-1,1-\exp(-\lambda
I))$; see Section 5.1. Elementary calculation yields, for $i=1,2,\dots,d_{\rm
max}$, that

$\sum_{k=0}^{i-1}(i-k-1)p_{i,k}=(i-1)q_{I}$ (A.10)

and

$\sum_{k=0}^{i-1}(i-k-1)^{2}p_{i,k}=(i-1)(i-2)q_{I}^{(2)}+(i-1)q_{I}.$ (A.11)

Using (A.10) and (A.11), it follows from (A.9) and some algebra that, for
$i=1,2,\dots,d_{\rm max}$

$\displaystyle\sum_{k=0}^{i-1}$
$\displaystyle\left[\boldsymbol{c}(\tilde{\tau},u)(\boldsymbol{l}_{ik}^{(1)})^{\top}\right]^{2}p_{i,k}=[2c_{I}(\tilde{\tau},u)-q_{I}c_{J}(\tilde{\tau},u)]^{2}$
(A.12)
$\displaystyle\qquad+q_{I}p_{I}c_{J}(\tilde{\tau},u)^{2}(i-1)+(q_{I}^{(2)}-q_{I}^{2})c_{J}(\tilde{\tau},u)^{2}(i-1)(i-2)$
$\displaystyle\qquad+2[2c_{I}(\tilde{\tau},u)-q_{I}c_{J}(\tilde{\tau},u)]{\rm
e}^{-i(\tilde{\tau}-u)}[1-ib(\tilde{\tau})p_{I}]$ $\displaystyle\qquad+{\rm