Takeover times for a simple model of network infection

# Takeover times for a simple model of network infection

Bertrand Ottino-Löffler Center for Applied Mathematics, Cornell University, Ithaca, New York 14853    Jacob G. Scott Cleveland Clinic, Departments of Translational Hematology and Oncology Research and Radiation Oncology    Steven H. Strogatz Center for Applied Mathematics, Cornell University, Ithaca, New York 14853
August 21, 2019
###### 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 does it take to completely infect a network of nodes, starting from a single infected node? An analogy to the classic “coupon collector” problem of probability theory reveals that the takeover time 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 , the takeover time is distributed as a Gumbel for the star graph; as the sum of two Gumbels for a complete graph and an Erdős-Ré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 -dimensional lattices with (these distributions approach the sum of two Gumbels as 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.

## 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 , how long does it take for the entire network to become infected? We call this the takeover time . 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 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 , 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 , the takeover time 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 -dimensional cubic lattices, the dependence on is intriguing: we find that is normally distributed for and , then becomes skewed for and approaches the sum of two Gumbels as 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 nodes arranged into a ring.

Suppose that nodes are currently infected. Let denote the probability that a susceptible node gets infected in the next time step. Notice that for a more complicated graph, might not be a well-defined concept, because it could depend on more than 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 , 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 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

 pm= (Prob. of selecting node on boundary) ×(Prob. of selecting susceptible neighbor). (1)

Hence, for the ring, the probability that the infection spreads on the next time step reduces to for all .

Next, define the random variable as the number of time steps during which the network has exactly infected nodes. The probability that this state lasts for time steps is then given by

 P(Xm=k)=qk−1mpm,

for where . To see this, note that is the probability that no new infection occurs on the first steps, times the probability that infection does occur on step .

Thus, for any network where is well defined, the time spent with infected nodes is a geometric random variable, with mean and variance . In particular, since the ring has for all , we find that has mean and variance in this case.

The takeover time for any network is

 T=N−1∑m=1Xm,

the sum of all the individual times required to go from to infected nodes, for . (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 are identical. However, their means and variances depend on , which prevents us from invoking the usual Central Limit Theorem to deduce the limiting distribution of . However, we can invoke a generalization of it known as the Lindeberg-Feller theorem. See Appendix A for more details.

After normalizing by its mean, , and its standard deviation, , we find

 T−N(N−1)(N−1)√Nd→Normal(0,1), (2)

where the symbol means convergence in distribution as 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, 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 be the number of spoke nodes that are currently infected. As in the ring case, (the probability to go from to in the next time step) is a well-defined quantity that depends on alone, and not on any other details of the network state. Using the logic of Eq. (1), we get

 pm=1N+1⋅N−mN (3)

for . Here, is the probability of choosing the infected hub as the first node, and is the probability of selecting one of the currently susceptible spoke nodes, out of the spoke nodes in total, as its neighbor.

Now that is in hand for the star graph, we can define the random variable and the takeover time just as we did for the one-dimensional ring. The only difference is that the -dependence of is now controlled entirely by the factor .

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 distinct items, and you have of them so far. If you are given one of the items at random (with replacement), the probability it is new to your collection is , the same factor we saw above, and precisely analogous to the probability of adding a new node to the infected set. Likewise, the waiting time to collect all items is precisely analogous to the time needed to take over the whole star graph. The only difference is the constant factor 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 , 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 by exponential random variables , with density

 P(E(pm)=x)dx=pme−pmxdx;x≥0.

From here we define the random variable , which has mean .

It can be shown (see Appendix B) that for a large class of and normalizing factors that

 T−μL∼F−μL, (4)

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

In the traditional coupon collector’s problem we would take ; but because of that factor, what we want is . Thanks to the fact we are now using exponential variables, we now know

 F/L=N−1∑m=0E(pm/L)=N∑k=1E(k)

(using ). 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 we can simply write the relevant distribution function for as

 gN(x)=N∑k=1ke−kxN∏r≠krr−k,

which can be manipulated into

 gN(x)=Ne−x(1−e−x)N−1.

From here, we can find the distribution for , and by extension. Therefore, taking the limit of large and using the standard approximation of the harmonic sum for gives us

 f(x)=e−(x+γ)exp(−e−(x+γ)) (5)

as the density, where is the Euler-Mascheroni constant. The density in (5) is a special case of the Gumbel distribution, denoted and defined to have the density

 h(x)=β−1e−(x−α)/βexp(−e−(x−α)/β). (6)

Specifically, we find

 T−μN(N+1)d→G, (7)

where is a Gumbel random variable distributed according to .

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 nodes corresponds to a “well-mixed population” and is one of the most common topologies in infection models. This network consists of mutually connected nodes, so the location of the initial infection does not matter.

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

 pm=mN(1−m−1N−1)=mN⋅N−mN−1 (8)

for . For the sake of convenience, we will collect these probabilities into a vector . As in the case of the star graph, we can approximate the takeover time by summing exponential random variables instead of geometric ones. So

 T−μN∼N−1∑m=1E(pm)−1/pmN=:S(p). (9)

To compress notation, we defined to be the normalized sum of exponential random variables across the entries of the vector .

The specific in Eq. (8) has some helpful symmetry. Notice that if , then

 pk=kN⋅N−kN−1=N−mN⋅mN−1=pm.

This symmetry means that the second half of the takeover looks just like the first half played backwards. If we set to be the front half of the -vector and to be the back half of , then we know

 S(p)=S(p(f))+S(p(b)). (10)

Because we have a symmetry and the order we add the individual exponential variables will not matter, the random variables and should be equal in distribution. Although 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 and to , where . The sequence of represents the probabilities corresponding to a coupon collector’s problem. It is therefore known that

 S(r)d→G (11)

where is distributed as , as described in Section III.

On the complete graph, we can rewrite (the probability of going from infected nodes to on the next time step) as , where . So for , the ’s resemble the ’s quite closely. Therefore both the front tail of and the back tail of 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 is the sum of two Gumbel random variables.

There are a few hangups with this intuitive argument:

1. and are about half the length of .

2. doesn’t quite equal at small .

Addressing the first hangup involves, once again, the front tails of these vectors. Each has a standard deviation of , which tells us that the smallest values of 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 and at some point, we should expect the limiting distributions of or 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 . The expression simply means we round down to the nearest integer. If we define

then we have

 S(p(f))∼S(p(T)),S(p(b))∼S(p(T)) (12)

and

 S(r)∼S(r(T)). (13)

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

 S(r(T))∼S(p(T)). (14)

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

 T−μNd→G+G. (15)

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 , 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 -dimensional lattices we will assume periodic boundary conditions. The side length of the -dimensional cube of nodes is denoted by . We are also taking , 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 . 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 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 -dimensional lattice, we can guess the surface area of these blobs. For a shape with a length-scale of , we typically expect volume to scale as and surface area to go as . So given an infected cluster of nodes, we expect it to have a surface area proportional to . Assuming some uniformity, we should get that the typical probability of infecting a new node should be proportional to at the start of the dynamics, where the exponent is given by

 η=d−1d. (16)

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 , where

 pm≈mηN(1−mηN). (17)

The fact that we only got a grip on 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 , where and are empirically obtained values for the average and standard deviation of respectively. Then any proportionality constants just get absorbed by the anonymous .

This reasoning further suggests that for sufficiently large,

 T−μσ∼N∑m=1X(pm)−1/pmσX, (18)

where 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

 T−μσ∼1√2(F′+F′) (19)

where we define

 F′:=M∑m=1E(mη)−1/mη√H (20)

and

 H:=H(2η;M)=M∑m=11/m2η (21)

for the sum of variances.

The truncation point is some increasing function of , which can normally just be set to . In the limit of large , the distinction does not really matter. However, it seems frequently possible to tune to get a good fit on finite- 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- estimate for this distribution. However, we do not expect any of these distributions to have a large- 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 to always be something between a Gumbel and a normal distribution. After all, implies , which returns us to identical variables and the 1D ring, giving us the standard normal. Meanwhile, implies , 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 ’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 . Even though all the summands are distinct, they start to resemble each other once gets sufficiently large.

For , Eq. (16) gives , which means that in Eq. (21) is the harmonic series. This diverges with , giving each summand 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 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 in the specific case of the 2D lattice: as ,

 T−μσd→Normal(0,1) (22)

for . This prediction is borne out in simulation, as shown in Figure 7.

However, no dimension higher than can yield normally distributed takeover times. For each of , the distribution of will converge to a distinct limiting distribution between a normal and a Gumbel, as we initially suspected. The important distinction between and is that in the latter, always converges to a finite number. Because of that, 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 nodes, and add an edge between any two with some probability . 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 as a sum of geometric random variables. That train of logic only really involved the graph having a well-defined dimension . 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 nodes on an Erdős-Rényi graph. What is the surface area of said cluster? Well, in expectation, the nodes are externally connected to nodes, for . So as gets large, the number of external neighbors in any cluster gets large as well, regardless of . 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 . Defining and to be the empirical mean and standard deviation of , we find

 T−μσ∼G′+G′, (23)

where is a Gumbel random variable with a mean of zero and a variance of 1/2. One can check that the corresponding distribution for is .

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 and 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 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  Hindersin et al. (2016a). For networks of size , 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 -dimensional lattices.

Another natural question in first-passage percolation involves finding the long-time and large- 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 or less. This amounts to finding all the nodes that have been infected by the origin within time , a problem that percolation theorists have typically studied in -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 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 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 to be the time for out of members to be infected, with the interesting range of ’s being . For the sake of example, consider the complete graph as our network topology, so we have .

As in the analysis for the complete takeover times, we can split into a front and back part and , with the front covering up to about and the back covering the remaining. Then

 Tθ−μθN=T(f)−μ(f)N+σ(b)NT(b)−μ(b)σ(b),

where is the mean of , and and are the mean and standard deviation of . However, it is easy to show that converges to 0 as gets large, regardless of . So we expect the distribution of 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 .

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

provided .

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 have a dependence on , being Geo(). 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

 Y:=X−N(N−1)√N.

So , and , 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 ,

 limN→∞N−1∑m=1E[Y2;|Y|>ϵ]=0.

To verify this, first notice that means that we need either or . However, the minus case will not come up in the limit; as gets large it would require to be negative, which is not possible.

Letting , we get

 E[Y2;Y>ϵ] =∞∑k=c+1(k−N(N−1)√N)21N(1−1N)k−1 =1(N−1)N2(1−1N)cN[c2+N(N−1)].

And so,

 limN→∞N−1∑m=1E[Y2;|Y|>ϵ] =limN→∞1N(N−1)(1−1N)c[c2+N(N−1)] =limN→∞(1−1N)c+limN→∞c2N(N−1)(1−1N)c.

Here, the first limit looks like

 limN→∞(1−1N)N−ϵN1/2+ϵN3/2=0.

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

 limN→∞CN(1−1N)ϵN3/2=0.

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

 N−1∑m=1X−N(N−1)√Nd→Normal(0,1). (24)

## Appendix B Geometric variables converging to exponential variables

Proposition: Say we have a positive sequence , and some function such that and

 limM→∞M∑m=11pmL2=0.

Then if , , and , we have

 T−μL∼F−μL. (25)

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

Let us define

 Φ:=E[exp(itT−μL)].

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

 Φ=M∏m=1pmexp[(it/L)(1−1/pm)]1−qmexp(it/L). (26)

Similarly, if we set

 ϕ:=E[exp(itF−μL)],

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

 ϕ=M∏m=1exp[−it/(pmL)]1−it/(pmL). (27)

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

 ϕ/Φ=M∏m=1exp(−it/L)−qmpm[1−it/(pmL)].

We assumed that gets large, so there is some function that has vanishing magnitude with large such that

 exp(−it/L)=1+(−it/L)+R1t2/L2.

So then we have

 ϕ/Φ=M∏m=1(1+t2pmL2R11−it/(pmL)).

Notice that . In addition, we already know the sum of goes to 0, so it must be that each individual gets large for all . This ensures the second term is small, and therefore it can be rewritten exactly as an appropriate exponential.

So

 ϕ/Φ= M∏m=1exp[R2t2/(pmL)] = exp[t2R2M∑m=11pmL2]→1,

where the final limit comes from our assumption on . The limit converges to 1, which establishes the proposition.

## Appendix C Sum of exponentials

Proposition: If we have exponential random variables for , with distinct, then is distributed according to the density

 gn(x)=n∑k=1pke−pkxn∏r=1,r≠kprpr−pk (28)

on .

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

 gn+1(x)=∫x0pn+1e−pn+1(x−y)gn(y)dy.

After calculating for a bit, we find

 gn+1(x)= n∑k=1pke−pkxn+1∏r≠kprpr−pk +n∑k=1pkpn+1pk−pn+1e−pn+1xn∏r≠kprpr−pk.

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

 Second Term= e−pn+1x(n∏k=1pkpk−pn+1)b(pn+1),

where we define

 b(z):=n∑k=1n∏r≠kpr−zpr−pk.

We can interpret as a polynomial of at most degree in (a Lagrange polynomial, to be specific).

But notice that for , then . This means that 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 is a constant 0, so . Plugging this in and simplifying gives

 gn+1(x)=n+1∑k=1pke−pkxn+1∏r=1,r≠kprpr−pk,

which is the desired result.

## Appendix D Truncation of sequences

Proposition: Let , with . Further say that is integer valued with . Given a positive sequence , assume

 limM→∞M∑m=11(pmL)2=A<∞,

and given that . Then

 B∑m=1E(pm)−1/pmL∼M∑m=1E(pm)−1/pmL. (29)

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

 ϕ=M∏m=1exp[−it/(pmL)]1−it/(pmL),

and that the truncated series on the left has

 ^ϕ=B∏m=1exp[−it/(pmL)]1−it/(pmL).

So naturally, we fix a and get the ratio

 ϕ/^ϕ=M∏m=B+1exp[−it/(pmL)]1−it/(pmL).

Because of the our last condition, we know that is small for all in this range. So we can do a Taylor expansion and make a function which is small in magnitude so that

 ϕ/^ϕ=M∏m=B+1(1+t2(pmL)2R11−it/(pmL)).

Again, and is large, so we can again shift to an exponential to get

 ϕ/^ϕ=exp[M∑m=B+1R2t2(pmL)2],

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

 limM→∞M∑m=B+11(pmL)2 =limM→∞(M∑m=11(pmL)2−B∑m=11(pmL)2) =A−A=0.

And so .

## 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 with . Also recall that and are the truncated sequences up to . Let be the characteristic function associated with the normalized sum , and that is the characteristic function associated with . Then

 ϕ/^ϕ= B∏m=1exp[−itrmN+itpmN]1−it/(pmN)1−it/(rmN) = B∏m=1exp[−itm(1−11−ϵm)] ×1−it(1−ϵm)−1/m1−it/m.

Notice in the range of ’s in the product, and is therefore small. We go through some Taylor expansions and cancellations, and using to represent functions of small magnitude, we find

 ϕ/^ϕ=B∏m=1(1+itR2ϵmm)exp(itR1ϵmm).

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

 ϕ/^ϕ=exp(2itR3B∑m=1ϵmm).

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

 B∑m=1ϵmm= ⌊√N−1⌋∑m=11mm−1N−1 ≤ ⌊√N−1⌋∑m=11mmN−1 ≤(N−1)−1/2N→∞−−−−→0.

This means we get that , 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 . But once again, we can apply the Lindenberg-Feller theorem. We are going to focus on the 2D case, so we have . Let

 YN,m=E(√m)−1/√m√H.

Because we are only looking at the special 2D case, .

Notice for any higher dimension , then we would get , which quickly converges to a finite number as gets large, whereas with , we have gets large for large . This distinction is what lets us apply the theorem to the 2D case, but not the rest.

Anyway, it is easy to check that and . So then

 N∑m=1E[Y2N,m]=1HN∑m=11m=1.

So in order to apply the theorem, we only need to check if for any fixed , we have

 Lind.:=limN→∞N∑m=1E[Y2N,m;|YN,m|>ϵ]?=0. (30)

If this final condition holds, then we can cite the theorem and conclude that is distributed as a normal as gets large.

We need not care about the case, because this is equivalent to asking for . However, scales as in the limit of large whereas , so this quantity will always eventually become negative, whereas exponential variables are always positive.

Therefore, let us focus on the positive half. Letting and integrating, we get

 E[Y2N,m;YN,m>ϵ] =∫∞c√me−√mx(x−1/√m√H)2dx =1Hm(1+(√mc)2)e−√mc.

Substituting in gives us

 E[Y2N,m;YN,m>ϵ] =e−1(2Hm+2ϵ√m√H+ϵ2)exp(−ϵ√m√H).

That last term will be the dominant term as gets large, so we can choose some positive constant (which may depend on ) such that

 E[Y2N,m;YN,m>ϵ]≤C1exp(−ϵ√m√H).

So, we have

 Lind. =limN→∞N∑m=1E[Y2N,m;YN,m>ϵ] ≤limN→∞N∑m=1C1exp(−ϵ√m√H).

We can bound the harmonic sum from below with a constant times , so there is some positive such that

 Lind.≤limN→∞N∑m=1C1exp(−C2√log(N)√m).

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

 ∫N0C1exp(−C2√log(N)√x)dx = 2C1C2log(N)[1−exp(−C2√Nlog(N)) −C2√Nlog(N)exp(−C2√Nlog(N))] ≤ 2C1C2log(N).

Second moments are always nonnegative, which means

 0≤Lind.≤