Research Article

Punctuated mutagenesis promotes multi-step evolutionary adaptation in human cancers

Department of Data Science, Dana-Farber Cancer Institute, United States
Department of Biostatistics, Harvard T.H. Chan School of Public Health, United States
Department of Stem Cell and Regenerative Biology, Harvard University, United States
School of Computer Science, Carnegie Mellon University, United States
Department of Pathology, New York University School of Medicine, United States
Laura and Isaac Perlmutter Cancer Center, New York University, United States
Center for Cancer Evolution, Dana-Farber Cancer Institute, United States
The Eli and Edythe L. Broad Institute, United States
Ludwig Center at Harvard, Harvard Medical School, United States

Oct 9, 2026

https://doi.org/10.7554/eLife.108058.3

Open access
Copyright information

eLife Assessment

This valuable study presents a theoretical model of how punctuated mutations influence multistep adaptation, supported by empirical evidence from some TCGA cancer cohorts. This solid model points to the case of possible punctuated evolution rather than gradual genomic change. There was some disagreement amongst the reviewers in terms of how closely the theoretical results apply to the phenomena examined empirically, and alternative explanations should be considered in the future.

https://doi.org/10.7554/eLife.108058.3.sa0

Significance of the findings:

Valuable: Findings that have theoretical or practical implications for a subfield

Landmark
Fundamental
Important
Valuable
Useful

Strength of evidence:

Solid: Methods, data and analyses broadly support the claims with only minor weaknesses

Exceptional
Compelling
Convincing
Solid
Incomplete
Inadequate

During the peer-review process the editor and reviewers write an eLife Assessment that summarises the significance of the findings reported in the article (on a scale ranging from landmark to useful) and the strength of the evidence (on a scale ranging from exceptional to inadequate). Learn more about eLife Assessments

Abstract
Introduction
Results
Discussion
Methods
Appendix 1
Appendix 2
Appendix 3
Appendix 4
Appendix 5
Data availability
References
Article and author information
Metrics

Abstract

The rate of acquisition of genomic changes in cancer has been the topic of much discussion, with several recent investigations finding evidence of punctuated evolution instead of gradual accumulation of such changes. Despite forays into the description and quantification of these punctuated events, the effects of such changes on subsequent cancer evolution remain incompletely understood. Here we investigate how non-gradual mutagenesis affects the ability of tumor cells to acquire and retain fitness-enhancing adaptations. We find that punctuated mutagenesis significantly facilitates adaptation in scenarios where adaptation requires crossing a fitness valley, that is when multiple mutations are required which individually are maladaptive but jointly confer a fitness advantage. By increasing the probability that multiple mutations occur in close succession, punctuation increases the chance that mutants in a fitness valley mutate further to reach a fitness peak before going extinct. Analyzing data from The Cancer Genome Atlas, we find that tumors with signatures of APOBEC mutagenesis, which has been shown to proceed in episodic bursts, exhibit patterns consistent with higher rates of crossing fitness valleys. Lastly, we characterize how the interplay between this enhanced ability to cross fitness valleys and adaptation-limiting effects of clonal interference affects overall adaptability in complex fitness landscapes.

Introduction

Tumors evolve via the acquisition of randomly arising genetic and/or epigenetic alterations and their selection (Merlo et al., 2006). The rate at which tumor cells accumulate such genomic changes to produce potentially adaptive innovation is thus an important determinant of the capacity of tumor cells to adapt to diverse selection pressures (Yates and Campbell, 2012; Neinavaie et al., 2021), affecting their propensity to progress locally, to metastasize, or to evolve resistance to therapeutic intervention. Several methods to quantify mutation rates have been developed (Araten et al., 2005; Williams et al., 2018; Werner et al., 2020), and proxies for mutation rates such as tumor mutation burden often feature in prediction models of therapeutic outcomes (Galuppini et al., 2019; Sha et al., 2020).

Most existing methods deal with constant mutation rates (Werner et al., 2020; Altrock et al., 2015), assuming gradual evolutionary change over time. However, accumulating evidence suggests that mutagenesis in cancer is fluctuating. Indeed, a recent study Petljak et al., 2019 demonstrated in vitro that mutations associated with DNA-editing activity of APOBEC cytidine deaminases (Petljak et al., 2022) occur in episodic bursts, with more than 100-fold differences in the rate of APOBEC-associated mutations across otherwise identical cell culture replicates (Figure 1A). This episodicity observed in vitro aligned with patterns in previously published in vivo data investigating APOBEC mutagenesis in lung cancer (Jamal-Hanjani et al., 2017; Lee et al., 2017), and data from intestinal crypts (Wang et al., 2023) later confirmed explicitly that this episodic pattern also occurs in patients. Other recent studies using genomic data to reconstruct evolutionary histories of tumors have found evidence of punctuated evolution across several cancer types (Navin et al., 2011; Sottoriva et al., 2015; Gao et al., 2016). These studies have demonstrated the existence of distinct phases of mutation bursts (Gao et al., 2016; Minussi et al., 2021; Figure 1B), challenging the prevailing paradigm of gradual emergence of mutant lineages (Figure 1C). Various mutagenic processes such as chromothripsis (Stephens et al., 2011; Zhang et al., 2015), chromoplexy (Baca et al., 2013), and other drivers of chromosomal instability Drews et al., 2022 have been implicated in causing such punctuated patterns (Davis et al., 2017).

Figure 1

Download asset Open asset

Population dynamics under gradual vs punctuated mutagenesis.

(A) Evidence for punctuated APOBEC mutagenesis (Petljak et al., 2019) from successive cell culture expansions, seeded with single progenitors from the preceding expansion. Sequencing of the expanded populations revealed large fluctuations in APOBEC-associated mutagenic signatures SBS2 and SBS13. (B) Evidence for punctuated copy number evolution (Gao et al., 2016; Minussi et al., 2021). Patterns in branch lengths of reconstructed phylogenetic trees reveal fluctuations in the rates of copy number alterations. Dynamics prior to the most recent common ancestor (MRCA) are unidentifiable. (C) Gradual vs punctuated evolution. In gradual evolution, novelty emerges and spreads at a constant rate over time. In punctuated evolution, novelty emerges and spreads during distinct burst phases. (D) Fitness schematic of two one-step adaptations which each independently confer a fitness advantage. (E) Fitness schematic for a two-step adaptation, in which carrying one mutation confers a fitness disadvantage, but carrying two mutations confers a fitness advantage. Schematic of possible evolutionary dynamics for two one-step adaptations (F) and of modes of valley crossing with or without prior fixation of the first mutant (G). Sequential fixation occurs at low mutation rates where emerging mutant lineages are likely to have fixated or gone extinct before the next mutation occurs. At higher mutation rates, the second mutation can occur in a multi-clonal population. (H) Sketch of failing two-step adaptation under a uniform (time-invariant) mutation rate. (I) Sketch of a successful two-step adaptation via stochastic tunneling under a punctuated mutation process with distinct clusters of high mutation rates.

How punctuated mutagenesis affects the evolutionary dynamics of tumor cell adaptation across complex fitness landscapes remains incompletely understood. Addressing this question with existing experimental data is inconclusive, as only short time horizons are observed (single expansions [Petljak et al., 2019]; time to a most recent common ancestor [Gao et al., 2016; Minussi et al., 2021]), and as fitness differences between observed cell lineages are incompletely characterized. Mathematical modeling of human cancer genomics data, however, can elucidate the dynamics of adaptation during punctuated tumor evolution. Here, we set out to systematically investigate the evolutionary consequences of punctuated mutation acquisition during tumorigenesis using mathematical modeling and analysis of genomic data from The Cancer Genome Atlas (TCGA).

Results

Temporal clustering facilitates multi-step adaptation

We set out to investigate how temporal clustering of mutagenic events into distinct episodes affects the ability of a cell population to achieve multi-step evolutionary adaptations. We considered two scenarios: tumor evolution via a single advantageous step, such as an oncogenic adaptation (Chakravarty et al., 2017) (‘one-step adaptations’, Figure 1D), and accumulation of multiple mutations with synergistic fitness effects (‘multi-step adaptations’, Figure 1E). Indeed, the two hit-hypothesis (Knudson, 1971) for tumor suppressor genes (TGSs) suggests that inactivation of one copy of a TSG can be inconsequential or maladaptive, while bi-allelic inactivation confers a selective advantage (Iwasa et al., 2004). While some TSGs exhibit (context-dependent) haploinsufficiency (Paige, 2003; Wang et al., 2018; Park et al., 2021), synergistic fitness effects for successive (epigenetic; Paige, 2003; Wang et al., 2018; Issa, 2022) alterations in the same gene remain the prominent feature of TSGs. Analogously, there are ample examples of successive alterations across different genes acting synergistically (Gu et al., 2010; Hsu et al., 2010; Deming et al., 2014). High mutation rates are known to limit the rate of retained one-step adaptations per emerging one-step adaptation due to clonal interference (Gerrish and Lenski, 1998; Wilke, 2004; Fogle et al., 2008; Figure 1F), which has been highlighted as a potential consequence of punctuated cancer evolution (Sun et al., 2018). The dynamics of multi-step adaptation have been studied extensively (Iwasa et al., 2004; Komarova et al., 2003; Nowak et al., 2004; Weinreich and Chao, 2005; Weissman et al., 2009; Proulx, 2011; Ashcroft et al., 2015; Van Egeren et al., 2018), elucidating two modes of evolution (Figure 1G): ‘sequential fixation’ refers to the scenario in which the second mutation only emerges after cells harboring the first mutant have taken over the population, while ‘stochastic tunneling’ refers to situations in which the second mutation arises sooner. In large populations, mutants with a selective disadvantage become vanishingly unlikely to reach fixation, so that stochastic tunneling becomes the dominant mode of evolution for multi-step adaptations.

To investigate how the rate of successful stochastic tunneling events – enabling multi-step adaptation – is affected by punctuated mutagenesis, we simulated selection dynamics of a population of cells which proliferate according to a Wright–Fisher process (Wright, 1931; Fisher, 1923) and accumulate mutations according to a constant (Figure 1H) or temporally clustered (Figure 1I) rate. The Wright–Fisher process models evolution as successive, non-overlapping generations of constant size. Each new generation is populated by drawing with replacement from the cells in the previous generation, with probabilities proportional to their fitness (Figure 2A). Cells can accumulate mutations that change their fitness, tha is increase their probability to survive to the next generation. We considered fitness landscapes in which two-step adaptations are the only way by which cells can increase their fitness, with fixed fitness ratios across subsequent two-step adaptations (Figure 2B, C). We also confirmed that the findings observed for the Wright–Fisher process, described below, are qualitatively consistent for expanding populations that proliferate according to a branching process model (Methods, Figure 2—figure supplement 1).

Figure 2 with 4 supplements see all

Download asset Open asset

Simulation results: valley crossing under uniform vs temporally clustered mutation rates.

(A) Schematic of a Wright–Fisher process. (**B, C**) Fitness landscapes used for the simulations in panels (D) and (E), respectively. Mutations move a cell from left to right through the landscape. Each two-step adaptation corresponds to a fitness increase by a factor of 1.5. Having an odd number of mutations comes at a multiplicative fitness disadvantage of 0.5 in (B) and an advantage of 1.01 in panel (D). (**D, E**) Heatmaps of simulation results for a Wright–Fisher process with 50 cells. The mutation rate trajectories in each panel are chosen such that the total expected number of mutations under the uniform trajectory is identical to that under the temporally clustered trajectory. Each column of pixels represents a population after a specific number of generations. Homogeneity in pixel-color per column thus indicates that cells with a given number of adaptations have fixated in the population, whereas persistent diversity of pixel colors (see Figure 2—figure supplement 1) indicates coexistence of lineages with different numbers of adaptations.

Interestingly, for two-step adaptations for which the intermediate mutants that carry only one of the two required mutations have a strong selective disadvantage (Figure 2B), we found that the population undergoes substantially (3.55 times) more two-step adaptations when accumulating mutations in a temporally clustered rather than in a uniform way (Figure 2D). This effect arises because intermediate mutants in ‘fitness valleys’ have a high chance of going extinct before acquiring the next mutation. During mutation bursts, acquiring the next mutation in time before the disadvantageous mutant goes extinct becomes more likely. Similarly, intermediate mutants are more likely to emerge in a mutation burst. This temporal clustering of the likelihood of two succeeding mutation events increases the rate of successful two-step adaptations.

This effect also emerges for two-step adaptations that do not constitute proper fitness valleys, that is for which the intermediate mutant is not maladaptive (Figure 2E; 1.15-fold increase). Under stochastic selection dynamics, even mutants with a slight selective advantage (Figure 2C) have a high chance of going extinct due to random drift. Lineages acquiring sets of synergistic mutations, thus, often do so without prior fixation of each intermediate mutant – via stochastic tunneling (Iwasa et al., 2004). In those regimes, temporal clustering of mutations therefore also increases the speed of adaptation, as confirmed with simulations (Figure 2—figure supplement 2A–D). Moreover, the proportion of such sets of synergistic mutations acquired via stochastic tunneling rather than sequential fixation also increases with temporal clustering (Figure 2—figure supplement 2E).

Our results generalize to multi-step adaptations with arbitrarily many steps (Appendix 1). To demonstrate this result, we compared a uniform mutation process with mutation rate µ to temporally clustered mutation processes in which mutations emerge exclusively during burst phases that constitute a fraction $\frac{1}{k}$ of the time, but at a $k$ -fold increased rate, $k μ$ (Figure 2—figure supplement 2A). We showed analytically that in the limit of the average mutation rate µ going to zero, the rate $f_{k}$ of successful $(n + 1)$ -step adaptations is $k^{n}$ times higher in the temporally clustered mutation process than in the uniform mutation process. We validated this result with simulations of a Wright–Fisher process for two-step adaptations, $n = 1$ , and a range of small values for µ (Figure 2—figure supplement 2F, G), and with numerical calculations of $f_{k}$ in a branching process (see Appendices 2 and 3 and Figure 2—figure supplements 3 and 4 for complete characterizations of parameter dependencies). As the mutation rate increases, this fold-increase $\frac{f_{k}}{f_{1}}$ becomes smaller. However, since the absolute rates f_k and f₁ are proportional to µⁿ⁺¹ (Appendix 1), the absolute effect of temporal clustering on the rate of fixations of mutants with multi-step adaptations increases in µ (Figure 2—figure supplement 2F). Analogously, we found that $\frac{f_{k}}{f_{1}}$ increases with the selective disadvantage that cells in a fitness valley have, while $f_{k} - f_{1}$ decreases with this selective disadvantage (Appendix 2). In contrast, for low µ, and particularly for parameter choices in biologically relevant ranges ( $μ \leq 10^{- 7}$ ), the fold-change becomes insensitive to the magnitude of the selective disadvantage so that ever smaller fitness differences suffice to reach high fold-changes, $\frac{f_{k}}{f_{1}} \approx k^{n}$ (Appendix 2). For instance, in a population of constant size that evolves according to a branching process with a birth–death ratio of 1, if intermediate mutations in a multi-step adaptation sequence reduce the birth–death ratio to 0.98, we have $\frac{f_{50}}{f_{1}} = 49.4 = 0.988 k$ for two-step adaptations, and $\frac{f_{50}}{f_{1}} = 2440.6 = 0.976 k^{2}$ for three-step adaptations. In light of the large fold-changes in mutation rates during a burst and before or after a burst (Petljak et al., 2019; Navin et al., 2011; Gao et al., 2016; Minussi et al., 2021), our findings suggest that the resulting effect on the dynamics of multi-step adaptation is substantial.

Temporal clustering in an exploration–exploitation setting

The increased rate of valley crossing is driven by phases of high mutation rates, when the fold-increase in the chance of valley crossing is larger than the fold-increase in the mutation rate. This effect improves adaptability while the fitness landscape offers scope for multi-step adaptation. However, in fitness landscapes which have global maxima, as the population adapts there is ever less scope for (further) multi-step adaptations, and excessive exploration of the landscape becomes costly.

To investigate the effect of temporal clustering in such an exploration-exploitation setting, we constructed a two-dimensional fitness landscape (Figure 3A) which is randomly re-drawn at regular intervals, mimicking exposure to novel physical environments or drugs during cancer evolution and treatment (Methods). Between subsequent re-drawing, the population may move to a global fitness maximum where any further mutation decreases fitness. We considered mutation processes with different average mutation rates µ and different values of the clustering parameter k (Figure 3B) and investigated the average fitness of the cell population over long time horizons. As before, the clustering parameter k denotes the fold-increase of the mutation rate in burst phases relative to the average mutation rate µ. However, rather than fixing the out-of-burst mutation rate at zero, the burst duration is now held constant (Figure 3B) to control for differences in waiting time until a burst occurs.

Figure 3

Download asset Open asset

Exploration and exploitation with temporally clustered mutation rates.

(A) Cells move through two-dimensional fitness landscapes. These fitness landscapes are randomly re-drawn every 50,714 divisions. (B) Mutation rate trajectories are parameterized with a mean mutation rate μ, and a clustering parameter k. (C) Average fitness in simulations of a population of 20 cells in a Wright–Fisher process. Simulations were run for 10⁸ division events. (**D–F**) Average fitness trajectories for representative snippets of the simulations in (C).

We found that, in this scenario, relative to the uniform mutation process (k = 1) increasing k initially increases the average fitness of the population for any fixed µ (Figure 3C). However, at high $k$ , marginal increases of $k$ reduce fitness, since high in-burst mutation rates in small populations increase the chance of leaving fitness peaks due to stochastic drift, and low mutation rates outside of bursts delay the acquisition of one-step adaptations. Similarly, for a fixed k, increasing µ when starting at low µ-values increases the average fitness. However, moving further in either of these directions in the parameter space – toward increasing µ or toward increasing k – eventually brings about a decrease in average fitness reflecting a trade-off between exploration and exploitation (see Appendix 4 and Figure 3D–F for detailed descriptions and plots of the dynamics). For any fixed k, average fitness peaks at intermediate µ, and the highest average fitness on the k–µ plane is achieved by a k > 1. Maximizing average fitness, thus, involves temporal clustering even when there is a trade-off between exploration and exploitation.

Proxies for valley crossing and for temporal clustering found in patient data

We then sought evidence of the modeling prediction that temporal clustering facilitates multi-step adaptation in patient data, leveraging whole-exome sequencing (WES) data from TCGA. We reasoned that among tumors in which similar numbers of mutations had emerged, those tumors with more temporally clustered mutation processes in their evolutionary history would be more likely to have undergone a larger number of multi-step adaptations. Generalizing this idea to enable comparison of tumors which may differ in the number of mutations they have acquired, we reasoned further that tumors with high levels of temporal clustering would have a high fraction of successful multi-step adaptations per attempted multi-step adaptation, that is per emergence of a first mutant in a multi-step adaptation sequence. We would, thus, expect analogous ratios of detectable quantities that scale with the numbers of successful or attempted multi-step adaptations to correlate with indicators of fluctuations in the mutation rate history of a tumor such as those found for APOBEC mutagenesis. In particular, we considered the frequency of observing two deactivating mutations in a TSG as an indicator for successful multi-step adaptation, and the total number of observed mutations in a tumor as a quantity that scales with attempted multi-step adaptations.

Before considering the TCGA data, to confirm in silico that such a correlation would indeed emerge under biologically realistic parameters, we performed large-scale simulations of tumor growth from a single cell to up to realistically detectable cell numbers of 10⁶–10⁷ cells (Fischer et al., 2006) (Methods). We recorded the fraction of mutations that emerged during burst phases, constructed a TSG deactivation score (Figure 4A, B) by counting the number of TSGs (modeled as mutations with synergistic fitness effects) with at least two mutations across the population and divided by the total number of mutations across the genome (Methods). Our simulation results confirmed that these two readouts are correlated (Pearson correlation 0.42, p < 0.01; Figure 4C), showing that in silico predictions regarding the likelihood of acquiring two-step adaptations under mutation processes with more vs less temporal clustering are well reflected in our score for TSG deactivation.

Figure 4 with 4 supplements see all

Download asset Open asset

Proxies for valley crossing and temporal clustering in simulations and in TCGA data.

(A) Sketch of the simulation analysis workflow. (B) Distributions of the fraction of mutations acquired during mutation bursts and the TSG deactivation score in simulations. (C) Joint distribution of the quantities in (B) indicates strong correlation. (D) Schematic of analysis workflow for TCGA data (methods) and results for the four cancer types with highest APOBEC signature contribution. (E) Probabilities that single base substitutions (SBSs) in samples of the cancer types in (D) were caused by APOBEC-associated signatures SBS2 and SBS13. On the left, these probabilities are averaged over all SBSs. On the right, only those SBSs are considered which are classified as ‘nonsense’ or ‘missense’, and which appear in a TSG with at least two such deactivating SBSs. The average probabilities between both groups are significantly different (t-test, p < 0.001). (F) Average probabilities analogous to those in panel (D) were constructed for each individual sample and then averaged across samples. Samples without deactivated TSGs were excluded. The average probabilities between groups are significantly different for the four cancer types (*, **, and ***, respectively, indicate that the p-value in a t-test lies below 0.05, 0.01, and 0.001).

Next, to construct similar readouts from the TCGA data, we first performed SBS signature decomposition on each sample and computed the relative contribution of APOBEC-associated mutation signatures (SBS2 and SBS13) vs all other mutation signatures (Alexandrov et al., 2013) (Methods). Given the evidence for episodic APOBEC mutagenesis (Petljak et al., 2019; Wang et al., 2023), we used this relative contribution as a score for temporal clustering. As a proxy for multi-step adaptations, a list of known TSGs (the TSGene 2.0 database; Zhao et al., 2016) was considered, and for each TCGA sample, we counted the number of TSGs in this list which harbored at least two single-nucleotide substitutions classified as a missense or nonsense mutation (Methods) normalized by the total number of single-nucleotide substitutions in that sample. The TCGA WES data did not enable us to determine whether such mutations indeed deactivated different alleles of a TSG, rather than appearing on the same allele. If the likelihood of appearing on the same allele varied with the contribution of the APOBEC signature, this association could confound our results. Since mutations that are interdependently generated on the same copy of a gene would likely appear in close spatial proximity (Petljak et al., 2022; Nik-Zainal et al., 2012; Mas-Ponte and Supek, 2020), we thus verified in a robustness check that our results remain unchanged when filtering out mutations with close-by genomic coordinates (Methods).

We then ranked the different tumor types in TCGA by the average contribution of SBS2 and SBS13 relative to all detected mutation signatures, and for each tumor type, investigated the correlation between our proxies for multi-step adaptations and for temporal clustering. For the top four categories with the highest mean contribution of APOBEC-driven mutagenesis to all mutations, we observed a positive correlation between the two proxies. This correlation was significant at the 5% level (t-test on the Pearson correlation coefficient) for three out of four of these top four categories (Figure 4D, for results for categories with lower APOBEC signature contribution see Figure 4—figure supplement 1).

Consistent with the hypothesis that this correlation arises because TSG deactivation occurs more readily with APOBEC mutagenesis, we found that mutations which contribute to the TSG deactivation score on average have a strongly increased likelihood of being caused by APOBEC (Methods, Figure 4E). This finding is consistent across cancer types and independent of whether probabilities are averaged over samples or over individual mutations, suggesting that this result is not driven by outlier samples with many TSG inactivations and many APOBEC-associated mutations (Methods, Figure 4F).

An alternative score for multi-step adaptations, in which the list of TSGs is replaced with a list of pairs of genes with synergistic fitness effects (Gu et al., 2010) (Methods), again exhibits a positive correlation with the relative contribution of APOBEC-associated mutation signatures in all categories; this correlation is significant for two of these categories (Figure 4—figure supplement 2). Additionally, we verified that our results are robust to the choice of correlation metric: using Spearman instead of Pearson correlations yields significant positive correlations across all top four categories discussed above (Figure 4D, Figure 4—figure supplement 3).

Taken together, our analyses of TCGA data show that tumors with a larger share of APOBEC-associated mutations are enriched for pairs of mutations with synergistic fitness effects. While we cannot exclude the influence of unobserved confounders, given the limitations of the available data, we found that the patterns from our analysis are consistent with the theoretical predictions of our model.

Stochastic tunneling vs clonal interference during punctuated evolution

In tumors evolving at high mutation rates, clonal interference may limit the rate at which a cell population manages to acquire and retain one-step adaptations (Figure 1F). Previous literature has quantified this limiting effect as a function of the mutation rate and the distribution of fitness values of emerging mutants (Campos et al., 2004; Campos and de Oliveira, 2004; Park and Krug, 2007). However, the effects of non-constant mutation rates in this setting have not been elucidated. We thus set out to investigate clonal interference in the setting of temporally clustered mutation rates (Figure 5A).

Figure 5 with 1 supplement see all

Download asset Open asset

Stochastic tunneling vs clonal interference under temporally clustered mutation rates.

(A) Schematic of the simulation approach. (B) Sketch of the mutation process. (C) Simulation results: measuring fixations of higher-fitness mutations per generation as a function of the clustering parameter and of the fitness landscape. In the one-step adaptation setting, fitness is defined as ${1.5}^{# m u t a t i o n s}$ . Fitness in the two-step adaptation setting is defined as in Figure 2B. Results are averaged over simulation runs with 10⁷ generations. (D) Sketch of fitness landscape for simulations shown in (E). Cells start with an unmutated genome of 200 loci. Mutations on each locus have independent multiplicative effects on cell fitness. In the first 100 loci, a single mutation confers a multiplicative fitness change of 1.5. In the remaining 100, two mutations are required to reach this multiplicative fitness change of 1.5, with the first mutation conferring a multiplicative fitness change of 0.5. (E) Simulation results for adaptation with genome sketched in (D). Results are averaged over 100 simulation runs per value of k. Shaded regions indicate 10th and 90th percentiles. Smaller plots (left) are zoomed-in versions of the first 1200 generations of the larger plots (right), with identical color-coding.

We first considered a scenario in which cells can only acquire one-step adaptations and measured the rate at which novel one-step adaptations reach fixation in the population across simulations with different clustering parameters k (Figure 5B). We found that fixation rates quickly decrease as k increases (Figure 5C). This pattern emerges because the extent to which clonal interference reduces fixation rates disproportionately increases with the mutation rate; for a fixed average mutation rate, temporal clustering thus increases the effects of clonal interference.

These findings stand in contrast to our results for two-step adaptations whose rate of fixation increases when temporal clustering is introduced. However, if in-burst mutation rates are sufficiently large such that multiple clones in the population may acquire a two-step adaptation independently before one of these clones has fixated, effects of clonal interference also play a role for two-step adaptations. Performing simulations for two-step adaptations, we found that fixation rates are non-monotone in k. While at low k increasing k leads to a steep increase in the fixation rate, this trend eventually levels off and becomes negative, with further increases in k leading to a decrease in the fixation rate (Figure 5C).

Having observed that temporal clustering increases effects of clonal interference but also facilitates stochastic tunneling, we next set out to investigate the relative contribution of these two effects on overall adaptation rates. We considered a setting in which cells can increase their fitness through both one- and two-step adaptations (Figure 5D). We found that one-step adaptations are acquired more quickly under the uniform mutation process (k = 1), while two-step adaptations are acquired more quickly under a temporally clustered mutation process (k = 5; Figure 5E). The relative magnitude of both effects varies with time. As one-step adaptations spread more readily in the population, differences in the speed at which they fixate manifest early in the dynamics. Over time, cells gradually exhaust the possibilities for one-step adaptations and the difference between the adaptation rate in both processes becomes dominated by differences in the speed of acquiring two-step adaptations. In this latter phase we observe large differences between the two mutation processes in the expected time until any given proportion of the possible two-step adaptations have spread in the population. Analogous time differences for reaching fixed proportions of one-step adaptations in the initial phase of the dynamics are much smaller (Figure 5E).

In these analyses, parameters such as mutation rates, fitness effects, and population sizes were chosen so that simulations remained computationally manageable. Decreasing mutation rates toward more biologically realistic ranges decreases the likelihood that multiple lineages with adaptations emerge simultaneously and therefore weaken the effects of clonal interference. In contrast, larger population sizes and smaller fitness differences increase the likelihood of co-existence between competing lineages as each emerging lineage requires more time to fixate in the population, which strengthens the effects of clonal interference.

We also investigated how this pattern depends on the ruggedness of the fitness landscape (Figure 5—figure supplement 1A). To this end, we varied the relative proportions of loci with and without fitness valleys and quantified adaptability differences. We found that the duration of the initial phase in which the uniform process yields higher adaptability varies with the ruggedness of the fitness landscape: this phase is roughly three times as long when the proportion of loci with fitness valleys is 5% compared to when it is 95% (Figure 5—figure supplement 1B–D). However, over this range of proportions, we consistently observed that this initial phase only accounts for a small fraction of the time it takes to acquire all adaptations.

Discussion

Accumulating evidence suggests that mutagenesis in tumor cell populations proceeds in punctuated bursts rather than gradually; however, the effects of such punctuation on the ability of a tumor cell population to acquire and retain fitness-enhancing adaptations remain incompletely understood. Here, we set out to investigate the evolutionary dynamics of these processes using mathematical modeling and genomics data analysis of human tumors.

We found that when a population acquires multi-step adaptations via stochastic tunneling, the rate of adaptation is substantially enhanced by punctuation. Stochastic tunneling is the predominant mode of evolution in large populations that traverse fitness valleys. However, stochastic tunneling also occurs when the necessary mutation steps to reach a substantial fitness advantage are not individually maladaptive, that is if steps do not constitute a fitness valley but confer a slight fitness advantage. Punctuated evolution therefore facilitates the accumulation of sets of mutations that jointly and synergistically confer a fitness advantage, while the fitness effect of carrying only a subset of these mutations can range from being strongly maladaptive to being moderately adaptive. For the limit of low average mutation rates, we show analytically that the rate of stochastic tunneling under a temporally clustered process is kⁿ times larger than under the time-invariant process, where k is the fold-increase in the mutation rate relative to the average mutation rate during mutation bursts in the temporally clustered process, and n is the number of mutation steps that is required to exit the fitness valley. Moving toward higher average mutation rates, this fold-change in the stochastic tunneling rates decreases, but the absolute difference between the rates increases. Investigating these relative and absolute effects in a branching process as a function of the birth–death ratio elucidated a similar pattern: While the fold-change is highest if lineages in a fitness valley have a low birth–death ratio, we observed the largest difference between rates of stochastic tunneling at birth–death ratios close to one.

Applied to cancer evolution where cells’ birth–death ratios are influenced by a variety of internal and external factors, these findings highlight that for any given multi-step adaptation sequence, the magnitude of the facilitating effects described here likely varies depending on whether the tumor is in a phase of constant or growing population size.

While temporal clustering facilitates multi-step adaptations, it also impedes one-step adaptations due to clonal interference. We examined the interplay between these contrasting effects in a setting in which cells can acquire both types of adaptations. We showed that the relative importance of both effects varies over time. Since one-step adaptations tend to be acquired more readily, clonal interference matters most in the early phase of the process. As the share of possible two-step adaptations relative to one-step adaptations increases, differences in the tunneling rate become the more relevant determinant for the speed of adaptation. A uniform mutation rate, thus, makes a population more efficient at finding local fitness maxima. However, temporal clustering allows the population to more quickly move between local maxima, thus speeding up the search for a global maximum.

While our theoretical results are easily verified in simulation settings where we can choose a fitness landscape, applying these results to real data remains challenging because of the inherent complexities in determining fitness effects in biological systems. Estimating fitness effects of individual mutations requires either intricate experimental approaches (Salehi et al., 2021) or large patient cohorts to control for differences in the (epi-)genetic background against which these mutations emerge; such estimations become considerably harder when considering joint fitness effects of sets of multiple mutations.

To circumvent these challenges, we focused on TSGs as a representative set of gene pairs with synergistic effects and on one mutation process known to tend to fluctuate over time – APOBEC-driven mutagenesis. Moreover, we restricted attention to single-nucleotide substitutions as the principal mode of APOBEC-driven mutagenesis and the context in which mutation signatures are best characterized, and we did not consider other classes of mutations such as copy number changes or epigenomic modifications. Consequently, we observed only a subset of groups of mutations with synergistic fitness effects in the data and possibly only a small fraction of the variability in punctuation between different samples. Nevertheless, we found that these scores significantly correlate, which aligns with our model predictions. However, given this narrow scope, the generalization of our model toward a full quantitative estimate of the relevance of punctuation to adaptation dynamics in cancer requires additional investigations.

The observed empirical patterns may in part reflect the influence of unmeasured confounders. Rates of multi-step adaptation are determined not only by temporal clustering, but also by factors such as average mutation rates, fitness landscapes, and population-size trajectories. Structural differences in these factors between tumors with high or low APOBEC signature contributions therefore represent potential sources of confounding (for a review of known APOBEC-related biological mechanisms and clinical associations, see Butler and Banday, 2023). Additionally, we cannot measure rates of multi-step adaptation directly, and therefore need to consider a proxy; to this end, we chose a score that relates deactivating TSG mutations to the absolute mutation count. Structural differences affecting the relative rates of deactivating TSG mutations, such as APOBEC-related enrichment for nonsense mutations (Adler et al., 2023), compared to other mutations can therefore also influence our empirical findings.

Our results have implications for both mechanistic and statistical modeling of tumor evolution. A key quantity in mechanistic models is the rate at which different types of adaptations arise; our results show that this rate depends on the temporal dynamics of the mutation process. Accounting for this effect might be particularly important when modeling the emergence of treatment resistance in models used to optimize treatment schedules (McDonald et al., 2023). There is evidence linking APOBEC mutagenesis to resistance in breast cancer (Gupta et al., 2025). Our results suggest that APOBEC, and other punctuated mutation processes, are likely to be relevant more generally to settings in which developing resistance requires multi-step adaptation. Analogously, in statistical models, incorporating measures of punctuated mutagenesis may improve the ability to predict treatment outcomes by accounting for how these temporal dynamics shape adaptability.

Methods

Simulations of a Wright–Fisher and a branching process in an unbounded fitness landscape

In the simulations presented in Figure 2, Figure 2—figure supplement 1, mutations occur between selection steps and are independent of divisions. In the Wright–Fisher process simulations (Figure 2C) with a temporally clustered mutation rate, the out-of-burst mutation probability per cell and per update step (between two divisions) was set to 0.01 and was increased to 0.1 in bursts. For the corresponding uniform mutation rate trajectory, this probability was simply set equal to the total amount of mutations that occurred in the simulations under the temporally clustered mutation rate, divided by the total number of update steps.

Mutation rates in the branching model simulations (Figure 2—figure supplement 1) were chosen analogously. However, to achieve temporally equidistant mutation bursts in the branching process simulations, we scaled the burst duration and the time between bursts by the population size.

Simulations Wright–Fisher process in exploration/exploitation setting

To produce the results shown in Figure 3, we performed agent-based stochastic simulations of a Wright–Fisher process with 20 cells. Cells were characterized by their location on a discrete fitness landscape (a 30 by 30 two-dimensional lattice). The fitness values f associated with the positions on the lattice were drawn independently at random as f = 0.5 + x⁴ where x ∼ U[0,1] is distributed uniformly between zero and one. The lower bound of 0.5 ensures that cells at all positions have a non-negligible probability of dividing.

Raising x to the fourth power creates a landscape with few peak-positions surrounded by many valley-positions with little variation in fitness.

In each selection step, one cell was randomly chosen to divide with probability proportional to its fitness, and replaced a cell which was drawn uniformly at random. Mutations occurred between consecutive selection steps, and the mutating cell was chosen uniformly at random amongst all cells in the population (independently of the preceding selection step). If a cell mutated, it would move to a location in the fitness landscape chosen uniformly at random from the set of (at most 8) locations in the Moore neighborhood of its current location.

The rate at which mutation events occurred was governed by the parameters µ and k. Mutation bursts lasted 10² division events and started every 10³ division events. Fitness landscapes were re-drawn every 50,714 ( $\approx 50, 000 + \frac{10^{3}}{7}$ ) division events, to periodically vary the relative timing of mutation bursts and re-drawings of the fitness landscape. In this manner, we avoid artifacts in our results caused by phase-alignment of bursts and re-drawings.

Larger-scale simulations

We simulated tumor evolution as a branching process, starting from a single unmutated cell, up to a randomly drawn target cell number between 10⁶ and 10⁷ cells. We assume a constant death rate, so that over the course of a simulation the probability that the next event is a death event rather than a division event remains fixed (at a value of 0.35). In case of a death event, one cell picked uniformly at random is removed from the population. In case of a division event, one cell is chosen to divide with probabilities proportional to its fitness.

We assume that the number of mutations per division follows a binomial distribution Bin(10⁵, µ), and for each simulation we randomly draw a baseline mutation probability µ uniformly between 0 and 3 × 10⁻⁵. In mutation bursts, this probability gets multiplied by a factor of 20 or 50, again chosen randomly for each simulation. The population enters a burst phase with probability $\frac{n}{10}$ per division, where n is the current population size, and exits a burst phase with probability $\frac{n}{3}$ , so that in expectation bursts phases last 3 generations, and start every 10 generations. The first cell in a simulation is initialized to be in a burst with probability $\frac{3}{13}$ consistent with the expected time spent in bursts given those parameters.

The fitness effect of mutations is modeled as follows. The starting cell has a fitness of one. Each new mutation has an additive effect on the cell’s fitness. We assume that mutations occur uniformly at random across the genome and that multi-step adaptations can occur in a fraction θ = 0.01 of the genome, roughly aligning with the ratio of the number of TSGs vs the total number of genes in the human genome. We subdivide this part of the genome into 200 TSGs which are all hit by a mutation with equal probability. For each individual TSG, the first mutation reduces a cell’s fitness by 0.05. The second mutation increases fitness by 0.15, and all further mutations are fitness neutral. For the remaining (1 − θ) fraction of the genome, we assume that mutations are fitness-neutral with a probability of 0.9, and that fitness effects are otherwise drawn from a Gaussian with mean equal to −0.005 and standard deviation equal to 0.005.

To produce the results in Figure 4A, B, we keep track of all mutations that arise in a population, and remove all mutations with less than 1% variant allele frequency from the output, as those would be unlikely to be detected in WES. Moreover, we keep track of the fraction of mutations that a sample acquired during a mutation burst. For each simulation, we then count the number of TSGs for which there are at least two mutations found in the population and divide this by the total number of mutations in the sample.

Analysis of TCGA WES data

WES data was acquired from TCGA. We performed mutation signature decomposition using the cosmic fit() function in Python from the package SigProfilerAssignment (Díaz-Gay et al., 2023) version 0.1.8 with cosmic version 3.4. Moreover, we used the mutation-level signature probabilities generated via SigProfilerAssignment to compute mean probabilities of signatures SBS2 and SBS13 for individual SBSs (Figure 4D, E).

To compute our TSG deactivation score, we filtered the SNV data for missense and nonsense mutations in genes belonging to the TSGene 2.0 database (Zhao et al., 2016). For each sample we calculated the number of TSGs with at least two mutations and divided by the total number of SNVs.

Analogously, to construct our synergistic mutations score, we filtered for missense and nonsense mutations in frequently co-mutated gene pairs identified by Gu et al., 2010, and for each sample divided the number of pairs in this list with mutations in both genes by the number of SNVs (Park et al., 2021). This enrichment for co-occurrences of mutations in both genes in a pair suggests that mutations in the two genes tend to have a synergistic fitness effect. For each sample in the TCGA WES data we thus computed the number of gene pairs from this list for which there is at least one non-synonymous single-nucleotide substitution mutation in each of the two genes, and divided this number by the total number of single-nucleotide substitutions in the sample.

As a robustness check, we investigated whether the results change if we require a minimum distance between the genomic locations of the mutations in TSGs that we count to our TSG deactivation score. APOBEC has been linked to clustered mutagenesis through mutation processes of kataegis and omikli (Petljak et al., 2022; Nik-Zainal et al., 2012; Mas-Ponte and Supek, 2020). If samples with higher APOBEC activity have higher rates of having multiple close-by mutations on the same allele, and if such a sample has multiple deactivating mutations in a TSG, these inactivating mutations might be more likely to be localized on the same allele compared to samples with lower APOBEC activity. Since we use appearances of more than one inactivating mutation as a proxy for bi-allelic deactivation of TSGs, such a pattern would bias our results. To account for that, we explored filtering TSG mutations in the data based on different minimum-distance-thresholds ranging from 1 to 100 bp and found no impact on our results. We did not perform any analyses with allele-specific mutation information, as determining allele-specific TSG inactivation relies on informative co-occurring mutations. For any given mutation, tumors with lower mutation burden are less likely to harbor such co-occurring mutations, and preferentially filtering out observations in tumors with lower mutation burden would have biased our findings.

Appendix 1

Derivations of results on multi-step adaptation rates in the limit of infrequent mutations

In the main text, we discussed differences in the rate of multi-step adaptations f_k as a function of a temporal clustering parameter k. Here, we formalize this discussion and derive the result presented in the main text. We consider processes in which multi-step adaptations occur via stochastic tunneling rather than sequential fixation, that is processes in which mutants with a selective disadvantage are negligibly unlikely to reach fixation. We assume that mutations happen sufficiently infrequently, so that we can neglect valley crossing scenarios in which the same mutation in a multi-step adaptation sequence occurs multiple times. Additionally, we require that valley crossing happens solely via stochastic tunneling, that is that the probability that a maladaptive mutant fixates is vanishingly low.

We denote the average rate at which cells acquire mutations by µ, and compare dynamics under different mutation processes which we index by a temporal clustering parameter k ≥ 1. Starting every d units of time, a process with parameter k undergoes a mutation burst lasting $\frac{d}{k}$ units of time, during which mutations occur at a rate kµ. Outside of bursts, the mutation rate is 0. We assume that the time $\frac{(k - 1) d}{k}$ between subsequent burst phases is sufficiently long such that we can restrict attention to multi-step adaptations that occur during a single burst, and that the duration of each burst $\frac{d}{k}$ is sufficiently long relative to the time that it takes to cross a fitness valley, so that valley crossings, which fail because the burst phase ended but which would have been successful, have a negligible effect on f_k.

As stated in the main text, with this setup we can show that for any population dynamics process for which the above limits can be motivated, the fold-increase in the rate f_k at which (n + 1)-step adaptations occur in a process with clustering parameter k relative to the uniform process approaches kⁿ.

lim_{d \to \infty} lim_{μ \to \infty} \frac{f_{k}}{f_{1}} = lim_{μ \to \infty} lim_{d \to \infty} \frac{f_{k}}{f_{1}} = k^{n}

Deriving this result for any selection process indeed becomes straightforward, if it can be motivated that the mutation rate does not affect the fate of any mutant lineage once the first mutant in this lineage has emerged.

Somewhat more formally, we index the steps in an n-step adaptation sequence by i ∈ {1, …, n}, and denote the expected size of the lineage descending from mutant i by Y_i(t). This quantity Y_i(t) of course depends on the specificities of the selection process, but for our derivations it suffices to only require that Y_i(t) does not depend on the mutation rate.

Since all intermediate mutants have a selective disadvantage, and since we assume that we are in a parameter regime in which the likelihood that a disadvantageous mutant fixates is negligibly low, it follows that the integral $\int_{0}^{\infty} Y_{i} (t) d t$ converges.

\int_{0}^{\infty} Y_{i} (t) d t < \infty

In the limit of a low (and for now time-invariant) mutation rate µ, the probability that the i’th mutant produces a further mutant P(i → i + 1) is simply the product of the mutation rate and this integral.

P (i \to i + 1) = μ \int_{0}^{\infty} Y_{i} (t) d t

The probability that the first mutant spawns a sequence of n + 1 mutants can, thus, be written as follows:

\begin{array}{ll} P (1 \to n + 1) & = \prod_{i = 1}^{n} (μ \int_{0}^{\infty} Y_{i} (t) d t) \\ = μ^{n} \prod_{i = 1}^{n} (\int_{0}^{\infty} Y_{i} (t) d t) \end{array}

Lastly, we use $\tilde{μ}$ to denote the rate at which new mutants with only one mutation emerge. For a given selection process, this rate may vary with the population size. For the Wright–Fisher models of constant population size $N,$ this expression simplifies to $\tilde{μ} = N μ$ . We arrive at the following formulation for the rate at which new advantageous mutants emerge f₁.

\begin{array}{ll} f_{1} & = \tilde{μ} \cdot P (1 \to n + 1) \\ = \tilde{μ} \cdot μ^{n} \cdot \prod_{i = 1}^{n} (\int_{0}^{\infty} Y_{i} (t) d t) \end{array}

Finally, we can introduce our clustering parameter k to scale both µ and $\tilde{μ}$ , and multiply by a factor of $\frac{1}{k}$ to arrive at the average rate of valley crossing in the temporally clustered process with parameter k.

\begin{array}{ll} f_{k} & = \frac{1}{k} \cdot (k \cdot \tilde{μ}) \cdot {(k \cdot μ)}^{n} \cdot \prod_{i = 1}^{n} (\int_{0}^{\infty} Y_{i} (t) d t) \\ = k^{n} \cdot \tilde{μ} \cdot μ^{n} \cdot \prod_{i = 1}^{n} (\int_{0}^{\infty} Y_{i} (t) d t) \\ = k^{n} \cdot f_{1} \end{array}

One consequence of our assumption that none of the maladaptive intermediate mutants reaches fixation is that the composition of the population once the final mutant emerges in the limit of µ → 0 does not depend on the mutation rate. For population dynamics models with a constant population size in which the probability that an emerging mutant (absent further mutation events) reaches fixation only depends on this composition, such as the Wright–Fisher process, we can interpret $\frac{f_{k}}{f_{1}}$ therefore also as the ratio of the rates at which mutants with n + 1 mutations fixate.

Analogously, in certain models of branching evolution in which the prospects of one branch do not depend on other co-evolving branches, such as the Galton–Watson process (Watson and Galton, 1875), whether the lineage of an emerging mutant with (n + 1) mutations survives is independent of the mutation rate in the limit of µ → 0. The ratio $\frac{f_{k}}{f_{1}}$ therefore also reflects the ratio of the rates at which surviving lineages with (n + 1) adaptations arise in such models. We provide thorough treatment of dynamics in such processes below (Appendix 2).

The above result suggests that the effect of temporal clustering on relative rates of valley crossing is substantial, and gets exponentially stronger the wider the fitness valley is (n). We validated this result with simulations of a Wright–Fisher process for a range of small values for µ (Figure 2—figure supplement 2B, E). As we move away from the limit of rare mutations, the relative effect of temporal clustering $\frac{f_{k}}{f_{1}}$ gets smaller, as will be discussed below. However, as the absolute rates f_k and f₁ are proportional to µⁿ⁺¹, the absolute effect of temporal clustering on the rate of fixations of mutants with multi-step adaptations initially increases in µ (Figure 2—figure supplement 2D).

Appendix 2

Dynamics in a branching process

In the Galton–Watson process, a branching process with constant birth and death rates, it is possible to derive explicit formulas for the probability that a lineage will acquire a mutation before going extinct (often termed ‘evolutionary rescue’ Azevedo and Olofsson, 2021). Formulating this probability for the lineage of the first mutant in a two-step adaptation sequence thus allows us to explore numerically how the magnitude of the effects that we describe behaves across a wide range of parameter combinations, which is especially attractive for parameter regions for which simulations become prohibitively expensive.

We thus investigated this probability of evolutionary rescue, which translates to $f_{k} / μ$ in our model, as a function of $k, μ$ and the birth/death ratio (Figure 2—figure supplement 3A–C). Note that the birth/death ratio that we considered here is that of the first mutant, not that of the wild-type cell. Unsurprisingly, we found that as µ decreases, rescue probabilities for any k converge to zero for subcritical processes (birth/death ratios < 1). For supercritical processes (birth/death ratios > 1), rescue probabilities converge to one minus the extinction probability of the lineage – rescue only happens if the process expands indefinitely. Throughout, rescue probabilities of course increase in the expected time that the lineage persists, that is they increase with the birth/death ratio. However, this increase is most steep for large µ and flattens for the subcritical range as µ decreases.

For the ratio of these escape probabilities for different k relative to that of the uniform process, that is for $\frac{\frac{f_{k}}{μ}}{\frac{f_{1}}{μ}} = \frac{f_{k}}{f_{1}}$ , we found that in the subcritical regime, as µ decreases, $\frac{f_{k}}{f_{1}}$ approaches k, consistent with our derivations in the previous section (Appendix 1). In the supercritical regime, the condition $\int_{0}^{\infty} Y (t) d t < \infty$ used in our derivations above is violated – lineages have a non-zero chance of branching out indefinitely. Here we saw that as µ decreases toward zero, $\frac{f_{k}}{f_{1}}$ approaches 1 (Figure 2—figure supplement 3D–F).

Moreover, we found that $\frac{f_{k}}{f_{1}}$ consistently decreases with increasing birth–death ratio. As a result of the interplay between this decrease in $\frac{f_{k}}{f_{1}}$ and increase in the absolute rescue probability, we found that the absolute effect, f_k−f₁, is highest for birth/death ratios close to one (Figure 2—figure supplement 3G). These results remain as µ decreases (Figure 2—figure supplement 3H, I) but the magnitude of the absolute effect of course shrinks, as rescue of lineages that would otherwise go extinct (all subcritical lineages, and those lineages in the supercritical regime that go extinct by chance) generally becomes unlikely at small values of µ.

Appendix 3

Dynamics as a function of time between bursts

In our derivations, we assumed that bursts are sufficiently spaced out in time, and last sufficiently long so that we can restrict attention to mutant lineages that emerge and go extinct within the duration of a single burst ( $d \to \infty$ ). To investigate how effect sizes behave away from this limit, as a function of the time between bursts, we numerically approximated rescue probabilities for finite $d$ .

Specifically, we approximated the Galton–Watson process by a Markov process with finitely many states $s \in {0, 1, \dots, N, r e s c u e}$ , starting with $s (0) = 1$ , where $s \in {0, 1, \dots, N}$ corresponds to the number of viable cells. We considered 0, $N$ and $r e s c u e$ as absorbing states, and otherwise modeled transitions to states $s - 1$ or $s + 1$ as occurring at rates ${s \cdot r}_{d}$ and ${s \cdot r}_{b}$ , respectively, with constant birth and death rates $r_{d}$ and $r_{b}$ , and transitions to the $r e s c u e$ state as occurring at rate $s \cdot k \cdot μ$ within bursts and at rate zero otherwise. We iteratively solved the corresponding forward equations of this continuous time Markov process numerically for subsequent cycles of in-burst and out-of-burst periods, until the summed mass on states ${1, \dots, N - 1}$ fell below a tolerance threshold, and we considered the joint mass on states $N$ and $r e s c u e$ as approximation for the rescue probability. For our results in Figure 2—figure supplement 4, we used N = 300 and a tolerance threshold of 10⁻¹⁰.

Generally, we showed that differences in multi-step adaptation rates emerge because lineages with a first mutation experience higher mutation rates under a temporally clustered mutation process compared to a uniform mutation process (Appendix 1). As the burst duration and the time-interval between subsequent bursts decrease, this expected difference in mutation rates for mutant lineages diminishes because of the increased likelihood that a mutant lineage that emerged in a burst also experiences a non-burst phase. Consequently, we found that lowering $d$ diminishes the effect of temporal clustering on the probability of producing a second mutant (Figure 2—figure supplement 4A–C). This pattern manifests both in the relative effects, $f_{k} / f_{1}$ (Figure 2—figure supplement 4D–F) and in the absolute effects $f_{k} - f_{1}$ (Figure 2—figure supplement 4H–J). However, for the absolute differences $f_{k} - f_{1}$ , we found that the strength of this dependency varies non-monotonously with the birth–death ratio. While absolute differences under large $d$ peak close to a birth–death ratio of one (Figure 2—figure supplements 3–4J), we also observed the strongest decrease in absolute differences at these values for the birth–death ratio, so that under low $d$ this peak turns into a local minimum of the absolute-difference curve (Figure 2—figure supplement 4H, I). This effect arises because the expected lifetime of a mutant lineage before going extinct increases steeply when moving from a subcritical process (birth–death ratios <1) to a critical process (birth–death ratio = 1), making such lineages with higher birth–death ratios more sensitive to changes in mutation rates that occur long after the lineage first emerged.

Appendix 4

Dynamics in exploration–exploitation setting

Our simulations show that for a given µ the effect of increasing the clustering parameter k on the average fitness becomes negative at some k (Figure 3C). In the simulations in this section, this effect is driven by two main factors. First, after a re-drawing, the population often is not in a local maximum of the fitness landscape, and hence may be able to increase its fitness by single mutations (one-step adaptations), without tunneling. Waiting for a burst to occur and thereby delaying such local explorations decreases average fitness. Second, once the population has found a peak in the fitness landscape, having very high mutation rates in bursts temporarily scatters cells to lower points in the landscape, and may even cause drift to points of lower fitness. Both of these effects become apparent when considering a representative snapshot of the simulations (Figure 3D). The misalignment of re-drawing and burst causes the population to remain in a valley for several divisions until the first burst occurs. Moreover, relative to similar snapshots of simulations with lower k (Figure 3E, F), once the population has reached a peak, the scattering during bursts at higher k leads to sharper decreases in fitness, and the population may even remain in a lower fitness point after a burst (Figure 3D).

This increased scattering can also be seen when going from k = 1 to k = 5 (Figure 3E, F). However, at k = 1, the population spends much time in local maxima, as large jumps in fitness only occur right after re-drawings, whereas at k = 5, jumps occur also long after the landscape was re-drawn and the population found a local maximum. These later jumps to higher points in the landscape correspond to tunneling events.

Appendix 5

TSG deactivation and ROS-associated mutagenesis

Following a reviewer suggestion, we also explored the contribution of mutation signatures attributed to reactive oxygen species (ROS; Figure 4—figure supplement 4). However, we found that these signatures account for only small fractions of the total observed mutation signatures, suggesting that these mutagenic processes would not account for much of the fluctuation behavior of overall mutation rates, and we saw no consistent correlations with TSG inactivation.

Data availability

The results shown here are in part based upon data generated by the TCGA Research Network: https://www.cancer.gov/tcga. Simulation code is publicly available at https://github.com/Michorlab/Punctuated-mutagenesis (copy archived at Graser, 2026).

References

(2023) Mutational processes of tobacco smoking and APOBEC activity generate protein-truncating mutations in cancer genomes
Science Advances 9:eadh3083.

https://doi.org/10.1126/sciadv.adh3083
- PubMed
- Google Scholar
1. Alexandrov LB
2. Nik-Zainal S
3. Wedge DC
4. Aparicio SAJR
5. Behjati S
6. Biankin AV
7. Bignell GR
8. Bolli N
9. Borg A
10. Børresen-Dale A-L
11. Boyault S
12. Burkhardt B
13. Butler AP
14. Caldas C
15. Davies HR
16. Desmedt C
17. Eils R
18. Eyfjörd JE
19. Foekens JA
20. Greaves M
21. Hosoda F
22. Hutter B
23. Ilicic T
24. Imbeaud S
25. Imielinski M
26. Jäger N
27. Jones DTW
28. Jones D
29. Knappskog S
30. Kool M
31. Lakhani SR
32. López-Otín C
33. Martin S
34. Munshi NC
35. Nakamura H
36. Northcott PA
37. Pajic M
38. Papaemmanuil E
39. Paradiso A
40. Pearson JV
41. Puente XS
42. Raine K
43. Ramakrishna M
44. Richardson AL
45. Richter J
46. Rosenstiel P
47. Schlesner M
48. Schumacher TN
49. Span PN
50. Teague JW
51. Totoki Y
52. Tutt ANJ
53. Valdés-Mas R
54. van Buuren MM
55. van ’t Veer L
56. Vincent-Salomon A
57. Waddell N
58. Yates LR
59. Australian Pancreatic Cancer Genome Initiative
60. ICGC Breast Cancer Consortium
61. ICGC MMML-Seq Consortium
62. ICGC PedBrain
63. Zucman-Rossi J
64. Futreal PA
65. McDermott U
66. Lichter P
67. Meyerson M
68. Grimmond SM
69. Siebert R
70. Campo E
71. Shibata T
72. Pfister SM
73. Campbell PJ
74. Stratton MR
(2013) Signatures of mutational processes in human cancer
Nature 500:415–421.

https://doi.org/10.1038/nature12477
- PubMed
- Google Scholar
(2015) The mathematics of cancer: integrating quantitative models
Nature Reviews. Cancer 15:730–745.

https://doi.org/10.1038/nrc4029
- PubMed
- Google Scholar
1. Araten DJ
2. Golde DW
3. Zhang RH
4. Thaler HT
5. Gargiulo L
6. Notaro R
7. Luzzatto L
(2005) A quantitative measurement of the human somatic mutation rate
Cancer Research 65:8111–8117.

https://doi.org/10.1158/0008-5472.CAN-04-1198
- PubMed
- Google Scholar
(2015) Stochastic tunneling and metastable states during the somatic evolution of cancer
Genetics 199:1213–1228.

https://doi.org/10.1534/genetics.114.171553
- PubMed
- Google Scholar
1. Azevedo RBR
2. Olofsson P
(2021) A branching process model of evolutionary rescue
Mathematical Biosciences 341:108708.

https://doi.org/10.1016/j.mbs.2021.108708
- PubMed
- Google Scholar
1. Baca SC
2. Prandi D
3. Lawrence MS
4. Mosquera JM
5. Romanel A
6. Drier Y
7. Park K
8. Kitabayashi N
9. MacDonald TY
10. Ghandi M
11. Van Allen E
12. Kryukov GV
13. Sboner A
14. Theurillat J-P
15. Soong TD
16. Nickerson E
17. Auclair D
18. Tewari A
19. Beltran H
20. Onofrio RC
21. Boysen G
22. Guiducci C
23. Barbieri CE
24. Cibulskis K
25. Sivachenko A
26. Carter SL
27. Saksena G
28. Voet D
29. Ramos AH
30. Winckler W
31. Cipicchio M
32. Ardlie K
33. Kantoff PW
34. Berger MF
35. Gabriel SB
36. Golub TR
37. Meyerson M
38. Lander ES
39. Elemento O
40. Getz G
41. Demichelis F
42. Rubin MA
43. Garraway LA
(2013) Punctuated evolution of prostate cancer genomes
Cell 153:666–677.

https://doi.org/10.1016/j.cell.2013.03.021
- PubMed
- Google Scholar
1. Butler K
2. Banday AR
(2023) APOBEC3-mediated mutagenesis in cancer: causes, clinical significance and therapeutic potential
Journal of Hematology & Oncology 16:31.

https://doi.org/10.1186/s13045-023-01425-5
- PubMed
- Google Scholar
Book
(2004) Modelling stochastic clonal interference
In: Ciobanu G, Rozenberg G, editors. Modelling in Molecular Biology. Springer. pp. 21–38.

https://doi.org/10.1007/978-3-642-18734-6_2
- Google Scholar
1. Campos PRA
2. de Oliveira VM
(2004) Mutational effects on the clonal interference phenomenon
Evolution; International Journal of Organic Evolution 58:932–937.

https://doi.org/10.1111/j.0014-3820.2004.tb00427.x
- PubMed
- Google Scholar
1. Chakravarty D
2. Gao J
3. Phillips SM
4. Kundra R
5. Zhang H
6. Wang J
7. Rudolph JE
8. Yaeger R
9. Soumerai T
10. Nissan MH
11. Chang MT
12. Chandarlapaty S
13. Traina TA
14. Paik PK
15. Ho AL
16. Hantash FM
17. Grupe A
18. Baxi SS
19. Callahan MK
20. Snyder A
21. Chi P
22. Danila D
23. Gounder M
24. Harding JJ
25. Hellmann MD
26. Iyer G
27. Janjigian Y
28. Kaley T
29. Levine DA
30. Lowery M
31. Omuro A
32. Postow MA
33. Rathkopf D
34. Shoushtari AN
35. Shukla N
36. Voss M
37. Paraiso E
38. Zehir A
39. Berger MF
40. Taylor BS
41. Saltz LB
42. Riely GJ
43. Ladanyi M
44. Hyman DM
45. Baselga J
46. Sabbatini P
47. Solit DB
48. Schultz N
(2017) OncoKB: a precision oncology knowledge base
JCO Precision Oncology 2017:1–16.

https://doi.org/10.1200/PO.17.00011
- PubMed
- Google Scholar
1. Davis A
2. Gao R
3. Navin N
(2017) Tumor evolution: Linear, branching, neutral or punctuated? Biochimica et Biophysica Acta (BBA)
Reviews on Cancer 1867:151–161.

https://doi.org/10.1016/j.bbcan.2017.01.003
- PubMed
- Google Scholar
1. Deming DA
2. Leystra AA
3. Nettekoven L
4. Sievers C
5. Miller D
6. Middlebrooks M
7. Clipson L
8. Albrecht D
9. Bacher J
10. Washington MK
11. Weichert J
12. Halberg RB
(2014) PIK3CA and APC mutations are synergistic in the development of intestinal cancers
Oncogene 33:2245–2254.

https://doi.org/10.1038/onc.2013.167
- PubMed
- Google Scholar
1. Díaz-Gay M
2. Vangara R
3. Barnes M
4. Wang X
5. Islam SMA
6. Vermes I
7. Duke S
8. Narasimman NB
9. Yang T
10. Jiang Z
11. Moody S
12. Senkin S
13. Brennan P
14. Stratton MR
15. Alexandrov LB
(2023) Assigning mutational signatures to individual samples and individual somatic mutations with SigProfilerAssignment
Bioinformatics 39:btad756.

https://doi.org/10.1093/bioinformatics/btad756
- PubMed
- Google Scholar
1. Drews RM
2. Hernando B
3. Tarabichi M
4. Haase K
5. Lesluyes T
6. Smith PS
7. Morrill Gavarró L
8. Couturier D-L
9. Liu L
10. Schneider M
11. Brenton JD
12. Van Loo P
13. Macintyre G
14. Markowetz F
(2022) A pan-cancer compendium of chromosomal instability
Nature 606:976–983.

https://doi.org/10.1038/s41586-022-04789-9
- PubMed
- Google Scholar
(2006) How few cancer cells can be detected by positron emission tomography? A frequent question addressed by an in vitro study
European Journal of Nuclear Medicine and Molecular Imaging 33:697–702.

https://doi.org/10.1007/S00259-005-0038-6/METRICS
- PubMed
- Google Scholar
1. Fisher RA
(1923) XXI.—On the Dominance Ratio
Proceedings of the Royal Society of Edinburgh 42:321–341.

https://doi.org/10.1017/S0370164600023993
- Google Scholar
(2008) Clonal interference, multiple mutations and adaptation in large asexual populations
Genetics 180:2163–2173.

https://doi.org/10.1534/genetics.108.090019
- PubMed
- Google Scholar
(2019) Tumor mutation burden: from comprehensive mutational screening to the clinic
Cancer Cell International 19:209.

https://doi.org/10.1186/s12935-019-0929-4
- PubMed
- Google Scholar
1. Gao R
2. Davis A
3. McDonald TO
4. Sei E
5. Shi X
6. Wang Y
7. Tsai P-C
8. Casasent A
9. Waters J
10. Zhang H
11. Meric-Bernstam F
12. Michor F
13. Navin NE
(2016) Punctuated copy number evolution and clonal stasis in triple-negative breast cancer
Nature Genetics 48:1119–1130.

https://doi.org/10.1038/ng.3641
- PubMed
- Google Scholar
1. Gerrish PJ
2. Lenski RE
(1998) The fate of competing beneficial mutations in an asexual population
Genetica 102–103:127–144.

https://doi.org/10.1023/A:1017067816551
- PubMed
- Google Scholar
Software
1. Graser CJ
(2026) Punctuated-mutagenesis, version swh:1:rev:610ba8e16897dbc5216e0718b368da666d09fa23
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:89771b9bb260d6770ee2ce89430a0a96390e9541;origin=https://github.com/Michorlab/Punctuated-mutagenesis;visit=swh:1:snp:8e8f4b4d78e84c70b3362de2fa58bd603981009e;anchor=swh:1:rev:610ba8e16897dbc5216e0718b368da666d09fa23
1. Gu Y
2. Yang D
3. Zou J
4. Ma W
5. Wu R
6. Zhao W
7. Zhang Y
8. Xiao H
9. Gong X
10. Zhang M
11. Zhu J
12. Guo Z
(2010) Systematic interpretation of comutated genes in large-scale cancer mutation profiles
Molecular Cancer Therapeutics 9:2186–2195.

https://doi.org/10.1158/1535-7163.MCT-10-0022
- PubMed
- Google Scholar
1. Gupta A
2. Gazzo A
3. Selenica P
4. Safonov A
5. Pareja F
6. da Silva EM
7. Brown DN
8. Shao H
9. Zhu Y
10. Patel J
11. Blanco-Heredia J
12. Stefanovska B
13. Carpenter MA
14. Chen Y
15. Vegas I
16. Pei X
17. Frosina D
18. Jungbluth AA
19. Ladanyi M
20. Curigliano G
21. Weigelt B
22. Riaz N
23. Powell SN
24. Razavi P
25. Harris RS
26. Reis-Filho JS
27. Marra A
28. Chandarlapaty S
(2025) APOBEC3 mutagenesis drives therapy resistance in breast cancer
Nature Genetics 57:1452–1462.

https://doi.org/10.1038/s41588-025-02187-1
- PubMed
- Google Scholar
1. Hsu MS
2. Yu JC
3. Wang HW
4. Chen ST
5. Hsiung CN
6. Ding S
7. Wu PE
8. Shen CY
9. Cheng CW
(2010) Synergistic effects of polymorphisms in DNA repair genes and endogenous estrogen exposure on female breast cancer risk
Annals of Surgical Oncology 17:760–771.

https://doi.org/10.1245/s10434-009-0802-0
- PubMed
- Google Scholar
1. Issa JP
(2022) The two-hit hypothesis meets epigenetics
Cancer Research 82:1167–1169.

https://doi.org/10.1158/0008-5472.CAN-22-0405
- PubMed
- Google Scholar
(2004) Stochastic tunnels in evolutionary dynamics
Genetics 166:1571–1579.

https://doi.org/10.1534/genetics.166.3.1571
- PubMed
- Google Scholar
1. Jamal-Hanjani M
2. Wilson GA
3. McGranahan N
4. Birkbak NJ
5. Watkins TBK
6. Veeriah S
7. Shafi S
8. Johnson DH
9. Mitter R
10. Rosenthal R
11. Salm M
12. Horswell S
13. Escudero M
14. Matthews N
15. Rowan A
16. Chambers T
17. Moore DA
18. Turajlic S
19. Xu H
20. Lee S-M
21. Forster MD
22. Ahmad T
23. Hiley CT
24. Abbosh C
25. Falzon M
26. Borg E
27. Marafioti T
28. Lawrence D
29. Hayward M
30. Kolvekar S
31. Panagiotopoulos N
32. Janes SM
33. Thakrar R
34. Ahmed A
35. Blackhall F
36. Summers Y
37. Shah R
38. Joseph L
39. Quinn AM
40. Crosbie PA
41. Naidu B
42. Middleton G
43. Langman G
44. Trotter S
45. Nicolson M
46. Remmen H
47. Kerr K
48. Chetty M
49. Gomersall L
50. Fennell DA
51. Nakas A
52. Rathinam S
53. Anand G
54. Khan S
55. Russell P
56. Ezhil V
57. Ismail B
58. Irvin-Sellers M
59. Prakash V
60. Lester JF
61. Kornaszewska M
62. Attanoos R
63. Adams H
64. Davies H
65. Dentro S
66. Taniere P
67. O’Sullivan B
68. Lowe HL
69. Hartley JA
70. Iles N
71. Bell H
72. Ngai Y
73. Shaw JA
74. Herrero J
75. Szallasi Z
76. Schwarz RF
77. Stewart A
78. Quezada SA
79. Le Quesne J
80. Van Loo P
81. Dive C
82. Hackshaw A
83. Swanton C
84. TRACERx Consortium
(2017) Tracking the evolution of non-small-cell lung cancer
The New England Journal of Medicine 376:2109–2121.

https://doi.org/10.1056/NEJMoa1616288
- PubMed
- Google Scholar
1. Knudson AG
(1971) Mutation and cancer: statistical study of retinoblastoma
PNAS 68:820–823.

https://doi.org/10.1073/pnas.68.4.820
- PubMed
- Google Scholar
(2003) Mutation-selection networks of cancer initiation: tumor suppressor genes and chromosomal instability
Journal of Theoretical Biology 223:433–450.

https://doi.org/10.1016/s0022-5193(03)00120-6
- PubMed
- Google Scholar
1. Lee J-K
2. Lee J
3. Kim S
4. Kim S
5. Youk J
6. Park S
7. An Y
8. Keam B
9. Kim D-W
10. Heo DS
11. Kim YT
12. Kim J-S
13. Kim SH
14. Lee JS
15. Lee S-H
16. Park K
17. Ku J-L
18. Jeon YK
19. Chung DH
20. Park PJ
21. Kim J
22. Kim TM
23. Ju YS
(2017) Clonal history and genetic predictors of transformation into small-cell carcinomas from lung adenocarcinomas
Journal of Clinical Oncology 35:3065–3074.

https://doi.org/10.1200/JCO.2016.71.9096
- PubMed
- Google Scholar
1. Mas-Ponte D
2. Supek F
(2020) DNA mismatch repair promotes APOBEC3-mediated diffuse hypermutation in human cancers
Nature Genetics 52:958–968.

https://doi.org/10.1038/s41588-020-0674-6
- PubMed
- Google Scholar
1. McDonald TO
2. Cheng YC
3. Graser C
4. Nicol PB
5. Temko D
6. Michor F
(2023) Computational approaches to modelling and optimizing cancer treatment
Nature Reviews Bioengineering 1:695–711.

https://doi.org/10.1038/s44222-023-00089-7
- Google Scholar
1. Merlo LMF
2. Pepper JW
3. Reid BJ
4. Maley CC
(2006) Cancer as an evolutionary and ecological process
Nature Reviews. Cancer 6:924–935.

https://doi.org/10.1038/nrc2013
- PubMed
- Google Scholar
1. Minussi DC
2. Nicholson MD
3. Ye H
4. Davis A
5. Wang K
6. Baker T
7. Tarabichi M
8. Sei E
9. Du H
10. Rabbani M
11. Peng C
12. Hu M
13. Bai S
14. Lin Y-W
15. Schalck A
16. Multani A
17. Ma J
18. McDonald TO
19. Casasent A
20. Barrera A
21. Chen H
22. Lim B
23. Arun B
24. Meric-Bernstam F
25. Van Loo P
26. Michor F
27. Navin NE
(2021) Breast tumours maintain a reservoir of subclonal diversity during expansion
Nature 592:302–308.

https://doi.org/10.1038/s41586-021-03357-x
- PubMed
- Google Scholar
1. Navin N
2. Kendall J
3. Troge J
4. Andrews P
5. Rodgers L
6. McIndoo J
7. Cook K
8. Stepansky A
9. Levy D
10. Esposito D
11. Muthuswamy L
12. Krasnitz A
13. McCombie WR
14. Hicks J
15. Wigler M
(2011) Tumour evolution inferred by single-cell sequencing
Nature 472:90–94.

https://doi.org/10.1038/nature09807
- PubMed
- Google Scholar
(2021) The genomic processes of biological invasions: from invasive species to cancer metastases and back again
Frontiers in Ecology and Evolution 9:681100.

https://doi.org/10.3389/fevo.2021.681100
- Google Scholar
1. Nik-Zainal S
2. Alexandrov LB
3. Wedge DC
4. Van Loo P
5. Greenman CD
6. Raine K
7. Jones D
8. Hinton J
9. Marshall J
10. Stebbings LA
11. Menzies A
12. Martin S
13. Leung K
14. Chen L
15. Leroy C
16. Ramakrishna M
17. Rance R
18. Lau KW
19. Mudie LJ
20. Varela I
21. McBride DJ
22. Bignell GR
23. Cooke SL
24. Shlien A
25. Gamble J
26. Whitmore I
27. Maddison M
28. Tarpey PS
29. Davies HR
30. Papaemmanuil E
31. Stephens PJ
32. McLaren S
33. Butler AP
34. Teague JW
35. Jönsson G
36. Garber JE
37. Silver D
38. Miron P
39. Fatima A
40. Boyault S
41. Langerød A
42. Tutt A
43. Martens JWM
44. Aparicio SAJR
45. Borg Å
46. Salomon AV
47. Thomas G
48. Børresen-Dale A-L
49. Richardson AL
50. Neuberger MS
51. Futreal PA
52. Campbell PJ
53. Stratton MR
54. Breast Cancer Working Group of the International Cancer Genome Consortium
(2012) Mutational processes molding the genomes of 21 breast cancers
Cell 149:979–993.

https://doi.org/10.1016/j.cell.2012.04.024
- PubMed
- Google Scholar
1. Nowak MA
2. Michor F
3. Komarova NL
4. Iwasa Y
(2004) Evolutionary dynamics of tumor suppressor gene inactivation
PNAS 101:10635–10638.

https://doi.org/10.1073/pnas.0400747101
- PubMed
- Google Scholar
1. Paige AJW
(2003) Redefining tumour suppressor genes: exceptions to the two-hit hypothesis
Cellular and Molecular Life Sciences 60:2147–2163.

https://doi.org/10.1007/s00018-003-3027-6
- PubMed
- Google Scholar
1. Park SC
2. Krug J
(2007) Clonal interference in large populations
PNAS 104:18135–18140.

https://doi.org/10.1073/pnas.0705778104
- PubMed
- Google Scholar
1. Park S
2. Supek F
3. Lehner B
(2021) Higher order genetic interactions switch cancer genes from two-hit to one-hit drivers
Nature Communications 12:1–10.

https://doi.org/10.1038/s41467-021-27242-3
- PubMed
- Google Scholar
1. Petljak M
2. Alexandrov LB
3. Brammeld JS
4. Price S
5. Wedge DC
6. Grossmann S
7. Dawson KJ
8. Ju YS
9. Iorio F
10. Tubio JMC
11. Koh CC
12. Georgakopoulos-Soares I
13. Rodríguez-Martín B
14. Otlu B
15. O’Meara S
16. Butler AP
17. Menzies A
18. Bhosle SG
19. Raine K
20. Jones DR
21. Teague JW
22. Beal K
23. Latimer C
24. O’Neill L
25. Zamora J
26. Anderson E
27. Patel N
28. Maddison M
29. Ng BL
30. Graham J
31. Garnett MJ
32. McDermott U
33. Nik-Zainal S
34. Campbell PJ
35. Stratton MR
(2019) Characterizing mutational signatures in human cancer cell lines reveals episodic APOBEC mutagenesis
Cell 176:1282–1294.

https://doi.org/10.1016/j.cell.2019.02.012
- PubMed
- Google Scholar
1. Petljak M
2. Dananberg A
3. Chu K
4. Bergstrom EN
5. Striepen J
6. von Morgen P
7. Chen Y
8. Shah H
9. Sale JE
10. Alexandrov LB
11. Stratton MR
12. Maciejowski J
(2022) Mechanisms of APOBEC3 mutagenesis in human cancer cells
Nature 607:799–807.

https://doi.org/10.1038/s41586-022-04972-y
- PubMed
- Google Scholar
1. Proulx SR
(2011) The rate of multi-step evolution in Moran and Wright-Fisher populations
Theoretical Population Biology 80:197–207.

https://doi.org/10.1016/j.tpb.2011.07.003
- PubMed
- Google Scholar
1. Salehi S
2. Kabeer F
3. Ceglia N
4. Andronescu M
5. Williams MJ
6. Campbell KR
7. Masud T
8. Wang B
9. Biele J
10. Brimhall J
11. Gee D
12. Lee H
13. Ting J
14. Zhang AW
15. Tran H
16. O’Flanagan C
17. Dorri F
18. Rusk N
19. de Algara TR
20. Lee SR
21. Cheng BYC
22. Eirew P
23. Kono T
24. Pham J
25. Grewal D
26. Lai D
27. Moore R
28. Mungall AJ
29. Marra MA
30. McPherson A
31. Bouchard-Côté A
32. Aparicio S
33. Shah SP
34. IMAXT Consortium
(2021) Clonal fitness inferred from time-series modelling of single-cell cancer genomes
Nature 595:585–590.

https://doi.org/10.1038/s41586-021-03648-3
- PubMed
- Google Scholar
1. Sha D
2. Jin Z
3. Budczies J
4. Kluck K
5. Stenzinger A
6. Sinicrope FA
(2020) Tumor mutational burden as a predictive biomarker in solid tumors
Cancer Discovery 10:1808–1825.

https://doi.org/10.1158/2159-8290.CD-20-0522
- PubMed
- Google Scholar
1. Sottoriva A
2. Kang H
3. Ma Z
4. Graham TA
5. Salomon MP
6. Zhao J
7. Marjoram P
8. Siegmund K
9. Press MF
10. Shibata D
11. Curtis C
(2015) A Big Bang model of human colorectal tumor growth
Nature Genetics 47:209–216.

https://doi.org/10.1038/ng.3214
- PubMed
- Google Scholar
1. Stephens PJ
2. Greenman CD
3. Fu B
4. Yang F
5. Bignell GR
6. Mudie LJ
7. Pleasance ED
8. Lau KW
9. Beare D
10. Stebbings LA
11. McLaren S
12. Lin M-L
13. McBride DJ
14. Varela I
15. Nik-Zainal S
16. Leroy C
17. Jia M
18. Menzies A
19. Butler AP
20. Teague JW
21. Quail MA
22. Burton J
23. Swerdlow H
24. Carter NP
25. Morsberger LA
26. Iacobuzio-Donahue C
27. Follows GA
28. Green AR
29. Flanagan AM
30. Stratton MR
31. Futreal PA
32. Campbell PJ
(2011) Massive genomic rearrangement acquired in a single catastrophic event during cancer development
Cell 144:27–40.

https://doi.org/10.1016/j.cell.2010.11.055
- PubMed
- Google Scholar
1. Sun R
2. Hu Z
3. Curtis C
(2018) Big bang tumor growth and clonal evolution
Cold Spring Harbor Perspectives in Medicine 8:a028381.

https://doi.org/10.1101/cshperspect.a028381
- PubMed
- Google Scholar
(2018) Fitness variation in isogenic populations leads to a novel evolutionary mechanism for crossing fitness valleys
Communications Biology 1:151.

https://doi.org/10.1038/s42003-018-0160-1
- PubMed
- Google Scholar
1. Wang LH
2. Wu CF
3. Rajasekaran N
4. Shin YK
(2018) Loss of tumor suppressor gene function in human cancer: an overview
Cellular Physiology and Biochemistry 51:2647–2693.

https://doi.org/10.1159/000495956
- PubMed
- Google Scholar
1. Wang Y
2. Robinson PS
3. Coorens THH
4. Moore L
5. Lee-Six H
6. Noorani A
7. Sanders MA
8. Jung H
9. Katainen R
10. Heuschkel R
11. Brunton-Sim R
12. Weston R
13. Read D
14. Nobbs B
15. Fitzgerald RC
16. Saeb-Parsy K
17. Martincorena I
18. Campbell PJ
19. Rushbrook S
20. Zilbauer M
21. Buczacki SJA
22. Stratton MR
(2023) APOBEC mutagenesis is a common process in normal human small intestine
Nature Genetics 55:246–254.

https://doi.org/10.1038/s41588-022-01296-5
- PubMed
- Google Scholar
1. Watson HW
2. Galton F
(1875) On the Probability of the Extinction of Families
The Journal of the Anthropological Institute of Great Britain and Ireland 4:138.

https://doi.org/10.2307/2841222
- Google Scholar
1. Weinreich DM
2. Chao L
(2005) Rapid evolutionary escape by large populations from local fitness peaks is likely in nature
Evolution; International Journal of Organic Evolution 59:1175–1182.

https://doi.org/10.1111/J.0014-3820.2005.TB01769.X
- PubMed
- Google Scholar
(2009) The rate at which asexual populations cross fitness valleys
Theoretical Population Biology 75:286–300.

https://doi.org/10.1016/j.tpb.2009.02.006
- PubMed
- Google Scholar
1. Werner B
2. Case J
3. Williams MJ
4. Chkhaidze K
5. Temko D
6. Fernández-Mateos J
7. Cresswell GD
8. Nichol D
9. Cross W
10. Spiteri I
11. Huang W
12. Tomlinson IPM
13. Barnes CP
14. Graham TA
15. Sottoriva A
(2020) Measuring single cell divisions in human tissues from multi-region sequencing data
Nature Communications 11:1035.

https://doi.org/10.1038/s41467-020-14844-6
- PubMed
- Google Scholar
1. Wilke CO
(2004) The speed of adaptation in large asexual populations
Genetics 167:2045–2053.

https://doi.org/10.1534/genetics.104.027136
- PubMed
- Google Scholar
1. Williams MJ
2. Werner B
3. Heide T
4. Curtis C
5. Barnes CP
6. Sottoriva A
7. Graham TA
(2018) Quantification of subclonal selection in cancer from bulk sequencing data
Nature Genetics 50:895–903.

https://doi.org/10.1038/s41588-018-0128-6
- PubMed
- Google Scholar
1. Wright S
(1931) Evolution in mendelian populations
Genetics 16:97–159.

https://doi.org/10.1093/genetics/16.2.97
- PubMed
- Google Scholar
1. Yates LR
2. Campbell PJ
(2012) Evolution of the cancer genome
Nature Reviews. Genetics 13:795–806.

https://doi.org/10.1038/nrg3317
- PubMed
- Google Scholar
1. Zhang C-Z
2. Spektor A
3. Cornils H
4. Francis JM
5. Jackson EK
6. Liu S
7. Meyerson M
8. Pellman D
(2015) Chromothripsis from DNA damage in micronuclei
Nature 522:179–184.

https://doi.org/10.1038/nature14493
- PubMed
- Google Scholar
1. Zhao M
2. Kim P
3. Mitra R
4. Zhao J
5. Zhao Z
(2016) TSGene 2.0: an updated literature-based knowledgebase for tumor suppressor genes
Nucleic Acids Research 44:D1023–D1031.

https://doi.org/10.1093/nar/gkv1268
- PubMed
- Google Scholar

Article and author information

Author details

Christopher J Graser
1. Department of Data Science, Dana-Farber Cancer Institute, Boston, United States
2. Department of Biostatistics, Harvard T.H. Chan School of Public Health, Boston, United States
3. Department of Stem Cell and Regenerative Biology, Harvard University, Cambridge, United States
Contribution
Conceptualization, Data curation, Formal analysis, Supervision, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0009-0003-8552-7767
Wenbo Wu

Department of Data Science, Dana-Farber Cancer Institute, Boston, United States

Contribution
Investigation, Visualization, Methodology, Writing – original draft

Competing interests
No competing interests declared
Cole Christini
1. Department of Data Science, Dana-Farber Cancer Institute, Boston, United States
2. School of Computer Science, Carnegie Mellon University, Pittsburgh, United States
Contribution
Data curation, Investigation, Visualization, Methodology, Writing – original draft

Competing interests
No competing interests declared
Mia Peljak
1. Department of Pathology, New York University School of Medicine, New York, United States
2. Laura and Isaac Perlmutter Cancer Center, New York University, New York, United States
Contribution
Writing – review and editing

Competing interests
Shareholder in Vertex Pharmaceuticals and a compensated consultant for the GLG Network
Franziska Michor
1. Department of Data Science, Dana-Farber Cancer Institute, Boston, United States
2. Department of Biostatistics, Harvard T.H. Chan School of Public Health, Boston, United States
3. Department of Stem Cell and Regenerative Biology, Harvard University, Cambridge, United States
4. Center for Cancer Evolution, Dana-Farber Cancer Institute, Boston, United States
5. The Eli and Edythe L. Broad Institute, Cambridge, United States
6. Ludwig Center at Harvard, Harvard Medical School, Boston, United States
Contribution
Conceptualization, Supervision, Writing – original draft, Writing – review and editing

For correspondence
michor@jimmy.harvard.edu

Competing interests
Co-founder and consultant of Harbinger Health, a consultant for Zephyr AI, and on the board of directors of Recursion Pharmaceuticals; none of these relationships are directly or indirectly related to the content of this manuscript

"This ORCID iD identifies the author of this article:" 0000-0003-4869-8842

Funding

No external funding was received for this work.

Ethics

All data used in this study were obtained from The Cancer Genome Atlas (TCGA, https://www.cancer.gov/tcga), which provides de-identified data collected under established ethical guidelines, with informed consent from participants and approval by relevant institutional review boards.

Version history

Sent for peer review: June 30, 2025
Preprint posted: July 4, 2025
Reviewed Preprint version 1: October 2, 2025
Reviewed Preprint version 2: December 5, 2025
Version of Record published: October 9, 2026

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.108058. This DOI represents all versions, and will always resolve to the latest one.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.