Selective sweep probabilities in spatially expanding populations

Domain: neuroscience
Published: 2026-02-11
Editorial outcome: major_revision

Abstract
Evolution during range expansions shapes biological systems from microbial communities and tumours to invasive species. A fundamental question is whether, when a beneficial mutation arises during a range expansion, it will evade clonal interference and sweep to fixation. However, most theoretical investigations of range expansions have considered regimes in which selective sweeps are effectively impossible, while studies of selective sweeps have assumed constant population size or ignored spatial structure. Here we use mathematical modelling and analysis to investigate selective sweep probabilities and timings in biologically relevant scenarios, including the case in which mutants can displace a slowly spreading wild type. Assuming constant expansion speed, we find surprisingly simple approximate and exact expressions for sweep probabilities in one, two and three dimensions, which are independent of mutation rate. Agent-based simulations confirm that our predictions are accurate for the spatial Moran process and remain informative when mutation effects on fitness are random and multiplicative. We further compare and synthesise our results with those obtained for alternative growth laws. Parameterised for human tumours, our model predicts that selective sweeps are rare except during early solid tumour growth, thus providing a general, pan-cancer explanation for findings from recent sequencing studies.

Introduction
Range expansion—the spatial spread of populations into new regions—is ubiquitous across biological scales and alters the course of evolution in distinct, often profound ways that remain incompletely understood1. Among cell populations, evolution during range expansions determines the development and spatial heterogeneity of biofilms2, tumours3, mosaicism4 and normal tissue5. At the species level, range expansions influenced human evolution6 and are of growing importance as climate change forces organisms into new habitats7,8. Prior theoretical and experimental investigations of evolution during range expansion have considered the case in which the wild type population spreads into new territory much faster than any mutant can displace the wild type9–12. In this scenario, which is typical of microbial colonies growing in vitro, mutants essentially expand only at the population boundary and selective sweeps are precluded. The alternative case of slow range expansion and strong selection has been unexplored but is widely plausible. Consider, for example, an invasive species that is rapidly adapting to its new environment while gradually displacing a resident competitor13, or a bacterial colony whose growth is slowed by sub-inhibitory antibiotic treatment.

Cancer provides an especially strong motivation for investigating the likelihood of selective sweeps during range expansion. Having acquired the necessary driver mutations to escape homeostasis, solid tumours continue to accumulate neutral and driver (that is, cell-fitness-enhancing) mutations as they grow and invade surrounding tissue14,15. The mode of this evolutionary process has sparked debate16,17. An early, highly influential model based on colorectal cancer genetics posited a linear mode of evolution in which cancers acquire mutations through sequential selective sweeps18. Later studies demonstrated intra-tumour heterogeneity with respect to both neutral and driver mutations, suggesting that most cancers undergo branching evolution before they are detected19–23. In some cases, the extent of heterogeneity has been shown to predict clinical outcomes24–26. A possible explanation that reconciles these observations is that linear evolution is restricted to the very early stages of solid tumour evolution and mutations with very strong selective advantage, yet a systematic study of this general, pan-cancer hypothesis is still lacking16,17,27. Mathematical modelling offers a way to investigate the early stages of tumour evolution, which are typically impossible to observe in the clinic28,29. We26,30 and others31–34 have used models with alternative spatial structures and modes of cell dispersal to examine how spatial structure influences tumour evolutionary mode and the extent of intratumour heterogeneity. However, with notable exceptions35–38, theoretical investigations have ignored spatial structure39–41 or relied on agent-based models26,30–34 whose results can be difficult to interpret, provide only limited explanatory insights and are not readily generalisable.

Here we use mathematical analysis to explain why beneficial mutations typically fixate only—if at all—in the very early stages of range expansions, even when mutants can displace the wild type faster than the wild type expands its range. By solving our model in the canonical case of constant radial expansion speed, we derive exact and simple approximate expressions for sweep probabilities in one, two and three dimensions. We confirm the accuracy and robustness of our analytical results using extensive agent-based simulations of a spatial Moran process and we compare outcomes for alternative growth models. We discuss how these findings shed light on the nature of evolution in range expansions in general and cancer development in particular.

Results
Amacroscopicmodelofevolutionduringrangeexpansion
Our macroscopic model is designed to test whether clonal interference alone can prevent selective sweeps and to obtain upper bounds on selective sweep probabilities during range expansions. We consider a wildtype population that starts expanding spherically at time t = 0, such that its radius x grows at speed c_wt. Focusing on selective sweeps, we consider only advantageous mutations, which we assume spread within the wildtype at speed c_m > c_wt (Fig. 1). Mutations occur at per-capita rate μ_e with each surviving genetic drift with probability ρ. In our analytic model, it is enough to consider the compound parameter μ = ρμ_e. For brevity, we will refer to μ as the mutation rate unless otherwise mentioned.

For mathematical and biological reasons (see ‘Discussion’), we focus on a model with constant radial expansion speeds; later, we compare these results with those that pertain to alternative growth models. Various models link propagation speeds to fitness values and migration rates. The most prominent formula is obtained from a reaction-diffusion equation associated with Fisher42 and Kolmogorov43: c = 2\sqrt{aD}, where D is the diffusion coefficient carrying information of migration and a is the difference in the local proliferation rate. More recent studies have sought to refine and generalise this result44–47. Our main results apply to any model that generates approximately constant speeds.

The first surviving mutation achieves a selective sweep only if it reaches every part of the wildtype expansion front before a second mutant of equal or greater fitness arises within the wildtype (Fig. 1). Otherwise, the outcome is clonal interference or possibly a soft selective sweep if the competing mutations are sufficiently similar48. For simplicity, in our macroscopic model, we neglect mutants with fitness values between those of the wildtype and the first mutant, which would slow the expansion of the first mutant and so reduce rather than nullify the selective sweep probability. Neither do we investigate the case of a yet-fitter mutant that arises from the first mutant and achieves a selective sweep. Later, we will examine the effects of relaxing these assumptions.

The unconditional sweep probability is derived in four steps. First, we introduce the random variable X, the radius of the wildtype population when the first mutant arises and we compute its probability density f_X(x). Second, we introduce the random variable Y, the distance between the wildtype and mutant origins and we calculate its probability density conditioned on X, namely f_Y(y|X=x). Third, we derive an expression for the conditional sweep probability Pr(sweep|X=x, Y=y). Finally, we marginalise out X and Y to obtain the unconditional sweep probability Pr(sweep). We focus on the three-dimensional case; analogous results in one and two dimensions are presented in the SI Text Sections 5 and 6.

The following result (proved in SI Text Section 1) will be useful in various parts of subsequent derivations.

Claim 1. The probability that k mutations arise and survive during the time interval [0, t] is Poisson distributed,

P(k) = \frac{e^{-\lambda}\lambda^k}{k!} with \lambda = \mu \int_0^t N(s)ds, (1)

where N(s) is the population size at time s. In particular, the probability that no successful mutant arises is P_0 = e^{-\lambda}.

Although it is commonly assumed that mutations are coupled to divisions, it is straightforward to translate between per-capita and per-division mutation rates (see SI Text Section 2).

Arrival time of the first mutant
In the absence of mutants, the wildtype population in three dimensions grows as N_wt = 4\pi x^3 / 3 = 4\pi(c_wt t)^3 / 3. Applying Claim 1, we obtain the probability that no mutations arise in the time interval [0, t],

P_0 = e^{-\mu \int_0^t 4\pi(c_wt s)^3 ds} = e^{-(\kappa t)^4}, (2)

in terms of a characteristic duration \kappa = \sqrt[4]{\frac{4\pi \mu c_wt^3}{3}}. We identify 1 - P_0 as the cumulative distribution function for the arrival time T of the first surviving mutant. The probability density of T is then

f_T(t) = \frac{d(1 - P_0)}{dt} = 4\kappa^4 t^3 e^{-(\kappa t)^4}, (3)

which is the Weibull distribution with shape parameter 4 and scale parameter \kappa. Substituting x = c_wt t, we find that f_X(x), the probability density of the radius of the wildtype population X at the time the first surviving mutant arises, is the Weibull distribution with shape parameter 4 and scale parameter \theta = \sqrt[4]{\frac{4}{3\pi} \frac{c_wt}{\mu}}. It follows that E[X] \approx 0.91\theta and Var[X] \approx 0.065\theta^2. Analogous calculations in one and two dimensions yield similar Weibull distributions (SI Text Sections 5 and 6 and Fig. 2A–C).

Location of the first mutant
Next, we compute the probability density for the distance Y between the first surviving mutant and the centre of the wildtype population. Since mutants arise in proportion to the number of wildtype cells, we have f_Y(y|X=x)dy \propto D(y), where D(y) is the number of cells at distance y. D(y) corresponds to the infinitesimal shell, D(y) = 4\pi y^2 dy 1_{[0, x]}(y), where the last term is an indicator function that defines the boundary of the wildtype population. After normalisation, we obtain

f_Y(y|X=x) = \frac{3y^2}{x^3} 1_{[0, x]}(y). (4)

We calculate the unconditional probability density of Y by marginalising out X,

f_Y(y) = \int_0^\infty f_Y(y|X=x)f_X(x)dx = \frac{3}{\theta} \frac{y^2}{\theta^2} \Gamma\left(\frac{4}{1}, \frac{y^4}{\theta^4}\right), (5)

where \Gamma(a, z) = \int_z^\infty t^{a-1}e^{-t}dt is the incomplete gamma function. We then find E[Y] \approx 0.68\theta and Var[Y] \approx 0.070\theta^2. Similar results pertain in one and two dimensions (SI Text Sections 5 and 6; Fig. 2A–C).

Conditional sweep probability
To determine the sweep probability, we express the time required for the mutant to sweep, \tau, with \tau = 0 when the first mutant arises. Recall that t = 0 at the origin of the wildtype population, so t = \tau + x/c_wt. Once the mutant population starts expanding, we have N_wt(\tau) = N_wt(\tau) - \Delta(\tau), where \Delta(\tau) is the number of wildtype cells that the mutant has replaced and N_wt(\tau) = 4\pi x^3(\tau) / 3 is the wildtype population size had there been no mutant. While the mutant is surrounded by the wildtype, we have \Delta(\tau) = \Delta_1(\tau) = 4\pi x_m^3(\tau) / 3. After the mutant breaches the wildtype boundary, \Delta(\tau) is given by the intersection formula of two balls, which we denote \Delta_2(\tau) (see SI Text Section 3 for the explicit formula).

Together, we have

N_wt(\tau) = (N_wt(\tau) - \Delta_1(\tau)) 1_{[0, \tau_1]}(\tau) + (N_wt(\tau) - \Delta_2(\tau)) 1_{[\tau_1, \tau_2]}(\tau), (6)

where \tau_1 is the time at which the mutant reaches the wildtype boundary and \tau_2 is the time at which the mutant has entirely replaced the wildtype. Noting that

\tau_1 = \frac{x-y}{c_m - c_wt} (7)

and

\tau_2 = \frac{x+y}{c_m - c_wt}, (8)

we apply Claim 1 to obtain the conditional sweep probability

Pr(sweep|X=x, Y=y) = e^{-\mu \int_0^{\tau_2} N_wt(\tau)d\tau}. (9)

Although this integral can be solved analytically (SI Text Section 3), the resulting expression is complicated and not very enlightening and we therefore seek simpler approximations. An especially fruitful approach is to assume y = 0, so that f_Y(y|X=x) = \delta(y), where \delta is the Dirac delta function. This yields an upper bound on the conditional sweep probability because, due to geometrical symmetry, the time required for a mutant to sweep (\tau_2 in eqn. (8)) is minimal when y = 0. With this approximation, eqn. (9) simplifies drastically as \tau_1 = \tau_2 = \frac{x}{c_m - c_wt}. We do not need to integrate over \Delta_2(\tau). The analytic solution is

Pr(sweep|X=x, Y=0) = e^{-(\alpha x)^4}, (10)

where \alpha = \sqrt[4]{\frac{4\pi \mu (c_wt^2 - 3c_wt c_m + 3c_m^2)}{3(c_m - c_wt)^3}} is a characteristic length. Despite being based on a seemingly crude approximation, eqn. (10) is remarkably close to the exact conditional sweep probability (Fig. 2D). The corresponding approximations in one and two dimensions (SI Text Sections 5 and 6) are likewise simple and useful.

Unconditional sweep probability
We now have all the necessary ingredients to compute the unconditional sweep probability by marginalising out X and Y from the conditional sweep probability,

Pr(sweep) = \int_0^\infty \int_0^x Pr(sweep|X=x, Y=y)f_Y(y|X=x)f_X(x)dydx

= \int_0^\infty \int_0^x e^{-\mu \int_0^{\tau_2} N_wt(\tau)d\tau} \frac{3y^2}{x^3} \frac{4x^3}{\theta^4} e^{-x^4/\theta^4} dydx. (11)

As before, we use the approximation y = 0 to obtain the strikingly simple result

Pr(sweep) \le \int_0^\infty Pr(sweep|X=x, Y=0)f_X(x)dx

= \int_0^\infty e^{-(\alpha x)^4} \frac{4x^3}{\theta^4} e^{-x^4/\theta^4} dx = \frac{1}{1 + (\frac{c_m - c_wt}{c_m})^3}. (12)

This result generalises to analogous one- and two-dimensional models, with the exponent 3 replaced by the respective spatial dimension (SI Text Sections 5 and 6). The approximate expressions are close to the exact solution in one dimension and to numerical approximations of the exact solution in two and three dimensions (Fig. 2E). In the latter case, the approximation y^2 = 0.28125x^2 instead of y = 0 and hence shed light on why the upper bound is close to the exact solution (SI Text Section 5). The sweep probability is independent of the mutation rate, not only in these approximations but also in the general case of Eqn. (11). This follows from a more general result (SI Text Section 4). Eqn. (12) also leads to a simple upper bound on the probability of multiple sequential sweeps (SI Text Section 10).

Together, we have

Pr(sweep) \approx 1 - \left(\frac{r_wt - r_re}{r_m - r_re}\right)^d, (13)

where r_m, r_wt and r_re are the proliferation rates of the mutant invader, wildtype invader and resident populations, respectively and d is the spatial dimension.

Conditional arrival time of the first mutant
In biological systems where it is infeasible to track evolutionary dynamics, selective sweeps must be inferred from subsequent genetic data. For example, we might observe a public mutation in a tumour and ask when this mutation occurred. We can use our model to obtain the probability distribution of the radius X at the time the first mutant arose, given that we observe a selective sweep, by applying Bayes’ theorem,

f_X(x|sweep) = \frac{Pr(sweep|X=x)f_X(x)}{Pr(sweep)}. (14)

Using the approximation y = 0, we obtain the Weibull distribution

f_X(x|sweep) \approx \frac{4x^3}{\theta^4} \left(\frac{c_m}{c_m - c_wt}\right)^3 e^{-\frac{x^4}{\theta^4} (\frac{c_m}{c_m - c_wt})^3}, (15)

where \beta = \frac{c_m}{c_m - c_wt}. This approximation and the corresponding results in one and two dimensions are close to the exact solutions (Fig. 2F).

Sweepprobabilitiesinalternativegrowthmodels
To examine the robustness of our findings to variations in the model assumptions, we compare them to results obtained from alternative growth models. Instead of permitting mutants to grow within the wildtype population, one might instead assume that dispersal is restricted to the wildtype population boundary. In this case, a complete sweep is impossible and we instead ask whether the first arising mutant envelopes the wildtype. Interestingly, this envelopment probability is independent of the mutation rate and obtains values close to the sweep probabilities of our main model36 (Fig. 3A and SI Text Section 7). Although our focus is on expanding populations, we also applied our methods to compute the sweep probability in constant populations (SI Text Section 7), obtaining a result that depends on both mutation rate and population size (Fig. 3), in agreement with prior analyses50,51.

We have also investigated sweep probabilities in non-spatial models. In the case of exponential growth with no competition, we find that a single mutant is unlikely to become dominant unless its exponential growth rate is several times larger than that of the wildtype and that the sweep probability is relatively insensitive to the mutation rate (Fig. 3B and SI Text Section 7). In the case of constant population size and logistically growing mutant, the sweep probability is sensitive to both population size and mutation rate (Fig. 3B and SI Text Section 7). For a non-spatial population of constant size, the sweep probability decreases as the mutation rate increases, as in the spatial case. In summary, we find that the main conclusions drawn from our primary model also hold for alternative models of range expansions (exponential or boundary growth) but not for non-growing or tightly bounded populations.

Validationofourmacroscopicmodelusingagent-basedsimulations
To gauge the robustness of our theoretical predictions to the effects of stochastic growth and discretization, we measure the frequency of selective sweeps in a two-dimensional agent-based model. We suppose that the wildtype population invades a habitat initially occupied by a resident competitor, which is a plausible biological scenario for both invasive species13 and invasive tumours26,30. Our agent-based model has resident, wildtype invader and mutant invader populations with local proliferation rates r_re, r_wt and r_m,i. Localised competition between wildtype and resident slows the wildtype expansion and creates potential for selective sweeps. For consistency with our macroscopic model, we parametrised the simulations so that all mutations have equal effect (r_m,i = r_m for all i). We implemented this model using the demon agent-based modelling framework52 within the warlock computational workflow53, which facilitates running large numbers of simulations on a high-performance computing cluster. We have previously applied the same framework to studying cancer evolution30,53. Further model details are given in Methods.

The agent-based simulations provide useful insights in addition to the macroscopic model because, although the general setup is the same, they differ in several ways. Space in the simulations is divided into discrete patches; the times between birth and dispersal events are stochastic; population boundaries are not smooth; the expansion wavefront is typically not sharp and changes shape as the wave progresses; and the mutant will have increased propagation speed when competing with the resident population. Hence, we would not expect perfect agreement between the results of the two models.

Given the considerable differences between the models, the probability density functions resulting from the macroscopic and microscopic models are reassuringly consistent. The radius at the time the first surviving mutant arises is slightly lower in the simulations than in our analytical model (Fig. 4A). Similarly, the distribution for the location of the first surviving mutant coincides well except for a small offset of the mean (Fig. 4B). Such offsets are expected due to discretization effects and the fact that, in the simulations, the propagation front needs to be established before the expansion can proceed. The sweep probabilities in simulations are nevertheless very close to our analytical predictions and change very little when varying the mutation rate over orders of magnitude, confirming our prediction that the sweep probability is independent of the mutation rate (Figs. 4C and Supplementary Figs. S1–S4).

Randomfitnesseffects
Since the effects of mutations are, in reality, not equal but random, we next used our agent-based model to test whether our predictions remain informative when mutation fitness effects are drawn from an exponential distribution, while still assuming that no mutant can accumulate more than one mutation. Because the random fitness effect determines not only a mutant’s expansion speed but also its probability of evading stochastic extinction, it is more appropriate to compare results in terms of the mean fitness effects of the contending mutations that evade stochastic extinction, rather than the mean fitness of all generated mutations55,56 (see ‘Methods’). In terms of this e_s (Fig. 5C) or an alternative summary statistic (Supplementary Fig. S5; SI Text Section 9), sweep probabilities in random-fitness-effect simulations are close to the predictions obtained under the assumption of fixed effects. Although it is possible for many fitter second or subsequent mutants to arise in the wildtype population and achieve a selective sweep by replacing the wildtype and all less-fit competitor mutants, in our simulations this scenario contributes very little to the total sweep probability.

Multiplicativefitnesseffects
Next, we examined the consequences of allowing individuals to acquire multiple mutations with multiplicative fitness effects. The accumulation of beneficial mutations has two opposing effects on the sweep probability. By acquiring further mutations, a mutant lineage can become fitter, expand faster and hence increase its probability of achieving a selective sweep. On the other hand, competitor lineages can also become fitter and more likely to prevent a sweep. By grouping simulations that have the same distribution of mutation fitness effects, we find that the first factor outweighs the second, so that sweep probabilities are higher in the case of cumulative mutation effects (Fig. 5C). This result makes intuitive sense because the lineage closest to achieving a sweep represents the largest target for further fitness-enhancing mutations. Allowing the accumulation of beneficial mutations also increases the median time at which sweeps finish, in particular when fitness advantages are weak (Fig. 5D and Supplementary Figs. S11–S14). Nevertheless, the conclusion remains that, even when mutations have multiplicative effects, selective sweeps are infrequent unless contending mutations are on average strongly beneficial, in which cases sweeps occur early.

When the number of mutations that can accumulate is unlimited, we find that, for every distribution of fitness effects, a single selective sweep is a more common outcome than two sweeps and two sweeps occur more often than three (Fig. 5E). Only rarely do we observe more than three sweeps. This conclusion is robust to the caveat that, for small s values, we likely underestimate the frequency of multiple sweeps because some would complete after we terminate the simulations. By plotting sweep probability against the mean fitness increase over the course of the simulation, we find that sweeps are more often due to a single large-effect mutation than several mutations with small effects (Supplementary Fig. S5).

Applicationtocancer
Our findings have especially interesting implications for understanding cancer evolution, in which case we consider a wild-type tumour that evolves while invading a resident population of normal cells. The macroscopic model then has only three parameters: the mutation rate conditional on survival, \mu and the wildtype and mutant propagation speeds, c_wt and c_m. To estimate c_wt for human tumours, we follow a similar procedure to prior studies36,57. Consider a tumour that grows to a volume of V between 1 and 10 cm^3 in time T between 5 and 20 years. The propagation speed can then be estimated as

c_wt = \sqrt[3]{\frac{4\pi V}{3T}}, (16)

which equates to between 1 and 40 \mu m per day. Given that the diameter of a typical cancer cell l \approx 20 \mu m58 and the generation time (cell cycle time) is \tau_G \approx 4 days59, we can switch units to obtain c = c_wt \times \tau_G / l, which is between 0.15 and 7 cell diameters per generation. The rate of acquiring advantageous (driver) mutations is usually estimated at around \mu_e = 10^{-5} per cell per generation57,60. For the survival probability \rho, we assume values between 0.09 and 0.5, in agreement with inferred values in colorectal tumours61. Together, these parameter values lead to a typical length scale \theta = \sqrt[4]{\frac{4}{3\pi} \frac{c_wt}{\mu}} \approx 10 to 50 cells. The expected values of X and Y are then E[X] = 9 to 45 cell diameters and E[Y] = 7 to 35 cell diameters. It follows that a tumour, having acquired sufficient mutations to grow, will likely gain further driver mutations already during early development.

By eqn. (12), if a mutant expands 10% faster within a tumour than the tumour expands into surrounding tissue, then the sweep probability is predicted to be less than 1 - (1/1.1)^3 \approx 0.00075. If the mutant propagates twice as fast as the wildtype, we have Pr(sweep) < 0.13 and in the extreme case where c_m is ten times c_wt, we have Pr(sweep) < 0.73. In the latter case, the expected population radius when the sweep began, given that a sweep occurred, is approximately 50 cells, corresponding to a population size N of 400,000 cells. The time for the sweep to be completed is then \tau_2 \approx 40 generations, at a population size of approximately 800,000 cells. This result is relatively robust because the expected values E[X], E[Y] and E[X|sweep] are proportional to the characteristic length \theta, which varies with the fourth root of \mu and c_wt. Hence, in three dimensions, changing \mu or c_wt by a factor of 10 changes our estimates of the radius and time by only a factor of 10^{1/4} \approx 1.8 and estimates of the population size by a factor of 10^{3/4} \approx 5.6. In summary, our macroscopic model (assuming equal mutation effects and no accumulation of mutations) predicts that selective sweeps during a clonal expansion are rare unless mutations are very strongly beneficial, in which cases sweeps begin and end early in the expansion. This general conclusion also holds for alternative growth models (Fig. 3).

Methods
Numericalintegration
We performed numerical integration using the MATLAB function ‘trapz’. For values of x close to 0 and beyond 2\theta, f_X(x) is small and its contribution to the integrals is negligible. Hence, to avoid numerical errors without compromising precision, we set the lower and upper bounds of integration to 0.001\theta and 3\theta when integrating over f_X(x). Similarly, we set 0.001\theta as the lower bound of the integration over y. We set interval widths to 0.001\theta for x \le 0.5 and 0.01\theta for x > 0.5.

Agent-basedsimulations
We ran agent-based simulations using the Warlock automated computational workflow for the demon modelling framework30,53. Individuals in this agent-based model are subdivided into well-mixed demes on a regular two-dimensional grid. The demes have identical carrying capacities, K and are initially filled with residents, except that a single wildtype invader is introduced at the centre of the grid. At each time step, an individual is chosen at random, with probability weighted by fitness, to be replaced by two offspring. Each offspring then either migrates, with probability m, to a neighbouring deme in a randomly chosen direction or remains in its parent deme. Local density dependence is implemented by imposing a very high death rate whenever a deme is above the carrying capacity. Mutation is coupled to wildtype reproduction. Motivated by an invasive cancer cell population spreading in non-invasive healthy tissue, we do not permit resident individuals to disperse. We can account for this asymmetry by adapting the conversion from fitness advantage into expansion speeds (to be published in a later study). Further model details have been published previously30,53. All simulations were performed using City, St George’s, University of London’s Hyperion cluster. Example graphical outputs of the simulations, plotted using the ggmuller R package74, are included in Supplementary Figs. S7 and S8.

For all simulations, we set the deme size to K = 16 and the migration probability m = 0.05. For such a small m, the probability of surviving drift can be estimated as the probability of becoming locally fixed in one deme. Since the within-deme dynamics approximate a Moran process, this probability is \rho = 1 - (1 - r_wt/r_m)^K, where r_wt and r_m are the proliferation rates75. We measured propagation speeds by simulating the expansion of a wildtype into a large resident population, or of a mutant into a large wildtype population, in the absence of mutation and then applying linear regression to the growth curve of the effective radius, defined as the square root of the population size divided by \pi.

For each parameter set, we took the mean speed from ten replicates, for which the standard deviation was consistently below 2% of the mean. We set the proliferation rates to r_re = 0.91 and r_wt = 1.0 for the resident and wildtype, respectively, leading to a_wt = r_wt - r_re = 0.09. The measured wildtype speed was then c_wt \approx 0.15. In simulations with fixed mutation effects, we varied the mutant proliferation rate r_m from 1.1 to 3.0, so that a_m = r_m - r_wt ranged from 0.1 to 2.0. The measured mutant speeds are presented in Supplementary Table S3.

To implement random mutation effects, we multiplied the fitness effect of each mutation by a factor X drawn from an exponential distribution with mean 1. To prevent birth rates from being implausibly large, we assumed diminishing-returns epistasis associated with maximum birth rate M = 10. The birth rate after i mutations was then

r_m,i = min(r_m,i-1 + sX(1 - r_m,i-1/M), M), (17)

where parameter s > 0 determines the mean mutation effect.

Computationofsweepprobabilitiesinsimulationswithfixedmutationeffects
We ran a batch of 1000 simulations for each combination of mutation rate \mu_e \in {10^{-4}, 10^{-5}, 10^{-6}} and mutant proliferation rate r_m \in {1.1, 1.2, 1.3, 1.4, 1.5, 2.0, 2.5, 3.0}. Resident and wildtype proliferation rates were in all cases r_re = 0.91 and r_wt = 1.0. Simulations were stopped at a population size of 1,000,000 or 2000 generations (whichever was reached first). We called a sweep when all individuals at the end of the simulations shared a mutation.

To investigate the extent to which we failed to detect selective sweeps that would have completed after the simulation stop time, we examined the population sizes when the detected sweeps finished. For every set of parameter values for which sweeps were abundant enough to analyse, the population radius at the time of sweep completion exhibited a bell-shaped distribution that decayed steeply before the maximum possible radius \sqrt{10^6/\pi} \approx 564, indicating that we detected the vast majority of sweeps. In the remaining cases (\mu_e = 10^{-5} and r_m = 1.1; \mu_e = 10^{-6} and r_m = 1.1 or 1.2), sweeps were sufficiently rare that even relatively severe undercounting would not affect our conclusions.

Computationofsweepprobabilitiesinsimulationswithrandommutationeffects
For each variant of the random-mutation-effects model (permitting the accumulation of up to one, two, three or unlimited mutations), we ran eleven batches of 1000 simulations, with fitness effects calculated according to eqn. (17) and s \in {0.05, 0.075, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, 1, 1.5, 2}. As before, simulations were stopped when they reached a population size of 1,000,000 or 2000 generations.

The establishment probability \rho is lower for weakly beneficial than for highly beneficial mutations, which implies that the effective mutation rate \mu = \rho\mu_e is lower in the former case. For example, in a Moran process with K = 16, nearly 94% of nearly-neutral mutations succumb to stochastic extinction before they can contribute to the macroscopic evolutionary dynamics, compared to only 50% of mutations that confer fitness effects s = 1. To account for this difference in effective mutation rates, we calculated, in the at-most-one-mutation model, the mean effect e_s of mutations that reached an abundance of at least 10 individuals (see Supplementary Fig. S16 for the distributions). We then compared the sweep probability for each random-effects model with that of the equal-effects model with a_m = e_s. As before, we transformed a_m into c_m according to Supplementary Table S3.

To investigate whether our finite-time simulations omit sweeps that would occur if simulations were run for longer times, we investigated the population radii at sweep completion. Population size rather than number of generations was the dominant stopping criterion except when s = 0.05 in all models and s = 0.075 in the at-most-one-mutation model (Supplementary Fig. S9). In all cases, the radius distributions are bell-shaped and indicate that, for all but the smallest fitness effects, most selective sweeps finished long before the simulation stop time (Supplementary Figs. S11–S14). To estimate the extent of undercounting, we fitted right-truncated gamma distributions to the radius distributions, with the truncation at \sqrt{1,000,000/\pi} \approx 564 corresponding to the population size at which simulations were stopped. Using a maximum likelihood approach, we then estimated the scale and shape parameters of each gamma distribution. By evaluating the probability in the tail of the non-truncated gamma distribution beyond the truncation point, we estimated the fraction of sweeps omitted in each batch of simulations and adjusted the sweep probabilities in Fig. 5C accordingly. The adjustments are negligible except in the case of small fitness effects, for which sweep probabilities are in any case low. In tests using artificially truncated distributions from simulations with large s values (Supplementary Fig. S15), we find that this method typically overestimates the proportion of missing sweeps. The adjusted values can therefore be regarded as upper bounds.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Stein, Kizhuttil, Back, and Noble present a mathematical model for studying how an advantageous mutation spreads in a wildtype population and dominates it. The authors are interested in the sweep probability, which is the probability that the mutant and its descendants replace the wildtype population. The main scenario consists of wild type and mutant populations that expand radially in space. The authors derived the sweep probability and showed that it is independent of the mutation rate. Next, they studied a three-type agent-based model for verification of the theory. In addition, several other spatial and non-spatial scenarios are considered for comparison with the main scenario. Lastly, the authors applied the theory to cancer evolution.\nThe main text is well organized, and the mathematical derivations are easy to follow. The observation that the sweep probability is independent of the mutation rate is enlightening.\nHere are some other comments for the main text:\n1. Figure 2: Why is the shape of 1-d Y density so different from the 2-d and 3-d Y density? It seems that the formulas have a common form of the incomplete Gamma function.\n2. Line 157, display (8): tau1 before the double arrow should be tau2.\n3. Line 184: It would be nice if the authors give the full name of FKPP here. (There is one in line 223 but it might be better to put it here.)\n4. Line 184, display (13): It's hard to find the derivation of this approximation in the main text and the SI. Can the authors point it out explicitly here?\n5. Figure 3, panel B: The colors are different from panel (A) so there should be a legend here. In addition, could the author provide the number of simulations for generating these plots? It seems that the bin width of the histogram in (A) and (B) are different.\n6. Line 262: The location of the parenthesis is incorrect.\n7. Line 306: The discussion on the conditions that the sweep probability is independent of the mutation rate is great. I wonder if there are real biological cases in which these conditions are nearly satisfied and we can observe the independence. In addition, do the authors have an intuition on the approximation of locating the first mutant at the center? I think this may be due to the geometric symmetry of the expansion as well as the assumption that the mutant arises uniformly among the population.\nThere are some typos in the Supplementary:\n1. Line 14: The dummy variable of in the integral should not be t.\n2. Line 28, display [1]: Parentheses is missing.\n3. Line 150: There is a missing reference.\n4. Line 154, display [32]: The display formula should be corrected.\n5. Line 159, display [34]: The ‘wt’ in the last line should be a subscript.", "reviewer_2": "In this paper, Stein et al derive mathematical expressions for selective sweep probabilities in simple models of expanding populations. The authors are motivated by tumors, which, through recent multi-site and single cell sequencing efforts, have been revealed to be highly genomically heterogeneous. The authors place this at odds with older models from Vogelstein and others which hypothesize that tumors acquire driver mutations through sequential selective sweeps, which would reduce heterogeneity.\nThe authors' central claim is that range expansions restrict selective sweeps to small tumors (of the order tens of cells), which they show by analyzing a very simplified model of range expansion. In essence, they model a linearly expanding population (WT) as an initial tumor clone. They allow for a second driver (MT) to emerge at a prescribed mutation rate on this expanding WT background. In general, it is assumed the expansion speed of MT > expansion speed of WT. They define a fixation event as a case where MT completely displaces WT before another mutation arises on the WT background. This gives simple expressions for the fixation probability, since the problem, as laid out, can be cast in purely geometrical terms.\nGiven this relatively simple setup, the mathematical analysis looks sound to me. And I very much appreciate the attempt to place the model in context by estimating parameters in the \"Application to cancer\" section.\nHowever, I have broad concerns about the conceptual foundation of the model, i.e. simple models are good, but I am not convinced that this model is sufficiently complex to capture some of the more salient aspects of cancer evolution.\nMajor concern:\nThe most serious concern is about the motivation of the model itself -- is this the right model to study if cancer is the point of comparison? It is not clear to me that selective sweeps of the sort described are a relevant scenario in cancer. Rather, selective sweeps function as a sort of straw man in the manuscript. The model as laid out describes a highly restricted scenario, where only mutations on the WT background are accounted for. However, if the MT population is expanding faster or once the MT is larger than the WT (if mutations are per birth or per capita, respectively), it is the MT that is likely to get more subsequent beneficial mutations. These additional drivers on the MT background could abrogate the observed effect of reduced fixation of the MT allele (since drivers emerging on the WT background would have to catch up to carriers of the MT allele that have acquired additional beneficial mutations, assuming a fixed fitness effect size).\nThis is one relatively small tweak that can have a drastic impact on the conclusions that are drawn from the model. The authors themselves mention in the text and describe in the SI that the key result that the sweep probability is independent of the mutation rate is highly dependent on the details of the growth process. In fact, `either a positive or a negative effect of mutation rate on sweep probability' can be observed depending on the underlying growth dynamics. Cancers typically grow in highly heterogeneous environments, and exhibit different growth properties depending on the cell type of origin or metastatic site -- so this is likely a relevant perturbation to the model that can drastically change the results.\nIndeed, if we look to the literature, it has also been shown that, in spatially structured models parametrized by cancer contexts, local measures of heterogeneity are limited due to cell turnover (cf, Waclaw, et al. `A spatial model predicts that dispersal and cell turnover limit intratumour heterogeneity'). This is in stark contrast to what the authors conclude -- that spatial structure increases heterogeneity. My understanding is that this is due in part to the fact that the present authors primarily focus on the sweep probability (a global measure of loss of heterogeneity). However, in a spatial context, spatially continuous measures of heterogeneity might be more appropriate to use and to compare, at least over the relevant observation timescales.\nTaken collectively, these issues indicate to me that the results are of somewhat restricted applicability, and are highly dependent on the particular, relatively simple formulation of the model.\nMinor concern:\nPlease describe the \"warlock\" code in more detail, perhaps illustrating what its output looks like -- the manuscript would greatly benefit from having this be described in a self-contained way, even if the software was published elsewhere.\nThe code is contained in a rather large compressed file (>4GB) -- it would be much more preferable to store the code relevant to the paper (which is presumably somewhat simple, and significantly smaller than 4GB), in a Github repository."}, "author_response": "Response to the reviewers’ comments\nReviewer #1 (Remarks to the Author):\nStein, Kizhuttil, Back, and Noble present a mathematical model for studying how an\nadvantageous mutation spreads in a wildtype population and dominates it. The\nauthors are interested in the sweep probability, which is the probability that the\nmutant and its descendants replace the wildtype population. The main scenario\nconsists of wild type and mutant populations that expand radially in space. The\nauthors derived the sweep probability and showed that it is independent of the\nmutation rate. Next, they studied a three-type agent-based model for verification of\nthe theory. In addition, several other spatial and non-spatial scenarios are\nconsidered for comparison with the main scenario. Lastly, the authors applied the\ntheory to cancer evolution.\nThe main text is well organized, and the mathematical derivations are easy to follow.\nThe observation that the sweep probability is independent of the mutation rate is\nenlightening.\nWe are glad that the reviewer found our manuscript to be clear and enlightening.\nHere are some other comments for the main text:\n1. Figure 2: Why is the shape of 1-d Y density so different from the 2-d and 3-d Y\ndensity? It seems that the formulas have a common form of the incomplete Gamma\nfunction.\nThe different shapes can be understood from a mathematical perspective as follows:\nBy definition, the incomplete gamma function is a monotonically decreasing function\nof y>0. This explains the monotone behaviour in the 1-d case, for which the\nprobability density function (pdf) is simply proportional to the incomplete gamma\nfunction (Table S1). In the other cases, the pdf also contains a factor of y (for the 2-d\nmodel) or y^2 (3-d). A Taylor expansion shows that the factor of y or y^2 dominates\nwhen y is small, so that the pdf is zero when y=0 and increases approximately\nlinearly (2-d) or quadratically (3-d) when y is small. When y is large, the incomplete\ngamma function dominates and the curve decreases faster than exponentially.\n2. Line 157, display (8): tau1 before the double arrow should be tau2.\nWe corrected the typo.\n3. Line 184: It would be nice if the authors give the full name of FKPP here. (There is\none in line 223 but it might be better to put it here.)\nWe now do so.\n4. Line 184, display (13): It's hard to find the derivation of this approximation in the\nmain text and the SI. Can the authors point it out explicitly here?\nWe have added the necessary details (lines 195-198).\n5. Figure 3, panel B: The colors are different from panel (A) so there should be a\nlegend here. In addition, could the author provide the number of simulations for\ngenerating these plots? It seems that the bin width of the histogram in (A) and (B)\nare different.\nWe added a legend to panel B and made the bin widths identical. The number of\nsimulation replicates was already stated in the caption (we have put it in red text for\nease of reference).\n6. Line 262: The location of the parenthesis is incorrect.\nWe corrected the typo.\n7. Line 306: The discussion on the conditions that the sweep probability is\nindependent of the mutation rate is great. I wonder if there are real biological cases\nin which these conditions are nearly satisfied and we can observe the independence.\nIn addition, do the authors have an intuition on the approximation of locating the first\nmutant at the center? I think this may be due to the geometric symmetry of the\nexpansion as well as the assumption that the mutant arises uniformly among the\npopulation.\nDespite searching, we have been unable to find any empirical study of how sweep\nprobabilities in range expansions vary with mutation rate. Detecting selective sweeps\nin natural systems is generally difficult. One possibility would be to design lab\nexperiments like those of Hallatschek and colleagues using fluorescently labelled E.\ncoli (Hallatschek et al. PNAS, 2007), as we now mention in our Discussion (lines\n459-460).\nWe have added an explanation of why the upper bound on the sweep frequency,\nobtained by setting y = 0, is a good approximation to the exact result in the one-\ndimensional case. Linked to this explanation, we provide new, better approximations\nin the one-dimensional case (main text lines 185-188; supplementary text lines 119-\n128).\nThere are some typos in the Supplementary:\n1. Line 14: The dummy variable of in the integral should not be t.\n2. Line 28, display [1]: Parentheses is missing.\n3. Line 150: There is a missing reference.\n4. Line 154, display [32]: The display formula should be corrected.\n5. Line 159, display [34]: The ‘wt’ in the last line should be a subscript.\nWe have corrected the typos. Thank you for taking the time to check our manuscript\nso thoroughly.\nReviewer #2 (Remarks to the Author):\nIn this paper, Stein et al derive mathematical expressions for selective sweep\nprobabilities in simple models of expanding populations. The authors are motivated\nby tumors, which, through recent multi-site and single cell sequencing efforts, have\nbeen revealed to be highly genomically heterogeneous. The authors place this at\nodds with older models from Vogelstein and others which hypothesize that tumors\nacquire driver mutations through sequential selective sweeps, which would reduce\nheterogeneity.\nThe authors' central claim is that range expansions restrict selective sweeps to small\ntumors (of the order tens of cells), which they show by analyzing a very simplified\nmodel of range expansion. In essence, they model a linearly expanding population\n(WT) as an initial tumor clone. They allow for a second driver (MT) to emerge at a\nprescribed mutation rate on this expanding WT background. In general, it is\nassumed the expansion speed of MT > expansion speed of WT. They define a\nfixation event as a case where MT completely displaces WT before another mutation\narises on the WT background. This gives simple expressions for the fixation\nprobability, since the problem, as laid out, can be cast in purely geometrical terms.\nGiven this relatively simple setup, the mathematical analysis looks sound to me. And\nI very much appreciate the attempt to place the model in context by estimating\nparameters in the \"Application to cancer\" section.\nHowever, I have broad concerns about the conceptual foundation of the model, i.e.\nsimple models are good, but I am not convinced that this model is sufficiently\ncomplex to capture some of the more salient aspects of cancer evolution.\nWe are grateful to the reviewer for raising these concerns, which we have addressed\nby running extensive new simulations and analyses. We believe that the new results\nsubstantially strengthen our manuscript without changing its main conclusions.\nMajor concern:\nThe most serious concern is about the motivation of the model itself -- is this the right\nmodel to study if cancer is the point of comparison? It is not clear to me that\nselective sweeps of the sort described are a relevant scenario in cancer. Rather,\nselective sweeps function as a sort of straw man in the manuscript.\nWe have rewritten parts of the Introduction (lines 48-69) and Discussion (lines 425-\n436 and 451-455) and added more references to clarify how and why our model is\nrelevant to cancer.\nThe model as laid out describes a highly restricted scenario, where only mutations\non the WT background are accounted for. However, if the MT population is\nexpanding faster or once the MT is larger than the WT (if mutations are per birth or\nper capita, respectively), it is the MT that is likely to get more subsequent beneficial\nmutations. These additional drivers on the MT background could abrogate the\nobserved effect of reduced fixation of the MT allele (since drivers emerging on the\nWT background would have to catch up to carriers of the MT allele that have\nacquired additional beneficial mutations, assuming a fixed fitness effect size).\nWe agree that it is important to examine this scenario. We now report extensive new\nsimulations and analyses of models that permit mutants to accumulate beneficial\nmutations with random, multiplicative fitness effects (main text lines 258-309 and\n500-508; supplementary text Sections 9 and 10; Figures 4 and S5). We conclude\nfrom these new results that, “consistent with the predictions derived from our\nmacroscopic model, a selective sweep is highly probable only if a mutant with a large\nfitness advantage arises in the very early stages of a range expansion”.\nThis is one relatively small tweak that can have a drastic impact on the conclusions\nthat are drawn from the model. The authors themselves mention in the text and\ndescribe in the SI that the key result that the sweep probability is independent of the\nmutation rate is highly dependent on the details of the growth process. In fact, `either\na positive or a negative effect of mutation rate on sweep probability' can be observed\ndepending on the underlying growth dynamics. Cancers typically grow in highly\nheterogeneous environments, and exhibit different growth properties depending on\nthe cell type of origin or metastatic site -- so this is likely a relevant perturbation to\nthe model that can drastically change the results.\nThe sweep probability is independent of the mutation rate not only in our main model\nbut also in the case of boundary-driven growth. In the case of exponentially\nexpanding populations, we observe only a weak dependence on mutation rate\n(Figure 5, green curves). Only in the case of constant or tightly bounded population\nsize, which is outside the scope of our study, do we observe a strong dependence\non mutation rate. We have revised the relevant text (lines 315-317 and 325-332) to\nclarify these points. We conclude, “In summary, we find that the main conclusions\ndrawn from our primary model also hold for alternative models of range expansions\n(exponential or boundary growth) but not for non-growing or tightly bounded\npopulations.”\nIndeed, if we look to the literature, it has also been shown that, in spatially structured\nmodels parametrized by cancer contexts, local measures of heterogeneity are limited\ndue to cell turnover (cf, Waclaw, et al. `A spatial model predicts that dispersal and\ncell turnover limit intratumour heterogeneity'). This is in stark contrast to what the\nauthors conclude -- that spatial structure increases heterogeneity. My understanding\nis that this is due in part to the fact that the present authors primarily focus on the\nsweep probability (a global measure of loss of heterogeneity). However, in a spatial\ncontext, spatially continuous measures of heterogeneity might be more appropriate\nto use and to compare, at least over the relevant observation timescales.\nThe question of how spatial structure influences tumour evolutionary mode and\nintratumour heterogeneity was the focus of our previous study (Noble et al., Nature\nEcology & Evolution, 2022). In that paper we used extensive agent-based\nsimulations to examine not only a version of the model studied by Waclaw et al.\n(2015) but also models with alternative spatial structures, including the one we use in\nthe current study. Our 2022 paper specifically discusses how results obtained using\nthe agent-based model we use in the current study relate to the results of Waclaw et\nal. (2015). Consistent with Waclaw et al. (2015) and other prior studies, we broadly\nconcluded that spatial structure inhibits selection and thus increases intratumour\nheterogeneity, and that, among spatial models, those with higher cell turnover and\nhigher cell dispersal permit more selection and generate lower intratumour\nheterogeneity. We have also previously used spatial agent-based modelling to\ninvestigate how selective sweeps can purge intratumour heterogeneity (Noble &\nBurley et al., Evolutionary Applications, 2020). The scope of the current work is to\nobtain analytical results that provide general, mechanistic insights into one important\naspect of the evolutionary dynamics observed in these prior studies: the frequency of\nselective sweeps. We expect that our mathematical analysis could be extended to\nprovide insights into the dynamics of intratumour heterogeneity but this is beyond the\nscope of the current manuscript. We have added to the Introduction (lines 63-69)\nand the Discussion (lines 456-459) to clarify these points.\nTaken collectively, these issues indicate to me that the results are of somewhat\nrestricted applicability, and are highly dependent on the particular, relatively simple\nformulation of the model.\nAs above, we have added extensive new results and revised the text to show that\nour conclusions also hold for alternative models.\nMinor concern:\nPlease describe the \"warlock\" code in more detail, perhaps illustrating what its output\nlooks like -- the manuscript would greatly benefit from having this be described in a\nself-contained way, even if the software was published elsewhere.\nWe have added plots showing example model output (Figure 4A-B, S7 and S8) and\nmore details of the model (lines 500-508).\nReviewer #2 (Remarks on code availability):\nThe code is contained in a rather large compressed file (>4GB) -- it would be much\nmore preferable to store the code relevant to the paper (which is presumably\nsomewhat simple, and significantly smaller than 4GB), in a Github repository.\nWe have created such a GitHub repository and cited it in the Code Availability\nsection."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have addressed all my comments, and I recommend the manuscript for publication.", "reviewer_2": "Stein et al have amended their manuscript in two primary ways, adding more simulation results and expanding introductory and discussion points. I appreciate that effort has been put in to address my points of concern. However, the presentation of some of the additional results (extending the simulation model to random and accumulative fitness affects) come across as confusing, and potentially incorrect. Moreover, the relationship of these new simulation results to the analytical model is somewhat tenuous. Finally, the key predictions of the model wrt its relationship to cancer data seem to be demonstrably inconsistent with observations described in the literature.\nMajor points:\nI will first address the new individual-based simulations (those pertaining to Fig 4 in the revised text). There are two new simulations that are introduced:\n1. Simulations in which random effect size mutations are introduced on the WT background, but further mutations cannot accumulate.\n2. Simulations in which mutations can continue to accumulate “multiplicatively.”\nThe simulations in which random effect single mutations are introduced are problematic conceptually. The study aims to describe sweep probabilities, which are typically infinite time quantities, but individual-based models can only be run for finite time (presumably 1000 generations in this study, according to figure 4A, though this is never stated). The precise time to see a sweep (and therefore, the appropriate time to run a simulation) will sensitively depend on several underlying parameters like mutation rate, the mode of overall population growth (fixed, boundary-drive, exponential), and the distribution of fitness effects. Therefore, it is not clear that the “macroscopic” model and analytical results derived from it are the appropriate point of comparison for finite time data generated by individual-based simulations in the case where mutation effect sizes vary. Indeed, if the simulations were run for longer, we would eventually see an atypically large effect size beneficial mutation arise and eventually sweep in every simulation by chance (assuming an exponential distribution of effect sizes as the authors state).\nMoreover, the “anomaly” at low fitness observed in these simulations and described in Figure 4C is confusing, and I’m not sure the explanation provided in the text is satisfactory. Specifically, what is even being measured in Figure 4C? The probability of the first mutation sweeping or the probability of any mutation sweeping? What is the x-axis precisely? Figure 4D is more clear wrt the x-axis, but it shows that the sweep probability increases as more mutations are allowed to accumulate, which is a more realistic scenario than the analytical model where mutations do not accumulate.\nFor the simulations that allow for multiplicative accumulation of mutations, the problems are potentially more severe.\nAccording to the process by which the proliferation rate of mutations changes (Eq. 17), for some of the values stated (for instance, s=0.1) the proliferation rate decreases on average (this can be seen by simulating the process described by Eq. 17). Since further mutations are relatively deleterious on average in this model, it is perhaps not so surprising that the fixation of mutant lineages is somewhat limited. Moreover, simulations with different fitness scales s are apparently pooled to generate Figs 4C-E and may have contributed to the \"anomaly\" described in the text and above.\nFinally, I will highlight a key prediction that the authors bring up, that is seemingly untrue. The statement “we predict… detectable tumours will almost always contain multiple competing clones with different driver mutations” (line 360) is central to the message of the paper (and I appreciate that the authors were bold enough to go out on a limb and say this). Unfortunately, it is likely false. The authors cite a convenient figure from Dentro et al. (that 95% of tumors in PCAWG exhibit subclonal expansion, line 429) but neglect the key observation from that same paper that only 11% of subclones carry a driver. This suggests that subclonal dynamics are driven by considerations other than purely fitness on a fixed landscape per se (for instance, there could perhaps be ‘ecological’ effects or a deformed fitness landscape). In an ecological scenario, for instance, sweeps would be prevented by more “trivial” mechanisms like disjoint niche preferences. Beyond that, it is clear that there is evidence for \"multi-hit\" origins of certain cancers, which involve the acquisition of multiple clonal mutations and therefore imply selective sweeps. If the authors were to provide some hard evidence from cancer data to back up some of the key predictions of the model, I would love to see it.\nMinor points:\nI appreciate that a version of the code has been uploaded to Github. It appears that much of the underlying simulation code resides in a single >4000 line c++ file (https://github.com/robjohnnoble/demon_model/blob/master/src/demon.cpp). Considering that previous versions of this code have ostensibly been reviewed elsewhere, reviewing this is beyond my scope. However, this file is somewhat difficult to navigate and could stand to be more thoroughly commented or broken up into more easily parsed chunks. Moreover, I was not able to get it to run.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Reply to reviewer 1\nThe authors have addressed all my comments, and I recommend the manuscript for\npublication.\nWe thank the reviewer again for helping us to improve our manuscript.\nReviewer #1 (Remarks to the Author):\nThe GitHub and Zenodo repositories are well structured and contain Readme files to\nguide readers. I appreciate that some integrations are verified and documented\nusing Mathematica. Reproducing agent-based simulations is not easy, as it requires\ngenerating and analyzing data using the demon package from other publications.\nFor readers who want to re-run the simulations, we recommend using the automated\ncomputational workflow ‘warlock’ rather than the underlying demon C++ code. We\nnow highlight this workflow and cite the relevant preprint (Bak et al., arXiv, 2023) at\nthe start of the relevant Methods section (lines 487-489). As explained in that\npreprint, “The warlock repository encapsulates demon into an automated and\nreproducible computational workflow to simplify parallel simulations and make the\nsoftware accessible to a wider community. This workflow is implemented in\nSnakemake, one of the most popular workflow management systems adopted by the\nbioinformatics community. In addition to local execution, we use an inbuilt\nmechanism to provide an interface for parallel High Performance Computing via\nSLURM workload manager. Hundreds of demon instances with varied parameter\nvalues can thus be executed in parallel with a single operation. All dependencies are\nresolved by Anaconda. The analysis of model output can also be automated simply\nby adding further components to the pipeline.”\nReply to reviewer 2\nStein et al have amended their manuscript in two primary ways, adding more\nsimulation results and expanding introductory and discussion points. I appreciate\nthat effort has been put in to address my points of concern. However, the\npresentation of some of the additional results (extending the simulation model to\nrandom and accumulative fitness affects) come across as confusing, and potentially\nincorrect. Moreover, the relationship of these new simulation results to the analytical\nmodel is somewhat tenuous. Finally, the key predictions of the model wrt its\nrelationship to cancer data seem to be demonstrably inconsistent with observations\ndescribed in the literature.\nTo address the reviewer’s concern about confusing readers, we have thoroughly\nrevised and simplified the presentation of the new simulation results. To clarify the\nrelationship between the new simulation results and the analytical model, we have\nreanalyzed the data and obtained new analytical results. To explain better how our\nresults are consistent with empirical cancer data, we have rewritten the relevant\nparts of the Discussion, in consultation with authors of one of the most pertinent\npapers.\nMajor points:\nI will first address the new individual-based simulations (those pertaining to Fig 4 in\nthe revised text). There are two new simulations that are introduced:\n1. Simulations in which random effect size mutations are introduced on the WT\nbackground, but further mutations cannot accumulate.\n2. Simulations in which mutations can continue to accumulate “multiplicatively.”\nThe simulations in which random effect single mutations are introduced are\nproblematic conceptually. The study aims to describe sweep probabilities, which are\ntypically infinite time quantities, but individual-based models can only be run for finite\ntime (presumably 1000 generations in this study, according to figure 4A, though this\nis never stated). The precise time to see a sweep (and therefore, the appropriate\ntime to run a simulation) will sensitively depend on several underlying parameters\nlike mutation rate, the mode of overall population growth (fixed, boundary-drive,\nexponential), and the distribution of fitness effects. Therefore, it is not clear that the\n“macroscopic” model and analytical results derived from it are the appropriate point\nof comparison for finite time data generated by individual-based simulations in the\ncase where mutation effect sizes vary. Indeed, if the simulations were run for longer,\nwe would eventually see an atypically large effect size beneficial mutation arise and\neventually sweep in every simulation by chance (assuming an exponential\ndistribution of effect sizes as the authors state).\nThe reviewer is correct that, of necessity, we ran simulations for only a finite period\nwhereas our analytical predictions are based on integrating over all times. As we\nnow state in the Methods (lines 529-530), simulations were stopped at a population\nsize of 1 million or at 2,000 generations (whichever was reached first). However, the\nreviewer’s intuition that an exceptionally fit mutant is eventually certain to sweep is\nincorrect.\nFirst consider the models in which we restrict the number of mutations that can\naccumulate per individual. Typically, all individuals reach this maximum mutation\nlimit within the simulation time span, so no mutations (and in particular no mutations\nwith atypically large effect size) could occur after the simulation stop time.\nIn the model with unrestricted accumulation of mutations, if we were to use strictly\nmultiplicative fitness effects (so that a mutation with effect s increases proliferation\nrate by a factor of 1 + s) then evolution would be unconstrained. Fitter mutants would\nexpand faster, proliferate faster, and so acquire ever more mutations in a positive\nfeedback loop. The model would thus, within finite time, generate infinitely fit\n“Darwinian demons”. To prevent this physically impossible scenario, we instead\nassume in our simulations that evolution is constrained by diminishing returns\nepistasis, such that the proliferation rate of mutants cannot exceed ten times that of\nthe wildtype. We describe this simulation feature in the Methods (lines 520-523;\nEquation 17). Due to diminishing returns epistasis, the fitness effect of new\nmutations generally diminishes as the population size grows.\nAlthough the above arguments refute the claim that “we would eventually see an\natypically large effect size beneficial mutation arise and eventually sweep in every\nsimulation”, the reviewer is correct that we might have observed more sweeps had\nwe run the simulations for longer. To estimate the extent of such undercounting, we\nhave conducted a new analysis of the population sizes at which detected sweeps\noccurred (Methods lines 532-539 and 556-573; new Figures S10-S14). This analysis\nshows that the distribution of population sizes is bell-shaped and is well\napproximated by a Gamma distribution. By fitting truncated Gamma distributions to\nthe simulation results, we inferred the truncation points and hence the proportion of\nsweeps missing from each set of simulations. We have adjusted the data points in\nFigure 5C (previously Figure 4C) to account for the inferred numbers of missing\nsweeps. As we report in the Methods, “The adjustments are negligible except in the\ncase of small fitness effects, for which sweep probabilities are in any case low. In\ntests using artificially truncated distributions from simulations with large s values\n(new Figure S15), we find that this method typically overestimates the proportion of\nmissing sweeps. The adjusted values can therefore be regarded as upper bounds.”\nOur main results and conclusions are unaffected.\nMoreover, the “anomaly” at low fitness observed in these simulations and described\nin Figure 4C is confusing, and I’m not sure the explanation provided in the text is\nsatisfactory. Specifically, what is even being measured in Figure 4C? The probability\nof the first mutation sweeping or the probability of any mutation sweeping? What is\nthe x-axis precisely? Figure 4D is more clear wrt the x-axis, but it shows that the\nsweep probability increases as more mutations are allowed to accumulate, which is\na more realistic scenario than the analytical model where mutations do not\naccumulate.\nTo make a meaningful comparison between the results of the fixed-effect and\nrandom-effects models, we need to use an appropriate summary statistic of the\ndistribution of mutation fitness effects. As we now state in the Results (lines 291-\n293), the problem is that “the random fitness effect determines not only a mutant's\nexpansion speed but also its probability of evading stochastic extinction” and hence\ndetermines the effective mutation rate. As we now state in the Methods (lines 547-\n550), “For example, in a Moran process with K = 16, nearly 94% of nearly-neutral\nmutations succumb to stochastic extinction before they can contribute to the\nmacroscopic evolutionary dynamics, compared to only 50% of mutations that confer\nfitness effect s = 1”. Hence, it is inappropriate to compare results simply in terms of\nthe mean fitness of all generated mutations. We chose the x-axes in the previous\nFigures 4C and 4D after considering several less satisfactory alternatives. In\nresponse to the reviewer’s comment we have revisited this problem and arrived at a\nbetter solution.\nAs now explained in the Methods (lines 550-553), “To account for this difference in\neffective mutation rates, we calculated the mean effect of mutations that reached an\nabundance of at least 10 individuals in the at-most-one-mutation model (see Figure\nS16 for the distributions).” Considering in this way only the “contending” mutations\nthat escape stochastic extinction is conventional (see review by Bataillon & Bailey\n2014) and is consistent with previous studies (e.g. Gerrish & Lenski 1998; Ralph &\nCoop 2010). As recorded in the Acknowledgements, we discussed this approach\nwith a leading expert on distributions of mutant fitness effects, Guillaume Martin of\nthe Institut des Sciences de l'Évolution de Montpellier, who agreed it is appropriate to\nour study. We have rewritten the relevant sections of the Results (lines 286-323) and\nMethods (lines 540-573) accordingly.\nWith this simpler summary statistic on the x-axis, we are confident that the new\nFigure 5C is easier to understand than the previous Figures 4C and 4D. As we report\nin the Results (lines 296-298), when compared in terms of this better summary\nstatistic, “sweep probabilities in random-fitness-effect simulations are close to the\npredictions obtained under the assumption of fixed effects”. For a given distribution\nof fitness effects, the sweep probability increases as more mutations are allowed to\naccumulate (in line with intuition).\nAlthough we believe that the new Figure 5C is the best way to present our results,\nwe also think it is useful to show that the similarity between sweep probabilities in\nrandom-fitness-effect and fixed-effect simulations holds for an alternative choice of x-\naxis. We have therefore retained the previous Figure 4C as Supplementary Figure\nS5B and added a variant of this plot (without pooling results across parameter sets)\nin Figure S5A. As now mentioned in the Results (lines 323-326), Figure S5 also\nshows that “sweeps are more often due to a single large-effect mutation than several\nmutations with small effects”. We have rewritten the figure caption to answer the\nreviewer’s questions regarding the axes.\nWe have also retained Section 9 of the SI Text, which discusses the “anomaly” at\nlow fitness in Figure S5B. We believe SI Text Section 9 is useful because it sheds\nlight on the subtleties inherent to summarising and interpreting results when mutation\neffects are random; because it describes an effect that other researchers may\nobserve when analysing similar results; and because the transition matrix analysis\nlays the groundwork for further mathematical investigations. We have added a two-\nparagraph introduction to this section of the SI Text (lines 347-360).\nFor the simulations that allow for multiplicative accumulation of mutations, the\nproblems are potentially more severe. According to the process by which the\nproliferation rate of mutations changes (Eq. 17), for some of the values stated (for\ninstance, s=0.1) the proliferation rate decreases on average (this can be seen by\nsimulating the process described by Eq. 17). Since further mutations are relatively\ndeleterious on average in this model, it is perhaps not so surprising that the fixation\nof mutant lineages is somewhat limited.\nWe are grateful to the reviewer for spotting a typo in Eq. 17. A closing parenthesis\nwas misplaced so that the equation in the text differed from that in the simulation\ncode.\nWe wrote\nIt should be\nIn the correct equation, all mutations are beneficial. Importantly, all simulations were\nbased on the correct equation (as can be verified by reviewing the history of the\ncode in the GitHub repository).\nMoreover, simulations with different fitness scales s are apparently pooled to\ngenerate Figs 4C-E and may have contributed to the \"anomaly\" described in the text\nand above.\nAs we now explain in the SI Text (lines 359-360), the spike anomaly “is not simply\ndue to the pooling of results because it also appears in single-batch results”,\nas\nshown in new Figure S5C. We have revised the caption of the former Fig 4C (now\nFig S5B) to clarify that it shows pooled results.\nIn the new Figures 5C, 5D and S5A, each data point corresponds to one batch of\nsimulations with identical parameter values. There is no pooling of results and no\nspike anomaly at low fitness values.\nFinally, I will highlight a key prediction that the authors bring up, that is seemingly\nuntrue. The statement “we predict…detectable tumours will almost always contain\nmultiple competing clones with different driver mutations” (line 360) is central to the\nmessage of the paper (and I appreciate that the authors were bold enough to go out\non a limb and say this). Unfortunately, it is likely false. The authors cite a convenient\nfigure from Dentro et al. (that 95% of tumors in PCAWG exhibit subclonal expansion,\nline 429) but neglect the key observation from that same paper that only 11% of\nsubclones carry a driver. This suggests that subclonal dynamics are driven by\nconsiderations other than purely fitness on a fixed landscape per se (for instance,\nthere could perhaps be ‘ecological’ effects or a deformed fitness landscape). In an\necological scenario, for instance, sweeps would be prevented by more “trivial”\nmechanisms like disjoint niche preferences. Beyond that, it is clear that there is\nevidence for \"multi-hit\" origins of certain cancers, which involve the acquisition of\nmultiple clonal mutations and therefore imply selective sweeps. If the authors were to\nprovide some hard evidence from cancer data to back up some of the key\npredictions of the model, I would love to see it.\nWe appreciate that the quoted final sentence of the final section of our Results could\nreasonably be interpreted as going out on a limb to propose an unconventional\nhypothesis. This was not our intention. We have therefore rephrased the sentence\n(lines 357-360) to provide a more accurate and precise synopsis: “In summary, our\nmacroscopic model (assuming equal mutation effects and no accumulation of\nmutations) predicts that selective sweeps during a clonal expansion are rare unless\nmutations are very strongly beneficial, in which case sweeps begin and end early in\nthe expansion.”\nWe fully agree that “subclonal dynamics are driven by considerations other than\npurely fitness”. As we now clarify at the start of the Results (lines 82-84), “Our\nmacroscopic model is designed to test whether clonal interference alone can prevent\nselective sweeps, and to obtain upper bounds on selective sweep probabilities\nduring range expansions.” In a new Discussion paragraph (lines 424-440), we\nexplain the ways in which microenvironmental heterogeneity is likely to thwart\nselective sweeps. We also give other reasons why the general conclusion that\n“selective sweeps during clonal expansions are restricted to strong driver mutations\nduring the early stages of expansion” might be expected to hold even when driver\nmutation effects are random and multiplicative. We also now explain (on lines 473-\n476) that, “To obtain more precise predictions rather than upper bounds, it will be\nimportant to examine how microenvironmental heterogeneity, immune responses,\nand phenotypic plasticity inhibit clonal expansions in each cancer type.”\nAs the reviewer suggests that we might have misinterpreted or misrepresented the\nfindings of Dentro et al. (Cell, 2021), we sought advice from two of the co-first\nauthors of that paper: Jonas Demeulemeester and Maxime Tarabichi (who are also\nthe postdoc co-supervisors of Alex Stein, the first author of our paper). As now\nmentioned in the Acknowledgments, these two researchers read our manuscript,\nanswered our questions about their study, and offered suggestions to help us\nrephrase relevant passages. As the reviewer points out, not all clonal expansions in\nDentro et al. can be directly attributed to driver mutations. However, the 11%\nmentioned by the reviewer is a conservative lower bound for the following reasons:\n(i) The 11% refers only to known driver mutations. As shown by Gerstung et al.\n(Nature, 2020) and as also argued by Dentro et al., driver mutations that are typically\nfound to be subclonal arise from a more diverse pool of genes than driver mutations\nthat are typically found to be clonal. This suggests we have likely not yet identified all\npossible driver mutations. (ii) The 11% refers to SNVs and indels only. The analysis\nexcluded CNAs and SVs that can also act as driver mutations. (iii) Positive dN/dS\nratios suggest widespread positive selection of ‘functional’ mutations, further\nindicating the involvement of genetic mutations in clonal expansions. The complete\nquote from Dentro et al. acknowledges the first two limitations: “Our observations\nhighlight a gap in knowledge about drivers of subclonal expansions; only 11% of all\nsubclones carry a known SNV or indel driver mutation. Late tumor development\ncould be driven largely by different mechanisms (such as CNAs and SVs [Jamal-\nHanjani et al., 2017; Mamlouk et al., 2017] or epigenetic alterations), or the pool of\nlate driver mutations is too large to have been sampled extensively enough so far.”\nRegarding the third factor: “Analysis of dN/dS ratios revealed positive selection\nacross subclones and cancer types. Although our analyses do not exclude the\npossibility that a fraction of tumors evolve under weak or no subclonal selection, they\nsupport widespread selection.” We have revised and expanded the relevant parts of\nthe Discussion (lines 441-452) to clarify these points.\nOur results are entirely consistent with “multi-hit” cancer origins. In the Discussion,\nwe now refer explicitly to “multi-stage models” and explain how recent experimental\nstudies have tried to estimate the time of driver events with partial success. To\nlessen the risk of misunderstanding, we have revised the part of the Discussion that\nexplains how clonal driver mutations can arise through bottlenecks (lines 453-473).\nIn summary: With regard to cancer, our methods are designed to establish an upper\nbound on the selective sweep probability during a clonal expansion. Since we want\nan upper bound, we omit from the model factors that would impede selective\nsweeps, such as microenvironmental heterogeneity and accelerating propagation\nspeed. As we find that the upper bound is low, we conclude that multiple clonal driver\nmutations observed in tumours most likely arise not during the final clonal expansion\nbut instead through population bottlenecks in a stepwise expansion process. We are\nnot proposing an unorthodox hypothesis but rather determining which previously\nproposed explanation is the most likely, thus helping to resolve a debate in the field\n(as discussed in our Introduction, lines 53-70, and Discussion, lines 453-457).\nMinor points:\nI appreciate that a version of the code has been uploaded to Github. It appears that\nmuch of the underlying simulation code resides in a single >4000 line c++ file\n(https://github.com/robjohnnoble/demon_model/blob/master/src/demon.cpp).\nConsidering that previous versions of this code have ostensibly been reviewed\nelsewhere, reviewing this is beyond my scope. However, this file is somewhat\ndifficult to navigate and could stand to be more thoroughly commented or broken up\ninto more easily parsed chunks. Moreover, I was not able to get it to run.\nAs the reviewer points out, the simulation framework ‘demon’ was central to previous\npeer-reviewed papers (Noble & Burley et al., Evolutionary Applications, 2020; Noble\net al., Nature Ecology & Evolution, 2021) and has therefore been scrutinized before.\nWe appreciate that working with C++ code can be challenging, which is why we\ncreated the automated computational workflow ‘warlock’ to make our software more\naccessible. For readers who want to rerun the simulations, we recommend using\nwarlock rather than the underlying demon C++ code. We now mention this workflow\nat the start of the relevant Methods section (lines 487-489) and cite a preprint that\ndescribes how it works and how to use it (Bak et al., arXiv, 2023).\nWithout details of the reviewer’s software environment and the error messages they\nencountered, we are unable to comment on why they were unable to run our code.\nWe and others have successfully run the code on Windows, macOS and Unix\noperating systems on personal computers and on the computing clusters of three\nuniversities."}, {"round": 3, "reviewer_comments": {"reviewer_4": "As I have been brought in relatively late in the review process, I focused my attention on independently assessing how the previous reviewers' comments have been addressed. My assessment is that the authors have carefully resolved the confusion that arose when introducing new results from agent based simulations. They have amended and clarified some potentially problematic statements that were pointed out by the previous reviewers and addressed the concern about the apparent disagreement between some cancer data and their results, specifically in regards to the reference by Dentro et al.\nIn conclusion, I think all the points carefully raised by the previous reviewers have been fully addressed and the manuscript is ready for publication.\nThe available code has also been carefully amended according to the reviewers' suggestions."}, "author_response": ""}]

Decision letter