Drug resistance assays

Abstract

Methods of characterising an evolution of resistance in a population of proliferative cells exposed to an anti-proliferative treatment are disclosed, the method comprising: obtaining, by a processor, experimental data comprising total population size and lineage tracing data measured at one or more time points after the start of exposure of the cells to the anti-proliferative treatment, and fitting, by the processor, one or more first dynamic cell population growth models to the total population size data, and one or more second dynamic cell population growth models to the lineage tracing data Related methods and products are also described.

Description

[0001] DRUG RESISTANCE ASSAYS

[0002] Field of the Disclosure

[0003] The present invention relates to methods for characterising the evolution of resistance in proliferative cells subject to anti-proliferative treatment and particularly, although not exclusively, to methods that use total population size data and lineage tracing data to characterise such evolution, perform drug screening and drug design, and identify resistance biomarkers.

[0004] Background

[0005] Cancer treatment frequently fails due to the evolution of drug-resistant cell phenotypes caused by underlying genetic or non-genetic changes. The origin of these adaptations, their timing and rate of spread is key information for distinguishing the mechanism(s) of drug resistance, yet the dynamics cannot be observed directly.

[0006] The evolution of resistance remains the primary obstacle to the successful treatment of cancers. Initial studies focussed on resistance as a consequence of genetic alterations that produce resistant phenotypes that drive clonal evolution (Ding et al., 2012; Misale et al., 2012; Shi et al., 2014), but recent focus has shifted to non-genetic mechanisms of resistance (often termed ‘plasticity’) that enable cells to rapidly change phenotype (Emert et al., 2021; Hinohara et al., 2018; Rehman et al., 2021) promoting adaptive evolution following the change in selective pressure imposed by treatment. Efforts to tackle treatment resistance require knowledge of the molecular mechanism responsible so that it can be targeted. For example, new generations of targeted drugs use the resistance mechanisms of previous drugs as novel targets (Skoulidis et al., 2021). However, evolutionary-informed treatment strategies are also reliant on an understanding of the behaviour of resistant phenotypes: approaches to “steer” tumour evolution to forestall or prevent resistance emergence require an understanding of the heritability of resistance mechanism(s) through cell divisions (Acaret al., 2020), and the relative fitness of resistant cells compared to drug sensitive counterparts (West et al., 2018).

[0007] Resistance evolution offers a unique opportunity to study phenotypic evolution in cancer: treatment is one of the few environmental changes a tumour experiences with known timing that is under clinical control, and resistance is clearly defined as a phenotype that survives treatment. These features mean studies into resistance are well placed to distinguish evolutionary behaviours such as the stability of a resistance phenotype through cell divisions and the environmental dependence of phenotypic change (Luria and Delbruck, 1943; Russo et al., 2022; Shaffer et al., 2020). These behaviours dictate how resistance is distributed amongst cell lineages over time which can in turn determine the relatedness of surviving cells.

[0008] In a patient tumour, naturally occurring somatic mutations can be interpreted as genetic barcodes (Gabbutt et al., 2022b). In an experimental setting, lineage tracing technologies enable the tracking of cell relatedness (Kebschull and Zador, 2018): unique genetic sequences are incorporated into cells’ genomes via lentivirus infection, meaning all subsequent ancestors of the parental, barcoded population’s cells inherit this experimentally measurable tag. Quantitative assessment of barcode data, e.g. using the ClonTracertechnology (Bhang et al., 2015) enables the measurement of clonal dynamics (Blundell et al., 2019). When populations of cells are exposed to a change in extrinsic selection pressures, these genetic lineage tracing data provide a means to understand phenotype dynamics. For example, whether dominant lineages are shared between related experimental populations following exposure to drug treatment has been used as a qualitative indicator of whether drug resistance occurs via genetic or non-genetic mechanisms (Eyler et al., 2020).

[0009] However, there remain a lack of approaches that provide a quantitative readout of these phenotype behaviours.

[0010] Summary of the Disclosure

[0011] The present inventors have identified that there was a lack of approaches to characterise drug treatments in terms of resistance evolution in a pre-clinical context. They have developed a new approach that can infer the characteristics and dynamics of drug resistance without the need for direct measurement of the resistance phenotype using only genetic lineage tracing and population size data. Through extensive simulation experiments, they showed that the framework can recover ground-truth evolutionary dynamics from lineage tracing data. They further used experimental evolution to 5-Fu chemotherapy in two common colorectal cancer cell lines SW620 and HCT116 to provide empirical demonstration of the veracity of the framework. In SW620 cells, a stable pre-existing resistant subpopulation was inferred, whereas in HCT116 cells resistance emerged through phenotypic switching into a slow growing resistant state with stochastic exiting into a fully resistant phenotype. Extensive functional assays, including scRNA-seq and scDNA-seq were used to validate the distinct evolutionary routes and their molecular nature. The proposed approach relies on a new mathematical framework that can be applied to diverse experimental designs to infer the evolutionary dynamics of cancer cell therapy resistance evolution from readily obtained experimental data, enabling more rapid characterisation of resistance mechanisms. The approach as demonstrated relies on the use of both genetic lineage tracing and population size data to infer the temporal dynamics of cancer cell drug resistance phenotypes, without requiring specific measurement of cell phenotypes.

[0012] Thus, according to a first aspect, there is provided a method of characterising the evolution of resistance of a population of proliferative cells to an anti-proliferative treatment, the method comprising: obtaining, by a processor, experimental data comprising total population size and lineage tracing data measured at one or more time points after the start of exposure of the cells to the anti-proliferative treatment; fitting, by the processor, one or more first dynamic cell population growth models to the total population size data, thereby identifying a first set of values of each of one or more parameters of each of the one or more dynamic cell population growth models that meet a first predetermined model fit criterion, wherein the dynamic cell population growth models represent the growth of a plurality of subpopulations of cells having different responses to the anti-proliferative treatment comprising at least a first subpopulation of cells that is sensitive to the anti-proliferative treatment and a second population of cells that is resistant to the anti-proliferative treatment; and fitting, by the processor, one or more second dynamic cell population growth models to the lineage tracing data, each second dynamic cell population growth model representing the growth of the same plurality of subpopulations of cells and having the same parameters as a corresponding first dynamiccell population growth model, said fitting comprising evaluating the fit to data of the one or more second dynamic cell population growth models to the lineage tracing data parameterised using parameter values derived from the first set of values, thereby identifying a second set of values of each of one or more parameters of each of the one or more first and second dynamic cell population growth models that meet a second predetermined model fit criterion; wherein the second set of values of one or more parameters comprise: values that are indicative of the effect of the treatment on the growth of the plurality of subpopulations of cells, values that are indicative of the proportions of the one or more subpopulations of cells prior to exposure to the anti-proliferative treatment, and / or values that are indicative of the rate of transition of cells between the plurality of subpopulations of cells comprising a rate of transition between sensitive and resistant subpopulations of cells.

[0013] The method may have any one or more of the following optional features.

[0014] For each second dynamic cell population growth model there may be a corresponding first dynamic cell population growth model that represents the growth of the same plurality of subpopulations of cells using the same parameters, wherein the second dynamic cell population growth models represents the temporal dynamics of the numbers of cells of each of a plurality of cell lineages that belong to each of the plurality of subpopulations of cells (nis, niR, niE), and wherein the first dynamic cell population growth model represents the temporal dynamics of the total number of cells in each of the plurality of subpopulations of cells (ns, HR, nE).

[0015] In embodiments, the second set of values of one or more parameters comprise values that are indicative of one or more of: a proportion of cells that are resistant before the population of cells is exposed to the anti-proliferative treatment (p), a rate of transition of sensitive cells to resistant cells ( ), a probability of treatment-induced death for resistant cells (qj), a fitness penalty of resistant cells in the absence of treatment (5), a rate at which resistant cells escape a fitness penalty in the absence of treatment (a), a rate of change of a treatment induced effect on a death rate of resistant and / or sensitive cells upon exposure to the treatment(K), and a maximum strength of a treatment induced effect on a death rate of resistant and / or sensitive cells (De). In embodiments, the value of a rate of transition of sensitive cells to resistant cells (p) is indicative of whether transition between sensitive and resistant cells is likely to be due to one or more genetic mechanisms or one or more non-genetic mechanisms.

[0016] In embodiments, the method further comprises: obtaining, by a processor, experimental data comprising total population size and lineage tracing data measured at one or more time points after the start of exposure of the cells to a second anti-proliferative treatment; repeating the step of fitting the one or more first dynamic cell population growth models using the total population size data for the second antiproliferative treatment and the step of fitting the one or more second dynamic cell population growth models using the lineage tracing data for the second anti-proliferative treatment; and fitting one or more third dynamic cell population growth models jointly to lineage tracing data comprising the lineage tracing data for the two anti-proliferative treatments, wherein each third dynamic cell population growth model represents the growth of a plurality of subpopulations of cells comprising at least a first subpopulation of cells that is sensitive to both anti-proliferative treatments, a second population of cells that is resistant to the first anti-proliferative treatment, a third population of cells that is resistant to the second anti-proliferativetreatment, and a fourth population of cells that is resistant to both of the first and second anti-proliferative treatments. In some such embodiments, fitting the cross-resistance model comprise estimating, using both lineage tracing data and the parameters estimated by fitting the second dynamic cell population models independently for the two different anti-proliferative treatments, the value of a cross-resistance parameter (0), wherein the cross-resistance parameter (0) is a parameter that quantifies the strength of statistical association of inheritance of resistance phenotypes to the two anti-proliferative treatments. In embodiments, said fitting comprises evaluating the fit to data of the one or more third dynamic cell population growth models to the lineage tracing data parameterised using parameter values derived from the second sets of values for each of the two anti-proliferative treatments, thereby identifying a third set of values of the cross-resistance parameter that meet a third predetermined model fit criterion.

[0017] In embodiments, fitting, by the processor, the one or more first dynamic cell population growth models to the total population size data and / or fitting, by the processor, the one or more second dynamic cell population growth models to the lineage tracing data, and / or fitting, by the processor, the one or more third dynamic cell population growth models to the lineage tracing data, comprise using an approximate Bayesian computation inference approach, and / or using a simulation-based likelihood free approach and / or using a likelihood free method that relies on comparing simulated data to observed data using a distance measure and generating a posterior distribution of parameters that produce simulations that are within a predetermined distance of the empirical data and so satisfy the first and / or second predetermined model fit criteria and / or third predetermined model fit criteria.

[0018] In embodiments, fitting, by the processor, the one or more first dynamic cell population growth models to the total population size data comprises: simulating the one or more first dynamic cell population growth models using a first plurality of sets of parameter values for each of the one or more parameters of each first dynamic cell population growth model, the first plurality of sets of parameter values comprising the first set of values, thereby obtaining simulated total cell population size data for each of the first dynamic cell population growth models and each of the first plurality of sets of values for each of the one or more parameters; determining whether each set of values of parameters of each first dynamic cell population growth model meets the first predetermined model fit criterion by comparing the corresponding simulated total cell population size data to the measured total cell population size data; and selecting the first set of values as a subset of the plurality of sets of values that meet the first predetermined model fit criterion.

[0019] In embodiments, fitting, by the processor, the one or more first dynamic cell population growth models to the total population size data and / or fitting, by the processor, the one or more second dynamic cell population growth models to the lineage tracing data and / or fitting, by the processor, the one or more third dynamic cell population growth models to the lineage tracing data, comprise using a likelihood free method that relies on training a conditional density estimator to estimate a probability density function for the parameters of the model given simulated data, and using the trained model to predicta posterior distribution of parameters based on corresponding observed data instead of the simulated data, the posterior distribution thereby satisfying the respective predetermined model fit criterion.In embodiments, fitting, by the processor, the one or more second dynamic cell population growth models to the lineage tracing data comprises: simulating the one or more second dynamic cell population growth models using a second plurality of sets of parameter values for each of the one or more parameters of each second dynamic cell population growth model, the second plurality of sets of parameter values derived from the first set of values, thereby obtaining simulated lineage tracing data for each of the second dynamic cell population growth models and each of the second plurality of sets of values; determining whether each set of values of parameters of each second dynamic cell population growth model meets the second predetermined model fit criterion by comparing the corresponding simulated lineage tracing data to the measured lineage tracing data; and selecting the second set of values as a subset of the second plurality of sets of values that meet the second predetermined model fit criterion.

[0020] In embodiments, the first predetermined model fit criterion applies to the value of a distance metric between the measured total population size data and simulated total population size data corresponding to a first dynamic cell population growth model and a set of values for the one or more parameters of the first dynamic cell population growth model. In embodiments, the second predetermined model fit criterion applies to the value of a distance metric between the measured lineage tracing data and simulated lineage tracing data corresponding to a second dynamic cell population growth model and a set of values for the one or more parameters of the second dynamic cell population growth model.

[0021] In embodiments, fitting, by the processor, the one or more second dynamic cell population growth models to the lineage tracing data comprises: simulating the one or more second dynamic cell population growth models using a second plurality of sets of parameter values for each of the one or more parameters of each second dynamic cell population growth model, the second plurality of sets of parameter values derived from the first set of values thereby obtaining simulated lineage tracing data for each of the second dynamic cell population growth models and each of the second plurality of sets of values; training one or more machine learning models to predict conditional density distributions of the parameters given the simulated lineage tracing data, and using the trained one or more machine learning models to predict a posterior distribution of the parameters given the observed lineage tracing data any set of values drawn from these posterior distributions meeting the second predetermined model fit criterion.

[0022] In embodiments, fitting, by the processor, the one or more third dynamic cell population growth models to the lineage tracing data comprises: simulating the one or more third dynamic cell population growth models using a third plurality of sets of parameter values for each of the one or more parameters of each third dynamic cell population growth model, the third plurality of sets of parameter values derived from the second set of values for each of the first and second anti-proliferative treatments, thereby obtaining simulated lineage tracing data for each of the third dynamic cell population growth models and each of the third plurality of sets of values; determining whether each set of values of parameters of each third dynamic cell population growth model meets the third predetermined model fit criterion by comparing the corresponding simulated lineage tracing data to the measured lineage tracing data; and selecting the third set of values as a subset of the third plurality of sets of values that meet the third predetermined model fit criterion.In embodiments, fitting, by the processor, the one or more third dynamic cell population growth models to the lineage tracing data comprises: simulating the one or more third dynamic cell population growth models using a third plurality of sets of parameter values for each of the one or more parameters of each second dynamic cell population growth model, the third plurality of sets of parameter values derived from the second sets of values for each of the first and second anti-proliferative treatments, thereby obtaining simulated lineage tracing data for each of the third dynamic cell population growth models and each of the third plurality of sets of values; training one or more machine learning models to predict conditional density distributions of the parameters given the simulated lineage tracing data, and using the trained one or more machine learning models to predict a posterior distribution of the parameters given the observed lineage tracing data any set of values drawn from these posterior distributions meeting the third predetermined model fit criterion.

[0023] In embodiments, the one of more first dynamic cell population growth models and the one or more second dynamic cell population growth models, and optionally the third dynamic cell population growth models when used, are birth-death models that represent the number of cells alive in each of a plurality of subpopulations of cells as a function of time, the plurality of subpopulations of cells having different responses to the anti-proliferative treatment comprising at least a first subpopulation of cells that is sensitive to the anti-proliferative treatment and a second population of cells that is resistant to the anti-proliferative treatment, wherein the different subpopulations of cells differ by the birth rate in the presence of the treatment, their death rate in the presence of the treatment, their birth rate in the absence of treatment, and / or their death rate in the absence of treatment. In embodiments, the one of more first dynamic cell population growth models are hybrid differential equation and stochastic jump process models. In embodiments, the one of more second dynamic cell population growth models are stochastic agent-based models. In embodiments, a hybrid differential equation and stochastic jump process model is a model that represents the growth of each subpopulation of cells using a stochastic process when the total number of cells in the subpopulation is below a predetermined threshold, and using a deterministic differential equation model when the total number of cells in the subpopulation is at or above the predetermined threshold. In embodiments, the one of more third dynamic cell population growth models are stochastic agent based models. In embodiments, at least one first dynamic cell population growth model and corresponding second dynamic cell population growth model represent the growth of a subpopulation of cells that are resistant to the treatment and do not incur a fitness penalty in the absence of treatment.

[0024] In embodiments, the sensitive cells divide with a birth rate bs and die in the absence of treatment with a base death rate ds, and the resistant cells divide with a birth rate bR = bs (1 - 5) and die in the absence of treatment with a base rate d = ds (1 - 6), wherein 5 represents a fitness penalty of resistant cells in the absence of treatment. In embodiments, the sensitive cells die in the presence of treatment with an effective death rate equals to the sum of a base death rate (ds) and a sensitive treatment effect term (D(t)), and the resistant cells die in the presence of treatment with an effective death rate equals to the sum of a base death rate (dR) and a resistant treatment effect term (D(t)(1 - < / )), wherein the resistant treatment effect term is equal to the sensitive treatment effect term multiplied by the complement of a probability of treatment-induced death for resistant cells (< ). In embodiments, the sensitive treatment effect term is a time dependent parameter that is the product of a maximum strength of a treatment induced effect on a deathrate of cells (De) and a time-dependent factor (y(t)) that increases at a first predetermined rate (K) upon exposure to the treatment and decreases at a second predetermined rate (-K) upon withdrawing of the treatment. In embodiments, the first and second predetermined rate have the same absolute value.

[0025] In embodiments, the first and second dynamic cell population growth models are birth-date models that capture births and deaths of cells in each of the plurality of subpopulations of cells, wherein transitions between specific subpopulations of cells are associated with respective rates and coupled to birth events. In embodiments, the rate of transition between a first subpopulation and a second subpopulation is a parameter between 0 and 1 that represents the probability of a cell in the first subpopulation generating a cell in the second subpopulation upon cell division. In embodiments, the value of the rate of transition between a first subpopulation and a second subpopulation is indicative of the likely mechanism of action of the transition, wherein lower transition rate values are more likely to be associated with genetic or epigenetic mechanisms than higher transition rates, and higher transition rates are more likely to be associated with non-genetic mechanisms than lower transition rates.

[0026] In embodiments, the first and second dynamic cell population growth models each comprise a rate of transition of sensitive cells to resistant cells (p). In embodiments, at least one first dynamic cell population growth model and at least one corresponding second dynamic cell population growth model further comprise a rate of reversion of transition of resistant cells to sensitive cells (a). In embodiments, at least one first dynamic cell population growth model and corresponding second dynamic cell population growth model represent the growth of a subpopulation of cells that are resistant to the treatment and do not incur a fitness penalty in the absence of treatment, and a subpopulation of cells that are resistant to the treatment and do incur a fitness penalty in the absence of treatment, and said first and second dynamic cell population growth models comprise a rate of transition of resistant cells that do incur a fitness penalty in the absence of treatment escape to resistant cells that do not incur a fitness penalty in the absence of treatment. In embodiments, the rate of transition is obtained as the product of a base rate (a) and a time-dependent factor (y(t)) that increases at a first predetermined rate (K) upon exposure to the treatment and decreases at a second predetermined rate (-K) upon withdrawing of the treatment. In embodiments, the first and second predetermined rates have the same absolute value.

[0027] In embodiments, the total cell count data is data that has been obtained by cell counting of a cell population in cell culture, optionally using imaging data, or based on a tumour size determination. In embodiments, the lineage tracing data comprises counts of cells that have the same clonal identity or data derived therefrom, optionally wherein said counts are based on the number of cells that share each of a plurality of heritable barcodes. In embodiments, the lineage tracing data comprises diversity statistics derived from counts of cells that have the same clonal identity. In embodiments, the experimental data comprises total population size and lineage tracing data measured at one or more time points after the start of exposure of the cells in each of a plurality of replicates of the cell population exposed to the anti-proliferative treatment, and the diversity statistics comprise a first statistic that is representative of lineage diversity within the replicates and a second statistic that is representative of lineage diversity between the replicates. In embodiments, the diversity statistics are Hill Diversity statistics. In embodiments, the lineage tracing data comprises first lineage tracing data for a first anti-proliferative treatment and second lineage tracing datafor a second anti-proliferative treatment, and the diversity statistics comprise one or more diversity statics that quantify the degree to which shared lineages are enriched between the lineage tracing data associated with exposure to the first and second anti-proliferative treatments.

[0028] In embodiments, the experimental data comprises total population size and lineage tracing data measured at the start of exposure of the cells to the anti-proliferative treatment (tO) and at one or more time points (P1, P2,...) after the start of exposure of the cells to the anti-proliferative treatment; and / or wherein the experimental data comprises total population size data measured at a plurality of time points after the start of exposure of the cells to the anti-proliferative treatment (01, 02,..., P1, P2,...) and lineage tracing data measured at a subset of the same plurality of time points after the start of exposure of the cells to the antiproliferative treatment (P1, P2,...).

[0029] In embodiments, the experimental data comprises total population size and lineage tracing data measured for a cell population in vitro. In embodiments, the experimental data comprises total population size and lineage tracing data measured fora plurality of replicates of a cell population derived from the same parent population and exposed to the anti-proliferative treatment in vitro. In embodiments, the total population size data comprises total cell counts acquired using a non-disruptive technology, optionally imaging. In embodiments, the lineage tracing data comprises counts of cells in respective clonal lineages acquired using a sequencing technology. In embodiments, the lineage tracing data is data acquired at passaging steps applied to the in vitro cell culture. In embodiments, the method comprises fitting a plurality of first dynamic cell population growth models and a corresponding plurality of second dynamic cell population growth models, wherein the models represent different assumptions of how resistance occurs, and wherein the method further comprises comparing the plurality of fitted second dynamic cell population growth models using a third predetermined model fit criterion, and identifying a second dynamic cell population growth model that best fits the experimental data based on said comparison. In embodiments, the third predetermined model fit criterion is a model fit criterion that takes into account the numbers of fitted parameters of the models that are being compared. In embodiments, the third predetermined model fit criterion is a DIC score. In embodiments, growth in the first dynamic population growth models is modelled as density-dependent growth, optionally logistic growth with a fixed carrying capacity (K).

[0030] In embodiments, the anti-proliferative treatment comprises one or more anti-proliferative compounds or compositions, wherein the one or more anti-proliferative compounds or compositions comprise one or more active agents selected from: small molecules, large molecules, optionally antibodies, antigen binding molecules, peptides, nucleic acids, optionally miRNA, siRNAs and gene editing constructs, combinations or large and small molecules, optionally selected from antibody drug conjugates, bicycle peptide-drug conjugates, cells, optionally T cells or modified T cells, optionally CAR-T cells orTCR-T cells, and targeted protein degradation therapies, optionally selected from PROteolysis TArgeting Chimeras (PROTACs), molecular glues, Lysosome-Targeting Chimaeras (LYTACs), and Antibody-based PROTACs (AbTACs). In embodiments, the first dynamic population growth models comprise a model (A) that represents a first subpopulation of cells that is sensitive to the anti-proliferative treatment and a second population of cells that is resistant to the anti-proliferative treatment., and wherein the sensitive cells can transition toresistance cells but resistant cells cannot transition to sensitive cells. In some such embodiments, the resistant cells incur a fitness penalty in the absence of the anti-proliferative treatment. In embodiments, the first dynamic population growth models comprise a model (B) that represents a first subpopulation of cells that is sensitive to the anti-proliferative treatment and a second population of cells that is resistant to the anti-proliferative treatment, wherein the resistant cells incur a fitness penalty in the absence of the antiproliferative treatment, and wherein the sensitive cells can transition to resistance cells and resistant cells can transition to sensitive cells. In embodiments, the first dynamic population growth models comprise a model (C) that represents a first subpopulation of cells that is sensitive to the anti-proliferative treatment, a second population of cells that is resistant to the anti-proliferative treatment and that incur a fitness penalty in the absence of the anti-proliferative treatment, and a third population of escaped cells that is resistant to the anti-proliferative treatment and that do not incur a fitness penalty in the absence of the anti-proliferative treatment, wherein the sensitive cells can transition to resistance cells and resistant cells can transition to sensitive cells, and wherein the resistant cells can transition escaped cells but escaped cells cannot transition to resistant cells.

[0031] In embodiments, the first dynamic population growth models comprise a parameter that represents a proportion of cells that are resistant before the population of cells is exposed to the anti-proliferative treatment (p). In embodiments, there are no escaped cells before exposure to the anti-proliferative treatment. In embodiments, the first dynamic population growth models comprise subpopulation-specific birth and death rates. In embodiments, the first dynamic population growth models comprise a parameter (5) that represents a fitness penalty associated with the resistant phenotype in the absence of the antiproliferative treatment. In embodiments, the first dynamic population growth models comprise a rate of transition of sensitive cells to resistant cells (p). In embodiments, the first dynamic population growth models comprise a parameter that represents a probability of treatment-induced death for resistant cells (qj). In embodiments, qj = 1.0 denotes complete resistance, whereas when ip = 0.0 the sensitive and resistant cells experience the same level of drug-induced death. In embodiments, the first dynamic population growth models comprise a variable that represents the strength of a treatment induced effect on a death rate of resistant and / or sensitive cells. In embodiments, said variable is dependent on a rate of change of a treatment induced effect on a death rate of resistant and / or sensitive cells upon exposure to the treatment (K), and a maximum strength of a treatment induced effect on a death rate of resistant and / or sensitive cells (De).

[0032] In embodiments, the first dynamic population growth models comprise a model (A) that is a hybrid of a deterministic model specified by equations (7)-(8) below and a stochastic jump model specified by equation (9) below, and the second dynamic population growth models comprise a stochastic agent based model corresponding to the model specified by equation (9) below, where S is a sensitive cell, R is a resistant cell, nS is the number of sensitive cells, nR is the number of resistant cells, t is time, bs is a birth rate of sensitive cells, ds is a death rate of sensitive cells, bR is a birth rate of resistant cells, dR is a death rate of resistant cells, K is a carrying capacity, N is the total number of cells given by N=ns+nR, and p is a rate of transition of sensitive cells to resistant cells. In embodiments, the first dynamic population growth models comprise a model (B) that is a hybrid of a deterministic model specified by equations (10)-(11 ) below and a stochastic jump model specified by equation (12) below, and the second dynamic population growth models comprisea stochastic agent based model corresponding to the model specified by equation (12) below, where S is a sensitive cell, R is a resistant cell, nS is the number of sensitive cells, nR is the number of resistant cells, t is time, bs is a birth rate of sensitive cells, ds is a death rate of sensitive cells, bR is a birth rate of resistant cells, dR is a death rate of resistant cells, K is a carrying capacity, N is the total number of cells given by N=ns+nR, is a rate of transition of sensitive cells to resistant cells, and a is a rate of transition of resistant cells to sensitive cells. In embodiments, the first dynamic population growth models comprise a model (C) that is a hybrid of a deterministic model specified by equations (13)-(15) below and a stochastic jump model specified by equation (16) below, and the second dynamic population growth models comprise a stochastic agent based model corresponding to the model specified by equation (16) below, where S is a sensitive cell, R is a resistant cell, E is an escaped cell, nS is the number of sensitive cells, nR is the number of resistant cells, nE is the number of escaped cells t is time, bs is a birth rate of sensitive cells, ds is a death rate of sensitive cells, bR is a birth rate of resistant cells, dR is a death rate of resistant cells, bE is a birth rate of escaped cells, dE is a death rate of escaped cells, K is a carrying capacity, N is the total number of cells given by N=ns+nR+nE, is a rate of transition of sensitive cells to resistant cells, a is a rate of transition of resistant cells to sensitive cells, and ay(t) is a rate of transition of resistant cells to escaped cells, where y(t) is the time dependent factor that captures the increase or decrease in effective drug concentrations. In embodiments, y(t) is given by equation (4) and / or wherein o+a<1.0.

[0033] In embodiments, one or more of the first dynamic population growth models represent an effect of the antiproliferative treatment on cells in each subpopulation using an effective death rate that depends on a base death rate and a treatment effect term, wherein the treatment effect term is a time depend term that represents an exposure dependent effect of the anti-proliferative treatment on the cells death rate. In embodiments, the treatment effect term comprises a time dependent term that represents the effective concentration of the anti-proliferative treatment and a subpopulation specific factor that is a parameter that captures the subpopulation specific effect of the anti-proliferative treatment on the cell subpopulation. In embodiments, the method further comprises exposing one or more replicates of the cell population in an in vitro cell culture to the anti-proliferative treatment and measuring total population size and lineage tracing data at one or more time points after the start of exposure of the cells to the anti-proliferative treatment. In embodiments, said measuring comprises obtaining both total population size and lineage tracing data at one or more first time points after the start of exposure of the cells to the anti-proliferative treatment. In embodiments, said measuring comprises obtaining total population size data at one or more second time points after the start of exposure of the cells to the anti-proliferative treatment.

[0034] In embodiments, the obtained lineage tracing data and total population size data have been obtained from a cell population cultured in vitro, wherein the cell population was cultured using a process comprising: obtaining lineage tracing data for a parental population of cells, optionally wherein the parental population of cells has been obtained by barcoding a population of cells and expanding the barcoded population of cells during a predetermined period of time, separating the parental population of cells into a first plurality of replicates, and optionally a second plurality of replicates, exposing each of the first plurality of replicates to the anti-proliferative treatment for one or more predetermined periods of time, optionally separated by one or more predetermined periods of time in the absence of the anti-proliferative treatment, and passaging the cells in each of the first plurality of replicates, and optionally each of the second plurality of replicates,one or more times, optionally wherein lineage tracing data is acquired at passaging time points. In embodiments, the method further comprises exposing one or more replicates of the cell population in an in vitro cell culture to a first and second anti-proliferative treatments in respective replicates derived from the same parental population of barcoded cells, and / or receiving previously obtained data from these cells. In embodiments, fitting the first and second plurality of dynamic population growth models comprises sampling cells from a parental population comprising a proportion of resistant cells set by a fitted parameter p to obtain each of a first plurality of replicates comprising a number of cells matching the number of cells in the first plurality of replicates, and optionally each of a second plurality of replicates comprising a number of cells matching the number of cells in the second plurality of replicates, and sampling cells from each replicate at each passage time points to obtain a number of cells in each replicate corresponding to the number of cells in the replicate after passaging. In embodiments, the obtained lineage tracing data and total population size data have been obtained from a cell population cultured in vitro, wherein the cell population was cultured using a process comprising: obtaining lineage tracing data fora parental population of cells, wherein the parental population of cells has been obtained by barcoding a population of cells and expanding the barcoded population of cells during a predetermined period of time, and using a distribution of lineages obtained from said lineage tracing data from the parental population to estimate a parental population average birth rate and a parental population average death rate.

[0035] In embodiments, the method further comprises identifying one or more likely mechanisms of resistance using the parameters of (optionally selected) fitted second dynamic cell population growth models. In embodiments, the method further comprises selecting one or more resistance validation experiments based on the one or more likely mechanisms identified. In embodiments, the one or more mechanisms comprise a genetic mechanism and the selected one or more resistance validation experiments comprise a genome sequencing step, optionally whole genome sequencing or single cell sequencing. In embodiments, the one or more mechanisms comprise a non-genetic mechanism and the selected one or more resistance validation experiments comprise a transcriptome and / or chromatin state sequencing step, optionally single cell RNA sequencing step and / or a single cell ATAC-seq step. In embodiments, the one or more mechanisms include the presence of a subpopulation of cells prior to exposure to the treatment that is resistant to the treatment and experiences a fitness penalty in the absence of the treatment, and the selected one or more resistance validation experiments comprises isolating said subpopulation of cells and performing one or more of a genome sequencing step, a transcriptome sequencing step, and / or a chromatin state sequencing step.

[0036] In embodiments, the method further comprises outputting a report comprising one or more results of the method or information derived therefrom. In embodiments, the one or more results of information derived therefrom include one or more of: a value of one or more fitted parameters, a value of a fit criterion, an indication of an identified dynamic cell population growth model that best fits the experimental data, and an indication of a likely mechanism of resistance.

[0037] Also described according to a second aspect is a method of screening a plurality of candidate antiproliferative treatments, the method comprising: (i) characterising an evolution of resistance in a population of proliferative cells exposed to each anti-proliferative treatment using the methods of any embodiment ofthe first aspect; and (ii) comparing the results of said characterising to identify anti-proliferative treatments of the candidate anti-proliferative treatments that are less likely to be associated with an evolution of resistance.

[0038] The method according to the current aspect may have any of the features described in relation to the first aspect. The method may further comprise selecting one or more candidate anti-proliferative treatments for further characterisation, using the results of said comparing. The further characterisation may comprise pre-clinical and / or clinical characterisation.

[0039] Also described according to a third aspect is a method of identifying a combination of anti-proliferative treatments that has a reduced risk of evolution of resistance compared to a subset of the anti-proliferative treatments in the combination, the method comprising: (i) characterising an evolution of resistance in a population of proliferative cells exposed to the subset of the anti-proliferative treatments using the methods of any embodiment of the first aspect; and (ii) characterising an evolution of resistance in a population of cells exposed to one or more combinations of anti-proliferative treatments comprising respective additional anti-proliferative treatments in addition to the subset of the anti-proliferative treatments using the methods of any preceding claim and comparing the results of said characterising to identify combinations of antiproliferative treatments that are less likely to be associated with an evolution of resistance than said subset of anti-proliferative treatments; and / or (iii) identifying and performing one or more validation experiments based on the results of the characterising in (i), and identifying one or more additional anti-proliferative treatments that target a resistance evolution mechanism identified using said validation experiments. The method according to the current aspect may have any of the features described in relation to the first aspect.

[0040] Also described according to a fourth aspect is a method of identifying one or more biomarkers of resistance of a population of proliferative cells to an anti-proliferative treatment, the method comprising: (i) characterising an evolution of resistance in one or more different populations of proliferative cells exposed to the anti-proliferative treatment using the methods of any embodiment of the first aspect; and (ii) identifying one or more likely mechanisms of resistance using the parameters of (optionally selected) fitted second dynamic cell population growth models for each of the one or more different populations of proliferative cells and performing one or more resistance validation assays based on the identified likely mechanisms of resistance to identify one or more biomarkers of resistance; and / or (iii) comparing the results of said characterising between a plurality of the different populations of proliferative cells to identify one or more biomarkers associated with one or more of the plurality of populations of proliferative cells likely to be associated with an evolution of resistance.

[0041] According to a further aspect, there is provided one or more non-transitory computer readable media comprising instructions that, when executed by one or more processors, cause the one or more processors to perform the steps of any method described herein, such as a method according to any embodiment of any preceding aspect.

[0042] According to a further aspect, there is provided a computer program comprising code which, when the code is executed on a computer, causes the computer to perform the steps of any method described herein,such as a method according to any embodiment of any method described herein, such as a method according to any embodiment of any preceding aspect.

[0043] According to a further aspect, there is provided system comprising: a processor; and a computer readable medium comprising instructions that, when executed by the processor, cause the processor to perform the steps of any method described herein, such as a method according to any embodiment of any method described herein, such as a method according to any embodiment of any preceding aspect. The system may further comprise one or more of: an in vitro cell culture system, a cell count measuring means, optionally a microscope, and a lineage tracing data acquisition means, optionally a sequencer.

[0044] According to a further aspect, there is provided kit comprising: a lineage tracking composition and a computer readable medium comprising instructions that, when executed by a processor, cause the processor to perform the steps of any method described herein, method described herein, such as a method according to any embodiment of any preceding aspect.

[0045] The invention includes the combination of the aspects and preferred features described except where such a combination is clearly impermissible or expressly avoided.

[0046] Summary of the Figures

[0047] Embodiments and experiments illustrating the principles of the invention will now be discussed with reference to the accompanying figures in which:

[0048] Figure 1A is a flow diagram illustrating in schematic forms methods of the disclosure.

[0049] Figure 1B illustrates schematically an example of an approach to characterise resistance evolution experimentally.

[0050] Figure 2 shows an embodiment of a system for implementing methods of the disclosure.

[0051] Figure 3 shows an overview of the experimental simulation and sampling procedure used to generate synthetic data in examples of the disclosure. Unknown parameters govern the behaviour of cell phenotypes and their response to treatment throughout the simulation. Fixed parameters control key experimental steps, including the initial number of uniquely barcoded cells, their expansion time, sampling bottlenecks during splits into experimental replicates, and the timing of treatment windows. The drug-treatment replicates (DT1-4) and passages (P1-2) shown correspond to those used in the simulation modelling results Figure 4A shows a simple phenotypic compartment model with ‘sensitive’ and ‘resistant’ cells that divide and die with phenotype-specific birth and death rates, respectively. Parameters control the proportion of resistant cells when the simulation begins (p) and the reduction in birth and death rates of resistant cells relative to sensitive cells in the untreated environment (6).

[0052] Fig. 4B shows a hybrid model approach employing a stochastic jump process when nx< Nswitch and a deterministic ODE when nxs Nswitch, for phenotype x.

[0053] Fig. 5A-E illustrate an approach used to model the impact of treatment on population dynamics.

[0054] Fig. 5A shows that to model the effects of treatment, cells’ death rates may be modified according to the effective concentration of the drug at time t: D(t).Fig. 5B shows that the level of drug may be modelled via an uptake / decay model where De controls the maximum effective drug concentration and K controls the rate of uptake / decay experienced by cells. A set of pre-defined treatment on (Ton) and treatment off (Toff) times determine whether the drug is increasing or decreasing. The use of an experimental design with pre-defined treatment on (Ton) and treatment off (Toff) times allows to fit the parameters of the drug uptake / decay model as illustrated on Figs. 5A-5B. However, other experimental designs and models of drug effects are possible within the context of the present disclosure. For example, models that do not model the effect of the treatment using a dynamic model that describes accumulation and / or decay of treatment effect upon exposure and / or withdrawal, respectively, of the treatment can be fitted using experimental data that does not include data acquired using pre-defined treatment on (Ton) and treatment off (Toff) times, such as e.g. only including data at time point before exposure to the treatment (T=0) and at one or more time points during exposure to the treatment.

[0055] Fig. 5C illustrates schematically the components of model A: unidirectional transitions. Phenotype-specific birth and death rates are shown as intrinsic rates without accounting for logistic growth.

[0056] Fig. 5D illustrate schematically the components of model B: bidirectional transitions. Phenotype-specific birth and death rates are shown as intrinsic rates without accounting for logistic growth.

[0057] Fig. 5E illustrates schematically the components of model C: escape transitions. Phenotype-specific birth and death rates are shown as intrinsic rates without accounting for logistic growth.

[0058] Fig. 6A shows a schematic of the sampling scheme adopted to model a long-term resistance evolution experiment using barcoded cells. Following a mutual expansion period, cells are sampled without replacement into 4 drug treatment replicates (DT1, DT2, DT3, DT4). Cell population size measurements are made at two pre-defined observation timepoints (01, 02). Cells are sampled and grown again for a total of two Passages (P1, P2). Drug treatment is controlled via a set of pre-defined treatment on (t e Ton) and treatment off (t G Toff) windows. This figure illustrates a set up adopted in proof-of-concept work described in Examples 1-7. However, any other experimental set up may be used, depending e.g. on the specific models that are fitted and experimental constraints. For example, observation timepoints alone (where observation time points are time points at which cell population size measurements only are obtained) were included before the first passage due to experimental restrictions, with both cell population size and lineage tracing data collected thereafter. Further, the use of replicates enables the calculating of between replicate diversity statistics from the lineage tracing data, which were used in the present examples to evaluate model fit to lineage tracing data. However, the methods described herein can use any possible number of experimental replicates, and any number of passages (time points at which lineage tracing data and cell population size data are obtained), and is not reliant on the availability of observation time points (time points at which cell population size data alone are obtained).

[0059] Fig. 6B shows data that may be generated by this scheme: population size changes at a subset of times (Left), and lineage size distributions at each Passage per replicate (Middle). High dimensional lineage distributions are converted into two-dimensional diversity statistics per replicate Passage (DT, P) (Right).

[0060] Fig. 7A-D show quantitative models of resistance evolution, revealing distinct population dynamics during treatment. Figs. 7A-D (left) show schematics of the different evolutionary models, highlighting key parameters that control the behaviour of resistant, sensitive, and escaped phenotypes. Parameter valuesused for the simulations are shown. Figs. 7A-D (middle) show total cell population size trajectories (dashed lines) and the top 20 lineages (coloured lines) across 4 experimental drug -treatment replicates. Lineage colours are consistent across each set of 4 replicates. Four total cell size 'observations’ are made per replicate: two intermediate population size timepoints (01 and 02) and two passage size timepoints (P1 and P2), which are highlighted. Figs. 7A-D (right) show lineage diversity statistics for the 4 replicate subpopulations (rows) at two passage timepoints (P1 and P2 - columns), including within-replicate lineage diversity (x-axis) and between-replicate diversity (y-axis). Simulation fixed parameters: No= 1x10s, K = 1 x 107, birth rate b = 0.893 (days-1) and death rate d = 0.200 (days-1).

[0061] Fig. 7A shows schematics, trajectories and statistics for model A.

[0062] Fig. 7B shows schematics, trajectories and statistics for model B (unidirectional switching).

[0063] Fig. 7C shows schematics, trajectories and statistics for model B (bidirectional switching).

[0064] Fig. 7D shows schematics, trajectories and statistics for model C.

[0065] Fig. 8A-E illustrate that incorporating lineage data alongside population size changes increases the power to recover parameters governing resistance evolution.

[0066] Fig. 8E shows simulated population size (Left) and cell lineage (Middle) changes across 4 experimental drug-treatment replicates given set treatment windows and their corresponding lineage diversity statistics (Right).

[0067] Fig. 8A shows population size changes through treatment simulated using posterior predictive parameters from an early (top) and late (bottom) generation of the approximate Bayesian computation (ABC) inference first step that fits the total cell size changes at a set of given times to observed replicate timepoints (points). Fig. 8B shows posterior distributions for the resistance model parameters using solely the population data. Fig. 8C shows a subset of the simulated lineage diversity statistics for each sub-population’s two passage timepoints (P1 and P2) used for the parameter inference step, compared to the true diversity statistics for the given simulation (red points). This is highlighted for two of the model parameters: p (top: sensitive-to-resistant phenotype transition probability per cell division) and 5 (bottom: controls the fitness penalty of resistance in the absence of treatment). The adjacent panels (right) show the lineage diversity distance as a function of the two parameters, with the true value highlighted (red dashed line).

[0068] Fig. 8D shows posterior distributions of the inferred parameter values using the combined cell population size and cell lineage statistics (true parameter values shown by red dashed lines). Boxplots show the median, first quartile and third quartile, and whiskers at 1.5x the interquartile range.

[0069] Fig. 9A-C show simulated posterior distributions given different resistance models and evolutionary scenarios. Prior vs posterior distributions given a representative set of simulated synthetic data inferred using both population trajectories and lineage distributions for all three models considered. The following parameters are shown as -Iog10 values: p, a, a.

[0070] Fig. 9A shows results for Model A (uni-directional transitions).

[0071] Fig. 9B shows results for Model B (bi-directional transitions).

[0072] Fig. 9C shows results for Model C (escape transitions).

[0073] Fig. 10A-C show that the modelling framework struggles to recover the true parameter values when resistance is very common and / or the effect of treatment is weak.Fig. 10A shows population size changes through treatment simulated using parameters drawn from the posterior distributions from an early (left) and late (right) generation of the approximate Bayesian computation (ABC) process that fits the population size to observed replicate timepoints (points).

[0074] Fig. 10B shows a subset of the simulated lineage diversity statistics for each subpopulation's two timepoints used for the parameter inference step compared to the true diversity statistics for the given simulation (red points), highlighted for two of the model parameters: p (top: sensitive-to-resistant phenotype transition probability per cell division) and 5 (bottom: controls the fitness penalty of resistance in the absence of treatment). The adjacent panels (right) show the lineage diversity distance as a function of the two parameters, with the true value highlighted (red dashed line).

[0075] Fig. 10C shows posterior distributions of the inferred parameter values using the combined cell population size and cell lineage statistics (true parameter values shown by red dashed lines).

[0076] Fig. 11A-D shows that single re-barcoding can improve the recovery of true parameter values under certain conditions.

[0077] Fig. 11 A shows a schematic illustrating how cells were barcoded with a second lineage barcode tag at the time of P1. The counts of both barcodes were recorded at the end of the simulation. Within-replicate diversity and between-replicate diversity were calculated for barcode i, whereas only within-replicate diversity was calculated for barcode j.

[0078] Fig. 11B shows the lineage diversity statistics for each replicate's Passage (P1, P2) for a simulation using Model B (bidirectional transitions) calculated using the first barcode, comparing simulated to observed values as a function of the reversion transition parameter (a). The distances in diversity space are summarised in the middle column and the posterior estimates for all Model B parameters are shown on the right with red dashed lines indicating the true values.

[0079] Fig. 11C shows lineage diversity statistics as in Fig. 11B, but only showing the within-replicate diversity (qD) calculated using the second barcode (j).

[0080] Fig. 11 D shows lineage diversity statistics as in Fig. 11B, but for a simulation using Model C (escape transitions).

[0081] Fig. 12A-B illustrate model selection with DIC scores. Bar charts show the number of times the model selection step chose either of the two models, given the true underlying model was unidirectional switching or escape transitions is shown. Schematics and parameter values are shown for each model. A subset of diversity statistics are shown, given the posterior estimates of each model for the unidirectional model (top panels) and escape model (bottom panels). For each of the example outputs shown, the lower DIC score supports the correct model in each case.

[0082] Fig. 12A shows results when the true underlying model was unidirectional switching.

[0083] Fig. 12B shows results when the true underlying model was escape transitions.

[0084] Fig. 13 shows a schematic illustrating a stochastic agent-based lineage model that tracks the birth and death of cells with phenotypes and heritable cell lineage tags.

[0085] Fig. 14A-D illustrate a parameter inference framework.

[0086] Fig. 14A illustrates that simulations are performed given a sample from the parameters' prior distributions. Population sizes are recorded at the observation (O) and passage (P) times.Fig. 14B illustrates that the Euclidean distance is calculated between the simulated and observed datasets at each of these positions using ABC-SMC implemented in pyABC.

[0087] Fig. 14C illustrates for the final iteration of the ABC, the use of a stochastic agent-based model to simulate the lineage distributions and conversion of lineage distributions into diversity statistics per replicate Passage (DTx, Pn).

[0088] Fig. 14D shows calculation of the distance between the simulated and observed diversity statistics, only retaining the closest values to approximate the parameters’ posterior distribution.

[0089] Fig. 15 shows a comparison of posterior estimates for the pre-existing resistance fraction (p - top row) and sensitive-to-resistant transition probability (p - bottom row) parameters for a range of Nswitch values used in the hybrid phenotypic compartment model from Model A. Fits were performed on synthetic data where the true parameter values were known (red dashed lines). Posterior distributions are shown based only on the first step of the inference framework, prior to incorporating distance metrics from the cell lineages.

[0090] Fig. 16 shows the impact of differing diversity order q on parameter inference. Simulation output 2D diversity statistics are shown compared to the ground truth (red points) (LHS). This distance is shown as a function of the parameter p. The posterior distributions for all 6 parameters in model A are also shown (right).

[0091] Fig. 17A-F show quantification of the evolutionary dynamics during treatment in a barcoded colorectal cancer cell line (SW6bc).

[0092] Fig. 17A shows a simplified schematic of the experimental design for 4 drug-treatment replicates in a longterm evolutionary resistance assay. Cells were barcoded and expanded (POT) before being split into 4 replicate drug-treatment sub-populations that were exposed to periodic chemotherapy treatment (5-fluoruracil: 5-Fu) and sampled at two timepoints (Passage 1: P1 and Passage 2: P2). Other conditions including vehicle control treatment and later passages (P3-P4) are not shown.

[0093] Fig. 17B shows population size observations at four timepoints per-replicate (DT1-4), including early drugstop cell counts (not shown in Fig. 17A), image estimates and two passage timepoint cell counts (left panel); top 10 sequenced barcode lineage frequencies at the two passage timepoints, where area colour denotes lineage identity (middle panel); and lineage diversity statistics of the sequenced barcode distributions at the two passage timepoints (P1 and P2) for each replicate. Within replicate lineage diversity is shown on the x-axis) and between-replicate diversity dissimilarity is shown on the y-axis, with the range dictated by the number of sub-populations: [1,4] (right panel).

[0094] Fig. 17C shows posterior predictive distributions for Model A (unidirectional transitions - top) and Model C (escape transitions - bottom) of the within-replicate diversity (x-axis) and between-replicate diversity (y-axis) statistics (left) and the normalised diversity distance from the observed statistics (red points) to the simulated values (highlighted by the transition parameter p) (right).

[0095] Fig. 17D shows the Deviance Information Criterion (DIC), a measure of model fit, for Model A and Model C (lower values indicate higher model support) calculated using the posterior predictive distributions. Fig. 17E shows the posterior distribution for parameters in Model A (the model with the highest support). Boxplots show the median, the first and third quartiles, and whiskers at 1 5x the interquartile range.

[0096] Fig. 17F shows the data produced from a single simulation from the posterior predictive distribution using Model A, as in Fig. 17B.Fig. 18A-F show quantification of the evolutionary dynamics during treatment in a barcoded colorectal cancer cell line (HCTbc).

[0097] Fig. 18A shows a simplified schematic of the experimental design for 4 drug-treatment replicates in a longterm evolutionary resistance assay. Cells were barcoded and expanded (POT) before being split into 4 replicate drug-treatment sub-populations that were exposed to periodic chemotherapy treatment (5-fluoruracil: 5-Fu) and sampled at two timepoints (Passage 1: P1 and Passage 2: P2). Other conditions including vehicle control treatment and later passages (P3-P4) are not shown.

[0098] Fig. 18B shows population size observations at four timepoints per-replicate (DT1-4), including early drugstop cell counts (not shown in Fig. 18A), image estimates and two passage timepoint cell counts (left panel); top 10 sequenced barcode lineage frequencies at the two passage timepoints, where area colour denotes lineage identity (middle panel); and lineage diversity statistics of the sequenced barcode distributions at the two passage timepoints (P1 and P2) for each replicate. Within-replicate lineage diversity is shown on the x-axis) and between-replicate diversity dissimilarity is shown on the y-axis, with the range dictated by the number of sub-populations: [1,4].

[0099] Fig. 18C shows posterior predictive distributions for Model A (unidirectional transitions - top) and Model C (escape transitions - bottom) of the within-replicate diversity (x-axis) and between-replicate diversity (y-axis) statistics (left) and the normalised diversity distance from the observed statistics (red points) to the simulated values (highlighted by the transition parameter p) (right).

[0100] Fig. 18D shows the Deviance Information Criterion (DIC), a measure of model fit, for Model A and Model C (lower values indicate higher model support) calculated using the posterior predictive distributions. Fig. 18E shows the posterior distribution for parameters in Model C (the model with the highest support). Boxplots show the median, the first and third quartiles, and whiskers at 1.5x the interquartile range.

[0101] Fig. 18F shows the data produced from a single simulation from the posterior predictive distribution using Model C, as in Fig. 18B.

[0102] Fig. 19 shows the majority of cells contain a single barcode which adheres to the expected nucleotide sequence. The relative frequency of sequenced barcodes in expanded single cell colonies from each barcoded cell line (top: HCTbc, bottom: SW6bc) is shown, highlighted by whether or not the sequenced barcode adhered to the expected ClonTracer weak-strong (WS) 30bp nucleotide pattern. Colony 1 corresponds to a control well where 100 cells were seeded and expanded before sequencing. The lack of a dominant barcode and higher number of lower frequency barcodes that adhere to the WS pattern support this condition.

[0103] Fig. 20A-B shows long-term resistance evolution experiment dose-response curves. Each panel shows the response for the parental (POT) and a given experimental condition (CO: control, DS: drug-stop, DT: drug treatment). Adjacent panels show the IC50 values for all conditions.

[0104] Fig 20A shows dose-response curves for all experimental conditions from the P4 samples for HCTbc. Fig. 20B shows dose-response curves for all experimental conditions from the P4 samples for SW6bc.

[0105] Fig. 21 shows a schematic illustrating the sampling design of the expanded population of barcoded cells used to infer the population birth and death rates. Barcoded cells are expanded for a known At and thensampled (without replacement) K times for PCR amplification of the barcodes and sampled again (with replacement) J times given the allocation of reads to that sample when sequencing. The conditional probabilities correspond to the probability of seeing an individual lineage at each stage, with the unknown parameters highlighted in red.

[0106] Fig. 22 shows posterior distributions of the birth and death rates inferred from synthetic data over a range of combinations (denoted by simulation number) where the true values are known

[0107] (red points).

[0108] Fig. 23A shows posterior distributions of the birth and death rates inferred using sequenced lineage count data from two barcoded colorectal cancer cell lines, SW6bc (MSS) and HCTbc (MSI). Inference was performed on the combined lineage counts from three parental (‘POT’) replicates per cell line.

[0109] Fig. 23B shows posterior predictive distributions of lineage counts given 100 draws from each cell line’s birth and death rate posterior distributions.

[0110] Fig. 24A-J show the functional characterisation of resistance in SW6bc.

[0111] Fig.24A shows a schematic of sample design and naming schemes from the long-term resistance evolution experiment using barcoded SW620 colorectal cancer cells. Periods of on / off treatment (5-Fu and DMSO) are for illustrative purposes and are not to scale.

[0112] Fig. 24B shows UMAP dimensionality reduction of the scRNA-seq data from chosen sample replicates. Fig. 24C shows cluster assignment of UMAP results on a sample-by-sample basis highlighted by cluster identity.

[0113] Fig. 24D shows cluster assignment proportions by cell sample.

[0114] Fig. 24E shows cell cycle stage proportions by cell sample.

[0115] Fig. 24F shows genes associated with replicate-specific differential expression (DE) identified via specific comparisons (coloured points per panel): Resistance Treatment Response - DE in drug-treatment (DT) replicates but not drug-stop (DS) replicates; Acute Treatment Response - DE in DT replicates and DS replicates; Sensitive Treatment Response - DE in DS replicates but not DT replicates. Comparisons are repeated for each drug treatment replicate (DT1 and DT3) and the number of up / down DE genes in each group are also shown per comparison.

[0116] Fig. 24G shows gene set enrichment analysis (GSEA) for each drug-treatment replicate (DT1 and DT3) for genes identified as ‘Resistance’ genes in each replicate. Only gene sets that were found to be significant (adjusted p-value < 0.05) in one of the two replicates are shown. Point colour denotes normalised enrichment score (NES).

[0117] Fig. 24H shows the number of shared / unique resistance genes between each drug-treatment replicate. Fig. 24I shows ‘Resistance’ genes GSEA NES comparison between DT1 and DT3.

[0118] Fig. 24J shows a subset of single-cell copy number profiles from SW6bc experimental replicates. Each row is a single cell and chromosome bin colour denotes copy number. Cluster assignment and sample identity are also shown.

[0119] Fig. 25A-B show consensus copy number profiles from pooled scWGS of SW6bc experiment replicates. Fig. 25A shows consensus copy numbers for all chromosomes for the six condition replicates.Fig. 25B shows regions on chromosome 17 and 18 for the same condition replicates, highlighting an additional shared chromosomal gain unique to the two drug-treatment replicates.

[0120] Fig. 26A-J show functional characterisation of resistance in HCTbc.

[0121] Fig.26A shows a schematic of sample design and naming schemes from the long-term resistance evolution experiment using barcoded HCT 116 colorectal cancer cells. Periods of on / off treatment (5-Fu and DMSO) are for illustrative purposes and are not to scale.

[0122] Fig. 26B shows UMAP dimensionality reduction of the scRNA-seq data from chosen sample replicates. Fig. 26C shows cluster assignment of UMAP results on a sample-by-sample basis highlighted by cluster identity.

[0123] Fig. 26D shows cluster assignment proportions by cell sample.

[0124] Fig. 26E shows cell cycle stage proportions by cell sample.

[0125] Fig. 26F shows genes associated with replicate-specific differential expression (DE) identified via specific comparisons (coloured points per panel), as in Fig. 24F.

[0126] Fig. 26G shows gene set enrichment analysis (GSEA) for each drug-treatment replicate (DT3 and DT4) for genes identified as 'Resistance' genes in each replicate. Only gene sets that were found to be significant (adjusted p-value < 0.05) in one of the two replicates are shown. Point colour denotes normalised enrichment score (NES).

[0127] Fig. 26H shows the number of shared / unique resistance genes between each drug-treatment replicate. Fig. 26I shows ‘Resistance’ genes GSEA NES comparison between DT3 and DT4.

[0128] Fig. 26J shows mutations found in scRNA-seq data grouped by the number of replicates any mutation was found in and VAF distributions of mutations found in scRNA-seq data unique to the two drug-treatment replicates.

[0129] Fig. 27A-D show single cell fitness assays in HCTbc.

[0130] Fig. 27A shows a schematic of single cell isolation and clonal expansion, performed for two control (CO) and two drug-treatment (DT) replicates for the HCTbc cell line.

[0131] Fig.27B shows a snapshot of the Incucyte data used to generate the growth rates from confluence readings over time.

[0132] Fig. 27C shows log-confluence estimates and linear regression fits (coloured lines) for 12 chosen colonies per replicate.

[0133] Fig. 27D shows growth rates as estimated by the gradient fits of the log-confluence values per replicate sample. Pairwise comparisons were conducted using the Wilcoxon Rank Sum test. Significance levels are annotated as follows: ns p > 0.05, * p<=0.05, ** p<=0.01, *** p<=0.001, **** p <= 0.0001.

[0134] Fig. 28A-D show that sorting by size reveals distinct transcriptional types in HCTbc cells that survive multiple rounds of treatment.

[0135] Fig. 28A shows a comparison of brightfield images (top panels) of the large (left-hand panels) and small (right-hand panels) barcoded HCTbc cell phenotypes at the time of single-cell sorting following seven weeks of periodic chemotherapy (5-Fu) treatment. Illustrative single cell sorting images are also shown (bottom panels).Fig. 28B shows genes associated with replicate-specific differential expression (DE) identified via specific comparisons (coloured points per panel): Large / Small Treatment Response - DE in Large / Small cells, but not DS; Acute Treatment Response - DE in Large / Small cells and DS; Sensitive Treatment Response - DE in DS but not Large / Small cells. Comparisons are repeated for each size-sorted phenotype (Large - left panel, Small - right panel) and the number of up / down DE genes in each group are also shown per comparison.

[0136] Fig. 28C shows correlation coefficients between each replicate given the average expression across all genes and all cells.

[0137] Fig. 28D shows gene set enrichment analysis (GSEA) for each size-sorted phenotype (Large - left column, Small - right column). Only gene sets that were found to be significant (adjusted p-value < 0.05) in one of the two phenotypes are shown. Point colour denotes normalised enrichment score (NES).

[0138] Fig. 28E shows differences in NES from the genes analysed in Fig. 28D.

[0139] Fig. 29A shows an experimental design for testing a plasticity modulator. Two parallel experiments are shown: one with only cytotoxic treatment following the expansion step, and the other with the addition of a plasticity modulator throughout the experiment.

[0140] Fig. 29B shows results when a single phenotypic transition parameter was used to fit both experimental arms. Middle column: Comparison of within- and between-replicate diversity statistics between model fit and observed data (red points). Right column: Posterior estimate of the single transition rate p. Red dashed line indicates the true value.

[0141] Fig. 29C shows results as in Fig 29B, but for the case when separate transition rates were used to fit each experiment. Diversity statistics shown are for 'Experiment T in Fig. 29B and Fig. 29C.

[0142] Fig.30A-B illustrate the impact of over-estimating the barcoding efficiency on accurate parameter recovery.

[0143] Fig.30A shows a case where the estimated number of uniquely barcoded cells when the experiment begins (No) is accurate, producing simulations with similar diversity statistics to the observed data (LHS) and recovering the true parameter estimates.

[0144] Fig. 30B shows a case where the estimated number of barcoded cells (No = 105) is a significant overestimate, leading to incorrect parameter estimates (RHS).

[0145] Fig. 31 shows a conceptual model of cross-resistance linking molecular features, resistance mechanisms and cell phenotypes. Whether a cell exhibits a given phenotype depends on which mechanisms it expresses, which in turn determines the relative fitness of each phenotype in the two drug environments. Drug environments are defined as single concentrations of each (X and Y for Drug A and B, respectively). Sensitive cells are defined as those containing the ‘wild-type’ (WT) molecular feature, i.e. the absence of the feature(s) of interest.

[0146] Fig. 31A shows a schematic showing a single molecular feature controlling a single cross-resistance mechanism.

[0147] Fig. 31 B shows a schematic showing two independent molecular features controlling resistance mechanisms to two drugs.

[0148] Fig 31 C shows a schematic showing a group of co-inherited molecular features forming a single mechanism that confers cross-resistance.Fig. 32 shows the phenotype-specific transitions and drug-induced cell death in a two drug model of resistance. The possible transitions between phenotypic compartments and phenotype-specific birth and death rates (lefthand plot) and the elevated death rates of phenotypes in the two drug environments (middle and righthand plots) are shown. This model structure assumes there is no cross-resistance, and for simplicity assumes that the resistant phenotypes are ‘complete’ (ipA= 0.0 and

[0149]

[0150] = 0.0).

[0151] Fig. 33 shows a schematic representation of a model of cross-resistance between two drugs.

[0152] Fig. 33A shows an illustration of how varying the strength of cross-resistance (0) alters the statistical association between resistance to Drug A and Drug B. The area occupied by each coloured region represents the fraction of cells pre-existing in each phenotype: sensitive (S, blue), resistant to Drug A only (RA, red), resistant to Drug B only (RB, orange), or resistant to both (RAB, green). The total pre-existing fractions resistant to each drug, pxand pB, remain constant across panels; what changes with increasing 0 is the probability that these resistance traits co-occur within the same cells.

[0153] Fig. 33B shows corresponding schematics of the phenotypic transition structure. Arrow weights denote the relative probabilities of each transition. With no cross-resistance (0 = 0), independent mutation paths lead from the sensitive phenotype (S) to each single-resistant phenotype (RAor RB), and subsequently to double resistance (RAB). As 0 increases, transitions that generate RABdirectly from S become more likely, culminating in complete cross-resistance where a single mechanism confers resistance to both drugs simultaneously. For simplicity, the schematics assume symmetry in the two drugs’ pre-existing fractions and transition probabilities.

[0154] Fig. 34 shows an example cross-resistance inference workflow.

[0155] Fig. 34A shows an overview of the experimental design and the experimental measurements used for parameter inference (population sizes at passage events and diversity statistics derived from barcode counts). The experimental design measures the diversity of related cell lineages exposed to two separate drugs using 'DNA barcoding'. Lineage information can be converted into diversity statistics that capture the diversity within and between replicates, and within and between drug conditions.

[0156] Fig. 34B shows an inference workflow for the single-drug parameters (illustrated for Drug A). For clarity, the initial population-only ABC-SMC step (Stage 1a) is omitted, and only Stage 1b (training the neural density estimator using simulated diversity statistics) is shown. A neural density estimator is trained using data simulated from each parameter's prior distribution for each drug condition independently. The observed, experimental population size changes and lineage diversity information are then used to approximate the posterior distribution for parameters that capture the dynamics of resistance to a single drug using this optimized neural network.

[0157] Fig. 34C, D, E show an inference workflow for the cross-resistance strength 0. Example barcode distributions for each drug illustrate how 0 influences the between-drug lineage relationships used in the inference. The posterior for each single drug's resistant phenotype and the prior distribution of the strength of cross resistance (0), are used to jointly simulate the experiment across both drugs and train a second neural density estimator. The observed within and between drug diversity statistics are then used to approximate the posterior of the strength of cross resistance.

[0158] Fig. 35 shows example simulations for varying strengths of cross resistance in a two-drug model of resistance. The figure shows that the change in the phenotypic composition of the population when exposedto either drug in the scheme illustrated on Figs. 32-33 is a function of the strength of cross resistance (0), which also influences the relationship between the cell lineages that survive treatment.

[0159] Fig. 36 shows simulation results illustrating that the strength of cross resistance (0) determines the probability of the same cell lineages being enriched across both drug condition's experiments, which in turn determines the between-drug diversity statistics.

[0160] Fig. 37A-B show data from simulations demonstrating that the inventors’ cross resistance modelling framework (see Example 8) can recover the true value of cross-resistance strength (0 - red points / dashed line) from synthetic data generated under a range of evolutionary scenarios. On both figures, each row of plots shows for a particular evolutionary scenario: within and between-drug diversity statistics used to train the neural network and the ground truth value highlighted with a red dot (and black arrow); diversity distance as a function of theta (cross-resistance) for the synthetic data, with ground truth highlighted with a dashed vertical line; and strength of cross-resistance (0) posterior and ground-truth value (red dashed line) approximated using the optimised neural network.

[0161] Fig. 37A shows results for scenarios with a rare pre-existing resistance with a moderate S-> R transition probability. In each scenario, the drug’s resistant phenotypes have the same single-drug parameters, and the three rows compare scenarios with no, moderate, and complete cross resistance (from bottom to top).

[0162] Fig. 37B shows results for scenarios with a very high pre-existing resistance with low S-> R transition probability. In each scenario, the drug’s resistant phenotypes have the same single-drug parameters, and the three rows compare scenarios with no, moderate, and complete cross resistance (from bottom to top).

[0163] Fig. 38 shows a schematic illustrating an experimental design used for implementation of the crossresistance inference framework for three breast cancer cell lines, exposed to long-term treatment with two CDK4 / 6 inhibitors (Abemaciclib and Palbociclib).

[0164] Fig. 39A shows the experimental data obtained using the set up in Fig. 38, in terms of sequenced barcode distributions.

[0165] Fig.39B shows the strength of cross resistance estimates (densities) for each of the cell lines. The analysis predicts differing levels of cross-resistance across the three cell lines.

[0166] Fig. 39C shows the normalised time gained via treatment switching from Palbociclib to Abemaciclib compared to continued treatment with Palbociclib using the posterior distributions for the single drug parameters and strength of cross resistance estimate (0) for each of the ER+ breast cancer cell models (BT47, MCF7, T47D).

[0167] Fig. 39D shows the single drug posterior estimates for each of the drugs and cell lines.

[0168] Fig.39E schematically illustrates the process of estimating the benefit of treatment switching used to obtain the results in Fig. 39C. The phenotypes shown are purely illustrative.

[0169] Fig 39F shows the experimentally measured survival fraction of cells exposed to different concentrations of Abemaciclib or Palbociclib for each of the three ER+ breast cancer cell models. The three samples shown include the sensitive, ancestral barcoded population (POT) and the evolved palbociclib resistant line (Palbo.res). The concentration used during the long-term evolution experiment to generate resistance is highlighted (black dashed line).Detailed Description

[0170] Aspects and embodiments of the present invention will now be discussed with reference to the accompanying figures. Further aspects and embodiments will be apparent to those skilled in the art. All documents mentioned in this text are incorporated herein by reference.

[0171] The evolution of drug resistance remains the primary reason for treatment failure in cancer patients. Indeed, the evolution of resistance in cancer hampers effective treatment and limits the lifespan of newly developed drugs. Despite this, resistance evolution is not normally addressed until oncology drugs have been used to treat patients. Current pre-clinical assays (such as e.g. assays that assess cell viability across a range of drug concentrations) focus on the short-term response of cancer cells to treatment, and only look for the average response of a population of cells. These approaches fail to identify resistant sub-populations and therefore ignore the longer-term issue of resistance evolution. The present inventors have identified that there was a lack of pre-clinical assays, and specifically assays that can be performed in vitro, that characterise the evolution of resistance to treatments. They have developed a new platform, implemented in embodiments described herein as an in vitro platform, which they termed “Evolutionary-informed resistance assays” (EIRAs). The technology can provide a quantitative, pre-clinical readout of resistance evolution to a given drug. The data provided by the assay can enable the robust comparison of different candidate compounds, with a focus on identifying the sub-populations driving resistance to the treatment and measuring the ‘ease’ with which cells evolve resistance to either drug. This information is extremely value in the process of drug development.

[0172] The approach measures the behaviour of cells as resistance emerges over a longer period than traditional assays. Lineage tracing techniques (e.g. using lineage tracing tools in vitro, or naturally occurring somatic mutations) are used to quantify sub-clonal dynamics during treatment. Computational models simulate the cell population dynamics under different scenarios of how resistance can evolve, these simulations are fitted to measurements of cell population size changes through treatment and sequenced lineage distributions, finding the scenarios that best fit the observed data. The methods can predict whether resistance is the product of a stable genetic change or a transient non-genetic mechanism. They can also provide statistical readouts of the identity and phenotypic behaviours of the sub-populations driving resistance (for example, mutation rates from sensitive to resistant). Predictions generated with the models can be used to determine which approach to use for subsequent functional molecular characterisation of resistance (e.g. whether resistance is likely to be genetically caused or not). Further, the predictions can also be used to guide decision making at critical junctures of the drug development process such as: target discovery, triage of candidate compounds, treatment combination identification, and stratification of patients for clinical trials. For example, at the target discovery stage, the methods can be used to identify mechanisms to current treatments that provide novel targets for the generation of new targeted drugs. At the triage step, the methods can be used e.g. to identify which compounds in development minimise the propensity for resistance evolution. In the context of identifying treatment combinations, the methods can be used to test possible combination treatments to identify those that limit resistance evolution. In the context of patient stratification, the methods can be used to identify biomarkers that predict resistance, to identify sub-populations of responders for clinical trials (or indeed treatment selection).In Whiting 2022, the present inventors described models that explicitly capture both genetic and non-genetic sources of phenotypic variation in cell populations evolving resistance to therapy. While the approaches described therein bear some resemblance to the methods of embodiments of the present disclosure, there are also some important differences. Firstly, the present disclosure makes use of both total population size data and lineage tracing data to fit different models of resistance evolution. By contrast the work in Whiting 2022 only used lineage tracing data. The present inventors have realised that by using total population size data, it was possible to fit compartment models that track the dynamics of this data, and also corresponding models that track individual lineages and that model the same resistance evolution phenomena (i.e. sharing the same parameters representing the same phenomena), such that the former can be used to perform a first “rough” but computationally efficient fitting of the parameters of the models, while the latter can be used to more precisely and accurately identify model parameters. This not only improves the computational efficiency of the process of quantifying these key resistance parameters, it also results in fitted parameters that are more accurate and therefore more reliably quantify resistance evolution. Further, for the purpose of the first fitting step using total cell population size data, the inventors designed a new hybrid phenotypic compartment approach to simulate the population size of each of a plurality of phenotypes (i.e.subpopulations of cells that respond differently to treatment), treating the population as a stochastic jump process when below a predefined threshold and switching to a deterministic approximation when the population size reaches or exceeds this threshold. This advantageously balances accuracy of simulations at small population numbers, and computational efficiency. Further, for the purpose of the second fitting step, the inventors constructed agent-based models that preserves all the features of the first model. This second approach is more computationally expensive than the hybrid phenotypic compartment model, but also tracks the lineage identities of individual cells within the system. As such, this agent-based version can also generate cell lineage distributions. In summary, the use of a two-step method as proposed makes the model fitting process computationally tractable by leveraging first the total population size data then the lineage tracing data. This is, to the best of the inventor’s knowledge, a completely novel approach. Specific embodiments use in the first step a hybrid phenotypic compartment model (which is itself novel to the best of the inventors' knowledge) to exclude parameter combinations that could not account for observed changes in cell population sizes during treatment by comparing vectors of population size changes at predetermined times (07, 02, P1, P2). Specific embodiments use in the second step an agentbased lineage model to simulate the remaining non-identifiable parameter combinations and generate cell lineage distributions at each passage (P1, P2). Additionally, the inventors further designed a new approach to estimate posterior distributions of model parameters in the second fitting step that uses lineage tracing data, by retaining parameter combinations that produce lineage distributions with multiple diversity summary statistics (specifically, a within-replicate diversity and a between-replicate diversity, and in specific embodiments modelling cross-resistance, between drug-diversity) that closely match those derived from observed lineage tracing data. Specifically, the inventors designed a Bayesian approach to identify these parameters, resulting in a tractable but accurate and statistically grounded identification of important parameters that characterise evolution of resistance. By contrast, the work in Whiting 2022 only used the lineage tracing data to try and identify best parameter combinations in a non-Bayesian rudimentary ‘grid search’ approach. Finally, specific embodiments model cross-resistance between two drugs, in which anunderlying mechanism of resistance confers resistance of a population of cells to two drugs (i.e. providing a survival advantage in the presence of either drug), the strength of that cross-resistance being a quantifiable parameter of the model.

[0173] Methods of the present disclosure relate to the evolution of resistance in proliferative cells exposed to antiproliferative treatment. The proliferative cells may be selected from an immortalised cell line, a cancer cell line, or primary cancer cells. A cancer cell line may be a known pre-clinical model cell line. The cells may be cultured in vitro, or may be cells sampled from a tumour. The cells may be cultured in a cell culture environment as isolated cells, tissue or multicellular structures (e.g. organoids or spheroids). The cell culture environment may include a 2D or 3D cell culture environment. The term “cell culture” is used generally to refer to an environment in which cells (either isolated cells, tissues or multicellular cultures) can be maintained alive in an in vitro environment.

[0174] An anti-proliferative treatment refers to any compound, composition or condition (including e.g. radiation) that can have a cytotoxic or cytostatic effect on cells, including but not limited to proliferative cells. In embodiments, the anti-proliferative treatment is a cytotoxic treatment. An anti-proliferative treatment may comprise exposure to one or more chemicals (e.g. small molecules), one or more biologies (e.g. one or more active proteins, peptides, nucleic acids and combinations thereof), one or more genetic perturbations (e.g. knock out or knock down of one or more genes, using any technology known in the art such as e.g. CRISPR based KO, siRNA, etc.), one or more cell-based treatments (e.g. CAR-T, TCR-T, etc.), or any combination thereof. While an anti-proliferative treatment may generally be referred to herein as “drug”, this does not imply that only a single active compound, agent or perturbation is used, and instead a consistent combination of e.g. multiple active compounds may be referred to together as a “drug”. In some embodiments, a drug treatment (or anti-proliferative treatment) may refer to a treatment with a single active compound.

[0175] The term “lineage tracing data” refers to data that characterises the clonal composition of a cell population. The data may be obtained using a genetic lineage tracing system, also referred to as “cell barcodes or simply “barcode”) such as ClonTracer (Bhang et al. 2015), or any other experimental method that permits the tracking of cell lineages (i.e. tracking of individual clones). A genetic lineage tracing system can comprise a library of vectors (e.g. viral vectors) comprising barcodes that can be inserted into cells transduced with the vectors, such that each individual cell in the transduced population may be expected to include at least one unique barcode. Lineage tracing data may be obtained as sequencing data (e.g. a plurality of sequencing reads from which barcodes can be identified and counted or barcode counts derived therefrom, or a plurality of sequencing reads from which a clonal composition can be identified using naturally acquired somatic mutations in the cell population).

[0176] Methods of the present disclosure make use of dynamic population growth models that represent the growth of a plurality of subpopulations of cells having different responses to the anti-proliferative treatment comprising at least a first subpopulation of cells that is sensitive to the anti-proliferative treatment and a second population of cells that is resistant to the anti-proliferative treatment. The models may comprise one or more of: a deterministic model such as a population level differential equation based model, and a stochastic model such as an agent based model or stochastic branching process model. An agent basedmodel as used herein, refers to a stochastic model that tracks each individual cell in the population. Such a model can model lineage tracing information as each cell can be associated with a lineage identity. A stochastic branching process refers to a stochastic model that does not assign individual lineage identities to individual cells, where cells in each subpopulation collectively and stochastically undergo transitions and events with probabilities that are fitted parameters of the model. In embodiments, a pair of anti-proliferative treatments are analysed jointly to identify mechanisms of resistance to each of the therapies individually as well as mechanisms of resistance that confer cross-resistance (i.e. at least some survival advantage under treatment with either of the two anti-proliferative treatments).

[0177] The cell subpopulations may also be referred to as "compartments” or “phenotypes”. Cells in each subpopulation have a different response to the anti-proliferative treatment than cells in another subpopulation. The different responses can be modelled in terms of differences in birth and / or death rates in the presence and / or in the absence of treatment. Thus, the methods described herein can model the dynamics of a plurality of cell subpopulations that differ from each other by one or more of: their birth rate in the absence of treatment, their birth rate in the presence of treatment, their death rate in the absence of treatment, and their death rate in the presence of treatment. The plurality of subpopulations can be associated with transitions between specific subpopulations (also referred to as “allowed transitions”), where each allowed transition can be associated with a respective rate. Transitions that are reversible can be modelled as a pair of transitions (one in each direction), where each transition in the pair can be associated with a respective rate that can be the same or different. In embodiments, the plurality of cell subpopulations comprises a subpopulation of sensitive cells and a subpopulation of resistant cells. Sensitive cells are cells that are sensitive to the anti-proliferative treatment. For example, their death rate may be higher in the presence of the treatment than in its absence, and / or their birth rate may be lower in the presence of the treatment than in the absence of the treatment. Resistant cells are cells that are at least to some extent resistant to the treatment, in that at least one effect of the treatment is lower in these cells than in sensitive cells. For example, their death rate in the presence of the treatment may be lower than the death rate of sensitive cells in the presence of the treatment. Resistant cells may incur a fitness penalty in the absence of treatment. For example, resistant cells may have a lower birth rate than sensitive cells in the absence of the treatment. This is optional and in some embodiments the resistant cells and the sensitive cells may have the same birth and / or death rates in the absence of the treatment. In embodiments, the plurality of cell subpopulations comprises: (i) a subpopulation of sensitive cells, and (ii) one or more subpopulations of resistant cells and / or one or more subpopulations of escaped cells. Escaped cells may also be represented in embodiments of the disclosure. Escaped cells are resistant cells that do not incur a fitness penalty in the absence of treatment, or that have a lower penalty in the absence of treatment than other resistance cells. In embodiments, such as e.g. embodiments in which a single anti-proliferative treatment is under investigation or embodiments in which a pair of anti-proliferative treatments are under investigation but single anti-proliferative treatment parameters are being fitted, the plurality of cell subpopulations may comprise: (i) a single subpopulation of sensitive cells, (ii) a single subpopulation of resistant cells, and (iii) optionally, a single subpopulation of escaped cells. In embodiments, the plurality of cell subpopulations consists of between 2 and 4 subpopulations. Without wishing to be bound by theory, the present inventors have found through experimentation on different datasets that 2-4 types are oftenenough to capture observed population and lineage dynamics. In some embodiments in which a pair of anti-proliferative treatments (comprising a first and second anti-proliferative treatments) is under investigation, the plurality of cell subpopulations may comprise: (i) a single subpopulation of sensitive cells, (ii) a first subpopulation of resistant cells that are resistant to the first anti-proliferative treatment only, (iii) a second subpopulation of resistant cells that are resistant to the second anti-proliferative treatment only; and (iv) a third subpopulation of resistant cells that are resistant to the first anti-proliferative treatment and the second anti-proliferative treatment. Cells in the first, second and third subpopulations of resistant cells may be assumed to possess different mechanisms of resistance. For example, cells in the first subpopulation of resistant cells may be assumed to possess a first resistance mechanism conferring resistance to the first anti-proliferative treatment. Cells in the second subpopulation of resistant cells may be assumed to possess a second resistance mechanism conferring resistance to the second anti-proliferative treatment. Cells in the third subpopulation of resistant cells may be assumed to possess one or both of: (i) both the first and second resistance mechanisms, and (ii) a third resistance mechanism conferring resistance to both antiproliferative treatments. A fitted parameter (denoted as 0 in the examples of the disclosure) may moderate the rate of transition of sensitive cells to the third population of resistant cells, which quantifies the strength of departure of this rate of transition from a model in which resistance to both anti-proliferative treatments arises through independent acquisition (i.e. resistance to the first anti-proliferative treatment is acquired independently from resistance to the second anti-proliferative treatment) rather than joint acquisition (i.e. resistance to the first anti-proliferative treatment is not acquired independently from resistance to the second anti-proliferative treatment). In other embodiments in which a pair of anti-proliferative treatments (comprising a first and second anti-proliferative treatments) is under investigation, the plurality of cell subpopulations may comprise: (i) a single subpopulation of sensitive cells, and (ii) a single subpopulation of resistant cells that are resistant to the first anti-proliferative treatment with a first resistance parameter moderating the death rate of the population of resistant cells in the presence of the first anti-proliferative treatment and a second resistance parameter moderating the death rate of the population of resistant cells in the presence of the second anti-proliferative treatment.

[0178] The systems and methods described herein can be implemented in a computer system, in addition to the structural components and user interactions described. As used herein, the term “computer system” includes the hardware, software and data storage devices for embodying a system and carrying out a method according to the described embodiments. For example, a computer system can comprise one or more central processing units (CPU) and / or graphics processing units (GPU), input means, output means and data storage, which can be embodied as one or more connected computing devices. Preferably the computer system has a display or comprises a computing device that has a display to provide a visual output display. The data storage can comprise RAM, disk drives, solid-state disks or other computer readable media. The computer system can comprise a plurality of computing devices connected by a network and able to communicate with each other over that network. It is explicitly envisaged that computer system can consist of or comprise a cloud computer.

[0179] The methods described herein are computer implemented unless context indicates otherwise. Indeed, the features of lineage tracing data and the complexity of the simulations necessary to fit models as describedherein are such that the methods described herein are far beyond the capability of the human brain and cannot be performed as a mental act. The methods described herein can be provided as computer programs or as computer program products or computer readable media carrying a computer program which is arranged, when run on a computer, to perform the method(s) described herein. As used herein, the term ‘‘computer readable media” includes, without limitation, any non-transitory medium or media which can be read and accessed directly by a computer or computer system. The media can include, but are not limited to, magnetic storage media such as floppy discs, hard disc storage media, magnetic tape; optical storage media such as optical discs or CD-ROMs; electrical storage media such as memory, including RAM, ROM and flash memory; hybrids and combinations of the above such as magnetic / optical storage media.

[0180] Fig. 1A illustrates schematically methods of characterising an evolution of resistance in a population of proliferative cells exposed to an anti-proliferative treatment, methods of performing drug screening, and methods of designing an anti-proliferative treatment, e.g. that reduces the evolution of resistance in cells exposed to the anti-proliferative treatment. At optional step 10, one or more populations of proliferative cells are cultured under a defined scheme (e.g. experimental set up) that comprises exposing the cells to an anti-proliferative treatment. In embodiments, step 10 comprises exposing one or more replicates of the cell population in an in vitro cell culture to the anti-proliferative treatment. Step 10 may comprise culturing a cell population in vitro, using a process comprising: optionally obtaining a parental population of cells by barcoding a population of cells and expanding the barcoded population of cells during a predetermined period of time, obtaining lineage tracing data for a parental population of cells (part of step 12), separating the parental population of cells into a first plurality of replicates, and optionally a second plurality of replicates, exposing each of the first plurality of replicates to the anti-proliferative treatment for one or more predetermined periods of time, optionally separated by one or more predetermined periods of time in the absence of the anti-proliferative treatment, and passaging the cells in each of the first plurality of replicates, and optionally each of the second plurality of replicates, one or more times. In embodiments, lineage tracing data is acquired at passaging time points at step 12 as explained further below. In other words, an initial barcoded pool of cell may be split into replicate populations that are then independently evolved in parallel. The anti- proliferative treatment may be applied periodically, i.e. during multiple Ton periods separated by one or more Toff periods. This may be useful to parameterise the dynamic effect of the antiproliferative treatment on the cells. Barcoding may refer to labelling each cell with a stable genetic barcode system. Examples of such systems include the ClonTracer system (Bhang et al., 2015), the Clonmapper barcode system (Gutierrez et al. 2021), the CaTCH barcode system (Umkehrer et al. 2020) and many others. Growing cells fora predetermined period of time after the cells are barcoded advantageously ensures that most barcodes are represented multiple times. The second plurality of replicates may be control replicates not exposed to the anti-proliferative treatment. In embodiments, evolving barcodes may be used instead of stable genetic barcodes (see e.g. Kalhor et al., 2017). In embodiments, the cells are re-labelled at one or more passaging steps. This can be useful to identify clonal sweeps due to ‘escape transitions’ that occur following a complete loss of diversity before a passaging step.At optional step 12, experimental data comprising total population size and lineage tracing data are measured at one or more time points after the start of exposure of the cells to the anti-proliferative treatment. Alternatively, step 12 may simply comprise obtaining, by a processor, experimental data comprising total population size and lineage tracing data measured at one or more time points after the start of exposure of the cells to the anti-proliferative treatment. The data may have been previously measured and may be obtained by the processor from a computing device, data store, computer-readable memory or user interface. In embodiments, the total cell count data is data that has been obtained by cell counting of a cell population in cell culture, optionally using imaging data, or based on a tumour size determination. In embodiments, the total population size data comprises total cell counts acquired using a non-disruptive technology, e.g. imaging. The total cell count data may be selected from: cell counts that have been obtained from imaging data, cell counts that have been obtained from flow cytometry data, and cell counts that have been obtained from lineage tracing data. Lineage tracing data is also referred to simply as “lineage data”. Lineage tracing data refers to data that characterises a population of cells in terms of a distribution of clonal lineages. In embodiments, the lineage tracing data comprises counts of cells that have the same clonal identity or data derived therefrom. In embodiments, the lineage tracing data comprises counts of cells in respective clonal lineages acquired using a sequencing technology. In embodiments, the lineage tracing data is data acquired at passaging steps applied to an in vitro cell culture. For example, these counts are based on the number of cells that share each of a plurality of heritable barcodes. The lineage data may be selected from: data that has been obtained by extracting and sequencing heritable barcodes (e.g. heritable genetic barcodes) or lineage identities defined by a unique genetic and / or epigenetic mutation profile for each clone of a plurality of clones present at the start of exposure of the cells to the treatment. In embodiments, the experimental data comprises total population size and lineage tracing data measured for a cell population in vitro. In embodiments, the experimental data comprises total population size and lineage tracing data measured at one or more time points after the start of exposure of the cells in each of a plurality of replicates of the cell population exposed to the anti-proliferative treatment. In embodiments, the experimental data comprises total population size and lineage tracing data measured fora plurality of replicates of a cell population derived from the same parent population and exposed to the anti-proliferative treatment in vitro. The use of multiple replicate populations derived from the same parental population enables the calculation of lineage statistics from the lineage tracing data, specifically differences in lineage diversity between replicates. This in turn enables meaningful and computationally efficient comparison of simulated data from the stochastic models to experimental data. Indeed, the present inventors have identified that a lot of important information about resistance evolution can be captured in the between-replicate differences following prolonged exposure to drug. Note that the use of lineage statistics is not a requirement since e.g. full lineage distributions can be used instead. However, the present inventors have found that the use of lineage statistics advantageously enables accurate identification of relevant resistance evolution parameters with far greater efficiency than when using full lineage distribution information. In embodiments, step 12 comprises obtaining both total population size and lineage tracing data at one or more first time points after the start of exposure of the cells to the anti-proliferative treatment, and optionally obtaining total population size data at one or more second time points after the start of exposure of the cells to the anti-proliferative treatment. The one or more second time points may precedethe one or more first time points, or include one or more time points that precede the one or more first time points, and one or more time points that are after at least one of the one or more first time points. The one or more first time points may be time points at which the cell culture replicates are passaged. The one or more second time points may be any time point. For example, total population size data may be acquired using a technology that does not require sampling of the cells. This may be referred to as a non-disruptive (or non-destructive) technology. Thus, total population size measurements may be (but do not need to be) more frequent than lineage tracing data measurements, because the former may not require sampling of the cells. It is advantageous, although not necessary, for the time points at which lineage tracing measurements is obtained to also be associated with total cell population size measurements. This enables models fitted on the lineage tracing data and models fitted on the total population size data to learn from data at the same time point, thereby encouraging alignment between the models simulations of the same experiments. Lineage tracing data may be acquired at time points at which the cell cultures are passaged. For example, cell cultures may be harvested at one or more known times (which may be predetermined times or dynamically determined times, such as e.g. when the cell cultures reach a predetermined density, cell number or level of confluence). A predetermined number of cells may then be seeded in a new cell culture, and the remaining cells may be used to obtain lineage tracing data. In embodiments, the experimental data comprises total population size and lineage tracing data measured at the start of exposure of the cells to the anti-proliferative treatment (tO) and at one or more time points (P1, P2,...) after the start of exposure of the cells to the anti-proliferative treatment. In embodiments, the experimental data comprises total population size data measured at a plurality of time points after the start of exposure of the cells to the anti-proliferative treatment (01, 02,..., P1, P2,...) and lineage tracing data measured at a subset of the same plurality of time points after the start of exposure of the cells to the anti-proliferative treatment (P1, P2,...). Thus, data may be obtained at one or more time points at which only total cell population size data is obtained, in addition to one or more time points at which both total cell population size data and lineage tracing data is obtained. The former may be referred to as "observation” points (O) because these time points do not necessarily require interfering with the cells, which can instead simply be observed using non-disruptive technologies. The latter may be referred to as "passage” points (P) because lineage tracing data is typically acquired using technologies that require cell sampling and is therefore advantageously performed at passaging points, where cells are harvested and re-seeded anyway. Note that there is no requirement for the observation time points to be included. Instead, total population size data can be obtained solely at passaging time points (i.e. at the same time points where lineage tracing data is obtained). However, the present inventors identified that having these additional population size measurements, especially towards the start of the experiment (e.g. before the first passage), beneficially improved the model fitting process, i.e. leading to more accurate and confident model fits. Further, it is not a requirement that all of the total population size data is obtained from the same experiments as the lineage tracing data. For example, a parallel experiment with identical treatment schedule can be run solely for the purpose of obtaining total cell population size data at time points where the lineage tracing data is not available from a main experiment in which such data is collected.

[0181] In embodiments, the methods further comprise an optional step 13 of estimating respective average birth and death rates of the cell population prior to exposure to the treatment. This can be assumed to be equalto bs and ds, (i.e. the birth and death rate of the sensitive subpopulation of cells) e.g. when the proportion of resistant cells is expected to be very low prior to exposure to treatment. Alternatively, this can be assumed to represent a combination of the birth and death rates of the different subpopulations weighted by their respective proportions prior to exposure to treatment (i.e. at t=0), which can be used as a constraint to the model fitting process. These rates can be estimated using any experimental technology known in the art suitable for detecting live / dead (viable / nonviable) numbers or proportions of cells, such as e.g. live imaging, FACS sorting with live / dead cell markers and cell counting technologies that give a viable / nonviable fraction. These rates can also be estimated from lineage tracing data obtained after an expansion step in the absence of treatment for a predetermined period of time, assuming a simple birth / death model with a single population of cells (eq. (1)). The model can be fitted to data using Bayesian inference. The model can be fitted to lineage tracing data using a Bayesian inference framework to identify posterior distributions for the parameters b and d. The model can take into account noise introduced in the lineage tracing data by sampling of cells and sampling of reads to obtain the lineage tracing data. For example, this can use equations (56)-(59) below. In embodiments, the cell population was cultured using a process comprising: obtaining lineage tracing data fora parental population of cells, wherein the parental population of cells has been obtained by barcoding a population of cells and expanding the barcoded population of cells during a predetermined period of time, and step 13 comprises using a distribution of lineages obtained from said lineage tracing data from the parental population to estimate a parental population average birth rate and a parental population average death rate. In embodiments, the observed lineage distributions of the parental population are assumed to be a product of the birth-death process over a known period of time (the predetermined period of time) followed by cell sub-sampling, barcode extraction, amplification and sequencing. In embodiments, said estimating assumes that each cell contains a unique barcode prior to the predetermined period of time, and that therefore the distribution of barcode lineages in the expanded cells (Nt cells) are NO independent realisations of the birth-death process, where NO is the number of cells in the barcoded population prior to the predetermined period of time. In embodiments, cell subsampling and read sampling are both modelled with a Poisson distribution. In embodiments said estimating comprises calculate the variance introduced when sequencing J reads from a pool of K cells sampled from an expanded pool of Nt cells with birth and death rates b and d, and introducing this variance into the probability distribution for sampling K cells by replacing the Poisson distribution with a Negative Binomial distribution. In embodiments, said estimating comprises obtaining a negative binomial distribution giving the probability of sampling j reads from a particular lineage given that there are n cells of the particular lineage after the predetermined period of time, where the probability of there being n cells is given by equation (37)-(40) which depends only on the average birth and death rates in the population. This may then be used to determine the birth and death rates using a Bayesian approach (such as e.g. a full Bayesian statistical inference with MCMC sampling).

[0182] At step 14, the processor fits one or more first dynamic cell population growth models to the total population size data, thereby identifying a first set of values of each of one or more parameters of each of the one or more dynamic cell population growth models that meet a first predetermined model fit criterion, wherein the dynamic cell population growth models represent the growth of a plurality of subpopulations of cells having different responses to the anti-proliferative treatment comprising at least a first subpopulation of cells thatis sensitive to the anti-proliferative treatment and a second population of cells that is resistant to the antiproliferative treatment. Each cell subpopulation may be referred to as a "compartment”, and may be associated with a “phenotype” that characterises its response to the drug. Cells from a first subpopulation may acquire the phenotype of another subpopulation, which is referred to herein as a “transition”. In embodiments, transitions are assumed to be associated with cell births, i.e. when a cell divides it can generate two cells, one of which has the phenotype of another subpopulation. This is compatible with both genetic I epigenetic mechanisms of transition, as well as non-genetic mechanisms such as asymmetric segregation of transcripts during cell division. In other embodiments, transitions are assumed to be independent of cell births. This can be modelled simply by allowing transitions to occur at a transitionspecific rate. By contrast, transitions that are assumed to be associated with cell births can be modelled as occurring at an effective rate that is the product of the birth rate of cells in the subpopulation and a transition specific rate of transition per cell birth. This is also an accepted way to model transitions in cancer evolution, and is also compatible with genetic, epigenetic and non-genetic mechanisms.

[0183] At step 16, the processor fits one or more second dynamic cell population growth models to the lineage tracing data, each second dynamic cell population growth model representing the growth of the same plurality of subpopulations of cells and having the same parameters as a corresponding first dynamic cell population growth model, said fitting comprising evaluating the fit to data of the one or more second dynamic cell population growth models to the lineage tracing data parameterised using parameter values derived from the first set of values, thereby identifying a second set of values of each of one or more parameters of each of the one or more first and second dynamic cell population growth models that meet a second predetermined model fit criterion. The second set of values of one or more parameters comprise: values that are indicative of the effect of the treatment on the growth of the plurality of subpopulations of cells, values that are indicative of the proportions of the one or more subpopulations of cells prior to exposure to the anti-proliferative treatment, and / or values that are indicative of the rate of transition of cells between the plurality of subpopulations of cells comprising a rate of transition between sensitive and resistant subpopulations of cells. In embodiments, for each second dynamic cell population growth model there is a corresponding first dynamic cell population growth model that represents the growth of the same plurality of subpopulations of cells using the same parameters, wherein the second dynamic cell population growth models represents the temporal dynamics of the numbers of cells of each of a plurality of cell lineages that belong to each of the plurality of subpopulations of cells (n is, niR, niE), and wherein the first dynamic cell population growth model represents the temporal dynamics of the total number of cells in each of the plurality of subpopulations of cells (ns, HR, HE). Thus, for each second dynamic cell population growth model there is a corresponding first dynamic cell population growth model that represents the growth of the same plurality of subpopulations of cells using the same parameters. In other words, each pairof first and second dynamic cell population growth models represent the same subpopulations of cells and makes the same assumptions in relation to the dynamics of the growth of these subpopulations and the transitions that are possible between them, parameterised by the same parameters. However, the first and second dynamic cell population growth models in a pair may differ by their model formalism. For example, as described further below, the first dynamic cell population growth model may be at least partially deterministic (e.g. modelled with a system of differential equations, such as ordinary differential equations a least at cellsubpopulation sizes above predetermined thresholds), while the second dynamic cell population growth model may be fully stochastic. This advantageously enables the first fitting step to be used to perform a computationally efficient survey of parameter set values that fit experimental data, while the second fitting step may refine the parameter set using a fully stochastic, more accurate but more computationally expensive model that can make use of lineage tracing data.

[0184] The second set of values of one or more parameters may comprise values that are indicative of one or more of: a proportion of cells that are resistant before the population of cells is exposed to the antiproliferative treatment (p), a rate of transition of sensitive cells to resistant cells (p), a probability of treatment-induced death for resistant cells (qj), a fitness penalty of resistant cells in the absence of treatment (6), a rate at which resistant cells escape a fitness penalty in the absence of treatment (a), a rate of change of a treatment induced effect on a death rate of resistant and / or sensitive cells upon exposure to the treatment (K), and a maximum strength of a treatment induced effect on a death rate of resistant and / or sensitive cells (De, or Dcs and DCR). In embodiments, the second set of values of one or more parameters comprise values that are indicative of each of: (i) a proportion of cells that are resistant before the population of cells is exposed to the anti-proliferative treatment (p), (ii) a rate of transition of sensitive cells to resistant cells (p), (iii) a rate of change of a treatment induced effect on a death rate of resistant and / or sensitive cells upon exposure to the treatment (K), and (iv) a probability of treatment-induced death for resistant cells (qj) and a maximum strength of a treatment induced effect on a death rate of resistant and / or sensitive cells (De), or a maximum strength of a treatment induced effect on a death rate of resistant cells (DCR) and a maximum strength of a treatment induced effect on a death rate of sensitive cells (Dcs). In embodiments, the value of a rate of transition of sensitive cells to resistant cells (p) is indicative of whether transition between sensitive and resistant cells is likely to be due to one or more genetic mechanisms or one or more non-genetic mechanisms. For example, lower values (which can be expressed as probabilities when the rates are associated with cell birth events) are more likely to be associated with genetic mechanisms and higher values are more likely to be associated with non-genetic mechanisms. Thus, in embodiments, genetic mechanisms may be identified as unlikely (i.e. non-genetic mechanisms may be identified as likely) when transitions occur at rates that are too high to be biologically feasible. While there is no hard cut-off for this, when the rates are expressed as probabilities associated with cell birth events the inventors used 10_9to 10-6as a realistic genetic range, 10-3to 10° as a non-genetic range, and intermediate values were deemed not conclusive. In embodiments, a parameter indicative of a proportion of cells that are resistant before the population of cells is exposed to the anti-proliferative treatment (p) may be fitted using a prior distribution for a log transformed version of the variable. In embodiments, a parameter indicative of a rate of transition of sensitive cells to resistant cells (p) may be fitted using a prior distribution for a log transformed version of the variable.

[0185] In embodiments, the methods comprise fitting a plurality of first dynamic cell population growth models and a corresponding plurality of second dynamic cell population growth models. The models may represent different assumptions of how resistance occurs. In such embodiments, the method may further comprise optional step 18 of comparing the plurality of fitted second dynamic cell population growth models using a third predetermined model fit criterion, and identifying a second dynamic cell population growth model thatbest fits the experimental data based on said comparison. The third predetermined model fit criterion may be a model fit criterion that takes into account the numbers of fitted parameters of the models that are being compared. For example, the third model fit criterion may be a DIC score. This is optional as a single pair of first and second dynamic cell population growth models may be fitted.

[0186] In embodiments, the method comprises repeating steps 14 and 16 (and optionally any of steps 10, 12, 13, 18) for two different anti-proliferative treatments. In such embodiments, the method may further comprise step 19 of fitting a cross-resistance model jointly to lineage tracing data comprising lineage tracing data for cells cultured in the presence of a first of the two different anti-proliferative treatments (A) and lineage tracing data comprising lineage tracing data for cells cultured in the presence of a second of the two different anti-proliferative treatments (B). Fitting a cross-resistance model may comprise fitting one or more third dynamic cell population growth models to the lineage tracing data. Each third dynamic cell population growth model represents the growth of a plurality of subpopulations of cells comprising at least a first subpopulation of cells that is sensitive to both anti-proliferative treatments, a second population of cells that is resistant to the first anti-proliferative treatment, a third population of cells that is resistant to the second anti-proliferative treatment, and a fourth population of cells that is resistant to both of the first and second anti-proliferative treatments. Fitting the cross-resistance model may comprise estimating, based on the lineage tracing data and the parameters estimated at steps 14-16 independently for the two different antiproliferative treatments, the value of a cross-resistance parameter (0). The cross-resistance parameter (0) may be a parameter that quantifies the strength of statistical association of inheritance of resistance phenotypes to the two anti-proliferative treatments. The cross-resistance parameter (0) may be a parameter bounded between 0 and 1. The cross resistance parameter may be a parameter (0) that quantifies the departure of: (i) the proportion of cells that are resistant to both anti-proliferative treatments before the population of cells is exposed to any of the anti-proliferative treatments (PAB) and / or (ii) the rate of transition of sensitive cells to cells that are resistant to both anti-proliferative treatments (|JAB), from a model in which resistance to the two anti-proliferative treatments is acquired independently. For example, the parameter (0) may be the parameter 0 such that the proportion of cells that are resistant to both anti-proliferative treatments before the population of cells is exposed to any of the anti-proliferative treatments is given by equation (77), and / or the rate of transition of sensitive cells to cells that are resistant to both anti-proliferative treatments is given by / iAB= (1 - 0)( / z / zB) + ^min G< MB where pA and pB are respectively the rate of transition of sensitive cells to cells that are resistant to the first and second anti-proliferative treatments, and PA and PB are respectively the proportions of cells that are resistant to the first and second antiproliferative treatments before the population of cells is exposed to any of the anti-proliferative treatments. Thus, a third dynamic cell population growth model may be a model as defined in equations (69)-(75) and (77)-(81). Fitting the cross-resistance model may comprise drawing parameter values from posterior distributions obtained at steps 14-16 independently for the two different anti-proliferative treatments and values of the cross-resistance parameter from a prior distribution for the cross-resistance parameter, simulating the third dynamic cell population growth model using said parameter values thereby obtaining simulated lineage tracing data, and identifying a posterior distribution of the value of the cross-resistance parameter based on (i) a comparison between one or more observed diversity statistics computed from the lineage tracing data and corresponding diversity statistics computed from the simulated lineage tracing data, or (ii) a prediction of the posterior distributions from a machine learning model trained using parametervalues and corresponding diversity statistics computed from the simulated lineage tracing data. The one or more diversity statistics may be between drug diversity statistics computed at respective pairs of sampling time points (e.g. passage time points) in lineage tracing data associated with exposure to the first antiproliferative treatment and in lineage tracing data associated with exposure to the second anti-proliferative treatment. The between drug diversity statics may quantify the degree to which shared lineages are enriched between the lineage tracing data associated with exposure to the first and second anti-proliferative treatments. For example, the diversity statistics may be computed using equations (91)-(97). At optional step 20, one or more likely mechanisms of resistance are identified using the parameters of (optionally selected, e.g. in embodiments where a plurality of first dynamic population growth models are used) fitted second dynamic cell population growth models. At optional step 22, the results of any one or more of the preceding steps may be provided to a user, e.g. through a user interface. For example, step 22 may comprise outputting a report comprising one or more results of the method or information derived therefrom. The one or more results of information derived therefrom may include one or more of: a value of one or more fitted parameters, a value of a fit criterion, an indication of an identified dynamic cell population growth model that best fits the experimental data, and an indication of a likely mechanism of resistance. At optional step 24, one or more resistance validation experiments are selected and optionally performed based on the one or more likely mechanisms identified at step 20. For example, in embodiments in which the methods predict a pre-existing stable mechanism of resistance for an individual anti-proliferative treatment, assays that looked for genetic drivers of resistance (e.g. mutations or copy-number changes) can be prioritised. Such embodiments can be identified as embodiments in which a proportion of cells that are resistant before the population of cells is exposed to the anti-proliferative treatment (p) is non-zero and / or a rate of transition of sensitive cells to resistant cells (p) is lower than a predetermined threshold indicative of likely genetic transitions. The predetermined threshold may be e.g. between 105and 107, such as e.g. 106, when the rate is expressed per cell birth, i.e. the rate represents a probability of a sensitive cell birth event producing a resistant cell, or corresponding rates when expressed independently of cell births events. As another example, in embodiments in which the methods predict rapid phenotypic switching between sensitive and resistant cells for an individual anti-proliferative treatment, then scRNA-seq and / or scATAC-seq may be prioritised to identify the transcriptional and / or epigenetic regulators of switching. These assays may also be used to screen for therapeutics that modulate the rate of switching between sensitive and resistant cells. Such embodiments can be identified as embodiments in which a rate of transition of sensitive cells to resistant cells (p) is higher than a predetermined threshold indicative of likely non-genetic transitions. The predetermined threshold may be e.g. between 10-4and 10-2, such as e.g. 103, when the rate is expressed per cell birth, i.e. the rate represents a probability of a sensitive cell birth event producing a resistant cell, or corresponding rates when expressed independently of cell births events. As another example, in embodiments in which the methods predict that the initial response to an individual anti-proliferative treatment is driven by a population of slow-cycling / quiescent cells, this drug-tolerant persister phenotype may be isolated and distinguished from sensitive / fully resistant cells by first looking for DNAand then RNA changes. Such embodiments can be identified as embodiments in which a proportion of cells that are resistant before the population of cells is exposed to the anti-proliferative treatment (p) is non-zero and a fitness penalty of resistant cells in the absence of treatment (6) is non-zero, and / or a rate at which resistant cells escape a fitness penalty in the absence of treatment (a) is non-zero. Thus, in embodiments, the one or more mechanisms identified at step 20 comprise a genetic mechanism and the selected one or moreresistance validation experiments at step 24 comprise a genome sequencing step, optionally whole genome sequencing or single cell sequencing. In embodiments, the one or more mechanisms identified at step 20 comprise a non-genetic mechanism and the selected one or more resistance validation experiments at step 24 comprise a transcriptome and / or chromatin state sequencing step, optionally single cell RNA sequencing step and / or a single cell ATAC-seq step. In embodiments, the one or more mechanisms identified at step 20 include the presence of a subpopulation of cells prior to exposure to the treatment that is resistant to the treatment and experiences a fitness penalty in the absence of the treatment, and the selected one or more resistance validation experiments at step 24 comprise isolating said subpopulation of cells and performing one or more of a genome sequencing step, a transcriptome sequencing step, and / or a chromatin state sequencing step.

[0187] Results of the one or more resistance validation experiments can comprise the identification of one or more biomarkers of resistance. In other words, the methods described herein can be used to identify biomarkers that are associated with resistance to treatment having evolved in a population of cells. This identification is performed significantly more efficiently and at a much earlier stage of drug development than in current state of the art practice, in which resistance is typically only identified at the clinical stage rather than at the preclinical stage, and identification of mechanisms of resistance relies on trial and error and / or exhaustive characterisation in the absence of any information about the mechanisms that are likely to underlie the resistance. These biomarkers can in turn be used as patient stratification and / or monitoring biomarkers, e.g. to detect the emergence of resistance in a patient treated with the anti-proliferative treatment, and / or to detect the presence of a resistant population of cells in a patient prior to treatment of the patient with the anti-proliferative treatment.

[0188] In embodiments, the models are fitted using a Bayesian inference framework as described further below, to estimate the posterior probability distributions of each of the one or more models. In embodiments where multiple first and second dynamic cell population growth models are fitted, these can then be used to compare the different models to identify a model (i.e. a first and corresponding second dynamic cell population growth models) that best fits the data, as well as the most likely values of the parameters of this model given the observed data and assumed prior distributions for these parameters. Note that as will be explained further below, direct comparison between multiple models is optional and instead parameters that best fit the data for each model can be used to characterise evolution of resistance in the population of proliferative cells under the assumptions of the respective model. In embodiments, the prior distributions (also referred to herein simply as "prior”) are selected from those in equations (20) to (27) or (104a)-(104b) or distributions of the same type with different parameters. A prior distribution can be selected by the skilled person (e.g. a user) for each parameter based on prior biological knowledge and / or experimental design constraints, and / or assumptions associated with the parameter. In embodiments, the models comprise a parameter indicative of a proportion of cells that are resistant before the population of cells is exposed to the anti-proliferative treatment (p). Such a parameter may set an initial condition for the proportion of resistant cells at the start of a simulation of a model. In embodiments, this parameter is associated with a uniform prior in negative logarithmic space. The prior may be bounded between first and second values that are biologically grounded. For example, the prior for the parameter indicative of a proportion of cells that are resistant before the population of cells is exposed to the anti-proliferative treatment may bebounded between 0 and 7 in negative log 10 space (e.g. -log10(p)~Uniform(0, 7)). As the parameter p determines the proportion of resistant cells at the time of barcoding, the lower bound may be guided by the maximum number of cells barcoded. In embodiments, the models comprise one or more of: a rate of transition of sensitive cells to resistant cells (p), a rate at which resistant cells escape a fitness penalty in the absence of treatment (a), and a rate of reversion of transition of resistant cells to sensitive cells (a). In embodiment, one or more of said rate parameters are associated with a uniform prior in negative logarithmic space. Each prior may be bounded between first and second values that are biologically grounded. For example, the prior for one or more of said rates may be bounded between 0 and 9 in log 10 space (e.g. -log10(p)~Uniform(0, 9); -log10(o)~Uniform(0, 9); -log10(a)~Uniform(0, 9)), based on prior knowledge that probabilities of transitions of 10-9are lower than previously explored mutation rates in cancer resistance evolution studies. In embodiments, the models comprise a parameter indicative of a fitness penalty of resistant cells in the absence of treatment (5). In embodiments, this parameter is associated with a uniform prior in negative logarithmic space, such as e.g. -log10(5)~Uniform(0,3). This leads to high prior probability for low fitness penalties, such that the evidence in favour of higher penalties (5 > 0.1) must be large in the inference to confidently detect a 'cost' of resistance. This parameter is optional and may be omitted in particular in embodiments in which multiple anti-proliferative treatments are investigated (to identify crossresistance), e.g. for simplicity. Uniform priors may be selected as they are relatively uninformative priors that can be very simply bounded to biologically realistic values. In embodiments, the priors for one or more of the above transition rates parameters may be selected as non-uniform (i.e. more informed) priors, for example when previous evidence indicate that a particular transition is rare. In such cases, a prior that is associated with higher probabilities for lower rates may be selected (i.e. a probability mass skewed towards 0). In embodiments, a prior for the parameter indicative of a fitness penalty of resistant cells in the absence of treatment (5) may be selected as a non-uniform prior. The prior may have a probability mass skewed towards 1, such as a beta prior. This may be useful when a fitness penalty is suspected for biological reasons. The beta prior is naturally bounded between 0 and 1 and can therefore be applied on the untransformed parameter space. In embodiments, the models include a parameter indicative of a probability of treatment-induced death for resistant cells (qj). Such a parameter may indicate the strength of resistance of the resistant cell subpopulation, where ip = 1.0 denotes complete resistance. In such embodiments, the parameter may be associated with a p distribution which favours values closer to 1.0, such as e.g. ip ~P(3.0, 1.0). This may advantageously encourage the methods to identify cases where the ‘resistant’ phenotype (if present) confers some protection against treatment. Note that any prior distribution bounded between 0 and 1 with probability mass skewed to 1 may be used instead of a beta distribution. In embodiments, a separate short drug exposure assay may be used to help better calibrate the ip prior distribution and the proportion of resistance (p) prior distribution. In embodiments, a parameter indicative of a probability of treatment-induced death for resistant cells (qj) is not used and instead a parameter indicative of a maximum strength of a treatment induced effect on a death rate of resistant cells (DCR) is fitted separately from a parameter indicative of a maximum strength of a treatment induced effect on a death rate of sensitive cells (Dcs). In embodiments, the models comprise a parameter indicative of a rate of change of a treatment induced effect on a death rate of resistant and / or sensitive cells upon exposure to the treatment (K). In embodiments, the models comprise a parameter indicative of a maximum strength of atreatment induced effect on a death rate of resistant and / or sensitive cells (De) (or separate parameters for resistant and sensitive cells, as explained above). These parameters may control the behaviour of an effective drug concentration (effective drug concentration (De) and rate of drug accumulation / decay (kappa)). These parameters may be associated with a prior that is constrained to positive values. Any strictly positive distribution may be used, such as e.g. gamma, exponential or log-normal prior distributions. In embodiments, broad T distributions may be used, such as e.g. Dc~T(3.0, 1.0) and / or K~T(2.0, 1.0). In embodiments in which a parameter indicative of a maximum strength of a treatment induced effect on a death rate of resistant cells (DCR) is fitted separately from a parameter indicative of a maximum strength of a treatment induced effect on a death rate of sensitive cells (Dcs), both may use a uniform prior distribution constrained to positive values within a predetermined range. The predetermined range may be the same for both the sensitive cells and the resistant cells.

[0189] In embodiments, a single first dynamic cell population growth model and a single corresponding second dynamic cell population growth model may be fitted. Such embodiments may be useful when the different phenotypes of cells in the population are believed to be known, as are possible transitions between them, and the methods are used primarily to quantify parameters of these phenotypes and transitions. Such embodiments may also be useful when pairs of anti-proliferative treatments are investigated. Indeed, in such situations, where a simple model is believed to fit the data at hand satisfyingly, that model can be used as a basis to obtain a corresponding model that models cross-resistance, for simplicity given the added complexity of modelling multiple resistant subpopulations. In embodiments, a plurality of first dynamic cell population growth models and a corresponding plurality of second dynamic cell population growth models may be fitted. The models differ from each other by the assumptions made in relation to the effect of the treatment on the plurality of subpopulations of cells, and the transitions between plurality of subpopulations of cells. Thus, the models may differ from each other by one or more of: the plurality of subpopulations of cells that are represented, each subpopulation characterised by a different effect of the treatment on the cells birth and / or death rates, and the transitions that are assumed to be present between one or more pairs of subpopulations of cells. Such embodiments may be useful when the different phenotypes of cells in the population and / or the possible transitions between them are not believed to be known, and the methods are used identify these as well as the values of the parameters associated with these phenotypes and transitions. Note that it is possible for there to be more first dynamic cell population growth models than second dynamic cell population growth models, i.e. for some first dynamic cell population growth models to not be explored in the second fitting step. For example, this may be useful when a first dynamic cell population growth model amongst a set of candidate first dynamic cell growth models can be excluded at the first fitting step as not fitting the data as well as other candidate first dynamic cell growth models.

[0190] In embodiments, the one of more first dynamic cell population growth models and the one or more second dynamic cell population growth models are birth-death models that represent the number of cells alive in each of a plurality of subpopulations of cells as a function of time, the plurality of subpopulations of cells having different responses to the anti-proliferative treatment comprising at least a first subpopulation of cells that is sensitive to the anti-proliferative treatment and a second population of cells that is resistant tothe anti-proliferative treatment, wherein the different subpopulations of cells differ by the birth rate in the presence of the treatment, their death rate in the presence of the treatment, their birth rate in the absence of treatment, and / or their death rate in the absence of treatment. In embodiments, the sensitive cells divide with a birth rate bs and die in the absence of treatment with a base death rate ds, and the resistant cells divide with a birth rate bR = bs (1 - 6) and die in the absence of treatment with a base rate d = ds (1 - 5), wherein 6 represents a fitness penalty of resistant cells in the absence of treatment. In embodiments, 5 is not used, i.e. equivalent to 5=0, and bR = bs and dR = ds. In embodiments, the sensitive cells die in the presence of treatment with an effective death rate equals to the sum of a base death rate (ds) and a sensitive treatment effect term (D(t)), and the resistant cells die in the presence of treatment with an effective death rate equals to the sum of a base death rate (dR)and a resistant treatment effect term (D(t)(1 -( / )), wherein the resistant treatment effect term is equal to the sensitive treatment effect term multiplied by the complement of a probability of treatment-induced death for resistant cells (qj). In embodiments, the sensitive treatment effect term is a time dependent parameter that is the product of a maximum strength of a treatment induced effect on a death rate of cells (De) and a time-dependent factor (y(t)) that increases at a first predetermined rate (K) upon exposure to the treatment and decreases at a second predetermined rate (-K) upon withdrawing of the treatment, optionally wherein the first and second predetermined rate have the same absolute value. In other words, a term that can be interpreted as an effective concentration of the drug can impact the death rate of the sensitive cells, and that of resistant cells subject to a resistance parameter qj (which can be bounded between 0 and 1), e.g. using equations (5) and (6) below. In other embodiments, the sensitive cells die in the presence of treatment with an effective death rate equals to the sum of a base death rate (ds) and a sensitive treatment effect term (Ds(t)), and the resistant cells die in the presence of treatment with an effective death rate equals to the sum of a base death rate (dR) and a resistant treatment effect term (DR(t)), wherein the sensitive treatment effect term and the resistant treatment effect terms are time dependent parameters that are each the product of a maximum strength of a treatment induced effect on a death rate of cells respectively for the sensitive cells and the resistant cells (Dcs, DCR) and a time-dependent factor (y(t)) that increases at a first predetermined rate (K) upon exposure to the treatment and decreases at a second predetermined rate (-K) upon withdrawing of the treatment, optionally wherein the first and second predetermined rate have the same absolute value. In embodiments, the effective death rate of sensitive cells is given by equation (5) below, and the effective death rate of resistant cells is given by equation (6) below. In embodiments, the value of D(t) in equations (5) and (6) is given by equations (3) and (4) below. In embodiments, the effective death rate of sensitive cells is given by equation (65') below, and the effective death rate of resistant cells is given by equation (66’) below. In embodiments, the value of Ds(t) and DR(t) in equations (65’) and (66’) is given by equations (3’), (3”) and (4) below. In equation (4), Ton are periods in which the population of cells is exposed to the treatment, and Toff are periods following a Ton period, in which the population of cells is not exposed to the treatment. In other words, predefined treatment windows may be used (Ton), separated by windows of time without treatment (Toff), and the treatment windows may control a factor D(t) that captures the effective concentration of the drug treatment over time. This can effectively encode a lag in the effect of the treatment by expression the effective concentration of the drug treatment over time as a product of a maximum effective concentration of the drug and a time dependent parameter that is scaled between 0 and 1 and decreases linearly whenthe drug is absent and increases linearly when the drug is present, for example both at the same fixed rate (K).

[0191] In embodiments, the first and second dynamic cell population growth models are birth-date models that capture births and deaths of cells in each of the plurality of subpopulations of cells, wherein transitions between specific subpopulations of cells are associated with respective rates and coupled to birth events. In embodiments, the rate of transition between a first subpopulation and a second subpopulation is a parameter between 0 and 1 that represents the probability of a cell in the first subpopulation generating a cell in the second subpopulation upon cell division. In embodiments, the value of the rate of transition between a first subpopulation and a second subpopulation is indicative of the likely mechanism of action of the transition, wherein lower transition rate values are more likely to be associated with genetic or epigenetic mechanisms than higher transition rates, and higher transition rates are more likely to be associated with non-genetic mechanisms than lower transition rates. Non-genetic mechanisms of transition can include e.g. partitioning of gene transcripts during cell divisions. Non-genetic mechanisms of transition can include epigenetic modifications, such as DNA methylation changes or changes to chromatin organisation (e.g. chromatin accessibility changes). For example, the present inventors have demonstrated that the methods described herein can be applied in an embodiment to data from two different cancer cell lines, revealing that data from one cell line (SW6bc) were best explained by a model in which resistant cells do not transition back to sensitive cells, resistance was rare prior to exposure to treatment, and the transition rate from sensitive to resistant cells was low, consistent with a stable, pre-existing molecular driver of resistance. In contrast, data from the other cell line (HCTbc) were better explained by a model in which resistant cells can transition to resistant cells that do not incur a fitness penalty in the absence of treatment, and where fastergrowing, treatment-resistant cells emerged from a more common population of slower-growing, treatmentrefractory cells. In embodiments, the first and second dynamic cell population growth models each comprise a rate of transition of sensitive cells to resistant cells ( ). In embodiments, at least one first dynamic cell population growth model and at least one corresponding second dynamic cell population growth model further comprise a rate of reversion of transition of resistant cells to sensitive cells (a). In embodiments, at least one first dynamic cell population growth model and corresponding second dynamic cell population growth model represent the growth of a subpopulation of cells that are resistant to the treatment and do not incur a fitness penalty in the absence of treatment, and a subpopulation of cells that are resistant to the treatment and do incur a fitness penalty in the absence of treatment, and wherein said first and second dynamic cell population growth models comprise a rate of transition of resistant cells that do incur a fitness penalty in the absence of treatment escape to resistant cells that do not incur a fitness penalty in the absence of treatment, optionally wherein the rate of transition is obtained as the product of a base rate (a) and a time-dependent factor (y(t)) that increases at a first predetermined rate (K) upon exposure to the treatment and decreases at a second predetermined rate (-K) upon withdrawing of the treatment, optionally wherein the first and second predetermined rates have the same absolute value. In embodiments, a first dynamic population growth model comprises a parameter that represents a proportion of cells that are resistant before the population of cells is exposed to the anti-proliferative treatment (p). This may be referred to as a ‘pre-existing resistance fraction’ parameter (p) that controls the initial conditions of the model by setting the proportion of cells with the resistant phenotype before exposure to the anti-proliferative treatment(tO). In embodiments, there are no escaped cells (see below) before exposure to the anti-proliferative treatment. In embodiments, the first dynamic population growth models comprise subpopulation-specific birth and death rates. In embodiments, a first dynamic population growth models comprise a parameter (6) that represents a fitness penalty associated with the resistant phenotype in the absence of the antiproliferative treatment. This may be referred to as a cost parameter. In embodiments, a first dynamic population growth model comprises a rate of transition of sensitive cells to resistant cells (p). This may be referred to as a switching parameter ( ) that controls the probability of cells transitioning from the sensitive to resistant phenotype per cell division. In embodiments, a first dynamic population growth model comprises a parameter that represents a probability of treatment-induced death for resistant cells (qj), optionally wherein i = 1.0 denotes complete resistance, whereas when i = 0.0 the sensitive and resistant cells experience the same level of drug-induced death. This parameter may capture the strength of the resistance of the resistant population to the anti-proliferative treatment. In embodiments, a first dynamic population growth model comprises a variable that represents the strength of a treatment induced effect on a death rate of resistant and / or sensitive cells. In embodiments, said variable is dependent on a rate of change of a treatment induced effect on a death rate of resistant and / or sensitive cells upon exposure to the treatment (K), and a maximum strength of a treatment induced effect on a death rate of resistant and / or sensitive cells (De).

[0192] In embodiments, step 14 of fitting the one or more first dynamic cell population growth models to the total population size data comprises using an approximate Bayesian computation inference approach. The approximate Bayesian computation inference approach may be a simulation-based, likelihood free inference method, such as ABC-SMC (Approximate Bayesian Computation - sequential Monte Carlo) or Neural Posterior Estimation (NPE). In embodiments, step 14 comprises using ABC-SMC. In embodiments, step 14 comprises using a likelihood free method that relies on comparing simulated data to observed data using a distance measure and generating a posterior distribution of parameters that produce simulations that are within a predetermined distance of the empirical data and so satisfy the first predetermined model fit criterion. Specifically, in embodiments step 14 comprises: simulating, at step 14A, the one or more first dynamic cell population growth models using a first plurality of sets of parameter values for each of the one or more parameters of each first dynamic cell population growth model, the first plurality of sets of parameter values comprising the first set of values, thereby obtaining simulated total cell population size data for each of the first dynamic cell population growth models and each of the first plurality of sets of values for each of the one or more parameters; determining, at step 14B, whether each set of values of parameters of each first dynamic cell population growth model meets the first predetermined model fit criterion by comparing the corresponding simulated total cell population size data to the measured total cell population size data; and selecting, at step 14C, the first set of values as a subset of the plurality of sets of values that meet the first predetermined model fit criterion. In embodiments, the first predetermined model fit criterion applies to the value of a distance metric between the measured total population size data and simulated total population size data corresponding to a first dynamic cell population growth model and a set of values for the one or more parameters of the first dynamic cell population growth model. In embodiments, the first predetermined model fit criterion applies to the value of a distance metric between the measured total population size data and simulated total population size data corresponding to a first dynamic cellpopulation growth model and a set of values for the one or more parameters of the first dynamic cell population growth model. For example, the first predetermined model fit criterion may be whether the value of the distance metric is at or below a predetermined value. In embodiments, fitting, by the processor, the one or more first dynamic cell population growth models to the total population size data comprises iteratively performing: (i) at step 14A, simulating the one or more first dynamic cell population growth models using a first plurality of sets of parameter values for each of the one or more parameters of each first dynamic cell population growth model, the first plurality of sets of parameter values comprising the first set of values, thereby obtaining simulated total cell population size data for each of the first dynamic cell population growth models and each of the first plurality of sets of values for each of the one or more parameters; (ii) at step 14B, determining whether each set of values of parameters of each first dynamic cell population growth model meets the first predetermined model fit criterion by comparing the corresponding simulated total cell population size data to the measured total cell population size data; and (iii) at step 14C, selecting the first set of values as a subset of the plurality of sets of values that meet the first predetermined model fit criterion; wherein the first predetermined model fit criterion is updated to be more stringent at each iteration, and wherein the first plurality of sets of parameter values at each iteration is derived from the first set of values obtained at the preceding iteration. The first set of values selected at the last iteration may then be used to derive the parameter values used to fit the one or more second dynamic cell population growth models to the lineage tracing data. The first predetermined model fit criterion may be whether the value of a distance metric between the measured total population size data and simulated total population size data corresponding to a first dynamic cell population growth model and a set of values for the one or more parameters of the first dynamic cell population growth model is at or below a predetermined value. The predetermined value may be lowered at each iteration, thereby making the first predetermined model fit criterion more stringent. For example, in the first fitting step (step 14), the parameters of each said first dynamic cell population growth model can be sampled from a respective prior distribution TT(0) for each parameter, and the model can be simulated with the sampled parameters to generate synthetic total population size data x. This can then be compared to experimental population size data y using a distance function d(x, y). The distance function for use in step 14 and / or 16 can be any distance function known in the art, such as e.g. Euclidian distance (L2 metric) or Manhattan distance (L1 metric). Parameters that produce synthetic data within a tolerance Et can then be retained. This can be performed iteratively, as illustrated on Fig. 1A, with decreasing tolerance Et at each iteration, tightening the acceptance criterion until a target tolerance is reached, thereby approximating the posterior distribution TT(0 | y). In embodiments, the distance is the mean of the distances between all observations, across replicates. In embodiments, the distances are Euclidian distances.

[0193] In embodiments, step 16 of fitting the one or more second dynamic cell population growth models to the lineage tracing data comprises using an approximate Bayesian computation inference approach. The approximate Bayesian computation inference approach may be a simulation-based, likelihood free inference method, such as ABC-SMC (Approximate Bayesian Computation - sequential Monte Carlo) or Neural Posterior Estimation (NPE). In embodiments, step 16 comprises using NPE. In embodiments, step 16 comprises using a likelihood free method that relies on comparing simulated data to observed data using a distance measure and generating a posterior distribution of parameters that produce simulations that arewithin a predetermined distance of the empirical data and so satisfy the second predetermined model fit criterion. In embodiments, the second predetermined model fit criterion applies to the value of a distance metric between the measured lineage tracing data and simulated lineage tracing data corresponding to a second dynamic cell population growth model and a set of values for the one or more parameters of the second dynamic cell population growth model. In embodiments, the second predetermined model fit criterion applies to the value of a distance metric between the measured lineage tracing data and simulated lineage tracing data corresponding to a second dynamic cell population growth model and a set of values for the one or more parameters of the second dynamic cell population growth model. For example, the second predetermined model fit criterion may be whether the value of the distance metric is within a predetermined lowest percentage of values of the distance metric over the parameter values derived from the first set of values with which a second dynamic cell population growth model is simulated (i.e. second plurality of sets of parameter values). Thus, in embodiments, step 16 comprises: simulating, at step 16A, the one or more second dynamic cell population growth models using a second plurality of sets of parameter values for each of the one or more parameters of each second dynamic cell population growth model, the second plurality of sets of parameter values derived from the first set of values (e.g. sampled from posterior distributions obtained at step 14), thereby obtaining simulated lineage tracing data for each of the second dynamic cell population growth models and each of the second plurality of sets of values; determining, at step 16B, whether each set of values of parameters of each second dynamic cell population growth model meets the second predetermined model fit criterion by comparing the corresponding simulated lineage tracing data to the measured lineage tracing data; and selecting, at step 16C, the second set of values as a subset of the second plurality of sets of values that meet the second predetermined model fit criterion. For example, in the second fitting step (step 16), the parameters of each said second dynamic cell population growth model can be sampled from the approximated posterior distribution obtained in the first fitting step (e.g. obtained at the latest iteration) for the corresponding first dynamic cell population growth model. The sampled parameter values can be used to simulate the respective second dynamic cell population growth model. For example, for each replicate (1, 2,..., Z) in the experimental data, at each time point at which lineage data is available (e.g. passage time points (P1, P2,..., Pn)), the second dynamic cell population growth model can be used to generate a lineage distribution corresponding to the counts of lineages (1, 2,..., NO). In embodiments, these distributions were then converted into Z x n arrays for each of one or more lineage diversity statistics (e.g. 2 diversity statistics), as will be explained further below. The diversity statistics can then be normalised between 0 and 1 and the distance between these diversity statistics and the same diversity statistics calculated for the experimental data can be calculated. The closest p% (e.g. 10%) of simulations can be retained to generate a final posterior distribution TT(0 | y) for each second dynamic cell population growth model. In embodiments, step 16 comprises using a likelihood free method that relies on training a conditional density estimator to estimate a probability density function for the parameters of the model given simulated data, and using the trained model to predict a posterior distribution of parameters based on corresponding observed data instead of the simulated data, the posterior distribution thereby satisfying the second predetermined model fit criterion. Thus, in embodiments, step 16 comprises: simulating, at step 16A, the one or more second dynamic cell population growth models using a second plurality of sets of parameter values for each of the one or moreparameters of each second dynamic cell population growth model, the second plurality of sets of parameter values derived from the first set of values (e.g. sampled from posterior distributions obtained at step 14), thereby obtaining simulated lineage tracing data for each of the second dynamic cell population growth models and each of the second plurality of sets of values (e.g. see equation (101)); training, at step 16B, one or more machine learning models to predict a conditional density of the parameters given the simulated lineage tracing data (q(xIsxm)> where >x is the set of parameters and Ssimx is the simulated lineage tracing data), and using, at step 16C, the trained machine learning models to predict a posterior distribution of the parameters given the observed lineage tracing data (p(< P I sbs) where Sobsx is the observed lineage tracing data), any set of values drawn from these posterior distributions (including but not limited to a maximum density estimate or samples from the distributions) meeting the second predetermined model fit criterion. The same considerations described here in relation to step 16 apply equally to the step 19 of fitting the one or more third dynamic cell population growth models to the lineage tracing data, which can also comprise using an approximate Bayesian computation inference approach as described above. In embodiments, the one of more first dynamic cell population growth models are hybrid differential equation and stochastic jump process models. In embodiments, the one of more second dynamic cell population growth models are stochastic agent based models. For modelling the behaviour of individual cell lineages when exposed to treatment, embodiments of the disclosure model lineage tracing in silico using agentbased models that track each individual cell in a system and can assign unique cell population identities (also referred to herein as cell ‘barcodes’, by reference to the fact that they can be experimentally realised using barcoding libraries, even though this is not a necessary feature for example in the context of naturally genetically heterogeneous cell populations such as in tumours) that are inherited during cell divisions. Such models may assign and track the phenotype of each cell (‘agent’), i.e. recording which subpopulation of cells each cell belongs to. Modelling individual cell lineage evolution in such embodiments is performed using a model where all lineages are represented, which precludes the use of (deterministic) approximations to average population dynamics. By contrast, a purely stochastic agent-based model can track both the phenotype and the lineage relationship of all cells in the system. However, such models are extremely computationally expensive to run. Therefore, extensive simulation as is required for exhaustive parameterisation of the models (e.g. using approximate Bayesian inference as described further herein) is impractical. However, the present inventors have realised that the benefits of such models (which are both more exact than population models and fitted to detailed lineage tracing data) can still be realised if corresponding simpler models such as the hybrid models described herein are fitted to population level data (i.e. total cell counts) to identify regions of the parameter space that are more likely to fit the experimental data. These regions of the parameter space can then be characterised in further detail using the agent based models and the lineage tracing data. In embodiments of the second dynamic cell population growth models, each cell is characterised by a unique set of attributes: a phenotype (i.e. which compartment I subpopulation the cell belongs to - e.g. S, R or E), a lineage tag (i.e. each cell is assigned a unique lineage tag i at the beginning of the experiment, where i = 1, 2, 3,..., NO (NO is the number of uniquely barcoded cells when t=0. Such models therefore track the number of cells ns' with lineage tag i that have phenotype ‘S’, the number of cells HR' with lineage tag i that have phenotype ‘R’, etc. for each subpopulation of cells considered.In embodiments, a hybrid differential equation and stochastic jump process model is a model that represents the growth of each subpopulation of cells using a stochastic process when the total number of cells in the subpopulation is below a predetermined threshold, and using a deterministic differential equation model when the total number of cells in the subpopulation is at or above the predetermined threshold. In other words, the first dynamic cell population growth models may be hybrid models where each phenotypic compartment is modelled as a stochastic jump process at small population sizes (i.e. cells are born in each compartment at each iteration with a probability equal to a birth rate - e.g. S— > S+S with rate bs where S is a sensitive cell and bs is the birth rate of sensitive cells, cells die in each compartment with a probability equal to a death rate - e.g. S— >0 with rate ds where ds is the birth rate of sensitive cells, and cells transition to another compartment at each iteration with a birth event with a transition rate - e.g. S— > S+R with rate bS* where R is a resistant cell and is a rate of transition of sensitive cells to resistant cells) but a deterministic ordinary differential equation (ODE) at large population sizes (i.e. the number of cells in each compartment is given by an ODE that includes terms for each of the above events parameterised by the same rates, e.g. dns / dt, where ns is the number of sensitive cells at time t, is defined by an equation comprising a sum of terms capturing rates of births (ns*bs), deaths (-ns *ds) and transitions (e.g. transitions from sensitive to resistant: - ns*bs*p). The switch between a stochastic jump process model and a deterministic model may be Implemented as nx(t) > Nswitch, where x refers to any subpopulation being modelled, e.g. S, R, etc. (i.e. the switch may be evaluated separately for each compartment). In other words, the methods may switch to an ODE model for a compartment when the compartment size exceeds a predetermined threshold. The value of Nswitch is a predefined value that may be set by a user. Approaches for setting the values of Nswitch are described in the examples of the present disclosure. For example, Nswitch may be set at a value that balances computational efficiency and ability of the hybrid model to recover true parameter values in simulated data from a corresponding fully stochastic agent model. Without wishing to be bound by theory, the present inventors observed that values of Nswitch as low as 10 (i.e. switching to a deterministic model as soon as a subpopulation comprised at least 10 cells) was enough to recover true parameter values. Values of Nswitch=200 were found to reliably recover true parameter values and lead to computationally tractable simulations. Thus, the predetermined threshold (Nswitch) may be set to about 10, 20, 50, 100 or 200 cells.

[0194] In embodiments, the lineage tracing data (e.g. used at steps 16 and / or 18 and / or 19) comprises diversity statistics derived from counts of cells that have the same clonal identity. For example, the experimental data may comprise total population size and lineage tracing data measured at one or more time points after the start of exposure of the cells in each of a plurality of replicates of the cell population exposed to the antiproliferative treatment, and the diversity statistics may comprise a first statistic that is representative of lineage diversity within the replicates and a second statistic that is representative of lineage diversity between the replicates. In embodiments, the diversity statistics are Hill Diversity statistics. To enable the interpretable comparisons of the lineage distribution over time (clone size distribution over time, which is very high-dimensional data), embodiments of the disclosure condense the lineage relationship distributions from each of a plurality of simulation corresponding to replicate cell populations from the same parental population into 2-dimensional diversity statistics that capture the diversity within replicates and the diversity between replicates simultaneously. The diversity statistics may be calculated using equations (32) to (36).Indeed, through exploration of various models the present inventors identified that there were two important axes of information in the lineage distributions: first the number of cell lineages that survive within each of the replicate sub-populations, and second, the relationship of the cell lineages when comparing between the replicate sub-populations. They further identified that these aspects could be captured using respective ‘Hill Diversity Indices’. Thus, embodiments of the disclosure use as a diversity statistic the within-population diversity calculated using equation (32), or using equation (33). Instead or in addition to this, embodiments of the disclosure use as a diversity statistic the between-population diversity dissimilarity that compare, such as e.g. uses the ratio of (e.g. equation (34), the total diversity of order q when pooling all of a number I of replicate cell populations (qD(y) - eq. (35)), to the mean diversity of order q of all I replicate cell populations (qD(or) - eq. (36)). Thus, embodiments of the disclosure use as a diversity statistic a between-population diversity dissimilarity calculated using equations (34)-(36). Advantageously, Hill Diversity Indices capture a variety of diversity indexes used to model populations (e.g. Shannon, Simpson etc.) by simply modifying the q parameter (see below). Therefore, the use of Hill diversity indices encompasses use of any indices that are special cases of the Hill diversity indices and is flexible enough to use whichever index a user may want to use, by setting the value of q. For example, equation (33) is a special case of equation (32) in which q=1, resulting in the Shannon diversity. In embodiments, the value of the parameter q in the Hill within-population diversity and / or the Hill between-population diversity is set to a predetermined value. The predetermined value may be a user provided value or a default value. The predetermined value may be at least 1. Values of 1 and above advantageously reduce the effect of low frequency noise lineages, such as artifact lineages that appear due to preparation / sequencing artefacts and that artificially inflate the diversity. However, values of q that are too high (e.g. above 100) reduce the ability of moderate frequency lineages to impact the diversity statistics. Thus, the parameter q can be set to a predetermined value that balances i) reducing the impact of low frequency lineages and ii) allowing moderate frequency lineages to impact the diversity statistics. The present inventors have found values of q=2 to be suitable in the tested system, but an appropriate value may depend on the circumstances, such as e.g. the expected occurrence of noisy low frequency lineages. In embodiments, a value of 1 is selected between 1 and 100, between 2 and 50, between 2 and 20, between 2 and 10, or about 2. As explained above, in embodiments in which step 19 is performed, this can use between drug diversity statics may quantify the degree to which shared lineages are enriched between the lineage tracing data associated with exposure to the first and second anti-proliferative treatments. For example, the diversity statistics may be computed using equations (91)-(97).

[0195] In embodiments, growth in the first dynamic population growth models is modelled as density-dependent growth. Density dependent growth models are known in the art and include e.g. logistic growth, theta logistic growth, and the Richards model. In embodiments, growth in the first dynamic population growth models is modelled as density-dependent growth with a fixed carrying capacity (which is denoted herein as “K”). Such embodiments are particularly useful when modelling in vitro cultures, and specifically when the cell population density is not so low that cell number dependent growth inhibition can be ignored. Indeed, due to the limitation on population size imposed by the volume of tissue culture vessels in vitro, growth in all models may be advantageously modelled as logistic growth with a fixed carrying capacity, K. The carrying capacity may be set to a constant value guided by empirical data. For example, the number of cells atconfluence or at a predetermined confluence percentage (e.g. a predetermined confluence percentage at which passaging is performed) can be recorded, and used to derive an estimate of the cell number at confluence, which is equal to the carrying capacity K. In embodiments, cell populations may be grown when obtaining the experimental data until a maximum population size (Nmax) is reached (where Nmax < K). In other words, the experimental data may be from cell populations that have been grown such that no cell culture reached the carrying capacity. Logistic growth can be modelled by scaling the birth and death rates by the density dependence term (1 - N / K), where K is the carrying capacity and N is the total number of cells in the cell culture. Theta logistic growth can be modelled by scaling the birth and death rates by the density dependence term (1 - (N / K)e), where K is the carrying capacity, N is the total number of cells in the cell culture, and 0 is a parameter that changes the nature of the dependence between density and growth. For example, 0>1 weakens the density dependence compared to the logistic growth model at low N, and 0<1 strengthens the density dependence at N< K compared to the logistic growth model. The Richards model is a generalisation of the logistic model in which the growth curve is allowed to be asymmetrical about the point of inflection (see e.g. Tjorve & Tjorve, 2010).

[0196] In embodiments in which the experimental data has been measured from cell populations obtained by separating a parental population of cells into a first plurality of replicates and exposing each of the first plurality of replicates to the anti-proliferative treatment for one or more predetermined periods of time, fitting the first and second plurality of dynamic population growth models may comprise sampling cells from a parental population comprising a proportion of resistant cells set by a fitted parameter p to obtain each of a first plurality of replicates comprising a number of cells matching the number of cells in the first plurality of replicates (and similarly for each of a second plurality of replicates if used), and sampling cells from each replicate at each passage time points to obtain a number of cells in each replicate corresponding to the number of cells in the replicate after passaging. Cell cultures can be passaged one or more times during culture, as known in the art, to avoid reaching numbers of cells above a predetermined value. Passaging steps may be simulated in the first and second dynamic cell population growth models by sampling a predetermined number of cells at a passaging time. Sampling may be implemented in the first dynamic cell population growth models by weighting the probability of sampling each cell subpopulation by its frequency at the time of passaging. Sampling may be implemented in the second dynamic cell population growth models by drawing a predetermined number of times (referred to as Nseed by reference to the number of cells reseeded during the passaging step) from the individual cells represented in the models at the time of passaging, without replacement.

[0197] In embodiments, the anti-proliferative treatment comprises one or more anti-proliferative compounds or compositions, wherein the one or more anti-proliferative compounds or compositions comprise one or more active agents (also referred to herein as ‘‘anti-proliferative agents”, or “cytotoxic agents”) selected from: small molecules, large molecules (e.g. antibodies, antigen binding molecules, peptides, nucleic acids, optionally miRNA, siRNAs and gene editing constructs), combinations or large and small molecules (e.g. antibody drug conjugates, bicycle peptide-drug conjugates), cells (e.g. T cells or modified T cells, optionally CAR-T cells or TCR-T cells), and targeted protein degradation therapies (e.g. PROteolysis TArgetingChimeras (PROTACs), molecular glues, Lysosome-Targeting Chimaeras (LYTACs), and Antibody-based PROTACs (AbTACs)).

[0198] In embodiments, the first dynamic population growth models comprises a model (referred to herein as ‘‘model A”) that represents a first subpopulation of cells that is sensitive to the anti-proliferative treatment and a second population of cells that is resistant to the anti-proliferative treatment, wherein the resistant cells incur a fitness penalty in the absence of the anti-proliferative treatment, and wherein the sensitive cells can transition to resistance cells but resistant cells cannot transition to sensitive cells. The second dynamic population growth models may in such embodiments comprise a corresponding model (model A). Note that all of the specific models described herein are exemplary, and any dynamic cell population growth models that are suitable for modelling the number of cells and lineages of cells in a cell population including at least a sensitive subpopulation and a resistant subpopulation may be used. In embodiments, the firs dynamic population growth models comprise a model A that is a hybrid of a deterministic model specified by equations (7)-(8) below and a stochastic jump model specified by equation (9) below, and the second dynamic population growth models comprise a stochastic agent based model corresponding to the model specified by equation (9) below, where S is a sensitive cell, R is a resistant cell, nS is the number of sensitive cells, nR is the number of resistant cells, t is time, bs is a birth rate of sensitive cells, ds is a death rate of sensitive cells, bR is a birth rate of resistant cells, dR is a death rate of resistant cells, K is a carrying capacity, N is the total number of cells given by N=ns+nR, and is a rate of transition of sensitive cells to resistant cells

[0199] dnsN

[0200] — = (nsbs- nsds- nsbs\i) ■ (1 - -) (7) dnRN

[0201] = (.nRbR~ nRdR+ nsbs\i) ■ (1 - -) (8)

[0202]

[0203] S — * S + S, bs

[0204] S -» 0, ds

[0205] S -> S + R, bs\ (9)

[0206] R -> R + R, bR

[0207] R - * 0, dR

[0208] In embodiments in which step 19 is performed, this may comprise a derivative of model A as defined in equations (69)-(75) and (77)-(81).

[0209] In embodiments, the first dynamic population growth models comprises a model (referred to herein as “model B”) that represents a first subpopulation of cells that is sensitive to the anti-proliferative treatment and a second population of cells that is resistant to the anti-proliferative treatment, wherein the resistant cells incur a fitness penalty in the absence of the anti-proliferative treatment, and wherein the sensitive cells can transition to resistance cells and resistant cells can transition to sensitive cells. The second dynamic population growth models may in such embodiments comprise a corresponding model (model B). In embodiments, the first dynamic population growth models comprise a model B that is a hybrid of adeterministic model specified by equations (10)-(11) below and a stochastic jump model specified by equation (12) below, and the second dynamic population growth models comprise a stochastic agent based model corresponding to the model specified by equation (12) below, where S is a sensitive cell, R is a resistant cell, nS is the number of sensitive cells, nR is the number of resistant cells, t is time, bs is a birth rate of sensitive cells, ds is a death rate of sensitive cells, bR is a birth rate of resistant cells, dR is a death rate of resistant cells, K is a carrying capacity, N is the total number of cells given by N=ns+nR, p is a rate of transition of sensitive cells to resistant cells, and a is a rate of transition of resistant cells to sensitive cells

[0210] dnsN — = (nsbs- nsds- nsbs\k + nRbfio-) ■ (1 - -) (10) at dnRN — = (nRbR- nRdR+ - nRbRa) ■ (1 - -) (11)

[0211]

[0212] at K S — > S + S, bs

[0213] s 0, ds

[0214] S -» S + R, bs\i

[0215] (12) R -> R + R, bR

[0216] R - * 0, dR

[0217] R -> R + S, bRa

[0218] A similar extension of model B in which single and double resistant populations are present can be used in embodiments in which step 19 is performed, following the same principles set out in the model defined in equations (69)-(75) and (77)-(81), simply adding a rate of transition of each population of resistant cells to sensitive cells (i.e. OA, OB, OAB).

[0219] In embodiments, the first dynamic population growth models comprise a model (referred to herein as “model C”) that represents a first subpopulation of cells that is sensitive to the anti-proliferative treatment, a second population of cells that is resistant to the anti-proliferative treatment and that incur a fitness penalty in the absence of the anti-proliferative treatment, and a third population of escaped cells that is resistant to the anti-proliferative treatment and that do not incur a fitness penalty in the absence of the anti-proliferative treatment, wherein the sensitive cells can transition to resistance cells and resistant cells can transition to sensitive cells, and wherein the resistant cells can transition escaped cells but escaped cells cannot transition to resistant cells. The second dynamic population growth models may in such embodiments comprise a corresponding model (model C). In embodiments, the first dynamic population growth models comprise a model C that is a hybrid of a deterministic model specified by equations (13)-(15) below and a stochastic jump model specified by equation (16) below, and the second dynamic population growth models comprise a stochastic agent based model corresponding to the model specified by equation (16) below, where S is a sensitive cell, R is a resistant cell, E is an escaped cell, nS is the number of sensitive cells, nR is the number of resistant cells, nE is the number of escaped cells t is time, bs is a birth rate of sensitive cells, ds is a death rate of sensitive cells, bR is a birth rate of resistant cells, dR is a death rate of resistant cells, bE is a birth rate of escaped cells, dE is a death rate of escaped cells, K is a carrying capacity, N is the total number of cells given by N=ns+nR+nE, is a rate of transition of sensitive cells to resistant cells,a is a rate of transition of resistant cells to sensitive cells, and ay(t) is a rate of transition of resistant cells to escaped cells, where y(t) is the time dependent factor that captures the increase or decrease in effective drug concentrations, optionally wherein y(t) is given by equation (4) and / or wherein o+a<1.0:

[0220] dns— = (nsbs- nsds- nsbs\k + nRbfi<r) ■ (1 - (13) at dnR— = (nRbR- nRdR+ nsbs\k - nRbRa - nRbR(ay(tyf) ■ (1 - (14)

[0221] N

[0222] — = (nEbE- nEdE+ nR(ay(t))) ■ (1 - -) (15)

[0223]

[0224] S — > S + S, bs

[0225] S -> 0, ds

[0226] S -» S + R, bs\i

[0227] R R + R, bR

[0228] R - * 0, dR(16) R -> R + S, bRa

[0229] R -> R + E, bRay(t')

[0230] E -» E + E, bE

[0231] E ^ 0, dE

[0232] dy=( K t E Ton

[0233]

[0234] dt [ -K t G In embodiments, one or more of the first dynamic population growth models represent an effect of the antiproliferative treatment on cells in each subpopulation using an effective death rate that depends on a base death rate and a treatment effect term, wherein the treatment effect term is a time depend term that represents an exposure dependent effect of the anti-proliferative treatment on the cells death rate. In embodiments, the treatment effect term comprises a time dependent term that represents the effective concentration of the anti-proliferative treatment and a subpopulation specific factor that is a parameter that captures the subpopulation specific effect of the anti-proliferative treatment on the cell subpopulation. In other words, therapy can be modelled by modifying the death rate of cells. For example, for sensitive cells, an effective death rate can increase in the presence of the treatment as a function of D(t), with D(t) representing the effective treatment concentration at time t. The term D(t) can be expressed as a function that assumes a pharmacokinetic model where the focus is on the change in the realized effect of the drug on cells, rather than directly modelling the drug concentration. For example, D(t) can be the product of a parameter De and a function y(t) with a constant rate of change K that regulates the strength and rate of accumulation / reduction of this effect, respectively. In other embodiments, as explained elsewhere herein, D(t) is expressed separately for the sensitive and resistant cells as the product of a parameter DCR or Dcs, respectively for resistant and sensitive cells, and a function y(t) with a constant rate of change K that regulates the strength and rate of accumulation / reduction of this effect, respectively

[0235] The methods described herein find uses in a variety of contexts. The methods of the present disclosure find applications in the context of drug development, and specifically in the context of pre-clinical evaluationof candidate drugs, and drug combinations, for the treatment of diseases associated with abnormal cell proliferation, such as e.g. cancer. For example, the methods described herein can be used in comparing different disease models (e.g. different genetic backgrounds, for example - including but not limited to different genetic background obtained through genetic perturbation screens, such as e.g. CRISPR screens) and their propensity for resistance to a given drug / drug combination. This can be used for example to identify biomarkers for cells (and patients) that are likely to respond or not respond to the drug I drug combination (i.e. responders / non-responders) due to the evolution of resistance. This can also be used to identify the mechanisms that may be involved in resistance, and to design therapies that target these mechanisms. As another example, the methods described herein can be used in comparing the chronology of different drugs to quantify their propensity for resistance. For example, this can be used to determine whether drug A first followed by drug B leads to a lower chance of resistance than B then A. This can in turn be used to identify beneficial drug administration regimens. As another example, the methods described herein can be used in the identification of biomarkers for resistance by enabling more rapid identification of the subsequent functional assays necessary to characterise the molecular drivers of resistance. These can be used to design drugs that target the identified mechanisms of resistance, and / or to identify biomarkers of cells I patients that are unlikely to respond to a treatment, for example to set inclusion / exclusion criteria in clinical trials of the treatment and / or patient selection criteria in the clinic. A specific example is illustrated on Fig. 1B. As illustrated on Fig. 1B, a pre-clinical cancer cell model (e.g. a cell line representative of a cancer type to be treated with a candidate drug) can be labelled with a lineage tracing library, then expanded fora period of time in the absence of treatment to obtain a parental population (labelled “POT”). The parental population can be separated into a plurality of replicates here labelled DT1 to DT4. Each replicate can be subject to a predetermined drug treatment scheme, which can comprise periods of exposure to the drug (Ton) and periods in which the cells are cultured in the absence of the drug (Toff). The cells in each replicate can be observed in a non-disruptive manner to obtain total cell counts at a plurality of observation times (01, 02, 03, etc.). The cells can be sampled to read the barcodes from the lineage tracing library by sequencing, thereby obtaining lineage tracing data (and also optionally total cell counts) at a plurality of sampling times. Here the sampling times are illustrated as passaging times (P1, P2, P3, etc.) as the sampling is assumed to be performed at the same time as passaging the cells to avoid additional unnecessary disruption to the cells. The total cell count data (change in cell population size) and the lineage tracing data (sequenced barcode distribution, i.e. counts of each of a plurality of barcodes in the population) are used to fit a plurality of dynamic cell population growth models (here illustrated as candidate models A, B, C, although any number of models can be fitted, including a single model as there is also a benefit in identifying the parameters of a single model instead of comparing multiple different fitted models. Model fitting can be done in a Bayesian inference framework. Model fitting is advantageously done in a likelihood-free Bayesian approach, i.e. a Bayesian approach that does not require the use of a closed form likelihood function. In embodiments, model fitting is done using an approximate Bayesian computation (ABC) approach. In embodiments, model fitting is done in a Bayesian inference framework in which candidate sets of parameters are sampled from prior distributions, then the models are simulated with these distributions, and the counts of sets of candidate parameters that satisfy a model fit criterion (e.g. a criterion that applies to the distance between the simulated and observed data) are used to estimate a posteriordistribution for the parameters of the models. These estimated posterior distributions (one for each model that has been fitted, here illustrated as A, B, C) can then be used to compare the models and identify a model that is most likely given the observed experimental data and the chosen priors. The fitted parameters comprise parameters that quantify different behaviours within each model (such as the proportion of resistance, the stability of the resistant phenotype, the presence of quiescent persister cells, etc.), and different models include a different set of these behaviours. The selection of model and the values of the parameters obtained by fitting the selected model to the observed data represents a quantitative and qualitative characterisation of the evolution of resistance in the pre-clinical cell model. This can be used to guide subsequent functional characterisation (e.g. determining whether the resistant cells would be best characterised using genomic assays such as whole genome sequencing, or using non-genomic assays such as transcriptomics assays, e.g. single cell RNA sequencing, scRNA-seq). Instead or in addition to this, this information can be used to quantitatively compare different drugs, such as e.g. in terms of their propensity for resistance. This can be used to prioritise pre-clinical drug candidates for further assessment, e.g. to select drugs I drug combinations that are less likely to be associated with the development of resistance for further characterisation including clinical testing. By contrast, in current practice pre-clinical drug candidates are typically only assessed for safety and efficacy before being progressed to clinical trials, where evolution of resistance may be observed and may eventually lead to failure of the drug to achieve desired clinical outcomes. As such, many drugs are progressed to expensive and time-consuming clinical trials that eventually lead to resistance, where this could have been spared or another candidate prioritised (possibly despite lower efficacy) that is less likely to lead to resistance. This can also be used to select drugs for further assessment even when their pre-clinical efficacy may not be higher than an existing (e.g. standard of care) drug, when there is evidence that the selected drugs are less likely to be associated with resistance that has been observed clinically for the existing drug. By contrast, in current practice pre-clinical drug candidates may not be progressed if they do not demonstrate an increased efficacy compared to a comparative drug, even though the candidates may in fact have been associated with better overall outcomes clinically (where resistance evolution becomes apparent). This can also be used specifically to assess drug combinations and compare their propensity for resistance. Indeed, most treatments don’t enter clinical trials as a monotherapy, and characterising candidate combinations for their propensity for resistance can help position a new drug alongside existing drugs, prioritising combinations that are more likely to be successful in the clinic. As another example, the methods described herein can be used in investigating pairs of drugs to quantify their propensity for cross-resistance, i.e. the likelihood of a single mechanism of resistance evolving which confers resistance to both drugs.

[0236] Thus, also described herein is a method of screening a plurality of candidate anti-proliferative treatments, the method comprising: (i) characterising an evolution of resistance in a population of proliferative cells exposed to each anti-proliferative treatment using the methods described herein; and (ii) comparing the results of said characterising to identify anti-proliferative treatments of the candidate anti-proliferative treatments that are less likely to be associated with an evolution of resistance. The method may further comprise: (iii) selecting one or more candidate anti-proliferative treatments for further characterisation using the results of said comparing. The further characterisation may comprise pre-clinical and / or clinical characterisation. The plurality of candidate anti-proliferative treatments may differ by the identity of each ofone or more anti-proliferative agents that form part of the anti-proliferative treatments. Instead or in addition to this, the plurality of candidate anti-proliferative treatments may differ by the schedule of exposure of the cells to one or more anti-proliferative agents that form part of the anti-proliferative treatments. For example, a first candidate anti-proliferative treatment may comprise exposing the cells to drug A, then exposing the cells to drug B. Another candidate anti-proliferative treatment may comprise exposing the cells to drug B, then exposing the cells to drug A. Yet another candidate anti-proliferative treatment may comprise exposing the cells to drugs A and B simultaneously. One or more of these candidate anti-proliferative treatments may be compared using the methods described herein.

[0237] Also described herein is a method of identifying a combination of anti-proliferative treatments that has a reduced risk of evolution of resistance compared to a subset of the anti-proliferative treatments in the combination, the method comprising: (i) characterising an evolution of resistance in a population of proliferative cells exposed to the subset of the anti-proliferative treatments using the methods described herein; and (ii) characterising an evolution of resistance in a population of cells exposed to one or more combinations of anti-proliferative treatments comprising respective additional anti-proliferative treatments in addition to the subset of the anti-proliferative treatments using the methods of any preceding claim and comparing the results of said characterising to identify combinations of anti-proliferative treatments that are less likely to be associated with an evolution of resistance than said subset of anti-proliferative treatments; and / or identifying and performing one or more validation experiments based on the results of the characterising in i, and identifying one or more additional anti-proliferative treatments that target a resistance evolution mechanism identified using said validation experiments.

[0238] Also described herein is a method of identifying one or more biomarkers of resistance of a population of proliferative cells to an anti-proliferative treatment, the method comprising: (i) characterising an evolution of resistance in one or more different populations of proliferative cells exposed to the anti-proliferative treatment using the methods described herein; and (ii) identifying one or more likely mechanisms of resistance using the parameters of (optionally selected) fitted second dynamic cell population growth models for each of the one or more different populations of proliferative cells and performing one or more resistance validation assays based on the identified likely mechanisms of resistance to identify one or more biomarkers of resistance; and / or comparing the results of said characterising between a plurality of the different populations of proliferative cells to identify one or more biomarkers associated with one or more of the plurality of populations of proliferative cells likely to be associated with an evolution of resistance.

[0239] Figure 2 illustrates schematically an exemplary system according to the disclosure. The system comprises a computing device 1, which comprises a processor 101 and computer readable memory 102. In the embodiment shown, the computing device 1 also comprises a user interface 103, which is illustrated as a screen but may include any other means of conveying information to a user such as e.g. through audible or visual signals. In the illustrated embodiment, the computing device 1 is operably connected, such as e.g. through a network 6, to a cell culture system comprising a cell culture housing 2 and one or more sensors 3, and to a sequence data acquisition means 5, such as e.g. a sequencer. Alternatively, the computer 1 may be connected to a computer readable memory 4 (such as a database, other computing device or data store) from which data collected by the one or more sensors 3 and / or the sequencer 5 can be obtained bythe computer 1. The cell culture housing 2 may be an incubator or any other kind of housing suitable for live cell culture in a culture dish or vessel. The cell culture system may be an integrated system comprising a cell culture housing and at least one sensor, such as e.g. an Incucyte® live-cell analysis system. The computing device may be a smartphone, tablet, personal computer or other computing device. The computing device is configured to implement a method as described herein. In alternative embodiments, the computing device 1 is configured to communicate with a remote computing device (not shown), which is itself configured to implement a method as described herein. In such cases, the remote computing device may also be configured to send the result of the methods to the computing device 1. Communication between the computing device 1 and the remote computing device may be through a wired or wireless connection, and may occur over a local or public network such as e.g. over the public internet. Each of the sensor(s) 3 and sequencer 5 may be in wired or wireless connection with the computing device 1 or with another computing device (not shown) to which the computing device 1 is connected or from which the computing device 1 is able to acquire data (e.g. via a database 4 to which the computing device 1 is also connected). Thus, the connection between the computing device 1 and the sensor(s) 3 and sequencer 5 may be direct or indirect (such as e.g. through a remote computer). In alternative embodiments, the computing device 1 is configured to implement a method as described herein, using images received from a data store or remote computing device (such as e.g. a computing device associated with the cell culture system). Thus, the computing device 1 may not be directly connected to the cell culture system and / or the sequencer 5. The one or more sensors 3 may comprise at least one sensor configured to acquire cell count data of one or more cell population(s) in the cell culture housing, or data from which cell count data can be obtained. This sensor may be an imaging means, such as e.g. a microscope (e.g. a phase contrast microscope or bright-field microscope). The sequencer 5, also referred to as sequence data acquisition means, is configured to acquire sequence data associated with a sample of cells collected from the incubator 2. The data may be in the form of sequencing reads. The sample may be subject to one or more of cell separation, DNA purification, and library preparation prior to sequencing, as known in the art. The data from the sequencer 5 may be used by the computing device 1 or by another computing device (not shown) to obtain lineage tracing data in the form of counts of cells associated with each of a plurality of clonal identities (also referred to herein as lineages). The measurements from the sensors 3 and the sequencer 5 are communicated to the computing device 1 or the data store 4, which may store the data permanently or temporarily in memory 102. The computing device memory 102 may store one or more first dynamic population growth models and one or more second dynamic population growth models as described herein. The processor 101 may execute instructions to obtain lineage tracing data from the sequencing data and / or to obtain total cell count data from raw data from a sensor 3 (e.g. images). The processor 101 may execute instructions to fit the first and second dynamic population growth models as described herein using said data. The processor 101 may execute instructions to identify a likely mechanism of evolution of resistance in one or more cell populations that were cultured in the incubator 2 using the results of the fitting. The processor 101 may execute instructions to identify one or more follow up experiments to be performed on the cell populations in the incubator based on the identified likely mechanism of evolution of resistance. The processor 101 may execute instructions to display any of the results of any of the steps implemented by the processor to a user through the user interface 103.***

[0240] The features disclosed in the foregoing description, or in the following claims, or in the accompanying drawings, expressed in their specific forms or in terms of a means for performing the disclosed function, or a method or process for obtaining the disclosed results, as appropriate, may, separately, or in any combination of such features, be utilised for realising the invention in diverse forms thereof.

[0241] While the invention has been described in conjunction with the exemplary embodiments described above, many equivalent modifications and variations will be apparent to those skilled in the art when given this disclosure. Accordingly, the exemplary embodiments of the invention set forth above are considered to be illustrative and not limiting. Various changes to the described embodiments may be made without departing from the spirit and scope of the invention.

[0242] For the avoidance of any doubt, any theoretical explanations provided herein are provided for the purposes of improving the understanding of a reader. The inventors do not wish to be bound by any of these theoretical explanations.

[0243] Any section headings used herein are for organizational purposes only and are not to be construed as limiting the subject matter described.

[0244] Throughout this specification, including the claims which follow, unless the context requires otherwise, the word “comprise” and “include”, and variations such as “comprises”, “comprising”, and “including” will be understood to imply the inclusion of a stated integer or step or group of integers or steps but not the exclusion of any other integer or step or group of integers or steps.

[0245] It must be noted that, as used in the specification and the appended claims, the singular forms “a,” “an,” and “the” include plural referents unless the context clearly dictates otherwise. Ranges may be expressed herein as from “about” one particular value, and / or to “about” another particular value. When such a range is expressed, another embodiment includes from the one particular value and / or to the other particular value. Similarly, when values are expressed as approximations, by the use of the antecedent “about,” it will be understood that the particular value forms another embodiment. The term “about” in relation to a numerical value is optional and means for example + / - 10%.

[0246] Examples

[0247] INTRODUCTION - EXAMPLES 1-7: INVESTIGATING THE EVOLUTIONARY DYNAMICS OF DRUG RESISTANCE IN COLORECTAL CANCER

[0248] Cancer treatment frequently fails due to the evolution of drug-resistant cell phenotypes caused by underlying genetic or non-genetic changes. The origin of these adaptations, their timing and rate of spread is key information for distinguishing the mechanism(s) of drug resistance, yet the dynamics cannot be observed directly. In the work described in the following examples, the inventors construct a mathematical framework to infer the dynamics of drug resistance without the need for direct measurement of the resistance phenotype using only genetic lineage tracing and population size data. Extensive simulation experiments show that the framework can recover ground-truth evolutionary dynamics from lineage tracingdata. Experimental evolution to 5-Fu chemotherapy in two common colorectal cancer cell lines SW620 and HCT 116 provides empirical demonstration of the veracity of the framework. In SW620 cells, a stable preexisting resistant subpopulation was inferred, whereas in HCT116 cells resistance emerged through phenotypic switching into a slow growing resistant state with stochastic exiting into a fully resistant phenotype. Extensive functional assays, including scRNA-seq and scDNA-seq validate the distinct evolutionary routes and their molecular nature. The mathematical framework developed by the inventors can be extended to diverse experimental designs to infer the evolutionary dynamics of cancer cell therapy resistance evolution from readily obtained experimental data, enabling more rapid characterisation of resistance mechanisms.

[0249] METHODS - EXAMPLES 1-7: SIMULATIONSAND MODEL FITTING

[0250] Summary. To model resistance evolution, the inventors employed two separate simulation approaches. In the first, they designed a hybrid model that simulates the total population sizes in distinct phenotypic compartments during periodic treatment. This model switched between a stochastic (jump process) and deterministic (ordinary differential equation) model when a population threshold was crossed to maintain stochastic dynamics experienced at low population sizes. Depending on the model in question, the simulation tracks 2 (unidirectional (A) and bi-directional switching (B)) or 3 (escape transitions (C)) phenotypes that differ in their behaviour when exposed to treatment. However, this simulation is restricted to generating the total population of each phenotype and does not record lineage relationships.

[0251] The second simulation approach was a fully stochastic agent-based model simulated via a rejection-kinetic Monte Carlo algorithm. Here, individual cells were assigned lineage identities at the start of the simulation and the model output therefore included total population sizes per phenotype (as in the hybrid model) but also lineage size distributions, where each lineage was assigned at the beginning of each simulation. The inventors adopted a Bayesian approach for parameter inference: specifically, the more computational expedient hybrid model was used for the initial generations of approximate Bayesian computation (ABC), where observed data were compared to simulated population trajectories. Later generations simulated lineage distributions with the fully stochastic model and the distance was computed between observed and simulated lineage diversity statistics to generate the final posterior distribution. All computational simulations were written in Julia (version 1.7.2) (Bezanson et al., 2017).

[0252] Modelling approach. The inventors developed quantitative models that captured cell growth and death dynamics and phenotypic evolution in a population of cancer cells exposed to treatment that alters growth, death or phenotype dynamic parameters. Total cell population size trajectories and lineage information (as summarised by diversity statistics) were used to infer the model parameters that dictated a population’s response to treatment and to distinguish between different evolutionary models.

[0253] A birth-death model was constructed, where the birth and death rates of cells govern whether the population expands or contracts over time. Differences in phenotypic behaviours were encoded in the different birth and death rates cells experience during the periods on or off treatment. The treatment windows were predefined and controlled the effective concentration of the drug treatment overtime.

[0254] Other aspects of phenotype change were also described. For example, resistant cells can experience a fitness deficit or ‘cost’ in the absence of treatment, relative to sensitive cells. The rate at which cells cantransition between these phenotypes is also a key determinant of the population dynamics. At extremely low rates, transitions could be the product of rare (epi-)genetic mutations. Conversely, at higher values, more rapid transitions between phenotypes could emerge due to non-genetic phenomena, such as the stochastic partitioning of gene transcripts during cell divisions. To capture these different scenarios, transitions between phenotypic compartments with differing probabilities were modelled, with the number of different possible phenotypic transitions depending on the complexity of the specific model.

[0255] During the growth and treatment of cancer cells, population sizes are expected to vary across several orders of magnitude. At large cell population sizes, the random birth and death events average out and the system can be approximated deterministically. However, at small population sizes, dynamics are stochastic and deterministic approximations can fail to account for phenomena such as extinction or the stochastic emergence of phenotype populations following a rare transition event. These dynamics can be modelled with discrete time models using stochastic simulation algorithms. However, as the number of stochastic transitions per time step increases (e.g. when the population is large and the number of birth events are a function of the total population size), these approaches become computationally expensive, posing a barrier to parameter inference in simulation-dependent, likelihood-free approaches. As such, there have been attempts to implement hybrid approaches that adaptively switch from stochastic methods at low population sizes to deterministic approximations at high population sizes (Germane et al. 2024; Kreger, Komarova, and Wodarz 2021 ). In the work described in the present example, the inventors implemented a hybrid model where each phenotypic compartment is modelled as a stochastic jump process at small population sizes but a deterministic ordinary differential equation (ODE) at large population sizes.

[0256] The inventors also investigated the behaviour of individual cell lineages when exposed to treatment. Previous work has shown that the evolutionary dynamics experienced by a population are encoded in the lineage relationship between cells (Williams et al. 2018; Bollen et al. 2021; Karlsson et al. 2023; Gabbutt et al. 2022). These relationships can be measured retrospectively using heritable (epi)genetic marks, or prospectively using barcode-based approaches that incorporate a traceable lineage tag into cells (Kebschull and Zador 2018). Lineage tracing is easily implemented in silico: agent-based models track each individual cell in a system and can assign unique cell ‘barcodes’ that are inherited during cell divisions. These models can also easily accommodate different phenotypes by assigning and tracking the phenotype of each cell ‘agent’. However, modelling individual cell lineage evolution requires all lineages to be represented in a model, which precludes the use of (deterministic) approximations to average population dynamics.

[0257] Two modelling frameworks were therefore developed in parallel: a hybrid compartment model that tracks the total number of cells per phenotypic compartment and adaptively switches between an ODE and stochastic jump process, and a purely stochastic agent-based model that tracks both the phenotype and lineage relationship of all cells in the system.

[0258] Base Model. The base model described herein comprises features which are shared by all models. The base model includes two cell compartments demarcated by the distinct phenotype of cells that occupy each compartment: sensitive and resistant to drug treatment, with cell population sizes ns(t) and nR(t) at time t, respectively. The total number of cells at time t was denoted by N(f). Cells divided and died with rates specific to their phenotype: sensitive cells divide with rate bs and die with rate ds', resistant cells divide with rate 5R and die with rate C / R. The parameter 5 controlled a fitness deficit experienced by resistant cellsrelative to their sensitive counterparts, where resistant cells divide with rate bR= bs(l - 5) and die with rate dR= ds(l - S') (Fig. 4A). This formulation led to the interpretation of 5 as the relative reduction in net growth rate such that 2R=

[0259]

[0260] (1 - S), where A = (b - d). The inventors note that there are multiple possible modifications to bR and dRthat could lead to the same net fitness deficit. In particular, when a fitness deficit is included (which is not a requirement of the methods of the present disclosure, as cells that are resistant to the drug can also be assumed to not incur any fitness penalty in the absence of the drug, essentially modelling an escaped cell population instead of a resistant cell population, see below) it can be parameterised in different ways. For example, the net growth rate can be reduced by decreasing the birth rate, increasing the death rate, decreasing the birth rate and increasing the death rate, or decreasing both (as was done here). The specific combination chosen here reflected an assumption of the present model, which is the chosen parameterization in which a single parameter 5 denotes a reduction in both birth and death rates.

[0261] The inventors aimed to fit their models to long-term evolution experimental data, where the carrying capacity of the system is often known. For example, the maximum number of cells that a cell culture flask can hold is easily measurable. The carrying capacity was denoted as K cells. Growth was assumed to be logistic. Motivated by empirical observations that cell populations that approach the experimental carrying capacity enter a state of stasis where both cell birth and death rates are reduced (including treatment-mediated cell death), logistic growth was modelled by scaling the birth and death rates by the density dependence term (1 - -). This formulation led to a reduction in both birth and death rates as the total population size approaches the carrying capacity, K. This led to the pair of differential equations shown in Fig. 4B and reproduced below:

[0262] dn, A N\

[0263] -T- = (.bsns- dsns) [l --\ (1) at \ KJ dnR / N\ — = (bRnR- dRnR) (1 - -j (2)

[0264]

[0265] The inventors note that there are multiple ways to partition the density dependent term amongst the birth and death rates to produce net growth which is logistic (Huynh, Scott, and Thomas 2023) and that the true partitioning is often unknown. In other words, density dependent growth, including specifically logistic growth as used in the present examples, can be partitioned between the birth and death rates in any way that gives: dN / dt = r (1 - N / K) N (or a corresponding model for another density dependent growth model) where r is the net growth rate r=(b - d). Decreasing both the birth and death rates of cells is an assumption of the formulation of the model used here, but other implementations may only reduce the birth rate or increase the death rate, or increase the death rate and reduce the birth rates to different extents that together results in the assumed density dependent growth effect.

[0266] One common question facing investigations into resistance is whether resistant cells exist prior to treatment and if they do, what is the fraction of resistant cells. In the base model (and all derived models), the parameter p sets the initial conditions by controlling the proportion of resistant cells when t = 0, such thatPNW ■Modelling Treatment. To model the impact of treatment on the population dynamics, an assumption was made that the addition of drug increases the death rate according to a cell’s phenotype (Fig. 5A). This increase is a function of the effective drug concentration at time t, termed D(t). Based on the experimental design or treatment schedule, the simulation time is split into pre-defined ‘treatment on’ and 'treatment off’ windows, which begin at Ton and Toff, respectively (Fig. 5B).

[0267] Motivated by experimental observations that a lag time exists between the administration I removal of treatment and the cytotoxic effect of treatment on cancer cells, the change in drug effect was modelled as an uptake / decay model where

[0268] £>(t) = Dcy(t) (3)

[0269] -Y _ fKtE Ton

[0270] dt ~ l -K t E ToffWwhere Dcis the maximum effective concentration of the drug, y(t) is a dimensionless function representing the drug fraction over time, and K is a rate constant related to the drug uptake / decay dynamics. To represent the fraction of effective drug, y(t) is constrained to the interval [0, 1]. Unlike other pharmacokinetic models, the inventors only aimed to model the effective drug concentration over time, realised as the cytotoxic effect exerted on cells. For example, this could be the combined effect of a drug's metabolic stability and the activity of a cancer cell’s efflux pumps. The inventors assumed the rate at which this effect accumulates in cells after treatment commencement (t G Ton) was the same at which it diminishes when treatment ends (t G TOff) K. In sensitive cells, the effect of treatment was described as the increase in death rate such that ds(t) = ds+ D(t) (5) Whilst resistant cells were expected to have lower levels of drug-induced cytotoxicity, resistance might not be complete. The difference between drug-induced death experienced by resistance cells relative to sensitive cells is controlled by the parameter qj such that

[0271] dR(t) = dR+ (D(t) (1 - »))< 0 <. < 1 (6)

[0272] As such, in the models qj = 1.0 corresponds to cases where resistance is ‘complete’ and resistant cells are unaffected by treatment, whereas qj = 0.0 corresponds to cases where sensitive and resistant cells experience the same elevated death rate in the presence of treatment.

[0273] Hybrid model switching. To model the evolution of cell phenotype compartments over time, deterministic ODEs provide a computationally expedient approximation when populations are large. For the base model, this led to the pair of differential equations shown in Fig. 4B (left). However, deterministic ODEs fail to account for stochastic dynamics when cell population sizes are small. In the phenotypic compartment model described here, this might be the case if one of the two phenotypes is rare in the population, or when treatment drives both phenotypes to smaller cell numbers. To capture these dynamics, the system can be instead modelled as a discrete time stochastic jump process, shown in Fig. 4B (right), where birth and death events lead to an instantaneous, discrete change in the population size. The drawback of the stochastic simulation algorithms (SSAs) used in these models is their computational cost at large population sizes. One solution is to implement a hybrid approach where the model uses an SSA for small populationsizes of a given compartment, but switches to a deterministic ODE system when a population threshold has been crossed (Germane et al. 2024; Kreger, Komarova, and Wodarz 2021). The inventors employed a hybrid switching model where phenotype compartments are modelled as stochastic jump processes when nx(t) <= A / switch but switches to a deterministic ODE model when nx(t) > A / switch, for the phenotypic compartment x (Fig. 4B) (see Choosing NSWitch). This approach retains the stochastic forces experienced by phenotypic compartments when the cell population sizes are small whilst also retaining computational efficiency at large cell population sizes.

[0274] Model A: unidirectional transitions. The base model above describes the growth dynamics of sensitive and resistant populations during periods on / off treatment. To investigate evolutionary scenarios where cells could transition between phenotypic compartments, model A permitted switching between phenotypic compartments. Model A permits a single direction of switching, where sensitive cells can transition to resistant cells with probability p per cell division (Fig.5C).

[0275] As such, the ODEs controlling the change in each phenotypic compartment were:

[0276] dnsN

[0277] — = (nsbs- nsds- nsbs\i) ■ (1 - -) (7) dnHN — = (nRbR- nRdR+ nshsp) ■ (1 - -) (8)

[0278]

[0279] and the corresponding transition rates for the stochastic jump process were:

[0280] S — » S + S, bs

[0281] s 0, ds

[0282] S -» S + R, bsii (9)

[0283] R —> R + R, bR

[0284] R ~ * 0, dR

[0285] The stochastic jump process transitions were represented by intrinsic rates without accounting for logistic growth. When p is very low, the model can produce behaviours typical of genetic mutations that control a resistant phenotype, whereas high values can resemble rapid, non-genetic transitions into a resistant phenotype.

[0286] Model B: bidirectional transitions. Given observations that cancer cells can change phenotype rapidly (Emert et al. 2021; Hi-nohara et al. 2018; Rehman et al. 2021), model B was constructed to explore scenarios where these transitions could occur both rapidly and reversibly. Model B included an additional parameter that denotes the probability of a back transition: a (Fig. 5D).

[0287] The phenotypic compartment ODEs were:

[0288] dnsN

[0289] — = (nsbs- nsds- nsbsi + nRbRo) ■ (1 - -) (10) at K

[0290] dnRN — = (nRbR- nRdR+ nshsp - nRbRa) ■ (1 - -) (11)

[0291]

[0292] dt K and the corresponding transition rates for the stochastic jump process were:

[0293] S S + S, bs

[0294] (12) dsS -» S + R, bs\

[0295] R —> R + R, bR

[0296] R - * 0, dR

[0297] R -» R + S, bRo

[0298] This model can capture behaviours typical of both forwards and backwards mutations, as well as rapid, non-genetic, bi-directional phenotypic switching between sensitive and resistant phenotypes. The modified birth and death rates during treatment periods remained as described in model A (see Modelling Treatment).

[0299] Model C: escape transitions. Model C captures the dynamics of resistance evolution whereby cells can exist in a slow growing phenotype that is refractory to treatment before possibly entering a third phenotypic state which is resistant to treatment but does not exhibit a reduced proliferation rate, thereby escaping the fitness penalty associated with the resistant phenotype (‘escape transitions’).

[0300] An additional phenotype (‘escape’) was incorporated that performed identically to the previously described resistant phenotype, however did not incur the fitness penalty described by 5 (Fig. 5E). In addition to the transitions permitted in Model B, resistant cells could also transition to the escape phenotype with probability (a ■ y(t)) per cell division. The parameter y e [0, 1] controlled the relative drug concentration and y(t) was determined by the treatment timings and rate constant K (Fig. 5B). The transition from resistant to escape was therefore assumed to be ‘drug-induced’, similar to previous observations of drug-induced changes to mutation rates in cancer (Russo, Crisafulli, et al. 2019; Russo, Pompei, et al. 2022; Pisco et al.

[0301] 2013).

[0302] With the addition of the new ‘escape’ phenotype, the compartment ODEs were:

[0303] N

[0304] — = (nsbs- nsds- nsbs\i + nRbRo) ■ (1 - -) (13) dnRAf — = (nRbR- nRdR+ nsbs\i - nRbRo - nRbR(ay(tyj) ■ (1 - -) (14)

[0305] dnE,, N — = (nEbE- nEdE+ ■ (1 - -) (15)

[0306]

[0307] and the corresponding transition rates for the stochastic jump process were:

[0308] S — » S + S, bs

[0309] S -> 0, ds

[0310] S -> S

[0311]

[0312] + R,

[0313] R R + R, bR

[0314] R ^ 0, dR(16)

[0315] R -» R + S, bRa

[0316] R -> R + E, bRay(t)

[0317] E — > E + E, bE

[0318] E -» 0, dE

[0319] Given that resistant cells can now transition to either the sensitive or escape phenotype with probabilities o and a (where ay(t) < a), an additional constraint was included: ct + a < 1.0. This was introduced because both parameters ct and a are assumed to represent the probability of a transition during cell divisions.Alternative implementations could set these rates as rates rather than probabilities of transitions during cell divisions, obviating the need for the additional constraint.

[0320] The modified death rates during treatment for each phenotypic compartment were:

[0321] ds— ds+ D (t) (17)

[0322] dR= dR+ (D(tMi -ipy) (18)

[0323] dE= dE+ (D(t)(l - »)) (19) This model formulation assumes that the resistant and escaped population differ only by the fitness penalty in the absence of treatment, but that the effect of drugs on their death rates is the same. However, an alternative to this formulation would be to allow the parameter i to vary between the different resistant phenotypic compartments (i.e. resistant and escaped cells, in the present example). Therefore, equations (18) and (19), while written with a single i for simplicity, can also be expressed with respective ipR(equation (18)) and ipE(equation (19)).

[0324] Model Assumptions. The models outlined above rely on the following assumptions:

[0325] • Growth is logistic, realised as a decrease in both birth and death events as the population approached carrying capacity.

[0326] • Cell phenotypes can be described via discrete 'phenotypic compartments’.

[0327] • All cells in a phenotypic compartment have uniform birth and death rates.

[0328] • Phenotypic transitions are uniform for a given phenotype and are coupled to birth events.

[0329] • Transitions from the resistant to escape phenotype are drug-dependent.

[0330] • Any fitness penalty associated with resistance in the untreated environment (relative to the sensitive phenotype) is realised as a combination of lower birth and death rates, i.e. the fitness penalty is associated with a reduction in net growth rate. Note that this is juts one possible way of capturing a fitness penalty, which was designed to capture 'quiescent' dynamics where cells that escape treatment are slow cycling. This was thought to be best captured by a combined decrease in both the birth and death rates. However, a simpler fitness penalty could be represented as one that simply lowers the birth rate of the cells in the compartment that experiences the fitness penalty.

[0331] • The effect of treatment on cells is realised exclusively as an elevated death rate.

[0332] Agent-based lineage model. An agent-based stochastic branching process was designed that simulates the same process as the hybrid phenotypic compartment model, with the primary distinctions being that all individual cells are tracked along with their phenotype and lineage identity and the simulation is fully stochastic (Fig. 13).

[0333] Each cell was characterised by a unique set of attributes:

[0334] 1. Phenotype: As in the phenotypic compartment model (see Hybrid model switching), each cell can exhibit one of two or three phenotypes: S (Sensitive), R (Resistant), E (Escape), only in model C).

[0335] 2. Lineage Tag: Each cell is assigned a unique lineage tag i at the beginning of the experiment, where i = 1, 2, 3 No (No is the number of uniquely barcoded cells when t = 0).Therefore, for the agent-based lineage model, at any given time the model tracks the following quantities:

[0336] • The number of cells n with lineage tag i that have phenotype ‘S’.

[0337] • The number of cells n with lineage tag i that have phenotype ‘R’.

[0338] • The number of cells n with lineage tag i that have phenotype ‘E’.

[0339] When modelling long-term experimental evolution (see below, Modelling long-term experimental evolution), given the introduction of replicate sub-populations, the agent-based lineage model tracked the following quantities:

[0340] • The number of cells nzwith lineage tag i that have phenotype ‘S’ in replicate z.

[0341] • The number of cells nzwith lineage tag i that have phenotype ‘R’ in replicate z.

[0342] • The number of cells nzwith lineage tag i that have phenotype ‘E’ in replicate z.

[0343] Modelling long-term experimental evolution. To investigate the evolutionary dynamics of resistance using data collected from long-term evolution experiments, the inventors implemented the sampling of cells and treatment steps into the models to mimic common features of these experiments (Fig. 6). The simulations begin with No cells. Cells are then grown for a fixed period of time (texp) in the absence of drug treatment, producing a population of Nexpexpanded cells. Importantly, cells grow and transition between phenotypic compartments during this stage in the untreated environment according to the simulation’s evolutionary parameters. Cells were sampled without replacement into separate replicate sub-populations. For each replicate, denoted by z (where z = 1, 2,..., Z), Nseed cells were sampled, resulting in a total ofszeedcells per replicate. Each of the Z replicate sub-populations were exposed to periodic treatment as previously described (see Modelling treatment). When a pre-defined ’Passage time’ (tpass; or maximum population size (Nmax) (whichever occurred first) was reached, each of the Z replicate sub-populations were sampled without replacement for a second time and growth / treatment continued until either tmaxor Nmaxis reached. This step can be repeated multiple times to mimic the ‘Passaging’ steps common to in vitro experiments.

[0344] Parameter inference and model selection - overview. Due to the stochastic dynamics exhibited by the models, the resulting behaviours could not easily be described by probability distributions and calculating the likelihood was impractical. In such cases, ‘likelihood-free’ methods may be used, such as 'Approximate Bayesian Computation’ (ABC), that rely on comparing simulated data to observed data using a distance measure. By repeatedly comparing various input parameter values, as sampled from a prior distribution, a posterior distribution of parameters that produce simulations which closely resemble the empirical data can be generated and the distance measure can be minimised, thereby closely approximating the true posterior (Toni, Welch, et al. 2009). To infer parameters of interest across the three models of resistance evolution, the inventors implemented an ABC approach.

[0345] Two separate sources of information were used from the cancer cells exposed to periodic treatment to infer parameters of interest. These measurements were chosen to also capture features of the population commonly measured in long-term resistance evolution experiments:

[0346] the total number of cells in each replicate sub-population at Observation (O) and Passage (P) times, and

[0347] the cell lineage distribution of each replicate population at each Passage (P).Since measuring the total population size in these experiments can be performed more frequently than measuring the cell lineage distributions, which requires harvesting the cells, extracting their DNA and sequencing the lineage tags / barcodes, lineage measurements were only made at the ‘Passage’ events (Fig. 6A). At each passage event, a pre-determined number of cells are seeded into the subsequent flask. This value is fixed and determined by the known cell number seeded in the experimental design. The remainder are used for lineage tracing.

[0348] ABC relies on extensive simulations, but tracking each individual cell in large cancer populations (> 108cells) across many iterations is computationally prohibitive. To address this, a simplified model was first used, which tracked only the total size of each phenotypic compartment, allowing unlikely parameter values to be discarded. Then, the fully stochastic agent-based models were applied to infer parameters using ABC with simulated cell lineage distributions. In ABC technically one can only say that the final generation is an ‘approximation’ of the posterior. In theory, as Et — >0 and with sufficient statistics, ABC can approach the true posterior, but achieving this exactly is too computationally expensive in most practical applications. Therefore, the posterior distribution that is estimated is referred to as an “approximate” posterior distribution.

[0349] Parameter inference and model selection - step 1: population sizes with the phenotypic compartment model. The computationally cheaper phenotypic compartment model was first used to perform ABC-SMC, using the Python package pyABC (Schalte et al. 2022). Parameters were sampled from a prior distribution TT(0) to generate synthetic data x, which was then compared to observed data y using a distance function d(x, y). Parameters that produced synthetic data within a tolerance Et were retained. This tolerance Et decreased in subsequent iterations, tightening the acceptance criterion until a target tolerance was reached, thereby approximating the posterior distribution

[0350]

[0351] |y) (Toni and Stumpf 2010). In the present examples, the inventors used the pyABC package, which includes a built-in decreasing scheme for tolerance Et However, any implementation of the ABC algorithm may be used, as well as any tolerance Et updating strategy known in the art (including but not limited to those implemented in the pyABC package). The target tolerance was set by comparing the quality of parameter recovery in synthetic data and striking a balance between accuracy and computational run time.

[0352] For the first step, observed population size estimates (Oi, O2,..., Om) at predefined observation times {t01,t02, -^Om) and passage population size estimates (Pi, P2,..., Pn) at predefined passage times (tP1,tp2, -,tPn) were used These estimates were combined into a dataset of length m + n, consisting of matching times t and population sizes N. These datasets were generated for all Z replicates, resulting in a total of Z(m + n) observations {t, N}. The distance between observed and simulated data was calculated as the mean of the Euclidean distances between all {to, No} and {tp, Np} positions in this dataset (Fig. 14B).

[0353] Note that although the Euclidian distance was sued as a distance metric in these examples, any distance metric may be used that can quantify how “close” observed statistics are to simulated values. pyABC uses an ABC Sequential Monte Carlo (ABC-SMC) algorithm - sampled parameter values are propagated through generations, where each generation a sample is weighted according to the previous generation’s weight (this incorporates information from the prior distribution) and subject to a perturbation kernel to generate the next generations proposal distribution. This leads to faster convergence than a simple rejection algorithm (which was used in the second step below). However, the approach would still work with a simpler algorithm (such as a simple rejection algorithm as used in step 2), it would just be moderately slower.Parameter inference and model selection - step 2: lineage distributions with the agent-based model. For the second step of the inventors' ABC approach, parameter estimates were sampled from the latest pyABC SMC-ABC generation and the stochastic agent-based lineage model was simulated using these parameter values (i.e. exact values from the final generation of step 1 are sampled and used to simulate the full lineage version of the model). For each replicate (1, 2, Z), at each passage time point (Pi, P2, Pn), a lineage distribution was generated corresponding to the counts of lineages (1, 2 No) (Fig. 14C). These distributions were then converted into Z x n two-dimensional lineage diversity statistics (see Lineage diversity statistics). The diversity statistics were normalised to between 0 and 1 and the Euclidean distance between these diversity statistics and the observed data was calculated. The closest 10% of simulations were retained to generate the final posterior distribution (Fig. 14D). The 10% threshold was selected to balance accuracy and computational time using synthetic data - higher values often lead to worse parameter estimates with synthetic data, whereas lower values mean more simulations have to be run to retain the same number of samples to generate the final posterior distribution. The final posterior distribution is obtained based on the count of the number of samples that passed that target distance criterion in each region of the parameter space that was explored.

[0354] Parameter inference - prior distributions. To capture a broad number of possible scenarios concerning the level of pre-existing resistance and stability of the phenotypes of interest, uniform priors were set for the parameters controlling the initial conditions and the phenotypic transition probabilities such that each parameter was distributed in negative logarithmic space:

[0355] — loglO p) ~ Uniform(0, 7) (20)

[0356] — Zo^lO(.) ~ Uniform(0,9) (21)

[0357] —log 10 (cP) ~ Uniform(0, 9) (22)

[0358] — loglO(a) ~ Uniform(0,9~) (23)

[0359] Because No = 106, a value of p < 10~6corresponded to no resistant cells as an initial condition. The minimum probabilities in the three phenotypic transition rates (10~9) are all lower than previously explored mutation rates in cancer resistance evolution studies (Michor et al. 2005; Chmielecki et al. 2011; Diaz Jr et al. 2012). Uniform priors were also set on the fitness penalty associated with the resistant phenotype in negative logarithmic space:

[0360] — Iogl0(8~) ~ Uniform(0,3') (24)

[0361] This led to high prior probability for low fitness penalties; the evidence in favour of higher penalties (5 > 0.1 ) must be large in the inference to confidently detect a 'cost' of resistance.

[0362] The above prior were chosen to be flat (i.e. relatively uninformative) priors that are just bounded to biologically realistic values. The choice of priors may depend slightly on the experimental design and system under investigation. For example, p determines the proportion of cells at the time of barcoding, so the lower bound should be guided by the maximum number of cells barcoded. As another example, the transition rates might include more informed priors if previous evidence suggests phenotypic transitions between compartments are rare. As another example, when a fitness penalty is suspected for biologicalreasons, a non-uniform prior for 6 may be used, such as e.g. a beta prior (which is naturally bounded between 0 and 1 and can therefore be applied on the untransformed parameter space). The skilled person would be able to choose prior distributions that represent a given set of biological assumptions based on the teaching of the present disclosure.

[0363] Cases where the ‘resistant’ phenotype (if present) confers some protection against treatment are of particular interest. In the presently described model, this is controlled by the ‘strength of resistance’ parameter i, where i = 1.0 denotes complete resistance. Therefore, a distribution which favours values closer to 1.0 was used for the prior distribution of ip, resulting in a mean of 0.75 and a 95% confidence interval of approximately (0.38, 0.97):

[0364] ip ~ / ?(3.0, 1.0) (25)

[0365] Note that any prior distribution bounded between 0 and 1 with probability mass skewed to 1 may be used instead of a beta distribution. In embodiments, a separate short drug exposure assay may be used to help better calibrate the ip prior distribution and the proportion of resistance (p) prior distribution.

[0366] Parameters that controlled the behaviour of the effective drug concentration (effective drug concentration (De) and rate of drug accumulation / decay (kappa)) were constrained to positive values with broad T distributions:

[0367] De ~ T(3.0, 1.0) (26)

[0368] K ~ T(2.0, 1.0) (27)

[0369] For large values of K, the rate of change of the effective drug concentration was high enough that the drug-induced death became effectively instantaneous. The gamma distribution was chosen because it is strictly positive and flexible. Further, the gamma distribution is an appropriate distribution to model rates. However, alternative distributions may be used, such as e.g. exponential or log-normal prior distributions.

[0370] Model selection. For model selection, the deviance information criterion (DIC) was calculated, which is a measure of model fit that penalises model complexity (Francois and Laval 2011) by simulating from the approximate posterior distribution. DIC was calculated as follows:

[0371] DIC = 2D — D(0~) (28)

[0372] where

[0373] n

[0374] D = - 21og (- (V KE(sl- s0)) (29)

[0375] n _ <

[0376] i=l

[0377] given n replicates from the posterior distribution (i.e. n independent draws from the approximate posterior distribution, in other words n replicates each using a set of parameters sampled from the approximate posterior distribution obtained at step 2 above), where s' were the summary statistics obtained from the ithsample from the posterior predictive distribution (i = 1,..., n) and so were the same summary statistics given the observed data, and

[0378] n

[0379] D(0) = -21og KE(sl- s0)) (30)

[0380] n Z_i

[0381]

[0382] i=lwhere n replicates were instead simulated using a point estimate (the mean) of the posterior predictive distribution (i.e. taking the mean value of each parameter’s posterior distribution, and then simulating n times using those same values repeatedly), and

[0383] 1 u2

[0384] KE(U) = -=e~,u E R (31)

[0385]

[0386] ■V2TT

[0387] where KE(u) is a Gaussian error kernel which is used to model the distribution of distances around the observed summary statistics. This is a common kernel adopted when calculating the DIC.

[0388] Choosing NSWitch- In the phenotypic compartment model, the stochastic jump process is used when nx< Nswitch and the deterministic ODE is used when nxs Nswitch for phenotypic compartment x. The chosen value of Nswitch controls the trade-off between properly accounting for stochastic dynamics at low population sizes and computational efficiency at large population sizes. To assess how changing this switching threshold influenced the ability to recover true parameter values from synthetic data, the estimates given by different values of Nswitch were compared for the phenotypic compartment model.

[0389] Ground-truth data was simulated using the fully stochastic lineage model (see Agent-based lineage model) before performing the first step of the parameter inference process that fits the phenotypic compartment model to the total population size observations (see Parameter inference - Step 1: population sizes with the phenotypic compartment model) with Nswitch = [0, 1, 10, 100, 200, 500] (Fig. 15). By Nswitch = 10, the inference process could successfully identify a region within the parameter space that encompassed the true value. For some simulations where Nswitch = 500, the ABC failed to complete in the allocated time (the pyABC step had a maximum run time of 168hrs).

[0390] Nswitch was therefore set to 200 for all simulations. For validation steps that used synthetic data, ‘ground truth’ observed data was generated using the fully stochastic agent-based lineage model (see Agent-based lineage model), ensuring that assessment of the framework was not biased by the hybrid model approach and the chosen value of NSWitch-

[0391] Lineage diversity statistics. The population sizes simulated using the modelling framework meant that lineage relationships could be tracked in cell population sizes of the order 108. After also incorporating replicate sub-populations and multiple time-points, the lineage information produced from a single simulation could consist of millions of observations. To enable the interpretation of these distributions, summary statistics that captured salient features of the data were chosen.

[0392] Exploration of the models indicated that there were two important axes of information in the lineage distributions: first, the number of cell lineages that survive within each of the replicate sub-populations, and second, the relationship of the cell lineages when comparing between the replicate sub-populations. The results suggested that behaviours such as the stability of a resistant phenotype through time and difference in fitness between phenotypes could dictate these relationships following a period of selection (via treatment). To capture these two axes, the ‘Hill Diversity Indices’ were used to measure the lineage diversity within a population (Jasinska et al. 2020; Jost 2006; Roswell, Dushoff, and Winfree 2020). These diversity indices were also extended to measure the diversity between populations (and / or experimental replicates). Advantageously, Hill Diversity Indices capture a variety of diversity indexes used to model populations (e.g. Shannon, Simpson etc.) by simply modifying the q parameter (see below). Therefore, the use of Hilldiversity indices encompasses use of any indices that are special cases of the Hill diversity indices and is flexible enough to use whichever index a user may want to use, by setting the value of q.

[0393] Lineage diversity statistics - within-population diversity. To capture the diversity of lineages within a subpopulation, Hill diversity indices of o (rder q,qD, were used, defined as:

[0394] j > l / (l-q)

[0395] (32)

[0396]

[0397] 7=1 /

[0398] where pj > 0 was the relative frequency of the jthlineage, and q controls the weight that high frequency lineages contribute to the index. Relative frequencies are scaled according to q, where values of q < 1 and q > 1 preferentially leverage low- and high-frequency lineages, respectively. When q = 0, q=0D is simply the total number of lineages (often referred to as 'species richness’ in ecology). The exact solution for Equation (30) when q = 1 does not exist, however as q — 1,

[0399] j

[0400] q=1D = exp (- pj log (p7)) (33)

[0401]

[0402] 7 = 1

[0403] which is also referred to as the Shannon diversity (the natural exponential of the Shannon entropy). When all lineages are present in equal proportions (pi = p2 =... = pj) the values ofqD are insensitive to q and correspond to the total number of lineages.

[0404] The ability to over-leverage high-frequency lineages when calculating a diversity metric has two advantages. Firstly, the behaviour of successful lineages that survive treatment are of primary interest. Secondly, deriving lineage information of cells experimentally is noisy: errors in lineage tag amplification and sequencing can inflate the number of low-frequency lineages, giving a false impression of the true number. Higher values of q allow the weight with which putative low frequency noise lineages contribute to the diversity metrics to be decreased.

[0405] Lineage diversity statistics - between-population diversity dissimilarity. Hill-diversity indices of the order q allowed for the tuning of how strongly high-frequency lineages influenced diversity statistics. Lineage frequencies may also be transformed (xq) and back-transformed (x1 / <1~'i)) when quantifying the dissimilarity in diversity between populations. Specifically, this is captured inqD(0), the effective number of unique populations (or the true diversity of order q).qD(0) is the equivalent of the beta diversity, with the distinction that these transformations permitted adjustment of the contribution of high-frequency lineages, whilst also having an intuitive range of [1, I], where I is the total number of populations being compared (Jost 2006; Roswell, Dushoff, and Winfree 2020).qD(0) can be calculated by partitioning multiple sub-population’s diversity into two components:

[0406] qD(y)

[0407] = (34)

[0408] where (j,

[0409] Pq, j (35)

[0410]

[0411] 7 / =l /

[0412] where I is the total number of sub-populations being compared, Ji is the total numberof lineages amongst all I populations, and pjtis the relative frequency of lineage j amongst all I populations, and ( i J x Vfi-q)

[0413]

[0414] 7X i(2 iCPO)) / )(36)

[0415] where y is the relative frequency of lineage j in population / .

[0416] qD(y) captures the total diversity of order q when pooling all I sub-populations,qD(or) is the mean diversity of order q of all I sub-populations, andqD(jB) is the ratio between the two. When q = 0,q=0D(jB) = 1.0 when all lineages are shared amongst all I sub-populations - there is one ‘effective’ unique population; andq=0D(jB) = / when no lineages are shared amongst all I sub-populations - there are I effective unique populations.

[0417] Lineage diversity statistics - choosing q for lineage inference steps. The parameter q determines how preferentially high-frequency lineages are leveraged and directly influences the within- and between-replicate statistics used in the second step of the ABC parameter inference. The inventors explored how varying the value of q impacted the distribution of these summary statistics.

[0418] Specifically, the inventors investigated the effect of different q values on the ability to recovertrue parameter values from synthetic data where the ’ground truth’ was known. Fig. 16 (left-hand column) illustrates the relationship between within- and between-replicate diversity for a simulation using ’Model A’. When q is set to a low value (e.g., q = 0)), the statistics primarily reflect the total number of lineages without regard to their frequency. Given the dynamics are encoded in the success of a lineage as opposed to its presence / absence in a population, these statistics exhibit minimal variance between simulations. This, in turn, reduces the power to infer different evolutionary parameters.

[0419] As the value of q was increased (e.g., q = 1 and q = 2), the statistics began to capture more information from high-frequency lineages, leading to greater differentiation between evolutionary scenarios as described by varying parameter values. This improves the ability to accurately recover true parameter values, as shown in Fig. 16 (right-hand side). Further increases in q (e.g., q = 100) resulted in a situation where only the highest-frequency lineages dominated the statistics, significantly increasing the likelihood of high between-replicate dissimilarity across all parameter values simulated (measured on the y-axis as W)).

[0420] Based on these observations, Hill Diversity Indices of order q = 2 were selected when calculating the distance between observed and simulated data as part of the ABC inference framework. The failure to appropriately balance the contribution of high-frequency and low-frequency lineages can have significant consequences when dealing with real, sequenced cell lineage data. Errors introduced during PCR and sequencing can artificially inflate the number of low-frequency lineages, thereby distorting the diversity statistics if low-frequency lineages are over-leveraged. This value of q was chosen as a practical solution to balancing i) reducing the impact of low frequency lineages (we expect lots of these to be preparation / sequencing artefacts that artificially inflate the diversity) and ii) allowing moderate frequency lineages to impact the diversity statistics. In embodiments, a ‘calibration’ step can be included where the diversities are compared over a distribution of q values and this distribution is used to identify an appropriate value of q.Estimating birth and death rates with untreated samples. The net growth rate of a growing population (b -d) can be easily derived by comparing the change in population size, Nt - No, over a given time window, At. However, estimating individual birth and death rates (b and d) is more challenging. These estimates are essential for accurately simulating cell growth dynamics. Directly measuring the death rate often requires additional experiments, such as live-imaging or FACS (Johnson et al. 2019; Russo, Crisafulli, et al. 2019). Here, information from the initial expanded population of barcoded cells was used to estimate individual birth and death rates. Cells were uniquely barcoded and subjected to a shared expansion step where they grew to produce a population of Nexpcells (the same expanded pool that was used for the long-term resistance evolution experiment described in the examples below). Cells were subjected to the technical bottlenecks necessary to measure the cell lineage distributions. The observed lineage distributions are a product of the birth-death process over a known At followed by cell sub-sampling, barcode extraction, amplification, and sequencing. A Bayesian model that accounted for these bottlenecks was developed to estimate the average birth and death rates of each colorectal cancer cell line in vitro; HCT116 and SW620. These estimates can be used as a constraint in the parameterisation for the inference steps, e.g. where the values of d, b and 6 that are identified much be such that the population level observed birth and death rate (given the proportion of resistant cells at tO) in the absence of treatment is the one calculated here. This was not done in the present examples because the inventors were only interested in scenarios where resistance is rare so this initial average estimate of the rates is only slightly biased by the presence of rare resistant cells (if they exhibit a significant fitness penalty in the absence of treatment). Therefore, this was taken to represent an estimate of the birth and death rates for the sensitive population prior to exposure to treatment (i.e. bs and ds). Many other approaches can be used to infer the birth and death rates of the cell population in the untreated environment. These include e.g. live imaging, FACS sorting with live / dead cell markers and cell counting technologies that give a viable / non-viable fraction (to name a few).

[0421] Estimating birth and death rates with untreated samples - birth-death process. A schematic illustrating the sampling design of the expanded population of barcoded cells used to infer population birth and death rates is shown in Fig. 21. The probability mass function (p.m.f.) was used for the birth-death process, giving p(n), the probability of observing n cells at time t for given birth and death rates (b and d). Assuming that there is a single cell at t=0, this is given by equations (37)-(40) (see e.g. Durrett, R. 2015):

[0422] p(n) = (1 - a)(l — / ?) ■?n-1(37) for (n > 1), and

[0423] p(0) = a (38) for (n = 0), where

[0424] d(e(6-d)t _ i)

[0425] (39) b^e(.b-d)t >

[0426] (40)

[0427]

[0428] £e(&-d)t > It was assumed that each cell contained a unique barcode when the experiment began, and that therefore the distribution of barcode lineages in the Nt expanded cells were No independent realisations of the birthdeath process.The total number of cells after the mutual expansion step (Nt) was observed, which was a random variable depending on (b - d) (mean of the total population, Nt) and (b+d) (variance of the total population, cr^t).

[0429] Nt= No■ (41)

[0430] < =No ' ' ^b~d>t■ (e(b-d) t- 1) (42)

[0431]

[0432] 1b — a

[0433] where No is the number of cells at t = 0, assumed to be the number of uniquely barcoded lineages at t = 0. If barcode lineage frequencies could be observed directly, a Bayesian model could be used to infer b and d by using the pnp.m.f. for each lineage and assuming Nt ~ Normal(Nt,a^t). However, observed lineage sizes were subject to two technical bottlenecks (see sampling K cells and sampling J reads).

[0434] Estimating birth and death rates with untreated samples - sampling K cells. The first technical bottleneck was sampling a sub-population of cells before DNA extraction, where K cells were sampled from the Nt expanded pool, without replacement (Fig. 21 ). For a given lineage, because — JVt<< 1.0, this sampling step could be modelled as a Poisson distribution. The Poisson approximation allowed the sampling step to be modelled with a single parameter, simplifying the model by avoiding the need for the more complex parameters of the hypergeometric distribution. Given a lineage of size n, drawing K total cells from the Nt expanded cell pool gives the probability of sampling the lineage k times:

[0435] p(k\ri) ~ Pois(X) — — — — (43)

[0436]

[0437] K! where

[0438]

[0439] To determine the probability of seeing any barcode lineage k times after the birth-death process, the compound probability distribution was calculated:

[0440] JVtPW = ^P(fc|n) p(n) (45)

[0441]

[0442] n=0 Since lineages that are lost during sampling are not observed (p(k) = 0), this was conditioned on k > 0:

[0443] (Entop(fc|n) - p(n)) p(k\k > 0) =

[0444]

[0445] &n=oP(k> 0|n) p(n)) where

[0446] Nt p(k > 0) = 1 - p(k = 0|n) ■ p(n) (47)

[0447]

[0448] n=0 This gave a probability distribution for observing cells for a given lineage k times given K sampled cells from an expanded pool of Nt cells with birth and death rates b and d, respectively.

[0449] Estimating birth and death rates with untreated samples - sampling J reads. Another technical bottleneck is sampling reads during sequencing. After amplification and extraction, each cell’s barcode was sequencedon a flow cell. This was modelled as sampling each barcode lineage - K J times with replacement, where k was the number of times the barcode lineage appears in the K sampled cells and J is the total number of reads assigned to the sample (Fig. 21). Again, because - << 1.0, this step can be modelled as a second Poisson distribution:

[0450] Aje~* p(j\k) ~ Pois{A~) - (48)

[0451]

[0452] where

[0453] *

[0454]

[0455] = (49) Performing the nested marginalisation over both sampling distributions during parameter inference can be computationally challenging. However, the variance introduced during Poisson sampling is known and can be introduced into the p.m.f. for p(k) to account for the additional sampling variance from sequencing J total reads. This known variance was introduced into the probability distribution for observing K sampled cells.

[0456] Estimating birth and death rates with untreated samples - incorporating sequencing noise into sampled cells. The expectation and variance of a Poisson distribution where 4 = (J 7) are:

[0457] k

[0458] E[j\k] = var(j\k) = A = (-) J (50) To calculate the additional noise, the variance of the compound distribution of sampling J barcode reads from K barcoded cells was determined, assuming p(j\k~) ~ Pois - k

[0459] K - J~). The variance of a compound distribution, accounting for all probabilities of p(k), is:

[0460] var(f) = E[var( |k)] + var[E( / |k)] (51) where

[0461] K K K

[0462] var(f) = ^(yar(J\k) ■ p(k)) + (^(FQIk)2■ p(k)) - ^jE(J\E) ■ p(k)}2} (52)

[0463]

[0464] k=0 k=0 Lineages lost during sampling were not sequenced or observed in the final J reads. The compound variance conditioned on j > 0 and kwas given by:

[0465] var(j\j > 0,k~) = E\J2\j > 0, k] — (Ef I; > 0, k])2(53) Thus, the compound variance of j over all k was:

[0466] K / K, K \2\ var(f) = ^ ' (yar(j\j > 0,k) - p(k)) + I > 0,k]2■ p(E)) - E\j\j > 0, k] - p(fc) j j (54)

[0467]

[0468] where, because this applies to a lineage at size n in the expanded pool of Nt cells,

[0469] n

[0470] p(k) = p(k|n) ~ Pols (A = (— ) ■ K~) (55)This allowed for the calculation of the variance introduced when sequencing J reads from a pool of K cells sampled from an expanded pool of Nt cells with birth and death rates b and d. This extra known variance could be introduced into the probability distribution for sampling K cells by replacing the Poisson distribution with a Negative Binomial distribution.

[0471] The Negative Binomial was parameterised as:

[0472] p(k\ri) ~ NegBinom(\i, 0) = (k +J “ (— (— (56)

[0473]

[0474] \ k ' \\i + cp / \\i + cp / where the mean was:

[0475] n E[k\n] = IL = -- K (57)

[0476]

[0477] and the variance was:

[0478] u2var(k\ri) = \i + — (58)

[0479]

[0480] 0

[0481] „2

[0482] The additional variance in the Negative Binomial is —. For any n, the additional variance can be transformed into 0 using:

[0483] F

[0484] 0n=-2(59)

[0485]

[0486] var(j\ri) — p In summary, to account for the additional noise introduced during sequencing, the Poisson distribution p(k\ri)~ Pois (2 = (—

[0487] N) ■ K) was replaced with a Negative Binomial distribution, t

[0488] p(k\ri) ~ NegBinom (p. = K,<p - 0n, where <f>naccounts for the additional variance introduced during

[0489]

[0490] sequencing using Equation (59).

[0491] Estimating birth and death rates with untreated samples - recovery of rates from synthetic data. Lineage distributions were simulated using a modified version of the stochastic agent-based model (see above) with known birth and death rates. Cells were expanded for a known period of time and subjected to two rounds of sampling to simulate ‘sampling K cells’ and ‘sampling J reads’. A Bayesian framework for parameter inference was used to attempt to recover the true values, using the probability distributions described in the previous sections (see sampling K cells and sampling J reads). Bayesian statistical analyses were conducted using rstan (Stan Development Team. 2024. Stan Modelling Language Users Guide and Reference Manual, mc-stan.org), the R interface to the Stan modelling language, which facilitates full Bayesian statistical inference with MCMC sampling. Stan's programming language is specifically designed for specifying complex statistical models and employs advanced MCMC algorithms such as Hamiltonian Monte Carlo (HMC) to obtain posterior distributions. The inventors used STAN to calculate the posterior distribution (STAN uses a HMC-NUTS approach to sample from the posterior distribution). Due to correlations in the posteriors of b and d, they re-parameterised the priors such that: S_scaled ~ beta(1.0, 2.0); and L ~ gamma(1.0, 5.0) where S = (L*S_scaled); b = (L+S) / 2; d = (L- S) / 2, where L and S are reparametrised values of the birth and death rates where L is the 'cell turnover’ (L=b + d) and S is the 'net growth rate’ (S=b - d), and S_scaled is a beta distribution used to scale L. In other words, because strong correlations in the posterior distributions of b and d given the likelihood (since they both influence the net growth rate that is observed), the inventors postulated that the particular sampler used in STAN maystruggle to effectively sample these posteriors. They therefore reparameterised the problem in terms of S and L. Because it is known that (b + d) >= (b - d), L can just be scaled L a beta distribution (S_scaled) which is naturally bounded in the [0,1] interval to obtain S.

[0492] The Bayesian framework was able to closely recover the true parameters over a range of biologically feasible values from data generated using the stochastic agent-based model, validating the analytical approach to recovering average birth and death rates from a population using lineage distribution data (Fig.

[0493] 22).

[0494] Estimating birth and death rates with untreated samples - Bayesian inference of colorectal cancer cell line birth and death rates. The framework was implemented in the two barcoded CRC cell lines (HCT116 and SW620). For each cell-line, 3x ‘POT’ samples were harvested following the shared, mutual expansion stage. To infer the birth and death rates of each population, the Bayesian model outlined above was fit to all combined replicates. For barcoded MSI CRC cell line HCTbc, b = 0.693 days-1and d = 0.070 days-1. For barcoded MSS CRC cell line SW6bc, b = 0.656 days-1and d = 0.026 days-1(Fig. 23, top). These estimates were obtained using the simple birth / death model described above (equation (1), assuming a single population of S cells). This step comes before fitting the models A and B to the drug-treated experimental data and is used to estimate the birth / death rates of the treatment naive cells. Posterior predictive distributions of lineage counts given 100 draws from each cell line's birth and death rate posterior distributions are shown in Fig. 23 (bottom).

[0495] MATERIALS & METHODS - EXAMPLES 1-7: In vitro long-term experimental resistance evolution, phenotypic and molecular characterisation of phenotypes, bioinformatics analysis

[0496] Cell culture. Two colorectal cancer cell-lines HCT116 (ATCC CCL-247™) and SW620 (ATCC CCL-227™) were used for all data generating experiments in these examples. Cells were grown in 'standard conditions' which were as follows: growth medium, consisting of high glucose DMEM supplemented with 10% Foetal Bovine Serum (FBS) and 2% Penicillin-Streptomycin (herein 'full growth medium', unless stated otherwise). Cells were grown in T-175 vented flasks for the long-term drug-treatment experiment with 35mL of full growth medium. Flasks were grown at 37°C, 5% CO2 and 95% relative humidity. Cells were cryopreserved using a slow freezing technique. Briefly, cells were resuspended in freezing medium consisting of 90% FBS and 10% dimethyl sulfoxide (DMSO). The cell suspension was aliquoted into cryovials, which were then placed in freezing containers before being stored at -80°C overnight. After 24 hours, the vials were transferred to liquid nitrogen for long-term storage. Cell line identities were confirmed using STR analysis at the beginning of the experiment and were tested regularly for mycoplasma infection.

[0497] Lineage tracing. The ClonTracer library (Addgene #67267) was expanded in-house according to the original protocol (Bhang et al., 2015). The optimisation steps were also undertaken according to the original ClonTracer protocol such that the multiplicity of infection (m.o.i) of cells was estimated to be 0.1, minimising the number of multiple integrations. In addition, optimisation experiments were performed such that approximately 1x106cells (per cell line) could be confidently barcoded following infection and puromycin selection. The infection and selection steps were split over 3x150mm plates to ensure there was sufficientroom for population expansion following barcoding without the need for additional passaging events. Given the modelling framework used the lineage distributions enabled by these barcodes to estimate the evolutionary dynamics, single cell expansion and sequencing of representative clones was performed to make sure that the majority of cells did indeed contain a single barcode.

[0498] Long-term resistance evolution assay. Following barcoding of colorectal cancer cell lines with the ClonTracer genetic lineage tracing system, the initial population was expanded to ensure that each lineage was represented multiple times. Vials of this expanded population (POT) were frozen and 12 replicate subpopulations of 1x106cells each were seeded for 4 replicates across 3 conditions: vehicle control, CO, drugstop, DS, and drug-treatment, DT. This design was repeated in parallel for each barcoded cell line (HCTbc and SW6bc). Control replicates were treated with periodic vehicle control (DMSO) for 4 passages, drugtreatment replicates were treated with IC50 values of the chemotherapeutic 5-fluoruracil (5-Fu) for 4 passages and drug-stop were treated with a single period of IC505-Fu treatment then subjected to a single recovery period before harvesting and freezing (see Drug Inhibition Curves). All flasks were passaged before 80% confluency was reached. When passaging, 1x106cells were seeded into the subsequent replicate flask and the remaining sample’s cells were split into two vials. Frozen cells were kept for DNA extraction, barcode amplification and additional functional assays.

[0499] DNA extraction, barcode amplification and barcode sequencing. To extract genomic DNA (gDNA) for barcode amplification, cells were defrosted and a volume containing approximately 1x106cells was spun down into a pellet and re-suspended in 200pL of PBS. DNA was then extracted according to the DNeasy Blood & Tissue Kit (QIAGEN) protocol. DNA was eluted in DNase and RNase-Free PCR grade water. To quantify the concentration of DNA, the eluted volume was vortexed and then 1pL taken and measured using a Qubit 4 Fluorometer (Invitrogen) according to the manufacturer's protocol. DNA concentrations were recorded and the DNA frozen at -20°C until future processing. Following PCR steps, a magnetic bead based clean-up system was adopted for the purification of amplified barcode DNA (CleanNGS - CleanNA).

[0500] To amplify the barcodes from cell DNA, two universal amplification primers targeted the universal sequences flanking the semi-random barcode sequence in the ClonTracer construct directly. A bespoke adapter ligation protocol was optimised that used unique dual-indexes (IDT Illumina™ TruSeq UD Indexes 96) - this meant that each sample was identified via its own unique forward and reverse indexes. This approach minimised index hopping, a common problem for low complexity libraries sequenced on newer Illumina patterned flow cell systems (Costello et al., 2018). Final libraries were quantified using Agilent HS D1000 screentapes on an Agilent Tapestation 4200 system and sequenced on the Illumina Novaseq 6000 system using Novaseq S2 PE50 flow cells. Due to the low-diversity amplicon library, 15% of Phix control library was used.

[0501] Drug inhibition curves. The IC50 values of each cell line were measured as follows. Cells were trypsinised into a single-cell solution and cell viability was confirmed to be >80% using the Countess™ II Automated Cell Counter (Thermo Fisher Scientific). Cells were then seeded into a TC-treated 96-well adherent plate (8000 cells / well for HCT116 and 10000 cells / well for SW620) in 200pL of full growth medium. After 1 day (+24hrs from start), full growth medium + given drug concentrations were added to the respective wells.After three additional days (+96hrs from start), CellTiter-Glo® Reagent (Promega Ltd) was thawed and 50pL of reagent was added to all wells. Plates were placed in the incubator at standard conditions for 15 minutes before luminescent readings were taken using a PHERAstar FSX plate reader (BMG Labtech, Germany). An R script was used to calculate the dose-response curve from the results (Ritz et al., 2015). Single cell sorting for growth rate assays, scWGS and size-sorted scRNA. Barcoded HCT116 (HCTbc) and SW620 (SW6bc) cells from the long-term resistance evolution assay were single cell sorted for colony growth assays (HCTbc) and single-cell whole genome sequencing (scWGS, SW6bc). Vials were thawed and expanded briefly (<3 days) from the respective replicates. A cell count of >1x106and viability of >80% ensured there was no stringent bottleneck prior to characterisation. Single cell solutions were obtained by resuspending trypsinised cells in PBS and single cells were sorted into single wells using the cellenONE® single cell sorting platform (Cellenion). For the colony growth assays, cells were sorted into wells of a 96-well plate which contained conditioned growth media (50% standard media, 50% media from cells that were 80% confluent, filtered in 0.45pm size filter). For scWGS, cells were sorted into wells that contained a predistributed lysis buffer, as per the protocol in Laks et al. (2019. For scRNA-seq using the cellenCHIP 384-3’RNA Seq Kit, single cells were sorted directly into the chip: large cells were chosen as those with a diameter >36pm and an elongation factor <1.50 whilst small cells were chosen with a diameter filter of <20pm and an elongation factor <1.50.

[0502] Growth rate assays (HCTbc). To measure the growth rate of single cell colonies derived from 2 control (CO) and 2 drug -treatment (DT) HCTbc replicates, single cell sorted plates were kept in an Incucyte® S3 Live-Cell Analysis System, and confluence readings were recorded every 12 hours for 38 days. As most cells were sorted near the edge of the well, any wells that contained a cell sorted at the centre were excluded to control for well-position dependent differences in growth kinetics. After the full expansion duration, confluence readings were used to extract growth rates by fitting a linear regression model to the log-confluence values.

[0503] scl / l / GS (SW6bc). For whole genome sequencing of single cells from barcoded SW6bc cells, the protocol in Laks et al. (2019) was used. First, 1pl each of 384 unique i5 indexing primer at a concentration of 4pM were dispensed into empty 384-well plates and air-dried. 1 pl lysis buffer, consisting of 86.2% DirectPCR Lysis Reagent (Cell) (Viagen), 8.6% protease (Qiagen), 5.2% glycerol (Sigma) was added to each well and a single SW620 cell was sorted into each well using the CellenOne F1.4 system (Cellenion). Following centrifugation of the plate at 3000G for 5 min and overnight incubation at +4°C cell lysis was performed by incubating the plate at 50°C for 1h and 70°C for 15min. Tagmentation of genomic DNA was carried out by adding 1.8ul tagmentation mix, consisting of 1.775pl tagmentation buffer (20mM Tris-HCI pH 8, 10mM MgCI2, 20% dimethylformamide), 0.0165pl Tween, 0.00875pl Tn5 loaded (Diagenode), and incubating at 55°Cfor 10min. Reactions were neutralised by adding 1 pl neutralisation mix (0.5pl Qiagen protease, 0.01 pl Tween and 0.049pl ddH2O) and incubating at 50°C for 15min and 70°C for 10min. PCR amplification was performed after adding 3.8pl master mix (3.78pl 2x NEB Ultra II Q5 master mix (NEB), 1mM i7 indexing primer (500nM final, IDT)) with the following cycling parameters: Gap-filling at 72°C for 5min, initial denaturation at 98°C for 30s followed by 10 cycles of 98°C for 30s and 65°C for 75s and final elongation at 65°C for 5 min. PCR products were pooled and purified using a Zymo DCC5 spin column and then subjected to Exonuclease I digestion (NEB) and 1x Ampure bead purification. Final libraries were quantifiedusing Agilent HS D1000 screentapes on an Agilent Tapestation 4200 system and sequenced on the Illumina Novaseq 6000 system using Novaseq S2 PE50 flow cells.

[0504] scRNA sequencing (HCTbc and SW6bc). To perform scRNA-seq, cell samples were thawed (1x106-cells / sample vial) and cell viability was confirmed to be >80% using the Countess™ II Automated Cell (Thermo Fisher Scientific). Samples were treated for one additional round of drug-treatment (IC505-Fu for DS / DT conditions), vehicle control (CO condition) or grown in standard culture conditions (full growth medium - POT condition) for 3 days, followed by 3 days of standard culture conditions for all conditions. Cells were trypsinised and cell viability confirmed to be >80% and a fraction were re-suspended in PBS for an estimated number of 6000 captured cells per sample. The single cell suspension was loaded on a Chromium Single Cell 3' Chip C, followed by Chip D (10X Genomics). A single-cell gel bead-in-emulsion was generated using the Chromium Single Cell DNA kit and the Chromium Controller. The scRNA-seq library was prepared, and the concentration and quantity of the complete library was confirmed using the TapeStation High Sensitivity Screen Tape assay D1000 (Agilent) and a High-Sensitive Qubit™ dsDNA Kit (Life Technologies). Samples were normalised, pooled, and sequenced on an Illumina NovaSeq 6000 according to standard 10X Genomics recommendations at a median depth of at least 50K read pairs per cell. All samples from the same cell line were prepared on the same day and all samples from each cellline were sequenced on the same flow cell to minimise technical biases.

[0505] For the ‘size-sorted’ samples, experimental conditions (POT, CO, DS, DT) were re-sequenced from the original experiment (barcoded cell line HCTbc) alongside ‘large’ and 'small' cells following 7 weeks of treatment with IC50 values of 5- fluoruracil, and scRNA-seq was performed using the cellenCHIP 384-3’RNA-Seq Kit (Cellenion). This kit includes the cellenCHIP 384, using a nanowell array containing oligo-dT primers, unique cell barcodes (CB), and unique molecular identifiers (UMIs) for cDNA generation. The barcode oligos in the wells were rehydrated with Lysis and RT buffers followed by cell isolation and dispensing using the cellenONE (Cellenion). Reverse transcription was performed at 42°C for 90 minutes. The resulting cDNA was pooled by inverting the cellenCHIP 384 and centrifuging it into a recovery funnel for collection into microcentrifuge tubes. cDNA quantification was conducted using the High-Sensitive Qubit™ dsDNA Kit (Life Technologies), and amplification was carried out for up to 18 PCR cycles. The amplified cDNA was used to generate Illumina sequencing libraries, which were confirmed using the TapeStation High Sensitivity Screen Tape assay D1000 (Agilent).

[0506] Read merging, barcode extraction and clustering. Sequenced FASTQs from demultiplexed amplified barcode samples were merged using NGmerge (Gaspar, 2018) in ‘stitch mode’ with the following options: ‘-m 14, -p 0.2, -z’. The barcodes were then extracted and clustered using ‘Bartender’ (Zhao et al., 2018) with the following options: (extract) ‘-q? -p GACAG

[0030] AGCAG -m 2 -d b’; (cluster) ‘-c 10 -d 2 -z 5’. Due to issues with index-hopping encountered when previously sequencing lineage tracing marks without the use of non-redundant indexes, unique dual indexes (UDIs) were used for barcode sequencing. Despite this, patterns suggestive of index-hopping were observed in some samples (albeit to a much lesser degree). For example, barcodes that were found at high frequency in one cell line would often be found in low (but proportional) frequencies in the alternate cell line, despite independent barcoding steps and therefore putatively independent lineage relationships. These patterns also only appeared when samples were sequenced on the same flow-cell, confirming their technical nature. A list of ‘problem barcodes’ was createdbased on these patterns, which were excluded when calculating lineage diversity statistics. Of note, as diversity statistics that leveraged higher frequency lineages were chosen, these exclusions of low frequency ‘swaps’ had minimal impact when calculating the diversity statistics.

[0507] scRNA-seq expression and mutation analysis. For each barcoded cell line processed (HCTbc and SW6bc), the raw sequencing data were processed using the Cell Ranger software (10X Genomics) to perform sample de-multiplexing, barcode processing, and single-cell 3' gene counting. The cellranger count pipeline was used for each sample, aligning reads to the GRCh38-2020-A reference genome and generating feature- barcode matrices for further analysis. The resulting filtered feature-barcode matrices were imported into R for downstream analysis using the Seurat package (Hao et al., 2021). Cells expressing fewer than 500 genes and with fewer than 1000 RNA molecules detected, and cells with more than 20% mitochondrial genes expressed were excluded from further analysis to remove low-quality cells. To standardize cell numbers across samples and prevent biases due to variations in the number of single cells sequenced, 6000 cells per sample were used for scRNA-seq analysis (the sample with the lowest cell count had >6000 cells sequenced). The ‘DoubletFinder’ package was used to find and exclude putative doublets (McGinnis, Murrow and Gartner, 2019). Expression was normalised with the SCTransform function (Seurat v4.10) and the normalised expression matrix was used for all downstream analysis. The K-nearest neighbour graph based on the PCA-reduced dimensions was identified using the FindNeighbors function (Seurat v4.10). The first 30 components identified by PCAwere used for analyses, chosen based on the elbow plot method to capture the majority of the variance in the data. Clustering was performed using the FindClusters function in Seurat using resolution=0.2. The clusters were identified based on the Louvain algorithm applied to the K-nearest neighbour graph.

[0508] To call mutations from the scRNA-seq data, the output Cell Range BAMs with SComatic (Muyas et al., 2023) were used with the following settings - Splitting Alignment Files: n_trim = 5; max_nM = 5; max_NH = 1; Collecting Base Count Information: min_dp = 10; min_cc = 10; min_bq = 30; Detection of Somatic Mutations: min_cov = 10; min_ac_cells = 3; min_ac_reads = 4; max_cell_types = 7; min_cell_type = 4.

[0509] scl / l / GS analysis. Demultiplexed dual-index fastq files were obtained for single cells undergoing the DLP+protocol and the data were used as input for the workflow automation pipeline designed for the DLP+ method (https: / / github.com / shahcompbio / single_cell_pipeline) using default settings. Raw reads were adapter-trimmed using TrimGalore, mapped with BWAaln to the hg19 reference genome and deduplicated using picard MarkDuplicates. Subsequently, copy number calling was performed using the HMMcopy tool with reads segmented into non-overlapping 500kb genomic regions (Daniel Lai, 2017). Cells were required to have a minimum quality score of 95, a minimum of 250,000 reads and an S-phase probability of less than 0.1. Only bins considered “ideal” in all cells were used.

[0510] Differential Expression Analysis. Differential expression analysis was performed in R using the edgeR package (Robinson, McCarthy and Smyth, 2010). scRNA-seq data (raw RNA counts) was combined into pseudobulk expression count matrices for each replicate sample. Each barcoded cell line (HCTbc and SW6bc) was analysed separately. Negative-binomial dispersion was estimated using the ‘estimateDisp’ function, and then fit the GLM to the count data with the ‘glmQLFit’ and ‘glpQLFTest’ to derive log-fold changes and FDR values for each gene. When identifying differentially expressed genes that wereassociated with a specific sample comparison, those with a logFC > 2.0 and FDR < 0.05 were selected. When choosing genes that were not differentially expressed given the comparison between two samples, genes with logFC < 1.0 and FDR > 0.10 were selected. These subsets of genes were used to perform gene ontology (GO) or KEGG pathway analysis (Ritchie et al., 2015). For gene set enrichment analysis (GSEA) (Korotkevich et al., 2016), custom ranks were created based on the combinations of sample comparisons in question: if looking for genes associated with each comparison, the product of both logFC values was taken. If looking for genes whose expression was associated with one sample comparison but not the other, a score was generated with logFC_1 x (1 / (|logFC_2|+1 )), where logFC_1 is the log fold change in expression in comparison T, the comparison where of interest in differentially expressed genes. This custom score simultaneously penalises genes with a high differential expression in comparison ‘2’ (logFC_2).

[0511] EXAMPLE 1 - QUANTITATIVE MODELS OF RESISTANCE EVOLUTION

[0512] The inventors developed mathematical models of resistance evolution that described the response of a population of cancer cells to periodic drug treatment. They aimed to infer the dynamics solely using the change in total population size and relatedness between surviving cell lineages. The evolution of phenotypes within a cell population that control the response to treatment and the transitions between these phenotypes were modelled. The inventors began with simple models of resistance evolution and added additional complexity only when simple models failed to explain data. Three models of increasing complexity were designed (Models A, B and C), which encompassed behaviours previously observed during the emergence of cancer cell resistance.

[0513] Model A. The first and simplest model (herein 'unidirectional transitions’) consists of two phenotypes: sensitive and resistant. A 'pre-existing resistance fraction’ parameter (p) controls the initial conditions of the model by setting the proportion of cells with the resistant phenotype at the start of the simulation (to). Cells can divide with phenotype-specific birth and death rates (sensitive: bs,ds and resistant: bR,d ). A ’cost’ parameter (6) controls a fitness penalty associated with the resistant phenotype in the untreated environment. Including a fitness penalty allowed testing of the assumption that, in the absence of treatment, resistant cells exhibit a reduced net growth rate relative to sensitive cells. This behaviour has been observed in experimentally generated resistant cancer cell lines (Jensen et al., 2015), and is significant for phenomena including the existence of slowly dividing ‘persister’ cells (Oren et al., 2021) and drug ‘addiction’, whereby resistant cells fail to grow in the absence of treatment (Maltas et al., 2024). Fitness penalties associated with a resistant phenotype also influences tumour containment strategies (Viossatand Noble, 2021). A switching parameter (p) was included which controls the probability of cells transitioning from the sensitive to resistant phenotype per cell division.

[0514] Therapy was modelled by modifying the death rate of cells. For sensitive cells, the death rate was increased as a function of D(t), with D(t) representing the effective drug concentration at time t. Omitting the pharmacokinetics of the drug treatment may impair the model's ability to recover observed changes in cell population sizes, as the cytotoxic effects of a drug may not be realised immediately. A pharmacokinetic model was assumed where the focus was on the change in the realised effect of the drug on cells, rather than directly modelling the drug concentration. The parameters De and K were included and regulate thestrength and rate of accumulation / reduction of this effect, respectively. This model does not distinguish between intrinsic and extrinsic cellular factors; for instance, the increasing cytotoxicity experienced by cells following the addition of treatment could be affected by the drug’s chemistry, cellular efflux mechanisms, or some combination of the two. Furthermore, no assumption was made concerning the specific molecular change or mechanism that controls resistance in the models. Instead, the focus was on the phenotypes dictating the probability of survival during treatment. The parameter qj was used to modelle the survival probability of the resistant cells under treatment: qj = 1.0 denotes complete resistance, whereas when qj = 0.0, the sensitive and resistant cells experience the same level of drug-induced death.

[0515] To summarise, in Model A, the total cell population’s response to treatment is a product of the proportion of resistance (dictated by pre-existing resistance fraction parameter p, cost parameter 6 and switching parameter p), the effective strength of the drug at a given time (dictated by De and K, the parameters regulating the strength and rate of accumulation / reduction of the realised effect of the drug on the cells) and the strength of the resistant phenotype (dictated by qi, the survival probability of the resistant cells under treatment).

[0516] Model A allowed simpler evolutionary scenarios of interest to be explored: compared to sensitive cells, resistant cells can either be absent or pre-exist at varying frequencies (by...

Claims

1. Claims:

1. A method of characterising an evolution of resistance in a population of proliferative cells exposed to an anti-proliferative treatment, the method comprising:obtaining, by a processor, experimental data comprising total population size and lineage tracing data measured at one or more time points after the start of exposure of the cells to the antiproliferative treatment;fitting, by the processor, one or more first dynamic cell population growth models to the total population size data, thereby identifying a first set of values of each of one or more parameters of each of the one or more dynamic cell population growth models that meet a first predetermined model fit criterion, wherein the dynamic cell population growth models represent the growth of a plurality of subpopulations of cells having different responses to the anti-proliferative treatment comprising at least a first subpopulation of cells that is sensitive to the anti-proliferative treatment and a second population of cells that is resistant to the anti-proliferative treatment; andfitting, by the processor, one or more second dynamic cell population growth models to the lineage tracing data, each second dynamic cell population growth model representing the growth of the same plurality of subpopulations of cells and having the same parameters as a corresponding first dynamic cell population growth model, said fitting comprising evaluating the fit to data of the one or more second dynamic cell population growth models to the lineage tracing data parameterised using parameter values derived from the first set of values, thereby identifying a second set of values of each of one or more parameters of each of the one or more first and second dynamic cell population growth models that meet a second predetermined model fit criterion;wherein the second set of values of one or more parameters comprise: values that are indicative of the effect of the treatment on the growth of the plurality of subpopulations of cells, values that are indicative of the proportions of the one or more subpopulations of cells prior to exposure to the anti-proliferative treatment, and / or values that are indicative of the rate of transition of cells between the plurality of subpopulations of cells comprising a rate of transition between sensitive and resistant subpopulations of cells.

2. The method of claim 1, wherein for each second dynamic cell population growth model there is a corresponding first dynamic cell population growth model that represents the growth of the same plurality of subpopulations of cells using the same parameters, wherein the second dynamic cell population growth models represents the temporal dynamics of the numbers of cells of each of a plurality of cell lineages that belong to each of the plurality of subpopulations of cells (ms, niR, nT), and wherein the first dynamic cell population growth model represents the temporal dynamics of the total number of cells in each of the plurality of subpopulations of cells (ns, HR, HE).

3. The method of claim 1 or claim 2, wherein the second set of values of one or more parameters comprise values that are indicative of one or more of: a proportion of cells that are resistant before the population of cells is exposed to the anti-proliferative treatment (p), a rate of transition of sensitive cells to resistant cells (p), a probability of treatment-induced death for resistant cells (qj), a fitness penalty of resistant cells in the absence of treatment (6), a rate at which resistantcells escape a fitness penalty in the absence of treatment (a), a rate of change of a treatment induced effect on a death rate of resistant and / or sensitive cells upon exposure to the treatment(K), and a maximum strength of a treatment induced effect on a death rate of resistant and / or sensitive cells (De), optionally wherein the value of a rate of transition of sensitive cells to resistant cells (p) is indicative of whether transition between sensitive and resistant cells is likely to be due to one or more genetic mechanisms or one or more non-genetic mechanisms.

4. The method of any preceding claim, further comprising:obtaining, by a processor, experimental data comprising total population size and lineage tracing data measured at one or more time points after the start of exposure of the cells to a second anti-proliferative treatment;repeating the step of fitting the one or more first dynamic cell population growth models using the total population size data for the second anti-proliferative treatment and the step of fitting the one or more second dynamic cell population growth models using the lineage tracing data for the second anti-proliferative treatment; andfitting one or more third dynamic cell population growth models jointly to lineage tracing data comprising the lineage tracing data for the two anti-proliferative treatments, wherein each third dynamic cell population growth model represents the growth of a plurality of subpopulations of cells comprising at least a first subpopulation of cells that is sensitive to both anti-proliferative treatments, a second population of cells that is resistant to the first anti-proliferative treatment, a third population of cells that is resistant to the second anti-proliferative treatment, and a fourth population of cells that is resistant to both of the first and second anti-proliferative treatments; optionally wherein fitting the cross-resistance model comprise estimating, using both lineage tracing data and the parameters estimated by fitting the second dynamic cell population models independently for the two different anti-proliferative treatments, the value of a crossresistance parameter (0), wherein the cross-resistance parameter (0) is a parameter that quantifies the strength of statistical association of inheritance of resistance phenotypes to the two anti-proliferative treatments; and / or said fitting comprising evaluating the fit to data of the one or more third dynamic cell population growth models to the lineage tracing data parameterised using parameter values derived from the second sets of values for each of the two anti-proliferative treatments, thereby identifying a third set of values of the cross-resistance parameter that meet a third predetermined model fit criterion.

5. The method of any preceding claim, wherein: (i) fitting, by the processor, the one or more first dynamic cell population growth models to the total population size data and / or fitting, by the processor, the one or more second dynamic cell population growth models to the lineage tracing data and / or fitting, by the processor, the one or more third dynamic cell population growth models to the lineage tracing data, comprise using an approximate Bayesian computation inference approach, and / or using a simulation-based likelihood free method, and / or using a likelihood-free method that relies on comparing simulated data to observed data using a distance measure and generating a posterior distribution of parameters that produce simulations that are within a predetermined distance of the empirical data and so satisfy the first and / or second and / or thirdpredetermined model fit criteria; and / or (ii) the first predetermined model fit criterion applies to the value of a distance metric between the measured total population size data and simulated total population size data corresponding to a first dynamic cell population growth model and a set of values for the one or more parameters of the first dynamic cell population growth model; and / or wherein the second predetermined model fit criterion applies to the value of a distance metric between the measured lineage tracing data and simulated lineage tracing data corresponding to a second dynamic cell population growth model and a set of values for the one or more parameters of the second dynamic cell population growth model; and / or (iii) fitting, by the processor, the one or more first dynamic cell population growth models to the total population size data and / or fitting, by the processor, the one or more second dynamic cell population growth models to the lineage tracing data and / or fitting, by the processor, the one or more third dynamic cell population growth models to the lineage tracing data, comprise using a likelihood free method that relies on training a conditional density estimator to estimate a probability density function for the parameters of the model given simulated data, and using the trained model to predict a posterior distribution of parameters based on corresponding observed data instead of the simulated data, the posterior distribution thereby satisfying the respective predetermined model fit criterion.

6. The method of any preceding claim, wherein:(a) fitting, by the processor, the one or more first dynamic cell population growth models to the total population size data comprises:simulating the one or more first dynamic cell population growth models using a first plurality of sets of parameter values for each of the one or more parameters of each first dynamic cell population growth model, the first plurality of sets of parameter values comprising the first set of values, thereby obtaining simulated total cell population size data for each of the first dynamic cell population growth models and each of the first plurality of sets of values for each of the one or more parameters;determining whether each set of values of parameters of each first dynamic cell population growth model meets the first predetermined model fit criterion by comparing the corresponding simulated total cell population size data to the measured total cell population size data; and selecting the first set of values as a subset of the plurality of sets of values that meet the first predetermined model fit criterion; and / or(b) fitting, by the processor, the one or more second dynamic cell population growth models to the lineage tracing data comprises:simulating the one or more second dynamic cell population growth models using a second plurality of sets of parameter values for each of the one or more parameters of each second dynamic cell population growth model, the second plurality of sets of parameter values derived from the first set of values, thereby obtaining simulated lineage tracing data for each of the second dynamic cell population growth models and each of the second plurality of sets of values;determining whether each set of values of parameters of each second dynamic cell population growth model meets the second predetermined model fit criterion by comparing the corresponding simulated lineage tracing data to the measured lineage tracing data; andselecting the second set of values as a subset of the second plurality of sets of values that meet the second predetermined model fit criterion; orsimulating the one or more second dynamic cell population growth models using a second plurality of sets of parameter values for each of the one or more parameters of each second dynamic cell population growth model, the second plurality of sets of parameter values derived from the first set of values thereby obtaining simulated lineage tracing data for each of the second dynamic cell population growth models and each of the second plurality of sets of values;training one or more machine learning models to predict conditional density distributions of the parameters given the simulated lineage tracing data, andusing the trained one or more machine learning models to predict a posterior distribution of the parameters given the observed lineage tracing data any set of values drawn from these posterior distributions meeting the second predetermined model fit criterion; and / or(c) fitting, by the processor, the one or more third dynamic cell population growth models to the lineage tracing data comprises:simulating the one or more third dynamic cell population growth models using a third plurality of sets of parameter values for each of the one or more parameters of each third dynamic cell population growth model, the third plurality of sets of parameter values derived from the second set of values for each of the first and second anti-proliferative treatments, thereby obtaining simulated lineage tracing data for each of the third dynamic cell population growth models and each of the third plurality of sets of values;determining whether each set of values of parameters of each third dynamic cell population growth model meets the third predetermined model fit criterion by comparing the corresponding simulated lineage tracing data to the measured lineage tracing data; andselecting the third set of values as a subset of the third plurality of sets of values that meet the third predetermined model fit criterion; orsimulating the one or more third dynamic cell population growth models using a third plurality of sets of parameter values for each of the one or more parameters of each second dynamic cell population growth model, the third plurality of sets of parameter values derived from the second sets of values for each of the first and second anti-proliferative treatments, thereby obtaining simulated lineage tracing data for each of the third dynamic cell population growth models and each of the third plurality of sets of values;training one or more machine learning models to predict conditional density distributions of the parameters given the simulated lineage tracing data, andusing the trained one or more machine learning models to predict a posterior distribution of the parameters given the observed lineage tracing data any set of values drawn from these posterior distributions meeting the third predetermined model fit criterion.

7. The method of any preceding claim, wherein the one of more first dynamic cell population growth models and the one or more second dynamic cell population growth models are birth-death models that represent the number of cells alive in each of a plurality of subpopulations of cells as a function of time, the plurality of subpopulations of cells having different responses to the antiproliferative treatment comprising at least a first subpopulation of cells that is sensitive to the anti-proliferative treatment and a second population of cells that is resistant to the anti-proliferative treatment, wherein the different subpopulations of cells differ by the birth rate in the presence of the treatment, their death rate in the presence of the treatment, their birth rate in the absence of treatment, and / or their death rate in the absence of treatment.

8. The method of any preceding claim, wherein the one of more first dynamic cell population growth models are hybrid differential equation and stochastic jump process models, and / or wherein the one of more second dynamic cell population growth models are stochastic agent based models and / or wherein the one of more third dynamic cell population growth models are stochastic agent based models.

9. The method of claim 8, wherein a hybrid differential equation and stochastic jump process model is a model that represents the growth of each subpopulation of cells using a stochastic process when the total number of cells in the subpopulation is below a predetermined threshold, and using a deterministic differential equation model when the total number of cells in the subpopulation is at or above the predetermined threshold.

10. The method of any preceding claim, wherein the first and second dynamic cell population growth models, and optionally the third dynamic cell population growth models when used, are birth-date models that capture births and deaths of cells in each of the plurality of subpopulations of cells, wherein transitions between specific subpopulations of cells are associated with respective rates and coupled to birth events,optionally wherein the rate of transition between a first subpopulation and a second subpopulation is a parameter between 0 and 1 that represents the probability of a cell in the first subpopulation generating a cell in the second subpopulation upon cell division, and / or wherein the value of the rate of transition between a first subpopulation and a second subpopulation is indicative of the likely mechanism of action of the transition, wherein lower transition rate values are more likely to be associated with genetic or epigenetic mechanisms than higher transition rates, and higher transition rates are more likely to be associated with non-genetic mechanisms than lower transition rates;optionally wherein the first and second dynamic cell population growth models each comprise a rate of transition of sensitive cells to resistant cells (pi), and optionally wherein at least one first dynamic cell population growth model and at least one corresponding second dynamic cell population growth model further comprise a rate of reversion of transition of resistant cells to sensitive cells (a).

11. The method of any preceding claim, wherein at least one first dynamic cell population growth model and corresponding second dynamic cell population growth model represent the growth of a subpopulation of cells that are resistant to the treatment and do not incur a fitness penalty in the absence of treatment; and / or wherein at least one first dynamic cell population growth model and corresponding second dynamic cell population growth model represent the growth of a subpopulation of cells that are resistant to the treatment and do not incur a fitness penalty in theabsence of treatment, and a subpopulation of cells that are resistant to the treatment and do incur a fitness penalty in the absence of treatment, and wherein said first and second dynamic cell population growth models comprise a rate of transition of resistant cells that do incur a fitness penalty in the absence of treatment escape to resistant cells that do not incur a fitness penalty in the absence of treatment, optionally wherein the rate of transition is obtained as the product of a base rate (a) and a time-dependent factor (y(t)) that increases at a first predetermined rate (K) upon exposure to the treatment and decreases at a second predetermined rate (-K) upon withdrawing of the treatment, optionally wherein the first and second predetermined rates have the same absolute value.

12. The method of any preceding claim, wherein the lineage tracing data comprises counts of cells that have the same clonal identity or data derived therefrom, optionally wherein said counts are based on the number of cells that share each of a plurality of heritable barcodes, and / or wherein the lineage tracing data comprises diversity statistics derived from counts of cells that have the same clonal identity,optionally wherein the experimental data comprises total population size and lineage tracing data measured at one or more time points after the start of exposure of the cells in each of a plurality of replicates of the cell population exposed to the anti-proliferative treatment, and wherein the diversity statistics comprise a first statistic that is representative of lineage diversity within the replicates and a second statistic that is representative of lineage diversity between the replicates; and / or wherein the diversity statistics are Hill Diversity statistics; and / or wherein the lineage tracing data comprises first lineage tracing data for a first anti-proliferative treatment and second lineage tracing data for a second anti-proliferative treatment, and wherein the diversity statistics comprise diversity statics that quantify the degree to which shared lineages are enriched between the lineage tracing data associated with exposure to the first and second anti-proliferative treatments.

13. The method of any preceding claim, wherein the experimental data comprises total population size and lineage tracing data measured for a cell population in vitro; optionally wherein the experimental data comprises total population size and lineage tracing data measured for a plurality of replicates of a cell population derived from the same parent population and exposed to the anti-proliferative treatment in vitro; and / orwherein the total population size data comprises total cell counts acquired using a non-disruptive technology, optionally imaging, and / or wherein the lineage tracing data comprises counts of cells in respective clonal lineages acquired using a sequencing technology, optionally wherein the lineage tracing data is data acquired at passaging steps applied to the in vitro cell culture.

14. The method of any preceding claim, wherein the method comprises fitting a plurality of first dynamic cell population growth models and a corresponding plurality of second dynamic cell population growth models, wherein the models represent different assumptions of how resistance occurs, and wherein the method further comprises comparing the plurality of fitted second dynamic cell population growth models using a third predetermined model fit criterion, andidentifying a second dynamic cell population growth model that best fits the experimental data based on said comparison.

15. The method claim 14, wherein the third predetermined model fit criterion is a model fit criterion that takes into account the numbers of fitted parameters of the models that are being compared, and / or wherein the third predetermined model fit criterion is a DIC score.

16. The method of any preceding claim, wherein growth in the first dynamic population growth models is modelled as density-dependent growth, optionally logistic growth with a fixed carrying capacity (K).

17. The methods of any preceding claim, wherein the anti-proliferative treatment comprises one or more anti-proliferative compounds or compositions, wherein the one or more anti-proliferative compounds or compositions comprise one or more active agents selected from: small molecules, large molecules, optionally antibodies, antigen binding molecules, peptides, nucleic acids, optionally miRNA, siRNAs and gene editing constructs, combinations or large and small molecules, optionally selected from antibody drug conjugates, bicycle peptide-drug conjugates, cells, optionally T cells or modified T cells, optionally CAR-T cells orTCR-T cells, and targeted protein degradation therapies, optionally selected from PROteolysis TArgeting Chimeras (PROTACs), molecular glues, Lysosome-Targeting Chimaeras (LYTACs), and Antibody-based PROTACs (AbTACs).

18. The methods of any preceding claim, wherein the first dynamic population growth models comprise:a) a model (A) that represents a first subpopulation of cells that is sensitive to the antiproliferative treatment and a second population of cells that is resistant to the antiproliferative treatment, wherein the sensitive cells can transition to resistance cells but resistant cells cannot transition to sensitive cells, optionally wherein the resistant cells incur a fitness penalty in the absence of the anti-proliferative treatment; and / orb) a model (B) that represents a first subpopulation of cells that is sensitive to the antiproliferative treatment and a second population of cells that is resistant to the antiproliferative treatment, wherein the resistant cells incur a fitness penalty in the absence of the anti-proliferative treatment, and wherein the sensitive cells can transition to resistance cells and resistant cells can transition to sensitive cells; and / orc) a model (C) that represents a first subpopulation of cells that is sensitive to the antiproliferative treatment, a second population of cells that is resistant to the anti-proliferative treatment and that incur a fitness penalty in the absence of the anti-proliferative treatment, and a third population of escaped cells that is resistant to the anti-proliferative treatment and that do not incur a fitness penalty in the absence of the anti-proliferative treatment, wherein the sensitive cells can transition to resistance cells and resistant cells can transition to sensitive cells, and wherein the resistant cells can transition escaped cells but escaped cells cannot transition to resistant cells.

19. The method of any preceding claim, wherein one or more of the first dynamic population growth models represent an effect of the anti-proliferative treatment on cells in each subpopulation using an effective death rate that depends on a base death rate and a treatment effect term, wherein the treatment effect term is a time depend term that represents an exposure dependent effect of the anti-proliferative treatment on the cells death rate, optionally wherein the treatment effect term comprises a time dependent term that represents the effective concentration of the antiproliferative treatment and a subpopulation specific factor that is a parameter that captures the subpopulation specific effect of the anti-proliferative treatment on the cell subpopulation.

20. The method of any preceding claim, wherein the method further comprises exposing one or more replicates of the cell population in an in vitro cell culture to the anti-proliferative treatment and measuring total population size and lineage tracing data at one or more time points after the start of exposure of the cells to the anti-proliferative treatment,optionally wherein said measuring comprises obtaining both total population size and lineage tracing data at one or more first time points after the start of exposure of the cells to the antiproliferative treatment, and optionally obtaining total population size data at one or more second time points after the start of exposure of the cells to the anti-proliferative treatment; and / or wherein the method further comprises exposing one or more replicates of the cell population in an in vitro cell culture to a first and second anti-proliferative treatments in respective replicates derived from the same parental population of barcoded cells.

21. The method of any preceding claim, wherein the obtained lineage tracing data and total population size data have been obtained from a cell population cultured in vitro, wherein the cell population was cultured using a process comprising:obtaining lineage tracing data for a parental population of cells, optionally wherein the parental population of cells has been obtained by barcoding a population of cells and expanding the barcoded population of cells during a predetermined period of time,separating the parental population of cells into a first plurality of replicates, and optionally a second plurality of replicates.exposing each of the first plurality of replicates to the anti-proliferative treatment for one or more predetermined periods of time, optionally separated by one or more predetermined periods of time in the absence of the anti-proliferative treatment, andpassaging the cells in each of the first plurality of replicates, and optionally each of the second plurality of replicates, one or more times, optionally wherein lineage tracing data is acquired at passaging time points.

22. The method of any preceding claim, further comprising identifying one or more likely mechanisms of resistance using the parameters of (optionally selected) fitted second dynamic cell population growth models, optionally further comprising selecting one or more resistance validation experiments based on the one or more likely mechanisms identified, and / or outputting a report comprising one or more results of the method or information derived therefrom;optionally wherein the one or more mechanisms comprise a genetic mechanism and the selected one or more resistance validation experiments comprise a genome sequencing step, optionally whole genome sequencing or single cell sequencing; and / orwherein the one or more mechanisms comprise a non-genetic mechanism and the selected one or more resistance validation experiments comprise a transcriptome and / or chromatin state sequencing step, optionally single cell RNA sequencing step and / or a single cell ATAC-seq step; and / orwherein the one or more mechanisms include the presence of a subpopulation of cells prior to exposure to the treatment that is resistant to the treatment and experiences a fitness penalty in the absence of the treatment, and the selected one or more resistance validation experiments comprises isolating said subpopulation of cells and performing one or more of a genome sequencing step, a transcriptome sequencing step, and / or a chromatin state sequencing step; and / orwherein the one or more results of information derived therefrom include one or more of: a value of one or more fitted parameters, a value of a fit criterion, an indication of an identified dynamic cell population growth model that best fits the experimental data, and an indication of a likely mechanism of resistance.

23. A method of screening a plurality of candidate anti-proliferative treatments, the method comprising:(i) characterising an evolution of resistance in a population of proliferative cells exposed to each anti-proliferative treatment using the methods of any preceding claim; and (ii) comparing the results of said characterising to identify anti-proliferative treatments of the candidate anti-proliferative treatments that are less likely to be associated with an evolution of resistance.

24. The method of claim 23, further comprising selecting one or more candidate anti-proliferative treatments for further characterisation, optionally comprising pre-clinical and / or clinical characterisation, using the results of said comparing.

25. A method of identifying a combination of anti-proliferative treatments that has a reduced risk of evolution of resistance compared to a subset of the anti-proliferative treatments in the combination, the method comprising:(i) characterising an evolution of resistance in a population of proliferative cells exposed to the subset of the anti-proliferative treatments using the methods of any of claims 1 to 22; and(ii) characterising an evolution of resistance in a population of cells exposed to one or more combinations of anti-proliferative treatments comprising respective additional antiproliferative treatments in addition to the subset of the anti-proliferative treatments using the methods of any preceding claim and comparing the results of said characterising to identify combinations of anti-proliferative treatments that are less likely to be associated with an evolution of resistance than said subset of anti-proliferative treatments; and / or(iii) identifying and performing one or more validation experiments based on the results of the characterising in i, and identifying one or more additional anti-proliferative treatments that target a resistance evolution mechanism identified using said validation experiments.

26. A method of identifying one or more biomarkers of resistance of a population of proliferative cells to an anti-proliferative treatment, the method comprising:(i) characterising an evolution of resistance in one or more different populations of proliferative cells exposed to the anti-proliferative treatment using the methods of any of claims 1 to 22; and(ii) identifying one or more likely mechanisms of resistance using the parameters of (optionally selected) fitted second dynamic cell population growth models for each of the one or more different populations of proliferative cells and performing one or more resistance validation assays based on the identified likely mechanisms of resistance to identify one or more biomarkers of resistance; and / or(iii) comparing the results of said characterising between a plurality of the different populations of proliferative cells to identify one or more biomarkers associated with one or more of the plurality of populations of proliferative cells likely to be associated with an evolution of resistance.

27. One or more non-transitory computer readable media comprising instructions that, when executed by one or more processors, cause the one or more processors to perform the steps of any of claims 1 to 26.

28. A system comprising: a processor; and a computer readable medium comprising instructions that, when executed by the processor, cause the processor to perform the steps of any of claims 1 to 26.

Images