Takeover times for a simple model of network infection

TL;DR

This paper analyzes the time it takes for a simple stochastic infection model to spread across different network types, revealing how the distribution of takeover times varies with network structure.

Contribution

It provides a detailed characterization of the distribution of infection takeover times across various network topologies, linking extremal behavior to classical probability distributions.

Findings

Gumbel distribution for star graphs

Sum of Gumbels for complete and Erdős-Rényi graphs

Normal distribution for 1D and 2D lattices

Abstract

We study a stochastic model of infection spreading on a network. At each time step a node is chosen at random, along with one of its neighbors. If the node is infected and the neighbor is susceptible, the neighbor becomes infected. How many time steps $T$ does it take to completely infect a network of $N$ nodes, starting from a single infected node? An analogy to the classic "coupon collector" problem of probability theory reveals that the takeover time $T$ is dominated by extremal behavior, either when there are only a few infected nodes near the start of the process or a few susceptible nodes near the end. We show that for $N \gg 1$, the takeover time $T$ is distributed as a Gumbel for the star graph; as the sum of two Gumbels for a complete graph and an Erd\H{o}s-R\'{e}nyi random graph; as a normal for a one-dimensional ring and a two-dimensional lattice; and as a family of intermediate skewed distributions for $d$-dimensional lattices with $d \ge 3$ (these distributions approach the sum of two Gumbels as $d$ approaches infinity). Connections to evolutionary dynamics, cancer, incubation periods of infectious diseases, first-passage percolation, and other spreading phenomena in biology and physics are discussed.

Equations

$$p_m=$$

$$\times \left(\text{mbox}Prob.ofselectingsusceptibleneighbor\right).$$

$$P(X_m=k)=q_m^{k-1}{p}_m,$$

$$T=\sum _{m=1}^{N-1}{X}_m,$$

$$\frac{T-N(N-1)}{(N-1)\sqrt{N}}\stackrel{d}{\to }\text{mbox}Normal(0,1),$$

$$p_m=\frac{1}{N+1}\cdot \frac{N-m}{N}$$

$$P(\mathcal{E}(p_m)=x)dx=p_m{e}^{-p_{m}x}dx;x\ge 0.$$

$$\frac{T-\mu }{L}\sim \frac{F-\mu }{L},$$

$$F/L=\sum _{m=0}^{N-1}\mathcal{E}(p_m/L)=\sum _{k=1}^N\mathcal{E}(k)$$

$$g_N(x)=\sum _{k=1}^{N}k{e}^{-kx}\prod _{r\ne k}^N\frac{r}{r-k},$$

$$g_N(x)=N{e}^{-x}(1-e^{-x}{)}^{N-1}.$$

$$f(x)=e^{-(x+\gamma )}\exp \left(-e^{-(x+\gamma )}\right)$$

$$h(x)={\beta }^{-1}{e}^{-(x-\alpha )/\beta }\exp \left(-e^{-(x-\alpha )/\beta }\right).$$

$$\frac{T-\mu }{N(N+1)}\stackrel{d}{\to }G,$$

$$p_m=\frac{m}{N}\left(1-\frac{m-1}{N-1}\right)=\frac{m}{N}\cdot \frac{N-m}{N-1}$$

$$\frac{T-\mu }{N}\sim \sum _{m=1}^{N-1}\frac{\mathcal{E}(p_m)-1/p_m}{N}=:S(p).$$

$$p_k=\frac{k}{N}\cdot \frac{N-k}{N-1}=\frac{N-m}{N}\cdot \frac{m}{N-1}=p_m.$$

$$S(p)=S\left(p^{(f)}\right)+S\left(p^{(b)}\right).$$

$$S(r)\stackrel{d}{\to }G$$

$$p^{(T)}=(p_m{)}_{m=1}^B,\text{mbox}and{r}^{(T)}=(r_m{)}_{m=1}^B,$$

$$S(p^{(f)})\sim S(p^{(T)}),S(p^{(b)})\sim S(p^{(T)})$$

$$S(r)\sim S(r^{(T)}).$$

$$S(r^{(T)})\sim S(p^{(T)}).$$

$$\frac{T-\mu }{N}\stackrel{d}{\to }G+G.$$

$$\eta =\frac{d-1}{d}.$$

$$p_m\approx \frac{m^{\eta }}{N}\left(1-\frac{m^{\eta }}{N}\right).$$

$$\frac{T-\mu }{\sigma }\sim \sum _{m=1}^N\frac{X(p_m)-1/p_m}{{\sigma }_X},$$

$$\frac{T-\mu }{\sigma }\sim \frac{1}{\sqrt{2}}(F^{\prime }+F^{\prime })$$

$$F^{\prime }:=\sum _{m=1}^M\frac{\mathcal{E}(m^{\eta })-1/m^{\eta }}{\sqrt{H}}$$

$$H:=H(2\eta ;M)=\sum _{m=1}^{M}1/m^{2\eta }$$

$$\frac{T-\mu }{\sigma }\stackrel{d}{\to }\text{mbox}Normal(0,1)$$

Full text

I Introduction

Contagion is a topic of broad interdisciplinary interest. Originally studied in the context of infectious diseases Anderson and May (1991); Keeling and Eames (2005); Diekmann et al. (2012); Pastor-Satorras et al. (2015), contagion has now been used as a metaphor for diverse processes that spread by contact between neighbors. Examples include the spread of fads and fashions Bikhchandani et al. (1992); Watts (2002), scientific ideas Bettencourt et al. (2006), bank failures Aharony and Swary (1983); Allen and Gale (2000); May et al. (2008); May and Arinaminpathy (2010); Haldane and May (2011), computer viruses Kephart and White (1991), gossip Haas et al. (2006), rumors Daley and Kendall (1965); Draief and Massouli (2010), and yawning Provine (2005). Closely related phenomena arise in probability theory and statistical physics in the setting of first-passage percolation Aldous (2013); Auffinger et al. (2015), and in evolutionary dynamics in connection with the spread of mutations through a resident population Lieberman et al. (2005); Antal et al. (2006); Hindersin and Traulsen (2014, 2015); Ashcroft et al. (2015). We will use the language of contagion throughout, but bear in mind that everything could be reformulated in the language of the other fields mentioned above.

In the simplest mathematical model of contagion, the members of the population can be in one of two states: susceptible or permanently infected. When a susceptible individual meets an infected one, the susceptible immediately becomes infected. Even in this idealized setting, interesting theoretical questions remain, whose answers could have significant real-world implications, as we will argue below.

For example, consider the following model, motivated by cancer biology. Imagine a two-dimensional lattice of cells in a tissue, where each cell is either normal or mutated. At each time step a random cell is chosen, along with one of its neighbors, also chosen uniformly at random. If the first cell is mutated and its neighbor is normal, the mutated cell (which is assumed to reproduce much faster than its normal neighbor) makes a copy of itself that replaces the normal cell. In effect, the mutation has spread; it behaves as if it were an infection. This deliberately simplified model was introduced in 1972 to shed light on the growth and geometry of cancerous tumors Williams and Bjerknes (1972).

Here, we study this model on a variety of networks. Our question is: given a single infected node in a network of size $N$, how long does it take for the entire network to become infected? We call this the takeover time $T$. It is conceptually related to the fixation time in population genetics, defined as the time for a fitter mutant to sweep through a resident population. It is also reminiscent of the incubation period of an infectious disease, defined as the time lag between exposure to the pathogen and the appearance of symptoms; this lag presumably reflects the time needed for infection to sweep through a large fraction of the resident healthy cells.

For the model studied here, the calculation of the network takeover time is inherently statistical because the dynamics are random. At each time step, we choose a random node in the network, along with one of its neighbors, also at random. If neither of the nodes are infected, nothing happens and the time step is wasted. Likewise, if both are infected, the state of the network again doesn’t change and the time step is wasted. Only if the first node is infected and its neighbor is susceptible does the infection progress, as shown in Figure 1.

The time course of the infection is interesting to contemplate. Intuitively, when the network is large, it seems that the dynamics should be very stochastic at first and take a long time to get rolling, because it is exceedingly unlikely that we will randomly pick the one infected node, given that there are so many other nodes to choose from. Similarly we expect a dramatic slowing down and enhancement of fluctuations in the endgame. When a big network is almost fully infected, it becomes increasingly difficult to find the last few susceptible individuals to infect.

These intuitions led us to suspect that the problem of calculating the distribution of takeover times might be amenable to the techniques used to study the classic “coupon collector” problem in probability theory Feller (1968); Pósfai (2010). If you want to collect $N$ distinct coupons, and at each time step you are given one coupon at random (with replacement), what is the distribution of the time required to collect all the coupons? Like the endgame of the infection process, the coupon collection process slows down and suffers large fluctuations when almost all the coupons are in hand and one is waiting in exasperation for that last coupon. Erdős and Rényi proved that for large $N$, the distribution of waiting times for the coupon collection problem approaches a Gumbel distribution Erdős and Rényi (1961). This type of distribution is right skewed, and is one of the three universal extreme value distributions Fisher and Tippett (1928); Kotz and Nadarajah (2000).

In what follows, we will show that for $N\gg 1$, the takeover time $T$ is distributed as a Gumbel for the star graph, and as the sum of two Gumbels for a complete graph and an Erdős-Rényi random graph. For $d$-dimensional cubic lattices, the dependence on $d$ is intriguing: we find that $T$ is normally distributed for $d=1$ and $d=2$, then becomes skewed for $d\geq 3$ and approaches the sum of two Gumbels as $d$ approaches infinity. We conclude by discussing the many simplifications in our model, with the aim of showing how the model relates to more realistic models. We also discuss the possible relevance of our results to fixation times in evolutionary dynamics, population genetics, and cancer biology, and to the longstanding (yet theoretically unexplained) clinical observation that incubation periods for infectious diseases frequently have right-skewed distributions.

II One-dimensional Lattice

We start with a one-dimensional (1D) lattice. In this paper, we will always take lattices to have periodic boundary conditions, so imagine $N$ nodes arranged into a ring.

Suppose that $m$ nodes are currently infected. Let $p_{m}$ denote the probability that a susceptible node gets infected in the next time step. Notice that for a more complicated graph, $p_{m}$ might not be a well-defined concept, because it could depend on more than $m$ alone: the probability of infecting a new node could depend on the positions of the currently infected nodes, as well as on the susceptible node being considered. In such cases, we would need to know the entire current state of the network, not just the value of $m$, to calculate the probability that the infection will spread.

The 1D lattice, however, is especially tractable. Assuming that only one node is infected initially, at later times the infected nodes are guaranteed to form a contiguous chain. So for this simple case the graph state is indeed determined by $m$ alone. The only places where the infection can spread are from the two ends of the infected chain. (Even on more complicated networks, the dynamics of our model imply that the infected nodes always form contiguous regions, but few are as simple as this.)

The spread of infection involves two events. First, the node chosen at random must lie on on the boundary of the infected cluster. Then, one of its neighbors that happens to be susceptible must be picked. So

[TABLE]

Hence, for the ring, the probability that the infection spreads on the next time step reduces to $p_{m}=(2/N)\times(1/2)=1/N$ for all $m$.

Next, define the random variable $X_{m}=X(p_{m})$ as the number of time steps during which the network has exactly $m$ infected nodes. The probability that this state lasts for $k$ time steps is then given by

[TABLE]

for $k=1,2,\ldots,$ where $q_{m}:=1-p_{m}$. To see this, note that $P(X_{m}=k)$ is the probability that no new infection occurs on the first $k-1$ steps, times the probability that infection does occur on step $k$.

Thus, for any network where $p_{m}$ is well defined, the time spent with $m$ infected nodes is a geometric random variable, with mean $1/p_{m}$ and variance $1/p_{m}^{2}-1/p_{m}$. In particular, since the ring has $p_{m}=1/N$ for all $m$, we find that $X_{m}$ has mean $N$ and variance $N^{2}-N$ in this case.

The takeover time for any network is

[TABLE]

the sum of all the individual times required to go from $m$ to $m+1$ infected nodes, for $m=1,\ldots,N-1$. (Equality, in this case, means equality in distribution, as it will for all the other random variables considered throughout this paper.)

In the case of the 1D lattice, all the $X_{m}$ are identical. However, their means and variances depend on $N$, which prevents us from invoking the usual Central Limit Theorem to deduce the limiting distribution of $T$. However, we can invoke a generalization of it known as the Lindeberg-Feller theorem. See Appendix A for more details.

After normalizing $T$ by its mean, $\mu=N(N-1)$, and its standard deviation, $\sigma=(N-1)\sqrt{N}$, we find

[TABLE]

where the symbol $\xrightarrow{d}$ means convergence in distribution as $N$ gets large. Figure 2 confirms that the takeover times are normally distributed in the limit of large rings.

III Star graph

A star graph is another common example for infection models. Here, $N$ separate “spoke” nodes all connect to a single “hub” node and to no others, as illustrated in the upper left of Figure 3. We will assume the initial infection starts at the hub, since starting it in a spoke node would require only a trivial adjustment to the calculations below.

Let $m$ be the number of spoke nodes that are currently infected. As in the ring case, $p_{m}$ (the probability to go from $m$ to $m+1$ in the next time step) is a well-defined quantity that depends on $m$ alone, and not on any other details of the network state. Using the logic of Eq. (1), we get

[TABLE]

for $m=0,1,\ldots,N-1$. Here, $1/(N+1)$ is the probability of choosing the infected hub as the first node, and $(N-m)/N$ is the probability of selecting one of the $N-m$ currently susceptible spoke nodes, out of the $N$ spoke nodes in total, as its neighbor.

Now that $p_{m}$ is in hand for the star graph, we can define the random variable $X_{m}$ and the takeover time $T$ just as we did for the one-dimensional ring. The only difference is that the $m$-dependence of $p_{m}$ is now controlled entirely by the factor $(N-m)/N$.

That same factor turns up in a classic probability puzzle called the coupon collector’s problem Feller (1968); Pósfai (2010). At the time of this writing, millions of children are experiencing it firsthand as they desperately try to complete their collection of pocket-monsters in the Pokémon videogame series.

To see the connection, suppose you are trying to collect $N$ distinct items, and you have $m$ of them so far. If you are given one of the $N$ items at random (with replacement), the probability it is new to your collection is $(N-m)/N$, the same factor we saw above, and precisely analogous to the probability $p_{m}$ of adding a new node to the infected set. Likewise, the waiting time to collect all $N$ items is precisely analogous to the time $T$ needed to take over the whole star graph. The only difference is the constant factor $1/(N+1)$ in Eq. (3).

The limiting distribution of the waiting time for the coupon collector’s problem is well known. Although it resembles a lognormal distribution Read (1998), in fact it is a Gumbel distribution in the limit of large $N$, given the right scaling  Erdős and Rényi (1961); Pósfai (2010); Rubin and Zidek (1965); Baum and Billingsley (1965). We will now show that the same is true for our problem.

The first move is to approximate the geometric random variables $X(p_{m})$ by exponential random variables $\mathcal{E}(p_{m})$, with density

[TABLE]

From here we define the random variable $F:=\sum_{m=0}^{N-1}\mathcal{E}(p_{m})$, which has mean $\mu=\sum_{m=0}^{N-1}1/p_{m}$.

It can be shown (see Appendix B) that for a large class of $p_{m}$ and normalizing factors $L:=L(N)$ that

[TABLE]

where the symbol “$\sim$” means the ratio of characteristic functions goes to 1 as $N$ gets large. That is, the random variables on both sides converge to each other in distribution as $N$ gets large.

In the traditional coupon collector’s problem we would take $L=N$; but because of that $(N+1)$ factor, what we want is $L=N(N+1)$. Thanks to the fact we are now using exponential variables, we now know

[TABLE]

(using $k=N-m=1,2,\ldots,N$). A nice closed form for the distribution of the sum of a collection of distinct exponential variables is known Ashcroft et al. (2015); Baum and Billingsley (1965), but for convenience’s sake we rederive it in Appendix C. When we have any finite $N$ we can simply write the relevant distribution function $g_{N}(x)$ for $x\geq 0$ as

[TABLE]

which can be manipulated into

[TABLE]

From here, we can find the distribution for $(F-\mu)/L$, and $(T-\mu)/L$ by extension. Therefore, taking the limit of large $N$ and using the standard approximation of the harmonic sum for $\mu$ gives us

[TABLE]

as the density, where $\gamma\approx 0.5772$ is the Euler-Mascheroni constant. The density in (5) is a special case of the Gumbel distribution, denoted $\mathrm{Gumbel}(\alpha,\beta)$ and defined to have the density

[TABLE]

Specifically, we find

[TABLE]

where $G$ is a Gumbel random variable distributed according to $\mathrm{Gumbel}(-\gamma,1)$.

This distribution can be tested against simulation, and it works nicely as seen in Figure 3. Gumbel distributions have arisen previously in infection and birth-death models Ashcroft et al. (2015); Gautreau et al. (2007); Williams (1965), and are well known in extreme-value theory Feller (1968); Fisher and Tippett (1928); Kotz and Nadarajah (2000), but the fact that they show up here as a result of a network topology is unexpected.

IV Complete graph

The complete graph on $N$ nodes corresponds to a “well-mixed population” and is one of the most common topologies in infection models. This network consists of $N$ mutually connected nodes, so the location of the initial infection does not matter.

Given $m$ infected nodes, we once again have a well-defined $p_{m}$. Using the concept behind Eq. (1), we find

[TABLE]

for $m=1,\ldots,N-1$. For the sake of convenience, we will collect these probabilities into a vector $p=(p_{m})_{m=1}^{N}$. As in the case of the star graph, we can approximate the takeover time $T$ by summing exponential random variables instead of geometric ones. So

[TABLE]

To compress notation, we defined $S(p)$ to be the normalized sum of exponential random variables across the entries of the vector $p$.

The specific $p$ in Eq. (8) has some helpful symmetry. Notice that if $k=N-m$, then

[TABLE]

This symmetry means that the second half of the takeover looks just like the first half played backwards. If we set $p^{(f)}$ to be the front half of the $p$-vector and $p^{(b)}$ to be the back half of $p$, then we know

[TABLE]

Because we have a symmetry and the order we add the individual exponential variables will not matter, the random variables $S(p^{(f)})$ and $S(p^{(b)})$ should be equal in distribution. Although $N$ being odd or even may seem to be distinct cases, we will find that the distinction does not matter.

The basic concept here is to compare $p^{(f)}$ and $p^{(b)}$ to $r$, where $r=(m/N)_{m=1}^{N}$. The sequence of $r_{m}$ represents the probabilities corresponding to a coupon collector’s problem. It is therefore known that

[TABLE]

where $G$ is distributed as $\mbox{ Gumb}(-\gamma,1)$, as described in Section III.

On the complete graph, we can rewrite $p_{m}$ (the probability of going from $m$ infected nodes to $m+1$ on the next time step) as $p_{m}=r_{m}(1-\epsilon_{m})$, where $\epsilon_{m}=(m-1)/(N-1)$. So for $N\gg m$, the $p_{m}$’s resemble the $r_{m}$’s quite closely. Therefore both the front tail of $p^{(f)}$ and the back tail of $p^{(b)}$ look suspiciously like coupon collector’s processes. By this logic, we expect the total time to take over the complete graph should be just the sum of two coupon collector’s times. That is, we suspect that $T$ is the sum of two Gumbel random variables.

There are a few hangups with this intuitive argument:

$p^{(f)}$ and $p^{(b)}$ are about half the length of $r$. 2. 2.

$p_{m}$ doesn’t quite equal $r_{m}$ at small $m$.

Addressing the first hangup involves, once again, the front tails of these vectors. Each $\mathcal{E}(p_{m})$ has a standard deviation of $1/p_{m}$, which tells us that the smallest values of $p_{m}$ are the strongest drivers of the final distribution.

The fact that the events at low populations (of either infected or susceptible types) strongly determine most of the random fluctuations is something that has shown up in other evolutionary models, especially with selective sweeps Durrett and Schweinsberg (2004). So if we were to just truncate both $p^{(f)}$ and $r$ at some point, we should expect the limiting distributions of $S(p^{(f)})$ or $S(r)$ to not substantially change. We formalize this idea in Appendix  D, and find that it works out nicely.

We have a lot of options about where to truncate, but a useful truncation point is $B:=B(N)=\lfloor\sqrt{N-1}\rfloor$. The expression $\lfloor z\rfloor$ simply means we round $z$ down to the nearest integer. If we define

[TABLE]

then we have

[TABLE]

and

[TABLE]

Addressing the second hangup mostly involves formalizing $\epsilon_{m}$ as a rather small number. The details are outlined in Appendix E, where we find that

[TABLE]

From here we can daisy-chain the previous numbered equations in this section together and find that

[TABLE]

This means that we successfully piggy-backed on the result for star graphs to find that the resulting distribution for the complete graph is just a sum of two Gumbel random variables. The sum of two Gumbels has appeared previously in mathematically analogous places Aldous (2013); Van Der Hofstad et al. (2002). However our use of the coupon collector’s problem makes for a quick conceptual justification.

Figure 4 compares the takeover time distribution seen in simulations against the predicted distribution $\mathrm{Gumbel}(-\gamma,1)\star\mathrm{Gumbel}(-\gamma,1)$, and we see that this double-coupon logic works out well.

V d-dimensional lattice

As we did with the one-dimensional lattice, in our analysis of $d$-dimensional lattices we will assume periodic boundary conditions. The side length of the $d$-dimensional cube of $N$ nodes is denoted by $n=N^{1/d}$. We are also taking $1<d<\infty$, since we have already covered the 1D lattice and the infinite-dimensional lattice is a somewhat special case.

Unlike every previous case we have examined, we cannot consistently define $p_{m}$. The probability of infecting a new node will almost always depend on the specific location of all currently infected nodes. This means that all our previous approaches will not work well here. However, this does not bar us from making guesses based on reasonable approximations.

Although we could potentially get all kinds of weirdly-shaped clusters of infected nodes, that should not happen in expectation. Think back to the definition of our infection dynamics and Eq. (1). New infectees are added when a node on the boundary of the infected cluster gets randomly selected, and then one of its susceptible neighbors gets randomly selected and catches the infection.

Intuitively, it sounds like we have an expanding blob of infected nodes, with the expansion happening uniformly outward on every unit of surface area. This is a recipe for making sphere-like blobs in $d$ dimensions, at least at the start of the dynamics. As seen from the top half of Figure 5, this looks plausible in two dimensions.

The exact nature of this shape is actually a notoriously hard unsolved question. As we pointed out, there is a link between our infection model and first passage percolation on a lattice Auffinger et al. (2015). In that context, there is a rich literature surrounding questions about the nature of this cluster, but formal proofs of many of its properties have turned out to be difficult. However, convexity appears to be typical in the large size limit, and surface fluctuations should be relatively small Auffinger et al. (2015). Moreover, there is good reason to believe that on the two-dimensional (2D) lattice, the boundary of the expanding cluster is a one-dimensional curve, which will come in handy later Bramson and Griffeath (1980).

In any case, since the lattice is periodic, this infected cluster will keep expanding. This means that at the end of the dynamics, we should expect for the majority of susceptible nodes to also be in a single cluster, with insignificant enclaves elsewhere. This is borne out in simulations, as shown in the bottom half of Figure 5. If we focus on this majority susceptible cluster, we see that the end of the dynamics looks like a uniformly shrinking cluster of susceptible nodes, which is approximately the reverse of the uniformly growing infected cluster at the start. So, the beginning and end of the dynamics look similar once again, as they did for the complete graph.

More importantly, since this is a $d$-dimensional lattice, we can guess the surface area of these blobs. For a shape with a length-scale of $R$, we typically expect volume to scale as $R^{d}$ and surface area to go as $R^{d-1}$. So given an infected cluster of $m$ nodes, we expect it to have a surface area proportional to $m^{(d-1)/d}$. Assuming some uniformity, we should get that the typical probability of infecting a new node should be proportional to $m^{\eta}/N$ at the start of the dynamics, where the exponent $\eta$ is given by

[TABLE]

And just as in the case of the complete graph, this process at the start gets repeated backwards at the end.

This heuristic argument suggests that the total time to takeover should look like the sum of geometric variables $X(p_{m})$, where

[TABLE]

The fact that we only got a grip on $p_{m}$ up to a proportionality should not worry us. After all, that did not stop us when we worked through the star graph case earlier; back then we argued that such a proportionality constant would simply show up in the scaling factor in the denominator. If we treat this as a numerical problem, we do not need to explicitly find the scaling factor. Instead, we can examine $(T-\mu)/\sigma$, where $\mu$ and $\sigma$ are empirically obtained values for the average and standard deviation of $T$ respectively. Then any proportionality constants just get absorbed by the anonymous $\sigma$.

This reasoning further suggests that for $N$ sufficiently large,

[TABLE]

where $\sigma_{X}^{2}:=\sum_{m=1}^{N}(p_{m}^{-2}-p_{m}^{-1})$ is just the variance of the sum of geometric variables. But we already know how to approximate sums of geometric random variables. We can follow a similar procedure of truncation and perturbation as in the case of the complete graph. Assuming Eq. (18) is correct, we get

[TABLE]

where we define

[TABLE]

and

[TABLE]

for the sum of variances.

The truncation point $M$ is some increasing function of $N$, which can normally just be set to $M=N$. In the limit of large $N$, the distinction does not really matter. However, it seems frequently possible to tune $M$ to get a good fit on finite-$N$ cases, as the simulations of the 3D lattices in Figure 6 suggest.

In principle, we could try to use Eq. (28) to get a finite-$N$ estimate for this distribution. However, we do not expect any of these distributions to have a large-$N$ limit as easy as in the case of the star graph, nor for any of these distributions to have a name. For practical purposes, we can just simulate the right hand side of Eq. (20) directly, since generating and adding a large number of exponential variables is rather fast.

The Critical Dimension

Naively, we might expect the limiting distribution of $F^{\prime}$ to always be something between a Gumbel and a normal distribution. After all, $d=1$ implies $\eta=0$, which returns us to identical variables and the 1D ring, giving us the standard normal. Meanwhile, $d\to\infty$ implies $\eta\to 1$, which returns us to the coupon collector’s problem and the star graph, giving us the Gumbel. Incidentally, this argument suggests that the infinite-dimensional lattice has similar behavior as the complete graph under these dynamics. In between these extreme cases, we might expect the intermediate $d$’s to correspond to a family of intermediate distributions.

While this is generally true, there is a surprising caveat to be made about the case of $d=2$. Even though all the summands $(\mathcal{E}(m^{\eta})-1/m^{\eta})/\sqrt{H}$ are distinct, they start to resemble each other once $N$ gets sufficiently large.

For $d=2$, Eq. (16) gives $\eta=1/2$, which means that $H$ in Eq. (21) is the harmonic series. This $H$ diverges with $N$, giving each summand $(\mathcal{E}(m^{\eta})-1/m^{\eta})/\sqrt{H}$ a large denominator, and thus a small variance about a mean of zero. So, even though the summands are not identical random variables, they will become rather similar as we take $N$ to be large, suggesting that an improved version of the central limit theorem may apply. This intuition is confirmed by a careful analysis in Appendix F, showing that the Lindeberg-Feller theorem applies in this case.

Thus we predict a normal limiting distribution of $F^{\prime}$ in the specific case of the 2D lattice: as $N\rightarrow\infty$,

[TABLE]

for $d=2$. This prediction is borne out in simulation, as shown in Figure 7.

However, no dimension higher than $d=2$ can yield normally distributed takeover times. For each of $d=3,4,5,\ldots$, the distribution of $F^{\prime}$ will converge to a distinct limiting distribution between a normal and a Gumbel, as we initially suspected. The important distinction between $d=2$ and $d>2$ is that in the latter, $H$ always converges to a finite number. Because of that, $F^{\prime}$ will always have a nonzero third moment, preventing it from converging to a standard normal. For more details, see Appendix G.

VI Erdős-Rényi random graph

Unlike all the previous graphs we have seen, an Erdős-Rényi graph is randomly constructed. We start off with $N$ nodes, and add an edge between any two with some probability $0<\rho\leq 1$. In this section, we will condition on the graph being connected, so that complete takeover is always possible.

There is a good history of using generating functions to analyze desired properties on a random graph, including for various infection models Newman et al. (2001); Newman (2002); Newman and Watts (1999). But since we just finished analyzing the general lattice case, we can take another road.

Recall the central observation that let us recast $T$ as a sum of geometric random variables. That train of logic only really involved the graph having a well-defined dimension $d$. If we could define the dimension for other kinds of graphs, then all our observations from the previous section would simply carry over.

Imagine taking a cluster of $m$ nodes on an Erdős-Rényi graph. What is the surface area of said cluster? Well, in expectation, the $m$ nodes are externally connected to $O(pN)$ nodes, for $m\ll N$. So as $N$ gets large, the number of external neighbors in any cluster gets large as well, regardless of $m$. This is suggestive of an infinite-dimensional topology.

So, by collecting results from Eqs. (19) and (15), we can guess the limiting distribution of the takeover times $T$. Defining $\mu$ and $\sigma$ to be the empirical mean and standard deviation of $T$, we find

[TABLE]

where $G^{\prime}$ is a Gumbel random variable with a mean of zero and a variance of 1/2. One can check that the corresponding distribution for $G^{\prime}$ is $\mathrm{Gumbel}(-\gamma\sqrt{3}/\pi,\sqrt{3}/\pi)$.

We experimentally tested Eq. (23) by fixing a randomly generated Erdős-Rényi random graph, along with a seed at which the infection always started. Then we ran a million simulations of the stochastic infection process and compiled the observed distribution of takeover times. (The reason we fixed the graph beforehand was to avoid sampling multiple different values of $\mu$ and $\sigma$ over different realizations of the random graph.) The results of the experiment were consistent with our prediction, as shown in Figure 8.

VII Discussion

VII.1 Relation to other models

VII.1.1 Infection models

The model studied in this paper is intentionally simplified in several ways, compared to the most commonly studied models of infection. The purpose of the simplifications is to highlight how one aspect of the infection process – its network topology – affects the distribution of takeover times. However, the update rule also plays an important role. The assumptions we have made about it therefore deserve further comment.

Assumption 1: The infection is infinitely transmissible. When an infected node interacts with a susceptible node, the infection spreads with probability one. In a more realistic model, infection would be transmitted with a probability less than one.

Assumption 2: The infection lasts forever. Once infected, a node never goes back to being susceptible, or converts to an immune state, or gets removed from the network by dying. The dynamics of these more complicated models, known as SIS or SIR, have been studied on lattices and networks by many authors; for reviews, see Diekmann et al. (2012); Pastor-Satorras et al. (2015).

Assumption 3: The update rule is asynchronous. In other words, only one link is considered at a time. By contrast, in a model with synchronous updating, every link is considered simultaneously.

If the infection is further assumed to be infinitely transmissible, then at each time step every infected node passes the infection to every one of its susceptible neighbors. Such a infection, akin to the spreading of a flood or a wildfire, would behave even more simply than the process studied here. In fact, it would be too simple. The calculation of the network takeover time would reduce to a breadth-first search and its value would be bounded above by the network’s diameter. Note, however, that if the infection has a probability less than one of being transmitted to susceptible neighbors (such as in the original 1-type Richardson model Richardson (1973)), the system becomes nontrivial to analyze Auffinger et al. (2015); Richardson (1973).

VII.1.2 Models of evolutionary dynamics

More recently, the field of evolutionary dynamics Nowak (2006) has been extended to networks, and the field of evolutionary graph theory was born Lieberman et al. (2005). In general, the results in this field depend on modeling the spread of a mutant population using the Moran process Nowak (2006); Moran (1958). (Our model can be viewed as a limiting case of the Moran birth-death process, in the limit as the mutant fitness tends to infinity.) A number of important and interesting results have come from these studies of Moran dynamics, including the existence of network topologies that act as amplifiers of selection Adlam et al. (2015), increasing the probability of takeover, and also topologies that shift the takeover times we are considering Frean et al. (2013).

For example, working in the framework of evolutionary graph theory, Ashcroft, Traulsen, and Galla recently explored how network structure affects the distribution of “fixation times” for a population of $N$ individuals evolving by birth-death dynamics Ashcroft et al. (2015). The fixation time is defined as the time required for a fitter mutant (think of a precancerous cell in a tissue) to sweep through a population of less fit wild-type individuals (normal cells). Initially, a single mutant is introduced at a random node of the network. At each time step, one individual is randomly chosen to reproduce. With probability proportional to its fitness, it gives birth to one offspring, and one of its network neighbors is randomly chosen to die and be replaced by that offspring. The natural questions are: What is the probability that the lineage of the mutant will eventually take over the whole network? And if it does, how long does it take for this fixation to occur?

The calculations are difficult because there is no guarantee of mutant fixation (in contrast to our model, where the network is certain to become completely infected eventually). In the birth-death model, sometimes by chance a normal individual will be chosen to give birth, and its offspring will replace a neighboring mutant. If this happens often enough, the mutant population can go extinct and wild-type fixation will occur. Using Markov chains, Hindersin and colleagues provided exact calculations of the fixation probability and average fixation times for a wide family of graphs, as well as an investigation of the dependence on microscopic dynamics Hindersin and Traulsen (2015, 2014); Hindersin et al. (2016a, b). A challenge for this approach is that the size of the state space becomes intractable quickly: even with sparse matrix methods, it grows like $N2^{N}$ Hindersin et al. (2016a). For networks of size $N<23$, their computations showed that the distributions of mutant fixation times were skewed to the right, much like the Gumbels, sums of Gumbels, and intermediate distributions found analytically and discussed here in Sections III to VI.

VII.1.3 First-passage percolation

Our infection model is also closely related to first-passage percolation Aldous (2013); Auffinger et al. (2015). The premise behind this family of models can be described as follows: given a network, assign a random weight to each edge. By interpreting that weight as the time for an infection to be transmitted across that edge, and by choosing properly tuned geometric (or, more commonly, exponential) random variables as the edge weights, we can recreate our infection model.

Notice that percolation defines a random metric on the network, meaning that inter-node distances change from one realization to another. This leads to a number of natural questions. The most extensively studied is the “typical distance,” quantified by the total weight and number of edges on the shortest path between a pair of random nodes Bhamidi et al. (2013); Eckhoff et al. (2015a, b). It is also possible to analyze the “flooding time” Janson (1999); Van Der Hofstad et al. (2002), defined as the time to reach the last node from a given source node chosen at random. This quantity is the closest analogue, within first-passage percolation, of our takeover time. Indeed, a counterpart of our result for two Gumbels in the Erdős-Rényi random graph was obtained previously using these techniques Van Der Hofstad et al. (2002). However, we are unaware of flooding-time counterparts of our results about the takeover times for $d$-dimensional lattices.

Another natural question in first-passage percolation involves finding the long-time and large-$N$ limiting shape of the infected cluster. More precisely, given a fixed origin node, we can identify all nodes that can be reached from the the origin within a total path-weight of $t$ or less. This amounts to finding all the nodes that have been infected by the origin within time $t$, a problem that percolation theorists have typically studied in $d$-dimensional lattices. We saw an instance of such an expanding cluster in Section V. In a general number of dimensions, the provable nature of this shape may be complicated; the fluctuations of its boundary are thought to depend on the KPZ equations Auffinger et al. (2015); Kardar et al. (1986); Cieplak et al. (1996). The limiting shape is not typically a Euclidean ball, but it has been proven to be convex; see Aldous (2013); Auffinger et al. (2015) for an introductory discussion of these issues. In Figure 5, the nature of this cluster’s complement in a large torus was of concern to us, but that issue has not yet attracted mathematical attention, as far as we know.

VII.2 Applications to medicine: epidemic and disease incubation times and cancer mortality

For more than $100$ years, there have been intriguing empirical observations of “right-skewed” distributions in a remarkably wide range of phenomena related to disease Sawyer (1914); Miner (1922); Sartwell (1950, 1952, 1966). Examples include within-patient incubation periods for infectious diseases like typhoid fever Sawyer (1914); Miner (1922), polio Sartwell (1952), measles Goodall (1931) and acute respiratory viruses Lessler et al. (2009); exposure-based outbreaks like anthrax Brookmeyer et al. (2005) (see Nishiura (2007); Lessler et al. (2009) for more recent reviews); rates of cancer incidence after exposure to carcinogens Armenian and Lilienfeld (1974); and times from diagnosis to death for patients with various cancers Boag (1949) or leukemias Feinleib and MacMahon (1960).

The relationship between these phenomena and our model is intuitive: most of these processes depend on some sort of agent (a mutant cell, a virus, or a bacterium) invading and taking over a population, something which typically proceeds one “interaction” at a time. And as we have seen, our simple infection model automatically generates right-skewed distributions like Gumbels, sums of Gumbels, and intermediate distributions via a coupon-collection mechanism, for many kinds of population structures. So could it be that the right-skewed distributions so often seen clinically are, at bottom, a reflection of this same mathematical mechanism, a manifestation of an invasive, pathogenic agent spreading through a network of cells or people?

To test the plausibility of this idea, we need to amend our model slightly. Until now we have focused exclusively on the time $T$ to total takeover of a network. But in most scenarios related to disease, total takeover is not the relevant consideration. Sufficient takeover is what matters. For example, a patient need not have every single one of their bone marrow stem cells replaced by leukemic cells before they die from leukemia. Death presumably occurs as soon as some critical threshold is crossed – which is probably the case for diseases with infectious etiologies as well. So let us now check whether changing the criterion from total takeover to partial takeover changes our results, or not.

Times to partial takeover: truncation

Define $T_{\theta}$ to be the time for $\lfloor\theta N\rfloor$ out of $N$ members to be infected, with the interesting range of $\theta$’s being $0.5\leq\theta<1.0$. For the sake of example, consider the complete graph as our network topology, so we have $p_{m}=(m/N)(1-(m-1)/(N-1))$.

As in the analysis for the complete takeover times, we can split $T_{\theta}$ into a front and back part $T^{(f)}$ and $T^{(b)}$, with the front covering up to about $N/2$ and the back covering the remaining. Then

[TABLE]

where $\mu_{(f)}$ is the mean of $T^{(f)}$, and $\mu_{(b)}$ and $\sigma_{(b)}$ are the mean and standard deviation of $T^{(b)}$. However, it is easy to show that $\sigma_{(b)}^{2}/N^{2}$ converges to 0 as $N$ gets large, regardless of $\theta$. So we expect the distribution of $T_{\theta}$ in this case to asymptotically approach a Gumbel. As seen in Figure 9, similar results hold for Erdős-Rényi random graphs, even for $\theta=0.90$.

Thus, for complete graphs and Erdős-Rényi random graphs, the right-skewed distributions for complete takeover persist when we relax the criterion to partial takeover. In that respect our results seem to be robust.

The resilience of the Gumbel distribution is important to appreciate. As pointed out by Read Read (1998), a Gumbel distribution can be impersonated by a properly tuned 3-parameter lognormal distribution; see Appendix H for further details. A 3-parameter lognormal distribution has a density function

[TABLE]

provided $x>c$.

It is this 3-parameter lognormal distribution that has been frequently noted in empirical studies of disease incubation times. Originally proposed and elaborated by Sartwell Sartwell (1950, 1952, 1966) as a curve-fitting model, its seeming generality has led to it being called “Sartwell’s Law.” But it has always lacked a theoretical underpinning. Even recent reviews consider the origin of lognormal incubation times to be unresolved Nishiura (2007). In contrast, Gumbel and related distributions arise very naturally from the model studied here and from other infection models Williams (1965); Gautreau et al. (2007), and may provide a more suitable theoretical foundation than lognormals in that sense.

he spreading dynamics of processes on realistic, spatial topologies, Lieberman et al. introduced the field of ‘evolutionary graph theory’.

Appendix A The Lindeberg condition (1D)

When we were analyzing the 1D lattice, we could not quite make use of the classical central limit theorem, since our variables $X$ have a dependence on $N$, being Geo($1/N$). Fortunately, there is the Lindenberg-Feller variant of the central limit theorem, which lets us get the desired normal convergence. See reference Durrett (1991), page 98 for more details. However, we must first satisfy some special conditions before we can cite it.

For convenience’s sake, let us define

[TABLE]

So $E[Y]=0$, and $\sum_{m=1}^{N-1}E[Y^{2}]=1$, which satisfies two of the three conditions. However, there is still the matter of the titular Lindeberg condition on the restricted second moments, which says for any fixed $\epsilon>0$,

[TABLE]

To verify this, first notice that $|Y|>\epsilon$ means that we need either $X>N+\epsilon(N-1)\sqrt{N}$ or $X<N-\epsilon(N-1)\sqrt{N}$. However, the minus case will not come up in the limit; as $N$ gets large it would require $X$ to be negative, which is not possible.

Letting $c:=\lfloor N+\sqrt{N}(N-1)\epsilon\rfloor$, we get

[TABLE]

And so,

[TABLE]

Here, the first limit looks like

[TABLE]

Meanwhile the second limit can be bounded above (with some constant $C$) by

[TABLE]

So the total sum of conditional expectations converges to 0 as $N$ gets large, and so the Lindeberg condition is satisfied. This allows us to cite the theorem, and confirms that in the limit we get

[TABLE]

Appendix B Geometric variables converging to exponential variables

**Proposition: ** Say we have a positive sequence $(p_{m})_{m=1}^{M}$, and some function $L:=L(M)$ such that $\lim_{M\to\infty}L=\infty$ and

[TABLE]

Then if $T:=\sum_{m=1}^{M}X(p_{m})$, $F:=\sum_{m=1}^{M}\mathcal{E}(p_{m})$, and $\mu:=\sum_{m=1}^{M}1/p_{m}$, we have

[TABLE]

Proof: This is proven by finding the characteristic functions for both sides, and showing that the ratio of these functions goes to 1 as $M$ gets large. The characteristic function of a random variable uniquely determines its distribution, so this is a rather powerful statement.

Let us define

[TABLE]

If we split $T$ into the sum of geometric random variables and rearrange, we eventually get

[TABLE]

Similarly, if we set

[TABLE]

then after we proceed through some more algebra, we find that

[TABLE]

Let us fix $t$ so that we can pointwise consider the ratio of the characteristic functions. After some manipulation, we find

[TABLE]

We assumed that $L$ gets large, so there is some function $R_{1}:=R_{1}(M)$ that has vanishing magnitude with large $M$ such that

[TABLE]

So then we have

[TABLE]

Notice that $|1-it/(p_{m}L)|\geq 1$. In addition, we already know the sum of $1/(p_{m}L^{2})$ goes to 0, so it must be that each individual $p_{m}L^{2}$ gets large for all $m$. This ensures the second term is small, and therefore it can be rewritten exactly as an appropriate exponential.

So

[TABLE]

where the final limit comes from our assumption on $\sum_{m=1}^{M}\frac{1}{p_{m}L^{2}}$. The limit converges to 1, which establishes the proposition.

Appendix C Sum of exponentials

**Proposition: ** If we have exponential random variables $\mathcal{E}(p_{m})$ for $m=1,\ldots,n$, with $p_{m}$ distinct, then $\sum_{m=1}^{n}\mathcal{E}(p_{m})$ is distributed according to the density

[TABLE]

*on $x\geq 0$. *

Proof: This is a straightforward induction for the most part. The base case is simply checked by plugging in $n=1$. To get the inductive step down, we just convolve the previous step with a new exponential distribution, so

[TABLE]

After calculating for a bit, we find

[TABLE]

The first term is in the desired form, but the second term requires some work. After some further manipulation, we can get

[TABLE]

where we define

[TABLE]

We can interpret $b(z)$ as a polynomial of at most degree $n-1$ in $z$ (a Lagrange polynomial, to be specific).

But notice that for $l=1,\ldots,n$, then $b(p_{l})=1$. This means that $b(z)-1$ is a polynomial with n distinct roots, which is more than what its maximum degree should normally allow. The only way that is possible is if $b(z)-1$ is a constant 0, so $b(z)\equiv 1$. Plugging this in and simplifying gives

[TABLE]

which is the desired result.

Appendix D Truncation of sequences

**Proposition: ** Let $L:=L(M)$, $B:=B(M)$ with $\lim_{M\to\infty}L(M)=\lim_{M\to\infty}B(M)=\infty$. Further say that $B$ is integer valued with $1\leq B\leq M$. Given a positive sequence $p(M)=(p_{m})_{m=1}^{M}$, assume

[TABLE]

and $\lim_{M\to\infty}1/[p_{k(M)}L]=0$ given that $M\geq k(M)>B$. Then

[TABLE]

Proof: As before, the proof involves showing the ratio of characteristic functions converges to 1. The full series on the right has the function

[TABLE]

and that the truncated series on the left has

[TABLE]

So naturally, we fix a $t$ and get the ratio

[TABLE]

Because of the our last condition, we know that $1/(p_{m}L)$ is small for all $m$ in this range. So we can do a Taylor expansion and make a function $R_{1}$ which is small in magnitude so that

[TABLE]

Again, $|1-it/(p_{m}L)|>1$ and $p_{m}L$ is large, so we can again shift to an exponential to get

[TABLE]

where $R_{2}$ is small in magnitude again. But notice that this is based on the tail of a convergent sum. So

[TABLE]

And so $\lim_{M\to\infty}\phi/\hat{\phi}=1$.

Appendix E Edge perturbations

In principle we could show a more general statement here, but we are only going to directly calculate the effect of a perturbation once in this paper. So, for the sake of readability, we are just going to do this specific example.

Recall that for the complete graph $Np_{m}=m(1-(m-1)/(N-1))=Nr_{m}(1-\epsilon_{m})$ with $r_{m}=m/N$. Also recall that $p^{(T)}$ and $r^{(T)}$ are the truncated sequences up to $B=\lfloor\sqrt{N-1}\rfloor$. Let $\hat{\phi}$ be the characteristic function associated with the normalized sum $S(p^{(T)})$, and that $\phi$ is the characteristic function associated with $S^{(T)}$. Then

[TABLE]

Notice $\epsilon_{m}=(m-1)/(N-1)=O\left(N^{-1/2}\right)$ in the range of $m$’s in the product, and is therefore small. We go through some Taylor expansions and cancellations, and using $R_{j}$ to represent functions of small magnitude, we find

[TABLE]

The first term can be once again turned into an exponential (thanks to the smallness of $\epsilon$), and so we get

[TABLE]

Therefore, we get convergence to 1 if the sum converges to 0 as $N$ gets large. But this sum is easy to bound from above. That is,

[TABLE]

This means we get that $\phi/\hat{\phi}\to 1$, implying that the truncated sum for the complete graph converges to the truncated distribution for the coupon collector’s problem.

Appendix F The Lindeberg condition (2D)

Much like in the 1D lattice case, we are unable to directly use the typical central limit theorem, because the variables are not identical and have a dependence on $N$. But once again, we can apply the Lindenberg-Feller theorem. We are going to focus on the 2D case, so we have $\eta=1-1/d=1/2$. Let

[TABLE]

Because we are only looking at the special 2D case, $H=H(N,2(1/2))=\sum_{k=1}^{N}1/k$.

Notice for any higher dimension $d>2$, then we would get $H=\sum_{k=1}^{N}1/k^{2-2/d}$, which quickly converges to a finite number as $N$ gets large, whereas with $d=2$, we have $H$ gets large for large $N$. This distinction is what lets us apply the theorem to the 2D case, but not the rest.

Anyway, it is easy to check that $E[Y_{N,m}]=0$ and $E[Y_{N,m}^{2}]=m^{-2\eta}/H$. So then

[TABLE]

So in order to apply the theorem, we only need to check if for any fixed $\epsilon>0$, we have

[TABLE]

If this final condition holds, then we can cite the theorem and conclude that $\sum_{m=1}^{N}Y_{N,m}$ is distributed as a normal as $N$ gets large.

We need not care about the $Y_{N,m}<-\epsilon$ case, because this is equivalent to asking for $\mathcal{E}(\sqrt{m})<1/\sqrt{m}-\epsilon\sqrt{H}$. However, $H$ scales as $\log(N)$ in the limit of large $N$ whereas $1/\sqrt{m}\leq 1$, so this quantity will always eventually become negative, whereas exponential variables are always positive.

Therefore, let us focus on the positive half. Letting $c:=1/\sqrt{m}+\epsilon\sqrt{H}$ and integrating, we get

[TABLE]

Substituting in gives us

[TABLE]

That last term will be the dominant term as $N$ gets large, so we can choose some positive constant $C_{1}$ (which may depend on $\epsilon$) such that

[TABLE]

So, we have

[TABLE]

We can bound the harmonic sum $H$ from below with a constant times $\log(N)$, so there is some positive $C_{2}$ such that

[TABLE]

We can approximate this sum from above by interpreting it as a Riemann sum. By taking the appropriate integral, we get

[TABLE]

Second moments are always nonnegative, which means

[TABLE]

So Eq. (30) is finally confirmed. As a consequence, we can finally cite the Lindeberg-Feller theorem, and know that

[TABLE]

Appendix G Non-Normality of $d>2$

There are a lot of possible ways to show a distribution does not converge to a normal in a limit. But to show that the distribution of $F^{\prime}$ for a $d>2$ dimensional lattice (as defined in Eq. (20)) is not normal, it will suffice to consider the moments. We already know that $F^{\prime}$ has a mean of 0 and a variance of 1; so if $F^{\prime}$ went like a normal, then we should expect that $\lim_{N\to\infty}E[F^{\prime 3}]=0$ by symmetry.

We can reuse Eq. (27) to find the characteristic function $\phi$ of $F^{\prime}$ by plugging in $p_{m}=m^{\eta}$ and $L^{2}=H(2\eta)=\sum_{m=1}^{N}m^{-2\eta}$. Because $d>2$, then $2\eta>1$. By the definition of the characteristic function, we know that if we expand $\phi$ in powers of $t$, then

[TABLE]

So we can get the third moment by just reading off the coefficient of the $t^{3}$ term.

Returning to equations (20) and (27), let $x_{m}=itm^{-\eta}H(2\eta)^{-1/2}$. Using the standard expansions for $e^{x}$ and $1/(1-x)$, we find

[TABLE]

If we do not care about high order terms in $t$, then this is an easy product to take. In fact, if we collect terms and plug in for $x_{m}$, we get

[TABLE]

This means for any finite $N$, the third moment of $F^{\prime}$ is simply

[TABLE]

Although $H$ is a function that depends on $N$, this quantity will never get large. In fact, since $1\geq\eta>1/2$, then we know in the limit of large N that

[TABLE]

where $\zeta$ is the Riemann zeta function. In the range of $\eta$’s presented, $\zeta$ neither diverges nor hits zero, so the above will never be zero. In fact, the right hand side of Eq. (32) is monotone, so each distinct $1\geq\eta>1/2$ will produce a distinct third moment, and therefore a distinct distribution. As a side note, if we take $\eta\to 1$, we get $-12\sqrt{6}\zeta(3)/\pi^{3}$, which is the correct value for the third moment of a normalized Gumbel distribution, as expected. However, since the rest are distinct, that means we only transition to an exact Gumbel in the extreme limit.

In summary: given $d>2$, we never expect $F^{\prime}$ to have a zero third moment in the limit of large $N$, and so $F^{\prime}$ can never converge to a normal distribution. Moreover, because their third moments depend on $\eta=1-1/d$, we expect $F^{\prime}$ to converge to a different distribution for each $d>2$. Hence we expect there is no upper critical dimension.

Appendix H Lognormal distributions can masquerade as Gumbels

In Section VII.2 we noted that the distribution of disease incubation periods and other times of medical interest have often been fit by a lognormal. However, given the noise in real data, it is entirely possible that the true distribution should have been be a Gumbel (or a sum of two Gumbels), and was impersonated by a similar-looking lognormal.

Moreover, since most studies used 3-parameter lognormals, it would always be possible to match the first three moments of the data. We do as such in Figure 10, producing a very close fit to both a Gumbel and a convolution of two Gumbels. We can compare these densities using the Kolmogorov metric, given by the maximum difference between their cumulative distribution functions. Using this, we find that the normalized Gumbel is $\approx 0.0034$ away from its corresponding lognormal in this metric. For the sum of two Gumbels, we can numerically estimate that its corresponding lognormal is $\lesssim 10^{-2}$ away in the Kolmogorov metric.

Acknowledgments

Thanks to David Aldous, Rick Durrett, Remco van der Hofstad, Lionel Levine, and Piet Van Mieghem for helpful conversations. This research was supported by a Sloan Fellowship and NSF Graduate Research Fellowship grant DGE-1650441 to Bertrand Ottino-Löffler in the Center for Applied Mathematics at Cornell, and by NSF grants DMS-1513179 and CCF-1522054 to Steven Strogatz. Jacob Scott acknowledges the NIH for their generous loan repayment grant.