Improving Cancer Treatment via Mathematical Modeling: Surmounting the Challenges Is Worth the Effort.

Drug delivery schedules are key factors in the efficacy of cancer therapies, and mathematical modeling of population dynamics and treatment responses can be applied to identify better drug administration regimes as well as provide mechanistic insights. To capitalize on the promise of this approach, the cancer field must meet the challenges of moving this type of work into clinics.

Mathematical modeling of PDGF-driven glioblastoma reveals optimized radiation dosing schedules.

Glioblastomas (GBMs) are the most common and malignant primary brain tumors and are aggressively treated with surgery, chemotherapy, and radiotherapy. Despite this treatment, recurrence is inevitable and survival has improved minimally over the last 50 years. Recent studies have suggested that GBMs exhibit both heterogeneity and instability of differentiation states and varying sensitivities of these states to radiation. Here, we employed an iterative combined theoretical and experimental strategy that takes into account tumor cellular heterogeneity and dynamically acquired radioresistance to predict the effectiveness of different radiation schedules. Using this model, we identified two delivery schedules predicted to significantly improve efficacy by taking advantage of the dynamic instability of radioresistance. These schedules led to superior survival in mice. Our interdisciplinary approach may also be applicable to other human cancer types treated with radiotherapy and, hence, may lay the foundation for significantly increasing the effectiveness of a mainstay of oncologic therapy. PAPERCLIP:

Dissecting dual roles of MyoD during lineage conversion to mature myocytes and myogenic stem cells.

The generation of myotubes from fibroblasts upon forced MyoD expression is a classic example of transcription factor-induced reprogramming. We recently discovered that additional modulation of signaling pathways with small molecules facilitates reprogramming to more primitive induced myogenic progenitor cells (iMPCs). Here, we dissected the transcriptional and epigenetic dynamics of mouse fibroblasts undergoing reprogramming to either myotubes or iMPCs using a MyoD-inducible transgenic model. Induction of MyoD in fibroblasts combined with small molecules generated Pax7+ iMPCs with high similarity to primary muscle stem cells. Analysis of intermediate stages of iMPC induction revealed that extinction of the fibroblast program preceded induction of the stem cell program. Moreover, key stem cell genes gained chromatin accessibility prior to their transcriptional activation, and these regions exhibited a marked loss of DNA methylation dependent on the Tet enzymes. In contrast, myotube generation was associated with few methylation changes, incomplete and unstable reprogramming, and an insensitivity to Tet depletion. Finally, we showed that MyoD's ability to bind to unique bHLH targets was crucial for generating iMPCs but dispensable for generating myotubes. Collectively, our analyses elucidate the role of MyoD in myogenic reprogramming and derive general principles by which transcription factors and signaling pathways cooperate to rewire cell identity.

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.

Computational modeling of pancreatic cancer reveals kinetics of metastasis suggesting optimum treatment strategies.

Pancreatic cancer is a leading cause of cancer-related death, largely due to metastatic dissemination. We investigated pancreatic cancer progression by utilizing a mathematical framework of metastasis formation together with comprehensive data of 228 patients, 101 of whom had autopsies. We found that pancreatic cancer growth is initially exponential. After estimating the rates of pancreatic cancer growth and dissemination, we determined that patients likely harbor metastases at diagnosis and predicted the number and size distribution of metastases as well as patient survival. These findings were validated in an independent database. Finally, we analyzed the effects of different treatment modalities, finding that therapies that efficiently reduce the growth rate of cells earlier in the course of treatment appear to be superior to upfront tumor resection. These predictions can be validated in the clinic. Our interdisciplinary approach provides insights into the dynamics of pancreatic cancer metastasis and identifies optimum therapeutic interventions.

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.

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.

Peripheral blood TCR clonotype diversity as an age-associated marker of breast cancer progression.

Immune escape is a prerequisite for tumor growth. We previously described a decline in intratumor activated cytotoxic T cells and T cell receptor (TCR) clonotype diversity in invasive breast carcinomas compared to ductal carcinoma in situ (DCIS), implying a central role of decreasing T cell responses in tumor progression. To determine potential associations between peripheral immunity and breast tumor progression, here, we assessed the peripheral blood TCR clonotype of 485 breast cancer patients diagnosed with either DCIS or de novo stage IV disease at younger (<45) or older (≥45) age. TCR clonotype diversity was significantly lower in older compared to younger breast cancer patients regardless of tumor stage at diagnosis. In the younger age group, TCR-α clonotype diversity was lower in patients diagnosed with de novo stage IV breast cancer compared to those diagnosed with DCIS. In the older age group, DCIS patients with higher TCR-α clonotype diversity were more likely to have a recurrence compared to those with lower diversity. Whole blood transcriptome profiles were distinct depending on the TCR-α Chao1 diversity score. There were more CD8+ T cells and a more active immune environment in DCIS tumors of young patients with higher peripheral blood TCR-α Chao1 diversity than in those with lower diversity. These results provide insights into the role that host immunity plays in breast cancer development across different age groups.

Increased HOXA5 expression provides a selective advantage for gain of whole chromosome 7 in IDH wild-type glioblastoma.

Glioblastoma is the most frequently occurring and invariably fatal primary brain tumor in adults. The vast majority of glioblastomas is characterized by chromosomal copy number alterations, including gain of whole chromosome 7 and loss of whole chromosome 10. Gain of whole chromosome 7 is an early event in gliomagenesis that occurs in proneural-like precursor cells, which give rise to all isocitrate dehydrogenase (IDH) wild-type glioblastoma transcriptional subtypes. Platelet-derived growth factor A (PDGFA) is one gene on chromosome 7 known to drive gliomagenesis, but, given its location near the end of 7p, there are likely several other genes located along chromosome 7 that select for its increased whole-chromosome copy number within glioblastoma cells. To identify other potential genes that could select for gain of whole chromosome 7, we developed an unbiased bioinformatics approach that identified homeobox A5 (HOXA5) as a gene whose expression correlated with gain of chromosome 7 and a more aggressive phenotype of the resulting glioma. High expression of HOXA5 in glioblastoma was associated with a proneural gene expression pattern and decreased overall survival in both human proneural and PDGF-driven mouse glioblastoma. Furthermore, HOXA5 overexpression promoted cellular proliferation and potentiated radioresistance. We also found enrichment of HOXA5 expression in recurrent human and mouse glioblastoma at first recurrence after radiotherapy. Overall, this study implicates HOXA5 as a chromosome 7-associated gene-level locus that promotes selection for gain of whole chromosome 7 and an aggressive phenotype in glioblastoma.

Mechanisms of response and resistance to combined decitabine and ipilimumab for advanced myeloid disease.

The challenge of eradicating leukemia in patients with acute myelogenous leukemia (AML) after initial cytoreduction has motivated modern efforts to combine synergistic active modalities including immunotherapy. Recently, the ETCTN/CTEP 10026 study tested the combination of the DNA methyltransferase inhibitor decitabine together with the immune checkpoint inhibitor ipilimumab for AML/myelodysplastic syndrome (MDS) either after allogeneic hematopoietic stem cell transplantation (HSCT) or in the HSCT-naïve setting. Integrative transcriptome-based analysis of 304 961 individual marrow-infiltrating cells for 18 of 48 subjects treated on study revealed the strong association of response with a high baseline ratio of T to AML cells. Clinical responses were predominantly driven by decitabine-induced cytoreduction. Evidence of immune activation was only apparent after ipilimumab exposure, which altered CD4+ T-cell gene expression, in line with ongoing T-cell differentiation and increased frequency of marrow-infiltrating regulatory T cells. For post-HSCT samples, relapse could be attributed to insufficient clearing of malignant clones in progenitor cell populations. In contrast to AML/MDS bone marrow, the transcriptomes of leukemia cutis samples from patients with durable remission after ipilimumab monotherapy showed evidence of increased infiltration with antigen-experienced resident memory T cells and higher expression of CTLA-4 and FOXP3. Altogether, activity of combined decitabine and ipilimumab is impacted by cellular expression states within the microenvironmental niche of leukemic cells. The inadequate elimination of leukemic progenitors mandates urgent development of novel approaches for targeting these cell populations to generate long-lasting responses. This trial was registered at www.clinicaltrials.gov as #NCT02890329.

Progression signature underlies clonal evolution and dissemination of multiple myeloma.

Clonal evolution drives tumor progression, dissemination, and relapse in multiple myeloma (MM), with most patients dying of relapsed disease. This multistage process requires tumor cells to enter the circulation, extravasate, and colonize distant bone marrow (BM) sites. Here, we developed a fluorescent or DNA-barcode clone-tracking system on MM PrEDiCT (progression through evolution and dissemination of clonal tumor cells) xenograft mouse model to study clonal behavior within the BM microenvironment. We showed that only the few clones that successfully adapt to the BM microenvironment can enter the circulation and colonize distant BM sites. RNA sequencing of primary and distant-site MM tumor cells revealed a progression signature sequentially activated along human MM progression and significantly associated with overall survival when evaluated against patient data sets. A total of 28 genes were then computationally predicted to be master regulators (MRs) of MM progression. HMGA1 and PA2G4 were validated in vivo using CRISPR-Cas9 in the PrEDiCT model and were shown to be significantly depleted in distant BM sites, indicating their role in MM progression and dissemination. Loss of HMGA1 and PA2G4 also compromised the proliferation, migration, and adhesion abilities of MM cells in vitro. Overall, our model successfully recapitulates key characteristics of human MM disease progression and identified potential new therapeutic targets for MM.

Designing optimal allocations for cancer screening using queuing network models.

Cancer is one of the leading causes of death, but mortality can be reduced by detecting tumors earlier so that treatment is initiated at a less aggressive stage. The tradeoff between costs associated with screening and its benefit makes the decision of whom to screen and when a challenge. To enable comparisons across screening strategies for any cancer type, we demonstrate a mathematical modeling platform based on the theory of queuing networks designed for quantifying the benefits of screening strategies. Our methodology can be used to design optimal screening protocols and to estimate their benefits for specific patient populations. Our method is amenable to exact analysis, thus circumventing the need for simulations, and is capable of exactly quantifying outcomes given variability in the age of diagnosis, rate of progression, and screening sensitivity and intervention outcomes. We demonstrate the power of this methodology by applying it to data from the Surveillance, Epidemiology and End Results (SEER) program. Our approach estimates the benefits that various novel screening programs would confer to different patient populations, thus enabling us to formulate an optimal screening allocation and quantify its potential effects for any cancer type and intervention.

Epigenetic restriction of extraembryonic lineages mirrors the somatic transition to cancer.

In mammals, the canonical somatic DNA methylation landscape is established upon specification of the embryo proper and subsequently disrupted within many cancer types. However, the underlying mechanisms that direct this genome-scale transformation remain elusive, with no clear model for its systematic acquisition or potential developmental utility. Here, we analysed global remethylation from the mouse preimplantation embryo into the early epiblast and extraembryonic ectoderm. We show that these two states acquire highly divergent genomic distributions with substantial disruption of bimodal, CpG density-dependent methylation in the placental progenitor. The extraembryonic epigenome includes specific de novo methylation at hundreds of embryonically protected CpG island promoters, particularly those that are associated with key developmental regulators and are orthologously methylated across most human cancer types. Our data suggest that the evolutionary innovation of extraembryonic tissues may have required co-option of DNA methylation-based suppression as an alternative to regulation by Polycomb-group proteins, which coordinate embryonic germ-layer formation in response to extraembryonic cues. Moreover, we establish that this decision is made deterministically, downstream of promiscuously used-and frequently oncogenic-signalling pathways, via a novel combination of epigenetic cofactors. Methylation of developmental gene promoters during tumorigenesis may therefore reflect the misappropriation of an innate trajectory and the spontaneous reacquisition of a latent, developmentally encoded epigenetic landscape.

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.

Polyploidy can drive rapid adaptation in yeast.

Polyploidy is observed across the tree of life, yet its influence on evolution remains incompletely understood. Polyploidy, usually whole-genome duplication, is proposed to alter the rate of evolutionary adaptation. This could occur through complex effects on the frequency or fitness of beneficial mutations. For example, in diverse cell types and organisms, immediately after a whole-genome duplication, newly formed polyploids missegregate chromosomes and undergo genetic instability. The instability following whole-genome duplications is thought to provide adaptive mutations in microorganisms and can promote tumorigenesis in mammalian cells. Polyploidy may also affect adaptation independently of beneficial mutations through ploidy-specific changes in cell physiology. Here we perform in vitro evolution experiments to test directly whether polyploidy can accelerate evolutionary adaptation. Compared with haploids and diploids, tetraploids undergo significantly faster adaptation. Mathematical modelling suggests that rapid adaptation of tetraploids is driven by higher rates of beneficial mutations with stronger fitness effects, which is supported by whole-genome sequencing and phenotypic analyses of evolved clones. Chromosome aneuploidy, concerted chromosome loss, and point mutations all provide large fitness gains. We identify several mutations whose beneficial effects are manifest specifically in the tetraploid strains. Together, these results provide direct quantitative evidence that in some environments polyploidy can accelerate evolutionary adaptation.

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.

Non-cell-autonomous driving of tumour growth supports sub-clonal heterogeneity.

Cancers arise through a process of somatic evolution that can result in substantial sub-clonal heterogeneity within tumours. The mechanisms responsible for the coexistence of distinct sub-clones and the biological consequences of this coexistence remain poorly understood. Here we used a mouse xenograft model to investigate the impact of sub-clonal heterogeneity on tumour phenotypes and the competitive expansion of individual clones. We found that tumour growth can be driven by a minor cell subpopulation, which enhances the proliferation of all cells within a tumour by overcoming environmental constraints and yet can be outcompeted by faster proliferating competitors, resulting in tumour collapse. We developed a mathematical modelling framework to identify the rules underlying the generation of intra-tumour clonal heterogeneity. We found that non-cell-autonomous driving of tumour growth, together with clonal interference, stabilizes sub-clonal heterogeneity, thereby enabling inter-clonal interactions that can lead to new phenotypic traits.

Clonal evolution in breast cancer revealed by single nucleus genome sequencing.

Sequencing studies of breast tumour cohorts have identified many prevalent mutations, but provide limited insight into the genomic diversity within tumours. Here we developed a whole-genome and exome single cell sequencing approach called nuc-seq that uses G2/M nuclei to achieve 91% mean coverage breadth. We applied this method to sequence single normal and tumour nuclei from an oestrogen-receptor-positive (ER(+)) breast cancer and a triple-negative ductal carcinoma. In parallel, we performed single nuclei copy number profiling. Our data show that aneuploid rearrangements occurred early in tumour evolution and remained highly stable as the tumour masses clonally expanded. In contrast, point mutations evolved gradually, generating extensive clonal diversity. Using targeted single-molecule sequencing, many of the diverse mutations were shown to occur at low frequencies (<10%) in the tumour mass. Using mathematical modelling we found that the triple-negative tumour cells had an increased mutation rate (13.3×), whereas the ER(+) tumour cells did not. These findings have important implications for the diagnosis, therapeutic treatment and evolution of chemoresistance in breast cancer.

Mathematical modeling identifies optimum lapatinib dosing schedules for the treatment of glioblastoma patients.

Human primary glioblastomas (GBM) often harbor mutations within the epidermal growth factor receptor (EGFR). Treatment of EGFR-mutant GBM cell lines with the EGFR/HER2 tyrosine kinase inhibitor lapatinib can effectively induce cell death in these models. However, EGFR inhibitors have shown little efficacy in the clinic, partly because of inappropriate dosing. Here, we developed a computational approach to model the in vitro cellular dynamics of the EGFR-mutant cell line SF268 in response to different lapatinib concentrations and dosing schedules. We then used this approach to identify an effective treatment strategy within the clinical toxicity limits of lapatinib, and developed a partial differential equation modeling approach to study the in vivo GBM treatment response by taking into account the heterogeneous and diffusive nature of the disease. Despite the inability of lapatinib to induce tumor regressions with a continuous daily schedule, our modeling approach consistently predicts that continuous dosing remains the best clinically feasible strategy for slowing down tumor growth and lowering overall tumor burden, compared to pulsatile schedules currently known to be tolerated, even when considering drug resistance, reduced lapatinib tumor concentrations due to the blood brain barrier, and the phenotypic switch from proliferative to migratory cell phenotypes that occurs in hypoxic microenvironments. Our mathematical modeling and statistical analysis platform provides a rational method for comparing treatment schedules in search for optimal dosing strategies for glioblastoma and other cancer types.

Bisulfite-independent analysis of CpG island methylation enables genome-scale stratification of single cells.

Conventional DNA bisulfite sequencing has been extended to single cell level, but the coverage consistency is insufficient for parallel comparison. Here we report a novel method for genome-wide CpG island (CGI) methylation sequencing for single cells (scCGI-seq), combining methylation-sensitive restriction enzyme digestion and multiple displacement amplification for selective detection of methylated CGIs. We applied this method to analyzing single cells from two types of hematopoietic cells, K562 and GM12878 and small populations of fibroblasts and induced pluripotent stem cells. The method detected 21 798 CGIs (76% of all CGIs) per cell, and the number of CGIs consistently detected from all 16 profiled single cells was 20 864 (72.7%), with 12 961 promoters covered. This coverage represents a substantial improvement over results obtained using single cell reduced representation bisulfite sequencing, with a 66-fold increase in the fraction of consistently profiled CGIs across individual cells. Single cells of the same type were more similar to each other than to other types, but also displayed epigenetic heterogeneity. The method was further validated by comparing the CpG methylation pattern, methylation profile of CGIs/promoters and repeat regions and 41 classes of known regulatory markers to the ENCODE data. Although not every minor methylation differences between cells are detectable, scCGI-seq provides a solid tool for unsupervised stratification of a heterogeneous cell population.