SIR epidemics and vaccination on random graphs with clustering

In this paper we consider Susceptible \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rightarrow $$\end{document}→ Infectious \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rightarrow $$\end{document}→ Recovered (SIR) epidemics on random graphs with clustering. To incorporate group structure of the underlying social network, we use a generalized version of the configuration model in which each node is a member of a specified number of triangles. SIR epidemics on this type of graph have earlier been investigated under the assumption of homogeneous infectivity and also under the assumption of Poisson transmission and recovery rates. We extend known results from literature by relaxing the assumption of homogeneous infectivity both in individual infectivity and between different kinds of neighbours. An important special case of the epidemic model analysed in this paper is epidemics in continuous time with arbitrary infectious period distribution. We use branching process approximations of the spread of the disease to provide expressions for the basic reproduction number \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0$$\end{document}R0, the probability of a major outbreak and the expected final size. In addition, the impact of random vaccination with a perfect vaccine on the final outcome of the epidemic is investigated. We find that, for this particular model, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0$$\end{document}R0 equals the perfect vaccine-associated reproduction number. Generalizations to groups larger than three are discussed briefly.

duration of the contacts between individuals typically depend on the nature of their relationship. For this reason, recent interest has focused on the impact of the underlying social network on the spread of the disease. The social network is typically represented by a random graph , in which the nodes or vertices represent individuals and the edges represent social contacts between the individuals. Two nodes that share an edge are called "neighbours".
A popular choice when generating random graphs with a specified degree distribution is the configuration model (CM). It was introduced by Bollobás (1980) for the special case where the degree distribution is degenerate (i.e. every node of the graph has the same degree) and extended to more general degree distributions by Reed (1995, 1998). There is a vast literature on epidemics on configuration model graphs [see e.g. Andersson (1999), Britton et al. (2007), Janson et al. (2014), Barbour and Reinert (2013), Bhamidi et al. (2014)].
An important feature of the configuration model is that, under mild regularity conditions on the degrees, this type of graph is asymptotically unclustered. That is to say, it contains virtually no groups and short circuits. Real world networks do, however, typically exhibit clustering (Newman 2003), and there are a number of graph models that do allow for group structure (Bollobás et al. 2011;Karoński et al. 1999;Newman 2002). Epidemics on graphs with group structure were studied by Trapman (2007); Ball et al. (2009Ball et al. ( , 2010Ball et al. ( , 2014; Coupechoux and Lelarge (2015); Britton et al. (2008).
In this paper, we use a generalized version of the configuration model to incorporate clustering of the social network in the analysis of the spread of an infectious disease. The configuration model with clustering (CMC) was independently introduced by Miller (2009) and Newman (2009). It is an extension of the CM in the sense that, for each node u, in addition to the degree of u one also specifies the number of pairs of neighbours of u that are in turn neighbours of each others. In other words, one specifies the number of triangles (with non-overlapping edges) of which u is a member [see Sect. 2.1 for a precise definition of the graph model]. This allows for graphs with non-negligible clustering and a specified degree distribution. That is to say, the CMC deviates from the classical Erdős-Rényi graph model (Erdős and Rényi 1959) in two fundamental ways: it allows for for a non-Poissonian degree distributions and is asymptotically clustered. Epidemics on this type of graph have previously been studied by Miller (2009) and Volz (2011). Miller (2009) investigated the impact of clustering on the epidemic threshold, formulated as a bond percolation problem. This means that the infectivity of infected individuals is assumed to be homogeneous; an infected individual transmits the disease to each of its neighbours independently with some fixed probability T . Volz (2011) investigated the time evolution and final size of epidemics on CMC graphs under the assumption of exponentially distributed infectious periods during which individuals contact neighbours at a constant rate.
The main contribution of our research is that we extend the results of Miller (2009) andVolz (2011) by allowing for heterogeneous infectivity, i.e. by allowing for some infected individuals to be more contagious than others or for individuals to exhibit different contact behaviors for different types of neighbours. Such heterogeneity may, for instance, reflect variability in the infectious period or contact preferences on the part of individuals. We provide expressions for the probability of a major outbreak and the final size of a major outbreak. A key tool in our analysis is the approximation of the epidemic seen from a "generation of infection" or "rank" perspective by a multitype Galton Watson branching process. This approximation, which is interesting in its own right, gives rise to the rank-based reproduction number R 0 [see e.g. Pellis et al. (2008Pellis et al. ( , 2012]. We note that especially allowing for heterogeneity in the infectivity of individuals requires a more intricate branching process approximation than a model with homogeneous infectivity [as analysed by Miller (2009)]. To see this, consider an individual v which is infected through a triangle Δ. The "local epidemics" in Δ and in other triangles v is part of all depend on the infectivity of v and are therefore in general not independent.
The second contribution of this paper concerns vaccination. We investigate the impact of uniform vaccination (i.e. vaccinated individuals are selected uniformly at random) with a perfect vaccine (i.e. a vaccine that provides full and permanent immunity to the disease). We find that it is necessary to vaccinate a fraction 1 − 1/R 0 of the population in order to prevent a major outbreak of the disease, as in the case of homogeneous mixing. We illustrate our findings with numerical examples. This paper is structured as follows. In Sect. 2 we provide the preliminaries for the model. In Sect. 2.1 we give a more detailed description of how graphs are generated in the CMC and investigate the asymptotic clustering of such graphs and in Sect. 2.2 the epidemic model is specified. Sections 2.3 and 2.4 contains an overview of the concept of reproduction numbers and the necessary branching process background. In Sect. 3, we derive expressions for the probability of a major outbreak and the expected final size under the assumption of an unvaccinated and fully susceptible population, and in Sect. 4 the analysis is repeated under the assumption of uniform vaccination with a perfect vaccine. We illustrate our findings with numerical examples presented in Sect. 5 and discuss possible extensions in Sect. 6.

The configuration model with clustering
A CMC graph is constructed as follows. Let { p(k s , k Δ )} k s ,k Δ ∈N 0 be a prescribed joint degree distribution, where k s denotes the number of single edges attached to a node, and k Δ denotes the number of pairs of triangle edges. Throughout, (S, Δ) is assumed to be a generic random vector distributed according to p.
be a sequence of independent copies of (S, Δ). Analogously to the CM, a graph G N = G N ( p) of size N is constructed by first assigning the single degree S i and the triangle degree Δ i to the node v i , i = 1, 2, . . . , N . One may think of this step in terms of half-edges; to each node v i , we attach S i single half-edges and Δ i pairs of triangle half-edges. The single half-edges are then matched in pairs and the triangle half-edge pairs in threes by choosing a matching uniformly at random among all possible such matchings. The process of joining half-edges is illustrated in Fig. 1. As described in Miller (2009), the matching may be carried out as follows. Two lists of nodes, one single degree list and one triangle degree list are created. A node with joint degree (k s , k Δ ) appears k s times in the single list and k Δ times in the triangle list. The lists are then shuffled uniformly, Fig. 1 Schematic illustration of the construction of a CMC graph. Triangle half-edges (marked with a triangle) and single half-edges (marked with a perpendicular line) are assigned to the nodes of the graph (left). The half-edges are then matched uniformly at random (right). Note that two of the half-edges attached to v 3 are paired with each other and so form a self-loop and the nodes on positions 2m + 1 and 2m + 2 in the single degree list and positions 3m + 1, 3m + 2 and 3m + 3 in the triangle degree list are matched, m ∈ N 0 .
We define the total single degree as If the total single degree (that is, the length of the single degree list) is not even or if the total triangle edge degree (the length of the triangle degree list) is not a multiple of three we erase a single half-edge and/or one or two triangle half-edge pairs chosen uniformly at random. Similarly, we erase self-loops and merge multiple edges, so that the resulting graph is simple. Under assumption A1 (stated below) on p it holds that the number of single self-loops and single double edges converge in distribution to independent Poisson random variables with finite means [cf. Van der Hofstad (2016, prop. 7.13)].
For this reason, self-loops and multiple edges are negligible in the limit as N → ∞. In the remainder of this paper, we ignore the small differences in the topology of the graph that arise from erasing multiple edges or self-loops. In addition, we ignore the small differences in effective degree distribution that arise from erasing half-edges so that the number of single and triangle half-edges are multiples of two and three, respectively.
We make the following assumptions on p.
Note that the assumption A1 implies E(ΔS) < ∞. Assumption A2 ensures that the mean matrices of the approximating branching processes (presented below) are positively regular (we say that an r × r matrix M is positively regular if it has finite non-negative entries and for some n ∈ N all entries of M n are strictly positive).

Clustering coefficient of G N
For any undirected graph we can measure the amount of clustering in the network using the so-called clustering coefficient, which is defined as follows. Let G = (V , E) be an undirected graph with node set V and edge set E. Define the set of all ordered wedges (i.e. directed paths consisting of precisely two edges) of G and the set of all ordered triangles of G. The clustering coefficient C(G) of G is a measure of the degree of clustering of G and is defined as the fraction of the ordered wedges of G that are also triangles: Here | · | denotes the cardinality of a set. As stated in the following proposition, CMC graphs have asymptotically non-zero clustering as N → ∞. An analogous result for fixed degree sequences was presented in Newman (2009). Let P −→ denote convergence in probability.
Proposition 1 Let {G N } N be a sequence of CMC graphs with independent degrees drawn from p. If p satisfies assumption A1 then The proof is presented in the Appendix.

Downshifted size-biased degrees
The graph G N may be constructed by joining the half-edges in a random order. In particular, G N may be constructed as the epidemic progresses; starting with the initial infected case we sequentially match the half-edges along which the disease is transmitted. Since half-edges are chosen uniformly at random in the matching procedure, the probability to choose a specific node is proportional to the number of free half-edges attached to the node in question. That is, if we pair a single half-edge, the probability of choosing a specific node with k s unpaired single half-edges is proportional to k s . For this reason, the degree distribution of a node explored by joining a single half-edge in the early phase of the epidemic can be approximated by the single size biased [cf. the concept of excess degree in Meyers (2007)] degree distribution p (s) Similarly, the degree distribution of the nodes explored by joining three triangle half-edge pairs in the early phase of the epidemic can be approximated by the triangle size biased degree distribution p (Δ) In the epidemic process, we need to account for the fact that an infected individual has at least one non-susceptible neighbour (namely the direct source of its infection). For this reason, we introduce the downshifted size biased degree distributions p (s) • and p (Δ) • , given by Throughout, we will make frequent reference to the following random vectors and the expected values

The epidemic model
We use an SIR model to investigate the dynamics of the spread of the disease. At any given time point, the population is divided into three groups, depending on health status. The groups are susceptible (S), infectious (I) and recovered (R) [see e.g. Britton (2010), Diekmann et al. (2013)]. Individuals of the population make contact with other individuals at (possibly random) points in time. If, at some time point, an infectious individual contacts a susceptible individual then the susceptible individual instantaneously becomes infectious. An infectious individual will cease to be contagious after a period of time, which we call the infectious period of the individual in question, and is then transferred to the recovered group. Recovered individuals are those that are immune to the disease. Individuals belonging to this group play no further role in the spread of the disease. Because of this last observation, we can treat individuals that die because of the disease as "recovered". In summary, we allow only the transitions S → I and I → R. Note that the population is assumed to be closed; we ignore births, deaths and migration. More specifically, we consider an SIR epidemic in a generation framework on the clustered static graph G N and assume possible heterogeneity in infectivity, both between different individuals (individual heterogeneity) and between different kinds of edges (edge heterogeneity). Individual heterogeneity means that some infected individuals are more contagious than others. Such heterogeneity may, for instance, arise from variability in the infectious period. Edge heterogeneity reflects that individuals may exhibit different contact behaviors for different types of neighbours ; an individual may for instance prefer to spend more time with its triangle neighbours at the expense of spending less time with its single neighbours.
To construct a model that captures such heterogeneities, let T = (T s , T Δ ) be a random variable with support in [0, 1] 2 , and let be a sequence of independent copies of T . We allow for any dependence structure between T s and T Δ . Each node v i of G N is equipped with a two-dimensional transmission weight T i . If v i gets infected, then each susceptible single neighbour (neighbour by virtue of a single edge) of v i gets infected by v i independently in the next generation with probability T (i) s , and each susceptible triangle neighbour (neighbour by virtue of a triangle edge) of v i gets infected by v i independently in the next generation with probability Node v i thereafter becomes recovered, playing no further role in the epidemic. An infected node transmits the disease independently of the transmissions from other infected nodes. An infected node does not, however, transmit the disease to its neighbours independently, unless the distribution of T is degenerate. Conditioned on the transmission weights {T i } i and the structure of G N , the number of single and triangle neighbours that an infected node v i makes (infectious) contact with while infectious has a binomial distribution with parameters , respectively. The spread of this epidemic can be fully captured by a directed graph [see e.g. Pellis et al. (2012), Kenah and Miller (2011)]. To construct such directed graph from an undirected CMC graph G N , we replace each undirected edge of G N by two parallel directed edges, pointing in the opposite direction. The weight of an edge (v i , v j ), which represents the (potential) transmission time from v i to v j , is taken to be 1 if v i would make infectious contact with v j if infected, and ∞ otherwise. The individuals ultimately infected are then the individuals that can be reached from an initial case by following a path consisting of directed edges with finite edge weights.

Reproduction numbers
A key quantity in the study of epidemics is the basic reproduction number, often denoted by R 0 . It is usually defined as the expected number of infected cases caused by a "typical" infected individual in an otherwise susceptible population (Diekmann et al. 1990). For most stochastic epidemic models [including SIR epidemics in homogeneous mixing populations (Britton 2010), populations with households (Ball et al. 2016) and epidemics on networks (Britton et al. 2007)] it has the threshold property that a major outbreak is possible if and only if R 0 > 1.
For models where a suitable generation based branching process approximation is available, R 0 is usually defined as the Perron root (the dominant eigenvalue, which exists and is real-valued by assumptions A1 and A2, see for instance Varga (2009, Chapter 2) of the mean matrix of the approximating Galton Watson branching process. This is the definition used in this article. By standard branching process theory, the interpretation of R 0 as the expected number of cases caused by the typical individual in the early phase of the epidemic and its threshold properties are retained by this definition. The threshold property of R 0 is made precise in Theorem 1 below.
In Sect. 4, we investigate the spread of an epidemic in a population with vaccination. To this end, in addition to the basic reproduction number R 0 , we consider the perfect vaccine-associated reproduction number R V (Goldstein et al. 2009). A vaccine is perfect if it provides full and permanent immunity. That is, an individual vaccinated with a perfect vaccine cannot contract the disease. The perfect vaccine-associated reproduction number R V is defined as (Ball et al. 2016) where the critical vaccination coverage f (c) v is the fraction of the population that has to be vaccinated with a perfect vaccine in order to reduce R 0 to unity, if the vaccinated individuals are chosen uniformly at random. That is to say, f the fraction necessary to vaccinate in order to be guaranteed to prevent a major outbreak (Britton 2010) For many models, including epidemics on graphs generated by the CM (Britton et al. 2007) and the standard stochastic SIR epidemic model [i.e. individuals mix homogeneously, see for instance Britton (2010)], R V = R 0 . That is, vaccinating a fraction 1−1/R 0 of the population with a perfect vaccine is sufficient to surely prevent a major outbreak. On the other hand, for the households and households-workplaces model with uniform vaccination, R V ≥ R 0 (Ball et al. 2016) with strict inequality possible. In Sect. 4.1 we show that for the model analysed in this report, R V = R 0 .

Epidemics in continuous time: the rank-based approach
As mentioned above, heterogeneity in infectivity might arise from heterogeneity in the infectious period; an important special case of the above described model is epidemics in continuous time with random infectious periods where contacts between individuals take place according to point processes on R ≥0 . Ignoring the real time-dynamics of an epidemic does not impact results that concern the final outcome of the epidemic. This result was first presented by Ludwig (1975), see also Pellis et al. (2008) or Kiss et al. (2017, section 6.2.3) for a more recent discussion. This leads us to the often more tractable rank-based approach.
In the rank-based approach, however, v 3 is attributed to v 1 . Right: The resulting rank generation tree In order to define the rank of a vertex, denote the initial case by v * . The rank of a node v in G N is the distance from v * to v, if every edge along which the disease would be transmitted is assigned the edge weight 1, and every other edge is assigned the edge weight ∞. That is, the rank of v is the smallest number of directed edges that have to be traversed in order to follow a path of (potential) transmission from v * to v. We may then analyse the spread of the disease by letting generation n of the epidemic process consist of the individuals of rank n. If, for instance, v 1 is the first node in a triangle consisting of the nodes v 1 , v 2 , v 3 to be infected, and v 1 infects v 2 and thereafter attempts to infect v 3 , then v 3 is attributed to v 1 regardless of whether v 1 or v 2 infected v 3 . This is illustrated in Fig. 2.
Consider a continuous time epidemic formulated as follows. Suppose that each infected individual remains infectious for a (random) period of time. The infectious periods are distributed as the random variable τ , τ ∼ F, and independent (but identically distributed) for different nodes. Suppose further that a node makes contact with each neighbour independently at a Poisson rate while infected, and that susceptible individuals are fully susceptible, so that each infectious-susceptible contact results in transmission. If the contact rate is given by β s for single edge neighbours and β Δ for triangle edge neighbours then the transmission weights T s and T Δ are distributed as 1 − e −β s τ and 1 − e −β Δ τ , respectively. We then have is the Laplace transform of the infectious period.

Branching process approximations
To analyse the spread of the disease in the early stages of the epidemic, we employ a multi-type branching process approximation. The graph G N may be constructed by joining the half-edges in any (possibly random) order, provided that the uniform matching is not violated. In particular, the graph G N may be constructed (or explored) as the epidemic progresses; starting with the initial infected case u * we sequentially match the half-edges along which the disease is transmitted. In the early phase of the epidemic, short cycles (except for the triangles formed by triangle edges) are unlikely to appear. For these reasons, the early spread of the disease is well approximated by a suitably chosen branching process.
Similarly, a branching process approximation can be used to approximate the expected final size of the epidemic (Ball et al. 2009(Ball et al. , 2010(Ball et al. , 2014. In the graph representation of an epidemic, an individual contracts the disease if and only if there is a path of directed edges with finite edge weights from the initial case to the node representing the individual in question. Define the susceptibility set S(v) = S N (v) of a node v as the collection of nodes of G N that can be reached from v by tracing a path of finite length backwards. That is, the individuals that contract the disease are precisely the individuals with susceptibility sets that contain an initial case. Hence, if the initial case is chosen uniformly at random then the probability that a node v contracts the disease is proportional to the size of its susceptibility set S(v) and this probability can be approximated by exploring S(v). Figure 3 shows a schematic illustration of a susceptibility set.
By reversing the direction of the edges of the graph representation of an epidemic, but keeping the weights, the expected final fraction of the population infected in a major outbreak and the probability of a major outbreak are interchanged (Miller 2008), provided that the initial case is chosen uniformly at random. The process so obtained is called the backward epidemic process of the node v. If the underlying epidemic model is such that the backward epidemic process can be well approximated by a branching process, then we can use this branching process to compute the asymptotic distribution of the proportion of the population that ultimately escapes infection. This is made precise in the following theorem, due to Ball et al. (2014, Theorem 3.5), who proved the theorem for the related model of random intersection graphs. The statement of Theorem 1 carries over to the forward and backward branching processes considered in this paper. We omit the proof, which is analogous to the proof presented in Ball et al. (2014), see also Ball et al. (2009).

Theorem 1 Let q and q b be the extinction probabilities of the forward and backward approximating branching processes respectively, and let S N be the proportion of the population that ultimately escapes an epidemic in a population of size N . Then
In other words, in the limit of large population sizes, the epidemic "takes off" with probability 1 − q, and if this happens a fraction 1 − q b of the population is ultimately infected (with probability converging to 1 as N → ∞). Note that since R 0 is defined as the Perron root of the mean matrix of the forward branching process, q < 1 if and only if R 0 > 1.

An epidemic in a fully susceptible population
We now have the tools to analyse the spread of an infectious disease on a graph generated by the CMC. In the present section, the population is assumed to be fully susceptible to the disease, apart from the initially infectious individual.

Forward process
Before analysing the forward process, we need to set some terminology. For a given triangle u, v, w, where u is the first individual to be infected in the triangle u, v, w, we refer to v and w as twins. We approximate the spread of the disease during the early phase by a multi-type branching process consisting of the following three types (except for the initial case): Type 1: A node infected along a triangle edge whose twin (in the same triangle) is infected at the same time step or earlier Type 2: A node infected along a triangle edge that is not of type 1 Type 3: A node infected along a single edge Figure 4 shows three examples of possible paths of transmission within a triangle giving rise to type 1 and 2 individuals in the approximating branching process.
Denote by m 1,1 m 1,2 m 1,3 m 2,1 m 2,2 m 2,3 m 3,1 m 3,2 m 3,3 ⎞ ⎠ the mean matrix of the above described branching process. Suppose that v 1 is the first individual to be infected in the triangle v 1 , v 2 , v 3 . The probability that v 1 transmits the disease both to v 2 and v 3 is E(T 2 ). Similarly, the probability that v 1 transmits the disease to either v 2 or v 3 , but not to both, is 2E (T (1 − T )). Thus, by linearity of expectation and because the distribution of the susceptible neighbours of infected nodes in the early phase of the epidemic is given by the downshifted degree distributions in (4), we obtain (5) have the downshifted size biased distributions). Note that all entries of M f are finite and that S and Δ both have finite second moments by assumption A1.
If M f is positively regular (see the last paragraph before Sect. 2.1.1) then R 0 is given by the Perron root of M f . With little effort, one can use the expected values provided in (6) to show that necessary and sufficient conditions for M f to be positively regular are that assumptions A1-A2 hold and that 0 < E(T s ) < 1 and 0 < E(T Δ ) < 1. If some of these conditions are not satisfied, we may analyse the spread of the disease by reducing the number of types of the approximating forward branching process. It is worth pointing out that R 0 only depends on the marginal distributions of T s and T Δ (via their moments), not on the dependence structure between them.

Probability of a major outbreak
For two s-dimensional vectorsā = (a 1 , . . . , a s Let f : [0, 1] 3 → R 3 be the probability generating function of the offspring distribution of the three types in the approximating branching process. That is, for whereξ i = (ξ i,1 , ξ i,2 , ξ i,3 ) is distributed as the offspring of a type i individual, i = 1, 2, 3.
Similarly, let f * : [0, 1] 3 → R be the probability generating function of the offspring distribution of the initial case. Ifξ = (ξ * ,1 , ξ * ,2 , ξ * ,3 ) T is distributed as the offspring of the initial case, then f * is given by For i = 1, 2, 3, let (S (i) , Δ (i) ) be the joint degree of a type i case with offspring (ξ i,1 , ξ i,2 , ξ i,3 ) and transmission weight T = (T s , T Δ ). That is, Here d = denotes equality in distribution. By conditional independence we have Conditioned on the transmission weight T and the single degree S (1) , ξ 1,3 has a binomial distribution with parameters S (1) and T s . Thus • ) is independent of T .
Since the conditional offspring distribution of a type 2 individual is identical to the offspring distribution of a type 1 individual except that a type 2 individual may give birth to one additional type 1 individual with probability T Δ , we have Similarly, . (13) Substituting (11)-(13) into (10) gives an expression for f . By standard branching process theory, if R 0 > 1 the extinction probability of a process descending from a type i individual, i = 1, 2, 3, is given by q i , wherē q = (q 1 , q 2 , q 3 ) T is the unique solution ofq = f (q) in [0, 1) 3 . We also havē where f •n is the composition of f with itself n times. Since the approximating branching process dies out if and only if each of the processes started by the children of the initial case die out, the probability of extinction is given by f * (q). After some calculations, analogous to the calculations that led to (11)-(13), we find that the probability of extinction is given by where (S, Δ) is independent of T . We conclude that, by Theorem 1, the probability of a major outbreak is given by 1 − f * (q), whereq is the limit in (14) and also the fixed point of f in [0, 1) 3 .

Backward process
Let w be a given node of G N , chosen uniformly at random. We use a backward branching process to approximate the probability that w contracts the disease, which by an exchangeability argument equals the expected final size of a major outbreak.
The offspring of an individual v in the backward process are the individuals that would potentially have infected v, if they were infected themselves. The members of the susceptibility set are divided into the following two groups (which give rise to a two-type approximating backward branching process).
Type 1: The vertex is included in the susceptibility set by virtue of potential transmission along a single edge We assign kinship as follows. The children of type 1 of an individual v 1 are the individuals included in the susceptibility set due to potential transmission along a single edge. The children of type 2 of v 1 are the individuals included in the susceptibility set due to potential transmission of the disease to v 1 , within a triangle of which v 1 is a member. We note that, given a triangle v 1 , v 2 , v 3 where v 1 is the primary case, both v 2 and v 3 will be members of the susceptibility set of v 1 by virtue of transmissions within the triangle if at least one of the following events happens: Here "infects" is conditional on the "infector" being infected during the epidemic.
The events E 1 -E 3 are illustrated in Fig. 5. Standard calculations give that the probability of the union of the events E 1 -E 3 is given by Similarly, the probability that neither v 1 nor v 2 will be members of the susceptibility set of v by transmissions within the triangle is given by p 0 = (1 − E(T Δ )) 2 . For later use, denote 1 − p 0 − p 2 by p 1 .

Expected final size of a major outbreak
Let b be the probability generating function of the offspring distribution of the two types of the approximating backward branching process. Furthermore, let b * be the probability generating function of the offspring distribution of the ancestor w. Analogously to the forward branching process, the probability that a branching population whose ancestor is of type i, i = 1, 2, will go extinct is given by . The probability of extinction is given by b * (q b ).
Proceeding in the same manner as in Sect. 3.1.1 yields where p 0 , p 1 and p 2 are as in Sect. 3.2. Similarly , and the probability of ultimate extinction of the backward process is given by We conclude that the expected final size of a major outbreak is given by 1 − b * (q b ).

Random vaccination with a perfect vaccine
Assume that a fraction f v < 1 of the population is vaccinated, and that the vaccinated individuals are chosen uniformly at random (without replacement) from the population. The vaccine is perfect, in the sense that a vaccinated individual gains full and lasting immunity to the disease. If the population size N is large, we may use a slightly different model, where each individual is vaccinated with probability f v , independently of the vaccination status of other individuals. By the law of large numbers, for our purposes the models are equivalent in the limit as the population size N → ∞.
As before, we may approximate the early phase of the epidemic by a multi-type branching process. The individuals of the approximating branching process are now of the following three types.
Type 1: Infected along a triangle edge and has a twin that is known not to be susceptible Type 2: Infected along a triangle edge and has a twin that might be susceptible Type 3: Infected along a single edge To clarify the types, assume that in the early phase of the epidemic v 1 is the primary case in the triangle v 1 , v 2 , v 3 . If v 1 attempts to transmit the disease both to v 2 and v 3 and succeeds (that is, none of v 2 and v 3 are vaccinated) then both v 2 and v 3 are represented by type 1 individuals in the approximating branching process. This happens with probability If v 1 attempts to transmit the disease both to v 2 and v 3 , but only succeeds to transmit the disease to v 3 (that is, v 2 is vaccinated and v 3 is not vaccinated), then in the approximating branching process the individual representing v 1 gives birth to one type 1 individual (representing v 3 ) within the triangle v 1 , v 2 , v 3 . This happens with probability If v 1 attempts to transmit the disease only to v 2 and succeeds (that is, v 2 is not vaccinated) then in the approximating branching process, the individual representing v 1 gives birth to one type 2 individual (representing v 2 ) within the triangle v 1 , v 2 , v 3 . Fig. 6 Three examples of transmission dynamics within a triangle v 1 , v 2 , v 3 . An attempted transmission of the disease is represented by an arrow, an attempted transmission to a vaccinated individual is represented by an arrow and a blue bar. Left: v 1 attempts to transmit the disease both to v 2 and v 3 , and succeeds. Both v 2 and v 3 are represented by type 1 individuals in the approximating branching process. Center: v 1 attempts to transmit the disease both to v 2 and v 3 , the transmission to v 2 is blocked since v 2 is vaccinated. Then v 3 is represented by a type 1 individual. Right: v 1 succeeds to transmit the disease to v 2 , but does not attempt to infect v 3 . Then v 2 is represented by a type 2 individual (colour figure online) This happens with probability The above described events are illustrated in Fig. 6. Denote the mean matrix of the approximating branching process by M Using the expressions in (15) and (16) gives the expected number of type 1 individuals produced by a type 1 individual where m 1,1 is an element of the mean matrix M f of the forward branching process presented in (9).
Proceeding in the same fashion, we obtain the elements of the mean matrix M i, j=1 of the branching process with random vaccination. It turns out that

It is readily verified that the Perron root of M
where r f is the Perron root of M f . Setting r (v) f to 1 in (18) and solving for f v yields the critical vaccination coverage f We conclude that, for this particular graph model, equality holds between the basic reproduction number R 0 and the perfect vaccine-associated reproduction number R V as defined in (7). Note that R 0 is based on a rank-based perspective of infection and not on "who-infected-whom.

Probability of a major outbreak
Let h be the probability generating function of the offspring distribution of the three types in our model including vaccination. As in Sect. 3.1.1, we use the probability generating function to approximate the probability of extinction of the epidemic. To this end, let (ζ i,1 , ζ i,2 , ζ i,3 ) be distributed as the offspring of a type i individual with transmission weight T , i = 1, 2, 3, and let (S (i) , Δ (i) ) be distributed as the joint degree of this individual. That is, and Note that (S (i) , Δ (i) ) and T are independent. By conditional independence E z Conditioned on the transmission weight T and the joint degree (S (1) , Δ (1) ), the number of attempted transmissions from a type 1 individual along single edges has a binomial distribution with parameters S (1) and T s , and each attempted transmission succeeds with probability (1 − f v ). Thus, Similarly, for a type 1 individual w with triangle degree Δ (1) , by conditioning on the number of attempted transmissions (in k i of the Δ (1) − 1 triangles that are not yet affected by the disease, w attempts to transmit the disease to i individuals, i = 0, 1, 2) and the vaccination status of the individuals contacted by w we obtain Combining (19) and (20) yields By noting that the offspring distribution of a type 2 individual is identical to the offspring distribution of a type 1 individual, except that a type 2 may give birth to one additional type 1 individual with probability Similarly, Combining these results yields the probability generating function h of the offspring distribution of a type 1, 2, 3 individual respectively. That is, h(z) 1 is given by (21), h(z) 2 is given by (22) and h(z) 3 is given by (23).
The probability generating function h * of the initial case is given by Δ) is distributed as the joint degree of the initial case and independent of T . The probability of extinction of the approximating branching process is given by h * (q (v) ), whereq (v) is given by the point in [0, 1] 3 closest to the origin that satisfiesq (v) = h(q (v) ). Thus, by Theorem 1 the probability of a major outbreak is 1 − h * (q (v) ).

The backward process
We now turn our attention to the backward process and final size of an epidemic in a population where a fraction f v is vaccinated with a perfect vaccine. To this end, we introduce the following three types, where individuals are classified by their vaccination status and the type of the edge along which they would transmit the disease if infected.
Type 1: Transmits along triangle edge, no information on vaccination status is available Type 2: Transmits along triangle edge and is known not to be vaccinated since it is successfully infected by its twin Type 3: Transmits along single edge, no information on vaccination status is available To clarify the types a bit more, let v 1 , v 2 , v 3 be a given triangle. At least one of v 2 and v 3 belongs to the susceptibility set of v 1 by virtue of potential transmissions within the triangle if some the following events, illustrated in Fig. 7, happens. Note that all cases infected by virtue of transmission within the triangle v 1 , v 2 , v 3 are attributed to v 1 .
(E 1 ) v 2 attempts to infect v 1 and v 3 attempts to infect v 2 , both succeed, and v 3 does not attempt to infect v 1 . Or the same thing might happen, with v 2 and v 3 interchanged. This results in one type 1 and one type 2 individual in the approximating branching process. If v 1 is represented by a type 1 or 3 individual this happens with probability if v 1 is represented by an individual of type 2 this happens with probability (E 2 ) Only one of v 2 and v 3 attempts to infect v 1 , and succeeds. The other node does not attempt to infect any node within the triangle. This results in one type 1 offspring. If v 1 is represented by an individual of type 1 or 3 this happens with probability if v 1 is represented by an individual of type 2 this happens with probability (E 3 ) v 2 and v 3 both attempt to infect v 1 and succeeds. This results in two type 1 individuals born in the approximating branching process. If v 1 is represented by an individual of type 1 or 3 this happens with probability if v 1 is represented by an individual of type 2 this happens with probability .
(E 4 ) v 2 attempts to infect v 1 and succeeds. The other node, v 3 , attempts to infect v 2 , but fails due to v 2 being vaccinated. The individual v 3 does not attempt to infect v 1 . In this scenario, v 2 belongs to the susceptibility set of v 1 . However, we do not include v 2 is the approximating branching process. This does not have any impact on the result of our analysis, since we are only interested in the probability of extinction of the backward process and v 2 does not produce any offspring in this process.

Expected final size
Let b (v) and b (v) * be the probability generating function of the offspring distribution of the three types of the approximating backward branching process and of the ancestor, respectively. Furthermore, letζ i = (ζ b i,1 , ζ b i,2 , ζ b i,3 ) be distributed as the offspring of a type i, i = 1, 2, 3, individual and denote by E s the conditional expectation given that the parent of ) be distributed as the offspring of the ancestor. Denote the extinction probability of a process descending from a type i individual by q b i , i = 1, 2, 3 and letq where, as before, (S (i) , Δ (i) ) is distributed as the joint degree of a type i individual, i = 1, 2, 3. Now By conditioning on the number of triangles k 2 in which an event of type E 3 occurs, the number of triangles k a 1 in which an event of type E 1 occurs, the number of triangles k b 1 in which an event of type E 4 occurs and the number of triangles k c 1 in which an event of type E 2 occurs we obtain Inserting the right hand sides of (26) and (27) in (25) gives 3 ) Similarly and E(z Combining these results yields the probability generating function of the offspring distribution of the three types; b (v) (z) 3 is given by (28) and b (v) (z) 2 is given by (29). By replacing (S (s) • , Δ (s) • ) in the right hand side of (28) by (S (Δ) Also by replacing (S (s) • ) in the right hand side of (28), but now by (S, Δ) we obtain the probability generating function b (v) * (z) of the offspring of the initial case. The expected final size of the epidemic, conditioned on that a major outbreak occurs, is given by

Numerical example
Consider for now the special case where T s = T Δ . With some abuse of notation we denote T s = T Δ by T , i.e. T = T s = T Δ is one-dimensional here. Under very general assumptions, increasing the heterogeneity in infectiousness leads to a decrease in the probability of a major outbreak, the expected final size and R 0 (Kuulasmaa 1982;Meester and Trapman 2011;Miller 2008), see also Ball (1985); Kenah and Robins (2007);Miller (2007). In particular, for a fixed (marginal) transmission probability E(T ), the probability of a major outbreak and the expected final size are maximized if T = E(T ) with probability 1 and minimized if P(T = 1) = E(T ) = 1 − P(T = 0). Similarly, for given E(T ), R 0 is maximized if T = E(T ) with probability 1 and minimized if P(T = 1) = E(T ) = 1 − P(T = 0). We illustrate this with the following two examples. In this first example we assume that T = T s = T Δ . Consider the three degree distributions 1. p(2, 1) = 1 2. p(4, 0) = 0.95 = 1 − p(2, 1) 3. p(0, 2) = 0.95 = 1 − p(2, 1).
That is, in all three degree distributions the total degree is 4 with probability 1. In addition, distribution 1 corresponds to a network where every node is member of exactly one triangle. Distribution 2 corresponds to a network where a node is not a member of any triangle with probability 0.95, while with probability 0.05 a node is member of one triangle. Finally, distribution 3 corresponds to a network where a node is a member of two triangles with probability 0.95, while with probability 0.05 a node is member of one triangle.
Furthermore, let T have distribution Beta(α, α) for some α > 0. That is, T has density, C α x α−1 (1− x) α−1 , on the interval (0, 1), where C α is a normalizing constant. Then E(T ) = 1/2 and we can tune the heterogeneity of the infectivity of infected individuals by varying α. In particular Note that α = 1 corresponds to T ∼ U (0, 1), with α → ∞ corresponding to T becoming a point mass at 1/2. Here U (0, 1) denotes the uniform distribution on (0, 1). Figure 8 shows the probability that a major outbreak does not occur, the expected final size, R 0 and the critical vaccination coverage f (c) v as functions of α or E(T 2 ). As can be seen in Fig. 8, ignoring actual heterogeneity of infectivity in this case leads to an overestimation of the probability of a major outbreak (Fig. 8a, b). This effect is particularly evident in the presence of high clustering; the steeper slope of the curve corresponding to distribution 3 (Fig. 8b) and the relatively low probability of a major outbreak when α is small can be explained by the fact that the approximating forward branching process is close to being critical when α is small. Figure 8c, d shows that heterogeneity of infectivity has virtually no impact on the expected final size of a major outbreak and R 0 in the near absence of clustering, which is in line with known results for unclustered networks [see, for instance, section 4 in Miller (2008)]. In the presence of clustering, on the other hand, ignoring heterogeneity of infectivity leads to an underestimation of the expected final size and a substantial overestimation of the critical vaccination coverage f (c) v . Note that R 0 and f (c) v depend on the distribution of T only through the first and second moment of T .
Next, we relax the assumption T s = T Δ and investigate the impact of the correlation ρ between T s and T Δ on the spread of the disease. To this end, we consider a model where as before T s and T Δ both have distribution Beta(α, α) and where the correlation ρ = ρ(t) may be tuned by varying t ∈ [−1, 1]. Here ρ(t) is increasing in t with ρ(−1) = −1 and ρ(1) = 1. The degree distribution of the underlying graph is given by distribution 1 above.
To construct such a model, let the joint distribution of T s and T Δ be as follows. Let N 1 , N 2 , N 3 be three independent standard normal random variables, and assume that the joint distribution of T s and T Δ is given by where F α and Φ are the CDF's of the distributions Beta(α, α) and N (0, 1), and N s and N Δ are the standard normal random variables With little effort one can show that ρ(t) is indeed increasing in t, that (by the symmetry of the distribution Beta(α, α)) ρ(−1) = −1 and ρ(1) = 1, and that T s and T Δ are independent for t = 0. As can be seen in Fig. 9 the probability 1 − f * (q) that a major outbreak occurs increases as the correlation ρ(t) decreases. This effect is substantial when heterogeneity in individual infectivity is high (i.e. α is small) but wanes as the heterogeneity decreases (i.e. α increases). This can be explained by the fact that for a fixed value of α the probability that an infectious individual will not transmit the disease to any of its susceptible neighbours decreases as t increases, and that this probability is more sensitive to changes in t if the heterogeneity in individual infectivity is higher [cf. Kuulasmaa (1982)]. It should be noted that R 0 and the critical vaccination coverage f (c) v do not depend on the correlation ρ(t), which can be seen from the expression for the mean matrix in (9) and (18).

Discussion
In this paper, we have incorporated clustering in the spread of an infectious disease by allowing for groups of size three with non-overlapping edges. It is, in principle, straightforward to extend the methods used in this paper to larger group sizes. The CMC may, for instance, be generalized to larger group sizes as follows. Let K = {k 1 , . . . , k r } ⊂ N ≥2 be the set of possible group sizes. In the matching procedure, each node is equipped with an r -dimensional degree in N r 0 . The ith component (the k i -degree) of a degree specifies the number of groups of size k i to which the node in question belongs. Analogously to the construction of a CMC graph, groups are then formed by creating one list for each group size; a node with k i -degree d i appears precisely d i times in the list corresponding to groups of size k i . The lists are then shuffled and half-edges of nodes in positions k + 1, . . . , k + k i in the k i -list are joined. The structure of a graph so obtained would be characterized by fully connected cliques, and similar to that of a random intersection graph (Ball et al. 2014). One possible approach to investigate epidemics on such graphs would be to approximate the spread Fig. 8 The impact of heterogeneity in infectivity for the three degree distributions. With some abuse of notation we write T = T s = T Δ . a The probability that a major outbreak does not occur as a function of α. b The probability that a major outbreak does not occur as a function of E(T 2 ). c The expected final size of a major outbreak as a function of α. d The expected final size of a major outbreak as a function of E(T 2 ). e The basic reproduction number R 0 as a function of E(T 2 ). f The critical vaccination coverage f (c) v as a function of E(T 2 ) of the disease by a multitype Galton Watson process where groups (or cliques) are represented by the particles of the branching process. The types of the approximating branching process would then be vectors in N 2 of the form (m, n), where m represents the size of the clique and n represents the number of members of the clique that the primary case of the clique attempts to infect. Another possible approach would be to use an infinite type branching process in the spirit of Ball et al. (2014). We believe that the result would be analogous to the results obtained in Ball et al. (2014).

Fig. 9
The probability that a major outbreak does not occur as a function of α and t Acknowledgements This research is supported by the Swedish Research Council (Vetenskapsrådet) Grant 2016-04566. We wish to thank the members of the journal club on infectious diseases at Stockholm University and Daniel Ahlberg for suggestions that lead to substantial improvements of the paper.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix: Proof of proposition 1
Letd = {(S i , Δ i )} i∈N be a given (i.e. non-random) degree sequence that satisfies the following regularity assumptions.
where (S, Δ) has distribution p, which is assumed to satisfy A1-A2 in Sect. 2.1. Let further G = {G N } N ∈N be a sequence of graphs generated by the CMC, where the degree sequence of G N is given byd Under the assumptions A1-A2 the expected number of self-loops and the expected number of multiple edges are both of order O(1) [cf. Van der Hofstad (2016, prop. 7.11)]. Denote by A N the number of wedges of G N that are "deleted" when merging multiple edges and erasing self-loops, that is From the definition of A N , we deduce that the total number of ordered triangles of G N is bounded from below by |W G N Δ | ≥ N i=1 2Δ i − A N and the total number of ordered wedges is bounded from above by Therefore, by the definition of C(G N ) and the assumptions above as N → ∞. This lower bound is tight in the limit as the number of nodes N → ∞. Indeed, denote by W G N s the set of ordered triangles of G N that consists solely of single edges, i.e.
where the sums run over the integers 1, . . . , N . Dividing by N in (32) and letting N approach infinity gives E(|W G N S |)/N → 0 as N → ∞. Thus |W G N S |/N → 0 in probability. Repeating this procedure for triangles formed by a combination of triangle and single edges gives The assertion now follows by bounded convergence and the law of large numbers.