A modified fluctuation-test framework characterizes the population dynamics and mutation rate of colorectal cancer persister cells

Compelling evidence shows that cancer persister cells represent a major limit to the long-term efficacy of targeted therapies. However, the phenotype and population dynamics of cancer persister cells remain unclear. We developed a quantitative framework to study persisters by combining experimental characterization and mathematical modeling. We found that, in colorectal cancer, a fraction of persisters slowly replicates. Clinically approved targeted therapies induce a switch to drug-tolerant persisters and a temporary 7- to 50-fold increase of their mutation rate, thus increasing the number of persister-derived resistant cells. These findings reveal that treatment may influence persistence and mutability in cancer cells and pinpoint inhibition of error-prone DNA polymerases as a strategy to restrict tumor recurrence.

sets of drug-response growth assays (Fig. 1a). The first, the doseresponse assay in which CRC clones were exposed to increasing concentrations of targeted therapies ( Fig. 1b and Extended Data Fig. 1a), was used to analyze growth curves defined as the number of live cells versus time and drug concentration ( Fig. 1c and Extended Data Fig. 1b) and to quantify the intertwined processes of growth and transition to persister state on drug treatment. Supplementary  Fig. 3 shows the normalization process of the dose-response assay data used to obtain growth curves (Supplementary Note). The second, the single-dose assay, consisting of 3 weeks of exposure to a constant drug concentration, highlighted a biphasic (two time-scale) killing curve ( Fig. 1d and Extended Data Fig. 1c), characterized by a rapid decline of sensitive cells followed by a slower decline 25,26 , a hallmark of the emergence of persisters in bacteria 25 . The fraction of surviving persisters displayed a slow but measurable decay in cell number (Extended Data Fig. 2), suggesting a tendency for persisters to slowly die over time.
CRC persister cells slowly replicate during drug treatment. We next sought to elucidate the dynamics of persister proliferation and death under drug treatment. Staining of CRC persisters with carboxyfluorescein succinimidyl ester (CFSE), a cell-permeable fluorescent dye allowing flow-cytometric monitoring of cell divisions 27 , and 5-ethynyl-2′-deoxyuridine (EdU), a modified thymidine analog that is efficiently incorporated into DNA during active DNA replication, revealed that a fraction between 0.2% and 2.5% of persisters slowly replicates during treatment (Extended Data Fig. 3a-c and Supplementary Note), in line with recent data 28 . We next used a live cell microscopy imaging assay (Supplementary Note). Although most CRC persisters were nonreplicating, cell division events were visible in all the CRC clones analyzed (Extended Data Fig. 3d and Supplementary Videos 1-4). The cell division was occasionally successful and viable (Extended Data Fig. 3d). We also detected cell death events after cell division and in nondividing cells (Extended Data Fig. 3e,f).

The persister phenotype is induced by targeted therapies in CRC.
To quantify cell population dynamics during drug treatment, we developed a mathematical model of the transition of CRC cells to the persister state, which we term the 'transition-to-persister' (TP) model. This model incorporates birth-death parameters and phenotypic switching in the deterministic limit (that is, neglecting fluctuations due to stochastic demographic effects, see Supplementary Note) 29,30 . Figure 2a summarizes the TP model dynamics. We exploited the model to quantify the transition rate and assess whether a subpopulation of persisters predated drug administration or whether the persister phenotype emerged on drug treatment. The TP model considers three possible fates for drug-treated cells: (1) death, (2) replication and (3) switching to persister state at a rate, λ, in the presence of the drug; it further considers the pre-existence of an arbitrary steady fraction f 0 of persisters (Fig. 2a).
The following equations define the dynamics of sensitive (X(t)) and persister (Z(t)) cells according to the TP model (under drug treatment): The initial condition that specifies the solution to equation (1) is key for quantifying to what extent the transition to persister state is induced by the drug treatment. Specifically, if the sensitive-to-persister transition is fully drug induced, then untreated populations would not contain any persisters, that is, N(t 0 ) = 0. Conversely, if some persisters pre-exist, then the initial fraction of persister cells has a finite positive value (f 0 > 0). If f 0 is very small, some persisters may pre-exist, but the transition is mainly drug induced. If f 0 is actually comparable to the fraction of residual persisters after weeks of treatment, then the transition to persistence is not drug induced.
To determine the most likely scenario, we used experimental data collected from drug-response assays (Fig. 1). Using results from the dose-response assay, we defined parameters governing the dynamics of the model over a short timeframe, such as the initial fraction of persisters f 0 and the effective growth rate of treated cells. Similarly, the single-dose assay was used to quantify model parameters that affect long-term dynamics, such as the transition rate of sensitive to persister cells (λ) and the effective death rate of persisters (D p ). By constraining model parameters from experimental data, we established which scenario would best describe the cell-based results. The inferred parameters are compatible with the values obtained by live cell microscopy assay, supporting a balance between proliferation and cell death skewed slightly toward the latter (Extended Data Fig. 3g and Supplementary Note).
On treatment, the number of cells started to decline within 1-3 d (t 0 ), depending on the initial seeding density ( Supplementary Fig. 4). The observed cell dynamics were coherent in experiments with different seeding densities once the growth curves had been scaled (both in time and in measured viability) to the maximum value reached at t = t 0 ( Supplementary Fig. 4).
The parameters of the TP model were inferred with a standard Bayesian inference framework for both cell lines (Supplementary  Table 2 and Supplementary Note). DiFi displayed slower 'dying' dynamics compared with WiDr. In light of this, in WiDr we performed a joint fit of both the dose-response and the single-dose datasets, whereas in DiFi we assessed growth curves in response to multiple doses of targeted therapies for up to 19 d, which allowed performing a model fit based on the dose-response dataset only (Supplementary Note).
We identified the best-fit TP model parameters given the experimental data, considering different values of the initial number of persisters (f 0 ). The best fit between the inferred TP model and experimental data occurs when f 0 = 0, whereas the concordance decreases when f 0 increases; we note that a value of f 0 of 10% already leads to substantial deviations from the data (Fig. 2b). Therefore, the TP model is consistent with the persister phenotype being predominantly drug induced. In addition, the model properly describes the dynamics of the single-dose assay (Fig. 2c).
To further confirm the validity of the TP model, we next focused on the Bayesian statistics of the two model parameters describing the dynamics of persisters: the transition rate λ and initial fraction of persisters f 0 . The joint posterior distribution of the Bayesian inference of these two parameters is shown in Fig. 2d and Supplementary  Fig. 5. We found that the transition rate to persistence λ estimated by the model fit does not vary when considering different values of the initial fraction of persisters f 0 (Fig. 2d). The marginalized posterior probability of f 0 peaked at zero (Fig. 2d, bottom panel), and its upper boundary is much smaller than the ratio between the persister population size (after all persister cells have emerged) and the total population size at the start of treatment. This implies that the inferred value of transition rate λ is independent from f 0 and the best concordance of the TP model to the experimental data is obtained for f 0 = 0.
Finally, to compare the scenarios f 0 = 0 and f 0 > 0, we used the Bayesian information criterion (BIC) and Akaike's information criterion (AIC), which are standard Bayesian criteria used for model selection. According to both, the TP model with f 0 = 0 is preferred over f 0 > 0 (Extended Data Fig. 4, Supplementary Table 3 and Supplementary Note). We summarize the best-fit TP model parameters in Supplementary Table 4. It is interesting that we found the transition rate of WiDr and DiFi cells to persisters to be drug dependent (Extended Data Fig. 4 and Supplementary Fig. 6). These results indicate that, even if few persisters exist in the population before drug treatment, most of them must have transitioned to the persister phenotype after drug exposure. Our finding that WiDr cells show a transition rate to persistence that increases with drug concentration could be applied to design innovative strategies to restrict persister evolution. Notably, our analysis predicts that a linear increase of drug concentration, compared with a constant dosage, might reduce the number of persisters (Extended Data Fig. 5).
Persisters distribution supports a drug-induced scenario. We then measured how the number of persisters varied across multiple wells, because the distribution of this parameter is expected to be different between a drug-induced and a pre-existing scenario 31 . We seeded DiFi cl.B6 and WiDr cl.B7 in multiple 96-well plates and quantified the distribution of persisters (residual cell viability) among >400 independent wells after 3 weeks of drug treatment (Extended Data Fig. 6a). The observed abundance distribution across wells was consistent with a Poisson distribution (Extended Data Fig. 6b), supporting a drug-induced scenario, as confirmed by computational simulations (see ref. 31 , where a similar method was used for mutational processes, and Extended Data Fig. 6c). These numerical simulations show that, provided that persisters do not exist before drug treatment, the number of persisters emerging from sensitive cells under treatment is Poisson distributed. Conversely, pre-existing persisters would be generated with a constant rate from an exponentially expanding population before treatment administration. Hence the number of pre-existing cells is not Poisson distributed, but is described by a Luria-Delbrück 15 distribution (with variance ≫ mean). We found that the final distribution of persisters across wells is a Poisson distribution, in line with emergence after drug treatment.
A fluctuation assay quantifies persisters' mutation rates. Measurement of mutation rates in the absence or presence of anticancer drugs required the development of a second model, hereafter the 'mammalian cells-Luria-Delbrück' or 'MC-LD' model. The MC-LD model is a fully stochastic birth-death branching process, describing the growth of resistant cells before and during drug treatment (Fig. 3a). We designate with μ the effective rate at which one individual (cell) develops resistance whereas μ s and μ p indicate mutation rates of sensitive (untreated) and persister cells, respectively ( Fig. 3b and Supplementary Note).
WiDr and DiFi cells were seeded in 20 96-multiwell plates each and allowed to grow for a fixed number of cell divisions in drug-free standard culture conditions (Fig. 3c); afterwards, a constant, clinically relevant drug concentration was applied (Fig. 3d). The number of wells, the initial population size in each well and the time of cell replication in the absence of drug treatment were set by theoretical considerations incorporating the population dynamics parameters that we previously measured (Supplementary Tables 1 and 4 and Supplementary Note).
In accordance with our previous work 6 , a small number of early emerging resistant colonies was detected after 3-4 weeks of treatment (Fig. 3d,e). Conversely, in the vast majority of the wells sensitive cells died, whereas drug-tolerant persisters survived, as detected by measurement of residual cell viability (Extended Data Fig. 6) 6 . After several weeks of constant treatment of the residual persister cells (Fig. 3d), late-emerging resistant colonies appeared in a subset of wells in which persisters had previously been detected (Fig. 3d,e).
We ran multiple MC-LD model simulations, with input parameters inferred with the TP model, and found that resistant clones emerging at late time points (>4 weeks of treatment) are unlikely to originate from pre-existing resistant cells ( Fig. 4a and Extended Data Fig. 7). In accordance with previous work 6 , we considered the resistant colonies that became microscopically visible within the first 4 weeks of drug treatment (early emerging resistant) as those representing pre-existing resistant cells, that is, mutant cells that emerged during the expansion phase by spontaneous mutation. We also reasoned that resistant colonies that slowly emerged after ≥10 weeks of drug treatment (late-emerging resistance) in persister-containing wells could have developed drug-resistance mutations through the adaptive mutability process that we and others have observed 12,13 ( Fig. 4a and Extended Data Fig. 7).
As in a standard fluctuation test, the mutation rate can be inferred from the observed fraction of wells containing resistant cells. In the model, this fraction corresponds to the expected probability of observing a resistant clone in a well in a given time interval [0,T]. To compute this probability in the MC-LD model, we assumed that resistant cells divide with rate b and die with rate d, just like untreated cells. Supplementary Fig. 7 supports the stability of the inferred values of the mutation rates against variation of the division rate of resistant cells. As a result of reproductive Consequently, the probability of having at least one mutant is given by: To quantify the spontaneous mutation rate of cancer cells before drug administration, we focused on the resistant cells established by

Fig. 3 | A modified Luria-Delbrück fluctuation test to measure mutation rates in mammalian cells. a,
The modified fluctuation test, based on the inferred population dynamics and the MC-LD model, allows estimation of spontaneous (μ s ) and persister (μ p ) mutation rates. b, Schematic representation of cell population dynamics of CRC cells during the fluctuation test. During the initial expansion in the absence of drug treatment, CRC cells mutate with the spontaneous mutation rate (μ s ). When cells are exposed to targeted therapies, pre-existing resistant cells are selected by the drug and give rise to early emerging resistant colonies (red dashed line), whereas sensitive cells start to decline in number (black solid line) and switch to the persister state (dark-yellow solid line). Resistant cells derive from persisters with a mutation rate μ p and give rise to late-emerging resistant colonies (blue dashed line). c, Schematic representation of the experimental design underlying the fluctuation assay. WiDr and DiFi cells were seeded in 20 96-multiwell plates, for a total of 1,920 wells, and allowed to expand in the absence of drug for about 8 generations (reaching ~20,000 cells per well). After the expansion, all the wells were treated with targeted therapy (100 μg ml −1 of cetuximab for DiFi and 1 μM dabrafenib + 50 μg ml −1 of cetuximab for WiDr). d, Two sets of resistant clones were identified during the MC-LD experimental assay: the early emerging resistant clones grown after 3-4 weeks (stage 1) and the late-emerging resistant clones arising after >10 weeks (stage 2) of constant drug treatment. Scale bar, 100 μm. e, Bar graphs listing the number of resistant clones counted at the indicated timepoints during the MC-LD experiment for each CRC clone. Red bars indicate early emerging resistant clones (appearing in the first 4 weeks of drug treatment); blue bars indicate late-emerging resistant clones (appearing after ≥10 weeks of drug treatment). Results of two independent biological replicates for each clone are shown.
the time T treat before treatment administration. The expected number of resistant cells that emerged from sensitive cells in this time interval reads: where μ s is the mutation rate of sensitive cells.
To quantify the mutation rate of persister cells μ p , we consider resistant cells that emerged by the time T since the beginning of the drug treatment. The expected number of resistant cells that emerged from persister cells reads: We emphasize that equations (2)-(4) are connected to the solution of the TP model, equation (1). Hence, the solution of the MC-LD model is defined in terms of the same parameters that were estimated with the TP model (Supplementary Note).
We used this solution of the MC-LD model to derive estimators of mutation rates of sensitive cells μ s (encompassing the fraction of wells with early emerging resistant cells) and of persister cells μ p (corresponding to the fraction of wells with late-emerging resistant clones).

Persisters show an increased mutation rate under treatment.
Data collected with two-step MC-LD fluctuation tests for each clone allowed inferring mutation rates of sensitive (μ s ) and persister (μ p ) cells. We conservatively evaluated the mutational processes as chronological (measured in mutations per day) rather than replicative (mutations per generation). This choice is safe, as the ratio between replicative mutation rates of cells displaying the two phenotypes must always be higher than for chronological rates, because (beyond any uncertainty) measured cell division in persister cells was very low compared with that of untreated cells.
We found that mutation rates were increased by a factor of 7-to 50-fold in cells that survived and tolerated for several weeks doses of targeted therapies that were lethal for most of the parental population ( Fig. 4b and Supplementary Table 5). This result was consistent across multiple biological replicates, both in DiFi and WiDr cells and in response to clinically relevant concentrations of either EGFR blockade or EGFR/BRAF concomitant inhibition, respectively ( Fig. 4b and Supplementary Table 5).
To further validate the consistency of the mutation rate inference based on the MC-LD model, we ran multiple simulated replicates of the experiment, using a set of sensitive (μ s ) and persister (μ p ) mutation rates, and we then used the MC-LD estimators on the synthetic data. Figure 4c compares boxplots of the estimated mutation rates across replicates of simulated experiments, with the actual values of mutation rates used as inputs to the simulations. The agreement between these values validates our estimates.
We next assessed whether and to what extent the inferred value of the mutation rate is affected by the presence of different numbers of pre-existing persister cells f 0 using our estimators within a Bayesian framework (Fig. 4d). This approach returns the mutation rate, considering a range of realistic values of f 0 and their probability. We obtained the fold increase of the mutation rate of persister cells as a function of f 0 , in the entire range of values that are compatible with the dynamics observed in the growth curve assays experimentally assessed in Figs. 1 and 2. Figure 4d summarizes the results of this inference. We found that, considering all representative values of f 0 that are compatible with our experimental data, the increase of mutation rate in persister cells remains strongly supported.
To corroborate these results we replicated the full set of experiments and ran the analysis pipeline for two additional clones, one for each cell model (WiDr cl.B5 and DiFi cl.B3), thereby confirm-

Fig. 4 | Quantification of mutation rates in persister cells. a,
Simulated data for the assay described in Fig. 3. The experimentally measured MC-LD model parameters and the model-derived estimators of mutation rate, for sensitive and persister cells (dark-yellow solid lines), were used to simulate the time of appearance of pre-existing (red) and persister-derived (blue) resistant cells. b, Quantification of mutation rates for sensitive (red) and persister (blue) cells in the MC-LD experiment. The indicated cell models were seeded and treated as described in Fig. 3. Mutation rates were calculated from the experimental data, based on population parameters and the number of pre-existing (early emerging) and persister-derived (late-emerging) resistant clones as described in Fig. 3. Results represent inferred mutation rates (for each clone of sensitive and persister cells) with bar plots showing the mean of the posterior distributions of the mutation rates (n = 2 biologically independent experiments). Here, the bar chart is used as a graphic representation of inferred mutation rates (see Supplementary Table  5 for the corresponding numerical values). c, Validation of mutation rates estimator with model simulations. The boxplots represent the distribution of the estimated mutation rates for n = 100 independent simulations of the entire experiment using the parameters reported in Supplementary Tables 1 and 4 ing our findings and excluding a clonal bias effect (Extended Data Figs. 8 and 9). Molecular profiling of persister-derived resistant clones isolated from the fluctuation assays revealed acquisition of single-nucleotide variations (SNVs) or copy-number alterations (CNAs) in genes involved in the RAS-MEK pathway, which are known drivers of resistance to anti-EGFR/anti-BRAF inhibitors in CRC 21,34,35 (Supplementary Fig. 8 and Supplementary Table 6).
We propose a quantitative model for the evolutionary dynamics of CRC cells exposed to targeted therapies (Fig. 5a). Untreated cancer cells replicate and spontaneously acquire mutations that can confer resistance to targeted therapies (pre-existing resistant mutations) at a replicative mutation rate μ s . However, when cells are exposed to targeted therapies, most of them quickly die, whereas a subset of parental cells switch to a long-lasting surviving persister state at a rate λ and in a drug-induced fashion. Previous and current findings indicate that persister cells, under constant exposure to lethal doses of drugs, initiate a stress response that affects DNA replication fidelity 12,13 , thus leading to a measurable increase of their mutation rate (μ p ), therefore raising the probability that alterations conferring drug resistance could occur.
Inhibition of mutagenic REV1 extends the efficacy of targeted therapy. We previously reported that, in response to drug treatment, cancer cells switch from high-to low-fidelity DNA replication through downregulation of DNA repair genes and upregulation of specialized error-prone DNA polymerases 13 . This, in turn, could foster the temporary increase of the mutation rate observed in persister cells as quantitatively measured in the present study. Among the DNA polymerases that are upregulated in cancer cells on targeted therapy 13 , REV1 carries out translesion synthesis (TLS), a mutagenic process that allows cells to tolerate DNA damage by bypassing lesions that block normal DNA replication, resulting in the introduction of mutations 36,37 . Interfering with TLS using a REV1 inhibitor has been shown to enhance chemotherapy efficacy and suppress tumor growth both in vitro and in vivo 38,39 . Based on the above, we hypothesized that inhibition of mutagenic TLS would probably increase the cytotoxic effects of targeted therapy-induced DNA damage, therefore delaying the acquisition of resistance during adaptive mutability.
To assess this possibility, we performed a time-to-progression (TTP) assay, an approach that we previously established 40 to monitor the development of secondary resistance in cancer cells. DiFi and WiDr CRC cells, as well as the BRAF V600E-mutated cell line (JVE207), were treated with a MAPK (mitogen-activated protein kinase) pathway inhibitor, the REV1 inhibitor or their combination. Pharmacological blockade of REV1 remarkably delayed or prevented the development of secondary resistance to EGFR/BRAF inhibitors ( Fig. 5b and Extended Data Fig. 10).

Discussion
We present an experimental framework that integrates biological assays and mathematical modeling to investigate population dynamics of cell lines exposed to environmental perturbations. The MC-LD assay allows quantitative comparisons of spontaneous and drug-induced mutation rates and could, in principle, be applied to measure whether and how a wide range of environmental conditions affects persister phenotype and mammalian cells' mutation rates in a considerable number of biological systems.
An important caveat of our technique is that fluctuation tests measure 'phenotypic' mutation rates, that is, rates of conversion to a phenotype (here, resistance to targeted therapies) that could result from different mutational routes, including SNVs and CNAs, both of which we found to drive resistance in our persister-derived resistant clones. Nevertheless, our approach has the advantage of bypassing several hurdles associated with sequencing-based measurement of mutation rates.
Although we and others have recently shown that adaptive mutability fosters the acquisition of secondary resistance by increasing genomic instability in surviving persister cells 12,13 , the lack of models to quantitatively characterize the behavior of persisters under treatment has so far prevented reliable quantification of persisters' mutation rates. Although limited to cell lines, our controlled two-step fluctuation assay overcomes these issues. The results could be used to infer features of more complex systems (such as patient samples), where mathematical models are postulated and cannot be analogously validated.
Our results indicate that drug-induced sensitive-to-persister transition is a predominant path to the development of this phenotype. This is in line with recent evidence of a chemotherapy-induced persister state in CRC 41 . Although our results are coherently explained by the existence of a phenotypic switch of sensitive to persister cells, no direct observation of the switching is yet available. Alternative models whereby slower-proliferating tolerant cells generate faster-proliferating nontolerant phenotypes would also give rise to a biphasic killing curve 42,43 . In our framework, this would correspond to the case where f 0 ≫ 0, which is ruled out by our analysis. Hence, the interpretation linking the phenotypic switch to persistence with treatment appears to be the most likely scenario.
Importantly, even a small subset of pre-existing persisters does not affect our findings that mutation rates of cancer cells remarkably increase under treatment. Persister-derived resistant clones keep emerging after several weeks of continuous drug treatment. In the absence of an increased mutation rate, it would typically take (in a conservative estimate) >100-1,000 weeks for the cells in a single well to develop resistance based on the mutation rate of sensitive (untreated) cells. Equivalently a fluctuation assay performed on 2 × 10 5 -2 × 10 6 wells would be required to observe a few resistant clones after 10 weeks. In addition, the contribution of sensitive cells to resistance is exhausted after a few weeks of treatment, because we show that they become extinct within a few days (Fig. 1d). The evidence of active cell cycle progression, alongside an increase of mutagenic rate in persister cells, further supports previous findings of ongoing adaptive mutagenesis fostering acquisition of resistance 13 .
A recent study concluded that the impact of adaptive mutability on the mutational load in the protein-coding genome is small 44 . This in line with previous findings showing that the tumor mutational burden is not strongly increased in cells that acquired therapy resistance 12,13 . There are many confounding factors in these not precisely controlled systems. Indeed, the adaptive mutability phenotype is probably restricted in time (that is, when the cells are maladapted to the new environment 13 ), confined to a small subpopulation of cells and masked by the outgrowth of pre-existing resistant cells. Bulk analysis on tumor samples at relapse cannot disentangle these factors. Our characterization of persisters' dynamics and increased mutation rate have potential clinical relevance. First, the finding that higher drug concentrations induce an increased death rate of sensitive cells and an increased transition to persistence, a reservoir for the emergence of resistance, provides a rationale for therapeutic strategies to impair the emergence of persistence. Second, the involvement of error-prone DNA polymerases during adaptive mutability 13 offers opportunities for nonobvious combinatorial strategies to restrict drug resistance. Indeed, we show that inhibition of mutagenic TLS effectively delays the acquisition of secondary resistance.
Our methodology infers that clinically approved anticancer therapies can induce a temporary increase in the mutation rate of CRC cells. Our framework can be used to systematically measure mutation rates in mammalian cells exposed to a wide range of environmental stressors and to define drug combinations to restrict the emergence of therapeutic resistance.

Online content
Any methods, additional references, Nature Research reporting summaries, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41588-022-01105-z. Single-dose growth curve assay. DiFi and WiDr CRC cell clones were seeded in multiple 96-multiwell plates at 1,000 and 500 cells per well, respectively. Cells were allowed to expand for a fixed number of generations until a population size of 10,000-20,000 cells per well was reached. At that point, treatment was added (100 µg ml −1 of cetuximab for DiFi and 1 µM dabrafenib + 50 µg ml −1 of cetuximab for WiDr). Plates were then incubated at 37 °C in 5% CO 2 and cell viability was assessed at the indicated timepoints by measuring ATP content with the Cell Titer-Glo Luminescent Cell Viability assay (Promega), using the Tecan Spark 10M plate reader and the Tecan SparkControl Magellan software (v.2.2), over 22 and 32 d of constant treatment (for WiDr and DiFi, respectively). Medium and treatment were renewed once a week. To test the effect of different seeding densities on the residual viability assayed, each clone was seeded at different densities ((3-20) × 10 3 cells per well) in complete medium. The next day, treatment was added (100 µg ml −1 of cetuximab for DiFi and 1 µM dabrafenib + 50 µg ml −1 of cetuximab for WiDr) and viability was assessed at the indicated timepoints by measuring ATP content.
Staining with CFSE. CRC clones were seeded at 2.5 × 10 5 (WiDr) and 6.5 × 10 5 (DiFi) cells in multiple 10-cm dishes. The next day, untreated cells were stained with CellTrace CFSE Cell Proliferation Kit (Invitrogen) according to the manufacturer's instructions. At the indicated timepoints, starting from the day after staining (T0), cells were collected and resuspended in 1 ml of phosphate-buffered saline (PBS) with Zombie Violet dye 1,000× (BioLegend) to exclude dead cells. Cells were then analyzed by flow cytometry. For persister proliferation analysis, CRC clones were seeded at 2 × 10 4 cells per well in several 24-multiwell plates. The next day, cells were treated with 100 µg ml −1 of cetuximab (for DiFi) or 1 µM dabrafenib + 50 µg ml Staining with EdU. DiFi and WiDr clones were plated on several glass coverslips at 5 × 10 4 and 4 × 10 4 cells per coverslip, respectively. The next day, cells were treated with 100 µg ml −1 of cetuximab (for DiFi) or 1 µM dabrafenib + 50 µg ml −1 of cetuximab (for WiDr) and incubated at 37 °C in 5% CO 2 for 14 d (renewing treatment after 1 week), until a population of persister cells emerged on each coverslip. Then, at indicated timepoints (renewing medium and treatment once every week throughout the experiment), residual cells were stained with the Click-iT EdU Cell Proliferation Kit for Imaging (Invitrogen), according to the manufacturer's instruction. Briefly, cells were incubated with 10 µM EdU for 4 h. After that, cells were fixed in 4% paraformaldehyde for 20 min at room temperature and permeabilized with 0.5% Triton X-100 in PBS for 20 min at room temperature. Coverslips were then incubated in Click-iT reaction cocktail for 30 min, followed by nuclei staining with DAPI and F-actin staining with Alexa Fluor-555 phalloidin (50 μg ml −1 ). Slides were then mounted using the fluorescence mounting medium (Dako  48 . For each cell line pre-and post-treatment, gene copy number (GCN) was computed as follows: first the median read depth of all genomic regions was calculated; next, for each gene the median read depth was obtained and then divided by the former value. For each gene, its GCN in the preand post-treatment samples and the corresponding CNV (ratio between matched GCNs) were reported. DNAcopy R module was performed to cluster CNV using a circular binary segmentation algorithm.
TTP assay. TTP assays were conducted as previously described 40 . Briefly, 5 million cells (for WiDr and DiFi cells) and 4.5 million cells (for JVE207 cells) were plated for each treatment condition. Then, cells were treated with MAPK pathway inhibitors (dabrafenib 1 µM + cetuximab 30 µg ml −1 for WiDr and JVE207, cetuximab 50 µg ml −1 for DiFi), REV1 inhibitor (2 µM) or their combination, in parallel. Medium and treatment(s) were renewed weekly. Cells were counted each week; counts as 0 represent timepoints in which cells were too few and only medium and drug refreshments were done.
Materials availability. The CRC cell clones generated in the present study are available through A. Bardelli (Department of Oncology, University of Torino) under a material transfer agreement.
Theoretical modeling and computational analyses. In the present study, we developed and used two distinct mathematical models to investigate the dynamics of cell populations. The first model describes the transition-to-persister state (TP model), and is a birth-death model with phenotypic switching, which we explored in the deterministic limit. We considered four different model variants and compared them with experimental data to infer the most likely scenario for the sensitive-to-persister transition. The second model, which we named the 'MC-LD' model, is a fully stochastic birth-death branching process that includes the mutational processes of sensitive (untreated) and persister cells (under treatment).
To measure the mutation rate, stochastic fluctuations cannot be neglected. We simulated individual trajectories of the Markov process underlying the evolution of the MC-LD model by a coarse-grained version of the Gillespie algorithm 49 , which groups together all stochastic events happening in discrete time intervals of fixed duration ∆t.
For the inference of the birth-death rates b and d, we used the data on growth rates of CRC clones before drug treatment. Our inference scheme is summarized in Supplementary Fig. 2a. The parameters of the TP model were inferred using a Bayesian framework and data from the single-dose and dose-response assays.
Posterior distributions of the model parameters were sampled using a Hamiltonian Monte Carlo algorithm (Python 3, package pymc3, NUTS sampler) 50 . TP model variants were compared by means of the standard BIC and the AIC. To infer mutation rates for the MC-LD model, we computed an approximate analytical expression for the probability of the emergence of one mutant in an expanding population of cells in a given time interval [0,T], and used it to derive estimators for the emergence of mutations before and during treatment administration (from persisters). The mutation rate of persister cells was inferred with a Bayesian framework, to account for the uncertainty of the value of the initial fraction of persister cells, f 0 .
All the details on the theoretical/computational protocols are provided in Supplementary Note.   Fig. 3 | A fraction of persister cells slowly replicate during drug treatment. a, Distribution of the CFSE signal (Fitc) measured by flow cytometry is reported for the indicated timepoints. In each distribution the number of cells retained after the gating process was ≥ 3000. WiDr and DiFi cells were grown in standard conditions (untreated) or treated with 1 μM dabrafenib + 50 μg/ml cetuximab or 100 μg/ml cetuximab, respectively, for 2 weeks until the emergence of surviving persister cells. Both untreated (sensitive) and persister cells were stained with CFSE to quantify cell division and fluorescent signal (Fitc) was analyzed by flow cytometry at indicated time points. One representative experiment (n = 2 biologically independent experiments for untreated sensitive cells; n = 3 biologically independent experiments for persisters) is reported. b, Quantification of EdU positive cells at indicated time points. Surviving persisters emerged after 2 weeks of treatment of DiFi and WiDr clones with targeted therapies were labeled with EdU for 4 hours, then fixed and analyzed by fluorescence microscopy. Results represents mean ± SD of two independent experiments, each represented by two technical replicates. c, Representative images of EdU positive persister cells for each CRC clone analyzed of 2 biologically independent experiments. Green: EdU; Blue: DAPI; Red: Phalloidin. Scale bar 50μm. d-f, Representative snapshots of cell division events and apoptotic events observed in the timelapse microscopy assay of 3 biologically independent experiments. Surviving persisters emerged after 2 weeks of treatment of DiFi and WiDr clones were stained with Nucblue® (a dye for nuclei, in blue) and Caspase Cell Event TM (a live marker for the activation of Caspase3/7, in green) while maintaining drug pressure, and monitored for 5 days under an inverted widefield microscope.  Supplementary Table 4). b, In a modified drug delivery strategy, the drug concentration increases linearly over time. This strategy is constrained to have the same average drug concentration as in panel a, meaning that the total amount of drug delivered in 10 days is the same in the two strategies. In this case, at the beginning of the in-silico treatment, the delivered drug is below the MIC value. c, Expected number of persister cells predicted by the TP model (equation 1) for WiDr cl. B7 cells, corresponding to the two alternative drug delivery strategies: constant value (gray line) and linear increase (orange line). d, Relative fold change of the number of persister cells emerged with the drug strategy b vs drug strategy a, evaluated as a function of the average drug concentration. Fig. 6 | Distribution of persisters abundance across wells is compatible with the tP model expectation of f0=0. a, CRC clones were seeded in multiple 96-multiwell plates and allowed to expand for about 8 cell divisions in the absence of drug, and then treated with 1 μM dabrafenib + 50 μg/ml cetuximab (WiDr) or 100 μg/ml cetuximab (DiFi) for 3 weeks. Cell viability was determined by ATP assay. Wells containing rapidly proliferating pre-existing resistant colonies are marked in red, while the remaining wells contained a small number of viable cells, which we identified as drug-tolerant persisters (indicated in light-blue). Each bar represents one well. One representative experiment of two biologically independent replicates is reported. b, Cumulative distribution of persisters cell viability across wells (light-blue bars) is compatible with a Poisson cumulative distribution (magenta solid line, see Supplementary methods for details on the fitting procedure). One representative experiment of two biologically independent replicates is reported. c, The simulation of a stochastic model for the transition to persistence shows that in the case of f 0 = 0 (that is, in the drug-induced scenario) the number of persister cells per well is expected to have a Poisson distribution. Simulations were performed using the TP model fit values as input parameters, and setting the initial population size to 15000 cells/well. Simulations were stopped at 21 days; the distribution was evaluated using n = 1000 independent simulations.