Mathematical modeling in radiotherapy for cancer: a comprehensive narrative review

Abstract

Mathematical modeling has long been a cornerstone of radiotherapy for cancer, guiding treatment prescription, planning, and delivery through versatile applications. As we enter the era of medical big data, where the integration of molecular, imaging, and clinical data at both the tumor and patient levels could promise more precise and personalized cancer treatment, the role of mathematical modeling has become even more critical. This comprehensive narrative review aims to summarize the main applications of mathematical modeling in radiotherapy, bridging the gap between classical models and the latest advancements. The review covers a wide range of applications, including radiobiology, clinical workflows, stereotactic radiosurgery/stereotactic body radiotherapy (SRS/SBRT), spatially fractionated radiotherapy (SFRT), FLASH radiotherapy (FLASH-RT), immune-radiotherapy, and the emerging concept of radiotherapy digital twins. Each of these areas is explored in depth, with a particular focus on how newer trends and innovations are shaping the future of radiation cancer treatment. By examining these diverse applications, this review provides a comprehensive overview of the current state of mathematical modeling in radiotherapy. It also highlights the growing importance of these models in the context of personalized medicine and multi-scale, multi-modal data integration, offering insights into how they can be leveraged to enhance treatment precision and patient outcomes. As radiotherapy continues to evolve, the insights gained from this review will help guide future research and clinical practice, ensuring that mathematical modeling continues to propel innovations in radiation cancer treatment.

Introduction

Radiotherapy is a crucial treatment modality in the fight against cancer, used for over half of all cancer patients [1]. The application of mathematical modeling has been essential in the progression and advancement of radiotherapy over time. Figure 1 depicts some key examples of applying mathematical models in various aspects of radiotherapy. For cancer biology, mathematical modeling is used to help to understand the mechanisms of tumor growth, invasion, angiogenesis, and metastasis and predict treatment response. For radiobiology, modeling has been used to quantify biologically effective dose (BED), primarily through the use of the linear-quadratic (LQ) model, and has been used to assess risks (i.e., normal tissue injury) and benefits (i.e., tumor control) of radiotherapy [2, 3]. These are extensively utilized in the planning and assessment of radiotherapy.

Fig. 1

Applications of mathematical modeling in radiation oncology. A list of abbreviations is provided for reference at the end of the article

The LQ model provides a fundamental framework for comprehending the impact of radiation on living cells, allowing for the prediction and evaluation of the reaction of both tumor and normal tissues to different radiation dosages [2]. The BED (generally employing LQ conversions) serves a quantitative assessment for comparing various radiation treatments, accounting for factors such as the total dose, fractionation, and length of treatment [3].

Moreover, mathematical modeling is essential for the development and applications of tumor control probability (TCP) and normal tissue complication probability (NTCP) models [4, 5]. These models are crucial in optimizing treatment strategies by achieving an optimal tradeoff between maximizing tumor eradication and limiting harm to nearby healthy tissues.

Mathematical models are also extensively used throughout the clinical workflow of radiotherapy treatment planning. Radiation dose calculation has long depended on algorithms of various levels of complexity and stochastic beam transport simulation methods like convolution/superposition algorithms and Monte Carlo-based simulations [6, 7]. Inverse treatment plan optimization, a mainstay of modern intensity-modulated radiation therapy (IMRT), utilizes dose-informed linear and nonlinear objective functions, genetic and population level modeling methods, gradient-based and stochastic search algorithms such as simulated annealing methods, etc., allowing efficient radiotherapy dose optimization maximizing dose to targets while minimizing dose exposures to normal tissues [8, 9]. Techniques in computational mathematics such as finite element method and B-spline approximation have been the key methods that facilitate deformable image registration, which is critical when using various diagnostic imaging modalities that need to be fused to simulation scans (generally CT-based but also MRI-based with some commercial systems) [10,11,12].

With stereotactic body radiation therapy (SBRT) and stereotactic radiosurgery (SRS), mathematical modeling has increasingly examined the intricate connections between dose-fractionation and treatment effectiveness in modern radiotherapy methods [13]. Moreover, the advent of ultra-high dose rate FLASH radiation, along with the growing popularity and sophistication of spatially-fractionated radiation therapy (SFRT) have created new avenues for research in the field, with mathematical models playing a vital role in elucidating the mechanisms, examining impacts, and optimizing these treatment paradigms [14, 15].

One of the potentially most exciting applications of mathematical modeling in radiotherapy is the investigation of the intricate interaction between radiotherapy and the body's immune response, particularly in the context of immunotherapy [16, 17]. This synergy and the intricate considerations required in combining and timing radiotherapy with chemotherapy and surgery create a complex and exciting terrain in which mathematical models can provide essential insight. Another equally intriguing area of research is the integration of multi-scale, multi-omics data using a systems approach to better understand and combat cancer [18,19,20].

To guide readers through the diverse applications of mathematical modeling in radiotherapy, this review is organized into distinct sections, each addressing a specific domain of mathematical modeling and its applications. Following a technical overview in "Mathematical modeling approaches in radiation oncology," the review explores various domains including "Cancer biology and classical radiobiology models," "Mathematical approaches and models for modern radiotherapy treatment workflow," "Mathematical approaches and models for temporal and spatial fractionation," "FLASH radiotherapy," and "Integrating radiotherapy with immunotherapy and immune-oncology." The review concludes with sections on "Mathematical modeling and digital twins," as well as "Challenges, opportunities, and outlook," and "Conclusion”. This structure aims to provide a clear progression from foundational concepts to cutting-edge innovations, facilitating a comprehensive understanding of the field.

This review aims to provide a comprehensive overview of the diverse applications of mathematical modeling in radiotherapy, placing particular emphasis on cutting-edge applications and emerging methodologies. It seeks to elucidate the current state and future potential of radiotherapy, highlighting the critical role of mathematical modeling in driving innovations and enhancing patient outcomes in radiation cancer treatment.

Mathematical modeling approaches in radiation oncology

Mathematical modeling approaches used in radiation oncology can range from simple arithmetic formulas such as the LQ model to differential equations, and machine learning methods. The models can be derived from first principles, based on fundamental scientific laws, empirical fitting, using data-driven approaches to match observational outcomes, or a hybrid approach that combines two or more methods. Often incorporating physical laws, differential equations are a powerful mathematical modeling tool in radiation oncology to effectively describe spatial–temporal dynamics, with frequently used variables of time, space, or size. They can be classified into ordinary differential equations (ODE), partial differential equations (PDE), and integro-differential equations (IDE). ODEs are equations involving functions and their derivatives with respect to only one variable, for example, describing cell number changes over a single variable, time, in cancer cell growth models and their response to radiation. PDEs are equations of functions with partial derivatives of one or multiple variables, for example, describing interacting variables and spatial dimensions in the oxygen-depletion and reoxygenation process following FLASH radiation. While every ODE is technically a PDE with only one independent variable, PDEs are generally used to model phenomena involving two or more independent variables, such as time and space. Interacting variables are frequently modeled with systems of ODEs or PDEs. The variables of PDEs can, but don’t necessarily have to, stay separated. IDE involves both integrals and derivatives of the solution function, combining both differential and integral calculus, for example, calculating radiation dose from modeling and integrating the spatial dose transport through tissue. Other differential equations used for mathematical modeling are stochastic differential equations, which also include terms that describe a stochastic process. SDE are often applied to model tumor responses to radiotherapy or other radiobiological processes considering both spatial dose distribution and temporal dynamics. Mathematical models can also be classified into deterministic, with predictable outcomes based on fixed rules, or stochastic, which incorporates randomness and variability. One widely-known example is the two types of normal tissue radiation effects, where conditions like cataract formation are deterministic with a dose threshold, and conditions like radiation-caused malignancies are stochastic without a dose threshold [21]. In addition to continuous differential equations, the models can be discrete, dealing with distinct and separate values, where each cell is considered as a discrete element and a set of rules regulating cellular activities and interactions, also known as cell-based models, agent-based models, or individual-based models, allowing exploration of single-cell responses to treatment that lead to observable dynamics at the tumor scale. Table 1 lists these main categories of mathematical modeling approaches and example applications in radiation oncology. In each application, hybrid approaches may also be used combining the strengths of different methods. Besides conventional mathematical modeling approaches, machine learning is also an increasingly popular approach to leverage classical models with large datasets for pattern recognition and making predictions [22]. Mathematics is at the core of machine learning methods, as those methods rely on mathematical algorithms to identify patterns and relationships in data, usually optimize a mathematical objective, and leverage mathematical methods such as linear algebra, calculus, and probability.

Table 1 Different types of mathematical modeling approaches and example applications in radiation oncology

Cancer biology and classical radiobiology models

Mathematical modeling has been employed to offer valuable insights into cancer biology, especially the complex dynamics of tumor growth, invasion, and metastasis. Over the past century, these models have evolved from simple theoretical frameworks to highly sophisticated approaches, reflecting advances in tumor biology and computational capabilities. To describe how cancer cells proliferate over time, models from the simple exponential growth model to more complex models accounting for other biological factors have been developed. The exponential model, introduced in the early twentieth century, assumes the tumor grows at a constant rate proportional to its current size and has no asymptotical value, i.e., no carrying capacity. In reality, factors like nutrient availability and accessibility and cell death inhibit infinite growth. Two of the most frequently used sigmoidal models, the logistic growth model and the Gompertz model, incorporate these growth-limiting factors, such as applying a carrying capacity factor [33, 34]. Such models are also often employed to simulate the interplay between tumor and normal cells [35, 36].

In the mid-twentieth century, spatial–temporal modeling approaches emerged to better understand tumor invasion and spread into surrounding tissues. To model how cancer cells spread from the primary tumor into surrounding healthy tissues, PDEs are often used to capture the spatial and temporal dynamics of cell movement, degradation of the extracellular matrix (ECM), and interactions with the microenvironment [37]. Among many examples, a hybrid model is developed to simulate cancer progression, within this framework. Reaction–diffusion models are used to model nutrient transport, cell proliferation, and cell-ECM mechanical interactions are modeled using PDEs as a viscous fluid through a porous medium, a combined continuum-discrete model is used for angiogenesis, the level set method is used to capture the topological changes like tumor splitting, and finite-element models are applied to simulate vascular endothelial growth factor (VEGF)-receptor interactions and vessel network formation [38]. The even more complex bioprocess of cancer metastasis has leveraged multi-scale mathematical modeling that encompasses the cellular, tissue, and systemic levels. In the early twenty-first century, multi-scale mathematical models were developed to address more complex bioprocesses such as cancer metastasis, encompassing the cellular, tissue, and systemic levels. These frameworks of tumor biology are further integrated with the effects of radiation for modeling the complex system of radiobiology and radiotherapy applications, as illustrated in Fig. 2.

Fig. 2

Radiation-cell interactions and mathematical modeling. Blue texts denote radiotherapy techniques, green texts denote mathematical models, and red texts denote biological elements and states. A list of abbreviations is provided for reference at the end of the article

Radiation-based cell killing is the basis for radiotherapy, which is triggered largely by radiation-induced double-strand DNA breaks. Radiation effects are classically described with 5 R’s of radiobiology: “Radiosensitivity” that different cells exhibit a wide range of sensitivity to radiation; “Repair” of radiation-induced DNA damage which may lead to cell cycle arrest, recovery from sublethal injuries, or death; “Reoxygenation” in which vascular changes within the treated tumor leads to increased blood flow that reverses hypoxia, increasing radiosensitivity; “Reassortment/Redistribution” in which the differential survival of cells in cell cycle phases of varying radiosensitivity get redistributed in the cell cycle; “Repopulation” with which survived tumor and healthy cells regenerate between dose fractions [39, 40]. Mathematical modeling has played a central role in describing the relationship between radiation dose and cell kill/survival and in simplifying the intricate interplays of the 5 R’s. The original 4Rs, Repair, Reassortment/Redistribution, Repopulation, and Reoxygenation, were introduced in the mid-twentieth century to describe key biological mechanisms that influence tissue response to fractionated radiotherapy, serving as a foundational framework for optimizing radiotherapy efficacy. The fifth R, Radiosensitivity, was introduced later in the 1980s based on experimental findings that highlighted significant differences in the inherent radiosensitivity of tumor cell lines, further refining our understanding of this intrinsic factor which plays a critical role in determining treatment outcomes (Fig. 3).

Fig. 3

Complex interplays between radiation and the immune system

Although some alternative models have been proposed, the prevailing framework in describing cell survival in response to radiation remains the LQ model (Eq. 1). This model characterizes the cell survival fraction after receiving a certain dose by an exponential decay with a linear term, α, and a quadratic term, β. The linear term represents the initial damage caused by radiation that is directly proportional to the dose, and the quadratic term represents double-strand DNA breaks caused by combining two sublethal damage events where the damage is proportional to the square of the dose. The LQ model was developed in the mid-twentieth century to provide a mathematical framework for understanding the dose-dependent effects of radiation on cell survival. It remains widely used due to its ability to capture the effects of dose fractionation and its simplicity, which allows empirical fitting to experimental and clinical data. The $\alpha /\beta$ ratio is determined through experimental and clinical studies by fitting survival curves of irradiated tissues or cells to the LQ model. For example, in vitro studies use cultured cell lines irradiated at different dose fractions to analyze the survival-dose relationship and characterize $\alpha \;{\text{and}}\;\beta$ parameters, while in vivo or clinical studies analyze animal or patient radiation outcomes to refine $\alpha \;{\text{and}}\;\beta$ estimates for specific tumors and normal tissues. As the classical radiobiology model, the LQ model serves as the mathematical foundation of radiotherapy, and the empirically fitted $\alpha /\beta$ ratio is widely used to describe the radiosensitivity of different cells or tissues. The LQ model also effectively illustrates the effects of dose fractionation and how it impacts tumor control and normal tissue toxicity [41]. Based on the LQ model and the $\alpha /\beta$ ratio, the BED calculation allows the comparison of varying fractionation schemes by normalizing their biological effects (Eq. 2). On the other hand, the LQ model does not consider some of the other R’s in radiobiology such as reoxygenation and redistribution, and is simplistic in the modeling of damages and repairs. While the LQ has been used across various dose-fractionation schedules, with very high fractional disease (such as over 8–10 Gy/fraction for SRS or SBRT), or for higher dose rates, the radiobiology is different than with conventionally fractionated radiation, and the utility of the LQ model has been questioned and debated [42,43,44]. At higher fractional doses or higher dose rates, cellular DNA repair pathways may become overwhelmed, leading to increased DNA damage accumulation, and the tumor microenvironment may get modulated, potentially altering enhancing immune system activation. High dose fractions and high dose rates therefore can lead to increased direct cell death and may induce vascular damage within tumors. Modifications on the LQ models or alternative models have been developed for these cases which will be discussed in further detail in the next section. Table 2 summarizes a list of equations that are discussed in this section for easy reference.

Table 2 A list of the discussed equations/radiobiology models for easy reference

LQ model

$$S\left( D \right) = e^{{ - \left( {\alpha D + \beta D^{2} } \right)}}$$

(1)

where S(D) is the cell surviving fraction after irradiated with a dose D, α is the linear component of cell killing per unit dose, and β is the quadratic component of cell killing per unit dose squared.

BED

$${\text{BED}} = nd \left. {\left( {1 + \frac{d}{{\frac{\alpha }{\beta }}}} \right.} \right)$$

(2)

where n is the number of fractions, and d is the dose per fraction.

While the LQ model functions at the cell level. TCP and NTCP are employed at the tumor and tissue level to quantify the probabilities of tumor control and normal tissue toxicity based on the received dose. These models are usually empirically fitted from clinical and experimental data, often combined with some theoretical components. For TCP (Eqs. 3a and 3b), the logistic (sigmoidal) model and Poisson model are among the most widely utilized approaches in research and sometimes made available in commercial treatment planning and inverse optimization software [45, 46]. For NTCP (Eqs. 4a and 4b), prevailing models include the Lyman-Kutcher-Burman (LKB) model, the critical volume (CV) model, etc. [47, 48]. Some of these NTCP models have also been made available in commercial treatment planning systems for biological plan optimization. These models are also often used to examine dose toxicity relationships and have yielded normal tissue dose constraints widely adopted in clinical practice [49, 50]. With concerns raised about model limitations on non-uniform dose distributions, models continue to be enhanced to address involved complexities and large pooled data are used to fit these models [47, 49, 51]. To account for the spatial distribution of dose on normal tissue, the equivalent uniform dose (EUD) model (Eq. 5) is often used [52].

Logistic TCP model

$${\text{TCP}}\left( D \right) = \frac{1}{{1 + e^{{ - \frac{{D - D_{50} }}{k}}} }}$$

(3a)

where D is the delivered radiation dose, D₅₀ is dose for a 50% probability of tumor control, and k is the slope of the sigmoidal curve at D_50.

Poisson TCP model

$${\text{TCP }}\left( D \right) = e^{ - N \times S\left( D \right)}$$

(3b)

where N is the initial number of clonogenic cells and S(D) is cell survival fraction after dose D.

LKB NTCP model (simplified)

$${\text{NTCP}} = \frac{1}{{\sqrt {2\pi } }}\mathop \smallint \limits_{ - \infty }^{t} e^{{ - \frac{{u^{2} }}{2}}} du$$

(4a)

where t is the normalized dose difference and u is a dummy variable used in the integration to calculate the cumulative normal distribution.

CV NTCP model

$${\text{NTCP}} = \left\{ {\begin{array}{*{20}l} {0,} \hfill & {v < v_{crit} } \hfill \\ {1,} \hfill & {v \ge v_{crit} } \hfill \\ \end{array} } \right.$$

(4b)

where v is the fractional volume of the organ receiving a dose above the threshold dose, and v_crit is the critical volume.

EUD model

$${\text{EUD}} = \left( {\mathop \sum \limits_{i} v_{i} D_{i}^{a} } \right)^{\frac{1}{a}}$$

(5)

where a is a tissue-specific parameter for dose–response relationship, and v_i is the fractional volume receiving dose D_i.

In addition to the classical LQ model, other physical or mathematical models were also proposed, especially in earlier times, to describe the radiation effects on cells. Some notable ones include the target theory models such as the single-hit model that emerged in the 1920s (Eq. 6), the multi-target, single-hit model that emerged in the 1950s (Eq. 7), and the lethal-potentially lethal model that emerged in the 1970s (Eq. 8) etc. [53,54,55,56]. The single-hit model was the earliest attempt to describe cell survival to radiation, focusing on the probabilistic nature of radiation interactions. The multi-target, single-hit model was later introduced to account for cell repairs and describe the “shoulder” observed in the cell survival curve. The LPL was a further development later that differentiate between immediately lethal damage and damage that could become lethal if not repaired, highlighting the dynamic nature of post-irradiation DNA repairs.

Single hit model

$$S\left( D \right) = e^{ - kD}$$

(6)

where S(D) is the survival fraction given dose D, and k is a constant representing the sensitivity of the target to radiation.

Multi-target, single-hit model

$$S\left( D \right) = 1 - \left( {1 - e^{{ - \frac{D}{{D_{0} }}}} } \right)^{n}$$

(7)

where D is the dose, D₀ is the dose at which there is an average of one lethal hit per target, and n is the number of targets.

Lethal-potentially lethal model

$$S\left( D \right) = e^{{\left\{ { - \left( {\alpha D + \beta D^{2} + \gamma D^{3} } \right)} \right\}}}$$

(8)

where γ represents higher-order terms accounting for potentially lethal damage repair dynamics.

Radiation effects on cell survival are also affected by other factors. For example, oxygen enhances the radiation damage due to free radical formation and fixation of radiation damage not allowing it for easy repair. Therefore, the radiation dose required to cause the same amount of cell kill for cells in a hypoxic condition is higher than in a well-oxygenated condition. Oxygen-enhancement ratio (OER), the ratio of the radiation dose required to achieve the same cell killing in anoxic conditions to the dose required in oxygenated conditions, is used to describe this effect (Eq. 9) [57]. Different cell types have different OERs, ranging around 2.5 to 3 for mammalian cells [58], with lower OER values for radiation with high linear energy transfer (LET) or at low doses [57, 59].

OER

$${\text{OER}} = \frac{Dose\;in\;hypoxic\;conditions}{{Dose\;in\;oxygenated\;conditions}}$$

(9)

Similarly, relative biological effectiveness (RBE) is used to describe the biological effectiveness of different types of ionizing radiation relative to reference radiation modality (x-rays) (Eq. 10) [60]. This biological effectiveness depends on the LET of the radiation [61]. Higher LET radiation such as proton and other charged particles has a higher RBE than low LET radiation such as x-rays.

RBE

$${\text{RBE}} = \frac{{Dose \;of\;standard \left( {x - ray} \right)\;radiation\;needed\;to\;produce\;the\;same\;biologic\;effect}}{Dose\;of\;a\;second\;radiation\;needed\;to\;produce\;the\;same\;biologic\;effect}$$

(10)

OER addresses the radiobiological impact of oxygen conditions and is often used for hypoxia and FLASH RT modeling. RBE describes the biological effectiveness of different radiation modalities or types relative to standard x-rays, and is often used in particle therapy modeling. Both these models provide insights into biological responses to radiation under specific conditions and complement other models like TCP, NTCP, and EUD by addressing unique aspects of radiobiological interactions.

Mathematical approaches and models for modern radiotherapy treatment workflow

Beyond classical radiobiology, mathematics plays a pivotal role in various aspects of modern radiotherapy. In this section, we briefly outline notable and commonly used mathematical approaches within the general radiotherapy treatment workflow. In subsequent sections, we delve deeper into important and emerging radiotherapy applications that leverage advanced mathematical modeling techniques. The mathematical models referenced in this section primarily pertain to computational imaging and radiotherapy processes.

Modern radiotherapy utilizes individual patients’ volumetric images to develop customized treatment plans. Mathematical techniques and computational tools have long been implemented in facilitating the radiotherapy treatment workflow. For radiation dose computation, the Monte Carlo simulation is a fundamental stochastic method for beam transport simulation and remains the gold standard [62]. Other dose calculation algorithms with varying levels of accuracy and simplifications have also been adopted in clinical treatment planning systems, such as pencil beam, convolution superposition, collapsed cone convolution, and analytical anisotropic algorithm [63,64,65,66,67]. For inverse plan optimization, methods like linear programming, multi-criteria Pareto optimization, greedy gradient search, and simulated annealing are widely applied, and machine learning methods are employed to devise processes that automate treatment planning with greatly reduced reliance on human interactions [68,69,70,71,72,73,74,75].

Modern radiation oncology relies upon 3-dimensional volumetric images such as computed tomography (CT), magnetic resonance imaging (MRI), and positron emission tomography (PET). Mathematical methods, such as the Fourier transform and various mathematical filtering algorithms, are fundamental to image generation, processing, and enhancement [76,77,78,79]. Target (the tumor or area intended for radiation treatment) and critical organ (healthy tissues near the target that must be protected to minimize radiation-induced damage) delineation is an important step of radiotherapy treatment planning, and multi-model images such as CT, MRI, and PET are usually used to optimize the visualization of these structures. Mathematical approaches are an integral part of image registration and involve modeling tissue deformations to ensure accurate alignment of anatomical structures across imaging modalities, such as finite element methods and B-spline methods used in deformable image registration [10, 11]. In the meantime, routine clinical practice is greatly benefiting from numerous works that implement mathematical or machine learning models for automating image segmentation, such as convolutional neural networks (CNNs) that are widely used for image autosegmentation by learning features directly from imaging datasets [80]. Currently, such machine learning and hybrid models have largely outperformed atlas-based methods on autosegmentation and are widely used in commercial and in-house clinical solutions for automated segmentation and RT treatment planning. For radiotherapy delivery, mathematical modeling has also been critical in ensuring or improving delivery accuracy, such as for respiratory motion prediction and management during the treatment [81]. Most commercial respiratory gating systems integrate such predictive filters though their utility is currently more sporadic in the clinic due to imperfect robustness against breathing irregularities. On another modern RT frontier, image-based tissue-scale models leverage quantitative imaging modalities, such as MRI and PET, to optimize radiotherapy by integrating biological insights with patient-specific data. Hormuth et al. utilized dynamic contrast-enhanced (DCE)- and diffusion-weighted (DW)-MRI to predict tumor and vascular responses to radiotherapy in glioma models, achieving spatially accurate forecasts of treatment outcomes [82,83,84]​. Rockne et al. incorporated hypoxia-informed resistance using 18F-FMISO PET in a computational model for glioblastoma, improving the accuracy of tumor radiosensitivity predictions [85]. Lipková et al. combined FET-PET and MRI in a Bayesian framework to estimate tumor densities and guide personalized dose-escalation strategies while minimizing radiation toxicity. These approaches demonstrate the potential of image-driven models to refine radiotherapy by enabling adaptive, patient-specific treatments [86, 87].

Mathematical approaches and models for temporal and spatial fractionation

Radiotherapy dose fractionation favors normal tissue repairs, for the reason that tumors with a high α/β ratio (e.g. 10 Gy) tend to be early-responding, i.e. show radiation damage soon after exposure due to rapid cell turnover, while most organs (excluding skin) and normal tissues (excluding mucosal surfaces and marrow) with a low α/β ratio (e.g. 2–3 Gy) tend to be late-responding, i.e. show radiation damage long after exposure due to slow cell turnover. Models for early-responding tissues usually prioritize total dose effects and acute cell killing, such as using the LQ model with higher α/β ratios, and those for late-responding tissues tend to emphasize repair kinetics and long-term damage accumulation, incorporating lower α/β ratios and repair-misrepair parameters. Hypofractionation, radiotherapy delivered in one or limited numbers of fractions of large fractional doses, has gained traction in radiotherapy for cancer in recent years. In addition to enhancing patient convenience and cost-effectiveness, hypofractionation can be more effective in terms of cancer control in some circumstances. This greater efficacy may reflect differences in the radiobiology of hypofractionation. The effectiveness of dose fractionation is generally well modeled by classical models such as LQ model and BED calculation. However, as the fractional dose gets above a certain threshold (8-10 Gy), such as in the situations of SRS/SBRT, it becomes controversial how accurately these classical models could model the dose–effect and if they would overestimate cell survival. Biologically, such high fractional doses may overwhelm or saturate the DNA repair mechanisms, damage vasculature, and trigger immune effects, thereby violating the LQ radiation cell survival model. The effect of hypoxia may be lessened with hypofractionation. Also, some cancers (such as prostate cancer) behave more like late-responding tissue than early-responding tissue. Modifications to classical models and new model development are therefore necessitated to more accurately capture cell response to hypofractionated treatments. For example, an additional linear term that becomes dominant at higher doses is introduced in the linear-quadratic-linear (LQL) model which transitions the cell killing from quadratic to linear surpassing a crossover dose threshold [88]. This modification accounts for the saturation of sublethal damage repair mechanisms at high dose levels, and therefore allows the LQL model to better represent the biological response in scenarios involving high-dose fractions, such as SBRT. The universal survival curve (USC) model is another model that uses LQ model for lower doses and a linear model of two components for higher doses, combining the strengths of the LQ and multitarget models [89]. The transition in the USC model ensures the model's applicability across a wide dose range, improving its accuracy in describing cell survival in hypofractionation regimens.

Watanabe et al. developed a simple mathematical model to simulate tumor growth and response to a single high-dose fraction of irradiation, using key biological parameters such as the radiation sensitivity indicator α from the LQ model and a new vascular growth retardation factor θ to predict post-treatment tumor volume changes in both experimental and clinical settings [90]. Matsuya et al. compared the LQ model and two other microdosimetric-kinetic models (MKM) using a Markov chain Monte Carlo simulation, and found that damage repair during irradiation plays a key role and models without a repair factor during a short dose-delivery time may overestimate cell killing in fractionated radiotherapy [91]. Nakano et al. introduced a mathematical model combining ODE and MKM to predict the tumor cell lethal effect in non-small cell lung cancer (NSCLC) during SBRT and validated their approach with in vitro experimental data [92]. Tariq et al. developed a model to describe tumor volume dynamics during SBRT, focusing on ODE models rather than more complex PDE models [93]. This decision balanced model accuracy with practical clinical applicability, as the available clinical data did not include spatial distribution information. They explored one-population and two-population tumor growth models, employing the Akaike information criterion to identify the most appropriate model by balancing accuracy and complexity. The study concluded that the two-population exponential growth model provided the best trade-off, making it a promising candidate for future clinical applications.

Mathematical modeling has been used to examine conventional, hypo-, and hyper-temporally-fractionated and spatially fractionated radiotherapy, providing support for the clinical and preclinical evidence of their effectiveness. Table 3 lists some examples of studies that apply mathematical modeling to study these effects. Mathematical models are continually evolving in tackling the complexities of dose distribution and biological response in temporally and spatially fractioned RT. Although the LQ model remains widely used for its simplicity and applicability, there is growing interest in incorporating more sophisticated models that account for factors such as the fractional dose effects, hypoxia, and cell heterogeneity. Ongoing research aims to refine the mathematical modeling techniques to optimize treatment plans, improve TCP, and personalize treatment schedules based on patient-specific factors and comprehensive radiobiological insights [42, 88, 94, 95].

Table 3 Example studies applying mathematical modeling to study the effectiveness of temporal and spatial fractionations of radiotherapy

Mathematical approaches and models for FLASH radiotherapy: the ultra-high dose rate

FLASH radiotherapy (FLASH-RT) is an emerging radiotherapy concept with potential to deliver therapeutic doses of radiation with minimal side effects [105,106,107]. Characterized as an often-single-dose radiotherapy of ultra-high dose rate, > 40 Gy/s as opposed to the 0.01–0.4 Gy/s in conventional radiotherapy (CONV-RT), FLASH-RT has been shown in many preclinical studies to have reduced risks of normal tissue damages and hence improved healthy tissue sparing compared with CONV-RT [105, 108,109,110]. The FLASH effect is not a new concept, first described in the late 1950’s with reduced radiosensitivity of bacteria when exposed to ultra-high dose rate radiation in an anaerobic condition [111]. The newly ignited interests were built on a series of in vivo animal experiments in which FLASH-RT substantially reduced the healthy tissue toxicity while maintaining tumor control. Currently, FLASH-RT is actively researched on using photon and particle therapy systems such as proton, electron, and heavy ions, with pilot human clinical trials currently underway such as the FAST-01 trial [109, 112,113,114,115,116].

Despite numerous in vivo and in vitro demonstrations of the FLASH effect, the mechanisms behind the FLASH-RT are still largely controversial [117,118,119,120]. There are multiple popular theories. The most widely-accepted is likely the oxygen depletion theory, which can be traced back to the original FLASH effect paper [111], in which FLASH-RT is believed to create a transient hypoxic state in normal tissues through oxygen depletion, causing the tissues to become more radioresistant and hence better protected against the radiation, versus tumor tissues that are already in a baseline hypoxic state where additional effects are not expected [121,122,123,124]. Extensions of this theory include the reoxygenation dynamics in tumor and tissue environments with varying levels of oxygenation, the intrinsically different levels and responses to reactive oxygen species (ROS) between tumor and normal tissues, etc. Another popular theory is based on different tumor and normal tissue immunological responses to FLASH-RT, including their different immunogenic characteristics, the immune microenvironment, pro-inflammatory and anti-inflammatory responses, and the immune system's reaction to ionizing radiation [124,125,126,127]. A third main theory relates to DNA damages and responses to these damages between tumor and normal tissues under the FLASH-RT vs. CONV-RT conditions [128,129,130,131]. In addition, there are many other speculations based on the unique temporal and spatial dose distribution of FLASH-RT, as well as its damages to the vasculature, etc.[117,118,119]. Current FLASH-RT studies have not conformed the Bragg peak to tumor volume and have not verified dosimetry, which elicits hesitation in its clinical use [120]. Nevertheless, there is still much unknown and many controversies regarding the exact mechanism of the FLASH effect. Therefore, mathematical modeling provides a useful method to study possible physical, chemical, and biological interactions to gain a deeper understanding and insight into this complex system.

Based on the oxygen depletion theory that FLASH-RT leads to transient hypoxia in normal tissue cells, Zhou et al. employed dimensional analysis to estimate the minimum dose rate required for achieving the FLASH effect [24]. In this work, they divide the FLASH irradiation process by ultra-short radiation pulses into three phases: radiation-induced species production, oxygen molecule depletion, and reoxygenation through diffusion. Using differential equations describing these processes and applying dimensional analysis, they suggest that the estimated minimum dose rate for pulsed FLASH-RT is proportional to the product of the oxygen diffusion coefficient and intracellular oxygen concentration and inversely proportional to the square of the oxygen diffusion distance and the reduction of intracellular oxygen per unit radiation dose. Their estimation of the minimal required dose rate aligns with the order of magnitude of dose rates observed in previous FLASH radiotherapy experiments, providing a theoretical framework and preliminary estimates for FLASH-RT.

In another work, Abolfath et al. studied the molecular dynamics in FLASH-RT using Monte Carlo simulation of physical reactions including ionizations, electronic and vibrational excitations, and ROS generation and reactions [121. Applying analytic and numerical methods for a system of dimensionless differential equations with analytical and numerical solutions, their time-dependent molecular dynamics simulations establish the time scale of oxygen depletion in FLASH-RT to be femto- to nanoseconds after irradiation, and suggest that the FLASH effect takes place in physoxic normal tissues (4–5% oxygen) compared to the hypoxic tumor tissues, but not at higher oxygen levels above 10–15% oxygen.

Hu et al. applied PDE-based mathematical modeling on oxygen distribution in tissue based on the oxygen consumption rate of tissue and the distance between capillaries [122]. Their findings suggest a plateau effect that once above a certain threshold, further increasing the dose rate will not significantly alter the oxygen concentration change. They also suggest that based on the oxygen depletion theory alone the model results contradict the observed FLASH effect in some settings, especially in brain tissue, necessitating alternative hypotheses to complement oxygen depletion to fully elucidate the FLASH effect.

Rothwell et al. also applied mathematical models similar to the others based on oxygen depletion to study the FLASH effect, focusing on identifying the conditions under which radiation may induce sufficient depletion of oxygen to cause a diffusion-limited hypoxic cellular response [123]. They use a typical reaction–diffusion model in a one-dimensional slab geometry, combining Fickian diffusion with terms for metabolic and radiation-induced oxygen consumption. The radiation-induced consumption of oxygen is modeled using a two-stage lumped reaction, reflecting the physicochemical stages of radiolysis of water. With the model, they extensively analyzed how different parameters, such as dose rate, metabolic reaction rate, and oxygen concentration etc., influence the FLASH effect. The simulations show how these parameters interact and their cumulative impact on oxygen depletion and reoxygenation in cells.

Zhu et al. also used Monte Carlo simulations and mathematical modeling to investigate the impact of radiolytic oxygen depletion (ROD) on cellular responses to FLASH-RT [132]. Monte Carlo simulations were used to model the ROD process, study the interactions between radiation and biological tissues, and calculate oxygen depletion and DNA damage. Oxygen distribution and cellular responses were analyzed using numerical simulations. Similarly, they suggest that the FLASH effect due to ROD happens only in hypoxic cells with lower oxygen levels (pO₂ < 30–40 mmHg). They also conclude that single pulse radiation and multi-pulse radiation with shorter pulse intervals are better at achieving the FLASH effect without suffering from oxygen tension recovery. In addition, they suggest low LET source particles are preferred for FLASH-RT.

Similarly, Song et al. employed analytical mathematics to model the effects of varying LET on the FLASH effect based on oxygen depletion [133]. Integrating the oxygen depletion framework and the oxygen enhancement ratio values according to LET, they suggest that the FLASH effect is maximized and observable at intermediate oxygen level (10–100 mmHg), and low LET radiation is more effective at inducing the FLASH effect.

Mathematical modeling with PDE, ODE and other approaches, including molecular dynamics and Monte Carlo simulations, has provided insights into the dose rates required for FLASH-RT and its effects on tissue oxygenation. Still, more research is yet ongoing and needed to fully understand the FLASH effect and refine its clinical application.

Mathematical approaches and models for integrating radiotherapy with immunotherapy and immune-oncology

While radiotherapy is primarily a local therapy, the circulatory system can allow for possible systemic effects of radiation. Radiation can have both a stimulatory and a suppressive effect on the immune system. Radiation-related fatigue is an example of a systemic effect thought to be mediated by radiation-induced cytokines that circulate through the body. Similarly, radiotherapy can modulate the immune system and responses via cytokines as well as circulating immune cells. Thus, radiation can contribute to systemic, immune-mediated anti-tumor effects. Specifically, with its cytotoxic effects, radiation can trigger the release of tumor antigens and inflammatory factors from the dead tumor cells, leading to subsequent immune responses [134, 135]. Moreover, radiation can also induce immunogenic cell death, directly activate or suppress immune cells, modulate tumor microenvironment, and alter its immune responsiveness. For instance, in some cases, radiotherapy has been found to lead to an abscopal effect, where radiation treatment of a tumor at one site leads to recession of metastatic cancer at different sites [134]. The effect is believed to be mediated by the immune response activated by the primary site radiation that also targets distant tumor cells. Perhaps even more important, crucial for the immunomodulatory effects of radiation is the radiation-resulted tumor DNA breaks and ROS, and the cascade of biological events they trigger.

Once tumor DNA and damage associated molecular patterns (DAMPs) are sensed in the immune system, cells such as phagocytes, B cells, and T cells can mature, activate, and respond [134]. For example, phagocytes help to clean up cellular debris while antigen-presenting cells and dendritic cell exposure to tumor DNA and DAMPs allow for antigen-presenting function and release of cytokines to stimulate other immune cells. Similarly, T-cell exposure to DNA and DAMPs can help with tumor infiltration and the production of inflammatory cytokines [134]. With increased chemotactic factors being released, increased T cells infiltrate the tumor potentially increasing (and also potentially decreasing) the immune response to radiation [134, 136]. It should be noted, however, that prolonged DNA and DAMPs can activate apoptosis in T cells mitigating their anti-tumor effects [137]. The abscopal effect is mediated by the same antigen-presenting mechanism which activates T cells that can travel throughout the body to regulate tumor metastasis cells. Supporting this, Wang et al. studied melanoma in mouse models and found evidence of the abscopal effect due to increased responsiveness of tumor cells to immunotherapy following radiation treatment [138]. However, interestingly, the abscopal effect has also been observed in immuno-deficient mice and not only with RT, but also with surgery, indicating that the immune system is not the only mediator of this effect, and it is also not a unique feature of RT [139, 140]. Additionally, as noted by Walker et al., the interconnectedness of metastatic tumors through systemic T-cell redistribution suggests that changes at one site, whether from surgery or radiotherapy, can provoke systemic responses, further supporting the complexity of the abscopal effect as a multifaceted phenomenon [141]​.

Mathematical modeling can play a pivotal role in exploring the integration of radiotherapy with immune-oncology and immunotherapy, particularly given the complexity of the interactions involved. Radiation has the potential to result in either immunogenic or immunosuppressive effects. Radiation can interact with the tumor cells, the tumor microenvironment such as tumor vasculature and endothelial cells, and nearby healthy cells. Various factors such as dose, fractionation (hyper-, conventional, or hypo-fractionation), timing (concurrent or sequential, different delay intervals if sequential), and dose rate (conventional or FLASH) of the radiation, oxygen effects, LET can all affect the interactions.

To address these complexities, multiple mathematical models combining radiation with immunotherapy or investigating the immune-modulating effects of radiation have been proposed. Mathematical modeling has also extensively studied the balance between preserving healthy cells and killing cancer cells in the context of immunoradiotherapy, for instance, the impact of dose distribution and dose rates on the level of circulating lymphocytes. Table 4 lists some examples of these studies. Collectively, they provide further evidence for combining immunotherapy and radiotherapy treatments.

Table 4 Example studies exploring the effect of radiation on the immune system and the relationship between radiotherapy and immunotherapy

Overall, mathematical modeling provides a necessary tool for studying the complex relationship and interplay between radiotherapy, immunotherapy, and immune-oncology.

Mathematical modeling and digital twins

A promising new direction of research and application is digital twinning. A digital twin usually refers to a virtual copy or model of a physical entity [155,156,157]. The digital twin mirrors the physical entity using real-time data to provide real-time monitoring, detailed analysis, design, planning, and optimization for the real-world entity or process. Central to any digital twin is a proper model, with good accuracy and incorporating reasonable uncertainty, allowing the digital twin to effectively mirror and interact with its physical counterpart. Mathematical models, therefore, play an integral role in digital twinning, by describing mechanism-based models, characterizing uncertainties, simulating varying levels of complexity, integrating multiple scales of data, etc.

In radiotherapy, digital twin research and applications are only starting, but the field is expected to grow rapidly in the next few years. As described in the previous sections, the interplay between radiation and tumor, healthy, and immune cells is complex, which is further complicated by the sophistication of patient biology on an individual level and the modulation of radiation techniques and variations. Mathematical models, both mechanistic and machine learning ones, can therefore serve as vital tools in digital twin research and applications. In a pioneering study, digital twinning was employed to assimilate individual patients’ MRI data to personalize the growth and response parameters of the biology-based mathematical tumor growth, disease progression, and patient survival models for high-grade gliomas [158]. The digital twin based on MRI data before and during radiotherapy time points was used to predict outcomes post-radiotherapy. Eventually, such a system would then be able to adapt radiotherapy regimens and plans on an individual level based on predicted responses. In radiopharmaceutical therapy and theragnostic, because radiation dosimetry depends both on radiation isotope decay and particle transport and on dynamic radiotracer accumulation and washout in various tissues and organs, digital twins based on Monte Carlo simulations and pharmacokinetic modeling have provided new directions for accurate and individualized radiation dosimetry [159,160,161]. These advancements are then poised to personalize and optimize such radiotherapy applications. The digital twin concept has also been introduced to FLASH-RT modeling as well as proton adaptative radiotherapy incorporating anatomical uncertainty and variability [160, 162].

Challenges, opportunities, and outlook

Mathematical modeling has been an essential tool in radiobiology and prescribing radiation dosimetry in cancer radiotherapy, as well as advancing the precision, personalization, and overall effectiveness of radiotherapy. Current challenges include the lingering uncertainties in the models, the difficulty in clinically translating more advanced models, and the complexity in effectively integrating different data scales to navigate tumor and patient heterogeneity.

Radiotherapy clinical practice is largely based on classical radiobiology, dose calculation, and plan optimization models where uncertainties are often overlooked or addressed using a margin concept, i.e. adding margins to the gross tumor volume in the RT target definition to account for uncertainty and ensure the tumor receives the intended dose. Along with the recent advancement of radiation delivery accuracy and sophistication, these uncertainties began to be acknowledged and explicitly addressed in the models [163,164,165]. Even newer models also continue to involve uncertainties and could lead to poorer performance on heterogeneous datasets or populations [166,167,168,169]. Partially owing to such uncertainties and the model complexity, the more sophisticated models have experienced a slow clinical translation to replace classical models in standard clinical practices. Lastly, while multi-scale, multi-modal data integration is highly desirable for next-generation radiotherapy research and applications and is amenable to mathematical modeling approaches, such endeavors are still in the infancy.

Nevertheless, this is an exciting time for research in mathematical modeling within the field of cancer radiotherapy. Besides using these approaches to better elucidate mechanisms of advanced radiotherapy modalities such as FLASH-RT, SRS/SBRT, SFRT, and immune-radiotherapy, they can also be applied to leverage the ever-expanding available medical data to further optimize radiotherapy and cancer care. Mathematical modeling can be used to optimize radiotherapy and maximize therapeutical ratio, integrating functional and molecular imaging data of the patient, radiomics and other deep learning extracted features from these and other intravital patient images, genomic, radiogenomic, and many immunotyping information obtained through single-cell sequencing or liquid biopsies, and information from various wearable or implantable devices and fitness trackers. Building on the mathematical models that are already widely used in various aspects of radiotherapy and cancer biology, further development and integration of these models using multi-scale, multi-modal, digital twins and other systems approaches is expected to usher radiotherapy into a new era of personalization and optimization.

Conclusion

Mathematical modeling plays a crucial role in radiotherapy for cancer, with applications spanning radiobiology, clinical workflows, SRS/SBRT, FLASH-RT, immune-radiotherapy, radiotherapy digital twins, and beyond. This narrative review explores these areas, focusing on emerging trends and innovations. As the field advances, the importance of mathematical modeling will only grow, driving further improvements and integrating multi-scale and multi-model medical data in precision and personalized cancer treatment approaches.