Living organisms, intricately crafted by nature, exist as remarkable ensembles of cells, harmoniously navigating the intricate interplay between their internal mechanisms and the external environment. Fundamentally, cell division proceeds in an orderly and tightly regulated manner, guided by precise sequencing rules essential for the harmonious functioning of the organism's physiological processes [1]. However, on certain occasions, this delicate balance falters, and cells deviate from their normal patterns of growth and departure. Instead, they embark on a path of uncontrolled proliferation, exhibiting an aggressive behavior that disrupts the equilibrium between cell birth and cell death [2]. Such aberrant cells are commonly known as tumor or cancer cells [2]. Among the spectrum of life-threatening diseases, cancer stands as the second most formidable adversary, claiming countless lives worldwide [3]. Archaeological investigations spanning back 1.7 million years have yielded compelling evidence indicating that our human ancestors were not immune to the ravages of bone cancer [4]. These findings shed light on the ancient origins of this disease and provide insights into its enduring presence throughout human history [4]. Remarkably, the footprint of cancer transcends human history, extending back millions of years. Evidence derived from fossils dating as far back as 75 million years ago, including those of dinosaurs, has unveiled instances of bone cancer in various ancient species [5]. Projections indicate that in the year 2024, a staggering number of approximately 2,001,140 new cases of cancer will be diagnosed in the United States alone [6]. Furthermore, it is estimated that around 611,720 individuals will succumb to cancer-related causes, translating to an average of 50,000 deaths per month [6]. These alarming statistics underscore the critical importance of advancing research, implementing effective prevention measures, and providing optimal care to combat this pervasive and devastating disease.

Biological systems possess an inherent dynamism that is both intricate and captivating. Whether we examine population dynamics, inheritance patterns, or the progression of diseases, these systems undergo constant change, modification, and transformation. Nevertheless, achieving a thorough comprehension of these systems poses significant challenges, primarily due to the complex nature of cancer biology. This complexity stems from various factors [7], including genetic heterogeneity, cell-cell interactions within the microenvironment, the process of angiogenesis, and nonlinear dynamics.

Genetic mutation

Cancer cells undergo rapid mutations and evolution, resulting in a diverse array of cells within a tumor. Each cell may harbor distinct genetic mutations, thereby influencing its behavior and response to therapeutic interventions. Studies, such as the one conducted by Haynes et al. [8], have demonstrated that tumor stem cells and progenitor cells within tumor ecosystems follow predetermined paths of differentiation. These paths contribute significantly to the establishment and perpetuation of heterogeneity within breast tumors.

Cell-cell interactions

Tumors exist within complex microenvironments comprising various cell types, including immune cells, stromal cells, and blood vessels. The interplay between cancer cells and their microenvironment exerts significant influence over tumor progression, invasion, and response to treatment. Galon and Bruni [9] discussed upon the multifaceted roles played by various immune cells at different stages of tumor development and advancement.

Angiogenesis

The process of forming new blood vessels from existing ones is pivotal in cancer progression. Tumors release signals that stimulate angiogenesis, facilitating their growth and dissemination by acquiring essential nutrients and oxygen from the bloodstream. Heredia-Soto et al. [10] elucidated the critical roles of angiogenesis and the tumor microenvironment (TME) in ovarian cancer.

Nonlinear dynamics

The process of tumor growth and progression is complicated, as it involves interactions among numerous elements such as tumor cells, the surrounding microenvironment, immune cells, blood vessels, and signaling molecules. These interactions frequently demonstrate nonlinear behaviors [11], wherein minor alterations in one component can trigger significant effects on the entire system.

These complexities make it difficult to model tumor dynamics accurately because traditional mathematical models often oversimplify the underlying biology. However, advanced mathematical modeling provides a powerful tool to capture the complexity of biological systems and make robust predictions about their behaviour, for instance, Hormuth et al. [12] pioneered the development of an image-based mathematical model aimed at investigating the dynamics of tumor growth and the vascular response to angiogenesis within the context of radiation therapy integration.

The hallmarks of cancer listed in [13,14], briefly explained in Table 1, are set of characteristic features or behaviors that cancer cells typically acquire during tumorigenesis (the process of tumor formation and growth). While these are useful, but not sufficient to provide a full appraisal of the many aspects of cancer biology in all dimensions. More information is required to understand cancer from all angles to furnish as detailed as possible a full perspective of how diverse cancer disease is. Unravelling the precise influence and extent to which these fundamental hallmarks impact tumor growth, dynamic behavior, and response to treatment poses a formidable challenge. However, elucidating these complex mechanisms required interdisciplinary collaboration, advanced mathematical modeling, and comprehensive experimental studies. Bull and Byrne [15], in their work, provide valuable insights into the hallmarks of mathematical oncology;

Single versus hybrid frameworks: In single frameworks, focusing on one aspect of cancer dynamics, while in hybrid frameworks that combine multiple factors. Single frameworks simplify analysis, while hybrid frameworks offer a more comprehensive understanding of cancer's complexity by integrating various biological processes.

Homogeneity versus heterogeneity: Cancer cells and their environment can be depicted as either homogeneous (uniform) or heterogeneous (diverse). Homogeneous models assume uniform behavior across all cells, while heterogeneous models account for differences among cells, such as genetic mutations or responses to treatment.

Spatially averaged versus spatially resolved: Spatially averaged models provide an overview of cancer dynamics across a tissue or tumor without considering spatial variations. In contrast, spatially resolved models capture spatial differences in cell behavior, interactions, and environmental factors within the tumor microenvironment.

Single-scale versus multi-scale: Tumor models can operate at a single scale, focusing on specific biological levels (e.g., molecular, cellular) or multi-scale, integrating processes across different scales simultaneously. Single-scale models simplify complexity, while multi-scale models capture interactions between scales for a more comprehensive understanding.

Deterministic versus stochastic: In deterministic models, cancer behaviors are predicted based on fixed parameters and initial conditions, assuming no randomness. Stochastic models incorporate randomness into predictions, considering probabilistic events such as mutations, cell division, or treatment response, which better reflects real-world variability.

Continuum versus discrete: Continuum models represent cancer dynamics using continuous mathematical functions, ideal for describing large populations of cells. Discrete models treat individual cells or entities as distinct entities, suitable for studying interactions at a granular level, such as cell-cell interactions or genetic mutations.

They outline various aspects of these hallmarks, such as the modeling decisions made, the incorporation of biological data, the level of complexity in finding solutions, and the quality and nature of the data used to generate the models. By blending these hallmarks, researchers create varied mathematical models to grasp the intricacies of biological systems comprehensively.

In this review, the necessity of mathematical models to capture the complexities of tumor microenvironments is emphasized. The comparison of basic tumor growth models, which serve as the foundation for contemporary advanced mathematical models, is discussed. These models are instrumental in understanding the dynamics of tumor growth and progression. Mathematical models, ranging from discrete to advanced multiscale hybrid models, are explored. The challenges arising from the complexity and vast heterogeneity within tumor microenvironments are highlighted. Different mathematical treatment models are presented as solutions to these challenges. These advancements include cutting-edge models developed to capture the complexities of both traditional and emerging cancer treatments. By synthesizing these findings, this review provides a comprehensive overview of the current landscape of mathematical modeling in the context of tumor dynamics and treatment. It underscores the significance of incorporating mathematical models to enhance our understanding of tumor microenvironments and drive advancements in cancer research and treatment strategies.

Why the need for mathematical models?

Cancer is widely acknowledged as a complex system by the oncologist. These systems have a biological, chemical, and mechanical structure that changes over time to support tumor growth [7,8]. However, mathematicians have traditionally focused on simplified aspects of these systems. Every single day, there have been efforts to create comprehensive descriptions of tumor growth and treatment, including complex phenomena like cell-cell interactions, interactions with the extracellular matrix, and cell mutations [7]. Mathematicians have proposed various solutions using different methods, making the field of mathematical oncology diverse [15-17]. Perhaps one of the earliest studies by Mayneord [18], the mathematical model equation describing the growth of Jensen's rat sarcoma with irradiation was formulated, revealing a linear increase in its linear dimension i.e., measurements of the tumor's length and breadth. However, it's challenging to capture the full complexity of cancer in mathematical models (illustrated in next section). Additionally, proposing new models requires caution because including more phenomena makes the models more complex and requires knowledge of many physical parameters that are hard to measure experimentally [7]. Generalizing the behavior of the entire tumor microenvironment through mathematical models leads to the formulation of nonlinear higher-dimensional dynamical systems [11]. This highlights the importance of collaboration between mathematicians and oncologists.

We must acknowledge the significant contributions made by mathematicians (we will be seeing in coming sections). Mathematicians have also created practical tools for therapy procedures. However, there's still much to do to accurately predict tumor growth rates, metastasis risks, and the effects of drugs or radiation. Using computers for simulations can make experiments faster and more efficient. It allows us to simulate treatment cycles virtually, helping decide the best way to administer them. For instance, what's the ideal duration, intensity, and frequency of radiation treatments? Even a rough model can provide valuable insights. Our motivation in examining the current state of mathematical cancer research critically is to inspire new ideas, especially among young researchers.

Mathematical models

A plethora of tumor growth models have been proposed by various researchers to gain a more accurate understanding of the dynamics and progression of tumor growth. These models usually use differential equations (DEs) to explain how the tumor grows or changes over time. For instance, they help to figure out how many tumor cells there are at different points in time [16,19]. Differential equations offer a powerful mathematical framework to capture the complex interactions and dynamic changes observed in tumor growth. Some widely used mathematical models for studying tumor growth dynamics include Linear, Logistic, Exponential, Exponential-Linear, Gompertzian, Power-Law, Von Bertalanffy, and Dynamic Carrying Capacity [16]. In tumor growth models, carrying capacity refers to the maximum tumor size or population density that can be sustained within a given environment or host [16]. It represents the point at which growth slows down or stabilizes due to limitations in resources such as nutrients, oxygen, and space [20]). In Table 2 we have compared all these model equations, their advantages and disadvantages, goodness of fit (best fit for), identifiability, and predictability.

Model validation and accuracy

Model validation involves assessing how well a mathematical model represents the observed dynamics of a system using experimental data, often requiring collaboration with experimentalists to ensure accuracy and reliability. This process typically includes comparing model predictions with measured data, either from the same experimental conditions used for parameter estimation or from independent experiments [16,17]. A reliable model is expected to show strong alignment with observed data, with minimal errors and high predictive capability.

The validation process often evaluates the model’s goodness of fit (best fit for), which determines how closely the model reproduces the observed data. Metrics such as root mean squared error (RMSE) are commonly used to quantify the difference between predicted and actual values, with lower RMSE values indicating better fit. Additionally, statistical methods, such as analyzing the residuals or using coefficients of determination (R2), help assess the accuracy and reliability of model predictions [21]. Another critical aspect of validation is the reliability and identifiability of model parameters. In reliability, a robust model produces parameter estimates that are consistent and supported by confidence intervals and standard errors [16]. Identifiability ensures that the parameters can be uniquely determined from the available data, which is crucial for avoiding over parameterization and ensuring that the model remains interpretable and predictive [22]. Predictability is another key factor, reflecting the model's ability to forecast system behavior over time [23]. High predictability indicates that the model can generalize well beyond the data used for parameter estimation. Furthermore, sensitivity analysis is often employed to evaluate how changes in individual parameters affect the model’s output. This analysis helps identify critical parameters that significantly influence the system’s behavior and can guide model refinement by eliminating redundant or insignificant parameters.

A validated model should demonstrate strong goodness of fit, reliable and identifiable parameters, accurate predictions, and robust sensitivity characteristics. Such a model provides valuable insights into tumor dynamics and can guide treatment decisions with greater confidence.

Several studies compared different tumor growth models using both clinical and experimental data [16,17]. These studies found that there's a lot of variation in the parameters for each model, which was mainly due to differences in individuals' growth rates caused by factors like genetics, environment, and tumor characteristics like vascularization (the process where new blood vessels form within the tumor). Figure 1 shows a comparison of these growth models considering the estimated parameters in [16]. For instance, Sarapata and Pillis [17], looked at ten types of cancers and fit each model to data for each cancer type. Similarly, in [16], it provided a detailed analysis of experimental mouse models focusing on lung and breast cancers. They used data to compare different growth models in terms of how well they fit the data, their predictability, and how easy it was to identify their parameters. After comparing the models listed in Tables 2 for their descriptive ability, strengths, and weaknesses, it's evident that there isn't one model that fits all types of cancers perfectly. However, based on the analysis in references [16] and [17] and the summaries in these tables, the Gompertz model appears to provide a reasonable fit for cancers in the lung, breast, head and neck, bladder, liver, and pancreas.

Similarly, the logistic model also fits well for different types of cancers. It shows good identifiability [16] and accurately predicts lung cancer progression in humans [24], based on data from breast and lung cancers. The power-law model is the best fit for many cancer types, including breast and lung cancers. However, it's highly sensitive to changes in its parameters, particularly  and . Increasing  even slightly leads to biologically unjustifiable tumor growth. Despite being the best fit for most cancers, it's not recommended due to its unstable nature. Comparing the Gompertz and power-law models, the power-law model fits most cancer data best. However, because of its sensitivity to parameters and biologically questionable results, it's better to use the Gompertz or logistic models instead.

As basic models have limitations in fully illustrating tumor dynamics, continuous advancements have led to the development of various types of tumor growth models. These models aim to address some, but not all, complexities involved in tumor growth dynamics. They are valuable tools for understanding tumor growth, progression, and treatment response. In this section, we will explore different types of tumor dynamics modeling approach, such as discrete, continuous, hybrid, agent-based, stochastic, and multi-scale. Each of these approaches offers a unique perspective on tumor dynamics:

Discrete modeling techniques offer a means to describe the spatial and temporal dynamics of individual cells within a population [26]. In this approach, mathematical equations are employed to capture the behavior of each cell individually, allowing for the tracking and updating of internal interaction mechanisms that govern their physiological processes. By defining a discrete model in a time-dependent manner, the system's state can be described at specific time points [26]. These models are particularly suited for in vitro cellular studies or smaller-scale modelling applications because individual cells operate independently under a specific set of rules [27], such as cell-cell and cell-matrix interactions, including chemotaxis and haptotaxis, allowing for an examination of tumor dynamics at the level of single cells and below [28]. However, this approach is restricted to a smaller scale due to the high computational expenses involved [29]. Despite their effectiveness, discrete models are constrained by their intensive computational requirements, typically restricting their domain size to sub-millimetre scales or even smaller [30]. Dehingia et al. [31] proposed a discrete mathematical model to study the tumor-immune interaction in presence of a discrete time-delay term, given as:

Where,  stands for time;  for number of tumor cells;  for number of hunting T-cells (T-cells are a type of white blood cell crucial for the immune system's defense against infections and abnormal cells like cancer); , number of resting T-cells; Δ, discrete time-delay term for the growth dynamics of hunting T-cells. In Equation 1, the rate of change of cancer cells depends on several factors. The first term represents a linear increase, where s is the constant conversion rate of normal cells to tumor cells. The second term accounts for the logistic nature of growth dynamics. The third term reflects exponential decay, with α₁ representing the killing rate of malignant tumor cells by haunting T-cells. Similarly, equations 2 and 3 depict the rates of changes of haunting and resting T-cells, respectively, where the parameter Δ indicate the discrete-time delay in the growth of hunting T-cells. Additionally, it's assumed that tumor-specific resting T-cells can proliferate in the presence of cancer cells. This mechanism is represented by the Michaelis-Menten term , where ρ gives proliferation rate and η gives the steepness coefficient [32]. The simulation findings suggest that as the growth of hunting T-cells gradually increases within the patient's body, the immune system becomes more effective in stabilizing the growth of tumor cells. Even if there's a minor delay in this growth process, the patient's immune system may still successfully regulate tumor growth. However, over time, all three cell types start competing with each other, eventually leading to a Hopf bifurcation (it describes a sudden change in the behavior of a system as it reaches a critical point). The theoretical analysis of the model was validated numerically with the findings of Kaur and Ahmad [32] study. For other parameter definitions and details reader can follow [31].

Continuous-based modeling approaches, on the other hand, study tumor tissue entities as continuous variables, such as concentration, density, and volume fraction [33]. Models formulated using ordinary differential equations (ODEs) [34], partial differential equations (PDEs) [34], or integro-differential equations (IDEs) [35], are used to describe population-averaged or spatially determined dynamics of tumor cells. ODE models focus on temporal dynamics, making them suitable for analyzing time-resolved data ranging from the cellular to tissue scale. PDE models, on the other hand, capture the spatial dynamics of tumor growth and its invasion into neighboring healthy tissues, considering various factors such as oxygen and nutrient concentrations, vascular angiogenesis, tumor angiogenic factor concentration, drug concentrations, and haptotaxis effects [34]. These models are computationally efficient and can encompass a larger length scale volume fraction of tumor dynamics at a time [36]. Additionally, continuum models allow for the implementation of fundamental physical principles of continuum mechanics [37]. Anderson et al. [33] presented a continuous mathematical model having system of non-linear partial differential equation to illustrate the tumor cells-ECM (extracellular matrix) interaction response through haptotaxis consideration, presented as follows:

Where  stands for tumor cell density;  , ECM density;   , matrix degradation enzymes (MDEs) concentration. This model focuses only on how the tumor interacts with the tissue around it, without considering cell proliferation. Tumor cells produce MDEs that break down the surrounding tissue [38], creating space for the tumor cells to move through by diffusion (first term right hand side eq. 4). The model also takes into account how tumor cells move, which is called haptotaxis (second term in eq. 4). This means they respond to gradients of certain molecules like fibronectin. The model assumes that when MDEs come into contact with tissue, they start breaking it down as reported in [38]. The amount of ECM degraded depends on the concentrations of both ECM and MDE, with δ representing the rate at which ECM is degraded by MDE (illustrated in eq. 5). The role of MDE is to break down the basement membrane, allowing cells to move into the surrounding ECM. However, MDE can also target proteins besides ECM, like growth factors and cytokines [38]. This makes MDE activity complex, and it can either help or hinder tumor growth [38]. But in this model, it focusses only on MDE's role in breaking down ECM produced by tumor cells. This degraded ECM then spreads into the surrounding ECM and eventually breaks down, where  diffusion coefficient,  production rate, and  decay rate of MDEs (presented in eq. 6). For other parameters definition and detail reader can follow [33].

However, it is important to note that continuum models may not capture the complete range of diverse cellular and subcellular dynamical features exhibited at the single-cell level when individual cell mechanisms dominate [39]. Moreover, they may not fully capture tumor heterogeneity when cell properties vary over small spatial scales, such as genetic features of tumor cell heterogeneity [40]. Discrete models, on the other hand, allow for a more detailed examination of individual cell behavior but can be computationally intensive and limited in their ability to describe larger-scale systems [29,30]. Hybrid models have emerged to address these limitations and harness the benefits of both approaches.

Hybrid modeling, integrates both discrete and continuous approaches [41]. In this, the tumor microenvironment's extracellular components, such as oxygen concentration, drug concentration, and nutrients concentration, are modeled as a continuous medium, while individual cells and sub cells (e.g., tumor cells) are treated as discrete entities [42]. By combining discrete and continuum elements, hybrid models bridge the gap between the cellular and tissue scales, enabling a more comprehensive understanding of tumor behavior. Jeon et al. [43] extended the PDE model by Anderson et al. [33] to a hybrid mathematical model to demonstrate the tumor growth dynamics and its invasion. In the continuous part of the model, it added the concentration of one more microenvironment nutrient component, oxygen. Whereas the discrete part of the model describes individual cell behaviors like cell cycle, interactions between cells and with the ECM, and how cells move. For instance, here, in the discrete part we only present the cell movement equation.

Continuous Part of the model (concentration of oxygen component): Oxygen is crucial for tumor cells to survive, grow, and spread. This assumes that blood vessels, which provide oxygen, have formed and are continuously supplying oxygen through the extracellular matrix ECM. This also assumes that the amount of oxygen supplied by the ECM is directly related to how dense the ECM is, following previous research by Anderson and others [44]. Oxygen then diffuses into the ECM where tumor cells use it, and it naturally breaks down over time. This process is described using a mathematical equation that accounts for both diffusion and chemical reactions as:

Where,  stands for nutrient concentration (oxygen),  for time.

Discrete Part of the model: Cell movement is often described as a persistent random walk, which can be represented by a mathematical equation called the Langevin equation. In this equation (8), the first term accounts for friction between cells and between cells and ECM, while the second term represents random forces that affect cell movement, and the third term represents other associated extra forces values acting on the cell.

Where,  stands for mass of a single cell,  for cell velocity. For other parameters definition and detail reader can see [33,43].

Agent-based models (ABMs) serve as a valuable tool in modeling by simulating the behavior and interactions of individual entities, referred to as agents, within a system. These agents typically represent various cell types found in the tumor microenvironment, such as cancer cells, immune cells, and stromal cells. ABMs possess several key characteristics: Firstly, they are individual-based, meaning each agent possesses its own unique properties and behaviors, which may evolve over time based on interactions with other agents and the environment [45]. Secondly, ABMs are spatially explicit, allowing for the modeling of the spatial arrangement of cells and their impact on their interactions, thereby aiding in the understanding of tumor growth and invasion into surrounding tissues [46]. Lastly, ABMs are stochastic in nature, incorporating randomness into agent behavior to reflect the inherent variability observed in biological systems [47]. Agent-based models can be broadly categorized into two main types [48]: on-lattice based models and off-lattice models.

On-lattice models employ a fixed grid system, which simplifies simulation and computation [48]. However, this approach confines cell movement and interactions to neighboring grid sites, limiting their ability to fully represent the complexity of real tumors. In contrast, off-lattice models offer a continuous space, allowing for more realistic cell shapes and movements [48]. While this flexibility enables better representation of real-world phenomena, it comes with increased computational demands and potential challenges in managing cell overlaps and movement patterns. On-lattice models are particularly efficient for studying large-scale phenomena, whereas off-lattice models excel in capturing the sophisticated details of microscopic processes and the tumor microenvironment. Radszuweit et al. [49], developed an on-lattice based models to study tumor growth kinetics and compared with the in vivo Xenograft experimental results of mouse fibroblast cell line. Rocha et al. [50] studied an off-lattice agent based model at the cell level to describe dynamics of normal and tumor cells, accounting for cell movement and phenotypic transitions driven by microenvironment stimuli. While ABMs provide valuable insights into complex systems such as tumor dynamics, they are not without limitations. One significant challenge is the high computational run time associated with simulating the behavior and interactions of numerous individual agents, particularly in large-scale models [51]. Many research groups have released specialized tools for agent-based modeling in biology, however many researchers prefer to write and maintain their own code for agent-based simulations. For agent-based simulations, these tools include Physicell [52], CompuCell3D [53], Chaste [54], Morpheus [55], CellSys [56], VirtualLeaf [57], SEM++ [58], etc. The reader can refer to section 3, for illustration of an example of agent-based model.

Stochastic modeling, in the field of mathematical oncology, there's been a shift from deterministic to more stochastic modeling [59,60]. This change is thanks to advancements in computing technology. Initially, models relied on deterministic ordinary differential equations (ODEs), which didn't fully capture the randomness in biological systems. Now, with better computing power, stochastic models are gaining popularity. These models incorporate randomness, giving a more realistic view of tumor behavior. This shift allows researchers to better understand complex processes like tumor heterogeneity and drug resistance. Kim and colleagues [59] investigated cancer virotherapy using mathematical models. They constructed both deterministic and stochastic models to understand how different oncolytic viruses, used in cancer treatment, interact with tumors. The stochastic model was developed based on the deterministic one, aiming to analyze tumor extinction probabilities under different parameter conditions. Here, in this, we only consider the model equation for the population dynamics tumor cells , which are virus infected.

Deterministic part:

Stochastic part:

Where,  virus infected tumor cells,  virus,  cytotoxic T lymphocyte (CTLs),  uninfected tumor cells, α infection strength,  cytotoxicity of virus, α₂ other constant,  are components of a Wiener process  (here, Wiener process represent the random fluctuations in tumor growth or response to treatment over time, accounting for the inherent uncertainty in biological systems). Equation 9 shows the rate of changes of infected tumor cells over time, where the first term on the right-hand side represents the process where uninfected tumor cells become infected by contacting the virus at a rate denoted by . The second term accounts for the virus's ability to kill infected tumor cells, while the third term describes how virus-infected tumor cells are eliminated by CTLs at rates indicated by . Equation 10 introduces a stochastic differential equation model, which extends the deterministic model from equation 9. This extension incorporates stochastic processes into the system equations. Specifically, in the model proposed by Dingli et al. [61], free-virus particles are considered relatively small in size, necessitating the inclusion of randomness. In this context, they define  as a continuous random variable representing , where T denotes the transpose of a matrix. The stochastic differential model is then formulated as follows:

Following eq. (10.1), a system of four equations formed for each of the variables, where equation (10) is the first component. For other parameters definition and detail reader can go through [59]. Moreover, readers can follow up the following articles to understand a wide variety modeling structure for tumor growth [33,62-65].

Given the complexity of tumor dynamics and complex interactions among cells in the heterogeneous microenvironment, a multiscale modeling approach has emerged as a valuable tool [66]. This approach combines the advantages of both discrete and continuous modeling approaches to capture the different scales and processes involved in tumor growth and response to treatment. Multiscale models typically employ a hybrid approach, where signaling pathways [67] and cellular processes [67] are simulated using ordinary differential equations (ODEs), while the heterogeneity of the tumor microenvironment, drug concentration, and delivery are described using partial differential equations (PDEs) in the continuum modeling framework [67]. Researches highlighted that cancer encompasses diverse features such as growth, tissue invasion, and metastasis, which operate across multiple spatial and temporal scales [68]. Attempting to comprehensively capture and predict the dynamics and clinical implications of cancer using a single-scale model is an arduous if not insurmountable task [66]. The limitations emphasize the need for advanced multiscale modeling approaches that can effectively elucidate the complexity and heterogeneity of tumors.

In recent years, mathematics in cancer studies has emerged as a pivotal field encompassing studies at various biological scales, bridging the gap between theoretical cancer biology research and clinical oncology [69]. Different fields within life sciences have developed various definitions for biological scales [70]. Among these definitions, some of the primary spatial scales in biology include atomic [71], molecular [72], microscopic [72], macroscopic [73], and organismal/clinical scales [74]. The multi-scale modeling approach employed in this context offers several distinct advantages, making it a powerful tool in cancer research and clinical applications:

Multi-scale models offer profound insights into cancer progression by quantitatively capturing bio-physical mechanisms

Sadhukhan and colleagues [75] devised a computational off-lattice agent based model to explore the complexities of avascular tumor growth within epithelial tissue. Their model, which operates across multiple scales, unveiled the dynamics at the intracellular, cellular, and extracellular levels. Every cell was viewed as an agent with the ability to multiply, give rise to daughter agents that were exactly like it, or change into a different phenotype in response to changes in the external microenvironment and internal protein activity. Furthermore, the model took into consideration a range of biophysical forces that influenced the dynamics of tumor growth. The model addressed key aspects associated with the hallmarks of cancer. This included factors such as oxygen and nutrient utilization, variations in the tumor microenvironment, the proliferation of tumor cells through evasion of growth-suppressing signals, enhanced replication of tumor cells, and their resistance to programmed cell death (apoptosis)

Multi-scale models can simultaneously incorporate multiple spatial and temporal scales

Shipley et al. [76] employed mathematical homogenization techniques to construct a novel multiscale continuum model that investigated the transport of blood and chemotherapeutic drugs within both the blood vessels and the surrounding tissue of a vascular tumor. This model integrated the analysis of three distinct length scales within the tumor, briefly outlined below. Using a simplified representation of a dorsal skinfold chamber as a case study, the model's behavior was scrutinized to elucidate how the architectural layout of the vascular network impacted the distribution of fluids and drugs at the chamber's length scale. Comparative evaluations of various chemotherapeutic strategies from the model simulations indicated that a single injection approach consistently yielded superior outcomes compared to continuous perfusion methods. The model equations given are as follows:

Capillaries scale

Arterioles and venules scale

Chamber scale

 are fluid pressure and  are fluid velocity in the arterioles, venules, capillaries, and interstitium respectively. The interstitium consisted of cells surrounded by extracellular space, the capillaries were considerably larger than the gaps between cells in the interstitium. As a result, the interstitium was treated as an isotropic porous medium, and fluid flow through it was described by the Darcy law [77]. It was also assumed that fluid flow in both spaces was incompressible. Therefore, the model on the interstitium given as in eq. (11), where,  is blood plasma viscosity and  is interstitial permeability.

In larger blood vessels, the study assumed blood behaves as a continuous fluid with constant viscosity [78], also it omitted non-Newtonian effects to reduce complexity. Hence, in the capillaries, they considered the Navier-Stokes equations for fluid with constant viscosity, presented in eq. (12), where ρ is the density of fluid.

The model for fluid transport given in eq. (11), and (12) were then homogenized to describe movement in the arterioles and venules scales. The simplification process involved several steps: (a) assuming a difference in length scale between the micro and arteriole/venule lengths, (b) expanding all variables in powers of a small parameter  (a length scale ratio), (c) matching coefficients of powers of  (d) separating fluid velocity into microscale and arteriole/venule scale components, (e) averaging over the microscale to derive relationships for the arterioles/venules. On the scale of arterioles and venules, the capillaries and interstitium were treated as continuous mediums, behaving like a double porous medium with connected Darcy flow between them. Thus, on this scale, the equation governing fluid movement in the capillaries and interstitium represents as in equations (13) and (14). Lastly, this homogenization extends from the arteriole/venule scale to the chamber length scale, resulting in an effective model on the length scale of the chamber. This process produces equations for pressure and velocities in the arterioles, venules, capillaries, and interstitium as outlined in equations (15-18). The numerical model was compared and validated with the experimental study [79] and observed a good consistency. For other parameters definition and detail reader can further see [76].

Multiscale modeling advances traditional hybrid modeling approaches into a more sophisticated system

Jafari Nivlouei et al. [80] developed a multiscale mathematical model with the aim of examining the biological and physical dynamics of tumor evolution in the presence of normal healthy tissue. This model has considered a range of events occurring across different times and spatial scales. The study focuses on two distinct phases of tumor development; the avascular and vascular phases and investigates scenarios with and without the presence of normal healthy cells. The findings indicated an increase in tumor growth rate with the establishment of closed vessel loops around tumor cells during tumor-induced vascularization. Additionally, by introducing a multiscale simulation that accounts for cellular interactions and different cancer cell behavior, the study demonstrates better efficacy. The outline of model given as:

Agent-based modeling- cellular scale

Continuum modeling- Extracellular scale

The study used an agent-based cellular Potts model to simulate cellular system dynamics, employing a discretized lattice Monte Carlo method developed by Glazier and Graner. This model accommodates spatially extended generalized cells, allowing for the consideration of various levels, from intracellular to tissue scale. Each type of cell is given a special number τ, and each individual cell is marked with a specific value σ = 1,2,... .This helped us keep track of different types of cells and their locations in the simulated space.

The model operates on the principle of energy minimization, with each configuration's energy referred to as the Hamiltonian value . This Hamiltonian comprises four terms governing cellular dynamics: cell-cell adhesion, cell growth, chemotaxis, and cell continuity assurance. Here in this, for understanding we present model equations for cell adhesion only. Cell adhesion, a vital biological and physical cell property, influences tissue integrity and cell motility, fostering interactions between adjacent cells and between cells and the ECM. It is cell-type-dependent and reflects the coupling strength between entities  and the model given as in eq. (19). The sum includes all neighbouring pixels, where σ and σ' represent the cells' ID, and δ Kronecker symbol.

During tumor growth, cells experiencing oxygen deprivation release angiogenic factors like VEGF to promote the growth of new capillaries for accessing nutrients. This process was modeled using a diffusion-reaction equation at extracellular scale, depicted in eq. (20), where  represents nutrient concentration,  is the diffusion coefficient of nutrients, and  denotes the release of nutrients from vessels. The function B  accounts for nutrient uptake by cancer cells. VEGF secretion creates a concentration gradient between the tumor and nearby blood vessels, triggering activation of endothelial cells. The distribution of VEGF is described by a partial differential equation in a similar manner [80]. The simulations of the multiscale model were validated against the standard experimental study [81,82] and observed a consistency. For other parameters definition and details reader may look into [80].

Multi-scale modeling can integrate diverse data and mechanistic processes to systematically understand and depict disease onset and progression [83].

This approach reconciles the complexities associated with different scales, enabling researchers to obtain a more accurate representation of cancer dynamics. Managing vast amounts of data across different scales poses significant challenges. A key challenge in multi-scale modeling lies in achieving accuracy and efficiency in capturing both macroscopic and microscopic features of cancer [84]. Complex system simulation at various scales can be computationally costly [85]. Extensive data and complex calculations demand substantial computing power, often posing limitations, particularly for large-scale models [86]. Establishing an interface between discrete and continuous modeling poses a notable challenge [87]. Maintaining mathematical consistency, preserving the rules and properties of the governing equations, and seamlessly transitioning between discrete and continuous approaches are vital aspects of this process. To date, numerous computational tools for multiscale modeling have been developed and made available by various research groups. These tools may be offered as open-source or commercial software under different licensing agreements. Some of the prominent tools used for multiscale modeling include MATLAB [25], COMSOL Multiphysics [88], ANSYS [89], OpenFOAM [90], LAMMPS [91], Simvascular [92] etc.

Cancer therapy has made notable progress in recent years, but there are still numerous challenges that persist and require attention. Some of the key challenges are:

Drug resistance, cancer cells have the ability to mutate and become resistant to chemotherapy and other drugs, posing a significant challenge in treating advanced or recurrent cancers [93].

Side effects, many cancer treatments can lead to severe side effects, such as nausea, vomiting, hair loss, and fatigue [94].

Tumor heterogeneity, tumors often consist of various cell types with distinct genetic mutations. This heterogeneity complicates the development of treatments that effectively target all cancer cells within a tumor [95].

Late diagnosis, a significant number of cancers are diagnosed at an advanced stage, when treatment becomes more challenging [96]. This is often because early-stage cancers may not present noticeable symptoms.

Cost, cancer treatment can impose substantial financial burdens on patients and healthcare systems alike. The high costs associated with treatment can create barriers to access for some patients, limiting their ability to receive necessary care [97].

To address these limitations, continuous advancements in clinical experiments and treatment strategies are necessary. Integration of mathematical modeling approaches, such as optimizing drug distribution mechanism, cell death effect due to different treatment modalities and exploring combination therapies, holds great promise in improving treatment efficacy and reducing costs [98]. There are many options for cancer treatment, a summary of all different types of treatment described subsequently with sub-sections is given in Figure 2 below. In this section, our focus is not on reviewing numerous mathematical models developed for various treatment modalities, but rather on understanding how mathematical modeling strategies can be applied to different treatments based on underlying assumptions and experimental studies information. By illustrating the fundamental concepts behind developing these mathematical models for different cancer treatment, we aim to foster greater interest in advanced research and enhance the efficacy of therapy.

Surgery

Surgery has traditionally been the primary treatment option for most of the cancers. The goal of surgical intervention is to remove the tumor completely and potentially cure the patient. However, one of the major limitations of surgery is that the tumor must reach a certain size before it can be effectively removed [99]. Moreover, cancer cells have the ability to invade and metastasize to distant sites or sometime risk of cancer recurrence [100]. This highlights the need for additional treatment modalities, such as chemotherapy or radiation therapy, to target any remaining cancer cells that may not be surgically removed [101]. Kohandel et al. [102] explored treatment strategies for ovarian cancer by developing and analyzing a simple mathematical model that integrated tumor growth dynamics with the effects of chemotherapy and surgery. Various growth models and cell-kill hypotheses are considered, with surgery assumed to eliminate a fixed fraction of tumor cells instantaneously. The model representation is as follows:

Where, N tumor size,  pharmacokinetic (PK) and pharmacodynamics (PD) effect function (PK deals with how the body processes a drug, including its absorption, distribution, metabolism, and excretion). PD, on the other hand focuses on how drugs exert their effects on the body. It involves studying the relationship between the concentration of a drug at its site of action and the resulting pharmacological response),  surgery indicator function. Equation (21) represents the rate of tumor size change over time, with  describing tumor cell growth dynamics. The study explored a general form for  as shown in equation (22), wherein parameters β₁, β₂ , and α allow for modeling shifts between exponential, Gompertz, and logistic tumor growth. To simulate chemotherapy effects, the study considers a single non-cell-cycle-specific drug agent and various cell-kill models, presented in eq. (23). The log-kill hypothesis [103] assumes cell-kill is proportionate to tumor population, while the Norton and Simon (NS) [104] model links cell-kill to growth rate. Another model, the  model, proposes cell-kill is proportional to a saturable function resembling the Michaelis-Menton form [105]. Additionally, surgery effects are modeled as instantaneous, killing a fixed fraction  of tumor cells, with  related to the fraction of removed tumor cells. The outcomes of the study shed light on treatment sequencing strategies for ovarian cancer, providing insights into the effectiveness of combining chemotherapy and surgery. By considering different growth models and cell-kill hypotheses, the study elucidates how treatment outcomes may vary based on the timing and order of interventions. The numerical model simulations were compared and validated with the benchmark literature [106]. Reader may refer to [102], for other parameters definition and details.

Radiotherapy

Radiation therapy, also known as radiotherapy, is a prevalent cancer treatment modality that employs high-energy particles or waves to eliminate cancer cells [107]. By disrupting the DNA within cancer cells, radiation therapy impedes their ability to grow and multiply [108]. Prolonged or excessive exposure to low-dose ionizing radiation has the potential to cause genetic defects and can lead to serious carcinogenic effects [109]. Therefore, ensuring precise and targeted delivery of radiation to cancerous cells while minimizing damage to healthy tissues is a critical aspect of radiation therapy. Farayola et al. [110] presented a mathematical model that simulates the radiotherapy skin cancer treatment process, incorporating important radiobiological factors such as repair and repopulation of cells, given as follows:

Where,  and  are the population of normal and cancer cells respectively,  are the proliferation coefficients, carrying capacities, and competition coefficients of the normal and cancer cells respectively,   radiation strategy. Equations (24) and (25) describe the rates of change in healthy and cancer cell populations over time, respectively. In both equations, the first term on the right-hand side represents the logistic nature of their growth, while the second term reflects the decay resulting from competition among the cells. Additionally, the third term accounts for the administration of radiation, which removes a significant number of cancer cells and a smaller number of healthy cells from the system [110]. The simulation results aided to determine how big the tumor and healthy cells might be after treatment. The study might be used to guess how the number of cancer cells changes during treatment and what size they might be when treatment is done. The simulated results of the model were compared with the clinical results [111] and found a coincidence. Reader may go through [110], for other parameters definition and details.

Chemotherapy

Chemotherapy involves the administration of drugs that target and kill cancer cells [112]. Unlike radiation therapy, which focuses on specific areas, chemotherapy drugs circulate through the bloodstream and reach various parts of the body, with the goal of eliminating cancer cells regardless of their location [112]. Despite its successes, chemotherapy is not without its limitations. One of the main drawbacks is the toxicity of chemotherapy drugs, which can harm healthy cells along with cancerous ones. Recently, Miller et al. [69] utilized a multiscale modeling approach to evaluate the first-line chemotherapy response of individual patient tumors in intermediate and advanced non-small cell lung cancer. The study aims to classify patients' response to chemotherapy as disease-control, including complete-response, partial-response, and stable-disease categories. The outline of chemotherapy model in this given as:

Where,   drug concentration,  drug extravasation,  drug uptake, and  drug infusion. Equation (26) considers bolus administration, simulating drug transport with diffusivity Dc based on vascular extravasation within tumor tissue. The net rate  accounts for drug uptake by tumor and normal cells, as well as washout from interstitial space, representing drug half-life. In equation (27), drug extravasation by  is assumed constant from the vasculature, with  representing the transfer rate from pre-existing and new vessels,  indicating interstitial fluid pressure (IFP), and  denoting effective pressure (IFP at which net volume flux from vasculature is null). Initially, vascular drug concentration is , and extravasation follows first-order kinetics for constant drug infusion, as described in equation (28). The study findings suggest that employing diverse modeling techniques for cancer and analyzing the variation of chemicals within tumor tissue over time could enhance predictions of a patient's response to initial chemotherapy. The mathematical model and its simulations were validated and found consistent with the literature [113]. Readers may see [69], for details and other definitions.

Adjuvant and neoadjuvant chemotherapies: When chemotherapy is administered to patients before or after surgery or radiation therapy, it is referred to as adjuvant or neoadjuvant chemotherapy [114]. In neoadjuvant chemotherapy, the treatment involves administering chemotherapy drugs prior to surgery or radiation therapy [114]. The purpose of this approach is to reduce the size of the tumor and prevent its spread to other areas, a process known as metastasis. By shrinking the tumor, neoadjuvant chemotherapy aims to facilitate a more accurate and efficient surgical or radiation treatment. On the other hand, adjuvant chemotherapy is given after surgery and radiation therapy [114]. Its primary objective is to eradicate any remaining cancer cells that may be present after the initial treatment. It plays a crucial role in reducing the risk of cancer recurrence. By targeting residual cancer cells, it aims to prevent the return of the disease and improve long-term survival rates [114]. Gaffney [115] studied an adjuvant mathematical model for chemotherapy to find the best time to switch between drugs. This helps stop tumors from getting resistant to the treatment. They also thought about how the body absorbs drugs and when to give them. In the modeling framework, several key assumptions were made. Firstly, tumor cells experience homogeneous exponential growth without chemotherapy, especially during the subclinical phase [116]. This is justified by the relatively consistent doubling time of growing tumors at this stage. Secondly, the model considers four cell types based on drug resistance, with two types of drugs. Cell types include those sensitive to both drugs, resistant to one drug, or resistant to both drugs, all with identical proliferation rates. Lastly, the effects of both drugs on sensitive cells are combined linearly, reflecting the commonly observed log cell kill effect [103], where the surviving tumor cell fraction decreases linearly with increasing chemotherapy dose.

Following the assumption, the model equations are given in eq. (29) and (30), as follows:

Where, , are the tumor cell type according to sensitive to different chemo drug, , are cell kill function. The results of the study showed big improvements in how well treatments work when the time between switching drugs is made shorter in models that don't consider different cell growth stages. The mathematical model equations and their numerical simulations were compared and validated with their previous study [117]. One can go through [115], to understand the different parameter definitions and model in detail.

Hormone Therapy

Hormone therapy aims to disrupt the uncontrolled growth of cancer cells that rely on hormones for their growth [118]. Medications have the capability to inhibit hormone receptors on cancer cells, thereby impeding the attachment of hormones and their stimulation of cell growth [118]. Following initial treatment, hormone therapy may be employed to reduce the likelihood of cancer recurrence. Wang et al. [119] examined a mathematical model that looks at how the body's immune system responds to breast cancer, along with the effects of certain chemicals and hormones like estrogen. The following equations represent the model outline:

Where, CTLs, cytotoxic T cells, , are the populations of tumor cells, CTLs, and helper T cells. The proposed model consists of three time-dependent differential equations describing interactions among CTLs, helper T cells, and tumor cells under certain assumptions taken from previous studies [120,121]: a) tumor cells grow logistically without immune response, b) CTLs can kill tumor cells, c) tumor cells can activate naive and noncytotoxic cells, d) activated CTLs grow logistically, e) CTLs become inactive after interacting with tumor cells, and f) helper T cells can stimulate CTLs. Initially, tumor cell growth follows a logistic law with proliferation rate  and carrying capacity . After immune response stimulation, tumor cells interact with CTLs, leading to tumor cell death at rate . This process is represented by the differential equation (31).

Next, CTLs play a crucial role in the immune system's defense against tumor cells. The study models the logistic growth of CTLs with a growth rate  and a maximum carrying capacity . The interaction between tumor cells and CTLs is modeled by a hill function with a dissociation constant . After encountering tumor cells, CTLs eventually become inactive, represented by the term . Helper T cells further stimulate CTLs through cytokine, with the proliferation rate denoted by  and the natural death rate by . These processes are combined in the differential equation (32).

Lastly, Helper T cells are vital for supporting the immune response. The study models their growth with a constant influx represented by  and a hill function of tumor cells with a proliferation rate  and a dissociation constant  The natural death rate is denoted by . These factors are combined and represented in the differential equation (33). The model simulation results showed that the developed model was significant for studying how estrogen works and how it affected breast cancer and the immune system. They also found that hormone therapy, like estrogen-blocking drugs, could reduce the chances of getting breast cancer, especially in the early stages when the cancer was sensitive to hormones. The simulated results were compared with the clinical studies outcome [122] and found a good consistency.

Immunotherapy

Immunotherapy harnesses the power of the body's natural defense system to combat various diseases, including cancer, allergies, and infections [123]. It utilizes components such as antibodies, cytokines, and T-cells to enhance the immune response against specific targets [124]. However, it is important to acknowledge that immunotherapy has its limitations, and not all patients respond positively to this treatment approach. While it has recorded high remission rates in some cases, there are instances where patients do not experience significant benefits [125]. Castiglione and Piccoli [126], constructed a mathematical model to describe the interaction between cancer and the immune system and then applied the theory of optimal control to determine the most effective timing and dosage for stimulating the immune system with the therapeutic agent. The model focuses on several key immune cell populations, including CD4 T-helper cells and CD8 cytotoxic T cells, along with cancer cells. Dendritic cells, are the source of tumor-associated antigens (TAA), are introduced externally and play a crucial role in the system, given as follows:

Where H represents tumor-specific CD4 T helper cells, C denotes tumor-specific CD8 T cells or CTLs cytotoxic cells, M signifies cancer cells exposing TAA, D stands for mature dendritic cells loaded with TAA, and I represent IL-2 secreted by H, responsible for T cell growth. In Eq. (34), the term  represents the bone marrow's production of a few tumor-specific cells. Similarly, in Eq. (35), an equivalent term is provided for tumor-specific cytotoxic cells. The second term (RHS) in Eq. (34) signifies the clone expansion of tumor-specific helper cells following the presentation of tumor-associated antigens by dendritic cells. In Eq. (35), the third term (RHS), signifies the expansion of tumor-specific cytotoxic cells, induced by interactions with either tumor or dendritic cells. In Eq. (36), first term (RHS) denoted, cancer cells experience limited growth, and are eliminated by tumor-specific cytotoxic cells, indicated by second term. Dendritic cells in Eq. (37) are targeted and killed by cytotoxic cells, represented by . Interleukin IL-2, produced by helper cells upon recognizing tumor-loaded dendritic cells , is utilized by cytotoxic cells during clonal growth  (Eq. (38)), while  denotes the degradation of free interleukin. The simulation results observed that, high drug vaccination doses were given to reduce the size of the tumor. The numerical model was validated with the outcome of [127] and observed a good alignment.

Stem cell therapy

Stem cells possess a remarkable ability to self-renew and differentiate into various specialized cell types [128]. This unique characteristic of stem cells makes them a valuable tool in regenerative medicine, as they can be directed to attach themselves to specific cells and contribute to the repair and regeneration of damaged or diseased tissues [129]. In cancer therapy, high doses of radiation and chemotherapy can have detrimental effects on natural stem cells in the body [130]. Consequently, stem cell transplants have emerged as a useful approach for restoring and replenishing these damaged stem cells after treatment. Sfakianakis and colleagues [131] proposed a model aimed at understanding the dynamics of cancer therapy in the presence of regular cancer cells, cancer stem cells, and fibroblast cancer cells. The differential equations representing each cell type are nearly identical. Here, we will focus solely on describing the model equation for the time rate dynamics of cancer stem cells, given as:

Where, , are representing differentiated cancer cells, cancer stem cells, fibroblast cancer cells respectively,  extracellular matrix density, and the proliferation part function is defined as . In equation (39) and (40), the study initially uses a basic haptotaxis system in the model, where, the system involves cellular diffusion and logistic proliferation of cancer cells. Regular and cancer stem cells (CSCs), both interact with the ECM in a similar manner and connected through the epithelial-mesenchymal transition (EMT) (it’s a biological process where epithelial cells, which are typically tightly connected and organized in layers, undergo changes to become mesenchymal cells), where  denote the rate of EMT. The EMT rate depends on the amount of epidermal growth factor (EGF). Additionally, it is assumed that  may increase and reach a constant value asymptotically. The study also discussed the transition of cancer stem cells (CSCs) into cancer-associated fibroblast cells within the tumor environment. This transformation is represented by a term denoted as , where  represents the production coefficient for their transdifferentiation. The numerical model validation of this study was done by comparing the simulation results with the outcome of the literature [132] and found a good agreement. Readers may look into the article [131], for all other model equations and parameters detail.

Combination cancer therapy

Combination therapy, which involves the use of two or more different therapeutic approaches, has emerged as a cornerstone of future cancer treatment [133]. This approach has the potential to overcome drug resistance and enhance anti-cancer activity by reducing tumor growth, metastasis, and cancer stem cell populations, as well as inducing apoptosis [134]. Additionally, combination therapy offers cost efficiency, ensuring that advanced medical regimens are accessible to economically disadvantaged populations [135]. Conventional monotherapies often target highly proliferating cells without proper selectivity, resulting in the death of healthy cells along with cancer cells [136]. In combination therapy, the potential for toxicity is minimized as different agents target distinct pathways, allowing for lower quantities of each therapeutic agent and improving cost efficiency [137]. Frei et al. in 1972, demonstrated that using chemotherapy after surgery could substantially increase the likelihood of curing cancer [138]. Despite the advantages of combination cancer therapy, it has limitations. Unfavorable drug interactions and pharmacokinetics among the combined drugs can impact the expected therapeutic outcomes [139]. The complexity of drug combinations requires careful attention and control, as empirical experiments and clinical tools are necessary to analyze the effects of drug combinations properly [140]. Additionally, determining and predicting the ideal combination approaches to minimize complications faced by patients during combination cancer therapy can be challenging [141]. To address these challenges, mathematical models and computational simulations are utilized alongside clinical experiments to accelerate predictive analysis [142]. Some important combination therapy strategies include chemoimmunotherapy [143,144], chemoradiotherapy [145], chemovirotherapy [146,147], immunovirotherapy [131], radioimmunotherapy [148], and radiovirotherapy [149].

Chemoimmunotherapy: Chemoimmunotherapy represents a promising approach that harnesses the synergistic potential of tumor-suppressive drugs and the immune system's intrinsic anti-tumor response to achieve enhanced and durable therapeutic outcomes [150]. In the field of chemoimmunotherapy research, mathematical modeling has emerged as a valuable tool for studying the complexities of drug-immune system interactions. Pioneering work by Kirschner and Panetta laid the foundation for mathematical modeling of the immune system in the context of immunotherapy [127]. Building upon this foundation, recent advancements have been made, such as the mathematical model proposed by Ophir Nave [143] for chemoimmunotherapy in brain cancer. Nave's model encompasses a system of first-order nonlinear ordinary differential equations, providing explicit and analytical representations of the underlying biological processes. This study aimed to find the best treatment plan by trying different schedules and checking how stable they were. It found that a treatment schedule with a 7-day gap between sessions worked best. Rodrigues and colleagues [144] created an ordinary differential equation model to scrutinize the efficacy of chemoimmunotherapy in addressing chronic lymphocytic leukemia, a form of blood cancer. The model's focus was on exploring how both the immune system and chemoimmunotherapy contribute to combating the disease and potentially achieving remission. The model assumed cancer cells grew logistically, immune cells were naturally replenished by the host but also decay exponentially, and cancer cells stimulate new immune cell production until saturation. Interaction between cancer and immune cells has negatively impacted each other, with the rate proportional to their encounters. Chemotherapeutic drug dynamics followed the log-kill hypothesis with a Michaelis-Menten saturation response and first-order elimination kinetics. The outline of the model combining all these assumptions given is as follows:

Where, N,I number of cancer and immune cells respectively, Q the amount of chemotherapeutic agent. The equation (41) illustrates the conventional logistic cancer growth through its first term (RHS), while the  term gives the interaction between cancer and immune cells within each population , respectively). Additionally, the immune system is characterized as , as elaborated as in [151]. Equation (43) depicts the first-order pharmacokinetics of a chemotherapeutic drug originating from an external source. Regarding pharmacodynamics, the last term in equations (41) and (42) accounts the log-kill hypothesis, featuring a Michaelis–Menten drug saturation response. A similar functional response, exhibiting a positive sign, governs the production of immune cells stimulated by cancer, denoted as . The simulated outcomes indicated that chemoimmunotherapy plans might work well in treating Chronic Lymphocytic Leukemia, as long as chemotherapy doesn't hinder immunotherapy's effectiveness.

Chemovirotherapy: Virotherapy, an innovative and advanced approach in cancer treatment, harnesses the potential of viruses to selectively target and destroy cancer cells while sparing healthy cells [152]. However, its application as a standalone therapy remains limited, and many clinical trials have shown only modest success [152]. Salim et al. [146] investigated the combined impact of oncolytic virotherapy and chemotherapy on tumor cells using an ordinary differential equation (ODE) based mathematical model. The study's model was built on several key assumptions. Firstly, tumor cells were treated as a single entity, subdivided into uninfected and infected cancer cells. Secondly, an enhanced immune response was anticipated with larger virus burst sizes. Thirdly, conflict between uninfected tumor cells and the virus led to a surge in infected tumor cells. Fourthly, the tumor microenvironment experienced alteration with increased drug infusion, impacting its concentration in body tissue. Fifthly, tumor load growth followed a logistic pattern in the absence of treatment, reaching a carrying capacity of 10⁶ cells. Sixthly, the response of tumor cells to drugs, immune cell production, and virus infection was modeled using the Michaelis-Menten constant, which considers saturation effects. Lastly, the density of infected tumor cells escalated as viruses proliferate in uninfected cells and burst. The model, then formed based on the above assumptions, comprised a system of seven nonlinear differential equations. Here, for understanding we considered the model equation for variable , virus infected tumor cell density:

Where, I  stands for virus-infected tumor cell density, U  for uninfected tumor density, V for free virus particles,  for tumor-specific immune response,  for virus-specific immune response, and C for drug concentration. The model describes the dynamics of uninfected tumor cells becoming infected by the virus at a rate . The drug's effects on both uninfected and infected tumor cells are represented by the term , which signifies the lysis induced rates by the drug. The parameters  and  denote the Michaelis-Menten constants, indicating killing rates when the virus and drug are at half-maximal levels. Additionally, both virus-specific and tumor-specific immune responses contribute to killing infected tumor cells, represented by the terms  and , where  and  are the respective lysis rates. Infected tumor cells die at a rate δI. The model simulation results were validated and found to be in line with the study [153]. For other variable equations and details please follow [146].

Chemoradiotherapy: The combination of chemotherapeutic drugs and radiation therapy is a widely used approach in the treatment of advanced cancer [154]. In a comprehensive review by Zhang et al. [155], the relationship between microRNAs and various microenvironmental parameters of colorectal cancer, including tumor cell proliferation, invasion into neighboring tissues, cell death, and response to chemoradiotherapy, was explored. Bashkirtseva et al. [145] developed a nonlinear ordinary differential equation model to explore how altering the sequence of administering cytotoxic drugs and radiation affected cancer therapy. Model equation represented as follows:

Where, P population of tumor cells, the equation (45) describes the tumor growth in its first positive term, while negative two terms representing cell destruction by cytotoxic drugs and radiation, respectively. The logistic growth term in Eq. (45) depicts tumor growth with a rate   and carrying capacity K. Cell kill is modeled by the function , where b denotes the maximum fractional cell kill, c is the drug concentration, and ρ represents tumor cell sensitivity to drugs. Additionally, a Holling type II functional response  accounts for the Norton-Simon hypothesis, reducing cell kill as tumor cell population approaches carrying capacity. Radiotherapy sessions are represented by Dirac delta distributions , with  indicating the session index and  its delivery time. The model simulations were compared and validated with the clinical study [156]. See [145], for other parameters definition and details.

Radioimmunotherapy: The use of targeted antibodies combined with radionuclides, a therapy known as radioimmunotherapy. Traditionally, radiation beams have been used for the treatment of localized diseases, but radioimmunotherapy represents a more advanced treatment approach for both local and diffuse tumors [157]. Ashrafizadeh et al. [158] conducted a review exploring how the Abscopal effect, which demonstrates how radiation can activate the immune system against cancer cells, can be exploited in radioimmunotherapy. Gonzalez-Crespo et al. [148] explored a study using a bio-mathematical model to investigate how tumors react to radioimmunotherapy, a treatment that combines radiotherapy with immunotherapy. The study observed that using typical doses of radiation might not be the ideal choice. However, giving very large doses in a short time might also not work well. This could happen because these high doses might not work as well, or they might harm the blood vessels in the tumor too much.

Immunovirotherapy: Oncolytic Virotherapy, which involves the induction of the immune response through the use of viruses, plays a crucial role in combating cancer [159]. In this treatment approach, viruses are genetically engineered to enhance the production of inflammatory cytokines and stimulatory molecules, thereby activating immune cells to target cancer cell growth [160]. Engeland et al. [161] conducted a recent review on the current state of immunovirotherapy, with a focus on the challenges of implementation and the roles of mathematical models and experimental approaches in improving cancer treatment. Sfakianakis et al. [131] presented a mathematical modeling framework to study tumor virotherapy with mediated immunity. The model was designed to capture the complex interactions among tumor cells, oncolytic viruses, and the immune system in a dynamic environment with a moving boundary. Senekal and colleagues [162] explored how the recruitment of natural killer cells to the tumor microenvironment affects oncolytic virotherapy using a mathematical model.

Radiovirotherapy: Radiovirotherapy, which involves the use of radiation-mediated isotopes and viruses as carriers to target tumor cells, synergistically combines anti-tumor properties in therapy regimens. Touchefeu et al. [163] conducted an extensive review on the principles and biological characteristics of radiovirotherapy relevant to oncology. Al-Tuwairqi et al. [149] developed a set of ordinary differential equations to simulate the dynamics of aggressive tumor growth during radiovirotherapy treatment.

First, in the context of virotherapy, model incorporates three variables:  the density of uninfected cancer cells,  infected cancer cells, and  free virus particles. It is considered that uninfected cells grow exponentially at a rate  and become infected when encountering free virus particles at a rate β, while the virus replicates within infected cells, causing cell lysis at rate δ and generating new virion particles with burst size b. Free virus particles are eliminated from the body at a rate α. In radiotherapy, the model integrates the density of tumor cells irreversibly damaged by radiation , assuming radiation  causes damage to both tumor cell types at rates α₁ and α₂ respectively, with damaged cells subsequently removed at rate . These elements are captured in a system of nonlinear ordinary differential equations (ODEs) outlined in equations (46-49). The study findings suggest that radiovirotherapy proves to be more effective than virotherapy, offering a promising alternative by completely eliminating tumors. Numerical simulations were conducted, validated against analytical results of [163] and demonstrated good agreement.

Targeted therapy

Targeted therapy utilizes specific drugs or monoclonal antibodies that selectively target and inhibit the proteins involved in cancer cell growth and metastasis [164]. In targeted therapy, various mechanisms are employed to disrupt the tumor microenvironment and impede cancer progression. One approach involves suppressing angiogenesis [165]. Abazari and colleagues [166] developed a computational framework to assess the effectiveness of targeted drug delivery to cancer cells using antibodies against oncogenic cell-surface receptors. With the help of mathematical modeling of angiogenesis, they studied three different stages of tumor progression [166]. The results showed that lowering the number of blood vessels in the tumor decreased the interstitial fluid pressure, making it easier for drugs to get into the tumor. With regular chemotherapy, the drug Doxorubicin gathered unevenly around the blood vessels. However, with the new treatment using immuno-liposomes, along with releasing drugs inside cells and keeping them in tumors for longer, the drugs spread more evenly and killed more cancer cells than regular chemotherapy alone. It's also noted that intracellular Doxorubicin tends to be deposited heterogeneously primarily in the perivascular areas due to its limited penetration into tumor tissue [166]. Moreover, compared to conventional chemotherapy, immunochemotherapy significantly enhances the intracellular Dox concentration [166]. Rezaeian et al. [167] conducted a study to examine the effectiveness of magnetically controlled drug targeting (MCDT) in intraperitoneal (IP) chemotherapy for treating peritoneal malignancies (PMs). They developed a mathematical model to compare the efficacy of drug delivery using magnetic nanoparticles (MNPs) with conventional IP chemotherapy. In this study, we are only focusing on those model equations from [167], which illustrate how magnetic force can drive targeted nanodrug transport by exploring the magnetically controlled mass transport equation. The investigation involved computing the free drug concentration in the interstitial fluid using the convection-diffusion equation, shown in equation (50). Here,  stands for the free drug concentration, while  and φ represent the concentration of cell surface receptors and the tumor volume fraction available for the drugs, respectively. Effective diffusion coefficient given by . Moreover,  denotes the constants for the drug association rate with cancer cells. We also examined the velocity  in tissue, it was determined by summing velocity of interstitial fluid  with the equilibrium velocity . This equilibrium velocity occurs when the Stokes drag force  equals the magnetic force , where,  stands for dynamic viscosity of the interstitial fluid [167], are presented in equation (51) and (52).

In the equation (50),  represents the net rate of free drug exchange from blood and lymphatic vessels, while  stands for the drug concentration bound to the cancer cell. In eq. (52), H denotes the magnetic field,  represents the volume of MNPs and  represents the magnetic susceptibility of MNPs. Reader may follow [167] for complete model equations, parameters definition and details. Outcome of the study observed that, during a 60-minute treatment, MCDT notably enhances both the depth and area of drug penetration compared to conventional IP chemotherapy. Moreover, MCDT results in a higher fraction of killed cells within the tumor, suggesting its potential to enhance treatment effectiveness. The numerical model [167] for concentration distribution was validated with the benchmark study by Au et al. [168].

Advancements in understanding cancer evolution and personalized treatment offer opportunities to manage cancer as a chronic condition, a concept increasingly explored in oncology. Model-based approaches integrating tumor growth and resistance evolution hold promise in achieving this goal. Mathematical models, drawing on clinical trial data, facilitate quantitative understanding of tumor growth dynamics, drug PK-PD profiles, and resistance development in cancer patients. The adoption of various mathematical models in cancer research has grown significantly in recent decades. This review showcases the development of tumor growth models ranging from basic differential equation model to advanced models including agent-based, stochastic, and multiscale model, rooted in biophysical phenomena and experimental data, demonstrating their utility in simulating complex biological processes and interactions. Feasible model structures used to describe and forecast tumor dynamics and resistance evolution in solid tumor treatments are discussed herein.

In tumor size modeling, different types of equations like ODEs, algebraic equations, and PDEs are commonly used. Notably, ODE models primarily examine temporal dynamics, analyzing time-resolved data across cellular to tissue scales effectively [34]. Conversely, PDE models focus on spatial dynamics, depicting tumor growth's spatial spread into surrounding healthy tissues. They consider factors like oxygen and nutrient levels, vascular angiogenesis, tumor angiogenic factor concentrations, drug distributions, and haptotaxis effects [34]. The choice of model varies because it's hard to accurately capture the long-term natural growth of tumors in patients [16]. While some models set limits on tumor growth, others, like exponential growth models, don't. Hybrid modeling techniques give comprehensive understanding of biological processes by combining both discrete and continuous approach. In the work by Jeon et al. [43], they added the concentration of oxygen as another nutrient component in the continuous part of the model. Meanwhile, the discrete part of the model describes individual cell behaviors. On-lattice models use a fixed grid system, which makes simulations and computations simpler. However, this method restricts cell movement and interactions to nearby grid sites, limiting their ability to fully represent real tumors. On the other hand, off-lattice models allow cells to move more realistically in a continuous space, reflecting actual tumor shapes and movements better. Stochastic modeling incorporates randomness into models, providing a more realistic view of tumor behavior. Multi-scale models can consider multiple spatial and temporal scales simultaneously. Additionally, multi-scale models maintain the ability to accurately capture biological mechanisms in quantitative terms.

Using mathematical modeling approaches, such as optimizing how drugs spread, understanding how different treatments kill cancer cells, and exploring combination therapies, can greatly improve treatment effectiveness. For example, the review discussed various treatment effects based on important cell-killing theories like log-kill, Norton and Simon, and Emax. By examining these theories, the study showed how treatment results might change depending on when and how treatments are given. For instance, in radiotherapy, controlled radiation doses can remove many cancer cells and a few healthy cells from the body. Additionally, studying drug concentration in chemotherapy using mathematical models helps analyze how chemicals distribute in tumor tissue over time, improving predictions of how patients will respond to treatment. Neoadjuvant chemotherapy aims to shrink tumors before surgery or radiation therapy, making these treatments more accurate and efficient. Adjuvant chemotherapy, on the other hand, is given after surgery and radiation to eliminate any remaining cancer cells, reducing the chances of cancer returning.

Originally, models used deterministic ordinary differential equations (ODEs), but they couldn't fully capture the randomness in biological systems [59]. However, a study has shown that ODE models can still give comparable results to those derived from stochastic models [169]. But, with improved computing power, stochastic models are becoming more popular because they include randomness, making them better at depicting tumor behavior realistically [59]. Moreover, to make these models accurate, we need precise parameter estimates from lab experiments or clinical trials [85]. But there's no single model that covers all the processes in a real tumor [16]. Instead of aiming for full complexity, it might be better to enhance existing models gradually, making them more flexible and capable of providing detailed information on specific aspects. Therapeutic treatments also need deeper exploration, including how drugs move within extracellular and intracellular space, and how they affect the tumor's overall mechanics by causing cell death. In summary, the biggest challenge is combining the various complexities associated with tumor growth and treatments that occur at different scales into one mathematical model.

It is crucial to note that the mathematical equations employed in the models discussed here are not exhaustively examined. Rather, the objective has been to provide a comprehensive overview of mathematical models elucidating tumor dynamics and treatment strategies, as these models hold significant value in guiding medically and economically optimal approaches to enhance the quality of life, adhering to the core objective of “cancer treatment for all”.

The field of tumor growth dynamic modeling has seen remarkable advancements over the past few decades, with further opportunities on the horizon. While existing models have evolved significantly [170,171], transitioning from classical functional forms to more complex joint models, challenges remain in advancing personalized tumor growth dynamic models that account for individual-level tumor heterogeneity and intra subject variability. Integrating tumor physiology [172], cellular and molecular heterogeneity, and genomic alterations [173] into models could deepen our understanding of the interplay between tumor growth and clinical factors.

Tumor growth dynamic modeling also holds a substantial promise in oncology drug development [174], from preclinical evaluations to late-stage clinical trials. Incorporating real-time tumor kinetics with prognostic outcomes could enable earlier identification of clinical benefits, guiding timely decisions on treatment optimization [175]. Additionally, joint modeling approaches that have been discussed here, may provide more accurate assessments of patient outcomes and facilitate the early advancement or termination of investigational treatments, improving both patient care and resource allocation. Overcoming current limitations, such as the limited application of individual tumor lesion dynamics, will further enhance the predictive power and clinical utility of TGD models.

In addition to summarizing existing mathematical models in tumor growth and treatment, this review suggests new directions for research. Understanding how CTLs, helper T cells, and tumor cells interact dynamically can improve hormone therapy, immunotherapy, and combination therapies. Exploring the roles of DNA and microRNAs can advance current tumor dynamics models. Furthermore, studying and quantifying the negative impacts of high doses of radiation and chemotherapy on natural stem cells through mathematical modeling can help optimize stem cell therapy. By highlighting these emerging areas, the review encourages further exploration and inspires researchers to delve into new aspects of cancer treatment optimization.

The paper provides a comprehensive review of mathematical models used to study tumor growth dynamics and treatment strategies in cancer research. It emphasizes the need for interdisciplinary collaboration to address the complexities of cancer. Mathematical models play a crucial role in elucidating tumor progression, resistance mechanisms, and microenvironmental influences, offering insights into treatment outcomes. The review highlights advancements in modeling, demonstrating their value in predicting tumor behavior and guiding therapeutic strategies. With expanding computational resources and clinical data, these models are increasingly capable of enabling personalized cancer treatments. While challenges remain in integrating all relevant biological and clinical factors, ongoing advancements will enhance treatment optimization and improve patient outcomes. In conclusion, this review contributes to the growing body of knowledge in mathematical oncology and highlights the vital role of mathematical modeling in comprehending tumor dynamics and guiding treatment decisions.

We sincerely thank the reviewers and editor for their valuable comments and suggestions, which greatly improved the quality of our manuscript. We also extend our gratitude to Motilal Nehru National Institute of Technology (MNNIT) Allahabad for providing the necessary support and resources for this study.

No funding is available.

The authors declare that they have no conflict of interest.

All data generated or analyzed during this study are included in the article.

All authors made equal and significant contributions to this study and approved the final version of the manuscript for publication.