article: Tumors consist of a hierarchical population of cells that differ in their phenotype and genotype. This hierarchical organization of cells means that a few clones (i. e. , cells and several generations of offspring) are abundant while most are rare, which is called clonal dominance. Such dominance also occurred in published in vitro iterated growth and passage experiments with tumor cells in which genetic barcodes were used for lineage tracing. A potential source for such heterogeneity is that dominant clones derive from cancer stem cells with an unlimited self-renewal capacity. Furthermore, ongoing evolution and selection within the growing population may also induce clonal dominance. To understand how clonal dominance developed in the iterated growth and passage experiments, we built a computational model that accurately simulates these experiments. The model simulations reproduced the clonal dominance that developed in in vitro iterated growth and passage experiments when the division rates vary between cells, due to a combination of initial variation and of ongoing mutational processes. In contrast, the experimental results can neither be reproduced with a model that considers random growth and passage, nor with a model based on cancer stem cells. Altogether, our model suggests that in vitro clonal dominance develops due to selection of fast-dividing clones.
  Intratumoral heterogeneity, the genotypic and phenotypic differences within a single tumor, is a well known feature of cancer [1] and strongly influences the effectiveness of cancer therapy [2]. Genotypic heterogeneity is the result of random mutations, and while most of these mutations are neutral passenger mutations, some are functional mutations that add to phenotypic heterogeneity. Phenotypic differences may also be caused by phenomena such as differential signaling from the local tumor micro-environment, epigenetic changes, and stochastic gene expression [3]. Another proposed source of intratumoral, phenotypic heterogeneity is the presence of so-called cancer stem cells (CSCs) that have an unlimited potential to renew and can give rise to differentiated cells (DCs) with a limited potential to renew [4]. The presence of CSCs would result in a tumor containing a mixture of CSCs, and DCs that all derive from a small number of CSCs. For a long time, evidence for the presence of CSCs was primarily based on xenograft models in which transplantation of tumor cells into immunodeficient mice resulted in tumor growth in only a small fraction of the mice [1,5], suggesting that only a subset of the tumor cells has the ability to sustain long-term growth. However, the lack of success of initiating tumor growth in immunodeficient mice may also be related to the incomplete inhibition of the immune response [6], or to the dramatic change in tumor micro-environment upon transplantation [5]. An alternative approach to identify the existence of CSCs is to perform lineage tracing by fluorescent marking of a subpopulation of cells [7,8]. For example, Schepers et al. [9] managed to trace the lineage of CSCs by fluorescently labeling cells expressing stem cell markers in mice developing intestinal adenomas and thereby showed that all cells in small adenomas descended from a single stem cell. However, fluorescent labeling of stem cells is not possible in all cancer types, and for those cancer types an alternative approach is taken by labeling a small fraction of the tumor cells in animal models. Studies employing this strategy showed that the number of colored patches reduced during tumor growth, while the size of these patches increased [10–12]. These observations are compatible with the hypothesis that tumor cells descend from a small number of CSCs. Another, high-resolution, approach to lineage tracing is the application of unique genetic tags, also called cellular barcodes, to a population of tumor cells [13–19]. Tumors grown in immunodeficient mice injected with barcoded tumor cells are dominated by cells that express only a small subset of the barcodes [15,16]. Serial implantation of barcoded leukemic cells showed that rare clones in one mouse can develop dominance after transplantation into a second mouse, indicating that clonal dominance is not a predetermined property of certain clones [18]. Porter et al. [13] used cellular barcoding to follow the development of clonal dominance over time in an in vitro setup, thereby controlling the external factors that could affect clonal dominance. Populations derived from several polyclonal cell lines were barcoded and grown for three days after which a fraction of the cells was passed on to the next generation. By repeating this process (Fig 1A) and analyzing the intermediate clone distributions, Porter et al. [13] showed that clonal dominance progressed over time (Fig 1B). Interestingly, similar experiments with a monoclonal K562 cell line resulted in a minimal progression of clonal dominance, hinting that an intrinsic variability in the cell population may cause the progression of clonal dominance. Moreover, random passage could also play a role in developing clonal dominance. Although the studies discussed above represent a strong base of evidence for the development of clonal dominance in vivo as well as in in vitro tumor cell populations, the mechanism that drives this dominance remains unknown. The presence of CSCs is consistent with the induction of clonal dominance, but only in some cancers direct evidence for the involvement of CSCs is available [5]. Alternatively, evolutionary selection pressure may cause clones with a higher division or survival rate to dominate the tumor. Hence, it is necessary to further investigate the role of these mechanism in the development of clonal dominance. One way to do this in a formal way is to construct computational models that incorporate different hypotheses and compare the outcome of computer simulations to spatio-temporal clonal dynamics observed experimentally. Such an approach has been used before in numerous studies addressing the temporal and spatial evolution of tumor cell populations with CSCs, which are thoroughly reviewed in [20]. Several of these modeling studies focused on the development of spatial tumor heterogeneity by employing cell-based models in which intratumoral heterogeneity is induced by CSCs [21,22], by CSCs and epigenetic changes [23], or by mutations [24]. Here, we built a computational model that simulates the iterated growth and passage experiments described in [13]. By incorporating different hypotheses for the development of clonal dominance in our model, we show that heritable variability of the rates of cell divisions amongst tumor cells, which is due to a combination of initial variation and mutational processes, is sufficient to induce clonal dominance in iterated growth and passage experiments.
  We started by re-analyzing the data for the in vitro experiments with the lentivirally barcoded cell lines, in order to subsequently make quantitative comparisons with simulation results (see Methods section for more details). In these previously published experiments, the barcode-transduced cells were grown and aliquots containing 3 ⋅ 105 cells were used to initialize three biological replicates of the iterated growth and passage experiment. The cells were then, iteratively, grown for 3 days after which 3 ⋅ 105 cells were passed on to the next generation (Fig 1A). In agreement with Porter et al. [13], iterated growth and passage with the polyclonal K562 cell line resulted in clone loss (Fig 1B) and progressing clonal dominance (Fig 1C). To evaluate clonal dominance, we plotted clone size, sorted from large to small, against the cumulative population fraction for clones of that size (Fig 1C). In the remainder of this paper, we quantify clonal dominance based on this graph (Fig 1D) by employing the Gini coefficient [25]: G=12C0∑i=1C0∑j=1C0| Ni−Nj |n with C0 the number of clones at initialization, Ni the size of clone i, and n the number of cells. As expected from the observed clone distributions, the Gini coefficient is already greater than zero at the start of the experiment (indicating mild clonal dominance) and further increases over time (Fig 1E). We also processed the data for the in vitro experiments with monoclonal K562 cells (Fig 1F) and with HeLa cells (Fig 1G) and clone loss was in good agreement with those published in [13] (S1 Fig). The monoclonal K562 cell line exhibits a more limited clone loss and development of clonal dominance (at time points beyond the initial measurement) compared to the polyclonal K562 cell line. The HeLa cells show a similar trend as the K562 cells: the number of unique clones declines while clonal dominance increases. However, clone loss occurred later for the HeLa cells than it did for the polyclonal K562 cells. The limited development of clonal dominance during passaging in monoclonal K562 cells compared to polyclonal K562 cells indicates that chance could play a role in the development of clonal dominance, hence the most straightforward explanation for the development of clonal dominance is that the iterated passages cause small clones to completely disappear while larger clones remain and grow. This hypothesis is supported by Porter’s in vivo experiments in which no passage occurred and no loss of clones nor clonal dominance was observed [13]. To test whether clone loss during passage explains clone loss and progressive clonal dominance, Porter et al. employed a computational model of iterated deterministic growth phases and passage steps in which a subset of the cells were selected at random, showing that this could partially explain the experimentally observed clone loss but not the development of clonal dominance (see Methods section of [13]). In these simulations, all clones grow according to N i (t + t growth) = N i (t) e r t growth, where Ni (t) represents the size of clone i at time t, tgrowth is the duration of the growth phase, and r is the division rate of the cells. At the end of the passage interval, npass cells are selected randomly and passed to the next generation (Fig 1A). As a first step, we confirmed the Porter simulation results using the code that was published with Porter’s paper [13] (https: //github. com/adaptivegenome/clusterseq/). We ran simulations starting with C = 14,000 clones uniformly distributed over 3 ⋅ 105 cells and r = 24 19 log (2) to match the experimentally observed doubling time of 19 hours. As expected, these simulations resulted in a moderate reduction of the number of clones (Fig 2A; blue bars). In contrast to the in vitro results, clonal dominance developed only slightly and hardly increased over time beyond passage 10 (Fig 2B; blue bars). To test if a realistic initial clone distribution improved the resemblance between the experimental observations and the simulation results, we extended the simulations by employing the initial clone distribution of the polyclonal K562 cell line (Fig 2C) in which clonal dominance is already slightly developed. In this setting the expected number of clones left still remained larger than that for the polyclonal K562 cells while the clone distribution was hardly affected (Fig 2A and 2B; green bars). Porter et al. [13] noted that in simulations in which more cells were passaged, fewer clones disappeared, indicating that the random process of passage causes some clones to become smaller and finally disappear. However, because cell division is modeled as a deterministic process, the clone sizes, relative to other clones, remain similar over time, and the probability to disappear for individual clones remains unchanged over time. The assumption of deterministic growth would be reasonable for large clones, but for small clones, probabilistic events could strongly affect the simulation outcome. As shown in Fig 2C, a large part of the clones contain only a few cells, e. g. , ∼ 25% of the clones have less than 10 cells. Therefore, we next explored the effect of stochastic cell division by replacing deterministic growth by a stochastic growth model using Gillespie’s Stochastic Simulation Algorithm (SSA) [26] (see Methods section for details and Table 1 for the model parameters). In contrast with Porter’s deterministic growth model, passage occurs when the population size reaches the critical population size ncrit = 4 ⋅ 106, which corresponds to approximately 3 days of growth with a population doubling time of 19 hours. Because growth continues until the critical population size is reached, the division rate no longer affects the simulation results and we therefore set it to the arbitrary value of 1. As before, npass cells are then chosen randomly and passed on to the next generation. We performed stochastic growth simulations initialized with either a uniform clone distribution (Fig 2A and 2B; purple bars) or the initially observed clone distribution of polyclonal K562 cells (orange bars). In both cases the results were similar to those with deterministic growth, albeit with a larger decrease in the number of clones. For simulations with stochastic growth that were initialized with the polyclonal K562 clone distribution, the clone loss resembles the average in vitro clone loss for passages 10 and 20, but for passage 30 the in vitro clone loss overtook the simulated clone loss. This suggests that the early clone loss is dominated by the effects of stochastic division and passage, and that at a later stage another, unknown, mechanism further increases clone loss. Altogether, these results indicate that passage and growth can cause clone loss during passage, but neither of these mechanisms can induce progressive clonal dominance. The presence of cancer stem cells is thought to induce tumor heterogeneity because they generate a, hierarchically organized, population of cancer stem cells and differentiated cells [27]. To introduce cancer stem cells in our growth model we replaced the unlimited growth in the stochastic model with a previously published model of cancer stem cell driven growth [28]. In this model cells are either cancer stem cells (CSCs) that can divide indefinitely, or differentiated cells (DCs) that divide a limited number of times. CSCs proliferate at a division rate of rCSC and division can result either in two stem cells (with probability p1), in a stem cell and a differentiated cell (with probability p2), or in two differentiated cells (with probability p3) (Fig 3A). Differentiated cells proliferate at a division rate rDC until they reach their maximum number of divisions M, after which, following Weekes et al. [28], they die with a rate rDC. We did not consider random cell death of cancer stem cells and differentiated cells, because this process only affects the population growth rate. We fine-tuned the parameters of the CSC growth model such that the population growth rate is consistent with the 19 hour doubling time reported by Porter et al. [13]. For this we exploited the analytical solution of the CSC model elegantly derived by Weekes et al. [28]. This solution predicts that the population of cells initially grows, and then develops according to one of three growth regimes determined by β = (p1 − p3) rCSC. When β > 0, the population continues to grow, when β = 0 the population reaches an equilibrium, and when β < 0 the population will eventually go extinct. Because the in vitro cells were reported to be in “log phase growth” [13], we limited the parameter space to β > 0 by setting p3 to 0, which ensures a monotonically growing cell population. A complicating factor is that the parameters determining β do not affect the growth in the first couple of days. Instead, during this time interval, growth is determined by the division rate of DCs (rDC) relative to that of CSCs (rDC), and the maximum age of DCs (M). Therefore, we explored how much the simulated population doubling time deviates from the experimental doubling time of 19 hours when M and rDC were varied while the other simulation parameters remained at their default value according to Table 1 (Fig 3B). This exploration resulted in a well-defined parameter range for which the in vitro doubling time could be reproduced in silico. To illustrate the clonal development of the CSC growth model, we tested the model for several parameter sets from the region that matches with the experimentally observed population growth rate (Fig 3B; crosses). The simulations showed that the maximum age of DCs or the CSC division rate hardly affect clonal development (Fig 3C and 3D). Furthermore, an at first sight surprisingly strong reduction in clone number occurred in the simulations while the clone distribution remained virtually unchanged (Fig 3C and 3D). The dramatic reduction in clone number can readily be explained: clones frequently loose all CSCs during passage, which takes away the long-term self-renewal capacity of such clones. We tested this explanation by varying the initial percentage of CSCs and the probability of a symmetric CSC differentiation while keeping the other parameters constant (Table 1 and white cross in Fig 3B). In line with our explanation for the massive clone loss, increasing the probability of symmetric CSC division strongly reduced clone loss (Fig 3E). Increasing the initial percentage of stem cells had a less pronounced effect on clone loss (Fig 3E), which fits with the observation in [28] that the proportion of CSCs in the population will become constant in the long term and does not depend on the initial percentage of CSCs. Altogether, increasing the probability of a symmetric CSC to 1 and the initial CSC percentage to 100% did not lead to clone loss matching that observed experimentally. Moreover, changing these parameters had only little effect on clonal dominance (Fig 3F). The lack of progression of clonal dominance over time can be explained by a lack in difference between the clones: although in the long run all cells descend from only a few clones, there is no difference in the speed at which each clone generates offspring. In conclusion, incorporating CSC growth into our passaging simulation resulted in a poor match to the experimental observations because far too many clones disappeared and the distribution of the remaining clones did not exhibit clonal dominance. Because neither simulations with random growth and passaging of tumor cells nor simulations incorporating CSC growth can explain the development of clonal dominance, we next considered intrinsic variability of growth characteristics as a source of clonal dominance development. To this purpose, we let go of the distinction between cell types and let all cells divide with their own division rate ri, that is subject to mutations. Note that the term mutation is used here to describe any change in a cell’s phenotype that can be inherited, which includes both genetic and epigenetic changes. Mutation of cellular properties requires the explicit representation of single cells, therefore we can no longer use Gillespie’s SSA because clones rather than single cells are explicitly described in that type of simulation. Therefore, we implemented an agent based model (ABM) in which each cell is represented explicitly, and both a barcode and division rate are associated to each cell. In the ABM there are two sources for division rate variability. First, the division rate of cells at initialization may vary due to prior mutations. To implement this, the division rate of each cell is initialized to ri = Xi r0, with ri the cell’s division rate, r0 the mean division rate of all cells and Xi taken from a normal distribution with mean 1 and SD σ 0 *. Second, mutations are included by setting the division rate of each offspring to ri = Yirp, with Yi taken from a normal distribution with mean 1 and SD σ m *, and rp the division rate of the parent cell. Further details on the implementation of the ABM can be found in the Methods. To test if the ABM can induce similar clone loss (Fig 4A) and clonal dominance (Fig 4B) as observed in vitro, we ran simulations for a range of initial division rate SDs (σ 0 *) and mutation SDs (σ m *), all initialized using the initial distribution of the polyclonal K562 cells. Note that we did not vary the initial mean division rate (r0) because this parameter does not affect the simulation results (S2 Fig). The resulting ranges of clone loss and Gini coefficient do indeed include those observed for the polyclonal K562 cells. Furthermore, simulations with division rates that vary due to either initial variation (σ 0 * > 0) or mutations (σ m 0 > 0) result in higher levels of clone loss and a higher Gini coefficient compared to simulations without any variation (σ 0 * = σ m * = 0). Thus, the ABM that describes evolution of division rate variability has the potential to match the in vitro results. To fine-tune the model for the polyclonal K562 cell line measurements, we perform a maximum likelihood estimation (MLE) for the simulations shown in Fig 4 to identify the best fitting values for σ 0 * and σ m *. For this we consider the errors of all metrics included in the MLE to be normally distributed, so we can use the following definition of the loglikelihood: ℓ (σ 0 *, σ m *) ∼ ∑ m ∈ metrics ∑ P ∈ passages - (m ¯ (P) - μ m, P) 2 2 σ m, P 2 with m ¯ (P) the mean of 10 simulations for metric m at passage P, and μm, P and σm, P respectively the mean and SD of the three polyclonal K562 replicates for metric m at passage P. Using passages 10,20, and 30, and either the fraction of clones left, the Gini coefficient, or both metrics, we obtained heatmaps of the likelihood (Fig 5A–5C). The simulated clone loss is quite similar to that for the polyclonal K562 cells (Fig 5A), independent of σ m * and σ 0 *. This fits with the previous observation that a model with stochastic growth, but with identical division rates for all cells, already closely fits clone loss (Fig 4A). In contrast, the Gini coefficient only fits well for a defined region of parameter values (Fig 5B) and as a result this metric dominates the MLE based on both parameters (Fig 5C). The area of best matching parameter sets also includes sets without mutation or initial division rate variation, suggesting that the source of the division rate variability cannot be determined. However, it seems most likely that both mechanisms are active because mutations before the experiment started would already induce initial division rate variation and these mutations continue to happen upon cell division. Consistent with this, the best fit corresponds to a model that includes both mechanisms (σ 0 * = 0. 036 and σ m * = 0. 0018) and the simulation results are close to the experimentally observed Gini coefficient (Fig 5D) and clone loss (Fig 5E). Furthermore, major clones, i. e. , clones representing more than 1% of the population, do develop in the ABM simulations (Fig 5F) while the mean division rate of the population only increases with 10% (Fig 5G). To test the effect of initially present heterogeneity on clone loss and the development of clonal dominance, Porter et al. [13] created a monoclonal K562 cell line and showed that clone loss was reduced and clonal dominance only slightly developed (Fig 1F). Because this cell line was derived from a single cell, it is expected to have a lower initial variation yet similar mutation dynamics as other K562 cells. We tested this hypothesis by running a parameter sweep for σ 0 * and σ m * in which each simulation starts with the initial distribution of the monoclonal K562 cells. Performing an MLE for either clone loss (Fig 6A) or the Gini coefficient (Fig 6B), we found that our model fails to produce a good fit for the clone loss (the maximum ℓ in Fig 6A is -780), while the fit for the Gini coefficient is much better (the maximum ℓ in Fig 6B is -40). In line with our expectations, the optimal fit of the Gini coefficient occurs at a σ 0 * that is clearly lower than the optimal fit for the polyclonal K562 data (cf. Fig 5B). Moreover, the best-fitting σ m * values for mono- and polyclonal K562 data are similar. Although the optimal fitting parameters for the Gini coefficient show a close fit for the Gini coefficient (Fig 6C), this same parameter set produces a clone loss that follows a similar pattern as the monoclonal K562 cells, but the simulated clone loss is much higher (Fig 6D). However, the overestimation of clone loss by our model may actually be an underestimation of in vitro clone loss that is quantified by the PCR and sequencing procedure. Because any clone that had a barcode in the reference library is counted, so-called spurious reads may cause false positives and thus underestimate actual clone loss. These reads are at least partly the result of systematic errors that are common with next-generation sequencing methods [29]. Porter et al. [13] elegantly showed that the sequencing results were reproducible by pair-wise comparing the clone frequencies in four samples of the plasmid library. Nevertheless, this analysis does not completely rule out spurious reads because read errors are correlated to specific base sequences [30] and such errors are expected to occur at similar frequencies in each replicate [31]. To test the effect of spurious reads we artificially ‘contaminated’ the clone sizes returned by the simulations with spurious reads (see details in S1 File). This analysis showed that such spurious reads can decrease clone loss to such an extent that the model matches in vitro clone loss. The effect of spurious reads on the Gini coefficient was less pronounced, arguing in favour of fitting to the Gini coefficient rather than to clone loss data or to a combination of the two metrics. Furthermore, results were similar for the HeLa cells, i. e. for that cell line clone loss was not matched well whereas the Gini coefficient and the number of major clones was (for details see S2 File). Altogether, these results indicate that variation in division rates explains the dynamics of the clone size distribution described by Porter et al. [13].
  Employing a simple model of stochastic cell division and passage, we showed that most of the clone loss taking place during passage, observed by Porter et al. [13], is due to the random loss of clones during passage. In contrast, the progressive clonal dominance that developed in the same experiments, cannot at all be explained by random clone loss. Extending the simulations with indefinitely dividing CSCs and DCs with limited division capacity did not induce progressive clonal dominance and led to a much larger clone loss than observed experimentally. However, a model in which variability of division rates was present and was subject to an evolutionary selection pressure resulted in a progressive clonal dominance matching the dominance observed experimentally. The level of correspondence between clone loss and the simulations varied per tested cell line, which we hypothesize to result from errors during sequencing as these are known to occur frequently and at a reproducible rate [29,31]. We showed that indeed such errors have a much larger effect on the measured clone loss than on the Gini coefficient, which we used to characterize clonal dominance. Altogether, our model provides strong evidence that heritable division rate variation rather than the presence of CSCs induces the changes in clonal distribution observed by Porter et al. [13]. Note that our model was set up to mimic in vitro growth and cannot be directly extrapolated to predict the development of heterogeneity within real tumors, because in that case spatial effects such as contact inhibition of growth and signals from the tumor micro-environment are likely to play a role. Nevertheless, our simulations elucidate which effect evolution of division rates can have on the clonality of a cell population, and such effects are likely to also be important in vivo. For example, several in vivo, microscopy-based, lineage tracing studies [10,11], ascribe the development of large monochromatic patches from a mosaic pattern to CSCs. Based on these observations, Driessens et al. [10] proposed a mathematical model with CSCs and DCs that closely fitted the clone sizes they observed in their experiments with skin papillomas. However, their model considered CSCs to divide twice as fast as DCs, while CSCs are typically thought to divide slower or at a similar rate as DCs [4]. It would therefore be interesting to investigate whether division rate variability could also contribute to such in vivo lineage tracing data. Although there is in general ample evidence for the existence of CSCs [1,3, 4], our model results do not point to a role for CSCs in the experiments of Porter et al. Indeed, in the CSC growth model with division rate heterogeneity, clonal dominance appeared in combination with a massive clone loss. This clone loss occurred because only clones that had a CSC at initialization had any chance of generating offspring in the long term, and even those clones could disappear when CSCs were accidentally lost during passage. Consistent with this explanation, removing the distinction between CSCs and DCs in our model, led to clone loss closely matching the experimental observations. Recent studies showed that the CSC fate is plastic, meaning that differentiated cancer cells sometimes can become CSCs [4,32]. This CSC plasticity might provide an alternative mechanism to prevent clone loss, by enabling clones to (re) acquire CSCs. Adding such plasticity to our CSC growth model would give all clones the potential to generate offspring indefinitely, making it similar to a model in which all cells divide indefinitely, and should thus be able to fit in vitro clone loss. To test if such a model with clonal evolution of division rates would develop clonal dominance we extended the CSC model with heritable heterogeneous division rates (see Methods). Simulation with this model indeed showed that clonal dominance develops (see S3 Fig). However, in order to limit clone loss, it was required to have a high probability of symmetric CSC division and a large fraction of initial CSCs in the simulations. Therefore, we expect that a model with CSC plasticity can only reproduce the in vitro findings when the transition from DC to CSC is common. Reproducing the in vitro results in our simulation was possible when we considered the division rates of tumor cells to vary between clones and to be fully and directly inherited from the parent cell (although the rate could be changed by mutation). Whereas classical studies have provided ample evidence for division rate heterogeneity among tumor cells [33–35], we are unaware of observations of its full and direct inheritance. Nevertheless, Gray et al. [33] showed that when melanoma cells are iteratively grown in mice, isolated, and transferred into new mice, the tumor growth speed increased every generation, which indicates that fast dividing cells generate offspring that also divide fast. More recent work further supports the assumption of heritable division rates by showing a strong, positive, correlation in the division rate of B-cell siblings [36]. However, other studies have shown that the division times of breast cancer cells correspond less between parents and offspring than between siblings [37], and that in lymphoblasts the division times of parents and offspring do not correlate at all [38]. These findings indicate that the child’s division rate is not a direct copy of the parent’s division rate. Sandler et al. [37] propose a kicked cell cycle model where the cell cycle length is determined by the level of an oscillating protein which is inherited from the parent and the phase of this protein determines the time between birth and division. Hence, such a model results in similar division times for siblings, while the correlation between division times of parents and offspring depends on the cell cycle duration [37]. These observations show the need for a better understanding of how the division rate of child cells depend on the parent, which can be achieved from lineage tracing studies employing imaging of multiple divisions over time. While our model explains the strong clonal dominance evolving over time in vitro, it does not perfectly match all in vitro observations of Porter et al. [13]. The main discrepancy is that our model overestimates the clone loss observed with monoclonal K562 cells or HeLa cells (and is also slightly off for polyclonal K562 cells). Clone loss is minimal when we omit all division rate variation (Fig 2A), but even in such simulations clone loss is overestimated. Clone loss can also be reduced by increasing the percentage of clones that is passed (S4A Fig), however this also reduces clonal dominance (S4B Fig). Alternatively, it could be that passaging occurs with a bias for small clones to survive passaging, for instance when clones closely attach to each other in vitro and small groups of attached cells are more likely to be selected for passaging than large groups. However, this explanation seems unlikely considering that trypsinization was sufficient to detach any adhering cells. Therefore, we consider it likely that in vitro clone loss was actually underestimated in the experimental data due to the presence of spurious reads. Our analysis (S1 File) confirmed that our model would better fit the in vitro data if such spurious reads are included in the simulation results. Altogether, in this work we used a computational approach to test different hypotheses for the development of clonal dominance and showed that only the presence and evolution of division rate heterogeneity, can reproduce the experimental observations. Hence, this study showcases the value of computational modeling in the interpretation of experimental results. In the future, the model could be further extended to improve its power, especially for comparison with in vivo data. For this, the model should be extended with an explicit representation of space and physical interactions between cells [39]. With such a model it becomes possible to explore the consequences of division rate variability while comparing with intra-vital images studies.
  The sequencing data were downloaded from the NIH Sequence Read Archive (https: //www. ncbi. nlm. nih. gov/sra/SRX535233) using the SRA Toolkit. The downloaded FASTQ files were processed with our own code (S3 File), which collects all barcodes for which the read quality is 56 or higher and counts the occurrence of each barcode. We used our own code rather than the ClusterSeq code (i. e. , the code developed and used in Porter et al. [13]), because ClusterSeq clusters barcodes for each dataset separately, which may result in identical barcodes ending up in a different cluster at different time steps of the same experiment. To generate a reference library for the cell lines barcoded with the lentiviral vector, we merged the four plasmid library samples, excluding any barcode with a frequency smaller than 0. 0002% or that appeared in only one biological replicate. The resulting reference library, containing 13325 barcodes, was then used to select only these known barcodes from the experimental data (S1 Dataset). Finally, we processed the FASTQ files containing the reads for the experiments with the polyclonal K562, monoclonal K562, and HeLa cell lines. The full processing of these files is outlined in S3 File. In order to understand the clone size dynamics of tumour cells, we employed Gillespie’s Stochastic Simulation Algorithm (SSA) [26] to simulate the growth and passage of a cell population. The SSA allows to follow the size Ni of each clone i over time rather than that of single cells and we applied this to the polyclonal K562 cell line. The simulations were initialized based on the read counts for this cell line at passage 0 (see S2 Dataset), which contains the counts for nlibrary = 13325 unique barcodes. Hence, we initialized Ni for 1 ≤ i ≤ nlibrary as N i (0) = n pass n library N K562, i, where NK562, i is the read count for clone i in the polyclonal K562 sample and npass = 3 ⋅ 105 is the number of cells that are selected during passaging (i. e. , also for the initial passage). The simulations closely follow the experimental procedure of iterated growth and passage: the initial population of 3 ⋅ 105 cells grows until the critical population size ncrit = 4 ⋅ 106 is reached, then npass cells are selected randomly and taken to the next generation. We implemented the SSA using the τ-leaping algorithm [40], in which size Ni of clone i is described by: N i (t + τ) = N i (t) + k i (1) with time step τ, and ki taken from a Poisson distribution with mean rNi (t) τ, where r denotes the division rate. These τ leaps are performed until the population size n = ∑ i = 0 C N i ≥ n pass, where C is the number of clones. We set τ = 0. 0005, which resulted in a simulation time of typically ∼100 seconds (for ∼12,000 clones with one species and r = 1) with similar results as for lower τ values (S5 Fig). The source code of the simulations can be found in S4 File. To test how CSC driven growth affects the clone size evolution, we incorporated a previously published model of CSC driven growth [28]. In this model CSCs proliferate at a division rate of rCSC and division can result either in two CSCs with probability p1, in a CSC and a DC with probability p2, or in DCs with probability p3 (Fig 3A). DCs proliferate at a division rate rDC until they reach their maximum number of divisions M and then, following Weekes et al. [28], they die with a rate rDC. For simplicity, we did not consider random cell death of CSCs and DCs, because this process only affects the population division rate. This division scheme is incorporated in the growth model by defining for each clone i the number of CSCs NCSC, i and the number of DCs NDCm, i of age m (where m can take values ranging from 0 to the maximum age M). Then, for each clone i there are 5 possible transitions: where ki is obtained as described above, except for transitions 1-3 where the CSC division rate is multiplied by the respective transition probability. At initialization, the number of CSCs in each clone is set to NCSC, i = CSC0Ni, rounded to the nearest integer. The remaining cells are distributed evenly over the M DC species, while rounding to the lowest integer: N DC, i = ⌊ 1 M (1 - CSC 0) N i ⌋ ∀ 0 ≤ i ≤ M. When, due to rounding, N CSC, i + ∑ i = 0 M N DC, i < N i, the remaining cells are randomly distributed over all species of clone i. Note that the actual fraction of CSCs and DCs of varying ages will settle to equilibrium values that depend on the CSC and DC division rate and the probabilities for the CSC division modes [28], and that the initial CSC fraction therefore does not strongly affect simulation outcome (see Fig 3E). To test the effect of division rate heterogeneity on the development of the clone size distribution, the division rates ri, CSC and ri, DC of each clone i are multiplied by a randomly chosen value Xi that is obtained from a normal distribution with mean 1 and SD σ r *, and Xi is set to 0 when Xi <0. Note that σ r * is the division rate SD relative to the mean division rate; the actual division rate SD is σ r * · r. Because Xi is linked to clone i, the division rates are inherited upon division. The agent based model (ABM) differs from the SSA model in the explicit representation of individual cells rather than of clones. In the ABM, each cell i has a barcode bi and a division rate ri. All simulations are initialized based on sequencing data of either the polyclonal K562 cells, the monoclonal K562 cells, or the (polyclonal) HeLa cells, which can all be found in S2 Dataset. First, we create a master population from which the initial population of each replicate is selected. The master population consists of 5 ⋅ 106 cells, each of which contains a barcode following read abundance in the sequencing data, such that the relative abundance of the barcodes is initially correct. Then, a division rate ri = r0 max (X, 0) is assigned to each cell in the master population, where X is randomly drawn from a normal distribution with mean 1 and SD σ 0 *, and r0 is the initial mean division rate. The seed of the random number generator used for division rate selection is kept constant across simulation replicates such that the same division rate is associated with the same cell and the same barcode in every replicate. Finally, ninit cells are taken from the initial population to be used in the simulation of iterated growth and passage. Note that the initialization omits the clonal dynamics during the 8 to 9 days before the growth & passage experiments started. During this time interval, a within-clone correlation of the division rates likely develops, which is missing in the simulated cell population. As a result the model does not reproduce the overlap between major clones that was observed by Porter et al. [13]. Division is modeled using the dynamic Monte Carlo method [41] implemented with a first reaction method. Whenever a cell is created, either during initialization or after division, the next division time for that cell is assigned according to t i = - ln R r i, where ri is the cell’s division rate and R is a random number from a uniform distribution between 0 and 1. Then, by creating an ordered list of next division times, population growth can be simulated efficiently. When division occurs, a new division rate rc is assigned to each newly created cell according to rc = rp max (Y, 0) with rp the parent’s division rate and Y taken from a normal distribution with mean 1 and SD σ m *. Note that σ 0 * is the initial division rate SD relative to r0; the actual SD of the initial division rates is σ 0 * · r 0. As in the SSA, division continues until the population size reaches its maximum ncrit after which npass cells are randomly selected and passed to the next generation. The source code of the simulations can be found in S5 File. Note that there are some implementation differences between the aforementioned SSA model and the ABM, which result in small quantitative differences in the simulation results. These differences occur for two reasons: First, growth is slightly underestimated in the SSA model, and the exact underestimation depends on the clone size. Second, a clone in the ABM can include cells with different division rates, which is not the case in Gillespie simulations.
summary: Tumors consist of numerous cell populations, i. e. , clones, that differ with respect to genotype, and potentially with respect to phenotype, and these populations strongly differ in their size. A limited number of clones tend to dominate tumors, but it remains unclear how this clonal dominance arises. Potential driving mechanisms are the presence of cancer stem cells, which either divide indefinitely of differentiate into cells with a limited division potential, and ongoing evolutionary processes within the tumor. Here we use a computational model to understand how clonal dominance developed during in vitro growth and passage experiments with cancer cells. Incorporating cancer stem cells in this model did not result in a match between simulations and in vitro data. In contrast, by considering all cells to divide indefinitely and division rates to evolve due to the combination of division rate variability and selection by passage, our model closely matches the in vitro data.
keywords: medicine and health sciences
  cell division analysis
  hela cells
  cell cycle and cell division
  cell processes
  biological cultures
  cloning
  cancer stem cells
  cell differentiation
  simulation and modeling
  oncology
  developmental biology
  stem cells
  cell cultures
  bioassays and physiological analysis
  molecular biology techniques
  research and analysis methods
  evolutionary rate
  cell analysis
  animal cells
  cell lines
  molecular biology
  carcinogenesis
  cell biology
  biology and life sciences
  cellular types
  cultured tumor cells
  evolutionary biology
  evolutionary processes
title: Heritable tumor cell division rate heterogeneity induces clonal dominance