A variety of tumor cell populations, such as treatment sensitive and resistant cells or supporting stroma and immune cells, have been modeled mathematically. In this section, we explored how researchers have studied the interactions between different cell populations using ecological modeling approaches such as Lotka-Volterra competition models. The most common approach to model the temporal dynamics of tumor populations is using ordinary differential equations (ODEs) or implemented as an agent-based model (ABM). ODEs describe changes in variables over a continuous time. The input are often population sizes, concentrations, or other continuous variables. On the other hand, ABMs simulate the discrete dynamics of a population of agents with specific rules regarding the behavior of interactions of the agents with other agents and their environment. Based on the discrete nature, ABMs often capture stochastic interactions in heterogenous systems and have often been applied to studying dynamics where spatial patterning is of interest.

One notable application of deterministic population dynamics modeling that applies ecological constructs to cancer research is the ODE model developed by Grassberger et al. that predicted the growth of drug-sensitive cells, drug-tolerant or drug-persistent cells, and drug-resistant cells to tyrosine kinase inhibitors (TKI) in EGFR-mutant lung cancer. The ODEs for the TKI-persister and TKI-resistant populations included a term that described the transition between the two phenotypes as well as an interaction term that explained the interplay between TKI-persistent and TKI-resistant cells. By using continuous lung tumor volume data abstracted from CT scans, they demonstrated that the growth rate of TKI-resistant cells was correlated to tumor progression (Grassberger et al., 2019). Cerasuolo et al. studied the evolution of TKI resistance using a more complex system of five differential equations describing the growth of sensitive and resistant cells, the changes in pharmacokinetics of a prostate cancer drug called enzalutamide, and the dynamics of the tumor microenvironment. In this model, the authors parameterized the ODE model using in vitro and in vivo data. Cell counts were collected in vitro from murine prostates adenocarcinoma cell lines cultures using flow cytometry. Tumor volume data was extracted from measuring in vivo mouse tumors. For both sets of data, their growth was assumed to follow logistic growth. Cerasuolo et al. also introduced a diffusion term to account for the random variation in cell number in the tumor microenvironment, also known as Brownian motion (Cerasuolo et al., 2020). The addition of stochasticity provided a more accurate representation of the behavior of the microenvironment in prostate cancer allowing the researchers to explore a second drug used in combination with enzalutamide. ODEs have also been used to model androgen-dependent and androgen-independent prostate cancer cells. Draghi et al. developed a system of equations that described the interaction between the health status of the patient including their ability to sustain treatment and develop resistance, the prostate-specific antigen (PSA) level produced by the androgen-dependent cells and androgen-independent cells, and the androgen level (Draghi et al., 2019). PSA levels collected from patients under intermittent therapy over time were used to parameterize this ODE model. The researchers were able to show that mathematically modeling PSA levels over time has the potential to predict treatment response and optimize treatment schedules.

When it comes to studying the interactions of different cancer populations that co-exist in the same microenvironment and compete for survival, Lotka-Volterra competition models are frequently applied. These models originate from ecology through studying predation in different animal populations, where one species acts as the prey whose population size is dependent on the rate of predation by the predator. These concepts have also been applied to model different cancers using ODEs that represent populations of cells interacting and competing for the same resources. The equations are usually modeled as logistic growth where total population growth is ultimately dampened by the availability of resources and some carrying capacity of the environment. However, different variations of the Lotka-Volterra models have been explored in order to fit the growth patterns of cell populations (Michelson, 1987; Michelson et al., 1987). Evolutionary-based therapy scheduling, also referred to as adaptive therapy, uses Lotka-Volterra interactions between cancer cell populations of interest to inform chemotherapeutic schedules. These schedules aim to reduce tumor sizes without allowing dominance of a treatment-resistant cell type by considering the dynamics between different cell types within the tumor. Adaptive therapy aims to move away from administrating maximum tolerated doses, which often result in the development of treatment resistance; to administering strategically timed chemotherapy doses. This aims to better control tumor size at lower cumulative treatment doses and allows for treatment ‘holidays’ for the patients, which may also reduce side effects patients experience (Belkhir et al., 2021). Lower cumulative drug doses determined by mathematical models allow for enough sensitive cells to still outcompete resistant cells for dominance in the tumor population. By delaying the development of resistance, sensitive clones can still be treated resulting on longer life expectancy. Lotka-Volterra based models are usually formatted as differential equations of pairs of treatment-resistant and treatment-sensitive populations or quiescent and proliferating populations and a term that captures the effects of the drug on the sensitive or proliferating population (Strobl et al., 2023). Adaptive therapy has shown to delay the onset of resistance of abiraterone in a clinical trial of metastatic castration-resistant prostate cancer (NCT02415621) where treatment decisions were made based on a mathematical model of androgen-dependent cells, androgen-independent cells, and testosterone producing cells with patient PSA level serving as a biomarker for tumor burden (Cunningham et al., 2020). This trial resulted in an increased of time to progression of 19.2 months with a 46% decrease in total drug dose compared to standard of care therapies (Zhang et al., 2017). Ongoing adaptive therapy clinical trials leveraging similar approaches have continue in prostate cancer and been applied to ovarian cancer patients (NCT05080556) (Draghi et al., 2019; Michelson, 1987; Gallaher et al., 2023). Outside of adenocarcinomas, adaptive therapy clinical trials have been developed for rhabdomyosarcoma (NCT04388839) (Reed et al., 2020), BRAF-mutant melanoma (NCT03543969), and basal cell carcinoma (NCT05651828) indicating that this treatment strategy could be applied to variety of cancer types.

Several ODE models have been used to study molecular changes within cancer cells exposed to different treatments as well as their effect on the evolution of resistance and overall changes in tumor size. These changes will cause the cells to adapt to new environments that lack certain biomolecules or compete with other populations for resources. This has been explored in both cervical and non-small cell lung cancer, where tumor volume has been analyzed in relation to the radiation fraction administered. Radiation can be modeled based on the number of oxygenated cells and hypoxic areas in the tumor, which are included in the ODEs as oxygenation rate constants (Chen et al., 2019; Belfatto et al., 2016). ODE models have also studied cancer progression regarding metabolic processes, in the absence of treatment. Voutouri et al. studied how osmotic pressures and mechanical changes in the microenvironment can affect tumor progression. They developed a model that incorporates the fluid and solid phase of the tumor, but also the transport of cation and anions on glycosaminoglycan chains. In this case, a deformation gradient tensor is used to quantify the three-dimensional changes in the tumor, which provides information about the degree of tumor deformation due to different elements (Voutouri and Stylianopoulos, 2014). A different research group developed a nutrient-depletion model, where they were able to study the interaction of leukemia and cervical cancer cells using population ecology concepts. Co-cultures experiments were developed to investigate whether proliferating cervical cancer epithelial cells can inhibit leukemia T-cells. The populations were modeled following prey-predation like interactions assuming that each cell type lived in a niche forming two different cell ecosystems that interact by exchanging molecules (Sega and Chignola, 2014). Studies in pancreatic and breast cancer have focused on tumor progression as it pertains to different metabolic conditions in the tumor microenvironment. These groups have developed systems of differential equations that explain changes on the amount of glucose, oxygen, pH, and other metabolite concentrations used in processes such as glycolysis, tricarboxylic acid cycle, or the pentose pathway (Roy and Finley, 2017; Damaghi et al., 2021). These studies allowed researchers to create in silico simulations of enzyme or element deprivation and evaluate the effects on cell proliferation. Specifically, oxygen levels in the tumor have been widely studied since the lack of this element in the tumor microenvironment leads to cancer invasion and metastasis. Oxygen levels, however, are most commonly modeled using linear models of diffusion and not differential equations (Milotti et al., 2020). Information obtained from using both approaches suggest different metabolic targets to slow the progression of certain cancers. Once again, mathematical modeling was able to predict cellular growth based on different metabolic processes.

While traditional ODEs have often been leveraged in population dynamics modeling of tumor ecology, some researchers hoping to capture the range of behaviors a tumor might exhibit under different selective pressures have turned to ABMs. This framework consists of individual agents, which are unique and autonomous cells or molecules that can interact with each other and their microenvironment following a set of specified rules and these frameworks often incorporate stochasticity and can include spatial resolved information. For instance, Deisboeck et al. developed an in silico ABM to simulate the possible effects of engineered cells in the progression of the tumor. In their model, the engineered cells were primary tumor cells from a patient modified to express higher proliferation rates than native tumor cells. They also carry an on-off switch that can be targeted therapeutically to induce apoptosis. The researchers hypothesized that the engineered cells could outgrow the native tumor cells and the final population could be eliminated using a therapeutic. Here, each cell (engineered or native) is an agent that can obtain ‘proliferation’ credits, which determine whether the cell will proceed to proliferation or remain in place allowing the other population to dominate in the environment. The simulations capture a 3D tumor microenvironment composed of an initial native tumor cell and an initial engineered cell, which has a higher proliferation rate compared to native tumor cells (Deisboeck and Wang, 2008). The simulations showed the efficacy of engineered tumor cells at controlling and potentially combating cancer once they reach a certain proliferation level. However, these engineered cells would have to be developed from primary tumor cells and subsequently implanted in close proximities to the tumor resulting in a highly invasive procedure for the patient. Maley et al. developed an ABM for the pre-cancerous condition Barrett’s esophagus, which can develop into esophageal adenocarcinoma. The researchers described cancer progression in three stages (normal, pre-cancerous, and malignant) based on the accumulation of neutral and selective mutations. The transition of these agents or cells between cancer stages was determined by the number of accumulated mutations and their ability to increase the mutation rates, as well as the time required for cell division. This research was able to provide insights on the number of different mutation types necessary for the development of malignancy (Maley and Forrest, 2000).

Population based modeling has been widely applied to adenocarcinomas given that it forms solid tumors, and temporal data such as tumor size, or population concentrations, can be obtained from CT scans, flow cytometry counts, and in vitro cultures. In this section, we have explored the usage two dynamics models: continuous models (e.g., ODE models) and discrete models incorporating stochasticity (e.g., ABMs) to obtain information about cancer progression. Deterministic models can provide continuous estimates of cancer behavior over time based distinct timepoint data. These models often use ODEs to describe systems that have well-defined parameters and known interactions. ODE modeling is a good platform to test out different mechanisms of action. In contrast, stochastic models include randomness or probabilistic elements into the system and thus need to be simulated many times to get the range of behaviors of the system of interest. Thus, these models are best suited to answer questions about behaviors such as the unpredictable nature of mutations, cell interactions, or treatment responses at low population sizes. Deterministic models are commonly used to predict average or bulk trends experienced by the populations of interest making their parameterization a computationally efficient process. Stochastic modeling requires additional input data and since it simulates random fluctuations it can be more computationally expensive. However, the use of stochasticity is not always appropriate when modeling cancer progression. Whenever populations are small, random mutations can greatly affect the development of that population and accuracy of the model. Mutations in large populations are not as impactful and could be left out of the model. Both deterministic and stochastic modeling approaches contribute to the understanding of the complex interactions in cancer progression, so their usage depends on the biological question, and available computational resources and data. Some researchers have also leveraged a hybrid approach that changes based on population size as demonstrated in modeling of immunotherapy treatment response in lymphoma (Kimmel et al., 2021).

While the interactions between cell populations is important for cancer development, their spatial distribution can also play a critical role. Similar to animal populations distributed within an ecosystem, cellular populations in the tumor and its microenvironment are often arranged in specific patterns, which impacts the progression of the disease. What differentiates cancer from the populations commonly studied in ecology is the greater degree of spatial heterogeneity due to the presence of different cell types, phenotypes, and mutations. Hence, it is crucial to study the spatial distribution and clonal evolution patterns of these populations. Spatial data can be collected using different modalities including imaging, transcriptomic, or genomic techniques. Often, this data is collected from different areas in the tumor or from the primary tumor and distant mets and the expression profiles are compared. For instance, single-cell RNA sequencing has often been used to offer spatial context of cell locations in addition to transcriptomic data. Histological tumor samples have also been used to extrapolate cellular locations and predict treatment effectiveness based on cell distributions. Spatial ecology modeling ensures that therapies can effectively target the entire tumor and its microenvironment, rather than just a subset of populations (González-García et al., 2002; Yuan, 2016).

Spatial modeling approaches used in cancer not only deal with cell-to-cell interactions, but also with changes in cell shapes and sizes. Sapi et al. created a mathematical model to understand the changes in colorectal tumor volume under bevacizumab therapy. Their 2D model used a digital caliper to collect different measurements of the subcutaneous tumor under different bevacizumab concentrations without the need of performing an MRI. Measurements including tumor length and width were recorded for each mouse at least six times during the duration of the experiment. The authors showed that daily low doses of bevacizumab were able to cause a greater tumor volume reduction than a singular high dose (Sápi et al., 2015). Yamamoto et al. have also investigated effective treatment schedules and drug concentrations for pancreatic cancer, particularly with a combination of FOLFIRINOX and nab-paclitaxel therapy. They used time series tumor volume data from primary and metastatic site MRIs in a three cell-type logistic branching growth model. The growth was modeled starting with an initial cell that gives rise to proliferating cells that can accumulate mutations and transform to metastatic cells and establish metastatic sites. Yamamoto et al. introduced a quantifiable metric called Local Advancement Index (LAI), as an indicator of metastasis where the growth rate of the cells decreases as the tumor size increases. They concluded that tumors with lower LAI have a higher potential to metastasize and develop more metastases, (Yamamoto et al., 2019).

Volume changes have also been applied to ground-glass opacities, which are areas in the lung tissue that appear opaque. These opacities are indicative of various conditions including non-small cell lung cancer, especially when large volume alterations are observed. To understand the spatial-temporal dynamics of these opacities, Farkas et al. developed a model using diffusive logistic partial differential equations (PDE), which incorporated random movement of particles (Brownian motion). In this case, CT histograms, a graphical representation of the pixel intensities in CT region, were used as inputs for the model, which outputs a spatial density at a specific time. This information was used to identify doubling times of the ground glass opacities volume, serving as a marker of lung cancer progression (Farkas et al., 2018).

While the spatial location and volume of the tumor are important when selecting treatment strategies, the cell-to-cell interactions within the tumor can also play a role in determining cancer progression. Research groups such as Yang et al. have used both bulk and single-cell transcriptomic data to explore the presence and impact of tumor-associated macrophages on pancreatic adenocarcinoma development. They used CellPhoneDB (Efremova et al., 2020), a database analyzing cell-cell communications, to construct a network describing the interactions between macrophages and cancer cells. Additionally, the authors identified the presence of inflammatory macrophages in the tumor microenvironment as well as specific ligand-receptor pairs between tumor and immune cells (Yang et al., 2023). Information obtained from these networks provided a better understanding of macrophage interactions with pancreatic cancer providing an opportunity to develop targeted immunotherapies. Shan et al. also utilized spatially resolved bulk and single-cell transcriptomics data to detect the presence of stroma, cancer, normal pancreatic, and duct epithelium cells and their interactions in pancreatic cancer. The authors used SpaMOD (spatial molecular patterns) (Andres et al., 2001), a spatial modeling tool based on a partial least squares regression approach, to integrate spatially resolved transcriptomic data derived from the same pancreatic cancer tissue. By combining the two types of data, researchers identified the location of certain cell types within the tissue and their function. Understanding the spatial and functional organization of tissue has the potential to provide information about cancer stages and progression (Shan et al., 2022). Casesent et al. introduced a novel approach to collect spatially resolved genomic data called topographic single-cell sequencing, to study the invasiveness of ductal carcinoma in the breast. This innovative technique integrates laser-capture microdissection, laser catapulting, whole-genome amplification, and single-cell DNA sequencing to provide spatially resolved genomic data. After studying in situ and invasive ductal carcinoma samples, topographic single-cell sequencing identified a genomic transition from in situ to invasive phenotypes due to mutations. Thus the presence of mutant clones in the duct could indicate the potential development of metastasis (Casasent et al., 2018).

Understanding clonal evolution within a spatial context offers valuable insights into the selection and metastasis of distinct clones. The strength of selection can be measured using next-generation sequencing data (NGS), quantified as variant allele frequencies (VAFs). These allele frequencies represent the proportion of alleles at a particular locus that differ from the reference genome sequence. In cancer, the distribution of VAFs represent the genetic variants within the population. These distributions can be used to simulate spatial tumor growth models aimed at identifying mechanisms driving cancer progression, including the influence of natural selection and genetic drift. This methodology has been applied to multiple types of adenocarcinomas such as colorectal, lung, and esophageal (Sun et al., 2017). Utilizing variant information obtained from NGS, Sievers et al. identified that cancerous mutations can arise in small colon polyps. Modeling the molecular pathways involved in the growth of colorectal polyps and considering their mutation profiles, could offer insights in the clonal evolution of polyps that may progress into cancer (Sievers et al., 2017). In lung adenocarcinomas, studies aimed to uncover the clonal evolution of histomorphological patterns used global DNA methylation data, whole exome sequencing (WES), and RNA sequencing (Dietz et al., 2019; Karasaki et al., 2023). Dietz et al. focused on studying DNA methylation sites and the change in gene copy number in a phylogenetic framework to identify the origin of metastasis in the lymph nodes. Karasaki et al. developed a phylogenetic approach using information from changes in genomic distances that exists within different regional growth patterns in the tumor. By combining morphological, genomic, and clinical data, both groups concluded that even though different growth patterns were present within the same tumor, the most aggressive type or high-grade tumor patterns were associated with a high metastatic potential.

Spatial tumor information allows clinicians to identify the location of malignant clones and researchers often study how the genetic content compares between primary tumor and metastatic sites. Transcriptomic data (e.g., single-cell RNA sequencing data), epigenetic data (e.g., DNA methylation sites), imaging data (e.g., MRI), and genomic data (e.g., WES) have been used to retrieve the location and function of different cell types within a tumor. This data is often collected from different sections in the tumor or metastatic sites that can be compared at the transcriptomic, epigenetic, or genomic level. Again, the biological question becomes important when choosing a data type. Transcriptomic data provides a snapshot of the RNA information expressed by the cells, while genomics data can provide long-term view of gene regulation. However, when identifying and comparing mutations present in different areas of the tumor, genomic data would be necessary. Collecting all these data types is clinically expensive and the results can be complex to interpret. Nevertheless, clinicians could leverage spatial models to identify aggressive sites with metastatic potential and aim to target those initially.

The deterministic modeling approaches presented thus far, neglect the stochasticity of random events including the well-known mutation accumulation that happens in many cancers. Therefore, some groups have focused on incorporating these stochastic elements when modeling cancer systems. While there are some known mechanistic progression pathways for certain cancer types, adenocarcinomas such as esophageal or prostate cancer can arise based on a random sequence of events. Aimed at improving screening techniques, researchers have modeled esophageal cancer as a Markov process. This random process delineates the order of different events in tumor progression, where the probability of each event happening only depends on the previous event. This model is suitable for studying esophageal cancer progression because it mirrors the typical patient experience: starting with gastroesophageal reflux disease, advancing to Barrett’s esophagus, followed by undetected cancer, detected cancer, and ultimately leading to death. The multistage Markov model that researchers have developed for esophageal cancer has been effective at determining the appropriate time for screening and starting treatment (Hur et al., 2010; Curtius et al., 2020; Curtius et al., 2021). Prostate cancer has also been modeled as a Markov process across different stages during disease progression. Peirolo et al. studied the transition between well, moderately, and poorly differentiated prostate cancer cells. These cell states were related to a lower and greater stage prostate cancer, respectively. This model estimated cancer progression and identified the grade of the disease using clinical data including survival probability, frequency of detection given a grade and a volume, and probability of metastasis (Peirolo and Scalerandi, 2004). Gastric cancer can also transition through increasingly malignant states from normal gastric mucosa, to chronic nonatrophic gastritis, gastric atrophy, intestinal neoplasia, dysplasia, and finally to gastric cancer. Most of these transitions are determined by a bacterium called H. pylori. Yeh et al. developed a Markov transition model, where an empirical calibration approach was used to estimate good-fitting parameters in line with epidemiologic data obtained from the literature. This model aimed to explore the relationship between H. pylori and gastric cancer to develop strategies for cancer prevention (Yeh et al., 2008).

While many studies have focused on the different stages of cancer as their basis for their Markov model, Veestraeten et al. developed a continuous-time stochastic process for prostate cancer progression based on the variations of the prostate-specific antigen (PSA). This stochastic variable reached a lower boundary after therapy, where progression of the disease stopped. However, during relapse, the PSA level increased. The differential equation model created by the researchers considers a Wiener process, also known as Brownian motion, able to describe the stochastic nature of cancer progression. The parameters included in the model explained how PSA levels were affected by the growth of the tumor cells, and cell death due to various treatment intensities. This approach can help create personalized treatment schedules and track the effects of therapy in the progression of the disease (Veestraeten, 2006).

Another open question in cancer research focuses on the stochastic differentiation trajectories that cells may have travelled along to develop into cancer. Lopez et al. investigated the effect of harmful mutations on cancer evolution by examining whole genome doubling. The authors analyzed sequencing data from lung and triple-negative breast cancers to identify the timing of mutations relative to whole genome doubling. The model used a Gillespie algorithm, a stochastic method leveraged here to created simulations of cancer progression based on mutation rates, the sizes of driver and passenger mutations as well as their effect on fitness (López et al., 2020). Lopez et al. aimed to investigate stochastic mutational events in cancer genes that could be suppressed. Similarly, Lakatos et al. developed a stochastic branching process of tumor evolution to evaluate the production of neoantigens, specific proteins produced by the tumor. A kinetic Monte Carlo algorithm was used to simulate the evolution of colon, stomach, and endometrial cancers and their associated antigen production. In this case, identifying these mutated proteins could be targets for immunotherapies. However, there are still challenges when it comes to identifying the appropriate therapy based only on tumor size and neoantigen presence (Lakatos et al., 2020).

The presence of certain proteins or biomarkers in a tumor microenvironment can inform effective treatment strategies. Several groups have focused on developing molecular functional networks to identify possible interactions between different molecules that exist in the microenvironment. For instance, Chen et al. investigated the presence of biological molecules in lung adenocarcinomas by analyzing a network of genes, microRNA (miRNA), and long non-coding RNA (lncRNA) data linked to various stages of lung cancer. They utilized random walk and Monte Carlo algorithms to group these biological molecules, selecting certain clusters as modules due to their significant association with biological functions. These selected modules, along with their primary biological functions, were integrated to form a core evolution network. This network depicted the progression of lung cancer based on the biological molecules and their functions across the four disease stages (Chen et al., 2022). Jiang et al. also studied lung adenocarcinomas and the point of transition from epithelial to mesenchymal phenotypes in those tumors. Using gene expression data from lung adenocarcinoma cells where the protein TGF-beta has been known to induce epithelial to mesenchymal transition (EMT), researchers developed a Dynamic Network Biomarker (DNB) model. This model identified the timing and components of the EMT transition by examining high variation in the expression of DNB genes (measured as a coefficient of variation), strong correlation between DNB members (measured as the Pearson correlation coefficient), and weak correlation between DNB members and non-members (also measured as the Pearson correlation coefficient). These metrics were combined to calculate a composite index indicating the stage of DNB that led to the transition tipping point. SMAD7 and SERPINE1 were found to be promoters of EMT in lung adenocarcinomas (Jiang et al., 2020). While these dynamic networks can be used to provide insights into cancer staging, they can also inform patient survival. Fang et al. developed a survival network called ‘dnet,’ which uses gene mutations and survival data from breast, colon, lung, and ovarian adenocarcinomas to identify which gene mutations were related to patient survival. Briefly, dnet used TCGA mutation and survival data from cancer patients. This data included mutation frequencies across patients for different genes, along with three clinical variables: age, gender and tumor type. Dnet conducts survival analysis based on the Cox proportional hazards model to evaluate the influence of mutations on patient survival. The Cox proportional hazard ratio (HR) and the associated p-values for each gene were calculated using a likelihood ratio test between the full Cox regression and the baseline regression, identifying potential predictors for various cancer types. Lastly, a patient-survival gene network was constructed from the highest-scoring subgraph determined by the Cox proportional HR p-values previously calculated. This network was used to identify different genes related to patient survival, which could be potential therapeutic targets in almost any type of cancer (Fang and Gough, 2014).

Phylogenetic analysis and the development of phylogenetic trees have been widely applied to understand cancer progression and evolution. Starting from genomic or transcriptomic data, phylogenetic trees can be reconstructed using various criteria, revealing different aspects of tumor evolution. Many of the papers reviewed used next-generation sequencing (e.g., WES, whole-genome sequencing (WGS), DNA sequencing, RNA sequencing, cytogenetic data) data to unravel copy number variations (CNV) or single nucleotide variants (SNVs) driving disease development, which inform the emergence of new branches in phylogenetic tree models (Turajlic et al., 2015). Both CNVs and SNVs are common types of genetic alterations that appear in cancer. CNVs are changes in the number of copies of a specific segment of DNA, which leads to the amplification or deletion of genetic material. On the other hand, SNVs represent changes in single nucleotides in the DNA such as substitutions, insertions, or deletions. Both variations affect gene expression, and protein function, leading to the development and progression of cancer. The detection of these genetic alterations in patients contributes to the construction of phylogenetic trees that can inform cancer pathogenesis and the identify potential therapeutic targets. To study these alterations and reconstruct tumor phylogenesis, researchers have used the following algorithms: Bayesian or stochastic approaches, tumor clonality, maximum parsimony algorithms, clustering methods, or a hybrid method. While established tree construction algorithms such as BEAST or CLONET have originated from using these approaches, other algorithms are still in development. The process of creating new phylogenetic algorithms is closely tied to availability of data for the cancer and treatment of interest. Here, we aim to explore how traditional and new phylogenetic methods have been used to study cancer evolution and identify gene targets that might slow cancer progression and prevent emergence of resistance.

Researchers have also leveraged different diversity indexes from ecology to quantify cellular and phenotypic heterogeneity in tumor samples. These can help assess the richness, evenness, and dominance of different cell types within a tumor sample. Richness refers to the amount or count of cell types present in a sample, while evenness represents the abundance of a certain type compared to the total. Dominance describes the type with the largest influence and control of the overall population. The Shannon index (Shannon, 1948) is the most commonly used metric to quantify diversity within a dataset, considering the uncertainty associated with the presence of specific cell types. A more general approach is leveraging the generalized diversity index (GDI) (Jost, 2006), which can identify diversity across a number of scales of diversity within a dataset (Roswell et al., 2021; Crupi et al., 2019). These indexes have been applied to various cancer types. For instance, in clear-cell renal cell carcinoma, GDI was used to evaluate immune cell infiltration from patient bulk RNA sequencing data. Additionally, in patients with chronic myelomonocytic leukemia (CMML), researchers use the Shannon index to quantify cytokine receptor diversity in different stem and myeloid progenitor populations from patient samples. Higher Shannon indices were associated with an increase in the cytokine receptor diversity in hematopoietic stem cells in patients with CMML compared to healthy patients (Ferrall-Fairbanks et al., 2022a). While the Shannon index has been informative, the GDI emerges as a more robust metric. GDI is constructed as a continuous and non-increasing function within a range of values that are determined by the scale parameter q, which is the order of diversity (Crupi et al., 2019). Lower order metrics favor clonal richness, whereas higher order metrics give greater importance to the most abundant population. Ferrall-Fairbanks et al. first applied GDI to quantify levels of heterogeneity in cancer from single-cell RNA-sequencing data, extracting insights into cancer evolution (Ferrall-Fairbanks et al., 2019). GDI was also used to study parasexual recombination in breast cancer, a process involving the exchange of genetic material between cells without meiosis, where 2 cell lines fused to produce hybrids. Miroshnychenko et al. showed that at early passages, hybrids have a higher GDI across all orders (values of q) compared to the parental cell lines. At later passages, the GDI decreased at low orders, meaning that the species richness decreased, while at high orders was still high, which suggests that the species evenness in the hybrid cells was higher than the parental cells (Miroshnychenko et al., 2021). GDI was also applied to quantify immune cell infiltration using CDR3 sequences recovered from bulk RNA-sequencing of clear cell renal cell carcinoma samples. The authors showed that larger and more advanced disease stage cancers have increased richness in tumor recovered CDR3 sequences, while the evenness of these distributions segregated patients based on survival (Ferrall-Fairbanks et al., 2022b). As seen in this section, the diversity indices can be useful to assess the heterogeneity in a patient population based on clinicopathologic features, which could predict disease progression and treatment response.

EGT has been used to study cancer progression by identifying different players and the payoffs in a cancer system, often the tumor microenvironment. Archetti et al. studied a group of engineered cell populations as a public goods game. In their framework, they engineered cancer cells by knocking-out genes that produce growth factors needed to live and proliferate. By combining non-engineered and engineered cancer cells in-vitro, the overall population growth decreased compared to non-engineered cancer cells alone. This occurred because the engineered cancer cells disrupted the cooperative nature of the population by consuming growth factors (public goods) produced at a cost by the non-engineered population. This cost led to complete or almost complete eradication of the main growth factor-producing clone, resulting in a reduction in population growth. This model was applied to different types of cancer cells such as non-small cell lung cancer and pancreatic neuroendocrine cancer (Archetti, 2021).

EGT has often been paired with some of the methodologies already described. Malekian et al., for instance, created a model that combines EGT concepts with an ABM to study the effects of gap junctions in ductal carcinoma in situ progression. The agent-based portion of the model considers each cell to exist in one of the following states: quiescent, proliferative, apoptotic, and necrotic. These cells can transition between states based on oxygen concentration and the phase of the cell cycle. The EGT portion of the model was designed to address the interactions between cells. Here each agent, or cell, represents a player that can interact with any of the 8 cells that surrounds it using one strategy at the time (quiescence, proliferation, apoptosis, and necrosis). In their model, the interaction between 2 cells through gap junctions created a pay-off that sends the cell a survival or death signal based on different factors in the microenvironment such as oxygen levels or growth factors concentrations (Malekian et al., 2016).

Zhang et al. combined EGT concepts with a Lotka-Volterra model. The authors aimed to explain the competition of different cells in metastatic castrate-resistant prostate cancer under adaptive therapy. Their players were the cancer cells, the strategies the different phenotypes that a cell can take, in this case androgen-dependent, testosterone-producing, or androgen-independent cells, and lastly the payoffs were proliferation and survival of the cells. The Lotka-Volterra model addressed the competition of the three different cell phenotypes while undergoing therapy (Zhang et al., 2017).

Spatial data can also be modeled using game theory. Kaznatcheer et al. defined cancer cells as players that can take one of two different strategies: those exhibiting autonomous growth that are capable of rapid proliferation (AG cells) and those with acquired motility and invasiveness conferred by mutations (INV cells). These phenotypes were tied to different payoffs: the cost of motility on INV cells, and the benefit of fitness of a cell that has every resource available. For example, if two INV cells meet, then one gets the benefit of using the surrounding resources and the other one migrates at a cost. However, if a INV cell encounters an AG cell, then the INV cell pays the cost of having to move to find resources, while the AG benefits from that interaction. Using data from prostate cancer, researchers were able to justify that the spatial location of the cells can affect tumor invasion and evolution based on these cellular interactions (Kaznatcheev et al., 2015).

The Gini index was originally developed as a measure of economic inequality. In the context of cancer, it can be used to quantify cellular heterogeneity within the tumor to predict disease progression and treatment response. Hinohara et al. used the Gini index to quantify heterogenic expression of specific markers in breast cancer. Specifically, this index was used to identify variations in the expression levels of a group of enzymes (KDM5) implicated in cell proliferation and cancer progression in breast cancer cell lines, both pre- and post-treatment with an inhibitor targeting that protein family. They concluded that inhibiting KDM5 enzymes decreased transcriptomic heterogeneity, which could increase treatment effectiveness (Hinohara et al., 2019). Also in breast cancer, Filho et al. used the Gini index to quantify heterogeneity of ERBB2 gene expression in HER2 positive patients. In this type of cancer, an intratumor heterogeneity expression in the ERBB2 gene is associated with cancer resistance. However, the researchers observed a relatively low Gini index, suggesting that the samples have an evenly distributed expression level of the ERBB2 gene. In this case, the Gini index was not a good predictor of drug response (Filho et al., 2021). This index can also be used to quantify expression of markers that suggest response or resistance to treatment. Gil Del Alcazar et al. used the Gini index to quantify changes in immune populations between stable and growing tumors in a breast cancer rat model. Stable tumors expressed a higher Gini index with respect to their B-cell and T-cell receptors. This implied a higher diversity in immune response making these tumors not suitable for immune therapies (Gil Del Alcazar et al., 2022).

The Gini index is usually used in combination with other diversity metrics. Miroshnychenko et al. used the GDI to quantify diversity between parental breast and epithelial cell types, and hybrids of the two, which could have been formed through cell fusions also known as somatic hybridization. While the GDI can differentiate cells into different phenotypes, the Gini index can capture cell-to-cell phenotypic variability. Therefore, the Gini dispersion index provided information about the gene expression variability of cells within the same phenotype. The researchers found that hybrid cell populations had higher Gini indexes compared to the parental populations indicating that somatic hybridizations increased phenotypic diversity (Miroshnychenko et al., 2021).

The application of economic mathematical principles has been useful to study cancer in several ways. Game theory has been able to explain different phenomena happening during cancer development such as resistance, cell interaction, and the emergence of different cell states. These processes were modeled as a game with players and pay-offs tied to different strategies. Game theory, in combination with ecology tools, has been able to capture additional cellular interactions creating a more robust framework to explain temporal data. Additionally, the Gini index originally developed to measure economic inequality, has been applied to quantify cellular diversity in tumors. Exploring the cellular composition indicates effectiveness of treatment. For instance, high diversity in immune populations within a tumor indicate that immunotherapy treatment would be ineffective. Data types such as population counts, or volume are necessary to utilize game theory models while sequencing data would be necessary to obtain a Gini index. Once again data availability becomes a challenge when using these tools to study cancer progression or treatment effectiveness.