In the field of network epidemiology, populations of hosts are modeled as a graph, where nodes are individuals and links are contacts along which infectious diseases can spread [1] . Here, we review important concepts from the study of random graphs to inform our understanding of epidemics. We explain how these concepts and results can be used on the one hand to inform interventions and guide policy, and on the other hand to solve common nonlinear models of epidemic spread based on differential equations. Our paper is structured such that the main text focuses on concepts and insights, with all mathematical details presented in complementary boxes and in the Appendix.
To model epidemic spreading with graphs, we define links or contacts as interactions that will transmit the disease from one individual to another should one individual becomes infectious. We can then analyze various compartmental epidemic models (e.g., SIR, SEIR, etc.) and how they depend on the structure of the resulting graphs [2] [3] [4] [5] [6] . For the simplest random graphs, epidemics end up being equivalent to classic, mass-action compartmentalized models in epidemiology, which can be analyzed using systems of differential equations. For instance, the SIR model on Erdős-Rényi random graphs, where every pair of nodes is independently connected with the same probability, behaves like the Kermack-McKendrick SIR solution in the limit where the number of nodes, i.e. the population size, is large [7] [8] [9] .
However, graphs allow us to model many kinds of heterogeneity and structure, and explore how and when heterogeneity affects epidemic thresholds and sizes. Since the degree of a node is the number of links or contacts it has, we can model superspreading events by having a heavy-tailed degree distribution, where some nodes have degree considerably larger than the mean [10] . A random contact is more likely to be with a high-degree individual, who has more secondary contacts in turn. Thus the average degree of a node hit by a random contact is greater than the average degree of a random node. This phenomenon is often called the "friendship paradox": on average, your friends have more friends than you do [11] . In the epidemic context it has two effects. First, the reproductive number is not simply the mean degree; it depends also on its second moment or variance (see Box I). In particular, even if we hold a. low-risk, but high-spread b. high-risk, but low-spread C. high-risk and high-spread the mean degree constant, increasing the variance lowers the critical transmission rate-i.e., the per contact probability of transmission needed to reach the epidemic threshold [2] . Second, as developed by [12] [13] [14] and recently pointed out as highly relevant for COVID-19 in [15, 16] , the "friendship paradox"can amplify the effectiveness of contact tracing, giving a better chance that a contact leads to a superspreading event, helping us prevent further cascades of infection. However, these results on superspreading events assume that the graph is undirected: that is, links are equally likely to spread the disease in either direction. While this may be true for individuals who attend a group event or work together in close quarters, some mechanisms of spread, such as surface contacts from handling and delivery, are inherently directed [17] . In network terms, an individual might have high out-degree (e.g. a high viral load if infected) but low in-degree (a low risk of receiving the infection from others), or vice versa. In addition, many non-pharmaceutical interventions to prevent disease spread have a directed effect. Wearing a simple mask is one such mechanism, as it is believed to reduce the probability that the wearer spreads a respiratory disease like COVID-19 to others, but not to strongly protect the wearer from receiving the disease from others [18, 19] .
For epidemics where directed effects are important, so that for some contacts transmission is likely in one direction but unlikely or impossible in the other, a better network model is a directed graph [4, 17, [20] [21] [22] . Here each node has an in-degree, i.e., the number of contacts from whom it could receive the disease, and an out-degree, the number of contacts to whom it could spread the disease if it becomes infected. We can think of the in-degree as the risk of this individual, and the out-degree as their spread. As illustrated in Fig. 1 , an individual might contribute to rapid spreading by having high risk, high spread, or both. For instance, an average immunocompromised individual might have high risk but low spread, while someone who fails to wear a mask even while those around them do might have low risk and high spread. Finally, someone who by choice or necessity interacts with many others in close quarters would have both high risk and high spread, as in the traditional undirected picture of superspreading events.
This distinction between risk and spread, or between incoming and outgoing links, has several consequences. First, the reproductive number is proportional to the correlation between risk and spread (see Box II). More precisely, the contribution of a given node to the reproductive number is proportional to the product of its in-and out-degrees. Nodes with both high risk and high spread are likely to become infected and then infect many others. But if they have high risk and low spread, or vice versa, then the number of new cases they are likely to generate is smaller. This has possible implications for non-pharmaceutical interventions: the effect of reducing an individual's spread is proportional to their risk, and vice versa.
Second, we now have two distinct types of friendship paradox (see Box III). Outgoing links connect to nodes with probability proportional to their in-degree, so following a link forward will tend to lead to individuals with high risk. Incoming links, in contrast, come from nodes with probability proportional to their out-degree, so going backward tends to lead to individuals with high spread.
If P (k) is the degree distribution, i.e., the fraction of nodes with degree k, then the average excess degree, or the expected number of secondary contacts resulting from a random contact, is
An epidemic occurs when R0 > 1. For the Poisson degree distribution we have Var k = k and thus R0 = k , but for distributions with larger variance we have R0 > k . This is the friendship paradox: a random friend of yours has more friends, on average, than you do (even not counting you).
In a directed graph, each node has an in-degree kin and an out-degree kout. We associate these respectively with their risk, the number of contacts from whom they could receive the infection, and their spread, the number of contacts to whom they could transmit it. The degree distribution is then a joint distribution, giving the fraction P (i, j) of nodes with kin = i and kout = j. Since the total number of incoming and outgoing links must match, the mean in-degree and the mean out-degree are equal:
In the directed configuration model of random graphs [23] , we assume that nodes are connected to each other randomly conditioned on their in-and out-degrees. The "heads" of the directed links are randomly matched with their "tails," and each outgoing link arrives at a given node with probability proportional to that node's in-degree. (This method of constructing a random graph can give rise to self-loops or parallel links, but in sparse networks these are rare.) In the limit of large population size, we can also say that for each pair of nodes a, b there is a link a → b with probability kout(a)kin(b)/m where m = k n is the total number of links and n is the total number of nodes. The reproductive number, or the expected number of secondary contacts resulting from a random contact, is
Just as R0 depends on the variance of the degree distribution in the undirected case, it now depends on the covariance, or correlation, between the in-and out-degrees. In particular, the contribution of each individual to the reproductive number is proportional to the product of their risk and their spread.
This suggests that to find superspreading events, and thus to prevent a large number of downstream cases, it may be important to perform backward contract tracing, determining the source of an index case, rather than focusing only on forward contact tracing, finding those who the index case may have infected. As shown in Box III, the average spread of the predecessor of an index case is larger than that of an average successor, especially if the spread has high variance, except in the unusual case where risk and spread are almost perfectly correlated. This insight also implies that the time window we use to define relevant contacts should include their likely predecessor as well as the likely successors of an index case. We can then trace the "siblings" of the index case, i.e., others who were infected by the same source, and prevent a larger number of future cases.
Here we review some results of Refs. [24, 25] from random graph theory which show how to calculate the epidemic size of SIR-like models with directed transmission, as a function of the joint distribution of risk and spread. We note that these results can be generalized [26, 27] and expanded to include semi-directed networks with a mixture of directed and undirected links [17, 28] .
As mentioned above, we define the graph as the links a → b that will transmit the disease from a to b if a becomes infected. More generally this transmission would occur with some probability, so these links form a randomly chosen subnetwork of the potentially infectious contacts. The probabilistic case can be handled with the same methods, and we give an example in Box VI in the Appendix.
Let us first compute the epidemic probability, i.e., the probability that a single infected individual generates a macroscopic epidemic as opposed to a small outbreak. It helps to think of a process where links, rather than nodes, propagate from one generation to the next. An outgoing link produces j "children," or secondary contacts, if the link arrives at a node with out-degree j. A classic question in the theory of branching processes is whether the resulting tree of contacts is finite or infinite. Since the total number of descendants of a link is finite if and only if this is true of all of its children, the probability of this event obeys a fixed-point equation which can be written compactly using probability generating functions (see Box IV). Now consider the epidemic size. Nodes receive the infection from their predecessors, who receive it from their predecessors in turn. This leads to a backwards branching process of incoming links, where a link has i "children" (or perhaps we should call them parents now) if it comes from a node with in-degree i. If a node has only a finite number of ancestors in this tree, then with high probability none of them will be the initial node and it will remain uninfected. But if it has a macroscopic (i.e. infinite) number of ancestors, then with high probability one of its ancestors will Box III: Forward and backward friendship paradoxes and contact tracing Suppose we are performing contact tracing from an index case. Since the probability a given individual is their predecessor (the source of their infection) is proportional to their out-degree or spread, the average excess spread of their predecessor-that is, the number of "siblings" with the same predecessor-is
Similarly, the probability a given individual is a successor of the index case is proportional to their in-degree or risk, so the average excess risk of a successor is
These are the backward and forward friendship paradoxes: the predecessor of a random node tends to have high spread, but a successor tends to have high risk, especially when both of these have high variance. On the other hand, the average spread of a successor is just the reproductive number calculated in (3) above,
Comparing this with (4), we see that if the spread kout has high variance then the predecessor has larger spread on average than a successor, unless the risk and spread are almost perfectly correlated, i.e., unless kinkout ≥ k 2 out − 1. This insight highlights the conclusion of Ref. [16] that contact tracing strategies should include tracing backwards from the index case as well as forwards. In fact, this recommendation becomes critical once we consider the directionality of contacts such that risk and spread might not be perfectly correlated.
become infected and it will too. Thus the epidemic size is equal to the epidemic probability in a time-reversed model, where the disease spreads backward from each node to its predecessors.
On undirected graphs, these forward and backward branching processes are the same, and the probability they die out obeys the same fixed-point equation. This leads to the observation in the SIR model that the epidemic size equals the epidemic probability. In the directed case, however, these are typically different, a fact which was pointed out in Ref. [17] . In the next section, we consider degree distributions of various kinds and with varying levels of correlation between risk and spread, and look at how these types of heterogeneity effect their epidemic size.
Let us first consider a simple, random, undirected graph, as illustrated in Fig. 2(a-b) . Its structure is fully specified by its degree distribution, which sets the number of contacts leading to and from a given node, such that k in = k out = k. As shown in Ref. [29] , the epidemic size obtained on these graphs corresponds to the classic Kermack-McKendrick solution if and only if the degree distribution is Poisson, i.e., P (k) = k k e − k /k!. What is so special about this distribution? As it turns out, the Poisson distribution is the only one where the excess degree has the same distribution as the degree itself: the probability a link leads to k secondary contacts equals the probability that a random node has k contacts. In terms of the friendship paradox, the probability that someone has k friends is the same as the probability their friend has k friends in addition to you. In particular, these two distributions have the same mean, so R 0 = k . This is similar to the mass-action hypothesis of Kermack and McKendrick since it implies that information about the contact structure does not matter. For any other degree distribution -even a more homogeneous one as in Fig. 2 (a) -the distribution of excess degrees is different from the degree distribution, and the Kermack-McKendrick solution for the epidemic size does not apply.
Consider now a directed version of the Poisson graph. As shown in Fig. 2(c) , we can think of undirected links as bidirectional contacts, i.e., pairs of directed links, and make this directionality explicit before randomizing the incoming and outgoing links. What we now have is a graph where k in = k out = k, still drawn from a Poisson distribution, but where incoming and outgoing neighbours are different. Even if the construction of the directed graph structure is very close to that of the undirected case, this directed version actually does not fall back on the Kermack-McKendrick solution. One key difference can be observed around transmission dead-ends, i.e., nodes with k out = 0, which are now impossible to reach. In other words, if a node has been reached via one if its incoming links, there is always at least one outgoing link available to reach the next node. As shown in the Appendix, equivalence Box IV: Generating functions and epidemic size Given the degree distribution P (i, j), we can define a probability generating function (PGF)
Many quantities can be written compactly in terms of G and its derivatives, including the average in-and out-degrees (2) and the reproductive number (3):
We can also use G to analyze branching processes on the network. In the forward branching process of an outbreak, the probability Qout(j) that a random contact results in j secondary contacts, i.e., that a random link arrives at a node of out-degree j, is
We can write a generating function for Qout in terms of G,
Now let uout be the probability that a link has a finite number of descendants, i.e., that it will lead to a finite outbreak rather than an epidemic. This is true if it is true for all j of its children. Since these are independent, averaging over j gives the fixed-point equation
of which uout is the smallest positive solution. Finally, since a node with out-degree j generates an epidemic with probability 1 − u j out , the probability that infecting a random initial node leads to an epidemic is
To compute the epidemic size, we consider the backwards branching process of incoming links. The probability Qin(i) that a link comes from a node with in-degree i, and the corresponding generating function, is
The probability uin that a link has a finite number of ancestors is then the smallest positive solution of
Finally, a node with in-degree i becomes infected with probability 1 − u i in , so the size of the epidemic is It turns out that some naive approaches to modeling heterogeneity do not, in fact, alter the epidemic size at all. Consider the following kind of simulation, which would be run until there are no new infectious individuals and where P out is an arbitrary distribution of out-degrees: This type of simulation can give transmission trees similar to the example shown in Fig. 2(d) . While the out-degree distribution P out can be arbitrarily heterogeneous, including various kinds of superspreading, the incoming links fall with equal probability on every susceptible individual because of the random selection of neighbours in step 3. As a result, the distribution of risk or in-degree is in fact Poisson, and the in-and out-degrees are uncorrelated. As shown in the Appendix, in the limit of large N , the typical epidemic size produced by this kind of simulation exactly corresponds to the Kermack-McKendrick solution regardless of the out-degree distribution P out .
These simple examples suggest that the Kermack-McKendrick solution relies on (i) a Poisson distribution of risk among individuals and (ii) a lack of correlation between risk and spread. It does not actually rely on homogeneous spread but only on the independence between the infection risk faced by an individual and the number of secondary infections caused by that individual if infected. The undirected Poisson case, where this independence does not hold since the in-and out-degrees are identical by definition, is somewhat unique since the lack of backtracking and the friendship paradox balance out. We formalize this intuition mathematically in the Appendix, and some important results are visually summarized in Fig. 3(left) .
The mathematical analysis presented thus far is completely general and exact on very large random directed graphs. To further assess the impact of risk heterogeneity on the epidemic size, we need to look into the mechanisms behind the joint distribution of both risk and spread in specific models. Here we show how one case of these distributions arises from age-structured models in epidemiology.
We consider a general case in which there are M types of nodes and where a fraction w of the nodes are of type . We assume that the in-degrees and the out-degrees are uncorrelated, and that the contribution to the in-degree from each node type is independent and distributed according to a Poisson distribution. Thus, a PGF G is defined for each node type and has the general form
where z in r is the average number of incoming links type-nodes received from nodes of type r. The average in-degree of type-nodes is therefore r z in r . Akin to Eq. (2), the average in-and out-degrees between any pair of node types are related via
Following Ref. [31] , the size of a macroscopic outbreak is
where u in = (u in 1 , . . . , u in M ) is the largest positive solution of
for = 1, . . . , M . The reproductive number, R0 = ρ(Z), is equal to the spectral radius, i.e., the largest singular value, of the matrix Z whose elements are z in r .
One particularly powerful and common approach in current models of the COVID-19 pandemic is age-structured SIR or SEIR models (e.g., Ref. [30, 32, 33] ) where individuals are separated in different types based on their age and where contacts between age groups are specified through a contact (or "mixing") matrix. This approach is powerful because it allows the modeller to include important sources of heterogeneity such as different susceptibilities and likelihood of symptoms based on age, and also allows an easy parametrization of interventions such as school closures and/or lockdowns through the contact matrix. Note that the different categories need not be related to age, and could be used to represent other distinguishing properties such as behavior, physiology, employement, etc.
For each group i, let us define S i (t), I i (t) and R i (t) as the numbers of Susceptible, Infectious, and Recovered individuals in that group at time t. These structured models follow the SIR dynamics through a coupled system of ordinary differential equations written asṠ
where β is the baseline transmission rate of the disease, σ i is the group-dependent susceptibility of individuals, M ij is the average number of contacts from someone in group i to people in group j, and γ is the recovery rate of the disease. Note that these equations are equivalent to the model described in the Supplementary Materials of Ref. [30] . They can be understood as setting a Poisson distribution of in-degree with average βσ i M ij /γ coming from group j for all individuals in group i. The model therefore corresponds to a mixture of directed Poisson processes, and it is unclear a priori how much we should expect these models to differ from the simple mass-action SIR, which is equivalent to an undirected, unstructured network. But note in particular that the matrix M ij is not generally symmetric, so these equations implicitly recognize that links in the social network are directed. We apply this framework to a multitype directed graph [34] (see Box V) where the type of a node corresponds to their age group and where directed links allow us to consider heterogeneous susceptibility. As shown in Fig. 3(right) , the multitype directed graph approach exactly solves for the final epidemic size produced by the system presented in Eq. (20) . This provided a general solution for the epidemic size in structured SIR-like models. Perhaps more importantly, this solution highlights how the solution of such systems is actually a simple mixture of Poisson processes and therefore a combination of Kermack-McKendrick-like terms. As can be observed in Fig. 3(right) , including this type of heterogeneity does not translate into epidemic sizes that are very different from the ones predicted by the more simpler Kermack-McKendrick solution.
One intuition for this is that, since the sum of independent Poisson variables is Poisson (with the sum of their means) the total risk of individuals in these models, i.e., the number of incoming links they have from all age groups, is Poisson-distributed. Thus, while these models are structured, they result in a sum of terms like that in the original Kermack-McKendrick model, where R 0 is given by a weighted average (specifically a matrix eigenvalue, see Box V). As a result, these models cannot truly embrace the full heterogeneity of risk and spread produced by biology, demographics, and human behaviour.
In this paper, we reviewed key results from random directed graphs to clarify the role of heterogeneity and correlations in risk and spread of an infectious disease. We noted that in directed networks the friendship paradox plays different roles in contact tracing when following links forward or backward: forward tracing leads to high risk individuals, and backward tracing leads to high-spread individuals.
We then focused on identifying when the final epidemic size deviates from the classic Kermack-McKendrick solution. In particular, we showed that implementing heterogeneity in differential equation models by having a small number of types of individuals (representing, for example, age, occupation, susceptibility, and/or activity level) does not significantly affect epidemic size. Risk in these models is a mixture of Poisson distributions which individually correspond to Kermack-McKendrick-like contributions and therefore result in a similar epidemic size. While contact patterns between groups may be important in assessing interventions (e.g. the effect of school closings) the epidemic size is not always as sensitive to these parameters as one might think.
The robustness of the Kermack-McKendrick result on directed networks highlights the fact that the critical assumption behind their analysis of the SIR model is not mass-action mixing per se, but instead an implicit assumption that risk is Poisson distributed. As we showed in Section III, this also implies that any model or simulation will fall back on the classic mass-action result if the secondary cases caused by an infectious individual are randomly selected from the population, no matter how heterogeneous the distribution of spread is. When there is genuine heterogeneity in both risk and spread, they play dual roles in the early and late stages of an epidemic: the probability that an outbreak becomes an epidemic is driven primarily by superspreading events [10, 35] , while the final epidemic size is driven primarily by high-risk individuals.
For emerging infectious diseases, especially those like COVID-19 where no proven pharmaceutical therapies nor prophylactics exist, our results highlight the critical need for detailed, contact tracing studies. These studies should quantify both the incoming and outgoing risk of transmission and quantify the relative contribution of these risks during the epidemic. Indeed, the efficacy of measures such as digital contact tracing [36] and "immune shielding" [37] will depend on how risk and contact heterogeneity are structured across populations. Only through the application of models able to capture relevant heterogeneity, and the collection of empirical data on contact networks and risk, can the efficacy of such non-pharmaceutical interventions be accurately determined. In this Appendix, we formally solve for the final epidemic size found in directed graphs representing different models of disease spread, including the examples used in Section III.
Let us first illustrate the formalism when in-and out-degree are independent. Denoting their respective PGFs by G in (x) and G out (y), Eqs. (7), (11) , (12) , (14) and (15) become
and u out = G out (u out ) ;
We conclude that the size of a macroscopic outbreak, S, is solely dictated by the distribution of the in-degrees, i.e. by the distribution of the risk of individuals, while the probability for an outbreak to become macroscopic depends only on the distribution of spread. Moreover, we see from Eq. (3) that the number of secondary infections from patient zero and from secondary cases are equal, that is
and conclude therefore that there is no "friendship paradox effect" (from the excess degree distribution) when in-and out-degree are not correlated.
In addition to being independent from the out-degrees, if we assume that the in-degrees are distributed according to a Poisson distribution, that is if G in (x) = e k ( We conclude that the predictions for the size of a macroscopic outbreak will coincide with the Kermack-McKendrick solution if the in-degrees and not correlated with the out-degrees and if they are distributed according to a Poisson distribution.
If the in-degree is Poisson-distributed with mean k , its PGF is
As discussed in the text, with R0 = k this leads to the same prediction for the epidemic size as Kermack-McKendrick. Another degree distribution that might result from an epidemic simulation is where each individual has d contacts, each of which transmits the disease with probability q. Then k = qd and the distribution of in-degrees is binomial, giving the PGF
When d is large, the identity (1 + w/d) d ≈ e w applies, giving the same PGF (A.6) as the Poisson case and again recovering the Kermack-McKendrick prediction for epidemic size.
Let us now investigate the effect of a perfect correlation between individual risk and spread, which can be encoded as
where P in (i) is the in-degree, or risk, distribution. Note that full correlation does not imply that links are reciprocal (i.e. links run in both directions between nodes, which would be equivalent to undirected graphs). Equations (11), (12) , (14) and (15) We conclude that in the presence of full correlation, the size and the probability of a macroscopic outbreak are equal. Also, Eq. (3) shows that
which differs from its undirected counterpart (i.e. R 0 = k(k − 1) / k [24] ). As shown in Ref. [29] , the undirected case, in which in-and out-degrees are perfectly correlated by construction, falls back on the result of Kermack-McKendrick when the degrees are distributed according to a Poisson distribution. However, this results does not extend to fully-correlated directed networks. Indeed, substituting G in (x) = e k (x−1) in Eq. (A.9b) yields only two solutions u in = 0, 1 for all k . The probability for a macroscopic outbreak and its size are therefore equal to the fraction of nodes with nonzero in-/out-degree S = P = 1 − e − k , (A.11) and Eq. (A.10) yields R 0 = k + 1, since a node that has been infected can always infect at least one other node.
To explore scenarios in-between the limiting correlation cases above, we consider the following PGF which allows to fix the marginal out-degree distribution and to interpolate between full correlation (a = 1) and independence (a = 0) G(x, y) = [G out (xy)] a [G in (x)G out (y)] 1−a .
(A.12)
Note that in the case of full correlations, the in-degrees are distributed according to G out . Equations (11), (12) , (14) and (15) 