Inheritable cell-states shape drug-persister correlations and population dynamics in cancer cells

Inheritable cell-states shape drug-persister correlations and1 population dynamics in cancer cells2 Anton Iyer1, Adri´ an E. Granada2 and Shaon Chakrabarti 1,* 1Simons Centre for the Study of Living Machines, National Centre for Biological Sciences, Bangalore 2 Charit´ e Comprehensive Cancer Center, Berlin *Email for correspondence: [email protected] 3 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint Abstract4 Drug tolerant persisters (DTPs) drive cancer therapy resistance by temporarily evading drug action,5 allowing multiple routes to eventual permanent resistance. Despite clear evidence for DTPs, the6 timing of their emergence, proliferative nature, and how their population dynamics arise from mea-7 sured single-cell kinetics remain poorly understood. Here we use time-lapse microscopy data from8 two cancer cell lines, integrating single-cell and population measurements, to develop a quantitative9 description of drug persistence. Contrary to the expectation that increasing levels of genotoxic stress10 should lead to slower times to division and faster times to death, we observe minor changes in the11 single-cell intermitotic and death time distributions upon increasing cisplatin concentration. Yet,12 population decay rates increase 3-fold, suggesting a surprising independence of the overall dynamics13 from the measured birth and death rates. To explain this phenomenon, we argue that the observed14 lineage correlations and concentration-dependent decay rates imply cell-state dependent fate choices15 made both pre and post-cisplatin as opposed to just post-drug birth/death rate-based competitive16 fate choices. We demonstrate that these cell-states, present in the ancestors of DTP and sensi-17 tive cells, exhibit no difference in cycling speed and are inherited across 2-3 cellular generations. A18 stochastic model implementing these rules simultaneously recapitulates the observed decay rates and19 cell-fate correlations, also explaining how pre-drug fate decisions are consistent with barcoding ex-20 periments where barcode diversity remains unchanged after drug administration. Our results provide21 a powerful perspective on drug tolerance based on general arguments, without requiring knowledge22 of the underlying molecular architecture of the heterogeneous cell states.23 Introduction24 Non-genetic heterogeneity amongst isogenic cells gives rise to drug-tolerant persisters or “DTP”s –25 cells that temporarily survive drug treatment, often providing a reservoir for future genetic changes26 driving permanent resistance [1, 2, 3, 4, 5, 6, 7, 8]. Originally discovered in the context of bacterial-27 cell survival in the presence of antibiotics [9], the existence of persisters was later suggested for cancer28 as well, where a fractional-killing effect was demonstrated to limit curability in leukemia [10]. A large29 body of literature has since accumulated, suggesting interesting parallels between the phenomenon30 of persistence in bacterial and cancer cells, including a characteristic bi-phasic exponential decay in31 1 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint the population dynamics after drug administration and the origin of persisters from a small fraction32 of quiescent or slow-cycling cells in drug-naive conditions [8, 11, 12, 13, 14]. Furthermore, there33 seems to be no unique non-genetic cell-state that gives rise to persistence either in bacteria or in34 cancer cells – a wide range of regulatory networks such as metabolic, extra-cellular matrix and EMT35 pathways are likely involved depending on the specifics of the cell-type, the drug being administered36 and cellular interactions with the micro-environment.37 While there have been a lot of recent advances in characterizing DTP’s and their origins, many38 questions remain unresolved and debated, particularly in the context of cancer drug tolerance [15].39 For instance, the timing of cellular decisions leading to DTPs remains unclear. Based on end-fate40 (survival or death) correlations on lineages, some previous studies concluded that heterogeneous41 cell states existing before drug administration (“pre-DTPs”) pre-determine cell fate outcomes post42 drug treatment [2, 16, 17, 18]. Conceptually similar conclusions were later drawn from elegant43 barcoding techniques that additionally gave insights into the gene expression states that might44 determine these pre-DTP states [19, 20, 21]. It was also demonstrated using theoretical models of45 transcriptional bursting that it is indeed possible to generate rare cell-states where fluctuations in46 gene expression persist for a few generations, leading to cell survival in the face of drug treatment [22].47 However, other studies have argued that fate decisions are drug-induced by combining mathematical48 modeling of cancer population dynamics, modified Luria-Delbruck fluctuation tests and barcoding49 approaches [23, 7, 24]. Similarly, the question of whether persister cells always arise from lineages50 with intrinsically slower division kinetics remains debated. Original studies in bacteria attributed51 persisters to the presence of slowly cycling sub-populations in drug naive conditions[12, 13]. Most52 studies in cancer cells have also suggested the same phenomenon, resulting in the terms “cancer53 persisters” and “quiescence/dormant” being used almost interchangeably [25, 26]. These studies54 suggest the presence of a fitness cost to non-genetic persistence in the same vein as is often ascribed55 to genetic resistance mechanisms [27]. However, emerging evidence suggests that the ancestors of56 cancer persisters before drug treatment might also be proficient in cell-division, therefore challenging57 the wide-spread notion that drug tolerant persisters necessarily emerge from quiescent subpopulations58 of cancer cells[20, 28, 29].59 In this study we explore the nature, timing and kinetics of these cell fate decisions by investigating60 2 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint how the observed population dynamics of cancer cells quantitatively emerges from the measured61 single-cell kinetics. Integrating single-cell and population dynamics measurements from time-lapse62 microscopy datasets, we first demonstrate the surprising result that increasing cisplatin concentration63 does not affect cellular birth or death rates, unlike what is widely assumed [30, 31, 32, 33, 34, 35]. As64 a consequence, the ubiquitous exponential growth models which assume stochastic competition of cell65 fates (birth-death process and age-structured population models), fail to explain experimental data66 from two cisplatin-treated cell types, HCT116 and U2OS. Instead, we argue based on observed lineage67 correlations that pre-existing cell states largely determine fate decisions, and are inherited across68 multiple cellular generations before drug addition. However, some degree of fate determination must69 also occur during the time of drug administration, which is necessary to explain drug concentration-70 dependent effects on the population dynamics. Furthermore, we demonstrate that cells in states that71 are primed for death or survival do not have fitness differences, and are both inherited across 2-372 generations before the drug is administered. Finally, we show how early fate decisions are consistent73 with recent barcoding experiments where no change in barcode diversity was observed after drug74 administration. Our arguments are based on general principles independent of the molecular details75 of the underlying cellular states, providing a powerful approach to deciphering the dynamics of drug76 persistence in cancer.77 Methods78 Mixture models to investigate the presence of multiple cellular popula-79 tions before cisplatin, with different cycling kinetics.80 While full lineages were tracked in the HCT116 dataset [18], only one cell per sister-pair was tracked81 in the U2OS dataset [36]. Hence for the latter, we resorted to using a mixture model to check whether82 multiple subpopulations with distinct proliferation rates could be identified from the pre-cisplatin83 cells. Since previous work has identified the Exponentially Modified Gaussian (EMG) as a good84 model of the IMT distribution [18], we used a mixture of EMGs for this purpose. The EMG is a85 convolution of an exponential and Gaussian distributions, the former parameterized by λ and the86 latter by µ and σ. We used a mixture of two EMGs to fit the IMT distribution data and compared87 3 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint with the fit to a single EMG using the Akaike Information Criterion (AIC). Details of the likelihood88 function and inferred parameters can be found in SI Section 3.89 Predicting population dynamics from single-cell measurements.90 Age structured population growth models link the net proliferation rate of the population to single-91 cell division and death time distributions [37, 38, 39, 40]. They relate the exponential proliferation92 rate of the population to the age dependent division and death rates (hazards) via the Euler-Lotka93 equation:94 1 = ˆ ∞ 0 ϕ(x)e−γxdx, (1) where ϕ(x) = b(x)e− ´ x 0 hnet(y)dy, hnet(y) = b(y) + µ(y) is the net hazard, and b(y) and µ(y) are95 the age dependent birth and death hazards respectively. γ is the unique root to Equation 1 that96 gives the population growth rate. The hazard functions are time-dependent generalizations of the97 rate parameter in exponential distributions. The parameters of the birth and death hazards can be98 inferred from the experimentally measured single-cell inter-mitotic time (IMT) and apoptosis time99 (AT) distributions respectively. These can then be used to solve the integral Equation 1 to obtain100 the population growth/decay rate γ. Detailed descriptions of parameter inference and calculation of101 the hazard and proliferation rates are provided in SI section 8.102 Stochastic simulations of age-structured population dynamics – Model103 M0.104 Model M0 corresponds to the conventional age-structured model of population dynamics which105 generalizes constant division and death rates to age-dependent rates, giving rise to exponential106 growth at long times. Stochastic simulations of this model were performed in R for both pre-drug107 (only division) and post-drug (division and death) scenarios[38]. Single-cell distributions of division108 and death times were taken from the experimental datasets and converted to age-dependent birth109 and death rates (hazards; see SI section 8 for details). Each new cell in a lineage was assigned an age110 0 at birth and for each following time-step, the probabilities of dividing versus nothing happening (or111 4 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint dividing, dying or nothing happening after drug) were computed based on the birth and death rates.112 This tracking was done for each extant cell. At each division event, two new cells were stored as113 nodes in the extant cells list and the graph that was uniquely defined for each lineage tree. Attributes114 like ‘birth time’, ‘state’, ‘fate’, etc. were stored in these newly created nodes for future reference.115 Further details can be found in SI Section 9.116 Competing risks analysis of exponentially distributed division and death117 time distributions.118 Since cells exposed to cisplatin can undergo either division or death, these events are mutually119 exclusive and hence give rise to ‘competing risks’ for each cell. A well known problem in statistics,120 the competing risks effect leads to significant biases in the experimentally measured times to division121 and death after cisplatin addition [18]. Since only the event with a shorter time will ultimately122 be observed for any single-cell, all the measured division and death times will be biased towards123 lower values, thereby hiding the actual underlying distributions of division and death times. These124 underlying distributions therefore are not directly measurable, but need to be inferred from the125 measured biased division and death times. A simple analytical calculation using exponential waiting126 times for the birth and death events illustrates the competing risks problem, and demonstrates how127 inference of the underlying distributions can be carried out (see SI Section 7 for details).128 A Markov Chain Monte Carlo approach to competing risks analysis for129 HCT116 and U2OS cells.130 While exponential distributions allow for an analytical solution of the competing risks inference prob-131 lem, non-exponential distributions along with the pre and post-drug scenarios in the HCT116/U2OS132 experiments complicate the inference problem and preclude analytic solutions. For the post-drug133 scenario we needed to find the underlying distributions of IMT and AT as mentioned in the com-134 peting risks section, which are modeled as EMGs having three parameters each ( µ, σ and λ). We135 performed the inference of these distributions from the biased IMT and AT values obtained from136 the experiment using a Metropolis-Hastings Markov Chain Monte Carlo approach. In this method137 we sampled a total of 8 parameters from their posteriors. Six of these parameters describe the two138 5 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint distributions, IMT and AT (µ, σ and λ for each distribution). The remaining 2 parameters describe139 the probability of the cell to enter a quiescent stage and the delay time for the drug to act. To find140 the posteriors of these parameters, we formulated likelihoods for all the division and death events141 that are possible post-drug and assumed a uniform prior to obtain the posterior using Metropolis142 Hastings MCMC. Details of the inference method are provided in SI section 8 .143 Beyond age-structured populations: stochastic simulations of Models144 M1, M2 and M3.145 Model M1 : In Model M1, we defined the cells to exist in one of two states, one being more susceptible146 to die (sensitive, S’) while the other is more likely to persist (pre-DTP, P’) in the face of drug147 treatment. The cells were allowed to interconvert between these states in drug naive conditions.148 Post cisplatin, sensitive cells and pre-DTP cells were assigned fates to die and divide respectively.149 The times taken for the division or death events were sampled from the respective distributions150 measured in the experiments. This model could recapitulate the experimentally observed lineage151 correlations, but not the drug-concentration dependent decay rates.152 Model M2 : In this model, the cells possessed no ‘states’ before cisplatin so that on its administration,153 the cells were randomly assigned fates with probabilities that depended on drug concentration. While154 this model could explain the population decay rates observed in the experiment, it failed in explaining155 the correlations seen in the final post-cisplatin fates of lineage related cells.156 Model M3 : This model is essentially a combination of Models M1 and M2, where we simulated cells157 in two states S’ (sensitive) and P’ (pre-DTP) prior to drug addition with inter-state transitions.158 After the drug, the cells were assigned fates from 3 possibilities – death, survival or division, with159 probabilities conditioned on their state at the time of drug treatment. The state-dependent fate160 probabilities were chosen from a fate-matrix F. If a cell was assigned ‘division’ as the fate, the two161 resultant daughters maintained the original state of the mother cell at time of drug addition, and162 the fate assignment process was repeated using the same matrix F. This model could quantitatively163 recapitulate both the lineage correlations as well as the drug concentration dependent decay rates164 measured in the experiments. Further details are provided in SI Sections 9 and 10.165 6 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint Results166 Subpopulations with distinct proliferation rates cannot be detected in167 HCT116 or U2OS cells.168 To obtain accurate single-cell intermitotic times (IMTs), apoptosis times (ATs) and cell lineages, we169 obtained datasets from two previous studies on colon cancer HCT116 [18] (Fig. 1a-c) and osteosar-170 coma U2OS [36] (Fig. 1d-f) cell lines. These studies performed five days of single-cell time lapse171 microscopy, tracking cells in two days of drug-free medium followed by three days of cisplatin treat-172 ment, with snapshots taken every 30 mins. While HCT116 were exposed to a single concentration of173 cisplatin approximately near its IC50 value [18], U2OS cells were treated with four concentrations of174 cisplatin, 0, 7, 10 and 13 µM corresponding to Control, Low, Medium and High labels [36]. Besides175 measurements of IMT and AT of single cells, these datasets allow simultaneous extraction of the176 population dynamics by counting the number of surviving cells in each time frame. Furthermore,177 the HCT116 dataset also provides information on lineage relationships between the various single178 cells. In both datasets there were a small fraction of cells which never divided from the beginning of179 tracking, and which did not show any bias towards end fate (survival or death post cisplatin). We180 excluded these cells from all further analysis. SI Section 1 provides further details of both datasets,181 number of cells included/excluded and the fraction of end-fates observed.182 7 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint Figure 1: Pre-cisplatin heterogeneity in inter-mitotic times arises from a single population of cycling cells. (a) HCT116 tracking records the trajectories of all lineages. The dashed vertical line shows time of administration of cisplatin. Blue lines track back to ancestors of surviving (persister) cells while the red lines track back to ancestors of cells that die (sensitive cells). (b) log of population versus time for HCT116 cells show a distinct bi-phasic exponential decay after cisplatin (denoted by the vertical line). (c) Boxplot comparing IMTs of distinct drug naive ancestor cells belonging to lineages of persister or sensitive phenotypes. (d) U2OS tracking records the trajectories of a randomly chosen daughter cell (forward lineages). (e) log of population versus time for U2OS cells with three different cisplatin concentrations. The biphasic decay observed in panel (b) is not observed here. (f) Pre-cisplatin single-cell IMT distribution is better explained with a single EMG (Exponentially Modified Gaussian) model than a mixture of 2 EMGs as verified by comparing AIC values for the two models. We first asked whether distinct sub-populations of slowly dividing cells could be detected in the183 drug naive HCT116 and U2OS cells. The HCT116 dataset had both the daughter cells of dividing184 8 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint cells tracked (Fig. 1a) allowing extraction of single-cell times and cell fates for all the cells in185 the population. To find the population-time curve (Fig. 1b), we summed the number of cells186 in each time frame. Starting with a population of 65 cells, the semi-log plot of the curve shows187 characteristic exponential growth before the addition of the drug (vertical line in Fig. 1b) and a188 bi-phasic exponential decay after cisplatin exposure. The first larger decay rate switches to a plateau189 after the 100th hour of the experiment, driven by the drug tolerant persisters in the population. We190 identified persisters as those cells surviving till the end of the experiment (resulting in the second191 phase of exponential population decay) and extracted the IMTs of their drug-naive ancestors by192 tracking back in time along their lineages (Fig. 1a, blue line). We compared these drug-naive193 ancestor IMTs to those of the sensitive cells that died after cisplatin treatment (Fig. 1a, red line),194 and could not detect any statistically significant difference between these two distributions (Fig. 1c;195 see SI section 3 for details).196 The U2OS dataset was generated by tracking a single randomly chosen daughter after each cell197 division (Fig. 1d; starting population ∼7 cells per field of view). The population-time curves,198 calculated by averaging over the number of cells across 51 fields of view per time frame, showed a199 characteristic exponential growth before cisplatin addition and concentration-dependent exponential200 decay post cisplatin exposure (Fig. 1e). However, five days of tracking was not sufficient to observe201 the second slower exponential decay unlike in the HCT116 cells. Therefore to investigate whether202 multiple sub-populations with distinct IMT distributions exist pre-cisplatin, we resorted to fitting the203 data with a mixture model (see SI section 3 for details). We thus compared the ability of single and204 two-component Exponentially Modified Gaussian (EMG) models to explain the IMT distribution205 of drug naive cells, using the Akaike Information Criterion (AIC). We found that the single EMG206 (AIC = 3687.4) fared better than the mixture model (AIC = 4137.649). Also, the two models207 generated best fit distributions that were essentially indistinguishable, suggesting the absence of208 multiple sub-populations with distinct proliferation rates within the U2OS population as well (Fig.209 1f).210 In summary, our analysis demonstrates that in both HCT116 and U2OS cell types, the heterogeneity211 in cell cycle times before addition of cisplatin is best explained by a single population of cycling cells,212 not by multiple sub-populations with distinctly different cycling kinetics. This suggests that if a pre-213 9 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint DTP state exists amongst these cells before drug administration, it does not have any fitness cost214 as compared to the sensitive cell state.215 Proliferation rate of cancer cells can be predicted from single-cell division216 times.217 Having established that a single IMT distribution best explains the datasets before cisplatin treat-218 ment, we next asked how accurately this IMT distribution can predict the population growth rates.219 For exponentially proliferating cells, the population growth rate is approximately given by the inverse220 of the mean of the IMT distribution. While this is the most popular model used to describe expo-221 nential cellular growth, the underlying assumption here is that the IMT distribution is exponential.222 However, this is clearly not true in our datasets (Fig. 1f) and hence a more general approach is223 required. An age-structured population model (Model M0 ; see Methods and SI section 9 for details)224 is a generalization of the simple exponential growth models, that allows for non-exponential IMT225 distributions, and connects the population growth rate to the full IMT distribution via the Euler-226 Lotka equation (see Methods and SI section 4 for more details)[38, 41, 37, 40]. From the measured227 single-cell IMTs we estimated the predicted population growth rates from both the inverse of the228 mean as well as the Euler-Lotka equation, and compared them with the true growth rate – the linear229 regression slope of the semi log plot of number of cells as a function of time (Table 1).230 For the HCT116 cells, we found that the inverse of the mean gave an error of around 62%, while the231 Euler-Lotka estimate gave a much reduced error of 13% (Table 1 and Fig. 2a). The large error in the232 former model was not surprising as the measured IMT distribution was neither exponential nor very233 narrow to be effectively represented by the mean. For the latter result, we further investigated the234 potential origins of the 13% error. To check whether this error could arise simply due to the small235 number of cells in our experiments, we performed stochastic simulations taking the measured IMT236 distribution as input to generate in silico proliferation trajectories with similar cell numbers.237 10 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint Figure 2: An age-structured model accurately predicts population growth rates from single-cell IMT distributions. (a)-(c) are plots for HCT116, (d)-(f) are for U2OS. (a) Semi-log plot of population versus time curve of HCT116 cells before cisplatin administration. Predicted growth rate using inverse of the mean (pink line; error = 62%) and the Euler-Lotka equation (blue line; error = 13%) are compared with the experimental data shown as black dots. (b) Maximum Likelihood inference of the age distribution of ancestor cells immediately after seeding. (c) Distribution of errors in Euler- Lotka predictions of growth rate, obtained by running 500 iterations of an age-structured population simulation with initial age distribution taken from panel (b). Vertical dashed line represents the Euler-Lotka prediction error in experimental data (corresponding to error of the blue line in panel (a)). (d) Semi-log plot of population versus time curve of U2OS cells before cisplatin, with colours as in panel (a). Predicted growth rate using inverse of the mean (pink line; error = 73%) and the Euler-Lotka equation (blue line; error = 22%) are compared to the experimental data shown in black dots. (e) Inferred age distribution of seeded population of U2OS cells, as in panel (b). (f) Distribution of errors in Euler-Lotka estimates, obtained by running 500 iterations of simulations as in panel (c). Initial age distribution was taken from panel (e) for these simulations. Vertical dashed line represents Euler-Lotka prediction error in experimental data (corresponding to error of the blue line in panel (d)). 11 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint For each simulation run we calculated the growth rate prediction errors using both the growth238 rate estimates, mimicking the procedure we carried out on the experimental data. Since the age239 distribution of the cells at the initial time is an input to the simulations, and it is known that this240 distribution can affect the emergent population dynamics [38], we estimated the age distribution241 of the initial cells directly from the experimental data (Fig. 2b). Based on the IMT distribution242 of the full pre-cisplatin population, we inferred the starting ages of the initially plated cells in the243 experiment using a Maximum Likelihood approach (Fig. 2b; details in SI). We then repeated the244 simulations using the inferred initial age distribution and measured IMT distributions and found245 that the Euler-Lotka errors peaked around 2% (Fig. 2c). We carried out a similar analysis for the246 U2OS dataset, where we once again found a large error of around 73% from the inverse of the mean247 prediction, and a much reduced error of 23% from the Euler-Lotka equation (Fig. 2d). We found248 that the inferred initial age distribution of the starting U2OS cells was broader than that of the249 HCT116 cells (Fig. 2e), and when used in the simulations, could account for about 4% error in the250 growth rate prediction (Fig. 2f). The unexplained errors of about 10 − 20% for HCT116 and U2OS251 cells could be due to artefacts of the image analysis pipeline, or from other unaccounted biological252 sources such as cell-size control of the cell cycle.253 Model for growth rate prediction HCT116 U2OS Inverse of mean IMT 62% 73% Euler-Lotka 13% 23% Expected Errors (Simulations) ∼ 2% ∼ 4% Table 1: An age structured model (Euler-Lotka equation) accurately predicts population growth rates of both HCT116 and U2OS cells before cisplatin treatment. Shown here are deviations (in %) of growth rates predicted by either the Euler-Lotka equation or the inverse of the mean IMT, from the experimentally measured growth rate (obtained via linear regression of the cell-population versus time curve in Fig. 2a,d). The growth rate prediction by the Euler-Lotka equation for HCT116 cells is within ∼ 11% of the error expected from age-structured population simulations with the inferred starting age-distributions. For the U2OS cells, the predictions are off by about ∼ 20%. Our results in this section demonstrate that in the absence of cisplatin, the age-structured model254 predictions of population growth rate using a single IMT distribution is accurate to within 10 − 20%255 in both HCT116 and U2OS cell types. This provides a benchmark for the expected error level when256 trying to predict population dynamics from our single-cell datasets. Additionally, the fact that a257 single IMT distribution is sufficient to predict the population growth rate is consistent with results258 12 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint from the previous section, which suggested the absence of sub-populations of cells with different259 cycling kinetics.260 First exponential decay of drug treated cancer cells is not established261 by single-cell division and death times262 Since the age-structured model (Model M0 ) could predict pre-cisplatin population growth rates from263 the measured cell division times within 20% error, we next asked if a similar prediction accuracy264 can be achieved in the post-cisplatin scenario, where an added dimension is cell death causing the265 populations to decay instead of grow. Throughout this section, we will be focusing only on the266 first phase of exponential decay after cisplatin treatment (for the U2OS cells the experiments were267 not long enough to observe the second exponential phase). To find the IMT and AT distributions268 for both HCT116 and U2OS cells, we extracted the times from cells that were either born before269 drug exposure but divided or died after drug exposure, or cells that were both born and reached a270 fate after drug exposure (Fig. 3a-d). For the HCT116 cells there was data from only one cisplatin271 concentration (Fig. 3a, b), but three different concentrations for U2OS cells allowed us to explore272 concentration effects on the population dynamics (Fig. 3c, d).273 Though the U2OS population decay rates were distinctly different between high (−25x10 −3h−1)274 , medium (−11.4x10 −3h−1) and low ( −7.2x10−3h−1) cisplatin, surprisingly there were negligible275 differences in the IMT and AT distributions across the three cisplatin concentrations (Figs. 3c and276 3d respectively). This directly showed, in a model-independent manner, that times to division or277 death post cisplatin (and consequently the division and death rates) do not determine the population278 decay rates. This result contradicts the widely used assumption that the effect of drugs can be well279 described either as increasing the cellular death rate or reducing the division rate [30, 31, 32, 33,280 34, 35]. To rigorously confirm that the minor differences in the measured IMT and AT distributions281 cannot explain the large differences in decay rates across cisplatin concentrations, we calculated the282 predicted rates using both inverse of the mean and Model M0 (see SI section 9 for details). The283 IMT and AT distributions were not sufficient to accurately predict the concentration-dependent284 decay rates of the U2OS cells, producing errors greater than 100% (‘Measured’ columns in Table285 2). Indeed, even for HCT116 cells where we had data from one drug concentration, the error in286 13 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint predicted287 Figure 3: Inferring the underlying post-drug IMT and AT distributions from the experimentally measured, biased distributions. (a)-(d) are plots of measured IMT and AT distributions of HCT116 and U2OS cells while (e)-(h) are the competing-risks corrected distributions inferred using a Markov Chain Monte Carlo approach (see Methods and SI section 8 for details). (a)-(b) IMT and AT distri- butions of HCT116 cells. All cells that divided or died during the treatment period were tracked and measured, including cells that were born before and after cisplatin addition. (c)-(d) Forward lineage distributions for U2OS cells that have division and death events that occurred during the treatment period. The three distributions in both panels (c) and (d) are almost identical, demonstrating that increasing drug concentration does not affect division or death rates post cisplatin. (e) Inferred IMT distribution of HCT116 cells. Compared to panel (a), the long tail is evident demonstrating the significant bias in the measured distribution. (f) Inferred AT distribution of HCT116 cells, which is similar to the measured distribution in panel (b). (g) Inferred IMT distributions for U2OS cells reveal increasingly long tails for increasing cisplatin concentrations. (h) Similar to panel (f), the inferred AT distributions for U2OS cells is almost identical to the measured AT distributions. population decay rates was greater than 100% (‘Measured’ column in Table 2). Taken together, these288 findings were remarkably surprising in light of our previous section results as well as the vast body289 of literature modeling population growth/decay as the difference in birth and death rates.290 We therefore wondered if our analysis using the measured IMT and AT distributions suffered from a291 competing-risks bias, a statistical effect that prevents direct measurement of the correct underlying292 14 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint distributions when the events are mutually exclusive. We have previously demonstrated that this293 effect leads to strong biases in the IMT and AT distributions measured after drug treatment [18],294 and hence can potentially lead to erroneous population growth rate predictions. A simple analytic295 calculation for exponentially distributed IMT and AT is presented in SI section 7 that provides296 intuition for this competing risks effect, and how it can be corrected using inference techniques.297 In the more complex case of non-exponentially distributed IMT and AT as is the case with the298 HCT116 and U2OS datasets, we used a Markov Chain Monte Carlo (MCMC) approach to infer299 the parameters of the unobserved distributions (see SI section 8 for details). The inferred IMT300 distributions of HCT116 cells (Fig. 3e) showed a significant right skew (slower division times)301 compared to the measured IMT distribution, which was not observed for the AT distribution (Fig.302 3f). For the U2OS cells, the degree of right skewness in the inferred IMT distributions increased with303 higher drug concentrations (Fig. 3g) while once again the inferred AT distributions (Fig. 3h) had304 no significant skew compared to the measured AT. Using these inferred IMT and AT distributions305 we then recalculated the predicted decay rates using the mean inverse and the Eula-Lotka equation.306 While the errors significantly reduced (‘Inferred’ columns in Table 2), they were still very high with307 the smallest error around 40% in both HCT116 and U2OS and the largest error greater than 100%308 for low and medium cisplatin doses in U2OS.309 HCT116 U2OS Measured Inferred Measured Inferred Model Low Med High Low Med High Mean Inverse > 100% 51.8% > 100% > 100% > 100% > 100% 85% 60.2% Euler-Lotka > 100% 40.3% > 100% > 100% > 100% > 100% > 100% 40.3% Table 2: Age structured models fail to predict population decay rates for cells treated with cis- platin. Shown here are deviations (in %) of the predicted decay rates from the measured decay rates (obtained via linear regression of post-cisplatin experimental data in Fig. 1b,e). Predictions were made using IMT and AT distributions obtained directly from the experiment (‘Measured’) and from competing risks-corrected unbiased distributions (‘Inferred’). In summary, we found that estimates of population decay rate using measured single-cell times to310 division and death post cisplatin treatment were erroneous by over 100%, predicting positive growth311 instead of decay. Even after correcting for competing risk biases, the population decay rates could312 be predicted at best with an error of 40%, often with errors as high as 100%. These counter-intuitive313 15 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint results suggest that measured single-cell IMT and AT are not sufficient to estimate population decay314 rates post addition of cisplatin. In turn these results suggest that post-drug stochastic competition315 of cell fates inherent to these exponential growth models (Model M0 ), is likely to be an incorrect316 model of cell-fate decisions in the development of cancer drug persistence.317 Cell-state inheritance along with fate decisions before and during drug318 treatment control cancer population dynamics.319 We next asked if the nature and timing of fate decisions can be determined by developing a series320 of models with varying assumptions and checking for compatibility with experimental data. The321 deviations from Model M0 predictions of decay rates suggest that the cell fates are likely decided322 either during or before drug addition, largely independent of the single-cell times. Therefore we323 developed three alternate models of the fate decision process, Models M1, M2 and M3 (Fig. 4a).324 In these models, cell fates are either (i) determined based on 2 cell-states before drug – a pre-DTP325 state P’ that is more likely to survive and generate persisters, and a sensitive state S’ that is more326 likely to lead to cell death after drug administration (M1 ), (ii) decided randomly at the time of drug327 addition with probabilities that depend on the drug concentration (M2 ), or (iii) a combination of328 both M1 and M2 where pre-existing cell-states determine the end fate with large probability while329 the drug concentration modulates this probability to some extent ( M3 ). Inter-conversion between330 states S’ and P’ in Models M1 and M3 is modeled using a Markov Chain described by a transition331 matrix T while the drug concentration-dependent fate probabilities in Model M3 are imposed using332 a fate matrix F. Further details and parameterizations of each model are provided in Methods and333 SI Sections 9 and 10.334 Interestingly, we found that though Models M1 and M2 could recapitulate certain aspects of the335 experimental datasets (see SI section 9 for details), only Model M3 could quantitatively explain all336 three major observations simultaneously: (1) decay rate for a fixed drug concentration,337 16 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint Figure 4: A two state Markov model with state-dependent fates (Model M3 ) quantitatively ex- plains the experimentally measured decay rate and correlations in lineage related HCT116 cells. (a) Schematics of models M0 to M3. Cells in gray indicate no role of cell-states in determination of fate. Coloured cells indicate pre-existing cell-states that pre-determine final fates to a large extent. The dice indicates probabilistic fate assignment at the time of drug administration. (b) Lineage corre- lations in fates of HCT116 cells measured from time lapse experiments. Correlations exist for first and second cousins (FC and SC), while third cousins (TC) are similar to randomly selected pairs of non-lineage related cells (horizontal dashed line). (c) An example run of Model M3 showing the rel- ative populations of the two states as a function of time before drug administration. (d) Population versus time trajectories from simulations of Model M3 using identical IMT and AT distributions, a fixed Transition matrix T, but different Fate matrices. Each colour corresponds to simulations with a different F. (e) Same as panel (d), but with different Transition matrices for a fixed Fate matrix F. (f) Distribution of decay rates from the simulations (red histogram). The vertical dashed line is the decay rate obtained from the HCT116 experimental dataset. (g) Comparison of lineage correlations obtained from the HCT116 experiment (red) with simulations of Model M3 (blue). Transition rates used in the simulation were symmetric for transition between the two states, 9 × 10−3 per hour. The fate probabilities for cells to die, divide and stay alive for the two states were: Sensitive state S’ = (0.9125, 0.05, 0.0375); pre-DTP state P’ = (0.35, 0.1, 0.55). 17 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint (2) lineage correlations in end-fate and (3) drug concentration dependence of the population decay338 rates. On simulating M3, we noticed that a broad spectrum of effective decay rates could be obtained339 after cisplatin exposure even with the IMT and AT distributions fixed. This arises either due340 to different Fate matrices for a given Transition matrix (Fig. 4d) or due to different Transition341 matrices for a given Fate matrix (Fig. 4e). This explains the surprising results from Fig. 3c-342 d where the IMT and AT distributions were found to be insensitive to drug concentration – the343 IMT and AT distributions have little contribution towards establishing the decay rate. Rather,344 the population dynamics is dominated by the relative proportion of the two states at time of drug345 exposure (determined by T), and the entries of F. Importantly, we found that for certain physically346 realistic values of the T and F matrices, we could approximately recapitulate the experimentally347 measured decay rates (Fig. 4f) as well as the lineage correlations (Fig. 4g) while maintaining the348 same IMT and AT distributions as measured in the HCT116 experiment.349 Our analysis in this section quantitatively explains why the birth and death rates (and more generally350 the single-cell IMT, AT distributions) do not determine the decay rates of drug treated cells at351 the population level. We demonstrate that the transition rates between cell states before drug352 treatment, the ability of these states to be passed down across cellular generations, and the cisplatin353 concentration-dependent propensity of these states to result in death or survival largely determine354 the population dynamics and lineage correlations of cells after drug treatment.355 Change in barcode diversity before and after drug cannot establish timing356 of fate decisions in the presence of cell-state switching.357 Our analyses thus far demonstrated the importance of pre-existing cellular states in the population358 that leads to persistence of cells in the face of cisplatin treatment. This conclusion was true for both359 colorectal cancer as well as osteosarcoma cell types. However, a recent study on colorectal cancer360 cells concluded that the persister fate arose after the addition of cisplatin, based upon analysis of361 a high-complexity barcode library that had been transduced into the cancer cells [24]. In short,362 this study found that the barcode diversity (measured by the Shannon Diversity Index) showed no363 reduction upon cisplatin treatment, which suggested that the persister fate choice was likely made364 post drug treatment. Since these results are conceptually different from our findings,365 18 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint Figure 5: Maintenance of barcode diversity before and after drug cannot be interpreted as a signa- ture of drug-induced cell fates. (a) Current barcoding methods relate fate determination to pre-drug cell-states if the barcode diversity, shown here by the presence of distinct barcodes (colored circles), significantly reduces after drug (i). However, if the diversity does not reduce then fate determina- tion is assumed to be in response to drug exposure (ii). Our model however, suggests that both the conditions (i & ii) can be obtained with fates largely determined before drug administration. Simulation results show that in conditions when the cell-states before drug transition on time scales much slower than the cell division time, barcode diversity shrinks significantly post drug treatment (iii). However, if the transition rates are of the order of a few cell cycles or faster, more sensitive cells transition to a pre-DTP state, resulting in insignificant change in barcode diversity post-drug (iv). (b) Barcode abundance plot post drug treatment obtained from simulations of Model M3. Circle fill colors represent distinct barcodes while their sizes show barcode frequency counts. Results of 3 simulations are compared. ‘Control’: simulations with cell-state independent fate assignment of cells; ‘No Trans’: simulations with cell-fate assignment conditional on cell-states that do not inter-convert; ‘SimExpt’: simulations with state dependent fate assignment with cell-state transitions allowed be- fore drug addition. (c) Diversity of barcodes can be measured by the Shannon Diversity Index. The boxplot here shows change in Shannon Diversity Index for the corresponding simulations described in panel (b). 19 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint Figure 5: (d) Variance over Mean Ratio (VMR) of persister cell counts as measured in the experiment (Expt) and the simulation matching experimental results (SimExpt). The horizontal dashed line shows the VMR for a Poisson distribution. (e) Comparison of Cumulative Distribution Function (CDF) for perister cell counts obtained at the end of simulation (green) and experiment (red). ‘Poisson fit’ (blue line) is the CDF expected from a Poisson distribution of persister cell counts. we next investigated potential reasons that may explain the apparently contradictory results.366 The barcoding approach to determine the timing of persister fate choice relies on comparing the367 barcode diversity, or abundance of unique barcodes, before and after drug exposure in many cellular368 replicates. A reduction in barcode diversity after drug treatment is interpreted as selection of a small369 set of pre-existing persister cells (Fig. 5a i). Maintenance of barcode diversity on the other hand370 suggests a scenario where persister cells arise randomly in cells after drug treatment, thereby leading371 to many surviving barcodes (Fig. 5a ii). However, these interpretations were originally derived in372 the context of genetic mutations [42] where reversal of mutations is rare. The same assumptions do373 not hold in case of non-genetic cell states (the pre-DTP state P ′ or sensitive state S′ in our case),374 and we argued in the previous sections that transitions between the two states is a fundamental375 aspect of the dynamics in both HCT116 and U2OS cell types.376 To quantitatively explore how cell-state switching before drug addition could affect interpretations377 of barcode diversity, we used Model M3 with the same parameters from the previous section, with378 the additional information of barcodes (Fig. 5a iii-iv; see SI sections 11-12 for details). In brief,379 each starting cell was given a barcode ID that was inherited by all descendants across cell divisions380 and maintained over cell-state changes. The fates were conditioned on cell-state as described before381 for Model M3. The results of the simulation are shown as an abundance bubble plot (Fig. 5b)382 for qualitative visualisation and as a Shannon Diversity Index (SDI) plot (Fig. 5c) for quantitative383 assessment of barcode enrichment or diversity. The drop in Shannon Diversity Index post drug was384 observed only when the simulation did not allow for state transitions before drug treatment, which385 is equivalent to very slow transitions compared to a typical cell cycle time (Fig. 5b-c, ‘No Trans’).386 However, if the transitions were allowed with rates that gave rise to the experimentally measured387 lineage correlations (Fig. 4g), it resulted in no change in the SDI value (Fig. 5b-c, ‘SimExpt’). This388 short analysis confirmed that the timing of development of pre-DTP states cannot be inferred using389 current barcode analysis techniques and might require development of novel methods to account for390 20 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint cell-state transitions.391 Finally, we used an orthogonal approach to analyse the timing of fate decisions for the HCT116392 cells – a Luria-Delbruck (LD) fluctuation test. This method is based on quantifying the persister393 count distribution after cisplatin, with a Poisson distribution suggesting drug-induced persistence,394 while a super-Poisson distribution implying pre-existence of a pre-DTP state before drug addition395 [43]. The comparison is made either by calculating the Index of Dispersion (variance over mean) or396 obtaining the cumulative distribution function of the persister cell counts. We calculated the Index397 of Dispersion for HCT116 data by comparing the persister cell counts across lineage trees and found398 it to be around 4.5, which meant that the distribution of persisters was not Poissonian and that399 the persister fate originated before drug addition (see SI section 13 for details). Note that for the400 U2OS cells this analysis was not possible since the second exponential decay phase was not observed,401 precluding identification of persisters (Fig. 1e). Finally, we also checked if our simulation results402 that matched the decay rates (within 5% error; Fig. 4f) and lineage correlations (Fig. 4g) seen in the403 experiment, would give similar results. Indeed, the simulation results also showed non-Poissonian404 distributions (Fig. 5d). By comparing the CDF of the persister counts in simulation and experiment,405 we further confirmed the distribution to be non-Poissonian (Fig. 5e).406 In summary, the results in this section show that change in barcode diversity before and after drug407 treatment is not always sufficient to establish the timing of fate decisions. This is true when fates408 are related to non-genetic cell-states which inter-convert between each other on time scales similar409 to the cell-cycle time. Maintenance of barcode diversity after drug treatment is therefore consistent410 with a pre-DTP state existing before drug treatment and driving fate outcomes, as our analyses411 demonstrates using both lineage correlations as well as the Luria-Delbruck fluctuation test.412 Discussion413 In this work we demonstrate how general principles underlying the dynamics of drug persistence can414 be discovered in the process of connecting measured single-cell kinetics to the emergent population415 dynamics of cancer cells. Quite unexpectedly, the single-cell measurements directly show that in-416 creasing drug concentration does not affect the birth or death rates (Fig. 3c-d), unlike the widespread417 21 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint assumption in the field [30, 32, 33, 34, 35]. These results suggest that stochastic competition be-418 tween fates based on drug-induced cellular birth and death rates is not the mode of decision-making419 underlying the emergence of non-genetic drug tolerance. Rather, we demonstrate that non-genetic420 cell-states existing before drug treatment largely bias the eventual fates towards drug-tolerance or421 susceptibility. These biases are further modulated at the time of addition of the drug, leading to a422 quantitative explanation of both the population dynamics as well as the lineage correlations observed423 in experimental datasets of HCT116 and U2OS cells (Fig. 4f-g). While earlier studies have discussed424 the role of pre-existent cell-state switching in cancer drug tolerance, they typically use markers of425 specific cell-states to demonstrate state-switching and emergence of drug tolerant persisters [3, 44,426 45, 46, 47]. Our work however is based on more general arguments for the existence of non-genetic427 states and transitions between them, requiring no a priori knowledge of the molecular architecture of428 the underlying states. This is particularly powerful since the molecular details of these states remain429 unknown in all generality, and are likely to comprise transcriptional, post-transcriptional as well as430 epigenetic elements. Additionally, to the best of our knowledge this work is the first to quantitatively431 explain the connection between measured single-cell kinetics and the emergent population dynamics432 of drug-treated cancer cells.433 In the widely used exponential growth models and their generalizations, the fraction of cells that die434 (establishing the first decay rate) depend entirely upon the parameters of the post-drug IMT and435 AT distributions. However, in the presence of non-genetic cellular states that pre-determine fate436 outcomes (sensitive and pre-DTP states), the fraction of cells that die no longer depends exclusively437 on the IMT and AT distributions, but also on the transition rates of these pre-existing cell states.438 Though this idea has been suggested in many previous theoretical works that typically analyze439 FACS datasets [3, 48, 49, 50, 51, 35, 52], it has not been directly demonstrated from single-cell440 measurements linking IMT and AT distributions to the population dynamics. Additionally, the idea441 that single-cell division and death kinetics may not be sufficient to predict drug-treated population442 dynamics has been suggested before in the context of lung cancer cells [53]. However, this study did443 not incorporate distributions of time to death or correct for competing-risks biases, thereby making444 it challenging to interpret the results.445 Furthermore, the implications of cell-fate lineage correlations have been poorly appreciated in the446 22 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint context of non-genetic heterogeneity and the timing of emergence of pre-DTP states. While early447 studies on TRAIL-induced apoptosis had already suggested that lineage correlations imply early fate448 decisions[2, 17], a careful study later explored the detailed implications of correlations in apoptosis449 and their decay across sisters, cousins and higher order lineage-related cells [16]. Based on careful450 time-lapse microscopy experiments, this study separated out the effects of spatial correlations po-451 tentially arising from inter-cellular signalling from correlations arising from inheritance of cell-states,452 to argue for decisions of apoptosis taken well before addition of the drug [16]. Recently we demon-453 strated that correlations also exist in the dynamics of signaling pathways that get upregulated only454 after drug treatment (rate of p53 increase in response to cisplatin), suggesting that even if cellu-455 lar states change upon drug addition, the propensity for these changes are largely pre-programmed456 into upstream regulatory networks of the ancestor cells existing before drug addition [18]. Taken457 together with our current study, the emerging evidence points towards a scenario where cell-states458 existing well before drug addition largely determine eventual cell fates after addition of drugs. In-459 terestingly, recent tour-de-force studies using clever combinations of high-complexity barcoding and460 single-molecule FISH have arrived at broadly similar conclusions. These newer studies have also461 begun to shed light on important and previously unavailable information on the genes whose expres-462 sion states pre-determine the ultimate fate of cancer cells in response to drugs[19, 21], and further463 demonstrate that transcriptional states post long-term drug treatment may be very different from464 those that existed prior to drug exposure [28, 54].465 An intriguing implication of our results is that both the pre-DTP as well as the drug-sensitive states466 exist in cells that are actively proliferating. This suggests that DTPs need not arise exclusively467 from quiescent (also called dormant) or slowly cycling cells (such as cancer stem cells), which is468 primarily the evidence that has been reported in the literature from bacteria, mammalian cell lines,469 primary cells as well as tissues [12, 55, 56, 57, 26, 58, 59]. Importantly, our inability to find multiple470 proliferative sub-populations in either HCT116 or U2OS cells before drug treatment (Fig. 1c,f),471 combined with the fact that we could accurately predict the population growth rates from the472 measured single IMT distribution (Fig. 2c,f), demonstrates that the pre-DTP state does not have473 any fitness cost. This result is intriguing since a fitness cost is typically associated with drug resistance474 arising from genetic mutations [34, 60, 61], and has recently been suggested for non-genetic persister475 23 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint states as well [27]. It is interesting to note however that though the IMT distributions do not476 distinguish the ancestors of persisters versus sensitive cells, the number of divisions undergone can477 be predictive of eventual cell fate after drug treatment [36].478 The use of single-cell time lapse microscopy datasets in this study provides an advantage over other479 barcoding studies in terms of the high resolution at which lineages can be reconstructed. Indeed,480 these results demonstrate the advantage of using single-cell resolved lineage correlations over barcode-481 diversity based methods to infer timing of fate decisions (Fig. 5b-c). However, these experiments482 and hence our accompanying analyses have the major limitation of exploring very short time-scales483 after drug treatment (on the order of a few days). As a result, we cannot comment on effects of484 long-term drug exposure on cell-state switching, or drug induced gradual cell-state reprogramming485 over long time scales [5, 6, 20, 7, 28]. Careful combinations of time-lapse microscopy along with486 barcoding studies will be required to delineate short-term versus long term dynamics and provide a487 full picture of non-genetic drug tolerance that eventually leads to genetic resistance.488 Conclusions489 We quantitatively demonstrate in this study that stochastic choices between survival and death based490 on post-drug division and death rates is not how cancer cells make fate decisions. Rather, the choices491 are largely pre-determined via cell-states that exist at least 2-3 generations before drug treatment,492 and are inherited by successive generations. This is true both for the cell-state that likely leads to493 drug-tolerant persisters as well as the state that makes cells susceptible, thereby suggesting that the494 notion of a fitness cost need not necessarily apply to non-genetic mechanisms driving persistence. To495 a smaller extent the fate decisions are also modulated by the drug, and together with the pre-existing496 cell-states, quantitatively explain the emergence of lineage correlations in end-fate and the population497 dynamics. We argue that the existence of lineage correlations is a better alternative to inferring the498 timing of fate decisions as opposed to barcode diversity analysis, which is hard to interpret in the499 presence of cell-state switching. These arguments do not require knowledge of the specific molecular500 details of the underlying cell-states, and hence provide a general and powerful approach to the study501 of drug-tolerance on short time scales driven by non-genetic heterogeneity.502 24 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint References503 [1] Sreenath V. Sharma et al. “A chromatin-mediated reversible drug-tolerant state in cancer cell504 subpopulations.” In: Cell 141.1 (Apr. 2, 2010). Place: United States, pp. 69–80.505 [2] Sabrina L. Spencer et al. “Non-genetic origins of cell-to-cell variability in TRAIL-induced506 apoptosis.” In: Nature 459.7245 (May 21, 2009). Place: England, pp. 428–432.507 [3] Piyush B. Gupta et al. “Stochastic State Transitions Give Rise to Phenotypic Equilibrium in508 Populations of Cancer Cells”. In: Cell 146.4 (Aug. 2011). Publisher: Elsevier, pp. 633–644.509 [4] Michael Ramirez et al. “Diverse drug-resistance mechanisms can emerge from drug-tolerant510 cancer persister cells”. In: Nature Communications 7.1 (Feb. 2016). Publisher: Nature Pub-511 lishing Group, p. 10690.512 [5] Sydney M. Shaffer et al. “Rare cell variability and drug-induced reprogramming as a mode of513 cancer drug resistance”. In: Nature 546.7658 (June 2017). Publisher: Nature Publishing Group,514 pp. 431–435.515 [6] Robert Vander Velde et al. “Resistance to targeted therapies as a multifactorial, gradual adap-516 tation to inhibitor specific selective pressures”. In: Nature Communications 11.1 (May 2020).517 Publisher: Nature Publishing Group, p. 2393.518 [7] Mariangela Russo et al. “A modified fluctuation-test framework characterizes the population519 dynamics and mutation rate of colorectal cancer persister cells”. In:Nature Genetics 54.7 (July520 2022). Number: 7 Publisher: Nature Publishing Group, pp. 976–984.521 [8] Yi Pu et al. “Drug-tolerant persister cells in cancer: the cutting edges and future directions”.522 In: Nature Reviews Clinical Oncology 20.11 (Nov. 1, 2023), pp. 799–813.523 [9] JosephW. Bigger. “TREATMENT OF STAPHYLOCOCCAL INFECTIONS WITH PENI-524 CILLIN BY INTERMITTENT STERILISATION”. In: The Lancet 244.6320 (Oct. 14, 1944).525 Publisher: Elsevier, pp. 497–500.526 [10] H. E. Skipper, F. M. Schabel, and W. S. Wilcox. “EXPERIMENTAL EVALUATION OF PO-527 TENTIAL ANTICANCER AGENTS. XIII. ON THE CRITERIA AND KINETICS ASSOCI-528 ATED WITH ”CURABILITY” OF EXPERIMENTAL LEUKEMIA”. In: Cancer Chemother-529 apy Reports 35 (Feb. 1964), pp. 1–111.530 25 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint [11] Karl Kochanowski et al. “Drug persistence – From antibiotics to cancer therapies”. In: Current531 Opinion in Systems Biology. Pharmacology and drug discovery 10 (Aug. 2018), pp. 1–8.532 [12] Nathalie Q. Balaban et al. “Bacterial Persistence as a Phenotypic Switch”. In: Science 305.5690533 (Sept. 10, 2004). Publisher: American Association for the Advancement of Science, pp. 1622–534 1625.535 [13] Nathalie Q. Balaban et al. “A problem of persistence: still more questions than answers?” In:536 Nature Reviews Microbiology 11.8 (Aug. 1, 2013), pp. 587–591.537 [14] Mariangela Russo, Alberto Sogari, and Alberto Bardelli. “Adaptive Evolution: How Bacteria538 and Cancer Cells Survive Stressful Conditions and Drug Treatment”. In: Cancer Discovery539 11.8 (Aug. 2021), pp. 1886–1895.540 [15] Charles C. Bell and Omer Gilan. “Principles and mechanisms of non-genetic resistance in541 cancer”. In: British Journal of Cancer 122.4 (Feb. 1, 2020), pp. 465–472.542 [16] Patrick D. Bhola and Sanford M. Simon. “Determinism and divergence of apoptosis suscepti-543 bility in mammalian cells”. In: Journal of Cell Science 122.23 (Dec. 1, 2009), pp. 4296–4302.544 [17] Fran¸ cois Bertaux et al. “Modeling Dynamics of Cell-to-Cell Variability in TRAIL-Induced545 Apoptosis Explains Fractional Killing and Predicts Reversible Resistance”. In: PLOS Compu-546 tational Biology 10.10 (Oct. 2014). Publisher: Public Library of Science, pp. 1–13.547 [18] Shaon Chakrabarti et al. “Hidden heterogeneity and circadian-controlled cell fate inferred from548 single cell lineages”. In: Nature Communications 9.1 (Dec. 18, 2018), p. 5372.549 [19] Sydney M. Shaffer et al. “Memory Sequencing Reveals Heritable Single-Cell Gene Expression550 Programs Associated with Distinct Cellular Behaviors.” In: Cell 182.4 (Aug. 20, 2020). Place:551 United States, 947–959.e17.552 [20] Yaara Oren et al. “Cycling cancer persister cells arise from lineages with distinct programs”.553 In: Nature 596.7873 (Aug. 1, 2021), pp. 576–582.554 [21] Benjamin L. Emert et al. “Variability within rare cell states enables multiple paths toward555 drug resistance”. In: Nature Biotechnology 39.7 (July 1, 2021), pp. 865–876.556 [22] Lea Schuh et al. “Gene Networks with Transcriptional Bursting Recapitulate Rare Transient557 Coordinated High Expression States in Cancer”. In: Cell Systems 10.4 (Apr. 2020). Publisher:558 Elsevier, 363–378.e12.559 26 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint [23] Angela Oliveira Pisco et al. “Non-Darwinian dynamics in therapy-induced cancer drug resis-560 tance”. In: Nature Communications 4.1 (Sept. 18, 2013), p. 2467.561 [24] Sumaiyah K. Rehman et al. “Colorectal Cancer Cells Enter a Diapause-like DTP State to562 Survive Chemotherapy.” In: Cell 184.1 (Jan. 7, 2021). Place: United States, 226–242.e21.563 [25] Robert J. Vanner et al. “Quiescent sox2(+) cells drive hierarchical growth and relapse in sonic564 hedgehog subgroup medulloblastoma.” In: Cancer cell 26.1 (July 14, 2014). Place: United565 States, pp. 33–47.566 [26] Fran¸ cois M. Vallette et al. “Dormant, quiescent, tolerant and persister cells: Four synonyms567 for the same target in cancer.” In: Biochemical pharmacology 162 (Apr. 2019). Place: England,568 pp. 169–176.569 [27] David A. Kessler and Herbert Levine. “Phenomenological Approach to Cancer Cell Per-570 sistence”. In: Phys. Rev. Lett. 129.10 (Aug. 2022). Publisher: American Physical Society,571 p. 108101.572 [28] Yogesh Goyal et al. “Diverse clonal fates emerge upon drug treatment of homogeneous cancer573 cells”. In: Nature 620.7974 (Aug. 2023). Publisher: Nature Publishing Group, pp. 651–659.574 [29] Corey E. Hayford et al. “A heterogeneous drug tolerant persister state in BRAF-mutant575 melanoma is characterized by ion channel dysregulation and susceptibility to ferroptosis”.576 In: bioRxiv (Jan. 1, 2022), p. 2022.02.03.479045.577 [30] Jasmine Foo and Franziska Michor. “Evolution of resistance to anti-cancer therapy during578 general dosing schedules”. In: Journal of Theoretical Biology 263.2 (Mar. 21, 2010), pp. 179–579 188.580 [31] Juliann Chmielecki et al. “Optimization of dosing for EGFR-mutant non-small cell lung cancer581 with evolutionary cancer modeling”. In: Science Translational Medicine 3.90 (July 6, 2011),582 90ra59.583 [32] Ivana Bozic et al. “Evolutionary dynamics of cancer in response to targeted combination ther-584 apy”. In: eLife 2 (June 25, 2013), e00747.585 [33] Marc Hafner et al. “Growth rate inhibition metrics correct for confounders in measuring sen-586 sitivity to cancer drugs”. In: Nature Methods 13.6 (June 1, 2016), pp. 521–527.587 27 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint [34] S Chakrabarti and F. Michor. “Pharmacokinetics and drug-interactions determine optimum588 combination strategies in computational models of cancer evolution”. In: Cancer Research in589 press (2017).590 [35] Tyler Cassidy et al. “The role of memory in non-genetic inheritance and its impact on cancer591 treatment resistance”. In: PLOS Computational Biology 17.8 (Aug. 2021). Publisher: Public592 Library of Science, e1009348.593 [36] Adri´ an E. Granada et al.The effects of proliferation status and cell cycle phase on the responses594 of single cells to chemotherapy. ISSN: 1939-4586 1059-1524 Issue: 8 Journal Abbreviation: Mol595 Biol Cell Pages: 845-857 Publication Title: Molecular biology of the cell. Apr. 1, 2020.596 [37] E. O. Powell. “Growth Rate and Generation Time of Bacteria, with Special Reference to Con-597 tinuous Culture”. In: Microbiology 15.3 (1956). Publisher: Microbiology Society Type: Journal598 Article, pp. 492–511.599 [38] Evgeny B. Stukalin et al. “Age-dependent stochastic models for understanding population600 fluctuations in continuously cultured cells.” In: Journal of the Royal Society, Interface 10.85601 (Aug. 6, 2013). Place: England, p. 20130325.602 [39] Farshid Jafarpour et al. “Bridging the Timescales of Single-Cell and Population Dynamics”.603 In: Phys. Rev. X 8.2 (Apr. 2018). Publisher: American Physical Society, p. 021007.604 [40] Jie Lin and Ariel Amir. “The Effects of Stochasticity at the Single-Cell Level and Cell Size605 Control on the Population Growth.” In: Cell systems 5.4 (Oct. 25, 2017). Place: United States,606 358–367.e4.607 [41] E. O. POWELL. “SOME FEATURES OF THE GENERATION TIMES OF INDIVIDUAL608 BACTERIA”. In: Biometrika 42.1 (June 1, 1955), pp. 16–44.609 [42] Hyo-eun C. Bhang et al. “Studying clonal dynamics in response to cancer therapy using high-610 complexity barcoding”. In: Nature Medicine 21.5 (May 2015). Publisher: Nature Publishing611 Group, pp. 440–448.612 [43] Chewei Anderson Chang et al. “Ontogeny and Vulnerabilities of Drug-Tolerant Persisters in613 HER2+ Breast Cancer”. In: Cancer Discovery 12.4 (Apr. 2022), pp. 1022–1045.614 [44] Adrian Biddle et al. “Phenotypic Plasticity Determines Cancer Stem Cell Therapeutic Resis-615 tance in Oral Squamous Cell Carcinoma”. In: eBioMedicine 4 (Feb. 2016). Publisher: Elsevier,616 pp. 138–145.617 28 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint [45] Sugandha Bhatia et al. “Interrogation of Phenotypic Plasticity between Epithelial and Mes-618 enchymal States in Breast Cancer”. In: Journal of Clinical Medicine 8.6 (June 2019). Number:619 6 Publisher: Multidisciplinary Digital Publishing Institute, p. 893.620 [46] Piyush B. Gupta et al. “Phenotypic Plasticity: Driver of Cancer Initiation, Progression, and621 Therapy Resistance”. In: Cell Stem Cell 24.1 (Jan. 2019), pp. 65–78.622 [47] Andrea Sacchetti et al. “Phenotypic plasticity underlies local invasion and distant metastasis623 in colon cancer”. In: eLife 10 (May 2021). Ed. by Margaret C Frame, Anna Akhmanova, and624 Luke Boulter. Publisher: eLife Sciences Publications, Ltd, e61461.625 [48] Stefano Zapperi and Caterina A. M. La Porta. “Do cancer cells undergo phenotypic switching?626 The case for imperfect cancer stem cell markers”. In: Scientific Reports 2 (June 2012), p. 441.627 [49] Joseph Xu Zhou et al. “Nonequilibrium Population Dynamics of Phenotype Conversion of628 Cancer Cells”. In: PLOS ONE 9.12 (Dec. 2014). Publisher: Public Library of Science, pp. 1–629 19.630 [50] Niraj Kumar et al. “Stochastic modeling of phenotypic switching and chemoresistance in cancer631 cell populations”. In: Scientific Reports 9.1 (July 2019). Publisher: Nature Publishing Group,632 p. 10845.633 [51] Einar Bjarki Gunnarsson et al. “Understanding the role of phenotypic switching in cancer drug634 resistance”. In: Journal of Theoretical Biology 490 (Apr. 2020), p. 110162.635 [52] Paras Jain et al. “Cell-state transitions and density-dependent interactions together explain636 the dynamics of spontaneous epithelial-mesenchymal heterogeneity”. In: iScience 27.7 (July637 2024). Publisher: Elsevier.638 [53] Darren R. Tyson et al. “Fractional proliferation: a method to deconvolve cell population dynam-639 ics from single-cell data.” In: Nature methods 9.9 (Sept. 2012). Place: United States, pp. 923–640 928.641 [54] Gustavo S. Fran¸ ca et al. “Cellular adaptation to cancer therapy along a resistance continuum”.642 In: Nature (July 2024). Publisher: Nature Publishing Group, pp. 1–8.643 [55] Nathan Moore, JeanMarie Houghton, and Stephen Lyle. “Slow-Cycling Therapy-Resistant644 Cancer Cells”. In: Stem Cells and Development 21.10 (July 1, 2012), pp. 1822–1830.645 29 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint [56] Antonio Ahn, Aniruddha Chatterjee, and Michael R. Eccles. “The Slow Cycling Phenotype: A646 Growing Problem for Treatment Resistance in Melanoma”. In: Molecular Cancer Therapeutics647 16.6 (June 1, 2017), pp. 1002–1009.648 [57] M. Perego et al. “A slow-cycling subpopulation of melanoma cells with highly invasive proper-649 ties”. In: Oncogene 37.3 (Jan. 2018). Number: 3 Publisher: Nature Publishing Group, pp. 302–650 312.651 [58] Tri Giang Phan and Peter I. Croucher. “The dormant cancer cell life cycle”. In: Nature Reviews652 Cancer 20.7 (July 2020). Number: 7 Publisher: Nature Publishing Group, pp. 398–411.653 [59] Akihisa Seita et al. “Intrinsic growth heterogeneity of mouse leukemia cells underlies differential654 susceptibility to a growth-inhibiting anticancer drug”. In: PLOS ONE 16.2 (Feb. 1, 2021).655 Publisher: Public Library of Science, e0236534.656 [60] Artem Kaznatcheev et al. “Fibroblasts and alectinib switch the evolutionary games played by657 non-small cell lung cancer”. In: Nature Ecology & Evolution 3.3 (Mar. 2019). Publisher: Nature658 Publishing Group, pp. 450–456.659 [61] Maximilian A.R. Strobl et al. “Turnover Modulates the Need for a Cost of Resistance in660 Adaptive Therapy”. In: Cancer Research 81.4 (Feb. 2021), pp. 1135–1147.661 30 .CC-BY-NC-ND 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 3, 2024. ; https://doi.org/10.1101/2024.10.30.621043doi: bioRxiv preprint