Abstract 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.

Author summary 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.

Citation: Craig M, Kaveh K, Woosley A, Brown AS, Goldman D, Eton E, et al. (2019) Cooperative adaptation to therapy (CAT) confers resistance in heterogeneous non-small cell lung cancer. PLoS Comput Biol 15(8): e1007278. https://doi.org/10.1371/journal.pcbi.1007278 Editor: Darragh Gerard McArt, Queen’s University Belfast, UNITED KINGDOM Received: January 8, 2019; Accepted: July 22, 2019; Published: August 26, 2019 Copyright: © 2019 Craig et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability: The raw sequence data are available from the NCBI’s Sequence Read Archive, accession number GSE129353 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc = GSE129353). Average growth and standard deviations of the mono- and co-culture growth assays to which the model was fit have been included as a supplementary Excel file in our resubmission (S1 File). Funding: MC was funded by Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grant RGPIN-2018-04546, an NSERC Postdoctoral Fellowship, and National Institutes of Health grant DP5OD019851. AG was funded by a Breast Cancer Alliance Young Investigator Award. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: I have read the journal’s policy and the authors of this manuscript have the following competing interests: Aaron Goldman is an employee of Mitra Biotech and holds equity. Kazuya Arai is employed by JSR Life Sciences Corporation. M. Mamunur Rahman is employed by MBLI Inc. All other authors have declared that no competing interests exist.

Introduction 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.

Methods Ethics statement Anonymous non-small cell lung cancer (NSCLC) tissue samples were collected under IRB approval with due written consent from each patient. Chemicals and reagents 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. Cell culture 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). Human explant studies 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. Flow cytometry analyses 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 sequencing 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). Library preparation and sequencing 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. RNAseq analysis 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. Modeling the population dynamics of adaptive resistance in mono- and co-cultures 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 x WT and mutant (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 l type (type = WT, M i ), where l type is the effective birthrate, or fitness, given by l type = b − d. Cell numbers eventually saturate at K type (the population’s effective carrying capacity, where with 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 (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 (2) where , , and are the carrying capacities of the WT and mutant in co-culture, and the terms a 12 and a 21 denote the interactions between cell types. In summary, we assumed the following: In absence of drugs, monoculture growth is initially exponential and eventually saturates at an effective carrying capacity specific to each cell type (logistic growth). Mathematically, we assume that Eq 1) are equal to 0.

Eq 1) are equal to 0. Under exposure to drugs, sensitive cells in monoculture can perform intra-species switches to drug tolerant type at rate ν type . The effective carrying capacity K type remains unchanged from the non-drug case, however

. The effective carrying capacity K remains unchanged from the non-drug case, however In the presence of the therapeutic stresses, sensitive subtypes phenotypically switch into more drug tolerant/resistant types. Further, each phenotypic trait is heritable, that is that daughter cells have the same phenotype.

Intra-species growth dynamics are weaker than inter-species co-operative adaptation to therapy, so that the inter-species interaction terms a 12 and a 21 dominate. This implies that the additional co-culture effects are not changing or modifying the intra-species phenotypic switching mechanism but rather inducing a game theoretic change in the fitness values due to genotype-genotype interactions. Mathematically, we hypothesized that, as x WT and Eq 2).

and a dominate. This implies that the additional co-culture effects are not changing or modifying the intra-species phenotypic switching mechanism but rather inducing a game theoretic change in the fitness values due to genotype-genotype interactions. Mathematically, we hypothesized that, as x and Eq 2). The proliferation potentials and apoptosis rates of the sensitive and drug-tolerant subtypes are dose-dependent (note, however, that since only a single dose was considered experimentally, we fixed the dose and did not elaborate the dose-response function of each genotype). Adaptive resistance model parameter optimization 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.

Discussion 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.

Acknowledgments The authors would like to thank Dr. Jeff Gore and Dr. Kirill Korolev for their helpful discussions during the concept development phase of the project and discussion during data acquisition. The authors would like to thank Dr. Feng Zhang for help with CRISPR/Cas9 technology and contribution of resources.