Homeostasis limits keratinocyte evolution.

Recent studies have revealed that normal human tissues accumulate many somatic mutations. In particular, human skin is riddled with mutations, with multiple subclones of variable sizes. Driver mutations are frequent and tend to have larger subclone sizes, suggesting selection. To begin to understand the histories encoded by these complex somatic mutations, we incorporated genomes into a simple agent-based skin-cell model whose prime directive is homeostasis. In this model, stem-cell survival is random and dependent on proximity to the basement membrane. This simple homeostatic skin model recapitulates the observed log-linear distributions of somatic mutations, where most mutations are found in increasingly smaller subclones that are typically lost with time. Hence, neutral mutations are "passengers" whose fates depend on the random survival of their stem cells, where a rarer larger subclone reflects the survival and spread of mutations acquired earlier in life. The model can also maintain homeostasis and accumulate more frequent and larger driver subclones if these mutations (NOTCH1 and TP53) confer relatively higher persistence in normal skin or during tissue damage (sunlight). Therefore, a relatively simple model of epithelial turnover indicates how observed passenger and driver somatic mutations could accumulate without violating the prime directive of homeostasis in normal human tissues.

Second-Order Effects of Chemotherapy Pharmacodynamics and Pharmacokinetics on Tumor Regression and Cachexia.

Drug dose response curves are ubiquitous in cancer biology, but these curves are often used to measure differential response in first-order effects: the effectiveness of increasing the cumulative dose delivered. In contrast, second-order effects (the variance of drug dose) are often ignored. Knowledge of second-order effects may improve the design of chemotherapy scheduling protocols, leading to improvements in tumor response without changing the total dose delivered. By considering treatment schedules with identical cumulative dose delivered, we characterize differential treatment outcomes resulting from high variance schedules (e.g. high dose, low dose) and low variance schedules (constant dose). We extend a previous framework used to quantify second-order effects, known as antifragility theory, to investigate the role of drug pharmacokinetics. Using a simple one-compartment model, we find that high variance schedules are effective for a wide range of cumulative dose values. Next, using a mouse-parameterized two-compartment model of 5-fluorouracil, we show that schedule viability depends on initial tumor volume. Finally, we illustrate the trade-off between tumor response and lean mass preservation. Mathematical modeling indicates that high variance dose schedules provide a potential path forward in mitigating the risk of chemotherapy-associated cachexia by preserving lean mass without sacrificing tumor response.

The harsh microenvironment in early breast cancer selects for a Warburg phenotype.

The harsh microenvironment of ductal carcinoma in situ (DCIS) exerts strong evolutionary selection pressures on cancer cells. We hypothesize that the poor metabolic conditions near the ductal center foment the emergence of a Warburg Effect (WE) phenotype, wherein cells rapidly ferment glucose to lactic acid, even in normoxia. To test this hypothesis, we subjected low-glycolytic breast cancer cells to different microenvironmental selection pressures using combinations of hypoxia, acidosis, low glucose, and starvation for many months and isolated single clones for metabolic and transcriptomic profiling. The two harshest conditions selected for constitutively expressed WE phenotypes. RNA sequencing analysis of WE clones identified the transcription factor KLF4 as potential inducer of the WE phenotype. In stained DCIS samples, KLF4 expression was enriched in the area with the harshest microenvironmental conditions. We simulated in vivo DCIS phenotypic evolution using a mathematical model calibrated from the in vitro results. The WE phenotype emerged in the poor metabolic conditions near the necrotic core. We propose that harsh microenvironments within DCIS select for a WE phenotype through constitutive transcriptional reprogramming, thus conferring a survival advantage and facilitating further growth and invasion.

Gattaca: Base-Pair Resolution Mutation Tracking for Somatic Evolution Studies using Agent-based Models.

Research over the past two decades has made substantial inroads into our understanding of somatic mutations. Recently, these studies have focused on understanding their presence in homeostatic tissue. In parallel, agent-based mechanistic models have emerged as an important tool for understanding somatic mutation in tissue; yet no common methodology currently exists to provide base-pair resolution data for these models. Here, we present Gattaca as the first method for introducing and tracking somatic mutations at the base-pair resolution within agent-based models that typically lack nuclei. With nuclei that incorporate human reference genomes, mutational context, and sequence coverage/error information, Gattaca is able to realistically evolve sequence data, facilitating comparisons between in silico cell tissue modeling with experimental human somatic mutation data. This user-friendly method, incorporated into each in silico cell, allows us to fully capture somatic mutation spectra and evolution.

IsoMaTrix: a framework to visualize the isoclines of matrix games and quantify uncertainty in structured populations.

SUMMARY: Evolutionary game theory describes frequency-dependent selection for fixed, heritable strategies in a population of competing individuals using a payoff matrix. We present a software package to aid in the construction, analysis and visualization of three-strategy matrix games. The IsoMaTrix package computes the isoclines (lines of zero growth) of matrix games, and facilitates direct comparison of well-mixed dynamics to structured populations on a lattice grid. IsoMaTrix computes fixed points, phase flow, trajectories, (sub)velocities and uncertainty quantification for stochastic effects in spatial matrix games. We describe a result obtained via IsoMaTrix's spatial games functionality, which shows that the timing of competitive release in a cancer model (under continuous treatment) critically depends on the initial spatial configuration of the tumor. AVAILABILITY AND IMPLEMENTATION: The code is available at: https://github.com/mathonco/isomatrix. SUPPLEMENTARY INFORMATION: Supplementary data are available at Bioinformatics online.

Roadmap on plasticity and epigenetics in cancer.

The role of plasticity and epigenetics in shaping cancer evolution and response to therapy has taken center stage with recent technological advances including single cell sequencing. This roadmap article is focused on state-of-the-art mathematical and experimental approaches to interrogate plasticity in cancer, and addresses the following themes and questions: is there a formal overarching framework that encompasses both non-genetic plasticity and mutation-driven somatic evolution? How do we measure and model the role of the microenvironment in influencing/controlling non-genetic plasticity? How can we experimentally study non-genetic plasticity? Which mathematical techniques are required or best suited? What are the clinical and practical applications and implications of these concepts?

The role of memory in non-genetic inheritance and its impact on cancer treatment resistance.

Intra-tumour heterogeneity is a leading cause of treatment failure and disease progression in cancer. While genetic mutations have long been accepted as a primary mechanism of generating this heterogeneity, the role of phenotypic plasticity is becoming increasingly apparent as a driver of intra-tumour heterogeneity. Consequently, understanding the role of this plasticity in treatment resistance and failure is a key component of improving cancer therapy. We develop a mathematical model of stochastic phenotype switching that tracks the evolution of drug-sensitive and drug-tolerant subpopulations to clarify the role of phenotype switching on population growth rates and tumour persistence. By including cytotoxic therapy in the model, we show that, depending on the strategy of the drug-tolerant subpopulation, stochastic phenotype switching can lead to either transient or permanent drug resistance. We study the role of phenotypic heterogeneity in a drug-resistant, genetically homogeneous population of non-small cell lung cancer cells to derive a rational treatment schedule that drives population extinction and avoids competitive release of the drug-tolerant sub-population. This model-informed therapeutic schedule results in increased treatment efficacy when compared against periodic therapy, and, most importantly, sustained tumour decay without the development of resistance.

Inferring Tumor Proliferative Organization from Phylogenetic Tree Measures in a Computational Model.

We use a computational modeling approach to explore whether it is possible to infer a solid tumor's cellular proliferative hierarchy under the assumptions of the cancer stem cell hypothesis and neutral evolution. We work towards inferring the symmetric division probability for cancer stem cells, since this is believed to be a key driver of progression and therapeutic response. Motivated by the advent of multiregion sampling and resulting opportunities to infer tumor evolutionary history, we focus on a suite of statistical measures of the phylogenetic trees resulting from the tumor's evolution in different regions of parameter space and through time. We find strikingly different patterns in these measures for changing symmetric division probability which hinge on the inclusion of spatial constraints. These results give us a starting point to begin stratifying tumors by this biological parameter and also generate a number of actionable clinical and biological hypotheses regarding changes during therapy, and through tumor evolutionary time. [Cancer; evolution; phylogenetics.].

Cell signaling heterogeneity is modulated by both cell-intrinsic and -extrinsic mechanisms: An integrated approach to understanding targeted therapy.

During the last decade, our understanding of cancer cell signaling networks has significantly improved, leading to the development of various targeted therapies that have elicited profound but, unfortunately, short-lived responses. This is, in part, due to the fact that these targeted therapies ignore context and average out heterogeneity. Here, we present a mathematical framework that addresses the impact of signaling heterogeneity on targeted therapy outcomes. We employ a simplified oncogenic rat sarcoma (RAS)-driven mitogen-activated protein kinase (MAPK) and phosphoinositide 3-kinase-protein kinase B (PI3K-AKT) signaling pathway in lung cancer as an experimental model system and develop a network model of the pathway. We measure how inhibition of the pathway modulates protein phosphorylation as well as cell viability under different microenvironmental conditions. Training the model on this data using Monte Carlo simulation results in a suite of in silico cells whose relative protein activities and cell viability match experimental observation. The calibrated model predicts distributional responses to kinase inhibitors and suggests drug resistance mechanisms that can be exploited in drug combination strategies. The suggested combination strategies are validated using in vitro experimental data. The validated in silico cells are further interrogated through an unsupervised clustering analysis and then integrated into a mathematical model of tumor growth in a homogeneous and resource-limited microenvironment. We assess posttreatment heterogeneity and predict vast differences across treatments with similar efficacy, further emphasizing that heterogeneity should modulate treatment strategies. The signaling model is also integrated into a hybrid cellular automata (HCA) model of tumor growth in a spatially heterogeneous microenvironment. As a proof of concept, we simulate tumor responses to targeted therapies in a spatially segregated tissue structure containing tumor and stroma (derived from patient tissue) and predict complex cell signaling responses that suggest a novel combination treatment strategy.

Hybrid Automata Library: A flexible platform for hybrid modeling with real-time visualization.

The Hybrid Automata Library (HAL) is a Java Library developed for use in mathematical oncology modeling. It is made of simple, efficient, generic components that can be used to model complex spatial systems. HAL's components can broadly be classified into: on- and off-lattice agent containers, finite difference diffusion fields, a GUI building system, and additional tools and utilities for computation and data collection. These components are designed to operate independently and are standardized to make them easy to interface with one another. As a demonstration of how modeling can be simplified using our approach, we have included a complete example of a hybrid model (a spatial model with interacting agent-based and PDE components). HAL is a useful asset for researchers who wish to build performant 1D, 2D and 3D hybrid models in Java, while not starting entirely from scratch. It is available on GitHub at https://github.com/MathOnco/HAL under the MIT License. HAL requires the Java JDK version 1.8 or later to compile and run the source code.

From cells to tissue: How cell scale heterogeneity impacts glioblastoma growth and treatment response.

Glioblastomas are aggressive primary brain tumors known for their inter- and intratumor heterogeneity. This disease is uniformly fatal, with intratumor heterogeneity the major reason for treatment failure and recurrence. Just like the nature vs nurture debate, heterogeneity can arise from intrinsic or environmental influences. Whilst it is impossible to clinically separate observed behavior of cells from their environmental context, using a mathematical framework combined with multiscale data gives us insight into the relative roles of variation from different sources. To better understand the implications of intratumor heterogeneity on therapeutic outcomes, we created a hybrid agent-based mathematical model that captures both the overall tumor kinetics and the individual cellular behavior. We track single cells as agents, cell density on a coarser scale, and growth factor diffusion and dynamics on a finer scale over time and space. Our model parameters were fit utilizing serial MRI imaging and cell tracking data from ex vivo tissue slices acquired from a growth-factor driven glioblastoma murine model. When fitting our model to serial imaging only, there was a spectrum of equally-good parameter fits corresponding to a wide range of phenotypic behaviors. When fitting our model using imaging and cell scale data, we determined that environmental heterogeneity alone is insufficient to match the single cell data, and intrinsic heterogeneity is required to fully capture the migration behavior. The wide spectrum of in silico tumors also had a wide variety of responses to an application of an anti-proliferative treatment. Recurrent tumors were generally less proliferative than pre-treatment tumors as measured via the model simulations and validated from human GBM patient histology. Further, we found that all tumors continued to grow with an anti-migratory treatment alone, but the anti-proliferative/anti-migratory combination generally showed improvement over an anti-proliferative treatment alone. Together our results emphasize the need to better understand the underlying phenotypes and tumor heterogeneity present in a tumor when designing therapeutic regimens.

Clonal hematopoiesis of indeterminate potential and its impact on patient trajectories after stem cell transplantation.

Clonal hematopoiesis of indeterminate potential (CHIP) is a recently identified process where older patients accumulate distinct subclones defined by recurring somatic mutations in hematopoietic stem cells. CHIP's implications for stem cell transplantation have been harder to identify due to the high degree of mutational heterogeneity that is present within the genetically distinct subclones. In order to gain a better understanding of CHIP and the impact of clonal dynamics on transplantation outcomes, we created a mathematical model of clonal competition dynamics. Our analyses highlight the importance of understanding competition intensity between healthy and mutant clones. Importantly, we highlight the risk that CHIP poses in leading to dominance of precancerous mutant clones and the risk of donor derived leukemia. Furthermore, we estimate the degree of competition intensity and bone marrow niche decline in mice during aging by using our modeling framework. Together, our work highlights the importance of better characterizing the ecological and clonal composition in hematopoietic donor populations at the time of stem cell transplantation.

Histoecology: Applying Ecological Principles and Approaches to Describe and Predict Tumor Ecosystem Dynamics Across Space and Time.

Cancer cells exist within a complex spatially structured ecosystem composed of resources and different cell types. As the selective pressures imposed by this environment determine the fate of cancer cells, an improved understanding of how this ecosystem evolves will better elucidate how tumors grow and respond to therapy. State of the art imaging methods can now provide highly resolved descriptions of the microenvironment, yielding the data required for a thorough study of its role in tumor growth and treatment resistance. The field of landscape ecology has been studying such species-environment relationship for decades, and offers many tools and perspectives that cancer researchers could greatly benefit from. Here, we discuss one such tool, species distribution modeling (SDM), that has the potential to, among other things, identify critical environmental factors that drive tumor evolution and predict response to therapy. SDMs only scratch the surface of how ecological theory and methods can be applied to cancer, and we believe further integration will take cancer research in exciting new and productive directions. Significance: Here we describe how species distribution modeling can be used to quantitatively describe the complex relationship between tumor cells and their microenvironment. Such a description facilitates a deeper understanding of cancers eco-evolutionary dynamics, which in turn sheds light on the factors that drive tumor growth and response to treatment.

A Mathematical Dissection of the Adaptation of Cell Populations to Fluctuating Oxygen Levels.

The disordered network of blood vessels that arises from tumour angiogenesis results in variations in the delivery of oxygen into the tumour tissue. This brings about regions of chronic hypoxia (i.e. sustained low oxygen levels) and regions with alternating periods of low and relatively higher oxygen levels, and makes it necessary for cancer cells to adapt to fluctuating environmental conditions. We use a phenotype-structured model to dissect the evolutionary dynamics of cell populations exposed to fluctuating oxygen levels. In this model, the phenotypic state of every cell is described by a continuous variable that provides a simple representation of its metabolic phenotype, ranging from fully oxidative to fully glycolytic, and cells are grouped into two competing populations that undergo heritable, spontaneous phenotypic variations at different rates. Model simulations indicate that, depending on the rate at which oxygen is consumed by the cells, dynamic nonlinear interactions between cells and oxygen can stimulate chronic hypoxia and cycling hypoxia. Moreover, the model supports the idea that under chronic-hypoxic conditions lower rates of phenotypic variation lead to a competitive advantage, whereas higher rates of phenotypic variation can confer a competitive advantage under cycling-hypoxic conditions. In the latter case, the numerical results obtained show that bet-hedging evolutionary strategies, whereby cells switch between oxidative and glycolytic phenotypes, can spontaneously emerge. We explain how these results can shed light on the evolutionary process that may underpin the emergence of phenotypic heterogeneity in vascularised tumours.

Mix and Match: Phenotypic Coexistence as a Key Facilitator of Cancer Invasion.

Invasion of healthy tissue is a defining feature of malignant tumours. Traditionally, invasion is thought to be driven by cells that have acquired all the necessary traits to overcome the range of biological and physical defences employed by the body. However, in light of the ever-increasing evidence for geno- and phenotypic intra-tumour heterogeneity, an alternative hypothesis presents itself: could invasion be driven by a collection of cells with distinct traits that together facilitate the invasion process? In this paper, we use a mathematical model to assess the feasibility of this hypothesis in the context of acid-mediated invasion. We assume tumour expansion is obstructed by stroma which inhibits growth and extra-cellular matrix (ECM) which blocks cancer cell movement. Further, we assume that there are two types of cancer cells: (i) a glycolytic phenotype which produces acid that kills stromal cells and (ii) a matrix-degrading phenotype that locally remodels the ECM. We extend the Gatenby-Gawlinski reaction-diffusion model to derive a system of five coupled reaction-diffusion equations to describe the resulting invasion process. We characterise the spatially homogeneous steady states and carry out a simulation study in one spatial dimension to determine how the tumour develops as we vary the strength of competition between the two phenotypes. We find that overall tumour growth is most extensive when both cell types can stably coexist, since this allows the cells to locally mix and benefit most from the combination of traits. In contrast, when inter-species competition exceeds intra-species competition the populations spatially separate and invasion arrests either: (i) rapidly (matrix-degraders dominate) or (ii) slowly (acid-producers dominate). Overall, our work demonstrates that the spatial and ecological relationship between a heterogeneous population of tumour cells is a key factor in determining their ability to cooperate. Specifically, we predict that tumours in which different phenotypes coexist stably are more invasive than tumours in which phenotypes are spatially separated.

High School Internship Program in Integrated Mathematical Oncology (HIP IMO): Five-Year Experience at Moffitt Cancer Center.

Modern cancer research, and the wealth of data across multiple spatial and temporal scales, has created the need for researchers that are well versed in the life sciences (cancer biology, developmental biology, immunology), medical sciences (oncology) and natural sciences (mathematics, physics, engineering, computer sciences). College undergraduate education traditionally occurs in disciplinary silos, which creates a steep learning curve at the graduate and postdoctoral levels that increasingly bridge multiple disciplines. Numerous colleges have begun to embrace interdisciplinary curricula, but students who double major in mathematics (or other quantitative sciences) and biology (or medicine) remain scarce. We identified the need to educate junior and senior high school students about integrating mathematical and biological skills, through the lens of mathematical oncology, to better prepare students for future careers at the interdisciplinary interface. The High school Internship Program in Integrated Mathematical Oncology (HIP IMO) at Moffitt Cancer Center has so far trained 59 students between 2015 and 2019. We report here on the program structure, training deliverables, curriculum and outcomes. We hope to promote interdisciplinary educational activities early in a student's career.

Evolutionary dynamics of competing phenotype-structured populations in periodically fluctuating environments.

Living species, ranging from bacteria to animals, exist in environmental conditions that exhibit spatial and temporal heterogeneity which requires them to adapt. Risk-spreading through spontaneous phenotypic variations is a known concept in ecology, which is used to explain how species may survive when faced with the evolutionary risks associated with temporally varying environments. In order to support a deeper understanding of the adaptive role of spontaneous phenotypic variations in fluctuating environments, we consider a system of non-local partial differential equations modelling the evolutionary dynamics of two competing phenotype-structured populations in the presence of periodically oscillating nutrient levels. The two populations undergo heritable, spontaneous phenotypic variations at different rates. The phenotypic state of each individual is represented by a continuous variable, and the phenotypic landscape of the populations evolves in time due to variations in the nutrient level. Exploiting the analytical tractability of our model, we study the long-time behaviour of the solutions to obtain a detailed mathematical depiction of the evolutionary dynamics. The results suggest that when nutrient levels undergo small and slow oscillations, it is evolutionarily more convenient to rarely undergo spontaneous phenotypic variations. Conversely, under relatively large and fast periodic oscillations in the nutrient levels, which bring about alternating cycles of starvation and nutrient abundance, higher rates of spontaneous phenotypic variations confer a competitive advantage. We discuss the implications of our results in the context of cancer metabolism.

NeoPredPipe: high-throughput neoantigen prediction and recognition potential pipeline.

BACKGROUND: Next generation sequencing has yielded an unparalleled means of quickly determining the molecular make-up of patient tumors. In conjunction with emerging, effective immunotherapeutics for a number of cancers, this rapid data generation necessitates a paired high-throughput means of predicting and assessing neoantigens from tumor variants that may stimulate immune response. RESULTS: Here we offer NeoPredPipe (Neoantigen Prediction Pipeline) as a contiguous means of predicting putative neoantigens and their corresponding recognition potentials for both single and multi-region tumor samples. NeoPredPipe is able to quickly provide summary information for researchers, and clinicians alike, on predicted neoantigen burdens while providing high-level insights into tumor heterogeneity given somatic mutation calls and, optionally, patient HLA haplotypes. Given an example dataset we show how NeoPredPipe is able to rapidly provide insights into neoantigen heterogeneity, burden, and immune stimulation potential. CONCLUSIONS: Through the integration of widely adopted tools for neoantigen discovery NeoPredPipe offers a contiguous means of processing single and multi-region sequence data. NeoPredPipe is user-friendly and adaptable for high-throughput performance. NeoPredPipe is freely available at https://github.com/MathOnco/NeoPredPipe .

EvoFreq: visualization of the Evolutionary Frequencies of sequence and model data.

BACKGROUND: High throughput sequence data has provided in depth means of molecular characterization of populations. When recorded at numerous time steps, such data can reveal the evolutionary dynamics of the population under study by tracking the changes in genotype frequencies over time. This necessitates a simple and flexible means of visualizing an increasingly complex set of data. RESULTS: Here we offer EvoFreq as a comprehensive tool set to visualize the evolutionary and population frequency dynamics of clones at a single point in time or as population frequencies over time using a variety of informative methods. EvoFreq expands substantially on previous means of visualizing the clonal, temporal dynamics and offers users a range of options for displaying their sequence or model data. CONCLUSIONS: EvoFreq, implemented in R with robust user options and few dependencies, offers a high-throughput means of quickly building, and interrogating the temporal dynamics of hereditary information across many systems. EvoFreq is freely available via https://github.com/MathOnco/EvoFreq.

The 2019 mathematical oncology roadmap.

Whether the nom de guerre is Mathematical Oncology, Computational or Systems Biology, Theoretical Biology, Evolutionary Oncology, Bioinformatics, or simply Basic Science, there is no denying that mathematics continues to play an increasingly prominent role in cancer research. Mathematical Oncology-defined here simply as the use of mathematics in cancer research-complements and overlaps with a number of other fields that rely on mathematics as a core methodology. As a result, Mathematical Oncology has a broad scope, ranging from theoretical studies to clinical trials designed with mathematical models. This Roadmap differentiates Mathematical Oncology from related fields and demonstrates specific areas of focus within this unique field of research. The dominant theme of this Roadmap is the personalization of medicine through mathematics, modelling, and simulation. This is achieved through the use of patient-specific clinical data to: develop individualized screening strategies to detect cancer earlier; make predictions of response to therapy; design adaptive, patient-specific treatment plans to overcome therapy resistance; and establish domain-specific standards to share model predictions and to make models and simulations reproducible. The cover art for this Roadmap was chosen as an apt metaphor for the beautiful, strange, and evolving relationship between mathematics and cancer.