Reconstruction of single-cell lineage trajectories and identification of diversity in fates during the epithelial-to-mesenchymal transition.

Exploring the complexity of the epithelial-to-mesenchymal transition (EMT) unveils a diversity of potential cell fates; however, the exact timing and mechanisms by which early cell states diverge into distinct EMT trajectories remain unclear. Studying these EMT trajectories through single-cell RNA sequencing is challenging due to the necessity of sacrificing cells for each measurement. In this study, we employed optimal-transport analysis to reconstruct the past trajectories of different cell fates during TGF-beta-induced EMT in the MCF10A cell line. Our analysis revealed three distinct trajectories leading to low EMT, partial EMT, and high EMT states. Cells along the partial EMT trajectory showed substantial variations in the EMT signature and exhibited pronounced stemness. Throughout this EMT trajectory, we observed a consistent downregulation of the EED and EZH2 genes. This finding was validated by recent inhibitor screens of EMT regulators and CRISPR screen studies. Moreover, we applied our analysis of early-phase differential gene expression to gene sets associated with stemness and proliferation, pinpointing ITGB4, LAMA3, and LAMB3 as genes differentially expressed in the initial stages of the partial versus high EMT trajectories. We also found that CENPF, CKS1B, and MKI67 showed significant upregulation in the high EMT trajectory. While the first group of genes aligns with findings from previous studies, our work uniquely pinpoints the precise timing of these upregulations. Finally, the identification of the latter group of genes sheds light on potential cell cycle targets for modulating EMT trajectories.

Mathematical Modeling Identifies Optimum Palbociclib-fulvestrant Dose Administration Schedules for the Treatment of Patients with Estrogen Receptor-positive Breast Cancer.

Cyclin-dependent kinases 4/6 (CDK4/6) inhibitors such as palbociclib are approved for the treatment of metastatic estrogen receptor-positive (ER+) breast cancer in combination with endocrine therapies and significantly improve outcomes in patients with this disease. However, given the large number of possible pairwise drug combinations and administration schedules, it remains unclear which clinical strategy would lead to best survival. Here, we developed a computational, cell cycle-explicit model to characterize the pharmacodynamic response to palbociclib-fulvestrant combination therapy. This pharmacodynamic model was parameterized, in a Bayesian statistical inference approach, using in vitro data from cells with wild-type estrogen receptor (WT-ER) and cells expressing the activating missense ER mutation, Y537S, which confers resistance to fulvestrant. We then incorporated pharmacokinetic models derived from clinical data into our computational modeling platform. To systematically compare dose administration schedules, we performed in silico clinical trials based on integrating our pharmacodynamic and pharmacokinetic models as well as considering clinical toxicity constraints. We found that continuous dosing of palbociclib is more effective for lowering overall tumor burden than the standard, pulsed-dose palbociclib treatment. Importantly, our mathematical modeling and statistical analysis platform provides a rational method for comparing treatment strategies in search of optimal combination dosing strategies of other cell-cycle inhibitors in ER+ breast cancer. SIGNIFICANCE: We created a computational modeling platform to predict the effects of fulvestrant/palbocilib treatment on WT-ER and Y537S-mutant breast cancer cells, and found that continuous treatment schedules are more effective than the standard, pulsed-dose palbociclib treatment schedule.

Single-cell heterogeneity of EGFR and CDK4 co-amplification is linked to immune infiltration in glioblastoma.

Glioblastoma (GBM) is the most aggressive brain tumor, with a median survival of ∼15 months. Targeted approaches have not been successful in this tumor type due to the large extent of intratumor heterogeneity. Mosaic amplification of oncogenes suggests that multiple genetically distinct clones are present in each tumor. To uncover the relationships between genetically diverse subpopulations of GBM cells and their native tumor microenvironment, we employ highly multiplexed spatial protein profiling coupled with single-cell spatial mapping of fluorescence in situ hybridization (FISH) for EGFR, CDK4, and PDGFRA. Single-cell FISH analysis of a total of 35,843 single nuclei reveals that tumors in which amplifications of EGFR and CDK4 more frequently co-occur in the same cell exhibit higher infiltration of CD163+ immunosuppressive macrophages. Our results suggest that high-throughput assessment of genomic alterations at the single-cell level could provide a measure for predicting the immune state of GBM.

Identification of optimal dosing schedules of dacomitinib and osimertinib for a phase I/II trial in advanced EGFR-mutant non-small cell lung cancer.

Despite the clinical success of the third-generation EGFR inhibitor osimertinib as a first-line treatment of EGFR-mutant non-small cell lung cancer (NSCLC), resistance arises due to the acquisition of EGFR second-site mutations and other mechanisms, which necessitates alternative therapies. Dacomitinib, a pan-HER inhibitor, is approved for first-line treatment and results in different acquired EGFR mutations than osimertinib that mediate on-target resistance. A combination of osimertinib and dacomitinib could therefore induce more durable responses by preventing the emergence of resistance. Here we present an integrated computational modeling and experimental approach to identify an optimal dosing schedule for osimertinib and dacomitinib combination therapy. We developed a predictive model that encompasses tumor heterogeneity and inter-subject pharmacokinetic variability to predict tumor evolution under different dosing schedules, parameterized using in vitro dose-response data. This model was validated using cell line data and used to identify an optimal combination dosing schedule. Our schedule was subsequently confirmed tolerable in an ongoing dose-escalation phase I clinical trial (NCT03810807), with some dose modifications, demonstrating that our rational modeling approach can be used to identify appropriate dosing for combination therapy in the clinical setting.

Breast tumours maintain a reservoir of subclonal diversity during expansion.

Our knowledge of copy number evolution during the expansion of primary breast tumours is limited1,2. Here, to investigate this process, we developed a single-cell, single-molecule DNA-sequencing method and performed copy number analysis of 16,178 single cells from 8 human triple-negative breast cancers and 4 cell lines. The results show that breast tumours and cell lines comprise a large milieu of subclones (7-22) that are organized into a few (3-5) major superclones. Evolutionary analysis suggests that after clonal TP53 mutations, multiple loss-of-heterozygosity events and genome doubling, there was a period of transient genomic instability followed by ongoing copy number evolution during the primary tumour expansion. By subcloning single daughter cells in culture, we show that tumour cells rediversify their genomes and do not retain isogenic properties. These data show that triple-negative breast cancers continue to evolve chromosome aberrations and maintain a reservoir of subclonal diversity during primary tumour growth.

Imaging dynamic mTORC1 pathway activity in vivo reveals marked shifts that support time-specific inhibitor therapy in AML.

Acute myeloid leukemia (AML) is a high remission, high relapse fatal blood cancer. Although mTORC1 is a master regulator of cell proliferation and survival, its inhibitors have not performed well as AML treatments. To uncover the dynamics of mTORC1 activity in vivo, fluorescent probes are developed to track single cell proliferation, apoptosis and mTORC1 activity of AML cells in the bone marrow of live animals and to quantify these activities in the context of microanatomical localization and intra-tumoral heterogeneity. When chemotherapy drugs commonly used clinically are given to mice with AML, apoptosis is rapid, diffuse and not preferentially restricted to anatomic sites. Dynamic measurement of mTORC1 activity indicated a decline in mTORC1 activity with AML progression. However, at the time of maximal chemotherapy response, mTORC1 signaling is high and positively correlated with a leukemia stemness transcriptional profile. Cell barcoding reveals the induction of mTORC1 activity rather than selection of mTORC1 high cells and timed inhibition of mTORC1 improved the killing of AML cells. These data define the real-time dynamics of AML and the mTORC1 pathway in association with AML growth, response to and relapse after chemotherapy. They provide guidance for timed intervention with pathway-specific inhibitors.

Beyond Deterministic Models in Drug Discovery and Development.

The model-informed drug discovery and development paradigm is now well established among the pharmaceutical industry and regulatory agencies. This success has been mainly due to the ability of pharmacometrics to bring together different modeling strategies, such as population pharmacokinetics/pharmacodynamics (PK/PD) and systems biology/pharmacology. However, there are promising quantitative approaches that are still seldom used by pharmacometricians and that deserve consideration. One such case is the stochastic modeling approach, which can be important when modeling small populations because random events can have a huge impact on these systems. In this review, we aim to raise awareness of stochastic models and how to combine them with existing modeling techniques, with the ultimate goal of making future drug-disease models more versatile and realistic.

Pharmacokinetic Profiles Determine Optimal Combination Treatment Schedules in Computational Models of Drug Resistance.

Identification of optimal schedules for combination drug administration relies on accurately estimating the correct pharmacokinetics, pharmacodynamics, and drug interaction effects. Misspecification of pharmacokinetics can lead to wrongly predicted timing or order of treatments, leading to schedules recommended on the basis of incorrect assumptions about absorption and elimination of a drug and its effect on tumor growth. Here, we developed a computational modeling platform and software package for combination treatment strategies with flexible pharmacokinetic profiles and multidrug interaction curves that are estimated from data. The software can be used to compare prespecified schedules on the basis of the number of resistant cells where drug interactions and pharmacokinetic curves can be estimated from user-provided data or models. We applied our approach to publicly available in vitro data of treatment with different tyrosine kinase inhibitors of BT-20 triple-negative breast cancer cells and of treatment with erlotinib of PC-9 non-small cell lung cancer cells. Our approach is publicly available in the form of an R package called ACESO (https://github.com/Michorlab/aceso) and can be used to investigate optimum dosing for any combination treatment. SIGNIFICANCE: These findings introduce a computational modeling platform and software package for combination treatment strategies with flexible pharmacokinetic profiles and multidrug interaction curves that are estimated from data.

Acquired resistance to combined BET and CDK4/6 inhibition in triple-negative breast cancer.

BET inhibitors are promising therapeutic agents for the treatment of triple-negative breast cancer (TNBC), but the rapid emergence of resistance necessitates investigation of combination therapies and their effects on tumor evolution. Here, we show that palbociclib, a CDK4/6 inhibitor, and paclitaxel, a microtubule inhibitor, synergize with the BET inhibitor JQ1 in TNBC lines. High-complexity DNA barcoding and mathematical modeling indicate a high rate of de novo acquired resistance to these drugs relative to pre-existing resistance. We demonstrate that the combination of JQ1 and palbociclib induces cell division errors, which can increase the chance of developing aneuploidy. Characterizing acquired resistance to combination treatment at a single cell level shows heterogeneous mechanisms including activation of G1-S and senescence pathways. Our results establish a rationale for further investigation of combined BET and CDK4/6 inhibition in TNBC and suggest novel mechanisms of action for these drugs and new vulnerabilities in cells after emergence of resistance.

ESTIpop: a computational tool to simulate and estimate parameters for continuous-time Markov branching processes.

SUMMARY: ESTIpop is an R package designed to simulate and estimate parameters for continuous-time Markov branching processes with constant or time-dependent rates, a common model for asexually reproducing cell populations. Analytical approaches to parameter estimation quickly become intractable in complex branching processes. In ESTIpop, parameter estimation is based on a likelihood function with respect to a time series of cell counts, approximated by the Central Limit Theorem for multitype branching processes. Additionally, simulation in ESTIpop via approximation can be performed many times faster than exact simulation methods with similar results. AVAILABILITY AND IMPLEMENTATION: ESTIpop is available as an R package on Github (https://github.com/michorlab/estipop). SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Synthetic Lethal and Resistance Interactions with BET Bromodomain Inhibitors in Triple-Negative Breast Cancer.

BET bromodomain inhibitors (BBDIs) are candidate therapeutic agents for triple-negative breast cancer (TNBC) and other cancer types, but inherent and acquired resistance to BBDIs limits their potential clinical use. Using CRISPR and small-molecule inhibitor screens combined with comprehensive molecular profiling of BBDI response and resistance, we identified synthetic lethal interactions with BBDIs and genes that, when deleted, confer resistance. We observed synergy with regulators of cell cycle progression, YAP, AXL, and SRC signaling, and chemotherapeutic agents. We also uncovered functional similarities and differences among BRD2, BRD4, and BRD7. Although deletion of BRD2 enhances sensitivity to BBDIs, BRD7 loss leads to gain of TEAD-YAP chromatin binding and luminal features associated with BBDI resistance. Single-cell RNA-seq, ATAC-seq, and cellular barcoding analysis of BBDI responses in sensitive and resistant cell lines highlight significant heterogeneity among samples and demonstrate that BBDI resistance can be pre-existing or acquired.

DIFFpop: a stochastic computational approach to simulate differentiation hierarchies with single cell barcoding.

SUMMARY: DIFFpop is an R package designed to simulate cellular differentiation hierarchies using either exponentially-expanding or fixed population sizes. The software includes functionalities to simulate clonal evolution due to the emergence of driver mutations under the infinite-allele assumption as well as options for simulation and analysis of single cell barcoding and labeling data. The software uses the Gillespie Stochastic Simulation Algorithm and a modification of expanding or fixed-size stochastic process models expanded to a large number of cell types and scenarios. AVAILABILITY AND IMPLEMENTATION: DIFFpop is available as an R-package along with vignettes on Github (https://github.com/ferlicjl/diffpop). SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

KDM5 Histone Demethylase Activity Links Cellular Transcriptomic Heterogeneity to Therapeutic Resistance.

Members of the KDM5 histone H3 lysine 4 demethylase family are associated with therapeutic resistance, including endocrine resistance in breast cancer, but the underlying mechanism is poorly defined. Here we show that genetic deletion of KDM5A/B or inhibition of KDM5 activity increases sensitivity to anti-estrogens by modulating estrogen receptor (ER) signaling and by decreasing cellular transcriptomic heterogeneity. Higher KDM5B expression levels are associated with higher transcriptomic heterogeneity and poor prognosis in ER+ breast tumors. Single-cell RNA sequencing, cellular barcoding, and mathematical modeling demonstrate that endocrine resistance is due to selection for pre-existing genetically distinct cells, while KDM5 inhibitor resistance is acquired. Our findings highlight the importance of cellular phenotypic heterogeneity in therapeutic resistance and identify KDM5A/B as key regulators of this process.

SIApopr: a computational method to simulate evolutionary branching trees for analysis of tumor clonal evolution.

SUMMARY: SIApopr (Simulating Infinite-Allele populations) is an R package to simulate time-homogeneous and inhomogeneous stochastic branching processes under a very flexible set of assumptions using the speed of C ++. The software simulates clonal evolution with the emergence of driver and passenger mutations under the infinite-allele assumption. The software is an application of the Gillespie Stochastic Simulation Algorithm expanded to a large number of cell types and scenarios, with the intention of allowing users to easily modify existing models or create their own. AVAILABILITY AND IMPLEMENTATION: SIApopr is available as an R library on Github ( https://github.com/olliemcdonald/siapopr ). SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online. CONTACT: michor@jimmy.harvard.edu.

Punctuated copy number evolution and clonal stasis in triple-negative breast cancer.

Aneuploidy is a hallmark of breast cancer; however, knowledge of how these complex genomic rearrangements evolve during tumorigenesis is limited. In this study, we developed a highly multiplexed single-nucleus sequencing method to investigate copy number evolution in patients with triple-negative breast cancer. We sequenced 1,000 single cells from tumors in 12 patients and identified 1-3 major clonal subpopulations in each tumor that shared a common evolutionary lineage. For each tumor, we also identified a minor subpopulation of non-clonal cells that were classified as metastable, pseudodiploid or chromazemic. Phylogenetic analysis and mathematical modeling suggest that these data are unlikely to be explained by the gradual accumulation of copy number events over time. In contrast, our data challenge the paradigm of gradual evolution, showing that the majority of copy number aberrations are acquired at the earliest stages of tumor evolution, in short punctuated bursts, followed by stable clonal expansions that form the tumor mass.

Evolutionary Action Score of TP53 Coding Variants Is Predictive of Platinum Response in Head and Neck Cancer Patients.

TP53 is the most frequently altered gene in head and neck squamous cell carcinoma (HNSCC), with mutations occurring in over two thirds of cases; however, the predictive response of these mutations to cisplatin-based therapy remains elusive. In the current study, we evaluate the ability of the Evolutionary Action score of TP53-coding variants (EAp53) to predict the impact of TP53 mutations on response to chemotherapy. The EAp53 approach clearly identifies a subset of high-risk TP53 mutations associated with decreased sensitivity to cisplatin both in vitro and in vivo in preclinical models of HNSCC. Furthermore, EAp53 can predict response to treatment and, more importantly, a survival benefit for a subset of head and neck cancer patients treated with platinum-based therapy. Prospective evaluation of this novel scoring system should enable more precise treatment selection for patients with HNSCC.

Wee-1 kinase inhibition overcomes cisplatin resistance associated with high-risk TP53 mutations in head and neck cancer through mitotic arrest followed by senescence.

Although cisplatin has played a role in "standard-of-care" multimodality therapy for patients with advanced squamous cell carcinoma of the head and neck (HNSCC), the rate of treatment failure remains particularly high for patients receiving cisplatin whose tumors have mutations in the TP53 gene. We found that cisplatin treatment of HNSCC cells with mutant TP53 leads to arrest of cells in the G2 phase of the cell cycle, leading us to hypothesize that the wee-1 kinase inhibitor MK-1775 would abrogate the cisplatin-induced G2 block and thereby sensitize isogenic HNSCC cells with mutant TP53 or lacking p53 expression to cisplatin. We tested this hypothesis using clonogenic survival assays, flow cytometry, and in vivo tumor growth delay experiments with an orthotopic nude mouse model of oral tongue cancer. We also used a novel TP53 mutation classification scheme to identify which TP53 mutations are associated with limited tumor responses to cisplatin treatment. Clonogenic survival analyses indicate that nanomolar concentration of MK-1775 sensitizes HNSCC cells with high-risk mutant p53 to cisplatin. Consistent with its ability to chemosensitize, MK-1775 abrogated the cisplatin-induced G2 block in p53-defective cells leading to mitotic arrest associated with a senescence-like phenotype. Furthermore, MK-1775 enhanced the efficacy of cisplatin in vivo in tumors harboring TP53 mutations. These results indicate that HNSCC cells expressing high-risk p53 mutations are significantly sensitized to cisplatin therapy by the selective wee-1 kinase inhibitor, supporting the clinical evaluation of MK-1775 in combination with cisplatin for the treatment of patients with TP53 mutant HNSCC.