Understanding intrinsic and acquired resistance is crucial to overcoming cancer chemotherapy failure. While it is well-established that intratumor, subclonal genetic and phenotypic heterogeneity significantly contribute to resistance, it is not fully understood how tumor sub-clones interact with each other to withstand therapy pressure. Here, we report a previously unrecognized behavior in heterogeneous tumors: cooperative adaptation to therapy (CAT), in which cancer cells induce co-resistant phenotypes in neighboring cancer cells when exposed to cancer therapy. Using a CRISPR/Cas9 toolkit we engineered phenotypically diverse non-small cell lung cancer (NSCLC) cells by conferring mutations in Dicer1, a type III cytoplasmic endoribonuclease involved in small non-coding RNA genesis. We monitored three-dimensional growth dynamics of fluorescently-labeled mutant and/or wild-type cells individually or in co-culture using a substrate-free NanoCulture system under unstimulated or drug pressure conditions. By integrating mathematical modeling with flow cytometry, we characterized the growth patterns of mono- and co-cultures using a mathematical model of intra- and interspecies competition. Leveraging the flow cytometry data, we estimated the model’s parameters to reveal that the combination of WT and mutants in co-cultures allowed for beneficial growth in previously drug sensitive cells despite drug pressure via induction of cell state transitions described by a cooperative game theoretic change in the fitness values. Finally, we used an ex vivo human tumor model that predicts clinical response through drug sensitivity analyses and determined that cellular and morphologic heterogeneity correlates to prognostic failure of multiple clinically-approved and off-label drugs in individual NSCLC patient samples. Together, these findings present a new paradox in drug resistance implicating non-genetic cooperation among tumor cells to thwart drug pressure, suggesting that profiling for druggable targets (i. e. mutations) alone may be insufficient to assign effective therapy.
Lung cancer is the leading cause of cancer-related death in the United States, and non-small cell lung cancer (NSCLC) accounts for the majority of lung cancer cases each year [1]. Recent advances in the molecular understanding of NSCLC progression have informed the development of new targeted therapies that are safer and more effective than standard chemotherapy: of the nearly two-thirds of patients who have an oncogenic driver mutation, about half have of these are druggable [2]. ATP-competitive small molecule inhibitors of mutant epidermal growth factor receptor (e. g. , erlotinib, gefinitib, afatinib, osimertinib), of mutant serine/threonine kinase b-raf (e. g. , vemurafenib, dabrafenib), and of mutant anaplastic lymphoma kinase or of ROS1 proto-oncogene receptor tyrosine kinase (e. g. , crizotinib, ceritinib, alectinib) have been approved for management of NSCLC [3]. Despite the increasing number of agents that have entered the clinic, successful treatment of NSCLC is hampered by drug resistance [4]. Two primary forms of resistance are predominantly studied: intrinsic and acquired [5]. Intrinsic resistance is largely considered to be due to aberrations including somatic mutations and DNA amplifications, which render primary therapy failures [6]. Acquired resistance, however, can emerge under selective pressures, and during treatment [7]. More recently, phenotypic plasticity and the role of ‘adaptive’ and drug-induced cell state transitions have introduced a new paradigm in acquired drug resistance [8,9]. It is increasingly clear that an underlying driver of acquired resistance is due to the dramatic range of genetic and phenotypic diversity that is conferred during tumor evolution, comprising both passenger and driver mutations [10]. This heterogeneity within tumors, termed intratumor heterogeneity (ITH), can significantly impact responses to therapy and sensitivities to certain targeted agents, such as those listed above [11]. In NSCLC, intrinsically resistant tumor subpopulations can expand, or drug-tolerant cells can acquire alterations that confer more robust resistance to therapy [4,12–15]. To improve the effectiveness of cancer therapeutics, it is therefore critical to understand whether, and to what extent, phenotypically distinct cells within the same tumor interact to promote resistance. In combination with experimental approaches, quantitative methods arising from evolutionary theory have contributed to our understanding of the role of genetic heterogeneity in tumor initiation, progression, and resistance [16–24]. Mathematical models have previously demonstrated that local cell density and cell cycle can contribute to spatiotemporal heterogeneity and differences in cell response to treatment [25]. However, applying computational approaches to examine phenotypic heterogeneity has been less well pursued [26]. Here, we explore how intratumor phenotypic heterogeneity in a model of NSCLC affects the acquisition of resistance and population growth. To generate phenotypically unique sub-clones, we use the CRISPR/Cas9 system to stochastically mutate Dicer1, a cytoplasmic riboendonuclease responsible for maturation of microRNA [27]. Indeed, dysregulation of microRNA is heavily implicated in NSCLC [28]. Using a three-dimensional in vitro culture platform, which most closely parallels in vivo growth, we show that intrinsically-sensitive, phenotypically distinct NSCLC sub-clones support each other to rapidly evolve a therapy-resistant phenotype via drug-induced cell state transitions, a behavior we term “cooperative adaptation to therapy” (CAT). We develop an in vitro-validated model of intra- and interspecies competition governed by replicator dynamics [29] to simulate how CAT affects population growth over time. These findings build on previous evidence that drug-induced phenotypic transitions can underpin resistance, and implicate a population-wide impact of this phenomenon.
Anonymous non-small cell lung cancer (NSCLC) tissue samples were collected under IRB approval with due written consent from each patient. Unless noted otherwise, all reagents, small molecule inhibitors and chemotherapies were of the highest grade purchased from Sigma-Aldrich (St. Louis, MO). The NCI Diversity Set VI was used to screen cells against clinically approved drugs for cancer therapy [30]. All chemotherapeutics and small molecule inhibitors were dissolved in DMSO to a stock concentration of 10mM and kept frozen. Parental cell lines were generated as previously described [31] as a clone LCC1. 11, with KrasG12D; p53-/-; Dicer1f/+, where Dicer 1 is heterozygous. Dicer1 mutant clones were generated by transfecting the parental clone with CRISPR plasmid targeting the Dicer1 locus, and selection and expansion of single colonies. Lentiviral particles expressing codon optimized fluorescent proteins under suCMV promoter were transfected into cell lines following manufacturer protocol (GenTarget, San Diego, CA). A blasticidin gene under RSV promoter was used to select positively transduced cells. Despite multiple rounds of selection, there were some noticeable cells ‘negative’ for the fluorescent protein expression, or populations of cells that remained at a low confident level of expression, which were subsequently excluded from analysis in flow cytometry (see flow cytometry section, below). Anonymous human NSCLC tissues were assessed by CANscript using fresh specimen. Fresh tumor tissues were collected immediately after surgical resection. The tumor samples were transported to the laboratory at 4°C, in appropriate transport buffer within 60 minutes post-resection, for ex vivo studies and molecular and pathological evaluations. Tissues were cut into thin sections and cultured in 48-well plate using optimized conditions. Tumors were treated with the indicated drugs at the clinical max concentration (Cmax) for 72 hours. DMSO was used as a vehicle control. Tissue was then formalin fixed and paraffin embedded (FFPE) for subsequent analyses. Quantifying intratumor phenotypic, morphologic, and cellular heterogeneity was performed by visual inspection of a clinical pathologist (Dr. David Goldman, MD, co-author on the present study) using the following methodology: FFPE tissue sections were stained with hematoxylin and eosin (H&E) to identify the respective nuclear DNA content and cytoplasm of cells in the tissue. 1) In an effort to quantify the outgrowth of different tumor ‘clones’, a clinical pathologist then counted the number of histologically distinct tumor ‘neighborhoods’, which were defined as grossly-distinct (clusters of tumor cells growing within a confined region and sharing unique distinguishing morphology from other clusters of tumor cells in the same visual field). 2) Within each tumor ‘neighborhood’, cells were scored based on nuclear density and uniformity as well as cellular morphology uniformity on a scale of 1-5, where 5 is the most distinct and 1 is the most similar. 3) A score for cellular and morphologic heterogeneity was developed by adding together the value attributed to each ‘neighborhood’ for nuclear content uniformity and cellular morphology uniformity (2), which was multiplied by the number of respective tumor neighborhoods (1) for that tissue. Predicting response to therapy was performed using a clinically trained algorithm that was previously described [32]. Briefly, multiple terminal and kinetic assay (tumor morphology, tumor cell proliferation, cell death, viability, cell growth, and metabolic status) inputs were trained in a proprietary machine learning algorithm [32]. The algorithm generates a single score (currently defined as M-Score, but previously published under the nomenclature S-Score) for each drug arm tested. An M-Score > 25 indicates positive response and a value of M-Score ≤ 25 is indicative of a negative response. Cells were cultured as indicated above using the tumor spheroid NanoCulture plates (MBLI, Woburn MA). To assess the total number of cells growing in the 3-D spheroid versus aberrant 2-D growth, cells were plated and counted every other day for 7 days. Cells were imaged by brightfield microscope and 2-D adherent were quantified. Number of 2-D growing cells were quantified as % of total cells growing in culture. Less than 2% of cells on any day were noted to grow in 2-D versus the majority (98%-99% of cells) growing in the 3-D tumor spheroids. Cells were cultured in various proportions with either the parent WT cell line or any one of a combination of the different Dicer1 mutants. Cells were removed from culture dishes with StemPro Accutase dissociation reagent (Invitrogen, Carlsbad CA) and fixed with 4% paraformaldehyde in PBS for 30 minutes at RT. Cells were then processed using the BD Fortessa 4 laser, 17 parameter flow cytometer. Measurements were made in the indicated fluorescent channel based on the fluorescent tag (typically blue fluorescent protein (BFP) and red fluorescent protein (RFP) ), and any overlap was compensated prior to analysis. A gating strategy was employed to select cells based on the side scatter (SSC) and forward scatter (FSC) from the vehicle-treated parental population, followed by a selection of singlets based on FSC width and height. Propidium iodide exclusion was used to validate viable cells are contained within the FSC: SSC gate. An equal volume of cells was analyzed for each experiment and events were recorded in the defined gates (described above) and in the correct fluorescent channels. Data analysis was performed using FlowJo software (Tree Star Inc. , Ashland OR). Experiments were performed a minimum of three times (biological replicates) on independent days. Notably, and despite blasticidin selection, there were identifiable populations of cells that expressed no detectable fluorescent tag. This was attributed to heterogeneity of both CMV infection, expression and blasticidin selection. During analysis in flow cytometry, there was no increase, decrease or change in the proportion of labeled to unlabeled cells in the vehicle treatment or drug treatment cohorts. RNA was isolated from cells that were cultured in 3-D culture after 48 hours of growth in spheroids. Cells were rinsed and then lysed followed by RNA extraction using manufacturer protocol (Qiagen, Hiden Germany). RNA libraries were prepared using Illumina TruSeq Stranded mRNA sample preparation kits from 500ng of purified total RNA according to the manufacturer’s protocol. The resultant dsDNA libraries were quantified by Qubit fluorometer, Agilent TapeStation 2200, and RT-qPCR using the Kapa Biosystems library quantification kit according to manufacturer’s protocols. Uniquely indexed libraries were pooled in equimolar ratios and sequenced on a single Illumina NextSeq500 run with single-end 75bp reads by the Dana-Farber Cancer Institute Molecular Biology Core Facilities. Sequenced reads were aligned to the UCSC mm9 reference genome assembly and gene counts were quantified using STAR (v2. 5. 1b) [33]. Differential gene expression testing was performed by DESeq2 (v1. 10. 1) [34] and normalized read counts (FPKM) were calculated using cufflinks [35]. RNASeq data was mapped to a reference obtained from NCBI GRCh38. p12. Indexing of the raw sequence data contained in GRCH38. p12 was performed by the STAR open source alignment tool. These services were run on the Amazon Web Service (AWS) Elastic Cloud Compute (EC2) infrastructure with 8 vCPU cores and 32 GB Ram. Mapped results were correlated to index files using STAR 2pass. Indels were realigned and bases recalibrated before variant calling by STAR. We utilized the HaplotypeCaller as outlined in GATK best practices with the Java implementation. Variant calls (. vcf) were generated and annotated back to the GRCh38. p12 genomic reference assembly in order to determine point mutations. The mathematical model of adaptive resistance consists of two ordinary differential equations describing intra- and inter-species phenotypic switching founded upon game theoretic interaction assumptions identical to the Lotka-Volterra model commonly used in ecology [36,37]. This choice is founded upon well-developed mathematical modelling approaches to understanding, for example, the dynamics of interacting microbes and cancer cell populations [38,39]. Full details of the model are provided in the Supplementary Information file S1 SuppInfo. Briefly, we consider a wild xWT and mutant x M i (i = 1,2, 3) type. Let b be the birth rate of cells and d be their rate of death. In absence of drugs, cells initially grow exponentially at rate ltype (type = WT, Mi), where ltype is the effective birthrate, or fitness, given by ltype = b − d. Cell numbers eventually saturate at Ktype (the population’s effective carrying capacity, where K = K ˜ (1 - d / b) with K ˜ being the standard carrying capacity–see S1 SuppInfo for extended details). Based on the observation of a temporary decline and eventual rebound in the monoculture growth assays after drug pressure was applied, we assumed each cell type is capable of developing an intra-species drug tolerant phenotype under drug pressure [40,41] that we termed ‘intra-species phenotypic switching’. The intra-species growth dynamics can then be summarized by d x t y p e s e n s d t = l t y p e s e n s x t y p e s e n s (1 - (x t y p e s e n s + x t y p e t o l) K t y p e) - ν t y p e x t y p e s e n s d x t y p e t o l d t = l t y p e t o l x t y p e t o l (1 - (x t y p e s e n s + x t y p e t o l) K t y p e) + ν t y p e x t y p e s e n s, (1) where the superscript sens describes the drug-sensitive phenotype, and tol the drug-tolerant phenotype; νtype measures the intra-species switching rate. Note that no backwards switching from tolerant to sensitive is considered. We assume that intra-species interactions are dominated by inter-species relationships, that is to say that the adaptive behaviors of cells in co-culture are dictated by CAT dynamics. Co-culture dynamics were modelled as d x WT d t = l 1 x WT (1 - x WT K W T c o + a 12 x M i) d x M i d t = l 2 x M i (1 - x M i K M i c o + a 21 x WT), (2) where x WT = x WT s e n s + x WT t o l, x M i = x M i s e n s + x M i t o l, K W T c o and K M i c o are the carrying capacities of the WT and mutant in co-culture, and the terms a12 and a21 denote the interactions between cell types. In summary, we assumed the following: The mathematical model’s parameter values were estimated from the experimental procedures described in S1 SuppInfo. To perform the optimization, we used the Matlab function fmincon that minimizes the objective (cost) function of a nonlinear model with constraints [42]. To reduce the influence of initial conditions on the outcome, goodness-of-fit was assessed by minimizing over the resulting set of cost function values obtained via a multi-step parameter estimation procedure, as described in the Supplementary Information file S1 SuppInfo.
To study the role of intratumor heterogeneity in therapy failure, we first derived and cloned a cell line [43]. This cell line possesses an oncogenic Kras in conjunction with homozygous p53 and heterozygous Dicer1 loss of function (Kras G12D/+; p53 /; Dicer1 +/) and is capable of inducing tumors when transplanted into immunocompromised mice. Because Dicer1 regulates the production of miRNA that can dramatically influence the phenotypic heterogeneity of a cell [44,45], we engineered clones using the CRISPR/Cas9 toolkit to introduce a variety of mutations in Dicer1 (Fig 1a). Confirmation of Dicer1 sequence mutations were determined using NCBI reference data, identifying insertions and deletions as well as 17 unique single nucleotide polymorphisms (SNPs) in the mutant cell lines compared to WT (Fig 1b and S1 Fig). The mutants and WT cell lines were cultured in a physiologically-relevant manner in vitro, using a three-dimensional nano-culture platform that provided a substrate-free growth surface containing uniform, nanofabricated imperfections, which forced spheroid growth (NanoCulture plates –Fig 1c). NanoCulture plates enabled better long term growth potential than ultra low adherence (ULA) plates, and growth of spheroids were also more uniform (S2 Fig). Using this approach, tumor spheroid growth is not confounded by external factors such as matrigel or collagen. We characterized transcriptional heterogeneity of the wild type (WT) and mutant (M1-3) cells via interrogation of the mRNA sequencing from 3D cell cultures (Fig 1c). Primarily, differential expression analysis confirmed Dicer1 target-gene expressions (i. e. miRNA) were significantly, statistically different among mutants compared to WT with varying expression profiles (Fig 1d). Differential expression profiling of the top and bottom 500 dysregulated genes in each cell line confirmed unique transcriptional heterogeneity was observed in each mutant cell line with many reaching statistical significance (p-adjusted < 0. 05) (see S1 Table). In order to study how clonal diversity might impact response to drugs (sensitivity or resistance) we grew Dicer1 mutant clones and parental wild-type in separate mono-cultures or mixed co-cultures in 3-D and exposed them to different drugs at lethal concentrations. Using lentiviral-transfected fluorescent ‘tags’, flow cytometry was deployed to quantify the abundance of the different mutant clones in dual co-cultures, as described in the methods section (Fig 2a). We cultured the pre-labeled clones alone or together with WT in 3-D culture at a ratio of 50: 50 in the presence or absence of multiple clinically-relevant anticancer drugs: docetaxel, bortezomib, and afatinib [46]. Primarily, we determined that, while all cell lines showed heterogeneous sensitivities to drugs in mono-cultures as indicated by the proportion of tumor cells in spheroids after 96 hours of drug exposure, a mutant-WT mixed co-culture resulted in an increased proportion of each clone relative to the vehicle control (Fig 2b). The data indicated that both cell lines ‘benefitted’ in a mixed co-culture (WT or mutant), increasing their relative ‘resistance’ to drug compared to the same cell line’s mono-culture. This finding was consistent in each drug tested and in every mixed culture, with the exception of the M3/WT co-culture exposed to afatinib. Next, we sought a method to study the growth dynamics of the phenotypically diverse mutant clones in mono- and co-culture with parental WT cells in a competitive growth assay. We cultured the pre-labeled clones alone or together with WT in 3-D culture at different population density ratios of 50: 50,90: 10, or 10: 90. Using flow cytometry, we then performed two measurements at 48 hour intervals over the course of 7 days: 1) the % and number of cells in the tumor sphere population that is either WT or Mu, and 2) the size of the tumor sphere based on the number of fluorescently-detected, gated cells (Fig 3), which is indicative of overall tumor sphere growth. These data were integrated into the subsequent math modeling experiments and computational analyses described in the following sections. We leveraged the mathematical model to better understand the growth dynamics in the in vitro competitive growth assays. We first characterized monoculture dynamics (with and without drug) for the WT and all three mutants. In the absence of drug pressure, cells initially grew exponentially in monoculture (Fig 4a and 4b, top left) with no suggestion of intra-species phenotypic switching (“No Drug” columns in Table 1). However, the ‘dip and rebound’ observed in the monoculture growth assays after the addition of docetaxel, bortezomib, or afatinib was successfully captured through the monoculture model in (Eq 1) (Fig 4a and 4b top right, bottom; Drug columns in Table 1—see S1 SuppInfo). Our estimates predicted decreased growth of the sensitive phenotype and increased growth of the tolerant subtype in the presence of drugs in monoculture. The complete results of the monoculture parameter estimation procedure are reported in Table 1. As previously described, we hypothesized that genotype-genotype interactions are the dominant growth mechanism in co-culture (see Eq (2) ). This assumption was borne out as the model successfully captures the WT and mutant growth dynamics in co-culture mixes (10: 90,50: 50,90: 10 WT: Mu), as shown for M1 in Fig 4, and for fixed proportion (50: 50) across all mutant co-culture mixes (S2 Table and S4 Fig), with the exception of the previously identified M3/WT 50: 50 co-culture exposed to afatinib. Next, we tested the hypothesis that cancer cells induce rapid adaptations within and among neighboring cells to improve fitness of the heterogeneous population under drug pressure and thwart destruction. That is to say, cancer cells behave in a ‘cooperative’ manner to promote drug resistance in neighboring cells that do not have similar genetic and phenotypic features, hereafter referred to as cooperative adaptation to therapy (CAT). We hypothesized that the form of ‘adaptive’ resistance describing monoculture growth does not impact on the overall co-culture dynamics (i. e. the additional co-culture effects induce a change in the fitness values due to genotype-genotype interactions in a cooperative game) [47]. Therefore, we leveraged the parameterized mathematical model for comparative analysis. For the four drug scenarios (no drug, docetaxel, bortezomib, and afatanib), we compared monoculture growth dynamics to those in all three co-culture proportions by simulating all scenarios with identical initial condition (xWT (0) = 1, x M i (0) = 1) for the WT and M1 mutant. In the absence of drug exposure, for both cell types, we found that growth in monoculture is demonstrably stronger than in co-culture (Fig 5 left, top and bottom). However, under drug pressure, growth in co-culture quickly outpaces that of the monoculture (Fig 5 docetaxel, bortezomib, and afatanib panels, top and bottom). These results demonstrate mathematically that our model, parameterized to our in vitro NanoCulture growth assays, predicts that an advantageous drug tolerant phenotype, due only to heterogeneity, is conferred in co-culture. Indeed, these findings are consistent with the CAT hypothesis. The data produced to this point indicate that heterogenous tumors can overcome the onslaught of anti-cancer drugs by rapidly inducing resistance through co-supportive effects (i. e. CAT). We tested the hypothesis that ITH directly correlates to therapy response of many drugs, not just those predisposed by mutations. Indeed, ITH is a consequence of genetic and phenotypic (i. e. mutations giving rise to unique cellular and morphologic characteristics) variability among tumor cells, which is tightly associated with treatment resistance [48]. To do this, we deployed a human tumor explant platform that replicates the native tumor microenvironment including stromal and immune cells and allows for testing multiple drug responses in a single tumor explant [32]. First, tumor explants were generated from fresh biopsies isolated from NSCLC patients, cultured ex vivo on matched tumor matrix proteins supplemented with autologous patient plasma, and exposed to various clinically-relevant drug regimens (Fig 6a and S3 Table). This platform uses a trained algorithm that incorporates phenotypic response assays, which enables clinical response prediction (given by the M-Score) with published clinical accuracy (Fig 6b) [32]. Next, we wanted to quantify phenotypic heterogeneity by assessing the morphologic and cellular diversity of tumor cells using an assessment of clinical pathology. i. e. phenotypic rather than genetic heterogeneity. To do this, we developed a ‘heterogeneity score’ for each patient tumor sample based on the cellular and morphologic diversity of the tumor as quantified by a clinical pathologist using the following approach: Fresh tumor biopsies were stained with hematoxylin and eosin to visualize and differentiate the nucleus and cytoplasm of each cell. Next, a clinical pathologist (see Methods) quantified the number of distinct tumor ‘neighborhoods’, which are defined by clusters of tumor cells growing with similar morphologies within a single visual field at 20X magnification. This was performed on three independent tumor fragments from the same patient. Within each ‘neighborhood’ we then quantified, on a scale of 1-5, the uniformity of cell morphology and nuclear content in each visual field wherein 1 is more uniform and 5 is highly diverse. We then calculated the ‘heterogeneity score’ by multiplying the number of histologically distinct tumor neighborhoods from a single field-of-view with cellular morphology and nuclear uniformity, each on a scale of 1-5 (Fig 6c). The resulting score was used to stratify NSCLC patient samples across a spectrum of heterogeneity from low to high (Fig 6d). Interestingly, and consistent with evidence that ITH associates with therapy failure [48], we determined that the degree of ITH (as defined by our ‘heterogeneity score’) significantly affects predicted antitumor effect of multiple anticancer therapeutic regimens and combinations in an individual patient sample, as determined by the M-Score (Fig 6e, p < 0. 001). While this finding is compatible with previous evidence that heterogeneity across a population of patients associates to therapy failure [49,50], these data provide evidence that in an ex vivo setting, within a single patient, numerous conventional treatment options are similarly ineffective when ITH is high, supporting the hypothesis that ITH and CAT are potentially linked.
Identifying the biomarkers and patterns that result in drug resistance is penultimate to eradicating therapy failure in cancer. To achieve this and to advance precision medicine, mutational profiling for ‘druggable’ targets has been an active area of investigation for oncology in the past several decades [51]. However, there is mounting evidence that non-mutational mechanisms result in drug resistance, regardless of the mutational load in tumors [52]. Indeed, we have outlined here a potentially new paradigm in drug resistance by showing how “cooperative adaptation to therapy” (CAT), due to ITH, can induce a drug tolerant phenotype, which has significant consequences for therapeutic response. We took an approach that combined experimental evidence with theoretical modeling, and constructed a replicator dynamics model of intra- and interspecies competition that can potentially explain how individual clones within a tumor lead to treatment resistance. We determined that the emergence of CAT is directly related to the presence of neighboring clonal subsets within a tumor that force de-novo drug resistance after drug exposure. We investigated how phenotypic switching can occur within heterogenous tumors by leveraging a three-dimensional nanoculture in vitro spherical growth platform and Dicer1 mutants derived from CRISPR/Cas9 gene editing. We characterized tumor spheroid growth over 7 days in monocultures and co-cultures (in proportions of 10: 90,50: 50, and 90: 10 WT: Mutant) in the absence of drugs and in the presence of docetaxel, bortezomib, and afatinib as unique drug classes with differing mechanisms of action. We observed persistent growth in co-cultures under drug pressure, despite previous drug sensitivity of both the WT and the mutant, demonstrating the ability for phenotypic switching in heterogenous tumors. To quantify this behavior, we developed a mathematical model of the growth dynamics in mono- and co-cultures. We assumed that previously sensitive types exhibit an increased drug tolerance in the presence of drugs through non-linear interactions. Our model successfully characterized the phenotypic switching and our results demonstrated how genotype/genotype interactions promote increased tolerance to drugs. The empirical data presented in Fig 2 demonstrated that the mixed, heterotypic co-culture conditions resulted in improved survival of both cell lines in the mixed culture (with the exception of M3 in afatinib). Indeed, in the case of docetaxel, such behavior could be attributed to decreased cell cycling or other mechanisms that would argue against CAT as a method of thwarting drug pressure. However, the evidence for mutually beneficial growth of both cell lines, relative to vehicle control, when exposed to multiple different classes of drugs with different mechanisms of action, some of which don’t rely on cell cycling or proliferation (e. g. bortezomib), support the hypothesis of CAT and indicates that it may explain a drug-agnostic, or more universal phenomenon. While CAT is one possible explanation for the behavior observed here, there are other potential mechanisms also at play. Given these surprising empirical data, mathematical modeling of game theoretic cooperation was employed to align theoretic and empirical evidences. Indeed, the numerous instances in which the in silico model fit the experimental data provide support for the overall hypothesis. There is more work to be done, and understanding the mechanisms that contribute to this phenomenon still need further investigation. Leveraging the ex vivo human-autologous explant platform, our results suggested that ITH affects a patient’s sensitivity to multiple anticancer drugs predominated by kinase specific as well as general cytotoxics. In the context of our findings, this means that, regardless of druggable targets, CAT may confer universal drug resistance behavior. Therefore, better-informed therapeutic interventions should be considered such as timing the sequence of drugs [9], or using combinations of rational agents based on computational modeling [53]. Indeed, many complex alternative therapeutic interventions exist, and are being tested, which could potentially address some of the challenges of CAT. While the present study focused on theoretical models, examining the transcriptional and proteomic profile of cells in co-culture could result in novel therapeutic targets to combine with therapy and thwart CAT. Future studies might focus on growth dynamics in models of complex heterogeneity that include the microenvironment, immune cells and stromal cells. This could paint a clear picture of how CAT influences drug response. Given the evidences presented here, using strategies that can inform mutational evolution and provide single cell transcriptional profiling should be applied in parallel to a computational effort to gain a complete picture of cell-cell interactions and help guide therapeutic options for patients receiving care.
Here, we provide mathematical and empirical evidence to support a potentially new paradigm in drug resistance, which we have termed “cooperative adaptation to therapy” (CAT). CAT is defined by a phenomenon wherein drug-sensitive cancer cells with different genetic and phenotypic features within a 3-dimensional heterogeneous tumor induce non-mutational resistance in their neighboring cells under pressure of cancer therapy. To develop this novel conclusion we deployed an interdisciplinary effort including an ex vivo human tumor model, a CRISPR/Cas9 platform with 3-dimensional in vitro experiments, and high throughput flow cytometry. Importantly, we wove these data together using a mathematical model of intra- and interspecies competition to understand how tumor heterogeneity influenced our observations. By estimating the model’s parameters, we determined that the combination of genetic clonal variants in co-cultures allowed for previously drug-sensitive cells to continue to grow despite drug pressure. We were thus able to characterize distinct growth regimens in mono- and co-cultures without and with drug pressure.
flow cytometry
sequencing techniques
medicine and health sciences
cancer treatment
cancers and neoplasms
social sciences
cloning
mathematical models
oncology
animal behavior
molecular biology techniques
rna sequencing
zoology
research and analysis methods
behavior
lung and intrathoracic tumors
mathematical and statistical techniques
molecular biology
spectrophotometry
cytophotometry
psychology
phenotypes
genetics
biology and life sciences
non-small cell lung cancer
spectrum analysis techniques
Cooperative adaptation to therapy (CAT) confers resistance in heterogeneous non-small cell lung cancer