  * Date: Oct 2021

# Vaccination with partial transmission and social distancing on contact
networks

Christian John Hurry1, Alexander Mozeika2, Alessia Annibale1,3 1Department of
Mathematics, King’s College London, Strand, WC2R 2LS
2School of Cancer and Pharmaceutical Sciences, King’s College London, London,
UK
3Institute for Mathematical & Molecular Biomedicine, Hodgkin Building, Guy’s
Campus, London, SE1 1UL<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
###### Abstract

We study the impact of vaccination on the risk of epidemics spreading through
structured networks using the cavity method of statistical physics. We relax
the assumption that vaccination prevents all transmission of a disease used in
previous studies, such that vaccinated nodes have a small probability of
transmission. To do so we extend the cavity method to study networks where
nodes have heterogeneous transmissibility. We find that vaccination with
partial transmission still provides herd immunity and show how the herd
immunity threshold depends upon the assortativity between nodes of different
transmissibility. In addition, we study the impact of social distancing via
bond percolation and show that percolation targeting links between nodes of
high transmissibility can reduce the risk of an epidemic greater than
targeting links between nodes of high degree. Finally, we extend recent
methods to compute the distributional equations of risk in populations with
heterogeneous transmissibility and show how targeted social distancing
measures may reduce overall risk greater than untargeted vaccination
campaigns, by comparing the effect of random and targeted strategies of node
and link deletion on the risk distribution.

††:

Keywords: Epidemic modelling, Cavity and replica method, Random graphs,
networks

##  1 Introduction

The study of compartmental models, where a population is split into, for
example, ‘Susceptible’, ‘Infected’ and ‘Recovered’ compartments in the SIR
model, have formed an important theoretical and computational basis for the
study of epidemics and design of vaccination campaigns. An increased
understanding of contact networks, which detail physical contacts of a
sustained duration between individuals in a population (see e.g. [1]), has
been complemented by mathematical results for compartmental epidemic models on
networks. One such contribution is the prediction of an epidemic threshold
dependent upon the rates of infection and recovery, as well as the contact
network topology [2, 3, 4]. Many important results concerning the SIR model on
networks have been derived by message-passing approaches, also known as the
cavity method, including the size and risk of epidemics on networks with
arbitrary degree distributions [5] and degree correlations [6, 7]. More
recently, these methods have been employed to model epidemic mitigation via
contact tracing apps [8] and competing strains of infectious diseases [9, 10,
11]. In addition to analytical results, the cavity method provides a set of
equations that allows for an efficient _parallel_ numerical implementation, as
opposed to direct simulations of the SIR model, which are usually _sequential_
and typically marred by long computation times, scaling with the size and
connectivity of the network and the infection rate relative to the recovery
rate.

One of the key problems that mathematical epidemiology addressed is the
optimisation of vaccination strategies, i.e given a finite supply of vaccines,
who should be vaccinated to mitigate an epidemic. Previous work has studied
vaccination strategies based upon network topology. Vaccination which
prioritises nodes of higher degree, i.e. people with a high number of social
contacts, is known to lead to better outcomes in comparison with a random
vaccination campaign [12]. Furthermore, it has been shown that vaccination
using information beyond the degree can improve upon degree-based strategies
[13, 14]. The recent discovery of the role that the Hashimoto non-backtracking
matrix [15, 16, 17] plays in epidemic processes has led to the use of
eigenvalues to rank nodes for prioritisation of vaccination [18]. The benefit
of modelling such strategies is that they provide principles upon which
vaccination campaigns can be based with only partial knowledge of the contact
network. For example, it can provide theoretical insight into prioritising
vaccination for parts of the population with higher than average social
connectivity.

Vaccination is usually modelled by placing an individual into a separate
vaccinated compartment, such that vaccinated nodes block incoming infections,
so that if enough people are vaccinated the infection can no longer spread
through the population, a phenomena referred to as herd immunity. This assumes
that a vaccine provides full protection against transmission of an infectious
disease for any individual. In general, it may be desirable to relax this
assumption. For example, although the efficacy of a vaccine in preventing
symptoms is ascertained before approval for public use, it is more difficult
to determine how a given vaccine prevents transmission until data is collected
during or after vaccine roll-out. In this case, it may be desirable to assume
that vaccination reduces the transmissibility of an individual to a small but
finite value. Furthermore, it is also important to recognise that vaccines are
usually given to priority groups first, such as people with underlying health
conditions, medical staff, or people above a certain age. Due to the
correlated nature of social contacts [19, 20], prioritising vaccinations may
lead to correlations between the vaccinated status of a node and its
topological properties. Variation in the transmissibility of individuals has
been explored in previous works [3, 21, 22, 23], however, analytical results
have been restricted to the average risk or size of the epidemic. It has
previously been shown that there is high variability in the risk of individual
nodes, due to differences in the node environment [17, 24], which is neglected
when only considering the average risk. In our work, we show that this
variability is more prominent when transmissibility varies between individuals
and focus on how this impacts the distribution of risk.

After reviewing node percolation approaches to vaccination, we focus our study
on vaccination with partial transmission by extending the cavity approach to
allow for _heterogeneity_ in the transmission between individuals. At first we
consider transmissibility and degree to be uncorrelated, and we derive the
herd immunity threshold in this scenario, showing that vaccination with
partial transmission will always require a greater proportion of the
population vaccinated to achieve herd immunity. We then relax the assumption
that transmissibility and degree are uncorrelated to show how these
correlations affect epidemic risk. Using this framework, we then study the
impact of social distancing between groups of different transmissibility, by
deriving equations for the risk of an epidemic under link percolation that
targets either links between nodes of high transmissibility, or nodes of high
degree. To go beyond the average risk of an epidemic, we follow techniques
developed in context of the micro-structure of percolation in the
configuration model [24] to investigate the distribution of risk in structured
and heterogeneous populations. Our analysis reveals a highly non-trivial
distribution of risk, even amongst nodes of the same degree and
transmissability. Finally, we extend the cavity approach for the
distributional equations of risk to account for node and link deletion. We
comment that this provides a succinct procedure to explore the impact of
targeted vaccination and social distancing strategies on the distribution of
risk.

The remaining sections of this article are organised as follows. In Sec. 2 we
review the network percolation approach to the study of the steady state
reached in the SIR model on networks with arbitrary degree distributions and
degree correlations. In Sec. 3 we extend the cavity method to account for a
network with nodes assigned to sub-types of different transmissibility, and
provide a closed set of equations for the average risk a node poses to the
network. We consider the cases where node transmissibility is and is not
correlated with social contact, separately. In Sec. 4 we derive equations for
the average risk under link percolation, and show how the choice of the links
which are deleted affects the risk. Sec. 5 demonstrates how the cavity method
may be used to derive the distributional equations of risk and how to solve
these equations via a population dynamics procedure. We conclude with a
discussion of the theoretical value of our results and potential avenues for
future work. Technical details are described in the appendices.

##  2 Impact of vaccination on the epidemic risk in contact networks

One of the main interests in the study of the risk of epidemics spreading on
contact networks is the exploration of vaccination strategies. A question that
arises, when assessing vaccination strategies, is what fraction of the
population needs to be vaccinated to achieve herd immunity. In this section we
consider the SIR epidemic model on contact networks and we show that this
question can be answered by studying a simple node percolation problem.

###  2.1 Herd immunity with perfect vaccination in the SIR model on structured
networks

We consider the SIR model on an undirected network with $N$ nodes and
adjacency matrix ${\bf A}$. The elements
$A_{ij}=A_{ji}\in\\{0,1\\}~{}\forall~{}i,j$ denote presence ($A_{ij}=1$) or
absence ($A_{ij}=0$) of link for each pair of nodes $(i,j)$ and
$A_{ii}=0~{}\forall~{}i$. Following the formulation in [17], we assume that an
infected node may pass an infection to a neighbouring node in a short time
interval $\rmd t$ with probability $\beta\rmd t$ and an infected node recovers
from an infection at a time drawn from an ‘infectious time’ distribution,
$\gamma(t)$, normalised such that $\int_{0}^{\infty}\rmd t\gamma(t)=1$, to
ensure that an infected node eventually recovers. In order to quantify the
fraction of the population that needs to be vaccinated to achieve herd
immunity, we note that if a vaccine provides perfect immunity, its action on
an individual is functionally equivalent to deleting a node in the network, as
this prevents the spread of infection through that node and its links. In
order to incorporate the effect of node deletion we define a binary variable
$\sigma_{i}\in\\{0,1\\}$ which describes whether a node has been vaccinated
($\sigma_{i}=0$) or not ($\sigma_{i}=1$). As detailed in [17], the probability
$r_{i}$ that an infection starting from a single node $i$ causes an outbreak
across the bulk of the network, also called the _risk_ of node $i$, can be
determined, via a message-passing approach. Using the same approach and
incorporating the effect of node deletion, we can write

$\displaystyle r_{i}({\bf A})$
$\displaystyle=\sigma_{i}\left[1-\int_{0}^{\infty}\rmd
t\,\gamma(t)\prod_{j\in\partial^{\bf A}_{i}}\left(1-\left(1-\rme^{-\beta
t}\right)r_{j}^{(i)}({\bf A})\right)\right]$ (1) $\displaystyle
r_{j}^{(i)}({\bf A})$ $\displaystyle=\sigma_{j}\left[1-\int_{0}^{\infty}\rmd
t\,\gamma(t)\prod_{\ell\in\partial^{\bf A}_{j}\setminus
i}\left(1-\left(1-\rme^{-\beta t}\right)r_{\ell}^{(j)}({\bf A})\right)\right]$
(2)

where $\partial^{\bf A}_{i}=\\{j:A_{ij}=1\\}$ is the set of the neighbours of
$i$ and $r_{j}^{(i)}({\bf A})$ is the risk that node $j$ causes an outbreak in
the ‘cavity’ graph where node $i$ is removed. These equations are exact on
trees, where loops are absent, and become exact on locally tree-like graphs,
in the limit of large system size, where the typical length of loops diverges
logarithmically [25]. For a given contact network, one can solve the cavity
equations (1) and (2) numerically. Alternatively, in absence of true knowledge
of the contact network one can assume that it is _random_ and use equations
(1) and (2) to derive equations for the global risk $g({\bf
A})=\frac{1}{N}\sum_{i=1}^{N}r_{i}({\bf A})$, which is expected to be _self-
averaging_ [26] when $N\rightarrow\infty$ and thus independent of its
microscopic details. Averaging (1) over all sites as shown in A, one obtains,
for a single graph instance ${\bf A}$

$\displaystyle g({\bf A})$ $\displaystyle=$
$\displaystyle\sum_{\sigma}\sigma\left[\mathrm{P}(\sigma|{\bf
A})-\sum_{k}\mathrm{P}(k,\sigma|{\bf A})\int_{0}^{\infty}\rmd
t\,\gamma(t)\left(1-\alpha(t)\hat{g}_{k}({\bf A})\right)^{k}\right]$ (3)
$\displaystyle\hat{g}_{k}({\bf A})$ $\displaystyle=$
$\displaystyle\sum_{\sigma^{\prime}}\sigma^{\prime}\left[\mathrm{P}(\sigma^{\prime}|{\bf
A})-\sum_{k^{\prime}\geq 1}\frac{\mathrm{W}(k;k^{\prime},\sigma^{\prime}|{\bf
A})}{\mathrm{W}(k|{\bf A})}\int_{0}^{\infty}\rmd
t\,\gamma(t)\left(1-\alpha(t)\hat{g}_{k^{\prime}}({\bf
A})\right)^{k^{\prime}-1}\right]$

where we have denoted for brevity $1-\rme^{-\beta t}=\alpha(t)$. We have
defined the likelihood to draw at random a node with degree $k$ and label
$\sigma$ for the graph instance ${\bf A}$ as $\mathrm{P}(k,\sigma|{\bf
A})=N^{-1}\sum_{i}\delta_{\sigma,\sigma_{i}({\bf A})}\delta_{k,k_{i}({\bf
A})}$, with $k_{i}({\bf A})=|\partial^{\bf A}_{i}|$ denoting the degree of
node $i$ and $\mathrm{P}(\sigma|{\bf A})=\sum_{k}\mathrm{P}(k,\sigma|{\bf A})$
the marginal distribution. Furthermore, we have defined
$\mathrm{W}(k;k^{\prime},\sigma^{\prime}|{\bf
A})=\sum_{ij}A_{ij}\delta_{k,k_{i}({\bf A})}\delta_{k^{\prime},k_{j}({\bf
A})}\delta_{\sigma^{\prime},\sigma_{j}({\bf A})}/N\bar{k}({\bf A})$ as the
likelihood that by drawing a link at random we choose a link with a node of
degree $k$ at one end and a node with degree $k^{\prime}$ and label
$\sigma^{\prime}$ at the other. We denote by $\bar{k}({\bf
A})=N^{-1}\sum_{i}k_{i}({\bf A})$ the mean degree. The distribution
$\mathrm{W}(k;k^{\prime}|{\bf
A})=\sum_{\sigma}\mathrm{W}(k;k^{\prime},\sigma|{\bf A})$ is known as the
degree correlations, and $\mathrm{W}(k|{\bf
A})=\sum_{k^{\prime}}\mathrm{W}(k;k^{\prime}|{\bf A})$ is its marginal
distribution. Here and below we adopt the convention to denote a joint
probability distribution of node quantities across connected node pairs with
$\mathrm{W}(\dots;\dots)$. We will consider vaccination strategies based upon
the degree of nodes, such that $\mathrm{P}(k,\sigma|{\bf A})=\mathrm{P}(k|{\bf
A})\mathrm{P}(\sigma|k)$ and $\mathrm{W}(k;k^{\prime},\sigma^{\prime}|{\bf
A})=\mathrm{W}(k;k^{\prime}|{\bf A})\mathrm{P}(\sigma^{\prime}|k^{\prime})$,
where $\mathrm{P}(\sigma|k)$ models the degree-dependent vaccination strategy.

In the limit of large networks, $N\to\infty$, we expect fluctuations of
intensive network observables to vanish and their value on a single graph
realisation to coincide with their ensemble averaged value. Assuming that the
degree distribution and degree correlations display such self-averaging
behavior in the limit of large $N$, we average the cavity equations over a
suitably defined random graph ensemble, obtaining

$\displaystyle g$
$\displaystyle=\sum_{\sigma}\sigma\left[\mathrm{P}(\sigma)-\sum_{k}\mathrm{P}(k,\sigma)\int_{0}^{\infty}\rmd
t\,\gamma(t)\left(1-\alpha(t)\hat{g}_{k}\right)^{k}\right]$ (4)
$\displaystyle\hat{g}_{k}$ $\displaystyle=\sum_{k^{\prime}\geq
1}\sum_{\sigma^{\prime}}\mathrm{W}(k^{\prime}|k)\mathrm{P}(\sigma^{\prime}|k^{\prime})\sigma^{\prime}\left[1-\int_{0}^{\infty}\rmd
t\,\gamma(t)\left(1-\alpha(t)\hat{g}_{k^{\prime}}\right)^{k^{\prime}-1}\right]$
(5)

where we have introduced the ensemble averages
$\mathrm{P}(\sigma)={\langle}\mathrm{P}(\sigma|{\bf A}){\rangle}_{\bf A}$,
$\mathrm{P}(k,\sigma)={\langle}\mathrm{P}(k,\sigma|{\bf A}){\rangle}_{\bf A}$,
$\mathrm{W}(k^{\prime}|k^{\prime})={\langle}\mathrm{W}(k^{\prime}|k,{\bf
A}){\rangle}_{\bf A}$ and ${\langle}k{\rangle}={\langle}\bar{k}({\bf
A}){\rangle}_{\bf A}$, with ${\langle}\cdot{\rangle}_{\bf A}=\sum_{\bf
A}\cdot\mathrm{P}({\bf A})$ and $\mathrm{P}({\bf A})$ the probability over the
set of symmetric adjacency matrices ${\bf A}\in\\{0,1\\}^{N(N-1)/2}$, which
defines the random graph ensemble.

The set of equations (5) may be solved numerically. It is easy to check that
they always have a trivial solution $\hat{g}_{k}=0\,\forall~{}k$,
corresponding to the absence of an epidemic outbreak, i.e. $g=0$. The
stability of this solution depends upon the Jacobian of the system of
equations (5), such that an epidemic will occur if $|\lambda_{1}^{M}|>1$ where
$\lambda_{1}^{M}$ is the largest eigenvalue of the matrix with entries

$\displaystyle
M_{k,k^{\prime}}=(k^{\prime}-1)\mathrm{W}(k^{\prime}|k)\sum_{\sigma^{\prime}}\sigma^{\prime}\mathrm{P}(\sigma^{\prime}|k^{\prime})T{\color[rgb]{0,0,1},}$
(6)

where $T=\int_{0}^{\infty}\rmd t\,\gamma(t)\alpha(t)$ is known as the
transmissability of the nodes [5]. The above result was first derived in [6],
in absence of node deletion. For uncorrelated networks,
$\mathrm{W}(k^{\prime}|k)=\mathrm{W}(k^{\prime})$ and one has
$\hat{g}_{k}=\hat{g}~{}\forall~{}k\geq 1$ with

$\displaystyle\hat{g}$
$\displaystyle=\sum_{\sigma}\sigma\left[\mathrm{P}(\sigma)-\sum_{k\geq
1}\mathrm{P}(\sigma|k)\frac{k\mathrm{P}(k)}{{{\langle}k{\rangle}}}\int_{0}^{\infty}\rmd
t\,\gamma(t)\left(1-\alpha(t)\hat{g}\right)^{k-1}\right]\equiv f(\hat{g})$

Due to the monotonic (i.e. $f^{\prime}(\hat{g})>0$) and concave (i.e.
$f^{\prime\prime}(\hat{g})<0$) nature of $f(\hat{g})$, a non-trivial solution
exists, and is stable, if $f^{\prime}(0)>1$. Evaluating this gives a condition
for an epidemic to occur,

$\displaystyle
T>\frac{{{\langle}k{\rangle}}}{\sum_{k,\sigma}\sigma\mathrm{P}(\sigma|k)k(k-1)\mathrm{P}(k)}.$
(7)

Focusing on the simplest case of a random vaccination campaign, where the
variables $\\{\sigma_{i}\\}_{i=1}^{N}$ do not depend on the network degrees
and are random i.i.d. with distribution

$\displaystyle\mathrm{P}(\sigma)=p\delta_{\sigma,1}+(1-p)\delta_{\sigma,0},$
(8)

where $(1-p)\in[0,1]$ describes the fraction of nodes that are vaccinated, one
finds a critical value of $p$ below which epidemics are prevented,

$\displaystyle
p<p_{c}=\frac{{{\langle}k{\rangle}}}{T{\langle}k(k-1){\rangle}}.$ (9)

This is the herd immunity threshold, a rearranged form of the percolation
threshold in the configuration model [27].

Figure 1:  Global risk as a function of the infection rate $\beta$, for nodes
with exponential infectious time distribution. Left: Results shown for an
Erdös-Rényi (ER) graph with average connectivity ${{\langle}k{\rangle}}=5$ and
nodes with mean infectious time $\xi^{-1}=1/2$. Solid lines indicate
theoretical predictions from the cavity equations (5) while symbols show
results from simulations on an ER graph of size $N=3000$, averaged over $100$
initial sites of infection, and $20$ runs starting from each site. The dashed
line indicating the epidemic threshold is computed from equation (9), which
simplifies to $p_{c}=1/T{{\langle}k{\rangle}}$ for ER graphs. Error bars show
the standard deviation of the risk across different initial sites of
infection. Right: Theoretical predictions from equation (5) are shown for
graphs with Poissonian degree distribution with ${{\langle}k{\rangle}}=3$ and
nodes with mean infectious time $\xi^{-1}=1.25$. Results are shown for graphs
with neutral, assortative and disassortative degree correlations. Dashed lines
indicate the predicted epidemic threshold from the largest eigenvalue of (6)

In Figure 1 (left panel) the global risk predicted from the cavity equation
(4) is plotted for an Erdös-Rényi (ER) graph with mean degree
${{\langle}k{\rangle}}=5$ and shown to be in excellent agreement with results
from simulations [28, 29, 30] of the stochastic SIR model on networks with
$N=3000$ nodes. We note that the wide error bars, showing the standard
deviation of the risk across different initial sites of infection, indicate
that there is large variation in the risk of individual nodes in a given ER
network. The impact of degree correlations is shown on the right panel of
Figure 1, where results from the cavity equations are shown for networks with
the same average connectivity and Poissonian degree distribution, but
assortative and disassortative degree mixing, respectively. To compute the
degree correlation function $\mathrm{W}(k,k^{\prime})$ for assortative and
disassortative Poissonian graphs we followed the methods described in [31].
These results are consistent with results for the size of the giant component
in networks with degree correlations [32]. In Figure 2 we plot the global risk
as a function of the fraction $1-p$ of vaccinated population in ER graphs with
average connectivity ${{\langle}k{\rangle}}=5$. In the left panel we see that
predictions from the cavity equations (5) are in agreement with results from
simulations on networks of size $N=1000$. The right panel shows that the herd
immunity threshold is lowest in Poissonian graphs with disassortative
correlations, and highest for ER graphs with assortative correlations.

Figure 2: Global risk plotted as a function of the fraction $1-p$ of
vaccinated nodes with infection rate $\beta=0.5$ and exponential distribution
of infectious times with parameter $\xi=0.4$. Left: For an Erdös-Rényi graph
with ${{\langle}k{\rangle}}=5$ we show predictions (solid line) from the
cavity equations (5) and results of simulations (circles) on graphs with size
$N=1000$, where nodes were selected with probability $1-p$ to be vaccinated,
averaged over $10$ configurations of vaccinated nodes, $25$ initial sites of
infection, and $25$ runs. The dashed line indicates the herd immunity as
computed from equation (9). Right: Theoretical results from the cavity
equations (5) are shown for graphs with Poissonian degree distribution, mean
degree ${{\langle}k{\rangle}}=3$, for neutral, assortative and disassortative
degree correlations

##  3 Impact of vaccination with partial transmission

In the previous section we have considered an idealised vaccine that prevented
all infection passing through paths containing vaccinated nodes. In general,
vaccinations may not prevent all transmission, and individuals may respond to
vaccines differently, leading to heterogeneities in individual infectious
potential. It is not possible to know _a priori_ how an individual will
respond to a vaccine, however it is feasible to gather data on the impact of
vaccination on transmissibility within a population and estimate differences
in transmission between different demographics. It is therefore of interest to
understand how variations in coverage and infectious potential of individuals
can affect vaccination strategies. In order to do so, we extend the cavity
approach to the SIR model to account for heterogeneity in the transmissibility
of individuals. We then go on to show how this affects the herd immunity
threshold.

###  3.1 Epidemic risk with heterogeneous transmissibility

In order to account for heterogeneity in the transmissibility of individuals,
we assume that individuals remain infectious for different times. Without loss
of generality we assume that the infectious time distribution of each
individual has the same functional form but its scale varies between
individuals of different sub-type. Therefore, we introduce $M$ node sub-types,
identified by labels $\xi\in\\{\xi^{1},...,\xi^{M}\\}$, and assign each node a
label, $\\{\xi_{i}\\}_{i=1}^{N}$, according to some probability distribution
$\mathrm{P}(\xi)$. The labels $\\{\xi^{1},...,\xi^{M}\\}$ take values in
$(0,\infty)$, such that $\xi$ parameterises the infectious time distribution
$\gamma(t|\xi)$ of nodes with label $\xi$. For the remainder of this paper we
choose $\gamma(t|\xi)=\xi\rme^{-\xi t}$ such that $\frac{1}{\xi}$ represents
the _mean infectious time_ of individuals with label $\xi$. This particular
choice for the recovery time distribution is known as Markovian recovery,
however, our approach also holds for non-Markovian recovery times, which could
be implemented by e.g. a Weibull distribution [33].

In this case, the risk of node $i$ depends on the label of node $i$, and all
labels downstream of node $i$. We denote the set of labels downstream of node
$i$, including node $i$, by $\boldsymbol{\xi}_{i}$. The equations for the
local risk for a given instance of a contact network ${\bf A}$ with prescribed
labels $\boldsymbol{\xi}=(\xi_{1},\ldots,\xi_{N})$ read as

$\displaystyle r_{i}({\bf A},\boldsymbol{\xi}_{i})$
$\displaystyle=1-\int_{0}^{\infty}\rmd
t\,\gamma(t|\xi_{i})\prod_{j\in\partial_{i}^{\bf
A}}\left(1-\alpha(t)r_{j}^{(i)}({\bf A},\boldsymbol{\xi}_{j})\right)$ (10)
$\displaystyle r_{j}^{(i)}({\bf A},\boldsymbol{\xi}_{j})$
$\displaystyle=1-\int_{0}^{\infty}\rmd
t\,\gamma(t|\xi_{j})\prod_{\ell\in\partial_{j}^{\bf A}\setminus
i}\left(1-\alpha(t)r_{\ell}^{(j)}({\bf A},\boldsymbol{\xi}_{\ell})\right).$
(11)

Averaging over the nodes, we find the equation for the global risk, which
depends on the whole label sequence $\boldsymbol{\xi}$ and the network ${\bf
A}$

$\displaystyle g({\bf A},\boldsymbol{\xi})=\frac{1}{N}\sum_{i=1}^{N}r_{i}({\bf
A},\boldsymbol{\xi}).$ (12)

In B we show that the assumption of a random tree-like structure allows us to
find a closed expression for the global risk in the $N\to\infty$ limit, that
once averaged over the graph and label ensemble, reads as

$\displaystyle g$
$\displaystyle=1-\sum_{k,\xi}\mathrm{P}(k,\xi)\int_{0}^{\infty}\rmd
t\gamma(t|\xi)\left(1-\alpha(t)\hat{g}_{k,\xi}\right)^{k}$ (13)
$\displaystyle\hat{g}_{k,\xi}$
$\displaystyle=\sum_{k^{\prime},\xi^{\prime}}\mathrm{W}(k^{\prime},\xi^{\prime}|k,\xi)\left[1-\int_{0}^{\infty}\rmd
t\gamma(t|\xi^{\prime})\left(1-\alpha(t)\hat{g}_{k^{\prime},\xi^{\prime}}\right)^{k^{\prime}-1}\right]$
(14)

where $\mathrm{P}(k,\xi)$ and $\mathrm{W}(k^{\prime},\xi^{\prime}|k,\xi)$ are
the ensemble-averaged values of $\mathrm{P}(k,\xi|{\bf
A},\boldsymbol{\xi})=\frac{1}{N}\sum_{i}\delta_{k,k_{i}({\bf
A})}\delta_{\xi,\xi_{i}}$ and $\mathrm{W}(k^{\prime},\xi^{\prime}|k,\xi,{\bf
A},\boldsymbol{\xi})=\mathrm{W}(k,\xi;k^{\prime},\xi^{\prime}|{\bf
A},\boldsymbol{\xi})/\mathrm{W}(k,\xi|{\bf A},\boldsymbol{\xi})$,
respectively, where $\mathrm{W}(k,\xi;k^{\prime},\xi^{\prime}|{\bf
A},\boldsymbol{\xi})=\frac{1}{N\bar{k}({\bf
A})}\sum_{ij}A_{ij}\delta_{k,k_{i}({\bf A})}\delta_{k^{\prime},k_{j}({\bf
A})}\delta_{\xi,\xi_{i}}\delta_{\xi^{\prime},\xi_{j}}$. For a given kernel
$\mathrm{W}(k^{\prime},\xi^{\prime}|k,\xi)$ of a network with maximum degree
$k_{\max}$ one may solve the system of $M\times k_{\max}$ equations (14)
numerically and substitute the result into (13) to find the global risk.

Alternatively, one may be able to simplify these equations for specific
choices of $\mathrm{W}(k,\xi;k^{\prime},\xi^{\prime})$. For example, if we
assume the transmissibility of an individual is uncorrelated with its degree,
the degree of its neighbours and the transmissibility of its neighbour, in
this case
$\mathrm{W}(k,\xi;k^{\prime},\xi^{\prime})=\mathrm{W}(k;k^{\prime})\mathrm{P}(\xi)\mathrm{P}(\xi^{\prime})$,
and upon substitution of this kernel into the r.h.s of (14) we see that
${\hat{g}}_{k,\xi}={\hat{g}}_{k}\forall\xi$, and the system of equations (14)
reduce to,

$\displaystyle\hat{g}_{k}$
$\displaystyle=\sum_{k^{\prime}}\mathrm{W}(k^{\prime}|k)\left[1-\int_{0}^{\infty}\rmd
t\left<\gamma(t|\xi)\right>_{\xi}\left(1-\alpha(t)\hat{g}_{k^{\prime}}\right)^{k^{\prime}-1}\right]$
(15)

where $\left<...\right>_{\xi}=\sum_{\xi}...\mathrm{P}(\xi)$. We see that
equation (15) takes the same form as (5) but now depends upon the ‘average’
infectious time distribution $\left<\gamma(t|\xi)\right>_{\xi}$. The stability
criteria is therefore identical to that of (5) but now depends upon the
average transmissibility, such that an epidemic will occur if

$\displaystyle|\lambda_{1}^{(\bf{J})}|>1$ (16)

where $\lambda_{1}^{(\bf{J})}$ is the largest eigenvalue of the Jacobian of
(15) with entries
$J_{k,k^{\prime}}=\left<T(\xi)\right>_{\xi}(k^{\prime}-1)\mathrm{W}(k^{\prime}|k)$
and $T(\xi)=\int_{0}^{\infty}\rmd t\gamma(t|\xi)\alpha(t)$. A further
simplification can be made if we assume the degrees are uncorrelated, in which
case equations (15) reduces to a single equation,

$\displaystyle\hat{g}$
$\displaystyle=1-\sum_{k}\frac{k\mathrm{P}(k)}{{{\langle}k{\rangle}}}\int_{0}^{\infty}\rmd
t\left<\gamma(t|\xi)\right>_{\xi}\left(1-\alpha(t)\hat{g}\right)^{k-1}$ (17)

which allows a non-zero solution when the epidemic threshold exceeds a value
dependent on the average transmissibility,

$\displaystyle\frac{{\langle}k^{2}{\rangle}-{\langle}k{\rangle}}{{{\langle}k{\rangle}}}\geq\frac{1}{\left<T(\xi)\right>_{\xi}}.$
(18)

Another interesting case is where the degrees of nodes at either end of a
randomly selected link are conditionally uncorrelated given their labels $\xi$
and $\xi^{\prime}$ i.e when
$\mathrm{W}(k,\xi;k^{\prime},\xi^{\prime})=\mathrm{W}(\xi,\xi^{\prime})\mathrm{W}(k|\xi)\mathrm{W}(k^{\prime}|\xi^{\prime})$.
For this choice of $\mathrm{W}(k,\xi;k^{\prime},\xi^{\prime})$ we find that
$\mathrm{W}(k^{\prime},\xi^{\prime}|k,\xi)=\frac{\mathrm{W}(k,\xi;k^{\prime},\xi^{\prime})}{\mathrm{W}(k,\xi)}=\frac{\mathrm{W}(\xi;\xi^{\prime})\mathrm{W}(k|\xi)\mathrm{W}(k^{\prime}|\xi^{\prime})}{\mathrm{W}(k|\xi)\mathrm{W}(\xi)}$.
As a result of this, the cavity equations (14) simplify such that
${\hat{g}}_{k,\xi}={\hat{g}}_{\xi}~{}\forall~{}k$ where,

$\displaystyle{\hat{g}}_{\xi}$
$\displaystyle=1-\sum_{k^{\prime},\xi^{\prime}}\mathrm{W}(\xi^{\prime}|\xi)\mathrm{W}(k^{\prime}|\xi^{\prime})\int_{0}^{\infty}\rmd
t\gamma(t|\xi^{\prime})\left(1-\alpha(t)\hat{g}_{\xi^{\prime}}\right)^{k^{\prime}-1}.$
(19)

By using this form of $\mathrm{W}(k,\xi;k^{\prime},\xi^{\prime})$ we have
reduced (14) from a system of $k_{\max}\times M$ equations to just $M$
equations. This is a significant reduction, particularly for networks with
large maximum degree $k_{\max}$. We note that although we have chosen $\xi$ to
label nodes of different transmissibility, we could choose
$\gamma(t|\xi)=\gamma(t)$ and allow $\xi$ to label some other node property
(e.g age). Equations (19) would therefore provide a quick numerical
implementation to find epidemic risk in networks where degree correlations are
generated by a node property, $\xi$, according to
$\mathrm{W}(k,\xi;k^{\prime},\xi^{\prime})=\mathrm{W}(\xi,\xi^{\prime})\mathrm{W}(k|\xi)\mathrm{W}(k^{\prime}|\xi^{\prime})$.
By analysis of the Jacobian of the system of equations (19), we find that an
epidemic will occur if $|\lambda_{1}^{(J_{\xi,\xi^{\prime}})}|>1$ where
$\lambda_{1}^{(J_{\xi,\xi^{\prime}})}$ is the largest eigenvalue of
$J_{\xi,\xi^{\prime}}=\sum_{k^{\prime}}\mathrm{W}(\xi^{\prime}|\xi)(k^{\prime}-1)\mathrm{W}(k^{\prime}|\xi^{\prime})T(\xi^{\prime})$.

###  3.2 Herd immunity for vaccines with partial transmission

Having extended the cavity approach to account for heterogeneity in individual
transmissibility, we can now consider vaccination with partial transmission.
In this case, vaccinated individuals can still catch and transmit the disease,
but we assume their transmissability is reduced, such that an individual $i$
who is vaccinated will have a lower mean infectious time $1/\xi_{i}$ than an
unvaccinated individual. To analyse this scenario, we split a population into
two clusters: an unvaccinated population with low infectious time decay rate
$\xi^{\ell}$ and a vaccinated population with high decay rate
$\xi^{h}>\xi^{\ell}$. The fraction of the population that is vaccinated is
given by $\mathrm{P}(\xi^{h})=1-\mathrm{P}(\xi^{\ell})$. To keep the analysis
simple, we focus on random vaccinations, where one can assume that the labels
$\xi$ are not correlated with the network degrees, and we will assume that
degrees are uncorrelated. The critical fraction $P_{c}(\xi^{h})$ of population
that needs to be vaccinated to prevent epidemics, assuming vaccination with
partial transmission, can be found from (18) by writing

$\displaystyle{\langle}T(\xi){\rangle}_{\xi}$ $\displaystyle=$
$\displaystyle\mathrm{P}(\xi^{h})T(\xi^{h})+[1-\mathrm{P}(\xi^{h})]T(\xi^{\ell})$
(20) $\displaystyle=$ $\displaystyle T(\xi^{\ell})-\mathrm{P}(\xi^{h})\Delta$

with $\Delta=T(\xi^{\ell})-T(\xi^{h})>0$, which, substituted in (18), gives

$\displaystyle\mathrm{P}(\xi^{h})>\frac{T(\xi^{\ell})}{T(\xi^{\ell})-T(\xi^{h})}\left[1-\frac{{{\langle}k{\rangle}}}{T(\xi^{\ell})\left<k(k-1)\right>}\right]\equiv
P_{c}(\xi^{h})$ (21)

Clearly, the fraction of the population that has to be vaccinated to achieve
herd immunity must be greater when the vaccination leads to partial
transmission, as opposed to no transmission. The herd immunity threshold under
vaccination without transmission is given by $1-p_{c}$ from (9), where $T$ is
the transmissibility of the unvaccinated population i.e $T=T(\xi^{\ell})$. One
can express the herd immunity threshold under vaccination with partial
transmission, $P_{c}(\xi^{h})$, in terms of $1-p_{c}$,

$\displaystyle
P_{c}(\xi^{h})=\frac{T(\xi^{\ell})}{T(\xi^{\ell})-T(\xi^{h})}(1-p_{c}).$ (22)

Since $T(\xi^{\ell})>T(\xi^{\ell})-T(\xi^{h})$, it is clear that vaccination
with partial transmission requires a greater number of people to be vaccinated
when compared with vaccination with no transmission i.e
$P_{c}(\xi^{h})>1-p_{c}$.

###  3.3 Networks with correlated structure and transmissibility

In general, social contacts will be correlated with the transmissibility of
individuals. For example, vaccinations are often prioritised for at risk
groups, and for people above a certain age. Furthermore, vaccine distribution
is dependent upon supply chains, and this can lead to higher levels of
vaccination in one geographic area to another. Aside from vaccinations,
transmissibility may correlate with social contacts, for example differences
in transmissibility between children and adults have been noted for some
infectious diseases, and hence if a social interaction network shows
correlation with age, this may affect the spread of epidemics.

To understand how this affects the epidemic risk and vaccination we need to
specify an ensemble of networks where links between individuals are based upon
their respective transmissibility. For simplicity we do not specify any hard
constraints on the degree sequence and instead define an ensemble of networks
with average degree ${{\langle}k{\rangle}}$ where links are drawn in a way
that allows for preferential attachment between nodes with specific
transmissibilities, on the basis of an arbitrary function
$\mathrm{W}(\xi;\xi^{\prime})$ of the labels of the two nodes concerned. More
specifically, we assume that for every pair of nodes $i$ and $j$, links are
drawn randomly and independently with probability

$\displaystyle\mathrm{P}(A_{ij})=\frac{{{\langle}k{\rangle}}}{N}\frac{\mathrm{W}(\xi_{i};\xi_{j})}{\mathrm{P}(\xi_{i})\mathrm{P}(\xi_{j})}\delta_{A_{ij},1}+\left(1-\frac{{{\langle}k{\rangle}}}{N}\frac{\mathrm{W}(\xi_{i};\xi_{j})}{\mathrm{P}(\xi_{i})\mathrm{P}(\xi_{j})}\right)\delta_{A_{ij},0}$
(23)

Here, $\mathrm{W}(\xi;\xi^{\prime})$ can be interpreted as the probability to
draw a link with node labels $\xi$ and $\xi^{\prime}$ at either end, hence it
must be non-negative and normalised to one
$\sum_{\xi,\xi^{\prime}}\mathrm{W}(\xi;\xi^{\prime})=1$. Due to the undirected
nature of the links, we also have
$\mathrm{W}(\xi;\xi^{\prime})=\mathrm{W}(\xi^{\prime};\xi)$. We will refer to
nodes with the same label, $\xi$, as a _cluster_ , hence $\mathrm{W}$ controls
the links between clusters. In C we show that the degree distribution of a
graph drawn from this ensemble is given by

$\displaystyle\mathrm{P}(k)=\sum_{\xi}\mathrm{P}(\xi)\mathrm{P}(k|\xi)=\sum_{\xi}\mathrm{P}(\xi)\rme^{-\bar{k}(\xi)}\frac{(\bar{k}(\xi))^{k}}{k!}$
(24)

where the average degree is
${{\langle}k{\rangle}}=\sum_{\xi}\mathrm{P}(\xi)\bar{k}(\xi)$. This is such
that each cluster of nodes in this ensemble has a Poissonian degree
distribution $\mathrm{P}(k|\xi)$ with mean degree $\bar{k}(\xi)$, where $\xi$
is the cluster label. Furthermore, we find that the degree correlations are of
the special form

$\displaystyle\mathrm{W}(k,\xi;k^{\prime},\xi^{\prime})=\mathrm{W}(\xi;\xi^{\prime})\mathrm{W}(k|\xi)\mathrm{W}(k^{\prime}|\xi^{\prime}),$
(25)

where $\mathrm{W}(k|\xi)=k\mathrm{P}(k|\xi)/\bar{k}(\xi)$, which allows us to
use the reduced form of the cavity equations (19) to study the risk on graphs
from this ensemble. Finally, we note that the marginal of
$\mathrm{W}(\xi;\xi^{\prime})$ is
$\mathrm{W}(\xi)=\sum_{\xi^{\prime}}\mathrm{W}(\xi;\xi^{\prime})=\mathrm{P}(\xi)\bar{k}(\xi)/{{\langle}k{\rangle}}$.

Figure 3: Left: Global risk plotted against infection rate $\beta$ for a
population where nodes have an infectious time distribution
$\gamma(t|\xi)=\xi\rme^{-\xi t}$ and are separated into two sub-types with
$\xi\in\\{\xi^{\ell},\xi^{h}\\}=\\{1.0,5.0\\}$. The network is drawn from the
ensemble (23) with ${{\langle}k{\rangle}}=3$ and the kernel $\mathrm{W}$
defined in (27) with $q=0.15$ and $\bar{k}(\xi^{\ell})={{\langle}k{\rangle}}$
(and $\bar{k}(\xi^{h})={{\langle}k{\rangle}}$, via (26)). Results from cavity
equations (lines) and simulations (symbols) are shown. Simulations were
performed on networks of size $N=2000$, averaged over 500 different sites of
initial infection and 100 repetitions. Dashed lines indicate epidemic
threshold predicted by the Jacobian of (19). Right: For the same population,
we plot the global risk against the fraction of fast recovery nodes
$\mathrm{P}(\xi^{h})$. We set the infection rate to $\beta=1$. Results from
the cavity equations (19) are shown for different inter-connectivity
$q=\\{0.1,0.15,0.2,0.3\\}$ (curves from top to bottom).

In the remainder of this section, we demonstrate results for networks with two
clusters, labelled by $\xi^{h}$ and $\xi^{\ell}<\xi^{h}$, so that the
transmissability kernel $\mathrm{W}(\xi;\xi^{\prime})$ is a $2\times 2$
matrix, $\mathrm{P}(\xi^{h})=1-\mathrm{P}(\xi^{\ell})$ and the average
connectivity is related to the intra-cluster connectivities by

$\displaystyle{{\langle}k{\rangle}}$ $\displaystyle=$
$\displaystyle\bar{k}(\xi^{\ell})+\mathrm{P}(\xi^{h})(\bar{k}(\xi^{h})-\bar{k}(\xi^{\ell}))$
(26) $\displaystyle=$
$\displaystyle\bar{k}(\xi^{h})+\mathrm{P}(\xi^{\ell})(\bar{k}(\xi^{\ell})-\bar{k}(\xi^{h})).$

Given the properties of symmetry, marginalization and normalization of
$\mathrm{W}$, we can parameterise $\mathrm{W}(\xi,\xi^{\prime})$ in terms of
two control parameters, the intra-cluster connectivity $\bar{k}(\xi^{\ell})$
and a parameter, $q$, that controls the inter-cluster connectivity, via the
matrix,

$\mathrm{W}=\begin{pmatrix}\frac{\bar{k}(\xi^{\ell})\mathrm{P}(\xi^{\ell})}{{{\langle}k{\rangle}}}-q&q\\\
q&1-\frac{\bar{k}(\xi^{\ell})\mathrm{P}(\xi^{\ell})}{{{\langle}k{\rangle}}}-q\end{pmatrix}$
(27)

Since the elements of $\mathrm{W}$ are probabilities, they must take values in
$[0,1]$, hence the free parameters can take values in the range
$q{{\langle}k{\rangle}}/\mathrm{P}(\xi^{\ell})<\bar{k}(\xi^{\ell})<(1-q){{\langle}k{\rangle}}/\mathrm{P}(\xi^{\ell})$
with $q\in[0,1/2]$.

In Figure 3 (left panel) the global risk, as predicted from the cavity
equations, is plotted for a population with fast and slow recovery nodes, as a
function of the infection rate $\beta$, for different values of the fast
recovery population density, $\mathrm{P}(\xi^{h})$, and it is found in good
agreement with simulations. Figure 3 (right panel) shows that as the
proportion of nodes of fast recovery increases, herd immunity is eventually
reached. Interestingly, as the inter-connectivity $q$ is increased, the global
risk decreases, as does the herd immunity threshold. This can intuitively be
explained by first considering the scenario where $q=0$ and the network is
split into two disconnected sets of nodes, one where an epidemic is unlikely
to be caused by any node, the other where nodes have a non-zero probability of
causing an epidemic. As $q$ is increased, the number of links between these
sets of nodes increases, with the fast recovery nodes acting as a blockade to
the path of infection, lowering the global risk and herd immunity threshold.

##  4 Social distancing in populations with heterogeneous transmissibility

We now show how social distancing affects epidemic risk in a population with
heterogeneous transmissibility. We model social distancing as a bond
percolation process. We introduce the random binary variables
$\tau_{ij}\in\\{0,1\\}$ which indicates whether a link between two individuals
is broken due to social distancing ($\tau_{ij}=0$) or not ($\tau_{ij}=1$).
These are drawn from the distribution,

$\displaystyle\mathrm{Q}(\boldsymbol{\tau}|\mathbf{y})=\prod_{i<j}\left[(1-y_{ij})\delta_{\tau_{ij},0}+y_{ij}\delta_{\tau_{ij},1}\right]$
(28)

where $1-y_{ij}$ is the probability that a link between $i$ and $j$ is broken.
Under bond percolation the equations for risk take the form,

$\displaystyle r_{i}({\bf A},\boldsymbol{\xi})$
$\displaystyle=1-\int_{0}^{\infty}\rmd
t\,\gamma(t|\xi_{i})\prod_{j\in\partial^{\bf
A}_{i}}\left(1-\tau_{ij}\alpha(t)r_{j}^{(i)}({\bf A},\boldsymbol{\xi})\right)$
(29) $\displaystyle r_{j}^{(i)}({\bf A},\boldsymbol{\xi})$
$\displaystyle=1-\int_{0}^{\infty}\rmd
t\,\gamma(t|\xi_{j})\prod_{\ell\in\partial^{\bf A}_{j}\setminus
i}\left(1-\tau_{jl}\alpha(t)r_{\ell}^{(j)}({\bf A},\boldsymbol{\xi})\right).$
(30)

We consider three different types of bond percolation: random percolation
$y_{ij}=y~{}\forall~{}i,j$, degree-based percolation $y_{ij}=y(k_{i},k_{j})$,
and sub-type percolation $y_{ij}=y(\xi_{i},\xi_{j})$. Following steps in B,
one can find the closed expressions for the global risk

$\displaystyle g$ $\displaystyle=1-\sum_{k\xi}\mathrm{P}(k,\xi)\int\rmd
t\gamma(t|\xi)\left(1-\alpha(t)\tilde{g}_{k,\xi}\right)^{k}$ (31)
$\displaystyle\tilde{g}_{k,\xi}$
$\displaystyle=\sum_{k^{\prime}\xi^{\prime}}y(k,\xi,k^{\prime},\xi^{\prime})\mathrm{W}(k^{\prime},\xi^{\prime}|k,\xi)\left[1-\int\rmd
t\gamma(t|\xi^{\prime})\left(1-\alpha(t)\tilde{g}_{k^{\prime},\xi^{\prime}}\right)^{k^{\prime}-1}\right]$
(32)

where $y(k,\xi,k^{\prime},\xi^{\prime})=y$ for random bond percolation,
$y(k,\xi,k^{\prime},\xi^{\prime})=y(k,k^{\prime})$ for degree-dependent bond
percolation and $y(k,\xi,k^{\prime},\xi^{\prime})=y(\xi,\xi^{\prime})$ for
sub-type dependent bond percolation. In figure 4 we consider an ER network,
with individuals randomly assigned one of two sub-types, high and low
transmissibility, $\xi\in\\{\xi^{\ell},\xi^{h}\\}$ with $\xi^{\ell}<\xi^{h}$,
such that degree and node sub-type are uncorrelated. We then compare the
global risk following random, degree-based, and sub-type bond percolation for
the same fraction of links removed. For degree-based percolation we choose
$y(k,k^{\prime})=\alpha\left(1-\frac{kk^{\prime}}{k_{\max}^{2}}\right)$ to
preferentially remove links attached to nodes of high degree. For sub-type
percolation, we consider the scenario where nodes of low transmissibility
(i.e. fast recovery $\xi=\xi^{h}$) do not reduce their social contact, i.e.
$y(\xi^{h},\xi^{h})=1$, and we set
$y(\xi^{\ell},\xi^{h})=y(\xi^{\ell},\xi^{\ell})$, such that
$y(\xi^{\ell},\xi^{\ell})$ is the only free parameter which controls the
overall fraction of links removed. We see that link percolation based upon
degree and node-type yields lower global risk relative to random link
percolation. Furthermore, for this parameterisation we see that percolation
based upon node-type is preferable to degree-based percolation. This is due to
the fact that degree-based percolation deletes a high concentration of links
between nodes of low transmissibility, as compared to percolation based upon
sub-types.

Figure 4:  Global risk plotted as a function of the fraction $y$ of links
removed from an ER graph with mean degree ${{\langle}k{\rangle}}=5$. Nodes
have infectious time distribution $\gamma(t|\xi)=\xi\rme^{-\xi t}$ and are
assigned to one of two sub-types of slow and fast recovery,
$\xi\in\\{\xi^{\ell},\xi^{h}\\}=\\{1,5\\}$ such that $\xi^{\ell}<\xi^{h}$,
with probability $\mathrm{P}(\xi^{h})=1-\mathrm{P}(\xi^{\ell})=0.65$. The
infection rate is set to $\beta=2$. Results are shown for the cases where
links are removed i) at random ii) with preference for links connecting nodes
of high degree
$y=\sum_{k,k^{\prime}}\mathrm{W}(k,k^{\prime})\alpha\left(1-\frac{kk^{\prime}}{k_{\max}}\right)$,
such that $\alpha$ set the value of $y$, and iii) with preference for nodes
with slow recovery
$y=\sum_{\xi,\xi^{\prime}}\mathrm{W}(\xi;\xi^{\prime})y(\xi,\xi^{\prime})$
where we set $y(\xi^{h},\xi^{h})=1$ and
$y(\xi^{\ell},\xi^{h})=y(\xi^{\ell},\xi^{\ell})$ such that
$y(\xi^{\ell},\xi^{\ell})$ is the free parameter that controls $y$.

##  5 Beyond the mean: distribution of risk in the SIR model

Figure 5: Results from population dynamics, with population size $S=2\times
10^{5}$ and $2\times 10^{6}$ samples to form the histogram, for an Erdös-Rényi
network with mean degree ${{\langle}k{\rangle}}=5$, infection rate $\beta=0.5$
and homogeneous transmission. Left: Cumulative distribution function of the
risk, $\mathrm{P}(X\leq r)$, for inverse mean infectious time $\xi=0.75$.
Results from population dynamics and simulations on a network of size $N=2000$
are shown. Inset shows same information for the risk distribution $\pi(r)$,
with simulations plotted as a histogram and population dynamics plotted as a
solid line. Right: Distribution of the risk distribution $\pi(r)$ (black line)
and degree conditional risk distributions $\pi_{k}(r)$, for
$k\in\\{0,1,2,3,4\\}$ (peaks from left to right) for inverse mean infectious
time $\xi=0.6$.

.

###  5.1 Distribution of risk in networks with degree correlations

In the previous sections we have used the cavity method to obtain a closed set
of equations for the average risk of epidemics in a network where nodes have
heterogeneous infectious time. Previous work has shown that epidemic risk
varies considerably from node to node, even with homogeneous transmission
probabilities, due to variations in the local environment [17]. The presence
of heterogeneities in node infectious times is expected to add an extra source
of fluctuations, therefore studying the distribution of the node risks is
expected to provide important information that is not captured by the mean
risk.

Recently, the cavity method has been used to assess the heterogeneous
behaviour of individual nodes in the context of percolation in sparse networks
[17, 24]. Here we extend this approach to account for degree correlations and
heterogeneities in the node transmissability and calculate the functional form
of the distribution of node risks. The latter is defined as

$\displaystyle\pi(r|{\bf A})=\frac{1}{N}\sum_{i=1}^{N}\delta\left(r-r_{i}({\bf
A})\right)$ (33)

and a closed set of equation can be derived for it by substituting the RHS of
equation (1) into equation (33) and following a similar procedure as in A for
the global risk, which leads, after averaging over the graph ensemble, to

$\displaystyle\pi(r)=\sum_{k\geq
0}\mathrm{P}(k)\left\\{\prod_{j=1}^{k}\int\rmd{\hat{r}}_{j}\hat{\pi}_{k}({\hat{r}}_{j})\right\\}\delta\left(r-\phi(k,\boldsymbol{{\hat{r}}})\right)$
(34) $\displaystyle\hat{\pi}_{k}({\hat{r}})=\sum_{k^{\prime}\geq
1}\mathrm{W}(k^{\prime}|k)\left\\{\prod_{j=1}^{k^{\prime}-1}\int\rmd{\hat{r}}^{\prime}_{j}\hat{\pi}_{k^{\prime}}({\hat{r}}^{\prime}_{j})\right\\}\delta\left({\hat{r}}-\phi(k^{\prime}-1,\boldsymbol{{\hat{r}}}^{\prime})\right)$
(35)

with $\boldsymbol{{\hat{r}}}=\\{{\hat{r}}_{1},...,{\hat{r}}_{k}\\}$,
${\bf{\hat{r}}}^{\prime}=({\hat{r}}^{\prime}_{1},\ldots,{\hat{r}}^{\prime}_{k^{\prime}-1})$
and

$\displaystyle\phi(k,\boldsymbol{\hat{r}})=1-\int_{0}^{\infty}\rmd
t\gamma(t)\prod_{j=1}^{k}\left(1-\alpha(t)\hat{r}_{j}\right).$ (36)

The set of equations (35) for the degree-dependent ‘cavity’ risk distributions
$\hat{\pi}_{k}(\hat{r})$ must be solved first to then solve equation (34) for
the distribution of risk $\pi(r)$. This can be accomplished numerically by
using a population dynamics algorithm [34]. This consists of approximating
each $\hat{\pi}_{k}({\hat{r}})$ by the empirical cavity risk frequencies
computed from a large number (i.e. a population) of cavity risks, which are
updated at each iteration of the algorithm following a stochastic map. We
start by initializing $k_{\max}$ populations of cavity risks, each of size
$S$, by drawing $k_{\max}\times S$ random variables ${\hat{r}}_{i}^{(k)}$,
$i=1,\ldots,S$ in the $[0,1]$ interval. Their empirical distribution
represents the zero-step approximation $\hat{\pi}_{k}^{0}({\hat{r}})$ of
$\hat{\pi}_{k}$. We then evolve the risks such that at each step $n$ the
empirical distribution $\hat{\pi}_{k}^{n}({\hat{r}})$ comes closer and closer
to the invariant distribution $\hat{\pi}_{k}({\hat{r}})$. Each iteration of
the algorithm consists of the following steps, which are repeated until
convergence

  1. 1.

Set $k=1$.

  2. 2.

Draw degree $k^{\prime}$ with probability $\mathrm{W}(k^{\prime}|k)$ .

  3. 3.

Randomly select $k^{\prime}-1$ risks from the $k^{\prime}$-th population of
risks.

  4. 4.

Compute the function $\phi(k^{\prime}-1,\boldsymbol{\hat{r}})$ using the
$k^{\prime}-1$ values of risk.

  5. 5.

Select a random risk from the $k$-th population and set its value to
$\phi(k^{\prime}-1,\boldsymbol{\hat{r}})$.

  6. 6.

Set $k=k+1$ and go back to step (2).

  7. 7.

Repeat all steps until all populations have converged.

Once the populations have converged, i.e. $\hat{\pi}_{k}(\hat{r})$ does not
change upon further iteration, the distribution of risk $\pi(r)$ can be
computed in a similar manner by taking samples from the population and
computing the histogram of $\phi(k,\boldsymbol{\hat{r}})$.

In Figure 5 (left panel) we show results from the population dynamics
algorithm for an ER graph with average connectivity ${{\langle}k{\rangle}}=5$.
Results are in very good agreement with simulations, carried out on networks
of size $N=2000$. As previously shown in [24], the distribution of risk can be
conditioned upon the degree of a node, leading to the degree-conditional risk
distributions

$\displaystyle\pi_{k}(r)=\left\\{\prod_{j=1}^{k}\int\rmd\hat{r}_{j}\hat{\pi}_{k}(\hat{r}_{j})\right\\}\delta\left(r-\phi(k,\boldsymbol{{\hat{r}}})\right).$
(37)

These are plotted in figure 5 (right panel) and show that the peaks and
troughs apparent in $\pi(r)$ are caused by nodes of different degree, with the
peak at $r=0$ arising from disconnected nodes. From this figure we see that
the mean risk increases with degree, showing the benefit of degree-based
vaccination strategies. However, the skew in the degree-conditional risk
distributions means that degrees which have a high _mean_ risk, also contain
nodes of relatively low risk, and vaccinating these nodes would reduce the
overall risk and epidemic size by little, highlighting the limitations of a
degree-based strategy. It can be checked analytically and verified numerically
that the risk distribution for random regular graphs are delta-peaked, with
the peak corresponding to the solution to equation (5) for the global risk,
which suggests that the random nature of the degree of the neighbours of a
node in graphs with heterogeneous structure are the source of variability in
the risk of nodes with the same degree.

Figure 6 (left panel) shows the cumulative distribution of risk for networks
with Poissonian degree distribution and assortative and disassortative degree
correlations, respectively. Results from the cavity method are found in good
agreement with simulations on networks with size $N=2000$. Degree correlations
are seen to have a significant impact on the risk distribution. This is
further shown on the right panel, where the survival function
$\mathrm{P}(X\geq r)$ of assortative graphs is seen to have larger tails than
disassortative or neutral graphs, due to a higher proportion of high risk
nodes.

Figure 6: Results from population dynamics, with population size $S=1\times
10^{4}$ and $1\times 10^{5}$ samples to form the histogram, for a network with
Poissonian degree distribution, mean degree ${{\langle}k{\rangle}}=5$,
infection rate $\beta=0.25$ and homogeneous recovery rate $\xi=0.75$. Left:
cumulative distribution function of the risk, $\mathrm{P}(X\leq r)$, computed
from the cavity method (dashed line) compared with simulations (solid line) on
a network of size $N=2000$, for assortative and disassortative degree
correlations. Right: same information from the cavity method shown for the
survival function, $\mathrm{P}(X\geq r)$, with the curve for a neutral graph
added for comparison

###  5.2 Distribution of risk in networks with heterogeneous transmission

The equations for the distribution of risk can be immediately generalised to
the case of heterogeneous transmission, giving

$\displaystyle\pi(r)$ $\displaystyle=\sum_{k\geq
0,\xi}\mathrm{P}(k,\xi)\left\\{\prod_{j=1}^{k}\int\rmd\hat{r}_{j}\hat{\pi}_{k,\xi}(\hat{r}_{j})\right\\}\delta\left(r-\phi(k,\xi,\boldsymbol{{\hat{r}}})\right)$
(38) $\displaystyle\hat{\pi}_{k,\xi}(\hat{r})$
$\displaystyle=\sum_{k^{\prime}\geq
1,\xi^{\prime}}\mathrm{W}(k^{\prime},\xi^{\prime}|k,\xi)\left\\{\prod_{j=1}^{k^{\prime}-1}\int\rmd\hat{r}^{\prime}_{j}\hat{\pi}_{k^{\prime},\xi^{\prime}}(\hat{r}^{\prime}_{j})\right\\}\delta\left(\hat{r}-\phi(k^{\prime}-1,\xi^{\prime},\boldsymbol{{\hat{r}^{\prime}}})\right)$
(39) $\displaystyle\phi(k,\xi,\boldsymbol{\hat{r}})$
$\displaystyle=1-\int_{0}^{\infty}\rmd
t\gamma(t|\xi)\prod_{j=1}^{k}\left(1-\alpha(t)\hat{r}_{j}\right).$ (40)

These equations can be solved using a generalisation of the population
dynamics algorithm above with $P\times k_{\max}$ populations. The risk
distribution $\pi(r)$ can be de-convoluted either in terms of degree, or node
label $\xi$ or both. To illustrate this, we first consider a random regular
graph with degree ${{\langle}k{\rangle}}=5$, where each node is independently
assigned a random label $\xi\in\\{\xi^{\ell},\xi^{h}\\}$. Figure 7 (left
panel) shows that the risk distribution neatly splits into the distribution of
risk for each cluster. The risk distribution for nodes of a given label $\xi$,
is seen to have several peaks, corresponding to the different neighbourhoods a
node may have: in a regular graph of degree $c$, with two clusters, there are
$c+1$ different labelling configurations for the neighbourhood of a node. The
height of the peaks in this distribution are therefore related to the
multiplicity of these possible labelling configurations. It is important to
note, that these peaks are not delta-peaked, but are broader, and overlap, due
to the risk being dependent not only upon the neighbourhood of a node, but on
the neighbours of neighbours as explored in [35]. Figure 7 (right panel) shows
the risk distribution and the degree conditional risk distributions for the
ensemble defined in (23), with average connectivity ${{\langle}k{\rangle}}=5$.
When compared to figure 5 we see a greater amount of heterogeneity for nodes
of a given degree. This highlights the importance of understanding differences
in transmissability in individuals and its interplay with the network
structure: the efficacy of sub-type based or degree based vaccination depends
on the level of links between each cluster and their relative
transmissability.

Figure 7: Left: Risk distribution $\pi(r)$ of regular graph with connectivity
${{\langle}k{\rangle}}=5$, rate of infection $\beta=0.5$, and two recovery
rates $\xi\in\\{0.5,0.75\\}$ with probability $\mathrm{P}(\xi)=0.5$ for all
$\xi$. The curves with larger peaks correspond to the risk distribution
conditioned upon $\xi$. Right: risk distribution $\pi(r)$ for a network from
the ensemble (23) with mean degree ${{\langle}k{\rangle}}=5$, rate of
infection $\beta=0.5$, and two recovery rates $\xi\in\\{0.6,1.0\\}$ with
probability $\mathrm{P}(\xi=1.0)=0.5$. The transmissability kernel is defined
in (27) with $\bar{k}(\xi=0.6)=\bar{k}(\xi=1.0)={{\langle}k{\rangle}}$. Larger
peaks correspond to the risk distribution conditioned upon $k$. In both plots
the population size is $S=10^{5}$ and $10^{6}$ samples are taken from the
converged cavity distribution to form $\pi(r)$.

###  5.3 Risk distribution in the limit of large connectivity

An exact expression for the risk distribution may be derived in the limit of
large connectivity, for graphs with homogeneous transmissibility, using
methods demonstrated in [24]. To extend these methods to graphs with
heterogeneous transmissibility we use the ansatz that as
${{\langle}k{\rangle}}\to\infty$ the cavity field distributions will become
delta-peaked i.e.
$\hat{\pi}_{k,\xi}(\hat{r})=\delta(\hat{r}-\hat{r}^{*}_{k,\xi})$. Inserting
this into equation (39), multiplying by $\hat{r}$ and integrating both sides
over $\hat{r}$, we find

$\displaystyle\hat{r}^{*}_{k,\xi}=\sum_{k^{\prime},\xi^{\prime}}\mathrm{W}(k^{\prime},\xi^{\prime}|k,\xi)\left[1-\int_{0}^{\infty}\rmd
t\gamma(t|\xi^{\prime})\left(1-\alpha(t)\hat{r}^{*}_{k^{\prime},\xi^{\prime}}\right)^{k^{\prime}-1}\right].$
(41)

Similarly, we insert
$\hat{\pi}_{k,\xi}(\hat{r})=\delta(\hat{r}-\hat{r}^{*}_{k,\xi})$ into equation
(38) which yields,

$\displaystyle\pi(r)$
$\displaystyle=\sum_{k,\xi}\mathrm{P}(k,\xi)\delta\left(r-\left[1-\int_{0}^{\infty}\rmd
t\gamma(t|\xi)\left(1-\alpha(t)\hat{r}^{*}_{k,\xi}\right)^{k}\right]\right).$
(42)

Given that $\pi(r)=\sum_{k,\xi}\pi(r|k,\xi)\mathrm{P}(k,\xi)$ we find the
conditional risk distribution to be delta-peaked under this ansatz,

$\displaystyle\pi(r|k,\xi)$ $\displaystyle=\delta\left(r-f_{\xi}(k)\right)$
(43)

where

$\displaystyle f_{\xi}(k)$ $\displaystyle=1-\int_{0}^{\infty}\rmd
t\gamma(t|\xi)\left(1-\alpha(t)\hat{r}^{*}_{k,\xi}\right)^{k}.$ (44)

We can simplify the above by choosing
$\gamma(t|\xi)=\delta\left(t+\frac{1}{\beta}\ln\left(1-T(\xi)\right)\right)$
such that,

$\displaystyle f_{\xi}(k)$
$\displaystyle=1-\left(1-T(\xi)\hat{r}^{*}_{k,\xi}\right)^{k}.$ (45)

In order to retrieve a non-trivial distribution of risk in the limit
$\bar{k}\to\infty$, we assume that $T(\xi)\ll 1~{}\forall~{}\xi$ (at finite
transmissibility, every node will be near-certain of causing an epidemic in
the large connectivity limit). This allows us to write
$f_{\xi}(k)=1-\rme^{-kT(\xi)\hat{r}^{*}_{k,\xi}}$. To proceed we assume that
degree correlations are of the form
$W(k^{\prime},\xi^{\prime}|k,\xi)=W(\xi^{\prime}|\xi)W(k^{\prime}|\xi^{\prime})$,
which by equation (41) yields $r^{*}_{k,\xi}=r^{*}_{\xi}~{}\forall~{}k$. To
find the risk distribution of nodes of a given label
$\pi_{\xi}(r)=\sum_{k}\pi(r|k,\xi)\mathrm{P}(k|\xi)$, we then use (43) and
properties of the Dirac delta-function to write,

$\displaystyle\pi_{\xi}(r)$
$\displaystyle=\sum_{k}\mathrm{P}(k|\xi)\frac{\delta\left(k-f^{-1}_{\xi}(r)\right)}{|f^{\prime}_{\xi}(f^{-1}_{\xi}(r))|}=\frac{1}{T(\xi)\hat{r}_{\xi}^{*}(1-r)}\sum_{k}\mathrm{P}(k|\xi)\delta\left(k+\frac{\ln(1-r)}{T(\xi)\hat{r}_{\xi}^{*}}\right).$
(46)

To proceed further we must specify the degree distribution, and so we consider
the distribution of risk in graphs drawn from the graph ensemble defined by
(23) and assume that as ${{\langle}k{\rangle}}\to\infty$ the conditional mean
degree $\bar{k}(\xi)\to\infty~{}\forall~{}\xi$. Noting that the degrees in
this ensemble are distributed according to
$\mathrm{P}(k|\xi)=\rme^{-\bar{k}(\xi)}\bar{k}^{k}(\xi)/k!$, which in the
limit of $\bar{k}(\xi)\to\infty$ is well approximated by
$\mathrm{P}(k|\xi)\sim\mathcal{N}(\bar{k}(\xi),\bar{k}(\xi))$, we find,

$\displaystyle\pi_{\xi}(r)=\frac{1}{T(\xi)r^{*}_{\xi}}\sqrt{\frac{1}{2\pi\bar{k}(\xi)}}\exp\\!\left[\\!-\frac{\bar{k}(\xi)}{2}\left(\\!1+\frac{\ln(1-r)}{T(\xi)r^{*}_{\xi}\bar{k}(\xi)}\\!\right)^{2}\\!\\!-\ln(1-r)\\!\right].$
(47)

The unconditional risk distribution is then given by
$\pi(r)=\sum_{\xi}\mathrm{P}(\xi)\pi_{\xi}(r)$. We can compare this expression
with the results of population dynamics for a system with two clusters with
transmissability $T(\xi)=\\{0.0225,0.03\\}$ of equal size
$\mathrm{P}(\xi)=0.5$. We see in figure 8 that for the case where
${{\langle}k{\rangle}}=50$ there is excellent agreement with the expression
(47). In figure 8 we also see the impact of changing the inter-connectivity in
this ensemble: for low $q$ we see two distinct peaks, corresponding to the
risk of nodes of different transmissability, but as $q$ increases the two
peaks begin to overlap, until the distribution becomes unimodal for high $q$.

Figure 8: Conditional risk distribution $\pi_{\xi}(r)$, for a network in the
ensemble (23) with ${{\langle}k{\rangle}}=50$ and two clusters, with recovery
time distribution $\gamma(t|\xi)=\delta(t+\frac{1}{\beta}\ln(1-T(\xi)))$,
parameterised by $T(\xi)\in\\{0.0225,0.03\\}$ with $\mathrm{P}(\xi)=0.5$ and
inter-connectivity $q$. Left: The dashed lines show the solutions to equations
(38) and (39) via population dynamics, with population size $S=10^{5}$ and
$10^{6}$ samples forming the final distribution. The solid curve shows the
large mean limit approximation (47). The inter-connectivity is set to
$q=0.01$. Right: Large mean limit approximation for different values of $q$,
as shown in the legend.

###  5.4 Distributional equations of risk with node and link percolation

Finally, we use the cavity method to investigate the effect of node and link
percolation on the distribution of risk in the SIR model. Node deletion models
perfect vaccination, as before, and link percolation models the loss of a
social contact between two individuals, due to, for example, social-distancing
measures. To do so we introduce a binary variable $\tau_{ij}$ which indicates
whether a link in a contact network has been deleted ($\tau_{ij}=0$) or not
($\tau_{ij}=1$). The local risk subject to node and link percolation is then
written,

$\displaystyle r_{i}({\bf A})$
$\displaystyle=\sigma_{i}\left[1-\int_{0}^{\infty}\rmd
t\,\gamma(t)\prod_{j\in\partial^{\bf
A}_{i}}\left(1-\tau_{ij}\alpha(t)r_{j}^{(i)}({\bf A})\right)\right]$ (48)
$\displaystyle r_{j}^{(i)}({\bf A})$
$\displaystyle=\sigma_{j}\left[1-\int_{0}^{\infty}\rmd
t\,\gamma(t)\prod_{\ell\in\partial^{\bf A}_{j}\setminus
i}\left(1-\tau_{jl}\alpha(t)r_{\ell}^{(j)}({\bf A})\right)\right].$ (49)

We consider the case where node and links are randomly and independently
deleted with probability dependent on their degree, so that

$\displaystyle\mathrm{P}(\boldsymbol{\sigma}|{\bf k})$
$\displaystyle=\prod_{i}\left[(1-x(k_{i}))\delta_{\sigma_{i},0}+x(k_{i})\delta_{\sigma_{i},1}\right]$
(50) $\displaystyle\mathrm{Q}(\boldsymbol{\tau}|{\bf k})$
$\displaystyle=\prod_{i<j}\left[(1-y(k_{i},k_{j}))\delta_{\tau_{ij},0}+y(k_{i},k_{j})\delta_{\tau_{ij},1}\right].$
(51)

By proceeding as in A, and averaging over the graph ensemble, we find that the
distribution of risk is given by,

$\displaystyle\pi(r)$ $\displaystyle=$
$\displaystyle\sum_{k}\mathrm{P}(k)\sum_{\sigma}P(\sigma|k)$
$\displaystyle\times\left\\{\prod_{j=1}^{k}\sum_{k_{j}^{\prime}}\mathrm{W}(k_{j}^{\prime}|k)\sum_{\hat{\tau}_{j}}\mathrm{Q}(\hat{\tau}_{j}|k,k_{j}^{\prime})\int\rmd{\hat{r}}_{j}\hat{\pi}_{k_{j}^{\prime}}({\hat{r}}_{j})\right\\}$
$\displaystyle\times\delta\left(r-\phi(k,\sigma,\boldsymbol{\hat{\tau}},\boldsymbol{{\hat{r}}})\right)$
$\displaystyle\hat{\pi}_{k}({\hat{r}})$ $\displaystyle=$
$\displaystyle\sum_{\sigma}\mathrm{P}(\sigma|k)\left\\{\prod_{\ell=1}^{k-1}\sum_{k_{\ell}}\mathrm{W}(k_{\ell}^{\prime}|k)\sum_{\hat{\tau}^{\prime}_{\ell}}\mathrm{Q}(\hat{\tau}^{\prime}_{\ell}|k,k_{\ell}^{\prime})\int\rmd{\hat{r}}_{\ell}^{\prime}\hat{\pi}_{k_{\ell}^{\prime}}({\hat{r}}_{\ell}^{\prime})\right\\}$
(52)
$\displaystyle\times\delta\left({\hat{r}}-\phi(k-1,\sigma,\boldsymbol{\hat{\tau}^{\prime}},\boldsymbol{{\hat{r}}^{\prime}})\right)$

with $\boldsymbol{\hat{\tau}}=\\{\hat{\tau}_{1},...,\hat{\tau}_{k}\\}$,
${\boldsymbol{\hat{\tau}}^{\prime}}=(\hat{\tau}^{\prime}_{1},\ldots,\hat{\tau}^{\prime}_{k^{\prime}-1})$
and

$\displaystyle\phi(k,\sigma,\boldsymbol{\hat{\tau}},\boldsymbol{\hat{r}})=\sigma\left[1-\int_{0}^{\infty}\rmd
t\gamma(t)\prod_{j=1}^{k}\left(1-\alpha(t)\hat{\tau}_{j}\hat{r}_{j}\right)\right].$
(53)

Equation (52) can be solved by a generalisation of the populations dynamics
algorithm above. To evaluate the cavity distribution $\hat{\pi}({\hat{r}})$:

  1. 1.

Set $k=1$.

  2. 2.

Draw $\sigma$ with probability $\mathrm{P}(\sigma|k)$

  3. 3.

Draw $k^{\prime}$ degree from the distribution $\mathrm{W}(k^{\prime}|k)$ .

  4. 4.

Select a risk from the $k^{\prime}$-th population of risks.

  5. 5.

Select $\hat{\tau}^{\prime}$ with probability
$\mathrm{Q}(\hat{\tau}^{\prime}|k,k^{\prime})$.

  6. 6.

Repeat steps (3-5) $k-1$ times, and store each $\tau$, and individual to form
the vectors $\boldsymbol{\hat{\tau}^{\prime}}$ and
$\boldsymbol{\hat{r}^{\prime}}$.

  7. 7.

Select a random risk from the $k$-th population and set its value to
$\phi(k-1,\sigma,\boldsymbol{\hat{\tau}^{\prime}},\boldsymbol{\hat{r}^{\prime}})$.

  8. 8.

Set $k=k+1$ and go back to step (2).

  9. 9.

Repeat all steps until all populations have converged.

Once the populations estimating $\hat{\pi}_{k}({\hat{r}})$ have converged, one
can compute the risk distribution in a similar manner. We note that, by
equations (52), the risk distribution $\pi(r)$ can be deconvoluted in terms of
$k$ and $\sigma$. By conditioning upon $\sigma=1$ it is possible to retrieve
the risk distribution of unvaccinated nodes, therefore providing a method to
see how vaccination strategies affect risk in the network of unvaccinated
nodes.

In figure 9 we show the effects of node and link deletion, representing
perfect vaccination and social distancing respectively. We consider two cases:
random and degree-based strategies. For random vaccination we have
$x(k)=x~{}\forall~{}k$ and $y(k,k^{\prime})=y~{}\forall~{}k,k^{\prime}$. For
degree based strategies we set $x(k)=1-\alpha\frac{k}{k_{\max}}$ and
$y(k,k^{\prime})=1-\alpha\frac{kk^{\prime}}{k_{\max}^{2}}$ such that nodes of
high degree and links connecting nodes of high degree are more likely to be
deleted. We show the case where $15\%$ of nodes/links are deleted in an ER
graph. We see that for the same fraction of nodes/links deleted, the overall
risk of the population is lower under random node deletion than random link
deletion. Furthermore, risk under degree-based node deletion is lower than
degree-based link deletion. This is consistent with the idea that social
distancing measures, which assume an imperfect break in the chain of social
contacts, are a less effective method of epidemic prevention than vaccination.
It is however, interesting to compare the effect of degree-based link
deletion, to random node deletion: in this case link deletion has the greater
reduction of risk in the population. This highlights the importance of
strategy in vaccination, although deleting a node blocks all paths of
infection through that node, a targeted strategy of link deletion can lead to
greater reductions in risk than random node deletion.

Figure 9: Survival function for the risk of nodes in a Erdös-Rényi graph with
mean degree ${{\langle}k{\rangle}}=3$. Infection rate is $\beta=0.7$ and nodes
have exponential infectious time distribution $\gamma(t)=\xi\rme^{-\xi t}$
with mean infectious time $\xi=0.75$. The dotted line indicates case where no
node or link deletion has occurred. In all other cases $15\%$ nodes/links are
deleted, subject to either random (blue curves) or degree-based (red curves)
deletion where high degree nodes/links are preferentially selected for
deletion. For degree based node deletion we have nodes deleted according to
(50) with $x(k)=1-\alpha\frac{k}{k_{\max}}$, and for degree based link
deletion we delete nodes according to (51) with
$y(k,k^{\prime})=1-\alpha\frac{kk^{\prime}}{k_{\max}^{2}}$

##  6 Discussion and conclusion

Predicting the effect of vaccination in epidemic models on networks continues
to be a source of interesting inquiry. In our work we have relaxed the common
assumption that vaccination blocks all transmission and instead assumed that
vaccinated individuals have reduced transmissibility. From a modelling
perspective this is useful when assessing the risk of epidemics and there is a
lack of data for the impact of a given vaccine on transmission. To explore
this, we have extended the cavity method to account for heterogeneity in the
transmissibility of groups of individuals. These groups could represent
differences in age, pre-existing medical conditions, or any other known
correlates with transmissibility. Here however we have focused on a population
of vaccinated and unvaccinated individuals, distinguished from each other by
their transmissibility. Our results reveal that herd immunity is still
attained from vaccination with partial transmission, if vaccination reduces
transmissibility to a sufficiently low level that depends on both the
infectious rate of the disease and contact network topology. We have shown
that vaccination with partial transmission requires a greater proportion of
the population to be vaccinated to achieve herd immunity, as intuitively
expected. This highlights the importance of estimating the transmissibility of
an infectious disease, and the impact that vaccination has upon individual
transmissibility. However, it also highlights the robustness of vaccination as
a strategy for the mitigation of epidemics, as even a partial reduction of
transmission can achieve the desired herd immunity effect.

The benefit of the cavity approach to the SIR model on networks is that it
allows one to evaluate the epidemic risk without simulation. Our equations for
the global risk provide a quick method to calculate epidemic risk in
populations with groups of different transmissibility. By solving the
distributional equations of epidemic risk we have revealed the stark impact of
heterogeneity in transmissibility. For graphs with homogeneous degree
structure, the distribution of risk is delta-peaked when transmissibility of
each individual is the same, but non-trivial when nodes take one of two values
of transmissibility. This reveals that it is a strong assumption to make that
the transmissibility of individuals is homogeneous, and that to ignore such
heterogeneity can potentially dampen epidemic mitigation strategies that focus
on network topology alone. Indeed our equations show that by introducing
groups of different transmissibility, the epidemic risk will depend upon the
size and transmissibility of each group, as well as the assortativity of links
between nodes in different groups.

While vaccination is an important part of epidemic mitigation, we have also
shown that the cavity method can model social distancing via bond percolation.
Our framework allows us to compare random bond percolation with percolation
that targets either links connecting nodes of high degree, or nodes of a
particular sub-type. This may be useful when studying variations in the social
distancing of different demographics and their impact on epidemic risk. One
example could be studying the effect of relaxed social distancing amongst the
vaccinated population who are perhaps less risk averse.

When studying populations with groups of different transmissibility we have
provided a method to reduce the complexity of the cavity equations for the
global risk by suitable choice of the graph ensemble. While we have used this
graph ensemble to study the effect of heterogeneity in transmissibility, we
can also choose that the transmissibility of each group is the same, such that
the groups of nodes now only vary in their mean degree and assortativity
between different groups. Therefore, this ensemble provides a low complexity
method to study the effect of degree correlations on epidemic risk, where
correlations are driven by some hidden variable, perhaps reflecting
differences in demographic.

The distributional equations of risk give insight into the impact of different
vaccination strategies and social distancing behaviours beyond their effect on
the global risk. This is particularly interesting with respect to the tails of
these distributions, corresponding to the probability of high risk nodes. We
have shown that assortative graphs have fat tailed distributions, and a
further avenue of work may be to elucidate more information about the tail of
the risk distribution in graphs with strong degree correlations from the
cavity method. Interestingly, the distributional equations allow one to study
the risk distribution in the sub-network of unvaccinated nodes with no further
technicalities.

There are several limitations of our work. Firstly, the equations we provide
only describe the steady state of the SIR model, but dynamics are an important
consideration in the modelling of epidemics, for example the timing of peak
infections is of particular interest. An interesting avenue of further work
could be to extend our investigations to the dynamics of the SIR model using
the dynamical cavity method, exploring how differences in transmission and
assortativity affect the peak of infections. Additionally, our modelling
approach is restricted to static networks, whereas realistically contact
networks vary with time, as discussed in [36, 37]. Lastly, it is important to
note that the SIR model does not allow for reinfection of individuals, making
it unsuitable for the modelling of some infectious diseases where reinfection
is a common occurrence. Despite these limitations our work shows that the
cavity approach to the epidemic risk in the SIR model reveals rich behaviour
governed by the heterogeneity in the transmissibility of individuals.

All authors would like to thank Reimer Kühn for useful discussions.

## References

  * [1] Salathé M, Kazandjieva M, Lee J W, Levis P, Feldman M W and Jones J H A high-resolution human contact network for infectious disease transmission 2010 Proc Natl Acad Sci U S A 107 22020–22025
  * [2] Dall’Asta L Inhomogeneous percolation models for spreading phenomena in random graphs 2005 J Stat Mech 2005 P08011
  * [3] Newman M E Spread of epidemic disease on networks 2002 Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics 66
  * [4] Gleeson J P Cascades on correlated and modular random networks 2008 Phys Rev E Stat Nonlin Soft Matter Phys 77 046117
  * [5] Karrer B and Newman M E Message passing approach for general epidemic models 2010 Phys Rev E 82
  * [6] Vázquez A and Moreno Y Resilience to damage of graphs with degree correlations 2003 Phys Rev E 67 015101
  * [7] Shiraki Y and Kabashima Y Cavity analysis on the robustness of random networks against targeted attacks: Influences of degree-degree correlations 2010 Phys Rev E 82 036101
  * [8] Bianconi G, Sun H, Rapisardi G and Arenas A Message-passing approach to epidemic tracing and mitigation with apps 2021 Phys Rev Res 3 L012014
  * [9] Min B and Castellano C Message-passing theory for cooperative epidemics 2020 Chaos 30 023131
  * [10] Sun H, Saad D and Lokhov A Y Competition, Collaboration, and Optimization in Multiple Interacting Spreading Processes 2021 Phys Rev X 11
  * [11] Moore S and Rogers T Heterogeneous node responses to multi-type epidemics on networks 2020 Proc R Soc Lond A Math Phys Sci 476
  * [12] Pastor-Satorras R and Vespignani A Immunization of complex networks 2002 Phys Rev E 65 036104
  * [13] Miller J C and Hyman J M Effective vaccination strategies for realistic social networks 2007 Physica A 386 780–785
  * [14] Ma J, Van Den Driessche P and Willeboordse F H The importance of contact network topology for the success of vaccination strategies 2013 J Theor Biol 325 12–21
  * [15] Hamilton K E and Pryadko L P Tight Lower Bound for Percolation Threshold on an Infinite Graph 2014 Phys Rev Lett 113 208701
  * [16] Karrer B, Newman M and Zdeborová L Percolation on Sparse Networks 2014 Phys Rev Lett 113 208702
  * [17] Rogers T Assessing node risk and vulnerability in epidemics on networks 2015 Europhys Lett 109
  * [18] Torres L, Chan K S, Tong H and Eliassi-Rad T Nonbacktracking eigenvalues under node removal: X-centrality and targeted immunization 2021 SIAM J Math Data Sci 3 656–675
  * [19] Mossong J, Hens N, Jit M, Beutels P, Auranen K, Mikolajczyk R, Massari M, Salmaso S, Tomba G S, Wallinga J et al. Social contacts and mixing patterns relevant to the spread of infectious diseases 2008 PLoS Med 5 e74
  * [20] Mistry D et al. Inferring high-resolution human mixing patterns for disease modeling 2021 Nat Commun 12
  * [21] Anderson R M and May R M Spatial, Temporal, and Genetic Heterogeneity in Host Populations And the Design of Immunization Programmes 1984 Math Med Biol 1 233–266
  * [22] Gou W and Jin Z How heterogeneous susceptibility and recovery rates affect the spread of epidemics on networks 2017 Infect Dis Model 2 353–367
  * [23] Baxter G and Timár G Degree dependent transmission probabilities in epidemic processes 2021 arXiv preprint arXiv:2102.11602
  * [24] Kühn R and Rogers T Heterogeneous micro-structure of percolation in sparse networks 2017 Europhys Lett 118 68003
  * [25] Melnik S, Hackett A, Porter M A, Mucha P J and Gleeson J P The unreasonable effectiveness of tree-based theory for networks with clustering 2011 Phys Rev E Stat Nonlin Soft Matter Phys 83 036112
  * [26] Mézard M, Parisi G and Virasoro M A 1987 Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications vol 9 (World Scientific Publishing Company)
  * [27] Newman M 2018 Networks (Oxford university press)
  * [28] Gillespie D T Exact stochastic simulation of coupled chemical reactions 1977 J Phys Chem 81 2340–2361
  * [29] Kiss I Z, Miller J C, Simon P L et al. Mathematics of epidemics on networks 2017 Cham: Springer 598
  * [30] Miller J C and Ting T EoN (Epidemics on Networks): a fast, flexible Python package for simulation, analytic approximation, and analysis of epidemics on networks 2019 J Open Source Softw 4 1731
  * [31] Annibale A, Coolen A C, Fernandes L P, Fraternali F and Kleinjung J Tailored graph ensembles as proxies or null models for real networks I: Tools for quantifying structure 2009 J Phys A Math Theor 42
  * [32] Newman M E J Assortative Mixing in Networks 2002 Phys Rev Lett 89 208701
  * [33] Kiss I Z, Röst G and Vizi Z Generalization of Pairwise Models to non-Markovian Epidemics on Networks 2015 Phys Rev Lett 115 078701
  * [34] Mézard M and Parisi G The bethe lattice spin glass revisited 2001 Eur Phys J B 20 217–233
  * [35] Tishby I, Biham O, Katzav E and Kühn R Revealing the microstructure of the giant component in random graph ensembles 2018 Physical Review E 97 042318
  * [36] Onaga T, Gleeson J P and Masuda N Concurrency-induced transitions in epidemic dynamics on temporal networks 2017 Phys Rev Lett 119 108301
  * [37] P Peixoto T and Gauvin L Change points, memory and epidemic spreading in temporal networks 2018 Sci Rep 8 1–10

##  Appendix A Risk with homogeneous transmission and node and link deletion

Here we show how to derive a closed set of equations for the average, or
‘global’, risk that a node causes an epidemic across a network that is subject
to node and link deletion, when nodes are homogeneous in transmissibility. We
introduce a binary random variable $\sigma_{i}\in\\{0,1\\}~{}\forall~{}i$ to
indicate if a node $i$ is deleted ($\sigma_{i}$=0) or not ($\sigma_{i}=1$) and
$\tau_{ij}\in\\{0,1\\}~{}\forall~{}i,j$ to indicate if the link between $i$
and $j$ is deleted ($\tau_{ij}=0$) or not ($\tau_{ij}=1$). We consider the
case where node and links are randomly and independently deleted with
probability dependent on their degree, according to,

$\displaystyle\mathrm{P}(\sigma_{i}|k_{i})=\left(1-x(k_{i})\right)\delta_{\sigma_{i},0}+x(k_{i})\delta_{\sigma_{i},1}$
(54)
$\displaystyle\mathrm{Q}(\tau_{ij}|k_{i},k_{j})=\left(1-y(k_{i},k_{j})\right)\delta_{\tau_{ij},0}+y(k_{i},k_{j})\delta_{\tau_{ij},1}.$
(55)

For a graph with adjacency matrix ${\bf A}$ the risk that node $i$ causes an
epidemic is given by,

$\displaystyle r_{i}({\bf A})$ $\displaystyle=$
$\displaystyle\sigma_{i}\left[1-\int_{0}^{\infty}\rmd
t\,\gamma(t)\prod_{j\in\partial^{\bf
A}_{i}}(1-\alpha(t)\tau_{ij}r_{j}^{(i)}({\bf A}))\right]$ (56) $\displaystyle
r_{j}^{(i)}({\bf A})$ $\displaystyle=$
$\displaystyle\sigma_{j}\left[1-\int_{0}^{\infty}\rmd
t\,\gamma(t)\prod_{\ell\in\partial^{\bf A}_{j}\setminus
i}\left(1-\alpha(t)\tau_{j\ell}r_{\ell}^{(j)}({\bf A})\right)\right].$ (57)

Summing (56) over $i$, dividing by $N$ and inserting unity in the form

$1=\sum_{\sigma}\delta_{\sigma,\sigma_{i}}\sum_{k\geq
0}\delta_{k,|\partial_{i}^{\bf A}|}\prod_{j\in\partial_{i}^{\bf
A}}\int\rmd\hat{r}_{j}\,\delta(\hat{r}_{j}-r_{j}^{(i)}({\bf
A}))\sum_{\hat{\tau}_{j}}\delta_{\hat{\tau}_{j},\tau_{ij}}$ (58)

we obtain

$\displaystyle g({\bf A})$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{\sigma}\sum_{i=1}^{N}\delta_{\sigma,\sigma_{i}}\sigma\Bigg{\\{}1-\sum_{k\geq
0}\delta_{k,k_{i}({\bf A})}\left[\prod_{j\in\partial_{i}^{\bf
A}}\int\rmd\hat{r}_{j}\,\delta(\hat{r}_{j}-r_{j}^{(i)}({\bf
A}))\sum_{\hat{\tau}_{j}}\delta_{\hat{\tau}_{j},\tau_{ij}}\right]$ (59)
$\displaystyle\times\int_{0}^{\infty}\rmd
t\gamma(t)\prod_{j\in\partial_{i}^{\bf
A}}(1-\alpha(t)\hat{\tau}_{j}\hat{r}_{j})\Bigg{\\}}$ $\displaystyle=$
$\displaystyle\sum_{\sigma}\sigma\Bigg{[}\mathrm{P}(\sigma)-\sum_{k\geq
0}\sum_{{\boldsymbol{\hat{\tau}}}}\int\rmd{{\boldsymbol{\hat{r}}}}\,\mathrm{W}_{c}(k,\sigma;{\boldsymbol{\hat{\tau}}},{{\boldsymbol{\hat{r}}}}|{\bf
A})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times\int_{0}^{\infty}\rmd
t\gamma(t)\prod_{j\in\partial_{i}^{\bf
A}}\left(1-\alpha(t)\hat{\tau}_{j}\hat{r}_{j}\right)\Bigg{]}$

where we have defined $k_{i}({\bf A})=|\partial^{\bf A}_{i}|$ the degree of
node $i$ in network ${\bf A}$,
${{\boldsymbol{\hat{r}}}}=(\hat{r}_{1},\ldots,\hat{r}_{k})$ the cavity fields
of the neighbours of a site with degree $k$,
${\boldsymbol{\hat{\tau}}}=(\hat{\tau}_{1},\ldots,\hat{\tau}_{k})$ the labels
that describe which links of a site with degree $k$ have been deleted,
$\mathrm{P}(\sigma)=N^{-1}\sum_{i}\delta_{\sigma,\sigma_{i}}$ the probability
that a node has label $\sigma$, and

$\displaystyle\mathrm{W}_{c}(k,\sigma;{\boldsymbol{\hat{\tau}}},{{\boldsymbol{\hat{r}}}}|{\bf
A})=\frac{1}{N}\sum_{i=1}^{N}\delta_{k,k_{i}({\bf
A})}\delta_{\sigma,\sigma_{i}}\prod_{j\in\partial_{i}^{\bf
A}}\delta(\hat{r}_{j}-r_{j}^{(i)}({\bf A}))\delta_{\hat{\tau}_{j},\tau_{ij}}$
(60)

the likelihood that a site drawn at random in network ${\bf A}$ has degree
$k$, label $\sigma$, and neighbours with cavity fields,
${{\boldsymbol{\hat{r}}}}$, with the links to these neighbours deleted
according to ${\boldsymbol{\hat{\tau}}}$. We denote probability distributions
of variables that are linked by an edge with $\mathrm{W}(~{};~{})$. We have
also denoted distributions which depend upon the cavity fields
${{\boldsymbol{\hat{r}}}}$ with a subscript $c$ and refer to them as ‘cavity
distributions’. We note that integrating a cavity distribution
$\mathrm{W}_{c}$ over its cavity fields leads to a normal (i.e. not cavity)
distribution $\mathrm{W}$. Furthermore, we note that integrating (60) over
${\boldsymbol{\hat{\tau}}}$ and ${{\boldsymbol{\hat{r}}}}$ leads to a
distribution of single site variables, that we denote with $P$. For later
reference, we will denote with $Q$ distributions of single link quantities, in
accordance with (55).

Using Bayes relation, we then have
$\mathrm{W}_{c}(k,\sigma;{\boldsymbol{\hat{\tau}}},{{\boldsymbol{\hat{r}}}}|{\bf
A})=\mathrm{P}(k,\sigma|{\bf
A})\mathrm{W}_{c}({\boldsymbol{\hat{\tau}}},{{\boldsymbol{\hat{r}}}}|k,{\bf
A})$ where we have used that, in the absence of correlations between node
labels, and between node labels and link labels, the distribution
$\mathrm{W}_{c}({\boldsymbol{\hat{\tau}}},{{\boldsymbol{\hat{r}}}}|k,{\bf A})$
in the cavity graph where node $i$ has been removed is independent of the
label $\sigma$ of node $i$. We proceed by writing
$\mathrm{W}_{c}({\boldsymbol{\hat{\tau}}},{{\boldsymbol{\hat{r}}}}|k,{\bf
A})=\sum_{\mathbf{q}}\mathrm{W}_{c}(\boldsymbol{q},{\boldsymbol{\hat{\tau}}},{{\boldsymbol{\hat{r}}}}|k,{\bf
A})$ where $\mathbf{q}=(q_{1},...,q_{k})$ are the degrees of the neighbours of
a node with degree $k$. Using Bayes relations, we can write
$\mathrm{W}_{c}(\boldsymbol{q},{\boldsymbol{\hat{\tau}}},{{\boldsymbol{\hat{r}}}}|k,{\bf
A})=\mathrm{W}(\mathbf{q}|k,{\bf
A})\mathrm{Q}({\boldsymbol{\hat{\tau}}}|k,\mathbf{q})\mathrm{W}_{c}({{\boldsymbol{\hat{r}}}}|{\bf
A},k,\mathbf{q})$ where we have used the conditional independence of
${{\boldsymbol{\hat{r}}}}$ (when conditioned on the degrees $k$ and ${{\bf
q}}$), on ${\boldsymbol{\hat{\tau}}}$, in the cavity graph where node $i$ and
all its links have been removed. Here
$\mathrm{Q}({\boldsymbol{\hat{\tau}}}|k,\mathbf{q})$ is the joint distribution
of deleting $k$ links according to ${\boldsymbol{\hat{\tau}}}$. Here
$\mathrm{W}_{c}({{\boldsymbol{\hat{r}}}}|{\bf A},k,\mathbf{q})$ is the joint
distribution of the cavity fields ${{\boldsymbol{\hat{r}}}}$ given that the
node removed from the cavity graph has degree $k$ with neighbours with degrees
${{\bf q}}$. By equation (55) we have that
$\mathrm{Q}({\boldsymbol{\hat{\tau}}}|k,\mathbf{q})=\prod_{j=1}^{k}\mathrm{Q}(\hat{\tau}_{j}|k,q_{j})$.
Furthermore, in the limit $N\to\infty$, by virtue of a locally tree-like
assumption, $\mathrm{W}_{c}({{\boldsymbol{\hat{r}}}}|{\bf A},k,\mathbf{q})$
also factorises $\mathrm{W}_{c}({{\boldsymbol{\hat{r}}}}|{\bf
A},k,\mathbf{q})=\prod_{j=1}^{k}\mathrm{W}_{c}(\hat{r}_{j}|{\bf A},k,q_{j})$,
for any $k>0$. For $k=0$, we can set
$\mathrm{W}_{c}({{\boldsymbol{\hat{r}}}}|{\bf
A},0,\mathbf{q})=\delta({{\boldsymbol{\hat{r}}}})$ as the product in (56) is
empty and evaluates to one, regardless of ${{\boldsymbol{\hat{r}}}}$. When we
insert this into (59) we find,

$\displaystyle g({\bf A})\\!=\\!\sum_{\sigma,k\geq
0,\mathbf{q}}\\!\sigma\mathrm{P}(k,\sigma|{\bf A})\mathrm{W}(\mathbf{q}|k,{\bf
A})\left[1\\!-\\!\int_{0}^{\infty}\rmd
t\,\gamma(t)\prod_{j=1}^{k}\left(1-\alpha(t)y(k,q_{j})\hat{g}_{k,q_{j}}({\bf
A})\right)\right]$ (61)

where $\hat{g}_{k,q}({\bf A})=\int\rmd\hat{r}\mathrm{W}_{c}(\hat{r}|{\bf
A},k,q)\,\hat{r}$ is the average contribution to the cavity field of a random
node with degree $k$ from a random neighbour with degree $q$. We now assume
that the degrees of the neighbours of a random site with degree $k$ are
independent when conditioned upon $k$ such that $\mathrm{W}({\bf q}|k,{\bf
A})=\prod_{j=1}^{k}\mathrm{W}(q_{j}|k,{\bf A})$. While this assumption is only
true for particular graph ensembles, it will act as a useful approximation for
ensembles where it is not true, as it leads to a closed set of equations.
Inserting this into (61) yields,

$\displaystyle g({\bf A})=\sum_{k\geq 0}\mathrm{P}(k|{\bf
A})x(k)\left[1-\int_{0}^{\infty}\rmd
t\,\gamma(t)\prod_{j=1}^{k}\left(1-\alpha(t)\tilde{g}_{k}({\bf
A})\right)\right]$ (62)

where

$\displaystyle\tilde{g}_{k}({\bf A})=\sum_{q}\mathrm{W}(q|k,{\bf
A})y(k,q){\hat{g}}_{k,q}.$ (63)

To make progress, we need an equation for $\tilde{g}_{k}({\bf A})$,

$\displaystyle\tilde{g}_{k}({\bf A})$
$\displaystyle=\sum_{q}\mathrm{W}(q|k,{\bf
A})y(k,q)\int\rmd\hat{r}~{}\mathrm{W}_{c}(\hat{r}|k,q,{\bf A})\hat{r}$ (64)
$\displaystyle=\sum_{q}\mathrm{W}(q|k,{\bf
A})y(k,q)\int\rmd{{\boldsymbol{\hat{r}}}}\frac{\mathrm{W}_{c}(k;{{\boldsymbol{\hat{r}}}},q|{\bf
A})}{\mathrm{P}(k|{\bf A})\mathrm{W}(q|k,{\bf
A})}\,\frac{1}{k}\sum_{j=1}^{k}\hat{r}_{j}$ (65)
$\displaystyle=\sum_{q}\mathrm{W}(q|k,{\bf A})y(k,q)\sum_{{{\bf
q}}}\int\rmd{{\boldsymbol{\hat{r}}}}\frac{\mathrm{W}_{c}(k;{{\boldsymbol{\hat{r}}}},{{\bf
q}}|{\bf A})}{\mathrm{P}(k|{\bf A})\mathrm{W}(q|k,{\bf
A})}\,\frac{1}{k}\sum_{j=1}^{k}\hat{r}_{j}\delta_{q,q_{j}}$
$\displaystyle=\sum_{q}\frac{y(k,q)}{k\mathrm{P}(k|{\bf
A})}\sum_{j=1}^{k}\sum_{{\bf
q}}\int\rmd{{\boldsymbol{\hat{r}}}}\frac{1}{N}\sum_{i}^{N}\delta_{k,k_{i}({\bf
A})}[\prod_{l\in\partial_{i}^{{\bf A}}}\delta_{q_{l},k_{l}({\bf
A})}\delta(\hat{r}_{l}-r_{l}^{(i)}({\bf A}))]\hat{r}_{j}\delta_{q,q_{j}}$
$\displaystyle=\sum_{q}\frac{y(k,q)}{k\mathrm{P}(k|{\bf
A})}\sum_{j=1}^{k}\int\rmd{{\boldsymbol{\hat{r}}}}\frac{1}{N}\sum_{i}^{N}\delta_{k,k_{i}({\bf
A})}[\prod_{l\in\partial_{i}^{{\bf A}}}\delta(\hat{r}_{l}-r_{l}^{(i)}({\bf
A}))]\hat{r}_{j}\delta_{q,k_{j}}$
$\displaystyle=\sum_{q}\frac{y(k,q)}{k\mathrm{P}(k|{\bf
A})}\frac{1}{N}\sum_{j=1}^{k}\sum_{i}^{N}A_{ij}\delta_{k,k_{i}({\bf
A})}\delta_{q,k_{j}({\bf A})}\hat{r}_{j}^{(i)}({\bf A})$ (67)

where we have used Bayes relations and the definition

$\mathrm{W}_{c}(k;{{\bf q}},{{\boldsymbol{\hat{r}}}}|{\bf
A})=\frac{1}{N}\sum_{i}^{N}\delta_{k,k_{i}({\bf
A})}[\prod_{j\in\partial_{i}^{{\bf A}}}\delta_{q_{j},k_{j}({\bf
A})}\delta(\hat{r}_{j}-r_{j}^{(i)}({\bf A}))]$

Hence, using the cavity equation (57) we find an expression for
$\tilde{g}_{k}$ as follows,

$\displaystyle\tilde{g}_{k}({\bf A})$ $\displaystyle\\!=\\!$
$\displaystyle\sum_{q}\frac{y(k,q)}{Nk\mathrm{P}(k|{\bf
A})}\sum_{i,j=1}^{N}A_{ij}\delta_{k,k_{i}({\bf A})}\delta_{q,k_{j}({\bf
A})}\sigma_{j}\Bigg{[}1\\!-\\!\int_{0}^{\infty}\\!\rmd
t\gamma(t)\\!\\!\prod_{\ell\in\partial^{\bf A}_{j}\setminus
i}\\!\\!\left(1\\!-\\!\alpha(t)\tau_{j\ell}r_{\ell}^{(j)}\right)\Bigg{]}$ (68)
$\displaystyle\\!=\\!$ $\displaystyle\sum_{q}\frac{y(k,q)}{Nk\mathrm{P}(k|{\bf
A})}\Bigg{[}\prod_{\ell\in\partial_{j}^{\bf A}\setminus
i}\int\rmd\hat{r}^{\prime}_{\ell}\delta(\hat{r}^{\prime}_{\ell}-r_{\ell}^{(j)})\sum_{\hat{\tau}^{\prime}_{\ell}}\delta_{\hat{\tau}^{\prime}_{\ell},\tau_{j\ell}}\Bigg{]}$
$\displaystyle\times\sum_{\sigma}\sum_{i,j=1}^{N}A_{ij}\delta_{k,k_{i}({\bf
A})}\delta_{q,k_{j}({\bf
A})}\delta_{\sigma,\sigma_{j}}\sigma\left[1-\int_{0}^{\infty}\rmd
t\,\gamma(t)\prod_{\ell=1}^{q-1}\left(1-\alpha(t)\hat{\tau}^{\prime}_{\ell}\hat{r}^{\prime}_{\ell}\right)\right]$
$\displaystyle\\!=\\!$ $\displaystyle\sum_{q}\frac{y(k,q)\bar{k}({\bf
A})}{k\mathrm{P}(k|{\bf
A})}\int\rmd{{\boldsymbol{\hat{r}}}}^{\prime}\sum_{{\boldsymbol{\hat{\tau}}}^{\prime}}\sum_{\sigma}\mathrm{W}_{c}(k;q,\sigma;{\boldsymbol{\hat{\tau}}}^{\prime},{{\boldsymbol{\hat{r}}}}^{\prime}|{\bf
A})\sigma$ $\displaystyle\times\left[1-\int_{0}^{\infty}\rmd
t\,\gamma(t)\prod_{\ell=1}^{q-1}\left(1-\alpha(t)\hat{\tau}^{\prime}_{\ell}\hat{r}^{\prime}_{\ell}\right)\right]$

where

$\displaystyle\mathrm{W}_{c}(k;q,\sigma;{\boldsymbol{\hat{\tau}}}^{\prime},{{\boldsymbol{\hat{r}}}}^{\prime}|{\bf
A})=\frac{\sum_{ij}A_{ij}\delta_{k,k_{i}({\bf
A})}\delta_{\sigma,\sigma_{j}}\,\delta_{q,k_{j}({\bf
A})}[\prod_{\ell\in\partial_{j}^{\bf A}\setminus
i}\delta(\hat{r}^{\prime}_{\ell}-r_{\ell}^{(j)}({\bf
A}))\delta_{\hat{\tau}_{\ell}^{\prime},\tau_{j\ell}}]}{N\bar{k}({\bf A})}$
(69)

is the likelihood that a randomly drawn link connects a node $i$ with degree
$k$ to a node $j$ with degree $q$, label $\sigma$ and neighbours ($i$
excluded) with cavity fields
${{\boldsymbol{\hat{r}}}}^{\prime}=(\hat{r}^{\prime}_{1},\ldots,\hat{r}^{\prime}_{q-1})$
and links deleted according to
${\boldsymbol{\hat{\tau}}}^{\prime}=(\hat{\tau}^{\prime}_{1},\ldots,\hat{\tau}^{\prime}_{q-1})$.
We next use Bayes and the independence of ${\boldsymbol{\hat{\tau}}}^{\prime}$
and ${{\boldsymbol{\hat{r}}}}^{\prime}$ on $\sigma$ when conditioned on $q$ to
write
$\mathrm{W}_{c}(k;q,\sigma;{\boldsymbol{\hat{\tau}}}^{\prime},{{\boldsymbol{\hat{r}}}}^{\prime}|{\bf
A})=\mathrm{W}(k;q,\sigma|{\bf
A})\mathrm{W}_{c}({\boldsymbol{\hat{\tau}}}^{\prime},{{\boldsymbol{\hat{r}}}}^{\prime}|k,q,{\bf
A})$ and write
$\mathrm{W}_{c}({\boldsymbol{\hat{\tau}}}^{\prime},{{\boldsymbol{\hat{r}}}}^{\prime}|k,q,{\bf
A})=\sum_{{{\bf
q}}^{\prime}}\mathrm{W}_{c}({\boldsymbol{\hat{\tau}}}^{\prime},{{\boldsymbol{\hat{r}}}}^{\prime},{{\bf
q}}^{\prime}|k,q,{\bf A})$ where ${{\bf
q}}^{\prime}=(q^{\prime}_{1},\dots,q^{\prime}_{q-1})$ are the degrees of the
neighbours of a node with degree $q$, excluding the neighbour with degree $k$.
Using again the independence of ${{\boldsymbol{\hat{r}}}}^{\prime}$ and
${\boldsymbol{\hat{\tau}}}^{\prime}$ when conditioned on the degrees $q$ and
${{\bf q}}^{\prime}$, we have
$\mathrm{W}_{c}({\boldsymbol{\hat{\tau}}}^{\prime},{{\boldsymbol{\hat{r}}}}^{\prime},{{\bf
q}}^{\prime}|k,q,{\bf A})=\mathrm{W}({{\bf q}}^{\prime}|q,{\bf
A})\mathrm{Q}({\boldsymbol{\hat{\tau}}}^{\prime}|q,{{\bf
q}}^{\prime})\mathrm{W}_{c}({{\boldsymbol{\hat{r}}}}^{\prime}|q^{\prime},{{\bf
q}}^{\prime},{\bf A})$ where we also note that ${{\bf q}}^{\prime}$,
${{\boldsymbol{\hat{r}}}}^{\prime}$ and ${\boldsymbol{\hat{\tau}}}^{\prime}$
are independent of $k$ when conditioned on $q$. By equation (55) the
${\boldsymbol{\hat{\tau}}}^{\prime}$ factorise, in addition to the
${{\boldsymbol{\hat{r}}}}^{\prime}$ which factorise under the tree-like
assumption, and we again assume that ${{\bf q}}^{\prime}$ factorise when
conditioned on $q$, hence we have

$\mathrm{W}_{c}({{\bf q}}^{\prime}|q,{\bf
A})\mathrm{Q}({\boldsymbol{\hat{\tau}}}^{\prime}|q,{{\bf
q}}^{\prime})\mathrm{W}_{c}({{\boldsymbol{\hat{r}}}}^{\prime}|q,{{\bf
q}}^{\prime},{\bf A})=\prod_{\ell=1}^{q-1}\mathrm{W}(q^{\prime}_{\ell}|q,{\bf
A})\mathrm{Q}(\hat{\tau}^{\prime}_{\ell}|q,q^{\prime}_{\ell})\mathrm{W}_{c}(\hat{r}^{\prime}_{\ell}|q,q^{\prime}_{\ell},{\bf
A}).$ (70)

Inserting this into (68) leads to the system of equations,

$\displaystyle\tilde{g}_{k}({\bf A})$ $\displaystyle=$
$\displaystyle\sum_{q\geq
1}\sum_{\sigma}y(k,q)\frac{\mathrm{W}(k;q,\sigma|{\bf A})}{\mathrm{W}(k|{\bf
A})}\sigma$ (71) $\displaystyle\times\left[1-\int_{0}^{\infty}\rmd
t\,\gamma(t)\left(1-\alpha(t)\sum_{q^{\prime}}\mathrm{W}(q^{\prime}|q,{\bf
A})y(q,q^{\prime})\hat{g}_{q,q^{\prime}}({\bf A})\right)^{q-1}\right]$

where $\mathrm{W}(k|{\bf A})=\sum_{q}\mathrm{W}(k;q|{\bf
A})=k\mathrm{P}(k|{\bf A})/\bar{k}({\bf A})$ and the sum over $q$ has been
restricted to $q\geq 1$, as $\mathrm{W}(k;q|{\bf A})$ vanishes for $q=0$.
Finally, by using (63) we have the following set of equations

$\displaystyle g({\bf A})$ $\displaystyle=$ $\displaystyle\sum_{k\geq
0}\mathrm{P}(k|{\bf A})x(k)\left[1\\!-\\!\int_{0}^{\infty}\\!\\!\\!\\!\rmd
t\,\gamma(t)\\!\left(1-\alpha(t)\tilde{g}_{k}({\bf A})\right)^{k}\right]$ (72)
$\displaystyle\tilde{g}_{k}({\bf A})$ $\displaystyle=$
$\displaystyle\sum_{q\geq 1}y(k,q)\mathrm{W}(q|k,{\bf
A})x(q)\left[1\\!-\\!\int_{0}^{\infty}\rmd
t\,\gamma(t)\left(1-\alpha(t)\tilde{g}_{q}({\bf A})\right)^{q-1}\right]$ (73)

where we have used (54) and denoted $\mathrm{W}(q|k,{\bf
A})=\mathrm{W}(k;q|{\bf A})/\mathrm{W}(k|{\bf A})$. These show that the global
risk $g({\bf A})$ has the same self-averaging properties of $\mathrm{P}(k|{\bf
A})$ and $\mathrm{W}(k;q|{\bf A})$. If the set of equations (73) is solved
numerically, the solution can be substituted into equation (72) to find the
global risk.

##  Appendix B Risk with heterogeneous transmission and link deletion

In this section we find closed equations for the global risk when nodes have
heterogeneous transmissibility. For brevity we do not consider node deletion
in this section. We consider link deletion but now based upon the label and
degree of nodes, such that,

$\displaystyle\mathrm{Q}(\tau_{ij})=\left(1-y(k_{i},\xi_{i},k_{j},\xi_{j})\right)\delta_{\tau_{ij},0}+y(k_{i},\xi_{i},k_{j},\xi_{j})\delta_{\tau_{ij},1}$
(74)

where $1-y(k,\xi,k^{\prime},\xi^{\prime})$ is the probability of deleting a
link with a node with degree $k$ and label $\xi$ at one end and a node with
degree $k^{\prime}$ and label $\xi^{\prime}$ at the other. Summing (29) over
$i$, dividing by $N$ and inserting unity in the form

$1=\sum_{\sigma}\delta_{\sigma,\sigma_{i}}\sum_{k\geq
0}\delta_{k,|\partial_{i}^{\bf A}|}\prod_{j\in\partial_{i}^{\bf
A}}\int\rmd\hat{r}_{j}\,\delta(\hat{r}_{j}-r_{j}^{(i)}({\bf
A}))\sum_{\hat{\tau}_{j}}\delta_{\hat{\tau}_{j},\tau_{ij}}$ (75)

we obtain

$\displaystyle g({\bf A},\boldsymbol{\xi})$ $\displaystyle=$ $\displaystyle
1-\frac{1}{N}\sum_{i=1}^{N}\sum_{k\geq 0}\delta_{k,k_{i}({\bf
A})}\sum_{\xi}\delta_{\xi,\xi_{i}}\left[\prod_{j\in\partial_{i}^{\bf
A}}\int\rmd\hat{r}_{j}\,\delta(\hat{r}_{j}-r_{j}^{(i)}({\bf
A},\boldsymbol{\xi}))\sum_{\hat{\tau}_{j}}\delta_{\hat{\tau}_{j},\tau_{ij}}\right]$
$\displaystyle\times\int_{0}^{\infty}\rmd
t\gamma(t|\xi)\prod_{j\in\partial_{i}^{\bf
A}}(1-\alpha(t)\hat{\tau}_{j}\hat{r}_{j})$ $\displaystyle=$ $\displaystyle
1-\sum_{k\geq
0,\xi,{\boldsymbol{\hat{\tau}}}}\int\rmd{{\boldsymbol{\hat{r}}}}\,\mathrm{W}_{c}(k,\xi;{{\boldsymbol{\hat{r}}}},{\boldsymbol{\hat{\tau}}}|{\bf
A},\boldsymbol{\xi})\int_{0}^{\infty}\rmd
t\,\gamma(t|\xi)\prod_{j\in\partial_{i}^{\bf
A}}\left(1-\alpha(t)\hat{\tau}_{j}\hat{r}_{j}\right)$

where we have defined
${{\boldsymbol{\hat{r}}}}=(\hat{r}_{1},\ldots,\hat{r}_{k})$ and

$\displaystyle\mathrm{W}_{c}(k,\xi;{{\boldsymbol{\hat{r}}}},{\boldsymbol{\hat{\tau}}}|{\bf
A},\boldsymbol{\xi})=\frac{1}{N}\sum_{i=1}^{N}\delta_{k,k_{i}({\bf
A})}\delta_{\xi,\xi_{i}}\prod_{j\in\partial_{i}^{\bf
A}}\delta(\hat{r}_{j}-r_{j}^{(i)}({\bf
A},\boldsymbol{\xi}))\delta_{\hat{\tau}_{j},\tau_{ij}}$ (77)

the likelihood that a site drawn at random in network ${\bf A}$ has degree
$k$, label $\xi$ and neighbours with cavity fields
${{\boldsymbol{\hat{r}}}}=(\hat{r}_{1},\ldots,\hat{r}_{k})$, and links deleted
according to
${\boldsymbol{\hat{\tau}}}=(\hat{\tau}_{1},\ldots,\hat{\tau}_{k})$. Similarly
to before, we use Bayes relation to write
$\mathrm{W}_{c}(k,\xi;{{\boldsymbol{\hat{r}}}},{\boldsymbol{\hat{\tau}}}|{\bf
A},\boldsymbol{\xi})=\mathrm{P}(k,\xi|{\bf
A},\boldsymbol{\xi})\mathrm{W}_{c}({{\boldsymbol{\hat{r}}}},{\boldsymbol{\hat{\tau}}}|k,\xi,{\bf
A},\boldsymbol{\xi})$ and we then write
$\mathrm{W}_{c}({{\boldsymbol{\hat{r}}}},{\boldsymbol{\hat{\tau}}}|k,\xi,{\bf
A},\boldsymbol{\xi})=\sum_{{\boldsymbol{\zeta}},{{\bf
q}}}\mathrm{W}_{c}({{\boldsymbol{\hat{r}}}},{\boldsymbol{\hat{\tau}}},{\boldsymbol{\zeta}},{{\bf
q}}|k,\xi,{\bf A},\boldsymbol{\xi})$, where
${\boldsymbol{\zeta}}=(\zeta_{1},\ldots,\zeta_{k})$ and ${{\bf
q}}=(q_{1},\ldots,q_{k})$ are the labels and degrees of the neighbours of a
random site with degree $k$ and label $\xi$. Using Bayes theorem this becomes
$\mathrm{W}_{c}({{\boldsymbol{\hat{r}}}},{\boldsymbol{\hat{\tau}}},{\boldsymbol{\zeta}},{{\bf
q}}|k,\xi,{\bf A},\boldsymbol{\xi})=\mathrm{W}({\boldsymbol{\zeta}},{{\bf
q}}|k,\xi,{\bf
A},\boldsymbol{\xi})\mathrm{Q}({\boldsymbol{\hat{\tau}}}|k,\xi,{\boldsymbol{\zeta}},{{\bf
q}})\mathrm{W}_{c}({{\boldsymbol{\hat{r}}}}|k,\xi,{\boldsymbol{\zeta}},{{\bf
q}},{\bf A},\boldsymbol{\xi})$ where we have also used the independence of
${\boldsymbol{\hat{\tau}}}$ and ${{\boldsymbol{\hat{r}}}}$ when conditioned on
${\boldsymbol{\zeta}}$ and ${{\bf q}}$. For $k=0$ we can set
$\mathrm{W}_{c}({{\boldsymbol{\hat{r}}}}|0,\xi,{\boldsymbol{\zeta}},{{\bf
q}},{\bf A},\boldsymbol{\xi})=\delta({{\boldsymbol{\hat{r}}}})$ as the
products in (LABEL:eq:hetero-risk) evaluate to one. In the limit $N\to\infty$,
by virtue of a locally tree-like assumption, the joint distribution of the
cavity fields ${{\boldsymbol{\hat{r}}}}$ factorises
$\mathrm{W}_{c}({{\boldsymbol{\hat{r}}}}|k,\xi,{\boldsymbol{\zeta}},{{\bf
q}},{\bf
A},\boldsymbol{\xi})=\prod_{j=1}^{k}\mathrm{W}_{c}(\hat{r}_{j}|k,\xi,q_{j},\zeta_{j},{\bf
A},\boldsymbol{\xi})$, as does
$\mathrm{Q}({\boldsymbol{\hat{\tau}}}|k,\xi,{\boldsymbol{\zeta}},{{\bf
q}})=\prod_{j=1}^{k}\mathrm{Q}(\hat{\tau}_{j}|k,\xi,q_{j},\zeta_{j})$ by
equation (74). Again, we assume the degrees and labels of the neighbours of a
random site with degree $k$ and label $\xi$, are independent, when conditioned
on the degree $k$ and label $\xi$, such that
$\mathrm{W}({\boldsymbol{\zeta}},{{\bf q}}|k,\xi,{\bf
A},\boldsymbol{\xi})=\prod_{j=1}^{k}\mathrm{W}(\zeta_{j},q_{j}|k,\xi,{\bf
A},\boldsymbol{\xi})$. As noted above, although this is only true for
particular graph ensembles, it acts as a useful approximation for ensembles
where it is not true, as it leads to a closed set of equations. When this is
inserted into (LABEL:eq:hetero-risk) we find,

$\displaystyle g({\bf A},\boldsymbol{\xi})=1-\sum_{k,\xi}\mathrm{P}(k,\xi|{\bf
A},\boldsymbol{\xi})\int_{0}^{\infty}\rmd
t\,\gamma(t|\xi)\left(1-\alpha(t)\tilde{g}_{k,\xi}({\bf
A},\boldsymbol{\xi})\right)^{k}$ (78)

where $\tilde{g}_{k,\xi}({\bf
A},\boldsymbol{\xi})=\sum_{q,\xi^{\prime}}\mathrm{W}(q,\zeta|k,\xi,{\bf
A},\boldsymbol{\xi})y(k,\xi,q,\zeta)\int\rmd\hat{r}\mathrm{W}_{c}(\hat{r}^{\prime}|k,\xi,q,\zeta,{\bf
A},\boldsymbol{\xi})\,\hat{r}$ is the average cavity field of a random
neighbour of a random node with degree $k$ and label $\xi$. To make progress,
we need an equation for $\tilde{g}_{k,\xi}({\bf A},\boldsymbol{\xi})$. To this
end we can write $\tilde{g}_{k,\xi}({\bf A},\boldsymbol{\xi})$ as

$\displaystyle\hat{g}_{k,\xi}({\bf A},\boldsymbol{\xi})$ $\displaystyle=$
$\displaystyle\sum_{q,\zeta}\mathrm{W}(q,\zeta|k,\xi,{\bf
A},\boldsymbol{\xi})y(k,\xi,q,\zeta)\int\rmd\hat{r}\mathrm{W}_{c}(\hat{r}|k,\xi,q,\zeta,{\bf
A},\boldsymbol{\xi})\,\hat{r}$ (79) $\displaystyle=$
$\displaystyle\sum_{q,\zeta}y(k,\xi,q,\zeta)\int\rmd{{\boldsymbol{\hat{r}}}}\,\frac{\mathrm{W}_{c}(k,\xi;q,\zeta,{{\boldsymbol{\hat{r}}}}|{\bf
A},\boldsymbol{\xi})}{\mathrm{P}(k,\xi|{\bf
A},\boldsymbol{\xi})}\,\frac{1}{k}\sum_{j=1}^{k}\hat{r}_{j}$ $\displaystyle=$
$\displaystyle\sum_{q,\zeta}y(k,\xi,q,\zeta)\sum_{{{\bf
q}},\boldsymbol{\zeta}}\int\rmd{{\boldsymbol{\hat{r}}}}\,\frac{\mathrm{W}_{c}(k,\xi;{{\bf
q}},\boldsymbol{\zeta},{{\boldsymbol{\hat{r}}}}|{\bf
A},\boldsymbol{\xi})}{\mathrm{P}(k,\xi|{\bf
A},\boldsymbol{\xi})}\,\frac{1}{k}\sum_{j=1}^{k}\hat{r}_{j}\delta_{q,q_{j}}\delta_{\zeta,\zeta_{j}}$
$\displaystyle=$
$\displaystyle\sum_{q,\zeta}\frac{y(k,\xi,q,\zeta)}{k\mathrm{P}(k,\xi|{\bf
A},\boldsymbol{\xi})}\sum_{j=1}^{k}\sum_{{{\bf
q}},\boldsymbol{\zeta}}\int\rmd{{\boldsymbol{\hat{r}}}}\,\frac{1}{N}\sum_{i}\delta_{k,k_{i}({\bf
A})}\delta_{\xi,\xi_{i}}$
$\displaystyle\times\left[\prod_{l\in\partial_{i}^{{\bf
A}}}\delta(\hat{r}_{l}-r_{l}^{(i)}({\bf
A},\boldsymbol{\xi}))\delta_{q_{l},k_{l}({\bf
A})}\delta_{\zeta_{l},\xi_{l}}\right]\,\hat{r}_{j}\delta_{q,q_{j}}\delta_{\zeta,\zeta_{j}}$
$\displaystyle=$
$\displaystyle\sum_{q,\zeta}\frac{y(k,\xi,q,\zeta)}{Nk\mathrm{P}(k,\xi|{\bf
A},\boldsymbol{\xi})}\sum_{i,j=1}^{N}A_{ij}\delta_{k,k_{i}({\bf
A})}\delta_{\xi,\xi_{i}}\delta_{q,k_{j}({\bf
A})}\delta_{\zeta,\xi_{j}}r_{j}^{(i)}({\bf A},\boldsymbol{\xi})$

Hence, using the cavity equation (30) we find an expression for
$\tilde{g}_{k,\xi}$ as follows,

$\displaystyle\tilde{g}_{k,\xi}({\bf A},\boldsymbol{\xi})$ $\displaystyle=$
$\displaystyle\sum_{q,\zeta}\frac{y(k,\xi,q,\zeta)}{Nk\mathrm{P}(k,\xi|{\bf
A},\boldsymbol{\xi})}\sum_{i,j=1}^{N}A_{ij}\delta_{k,k_{i}({\bf
A})}\delta_{\xi,\xi_{i}}\delta_{q,k_{j}({\bf
A})}\delta_{\zeta,\xi_{j}}r_{j}^{(i)}({\bf A},\boldsymbol{\xi})$
$\displaystyle=$
$\displaystyle\sum_{q,\zeta}\frac{y(k,\xi,q,\zeta)}{Nk\mathrm{P}(k,\xi|{\bf
A},\boldsymbol{\xi})}\sum_{i,j=1}^{N}A_{ij}\delta_{k,k_{i}({\bf
A})}\delta_{\xi,\xi_{i}}\delta_{q,k_{j}({\bf A})}\delta_{\zeta,\xi_{j}}$
$\displaystyle\times\left\\{\prod_{\ell\in\partial_{j}^{\bf A}\setminus
i}\int\rmd\hat{r}^{\prime}_{\ell}\,\delta(\hat{r}^{\prime}_{\ell}-r_{\ell}^{(j)}({\bf
A},\boldsymbol{\xi}))\sum_{\hat{\tau}^{\prime}_{\ell}}\delta_{\hat{\tau}^{\prime}_{\ell},\tau_{j\ell}}\right\\}$
$\displaystyle\times\left[1-\int_{0}^{\infty}\rmd
t\,\gamma(t|\zeta)\prod_{\ell\in\partial^{\bf A}_{j}\setminus
i}\left(1-\alpha(t)\hat{\tau}^{\prime}_{\ell}\hat{r}^{\prime}_{\ell}\right)\right]$
$\displaystyle=$
$\displaystyle\sum_{q,\zeta}\frac{y(k,\xi,q,\zeta)}{k\mathrm{P}(k,\xi|{\bf
A},\boldsymbol{\xi})}\sum_{{\boldsymbol{\hat{\tau}}}^{\prime}}\int\rmd{{\boldsymbol{\hat{r}}}}^{\prime}\mathrm{W}_{c}(k,\xi;q,\zeta;{\boldsymbol{\hat{\tau}}}^{\prime},{{\boldsymbol{\hat{r}}}}^{\prime}|{\bf
A},\boldsymbol{\xi})$ $\displaystyle\times\left[1-\int_{0}^{\infty}\rmd
t\,\gamma(t|\zeta)\prod_{\ell=1}^{q-1}\left(1-\alpha(t)\hat{\tau}^{\prime}_{\ell}\hat{r}^{\prime}_{\ell}\right)\right]$

where

$\displaystyle\mathrm{W}_{c}(k,\xi;q,\zeta;{\boldsymbol{\hat{\tau}}}^{\prime},{{\boldsymbol{\hat{r}}}}^{\prime}|{\bf
A},\boldsymbol{\xi})=\frac{\sum_{ij}A_{ij}\delta_{k,k_{i}({\bf
A})}\,\delta_{q,k_{j}({\bf
A})}\delta_{\xi,\xi_{i}}\delta_{\zeta,\xi_{j}}\prod_{\ell\in\partial_{j}^{\bf
A}\setminus i}\delta(\hat{r}^{\prime}_{\ell}-r_{\ell}^{(j)}({\bf
A},\boldsymbol{\xi}))\delta_{\hat{\tau}^{\prime}_{\ell},\tau_{j\ell}}}{N\bar{k}({\bf
A})}$

is the likelihood that a randomly drawn link connects a node $i$ with degree
$k$ and label $\xi$ to a node $j$ with degree $q$, label $\zeta$ and
neighbours (excluded $i$) with cavity fields
${{\boldsymbol{\hat{r}}}}^{\prime}=(\hat{r}^{\prime}_{1},\ldots,\hat{r}^{\prime}_{q-1})$
with links deleted according to
${\boldsymbol{\hat{\tau}}}^{\prime}=(\hat{\tau}^{\prime}_{1},\ldots,\hat{\tau}^{\prime}_{q-1})$.
By Bayes theorem
$\mathrm{W}_{c}(k,\xi;q,\zeta;{\boldsymbol{\hat{\tau}}}^{\prime},{{\boldsymbol{\hat{r}}}}^{\prime}|{\bf
A},\boldsymbol{\xi})=\mathrm{W}(k,\xi;q,\zeta|{\bf
A},\boldsymbol{\xi})\mathrm{W}_{c}({\boldsymbol{\hat{\tau}}}^{\prime},{{\boldsymbol{\hat{r}}}}^{\prime}|k,\xi,q,\zeta,{\bf
A},\boldsymbol{\xi})$ and we my then write
$\mathrm{W}_{c}({\boldsymbol{\hat{\tau}}}^{\prime},{{\boldsymbol{\hat{r}}}}^{\prime}|k,\xi,q,\zeta,{\bf
A},\boldsymbol{\xi})=\sum_{{{\bf
q}}^{\prime},{\boldsymbol{\zeta}}^{\prime}}\mathrm{W}_{c}({{\bf
q}}^{\prime},{\boldsymbol{\zeta}}^{\prime},{\boldsymbol{\hat{\tau}}}^{\prime},{{\boldsymbol{\hat{r}}}}^{\prime}|k,\xi,q,\zeta,{\bf
A},\boldsymbol{\xi})$ where ${{\bf
q}}^{\prime}=(q^{\prime}_{1},\ldots,q^{\prime}_{q-1})$ and
${\boldsymbol{\zeta}}^{\prime}=(\zeta^{\prime}_{1},\ldots,\zeta^{\prime}_{q-1})$
are the degrees and labels of the neighbours of a site with degree $q$ and
label $\zeta$ which is itself connected to a site with degree $k$ and label
$\xi$. By Bayes theorem we may write, $\mathrm{W}_{c}({{\bf
q}}^{\prime},{\boldsymbol{\zeta}}^{\prime},{\boldsymbol{\hat{\tau}}}^{\prime},{{\boldsymbol{\hat{r}}}}^{\prime}|k,\xi,q,\zeta,{\bf
A},\boldsymbol{\xi})=\mathrm{W}({{\bf
q}}^{\prime},{\boldsymbol{\zeta}}^{\prime}|q,\xi,{\bf
A},\boldsymbol{\xi})\mathrm{Q}({\boldsymbol{\hat{\tau}}}^{\prime},{{\boldsymbol{\hat{r}}}}^{\prime}|q,\xi,{{\bf
q}}^{\prime},{\boldsymbol{\zeta}}^{\prime})\mathrm{W}_{c}({{\boldsymbol{\hat{r}}}}^{\prime}|q,\xi,{{\bf
q}}^{\prime},{\boldsymbol{\zeta}}^{\prime},{\bf A},\boldsymbol{\xi})$, and
noted that ${\boldsymbol{\hat{\tau}}}^{\prime}$ and
${{\boldsymbol{\hat{r}}}}^{\prime}$ are independent when conditioned on the
labels and degrees. We have also noted that ${{\bf q}}^{\prime}$,
${\boldsymbol{\zeta}}^{\prime}$, ${\boldsymbol{\hat{\tau}}}^{\prime}$ and
${{\boldsymbol{\hat{r}}}}^{\prime}$ are independent of $k$ and $\xi$. By
equation (74) we have
$\mathrm{Q}({\boldsymbol{\hat{\tau}}}^{\prime}|q,\xi,{{\bf
q}}^{\prime},{\boldsymbol{\zeta}}^{\prime})=\prod_{l=1}^{q-1}\mathrm{Q}(\hat{\tau}^{\prime}_{\ell}|q,\xi,q^{\prime}_{\ell},\zeta^{\prime}_{\ell})$
and by the tree-like assumption
$\mathrm{W}_{c}({{\boldsymbol{\hat{r}}}}^{\prime}|q,\xi,{{\bf
q}}^{\prime},{\boldsymbol{\zeta}}^{\prime},{\bf
A},\boldsymbol{\xi})=\prod_{\ell=1}^{q-1}\mathrm{W}_{c}(\hat{r}^{\prime}_{\ell}|q,\xi,q^{\prime}_{\ell},\zeta^{\prime}_{\ell},{\bf
A},\boldsymbol{\xi})$. Furthermore, we again assume that degrees and labels of
neighbours factorise when conditioned upon the degree and labels of their
neighbours i.e. $\mathrm{W}({{\bf
q}}^{\prime},{\boldsymbol{\zeta}}^{\prime}|q,\xi,{\bf
A},\boldsymbol{\xi})=\prod_{\ell=1}^{q-1}\mathrm{W}(q^{\prime}_{\ell},\zeta^{\prime}_{\ell}|q,\xi,{\bf
A},\boldsymbol{\xi})$. Inserting this into (LABEL:app:hetero-g-tau-near) leads
to a closed set of equations for $\tilde{g}_{k,\xi}({\bf
A},\boldsymbol{\xi})$,

$\displaystyle\tilde{g}_{k,\xi}({\bf A},\boldsymbol{\xi})$ $\displaystyle=$
$\displaystyle\sum_{q,\zeta}y(k,\xi,q,\zeta)\mathrm{W}(q,\zeta|k,\xi,{\bf
A},\boldsymbol{\xi})$ (82) $\displaystyle\times\left[1-\int_{0}^{\infty}\rmd
t\,\gamma(t|\zeta)\left(1-\alpha(t)\,\tilde{g}_{q,\zeta}({\bf
A},\boldsymbol{\xi})\right)^{q-1}\right]$

Assuming that $\mathrm{P}(k,\xi|{\bf A},\boldsymbol{\xi})$ and
$\mathrm{W}(q,\zeta|k,\xi,{\bf A},\boldsymbol{\xi})$ are self-averaging over
the graph ensemble and the node label distribution, we finally average
equations (78) and (82) over the distribution $\mathrm{P}({\bf
A},\boldsymbol{\xi})$ of graphs and node labels, obtaining

$\displaystyle g$
$\displaystyle=1-\sum_{k,\xi}\mathrm{P}(k,\xi)\int_{0}^{\infty}\rmd
t\,\gamma(t|\xi)\left(1-\alpha(t)\,\tilde{g}_{k,\xi}\right)^{k}$ (83)
$\displaystyle\tilde{g}_{k,\xi}$
$\displaystyle=\sum_{q,\zeta}y(k,\xi,q,\zeta)\mathrm{W}(q,\zeta|k,\xi)\left[1-\int_{0}^{\infty}\rmd
t\,\gamma(t|\zeta)\left(1-\alpha(t)\,\tilde{g}_{q,\zeta}\right)^{q-1}\right]$
(84)

where $\mathrm{P}(k,\xi)={\langle}\mathrm{P}(k,\xi|{\bf
A},\boldsymbol{\xi}){\rangle}_{{\bf A},\boldsymbol{\xi}}$ and
$\mathrm{W}(q,\zeta|k,\xi)={\langle}\mathrm{W}(q,\zeta|k,\xi,{\bf
A},\boldsymbol{\xi}){\rangle}_{{\bf A},\boldsymbol{\xi}}$ with
${\langle}\cdot{\rangle}_{{\bf A},\boldsymbol{\xi}}=\sum_{{\bf
A},\boldsymbol{\xi}}\cdot\mathrm{P}({\bf A},\boldsymbol{\xi})$.

##  Appendix C An ensemble of networks linking nodes of similar or dissimilar
transmissability

In this section we derive relations for the average degree distribution
$\mathrm{P}(k|\boldsymbol{\xi})$ and average degree correlations
$\mathrm{W}(k,\xi;k^{\prime},\xi^{\prime}|\boldsymbol{\xi})$ for networks
drawn from the ensemble

$\displaystyle\mathrm{P}(\boldsymbol{A}|\boldsymbol{\xi})$
$\displaystyle=\prod_{i<j}\left[\frac{{\langle}k{\rangle}}{N}\frac{\mathrm{W}(\xi_{i};\xi_{j})}{\mathrm{P}(\xi_{i})\mathrm{P}(\xi_{j})}\delta_{A_{ij},1}+\left(1-\frac{{\langle}k{\rangle}}{N}\frac{\mathrm{W}(\xi_{i};\xi_{j})}{\mathrm{P}(\xi_{i})\mathrm{P}(\xi_{j})}\right)\delta_{A_{ij},0}\right]$
(85)

considered in (23), in the limit of large network size $N\to\infty$. We start
by computing the ensemble average

$\displaystyle\mathrm{P}(k|\boldsymbol{\xi})=\sum_{\boldsymbol{A}}\mathrm{P}(\boldsymbol{A}|\boldsymbol{\xi})\mathrm{P}(k|\boldsymbol{A})$
(86)

of the degree distribution for the single network instance

$\displaystyle\mathrm{P}(k|\boldsymbol{A},\boldsymbol{\xi})$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\delta_{k,\sum_{j(\neq i)}A_{ij}}$
(87)

Inserting (87) and (85) into (86), writing $\mathrm{P}({\bf
A}|\boldsymbol{\xi})=\prod_{i<j}\mathrm{P}(A_{ij}|\boldsymbol{\xi})$ and using
Fourier representations of Kronecker delta-functions, as well as the symmetric
property of the adjacency matrix $A_{ij}=A_{ji}$, we get

$\displaystyle\mathrm{P}(k|\boldsymbol{\xi})$
$\displaystyle=\frac{1}{N}\sum_{i}\sum_{\boldsymbol{A}}\int\frac{\rmd\omega}{2\pi}\rme^{\rmi\omega
k-\rmi\omega\sum_{j(\neq
i)}A_{ij}}\prod_{k<j}\mathrm{P}(A_{kj}|\boldsymbol{\xi})$ (88)
$\displaystyle=\frac{1}{N}\sum_{i}\sum_{\boldsymbol{A}}\int\frac{\rmd\omega}{2\pi}\rme^{\rmi\omega
k-\rmi\omega\sum_{k\neq
j}A_{kj}\delta_{ki}}\prod_{k<j}\mathrm{P}(A_{kj}|\boldsymbol{\xi})$
$\displaystyle=\frac{1}{N}\sum_{i}\sum_{\boldsymbol{A}}\int\frac{\rmd\omega}{2\pi}\rme^{i\omega
k-i\omega\sum_{k<j}A_{kj}(\delta_{ki}+\delta_{ji})}\prod_{k<j}\mathrm{P}(A_{kj}|\boldsymbol{\xi})$
$\displaystyle=\frac{1}{N}\sum_{i}\int\frac{d\omega}{2\pi}\rme^{i\omega
k}\prod_{k<j}\sum_{A_{kj}}\mathrm{P}(A_{kj}|\boldsymbol{\xi})\,\rme^{-i\omega
A_{kj}(\delta_{ki}+\delta_{ji})}.$

Taking the sum over $A_{kj}$ and employing the (asymptotic) identity
$1+x/N=\rme^{x/N+O(1/N^{2})}$, we have

$\displaystyle\mathrm{P}(k|\boldsymbol{\xi})$
$\displaystyle=\frac{1}{N}\sum_{i}\int\frac{\rmd\omega}{2\pi}\rme^{\rmi\omega
k}\prod_{k<j}\rme^{\frac{{\langle}k{\rangle}}{N}\frac{\mathrm{W}(\xi_{k};\xi_{j})}{\mathrm{P}(\xi_{k})\mathrm{P}(\xi_{j})}\left(\rme^{-\rmi\omega(\delta_{ki}+\delta_{ji})}-1\right)}$
(89)
$\displaystyle=\frac{1}{N}\sum_{i}\int\frac{\rmd\omega}{2\pi}\rme^{\rmi\omega
k}\prod_{k<j}\rme^{\frac{{\langle}k{\rangle}}{N}\frac{\mathrm{W}(\xi_{k};\xi_{j})}{\mathrm{P}(\xi_{k})\mathrm{P}(\xi_{j})}\left((\rme^{-\rmi\omega}-1)(\delta_{ki}+\delta_{ji})\right)}$
$\displaystyle=\frac{1}{N}\sum_{i}\int\frac{\rmd\omega}{2\pi}\rme^{\rmi\omega
k}\rme^{\frac{{\langle}k{\rangle}}{2N}\sum_{k\neq
j}\frac{\mathrm{W}(\xi_{k};\xi_{j})}{\mathrm{P}(\xi_{k})\mathrm{P}(\xi_{j})}\left((\rme^{-\rmi\omega}-1)(\delta_{ki}+\delta_{ji})\right)}$
$\displaystyle=\frac{1}{N}\sum_{i}\int\frac{\rmd\omega}{2\pi}\rme^{\rmi\omega
k}\rme^{(\rme^{-\rmi\omega}-1)\frac{{\langle}k{\rangle}}{N}\sum_{k(\neq
i)}\frac{\mathrm{W}(\xi_{k};\xi_{i})}{\mathrm{P}(\xi_{k})\mathrm{P}(\xi_{i})}}.$

At this point we make use of properties of the kernel
$\mathrm{W}(\xi;\xi^{\prime})$ as follows,

$\displaystyle\frac{{\langle}k{\rangle}}{N}\sum_{k(\neq
i)}\frac{\mathrm{W}(\xi_{k};\xi_{i})}{\mathrm{P}(\xi_{k})\mathrm{P}(\xi_{i})}$
$\displaystyle=\frac{{\langle}k{\rangle}}{N}\sum_{k(\neq
i)}\sum_{\xi}\delta_{\xi,\xi_{k}}\frac{\mathrm{W}(\xi;\xi_{i})}{\mathrm{P}(\xi)\mathrm{P}(\xi_{i})}$
(90)
$\displaystyle={\langle}k{\rangle}\sum_{\xi}\mathrm{P}(\xi)\frac{\mathrm{W}(\xi;\xi_{i})}{\mathrm{P}(\xi)\mathrm{P}(\xi_{i})}={\langle}k{\rangle}\frac{\mathrm{W}(\xi_{i})}{\mathrm{P}(\xi_{i})}$

where we have used the law of large numbers
$\mathrm{P}(\xi)=\lim_{N\to\infty}\frac{1}{N}\sum_{i=1}^{N}\delta_{\xi,\xi_{i}}$
and neglected ${\mathcal{O}}(N^{-1})$ terms. This allows us to derive the
following result

$\displaystyle\mathrm{P}(k|\boldsymbol{\xi})$
$\displaystyle=\frac{1}{N}\sum_{i}\int\frac{\rmd\omega}{2\pi}\rme^{\rmi\omega
k}\rme^{(\rme^{-\rmi\omega}-1){\langle}k{\rangle}\frac{\mathrm{W}(\xi_{i})}{\mathrm{P}(\xi_{i})}}$
(91)
$\displaystyle=\frac{1}{N}\sum_{i}\sum_{\xi}\delta_{\xi,\xi_{i}}\int\frac{\rmd\omega}{2\pi}\rme^{\rmi\omega
k}\rme^{(\rme^{-\rmi\omega}-1){\langle}k{\rangle}\frac{\mathrm{W}(\xi)}{\mathrm{P}(\xi)}}$
$\displaystyle=\sum_{\xi}\mathrm{P}(\xi)\,\rme^{-{\langle}k{\rangle}\frac{\mathrm{W}(\xi)}{\mathrm{P}(\xi)}}\int\frac{\rmd\omega}{2\pi}\rme^{\rmi\omega
k}\,\rme^{\rme^{-\rmi\omega}{\langle}k{\rangle}\frac{\mathrm{W}(\xi)}{\mathrm{P}(\xi)}}$
$\displaystyle=\sum_{\xi}\mathrm{P}(\xi)\,\rme^{-{\langle}k{\rangle}\frac{\mathrm{W}(\xi)}{\mathrm{P}(\xi)}}\int\frac{\rmd\omega}{2\pi}\rme^{\rmi\omega
k}\sum_{\ell\geq
0}\frac{({\langle}k{\rangle}\frac{\mathrm{W}(\xi)}{\mathrm{P}(\xi)})^{\ell}}{\ell!}\rme^{-\rmi\omega\ell}$
$\displaystyle=\sum_{\xi}\mathrm{P}(\xi)\,\rme^{-{\langle}k{\rangle}\frac{\mathrm{W}(\xi)}{\mathrm{P}(\xi)}}\frac{({\langle}k{\rangle}\frac{\mathrm{W}(\xi)}{\mathrm{P}(\xi)})^{k}}{k!}.$

The above can be simplified further. Starting from the definition of the
kernel $\mathrm{W}(\xi;\xi^{\prime})=\frac{\sum_{i\neq
j}A_{ij}\delta_{\xi,\xi_{i}}\delta_{\xi^{\prime},\xi_{j}}}{\sum_{i\neq
j}A_{ij}}$ we can write an expression for the marginal distribution
$\mathrm{W}(\xi)$,

$\displaystyle\mathrm{W}(\xi)$
$\displaystyle=\sum_{\xi^{\prime}}\frac{\sum_{i\neq
j}A_{ij}\delta_{\xi,\xi_{i}}\delta_{\xi^{\prime},\xi_{j}}}{\sum_{i\neq
j}A_{ij}}=\frac{\sum_{i\neq j}A_{ij}\delta_{\xi,\xi_{i}}}{\sum_{i\neq
j}A_{ij}}=\frac{\sum_{i}k_{i}\delta_{\xi,\xi_{i}}}{N{\langle}k{\rangle}}$ (92)
$\displaystyle=\frac{1}{N{\langle}k{\rangle}}\sum_{k}\sum_{i}k\delta_{\xi,\xi_{i}}\delta_{k,k_{i}}=\frac{1}{{\langle}k{\rangle}}\sum_{k}k\mathrm{P}(\xi,k)$
$\displaystyle=\frac{1}{{\langle}k{\rangle}}\mathrm{P}(\xi)\sum_{k}k\mathrm{P}(k|\xi)=\frac{\bar{k}(\xi)}{{\langle}k{\rangle}}\mathrm{P}(\xi)$

where we have defined $\bar{k}(\xi)=\sum_{k}\mathrm{P}(k|\xi)k$. We can now
substitute this back into our previous expression to find,

$\displaystyle\mathrm{P}(k|\boldsymbol{\xi})$
$\displaystyle=\sum_{\xi}\mathrm{P}(\xi)\rme^{-\bar{k}(\xi)}\frac{(\bar{k}(\xi))^{k}}{k!}.$
(93)

We note that the dependence on the specific realisation of $\boldsymbol{\xi}$
is lost and only the distribution $\mathrm{P}(\xi)$ enters the expression,
hence we will write $\mathrm{P}(k|\boldsymbol{\xi})=\mathrm{P}(k)$. We note
that by writing the degree distribution as,

$\displaystyle\mathrm{P}(k)=\sum_{\xi}\mathrm{P}(\xi)\mathrm{P}(k|\xi)$ (94)

comparison with (93) reveals the degree distribution of nodes with the same
label is the Poissonian distribution

$\displaystyle\mathrm{P}(k|\xi)=\rme^{-\bar{k}(\xi)}\frac{(\bar{k}(\xi))^{k}}{k!},$
(95)

with mean degree $\bar{k}(\xi)$. We also note that the joint distribution
$\mathrm{P}(k,\xi)$, which is required to solve (13), is obtained multiplying
the above times $\mathrm{P}(\xi)$

$\displaystyle\mathrm{P}(k,\xi)=\mathrm{P}(\xi)\rme^{-\bar{k}(\xi)}\frac{(\bar{k}(\xi))^{k}}{k!}.$
(96)

We can follow a similar procedure to find the average value of degree
correlations in the ensemble (85)

$\displaystyle\mathrm{W}(k,\xi;k^{\prime},\xi^{\prime}|\boldsymbol{\xi})$
$\displaystyle=\sum_{{\bf A}}\mathrm{P}({\bf
A}|\boldsymbol{\xi})\frac{1}{N{\langle}k{\rangle}}\sum_{ij}A_{ij}\delta_{k,k_{i}({\bf
A})}\delta_{k^{\prime},k_{j}({\bf
A})}\delta_{\xi,\xi_{i}}\delta_{\xi^{\prime},\xi_{j}}.$ (97)

By using Fourier representations of the Kronecker delta-function and taking
the sum over ${\bf A}$ we find, in the large $N$ limit,

$\displaystyle\mathrm{W}(k,\xi;k^{\prime},\xi^{\prime}|\boldsymbol{\xi})$
$\displaystyle=\mathrm{W}(\xi;\xi^{\prime})\mathrm{W}(k|\xi)\mathrm{W}(k^{\prime}|\xi^{\prime})$
(98)

with $\mathrm{W}(k|\xi)=k\mathrm{P}(k|\xi)/\bar{k}(\xi)$. This again reveals
that the degree correlations, when averaged over the graph ensemble, no longer
depend on the specific configuration $\boldsymbol{\xi}$, so we will write
$\mathrm{W}(k,\xi;k^{\prime},\xi^{\prime}|\boldsymbol{\xi})=\mathrm{W}(k,\xi;k^{\prime},\xi^{\prime})$.