State-dependent evolutionary models reveal modes of solid tumour growth

Cell Tree Rings: the structure of somatic evolution as a human aging timer

Biological age is typically estimated using biomarkers whose states have been observed to correlate with chronological age. A persistent limitation of such aging clocks is that it is difficult to establish how the biomarker states are related to the mechanisms of aging. Somatic mutations could potentially form the basis for a more fundamental aging clock since the mutations are both markers and drivers of aging and have a natural timescale. Cell lineage trees inferred from these mutations reflect the somatic evolutionary process, and thus, it has been conjectured, the aging status of the body. Such a timer has been impractical thus far, however, because detection of somatic variants in single cells presents a significant technological challenge. Here, we show that somatic mutations detected using single-cell RNA sequencing (scRNA-seq) from thousands of cells can be used to construct a cell lineage tree whose structure correlates with chronological age. De novo single-nucleotide variants (SNVs) are detected in human peripheral blood mononuclear cells using a modified protocol. A default model based on penalized multiple regression of chronological age on 31 metrics characterizing the phylogenetic tree gives a Pearson correlation of 0.81 and a median absolute error of ~4 years between predicted and chronological ages. Testing of the model on a public scRNA-seq dataset yields a Pearson correlation of 0.85. In addition, cell tree age predictions are found to be better predictors of certain clinical biomarkers than chronological age alone, for instance glucose, albumin levels, and leukocyte count. The geometry of the cell lineage tree records the structure of somatic evolution in the individual and represents a new modality of aging timer. In addition to providing a numerical estimate of "cell tree age," it unveils a temporal history of the aging process, revealing how clonal structure evolves over life span. Cell Tree Rings complements existing aging clocks and may help reduce the current uncertainty in the assessment of geroprotective trials.

A theory of evolutionary dynamics on any complex population structure reveals stem cell niche architecture as a spatial suppressor of selection

How the spatial arrangement of a population shapes its evolutionary dynamics has been of long-standing interest in population genetics. Most previous studies assume a small number of demes or symmetrical structures that, most often, act as well-mixed populations. Other studies use network theory to study more heterogeneous spatial structures, however they usually assume small, regular networks, or strong constraints on the strength of selection considered. Here we build network generation algorithms, conduct evolutionary simulations and derive general analytic approximations for probabilities of fixation in populations with complex spatial structure. We build a unifying evolutionary theory across network families and derive the relevant selective parameter, which is a combination of network statistics, predictive of evolutionary dynamics. We also illustrate how to link this theory with novel datasets of spatial organization and use recent imaging data to build the cellular spatial networks of the stem cell niches of the bone marrow. Across a wide variety of parameters, we find these networks to be strong suppressors of selection, delaying mutation accumulation in this tissue. We also find that decreases in stem cell population size also decrease the suppression strength of the tissue spatial structure.

The enhanced energy metabolism in the tumor margin mediated by RRAD promotes the progression of oral squamous cell carcinoma

The tumor margin as the invasive front has been proven to be closely related to the progression and metastasis of oral squamous cell carcinoma (OSCC). However, how tumor cells in the marginal region obtain the extra energy needed for tumor progression is still unknown. Here, we used spatial metabolomics and the spatial transcriptome to identify enhanced energy metabolism in the tumor margin of OSCC and identified that the downregulation of Ras-related glycolysis inhibitor and calcium channel regulator (RRAD) in tumor cells mediated this process. The absence of RRAD enhanced the ingestion of glucose and malignant behaviors of tumor cells both in vivo and in vitro. Mechanically, the downregulation of RRAD promoted the internal flow of Ca2+ and elevated its concentration in the nucleus, which resulted in the activation of the CAMKIV-CREB1 axis to induce the transcription of the glucose transporter GLUT3. GLUT inhibitor-1, as an inhibitor of GLUT3, could suppress this vigorous energy metabolism and malignant behaviors caused by the downregulation of RRAD. Taken together, our study revealed that enhanced energy metabolism in the tumor margin mediated by RRAD promotes the progression of OSCC and proved that GLUT3 is a potential target for future treatment of OSCC.

The molecular evolution of cancer associated genes in mammals

Cancer is a disease that many multicellular organisms have faced for millions of years, and species have evolved various tumour suppression mechanisms to control oncogenesis. Although cancer occurs across the tree of life, cancer related mortality risks vary across mammalian orders, with Carnivorans particularly affected. Evolutionary theory predicts different selection pressures on genes associated with cancer progression and suppression, including oncogenes, tumour suppressor genes and immune genes. Therefore, we investigated the evolutionary history of cancer associated gene sequences across 384 mammalian taxa, to detect signatures of selection across categories of oncogenes (GRB2, FGL2 and CDC42), tumour suppressors (LITAF, Casp8 and BRCA2) and immune genes (IL2, CD274 and B2M). This approach allowed us to conduct a fine scale analysis of gene wide and site-specific signatures of selection across mammalian lineages under the lens of cancer susceptibility. Phylogenetic analyses revealed that for most species the evolution of cancer associated genes follows the species' evolution. The gene wide selection analyses revealed oncogenes being the most conserved, tumour suppressor and immune genes having similar amounts of episodic diversifying selection. Despite BRCA2's status as a key caretaker gene, episodic diversifying selection was detected across mammals. The site-specific selection analyses revealed that the two apoptosis associated domains of the Casp8 gene of bats (Chiroptera) are under opposing forces of selection (positive and negative respectively), highlighting the importance of site-specific selection analyses to understand the evolution of highly complex gene families. Our results highlighted the need to critically assess different types of selection pressure on cancer associated genes when investigating evolutionary adaptations to cancer across the tree of life. This study provides an extensive assessment of cancer associated genes in mammals with highly representative, and substantially large sample size for a comparative genomic analysis in the field and identifies various avenues for future research into the mechanisms of cancer resistance and susceptibility in mammals.

Spatial interactions modulate tumor growth and immune infiltration

Lenia, a cellular automata framework used in artificial life, provides a natural setting to implement mathematical models of cancer incorporating features such as morphogenesis, homeostasis, motility, reproduction, growth, stimuli response, evolvability, and adaptation. Historically, agent-based models of cancer progression have been constructed with rules that govern birth, death and migration, with attempts to map local rules to emergent global growth dynamics. In contrast, Lenia provides a flexible framework for considering a spectrum of local (cell-scale) to global (tumor-scale) dynamics by defining an interaction kernel governing density-dependent growth dynamics. Lenia can recapitulate a range of cancer model classifications including local or global, deterministic or stochastic, non-spatial or spatial, single or multi-population, and off or on-lattice. Lenia is subsequently used to develop data-informed models of 1) single-population growth dynamics, 2) multi-population cell-cell competition models, and 3) cell migration or chemotaxis. Mathematical modeling provides important mechanistic insights. First, short-range interaction kernels provide a mechanism for tumor cell survival under conditions with strong Allee effects. Next, we find that asymmetric interaction tumor-immune kernels lead to poor immune response. Finally, modeling recapitulates immune-ECM interactions where patterns of collagen formation provide immune protection, indicated by an emergent inverse relationship between disease stage and immune coverage.

A seven-step guide to spatial, agent-based modelling of tumour evolution

Spatial agent-based models are frequently used to investigate the evolution of solid tumours subject to localized cell-cell interactions and microenvironmental heterogeneity. As spatial genomic, transcriptomic and proteomic technologies gain traction, spatial computational models are predicted to become ever more necessary for making sense of complex clinical and experimental data sets, for predicting clinical outcomes, and for optimizing treatment strategies. Here we present a non-technical step by step guide to developing such a model from first principles. Stressing the importance of tailoring the model structure to that of the biological system, we describe methods of increasing complexity, from the basic Eden growth model up to off-lattice simulations with diffusible factors. We examine choices that unavoidably arise in model design, such as implementation, parameterization, visualization and reproducibility. Each topic is illustrated with examples drawn from recent research studies and state of the art modelling platforms. We emphasize the benefits of simpler models that aim to match the complexity of the phenomena of interest, rather than that of the entire biological system. Our guide is aimed at both aspiring modellers and other biologists and oncologists who wish to understand the assumptions and limitations of the models on which major cancer studies now so often depend.

Spatial interactions modulate tumor growth and immune infiltration

Direct observation of immune cell trafficking patterns and tumor-immune interactions is unlikely in human tumors with currently available technology, but computational simulations based on clinical data can provide insight to test hypotheses. It is hypothesized that patterns of collagen formation evolve as a mechanism of immune escape, but the exact nature of the interaction between immune cells and collagen is poorly understood. Spatial data quantifying the degree of collagen fiber alignment in squamous cell carcinomas indicates that late stage disease is associated with highly aligned fibers. Here, we introduce a computational modeling framework (called Lenia) to discriminate between two hypotheses: immune cell migration that moves 1) parallel or 2) perpendicular to collagen fiber orientation. The modeling recapitulates immune-ECM interactions where collagen patterns provide immune protection, leading to an emergent inverse relationship between disease stage and immune coverage. We also illustrate the capabilities of Lenia to model the evolution of tumor progression and immune predation. Lenia provides a flexible framework for considering a spectrum of local (cell-scale) to global (tumor-scale) dynamics by defining a kernel cell-cell interaction function that governs tumor growth dynamics under immune predation with immune cell migration. Mathematical modeling provides important mechanistic insights into cell interactions. Short-range interaction kernels provide a mechanism for tumor cell survival under conditions with strong Allee effects, while asymmetric tumor-immune interaction kernels lead to poor immune response. Thus, the length scale of tumor-immune interactions drives tumor growth and infiltration.

Modeling the Effect of Spatial Structure on Solid Tumor Evolution and Circulating Tumor DNA Composition

Circulating tumor DNA (ctDNA) monitoring, while sufficiently advanced to reflect tumor evolution in real time and inform cancer diagnosis, treatment, and prognosis, mainly relies on DNA that originates from cell death via apoptosis or necrosis. In solid tumors, chemotherapy and immune infiltration can induce spatially variable rates of cell death, with the potential to bias and distort the clonal composition of ctDNA. Using a stochastic evolutionary model of boundary-driven growth, we study how elevated cell death on the edge of a tumor can simultaneously impact driver mutation accumulation and the representation of tumor clones and mutation detectability in ctDNA. We describe conditions in which invasive clones are over-represented in ctDNA, clonal diversity can appear elevated in the blood, and spatial bias in shedding can inflate subclonal variant allele frequencies (VAFs). Additionally, we find that tumors that are mostly quiescent can display similar biases but are far less detectable, and the extent of perceptible spatial bias strongly depends on sequence detection limits. Overall, we show that spatially structured shedding might cause liquid biopsies to provide highly biased profiles of tumor state. While this may enable more sensitive detection of expanding clones, it could also increase the risk of targeting a subclonal variant for treatment. Our results indicate that the effects and clinical consequences of spatially variable cell death on ctDNA composition present an important area for future work.

PhyloJunction: a computational framework for simulating, developing, and teaching evolutionary models

We introduce PhyloJunction, a computational framework designed to facilitate the prototyping, testing, and characterization of evolutionary models. PhyloJunction is distributed as an open-source Python library that can be used to implement a variety of models, through its flexible graphical modeling architecture and dedicated model specification language. Model design and use are exposed to users via command-line and graphical interfaces, which integrate the steps of simulating, summarizing, and visualizing data. This paper describes the features of PhyloJunction - which include, but are not limited to, a general implementation of a popular family of phylogenetic diversification models - and, moving forward, how it may be expanded to not only include new models, but to also become a platform for conducting and teaching statistical learning.

A new universal system of tree shape indices

The comparison and categorization of tree diagrams is fundamental to large parts of biology, linguistics, computer science, and other fields, yet the indices currently applied to describing tree shape have important flaws that complicate their interpretation and limit their scope. Here we introduce a new system of indices with no such shortcomings. Our indices account for node sizes and branch lengths and are robust to small changes in either attribute. Unlike currently popular phylogenetic diversity, phylogenetic entropy, and tree balance indices, our definitions assign interpretable values to all rooted trees and enable meaningful comparison of any pair of trees. Our self-consistent definitions further unite measures of diversity, richness, balance, symmetry, effective height, effective outdegree, and effective branch count in a coherent system, and we derive numerous simple relationships between these indices. The main practical advantages of our indices are in 1) quantifying diversity in non-ultrametric trees; 2) assessing the balance of trees that have non-uniform branch lengths or node sizes; 3) comparing the balance of trees with different leaf counts or outdegrees; 4) obtaining a coherent, generic, multidimensional quantification of tree shape that is robust to sampling error and inferential error. We illustrate these features by comparing the shapes of trees representing the evolution of HIV and of Uralic languages, and trees generated by computational models of tumour evolution. Given the ubiquity of tree structures, we identify a wide range of applications across diverse domains.

Modeling the effect of spatial structure on solid tumor evolution and ctDNA composition

Circulating tumor DNA (ctDNA) monitoring, while sufficiently advanced to reflect tumor evolution in real time and inform on cancer diagnosis, treatment, and prognosis, mainly relies on DNA that originates from cell death via apoptosis or necrosis. In solid tumors, chemotherapy and immune infiltration can induce spatially variable rates of cell death, with the potential to bias and distort the clonal composition of ctDNA. Using a stochastic evolutionary model of boundary-driven growth, we study how elevated cell death on the edge of a tumor can simultaneously impact driver mutation accumulation and the representation of tumor clones and mutation detectability in ctDNA. We describe conditions in which invasive clones end up over-represented in ctDNA, clonal diversity can appear elevated in the blood, and spatial bias in shedding can inflate subclonal variant allele frequencies (VAFs). Additionally, we find that tumors that are mostly quiescent can display similar biases, but are far less detectable, and the extent of perceptible spatial bias strongly depends on sequence detection limits. Overall, we show that spatially structured shedding might cause liquid biopsies to provide highly biased profiles of tumor state. While this may enable more sensitive detection of expanding clones, it could also increase the risk of targeting a subclonal variant for treatment. Our results indicate that the effects and clinical consequences of spatially variable cell death on ctDNA composition present an important area for future work.

cloneRate: fast estimation of single-cell clonal dynamics using coalescent theory

MOTIVATION

While evolutionary approaches to medicine show promise, measuring evolution itself is difficult due to experimental constraints and the dynamic nature of body systems. In cancer evolution, continuous observation of clonal architecture is impossible, and longitudinal samples from multiple timepoints are rare. Increasingly available DNA sequencing datasets at single-cell resolution enable the reconstruction of past evolution using mutational history, allowing for a better understanding of dynamics prior to detectable disease. There is an unmet need for an accurate, fast, and easy-to-use method to quantify clone growth dynamics from these datasets.

RESULTS

We derived methods based on coalescent theory for estimating the net growth rate of clones using either reconstructed phylogenies or the number of shared mutations. We applied and validated our analytical methods for estimating the net growth rate of clones, eliminating the need for complex simulations used in previous methods. When applied to hematopoietic data, we show that our estimates may have broad applications to improve mechanistic understanding and prognostic ability. Compared to clones with a single or unknown driver mutation, clones with multiple drivers have significantly increased growth rates (median 0.94 versus 0.25 per year; P = 1.6×10-6). Further, stratifying patients with a myeloproliferative neoplasm (MPN) by the growth rate of their fittest clone shows that higher growth rates are associated with shorter time to MPN diagnosis (median 13.9 versus 26.4 months; P = 0.0026).

AVAILABILITY AND IMPLEMENTATION

We developed a publicly available R package, cloneRate, to implement our methods (Package website: https://bdj34.github.io/cloneRate/). Source code: https://github.com/bdj34/cloneRate/.