Mathematical modeling of the synergistic interplay of radiotherapy and immunotherapy in anti-cancer treatments

Introduction While radiotherapy has long been recognized for its ability to directly ablate cancer cells through necrosis or apoptosis, radiotherapy-induced abscopal effect suggests that its impact extends beyond local tumor destruction thanks to immune response. Cellular proliferation and necrosis have been extensively studied using mathematical models that simulate tumor growth, such as Gompertz law, and the radiation effects, such as the linear-quadratic model. However, the effectiveness of radiotherapy-induced immune responses may vary among patients due to individual differences in radiation sensitivity and other factors. Methods We present a novel macroscopic approach designed to quantitatively analyze the intricate dynamics governing the interactions among the immune system, radiotherapy, and tumor progression. Building upon previous research demonstrating the synergistic effects of radiotherapy and immunotherapy in cancer treatment, we provide a comprehensive mathematical framework for understanding the underlying mechanisms driving these interactions. Results Our method leverages macroscopic observations and mathematical modeling to capture the overarching dynamics of this interplay, offering valuable insights for optimizing cancer treatment strategies. One shows that Gompertz law can describe therapy effects with two effective parameters. This result permits quantitative data analyses, which give useful indications for the disease progression and clinical decisions. Discussion Through validation against diverse data sets from the literature, we demonstrate the reliability and versatility of our approach in predicting the time evolution of the disease and assessing the potential efficacy of radiotherapy-immunotherapy combinations. This further supports the promising potential of the abscopal effect, suggesting that in select cases, depending on tumor size, it may confer full efficacy to radiotherapy.


Introduction
Immunological experiments during the last two decades have answered many important questions related to the causal relationship between chronic inflammation and carcinogenesis.The presence of inflammatory cells in the cancer milieu raises the question of the tumor progression despite a likely immune system reaction to tumor antigens.This aspect is particularly important since untreated tumors grow according to non-linear, macroscopic, laws as the Gompertz law (GL) Gompertz (1); Norton (2); Vaghi et al. (3) or the logistic one (LL) Verhulst (4); Vaghi et al. (3).Therefore those growth patterns emerge, at a larger level of magnification, from many microscopic biological factors, which turn out to be summarized by simple mathematical descriptions.
Since prolonged inflammation is a hallmark of cancer Hiam-Galvez et al. (5), initiating tumor genesis or supporting tumor growth, and the global immune response is significantly altered during tumor progression, immunotherapy is becoming a valid option in cancer treatment.However, the immune response can be detrimental rather than helpful [see, for example, Lin et al. (6)]: individual auto-antibodies play an antagonist role in cancer, but the agonist auto-antibody in some cancer patients turned out to be deleterious and harmful.Some preclinical and clinical evidence confirm the synergistic action of radiotherapy (RT) and immunotherapy against the tumor cells Zhao and Shao (7).Although The intrinsic sensitivity to radiation is patient-specific Puglisi et al. (8); Puglisi et al. (9) and may depend on different factors, RT is able to ablate cancer cells not only by directly induced necrosis or apoptosis but also by triggering an immune response that actively recruits immune cells within the tumor microenvironment.For example, RT promotes the release of tumor-associated antigens, which, once processed by antigenpresenting cells (APCs), prime CD8 + and CD4 + T cells in the draining lymph nodes.These lymphocytes attack both primary tumor and metastatic sites, posing the biological basis of the in situ vaccination driving the so-called abscopal effect Ngwa et al. (10); Mole (11); Demaria et al. (12): RT induces a systemic behavior that can activate the immune response against metastasis, i.e., in locations that are far from the RT-treated primary tumor.
The involvement of the immune system has been demonstrated in different experimental models such as melanoma Twyman-Saint Victor et al. (13), colorectal Dovedi et al. (14) and breast cancers Demaria et al. (15), but the clinical presentation where considered anecdotal or at least rare.
Its rarity in clinical practice is likely due to the simultaneous engagement of immune escape mechanisms, such as the recruitment of regulatory CD4 + (Treg) and myeloid-derived suppressor cells counterbalancing the anti-tumor CD8 + T cellmediated effects, and the tumor release of hypoxia-inducible factors with pro-survival activity Ji et al. (16).
More generally, susceptibility to the abscopal effect has been associated with several biological factors, such as tumor size or oxygen levels in tumor tissues.Indeed, the presence of hypoxic regions results in both increased resistance to the lethal effect of radiation mediated by reactive oxygen species (ROS) production and an immune suppressive tumor landscape ruled by Tregrecruiting chemokines and impaired APC function McNamee et al. (17); Castorina et al. (18) Mathematical modeling approaches have become increasingly abundant in describing immunotherapy and its synergy with RT Dewan et al.  (35) .
A large part of mathematical models on tumor growth and therapies are based on sets of coupled differential equations.The number of equations increases together with the details of the biological description and this implies a large number of parameters and initial conditions to be specified.The detailed analyses are often so complex to require surrogate models for a reliable determination of the parameters Browning and Simpson (36).
On the other hand, more economical models introduce a macroscopic evolution of tumor growth and therapy with fewer parameters and a coarse-grain dynamical evolution Norton (2); Wheldon ( 37 This is the key point to obtain an effective quantitative control of the tumor progression.Indeed, microscopic models, which have the advantage of a deep understanding of the biological dynamics, of the physiologically-based pharmacokinetic [see for example Maaß et al. (41)] and of the possible translation to different populations/diseases, require many parameters and, although some of them can be determined by previous analyses, the parametric error propagation will produce a large band of fluctuation in the prediction of the quantitative evolution of the disease.Therefore we prefer to apply macroscopic growth laws of the sigmoid family with two parameters.Our choice of the GL is due to the result that untreated tumor growth has been better described by it (see Vaghi et al. (3) for a recent study).Moreover, in a transplantable rat tumor, it was shown that control and regrowth curves after radiotherapy could be fitted by the same Gompertzian law, provided adjustments for the initial lag and the estimated number of clonogens immediately after irradiation were performed [Jung et al. (42)].Gompertzian growth has been assumed to describe human tumor repopulation during fractional radiotherapy also in Hansen et al. (43) and by O'Donoghue (44).
The main motivation of the proposed approach is, in our opinion, its complementary role in the clinical evaluation of disease progression, often based on macroscopic variables, and the better parameter identification, thus increasing model verifiability Braakman et al. (45).
Finally, this method offers a clear description of the complex interplay between radiotherapy, immunotherapy, and tumor progression, providing insights for advancing cancer treatment strategies that harness the abscopal effect.
The paper is structured as follows: the mathematical model, based on the GL and on the definition of the effective GL parameters, is recalled in the next section (Appendices A, B and C contain the corresponding calculations).Different macroscopic growth laws (as the LL) can be applied, without changing the underlying method.The emerging phenomenological approach, based on the suitable redefinition of the two GL parameters to describe the data, is reported in Section 3. The final sections are devoted to discussion and to the possible clinical use of the phenomenological model.

Methods
The proposed method is based on the mathematical model reported in detail in Appendices A, B and C in the Supplementary Materials.In what follows, the assumptions and some exact results are reported.Then, the phenomenological model is discussed.

Mathematical modeling
A general classification of macroscopic growth laws is reported in Castorina et al. (39) Castorina and Blanchard (46).For a population N(t) they are solutions of a general differential equation that can be written as where f(N) is the specific growth rate and its N dependence describes the feedback effects during the time evolution.If in Equation 1 f(N) becomes constant, the growth follows an exponential pattern.
The untreated tumor progression is described by the GL Norton (2); Vaghi et al. (3), solution of the previous equation (see Data Sheet 1) with where a,k,N 0 are constants that respectively indicate the exponential growth, the limiting factor, the initial cell number and N ∞ is the carrying capacity (N ∞ = N(0)exp(a/k)).
For untreated tumors, the GL emerges from microscopic, biological mechanisms where natural/adaptive immunity is taken into account Berendt and North (47); Gonzalez et al. (48); Castorina and Carco' (49) The further effects of immune therapy, I(t), can be described by a modification of the previous Equation 2as follows (see for example Wheldon (37)) where g is a constant.Notice that the sign of g indicates the agonist or antagonist effect of the immune response: a negative g increases the specific growth rate.The variable I(t) generically refers to the passive immunity resulting from the injection of anti-cancerspecific monoclonal antibodies i.e., the drug effects.However, it is, in general, unknown and requires a specific model.An example is the model of immunotherapeutic drug T11 target structure in the progression of malignant gliomas Khajanchi and Ghosh (50) Khajanchi and Banerjee (51).
The general solution of the previous equation and the IT effects on the tumor progression are discussed in Mathematical Formalism.For illustrative purposes, two specific cases are analyzed: I(t) = I(0) = constant and I(t) = I(0)exp(−rt).The role of the therapy can be assimilated to a redefinition of an effective carrying capacity (in the first of the two cases) and, in general, in the introduction of effective, time and therapy dependent, parameters a eff , k eff or N eff ∞ , k eff , whose quantitative relation with I(t) is given in Mathematical Formalism.The introduction of an effective carrying capacity is wellknown in population dynamics Royama (52).For example, the invention and diffusion of technologies lift the growth limit.Its possible time dependence is usually included by (at least) another differential equation, coupled with the growth equation.This will increase the number of parameters and initial conditions and, therefore, in our computational method the two effective parameters will be determined by data fits, giving a phenomenological indication about the disease progression.
There are different outcomes following an immune response to cancer Lin et al. (6): i) a surveillance role that inhibits the initiation and progression of the cancer; ii) the possibility that under certain conditions the immune response may nourish rather than curtail tumor growth.In other terms, monoclonal antibodies can exert antagonistic as well as sympathetic effects on tumor growth Lin et al. (6).
This possibility translates into a direct comparison among the parameters that describe the immunotherapy effects in Equation 3and its solutions and the available data.For example, the determination of the crucial sign of the constant g.
Let us now consider the combination of immune therapy (IT) and radiotherapy (RT).As a first step, we study the case of independent effects, i.e. no synergy between IT and RT.
The effect of radiotherapy is described by the linear quadratic model (LQM).Denoting by N(t − ), N(t + ) respectively the cell number before and after the single dose d at time t, the RT effect is given by where a and b are constants (numerically b ≃ a/10) and d is the dose, Van Leeuwen et al. (53).The result in Equation 4 assumes an instantaneous effect of the RT, which could be, in general, not strictly applicable.Also, In this case, the therapy and immune response effects can be translated in the definition of effective parameters of the GL (see Mathematical Formulation).
The number of tumor cells after n f treatments at time t n ,t n+1 ,t n +2 ,…t nf turns out to be (see Mathematical Formalism) with the functions D n f and W(t n f , t 0 ) as described in Mathematical Formalism.
According to Equation 5, if the tumor cell number decreases and the time evolution of the diseases moves toward complete recovery.
Therefore, although RT and IT are considered independent, if the effects of RT are such that then, a small impact of the immunotherapy, W, can produce a tumor volume regression.Moreover, the critical conditions in Equations (6, 7) depend on the fractioning of the radiotherapy, since different schedules give different values of D n f and W(t n f , t 0 ).Therefore, the previous conditions correspond to optimal control of the therapy effects, Khajanchi and Banerjee (54); Khajanchi (55); Khajanchi and Banerjee (56).

The synergy between immune and radio therapies
In the macroscopic framework, the description of the synergy between RT and IT requires a new term in the specific rate, which takes into account the immune response activated by RT, i.e.
where d is a constant, F(d,t) is a function of the dose d and of the time series of the treatments on the tumor, quantifying the cellkilling effect of the adaptive immune response Y (t), triggered by the RT.Y (t) is different from I(t) as it represents the outcome of the active immunization due to antigenic peptides coming from the disintegration of tissues hit by radiotherapy, and following the inflammation.To be more specific, Y (t) represents the immune response to tumor-associated antigens, promoted by the inflammation context due to the damage perpetrated by RT.Also, it must be specified that this immune response has a chance to exert an effect only before evasion mechanisms are established by the tumor (factors that are not counted in the present model).
If d = 0, F(d,t) = 0 there is no synergy.The specific form of the function F(d,t) requires a microscopic model, however one expects that the coupling Y (t)F(d,t), i.e. the immune activation due to RT, has a typical time decay, t, after the single dose radiotherapy described by the LQM.For the primary tumor, the parameter d is small, i.e., d << 1 and the synergy is small.The abscopal effect is described by considering a finite value of d for metastases that are far away from the primary tumor location.The result for the time evolution of the abscopal effect is given in Mathematical Formalism.

Results
According to the mathematical approach in the previous section a phenomenological, simplified, method of analysis of the experimental data emerges.Indeed, a large part of the therapy and immune response effects can be assimilated to a redefinition of the GL parameters, with time and treatment dependence, which turns out to be detailed enough to compare with data.This phenomenological, simpler, procedure facilitates the validation of the model [Braakman et al. (45)] with respect to more complex analytical approaches, as discussed later.
Moreover, the function I(t),W(t),Y (t) in the previous differential equations are largely unknown, therefore the data fits of the effective GL parameters (see the general formulas in Mathematical Formalism) give model-independent information about the role of IT and RT.

Analysis of experimental dataimmune therapy
In ref.Lin et al. (6) the different immune responses to cancer have been described.The authors highlight the agonist and antagonist effects of, respectively, AB93 and AB641 autoantibodies for the growth factor receptor TrkB in patients with breast cancer.After injection of MDA-MB-231 cells in immunodeficient mice they show the response to treatment, for different dosages of autoantibodies, measuring tumor growth.In particular, the data on the effects of AB641 and AB93 on tumor progression can be analyzed in the proposed macroscopic approach by Equation 3, by redefining GL parameters.Initially, one has to fit the untreated tumor progression data by GL to determine the corresponding parameters.Then, one repeats the analysis with immunotherapy.The results are shown in Figure 1, which clearly reveal the agonist or antagonist role of the different antibodies.
The experimental result in the case of therapy can be fitted by redefining the GL parameters with respect to the untreated ones (see Equation S2, Equation S23 and Mathematical Formalism in Supplementary Material).The values are reported in Table 1.The agonist effects increase both parameters, producing faster growth, corresponding to a negative g in Equation 3, whereas the antagonistic effect induces tumor depletion.Notice that the change in the GL parameters, i.e. of the therapy, implies a modification of the exponential rate and the carrying capacity with respect to the untreated tumor.

Analysis of experimental data -RT and abscopal effect
The abscopal effect has been experimentally studied in Nesseler et al. ( 58) by inoculation of undifferentiated fibrosarcoma cells (FSA1) into immunocompetent mice to simulate primary and metastatic conditions, successively divided in four treatment groups: no treatment, anti-PD-1 monoclonal antibody alone, RT alone, and combination of anti-PD-1 with RT.
Initially, only the effect of RT on the primary has been detected, showing a critical dose administration for tumor regression [Figure 2A   Comparison of the effective GL with data from the literature for untreated tumor and immunotherapy results.In the graph the symbols represent real data recorded and the dashed curves represent the Gompertzian fit, the colors are paired for both results.Fitted parameters are reported in Table 1.    2.

TABLE 1 GL effective parameters of the comparison in
to the untreated case, signal that the therapy produces a complete depletion of the size, analogously to the extinction in population dynamics.
The strong signal of the change of sign in the effective parameters can also be summarized by plotting the specific rate (1/NdN/dt), in Figures 4, 5, where its negative value indicates the complete tumor regression trend.
It should be stressed that the approach with GL effective parameters, suggested by the more rigorous previous mathematical model (see Mathematical Formalism), is a phenomenological one with the aim of a simplified clinical, but quantitative, understanding of the tumor progression at a more personalized level (see next discussion section).
The final analysis concerns the abscopal effect, according to Equation 8 with g = 0. Let us assume: a) an exponential decay of the immune response with a time delay t between the RT on the primary and the immunological effect on the secondary; b) the activated immune system continues its effect on the secondary with an exponential rate corresponding to the specific rate obtained at the end of the RT.
In other words, if t in is the starting day of RT on the primary, for t < t in the metastatic site evolves according to the GL progression with the untreated parameters.At t in the RT starts and the immune system targets the secondary (see Equation 8and its solution reported in Mathematical Formalism).At the end of the RT, the immunity response continues to reduce the metastasis with an exponential behavior if the specific rate turns out to be negative.
In Figure 6 the result is depicted, i.e., the abscopal effect, according to the previous approach, for different values of the parameter I(0), determining the response of the immune system on the secondary induced by RT (see Mathematical Formalism).

Discussion
The phenomenological approach, based on the previous mathematical formulation, consists in fitting the specific rate data by GL with effective parameters.Let us recall that the specific rate is much more reliable than the volume tumor variation, in determining its progression and the phase of growth or decrease.
We are aware that the coarse-grain proposed approach misses the detailed dynamics and can be considered an oversimplified description since the underlying pathways, Ng and Dai ( 22 8.However, one has to recall that, independently of the microscopic conditions, a large part of untreated tumors follow the GL (see Vaghi et al. (3) for a recent review) and that with a small number of parameters, one gets quantitative clinical indications for personalized treatments.If, for example, a patient gets a small specific rate by RT on the primary tumor, then a small contribution of IT might be able to result in a complete recovery.
The abscopal effect has been described by a macroscopic coupling between RT and immune system response, where the initial progression of the metastasis follows the GL with untreated primary tumor parameters.This is a reasonable assumption although the in-situ conditions can produce different results.
According to the specific conditions recalled in the introduction, choosing the best dose fractionation and timing with respect to immunotherapy is difficult.Notoriously, the use of protracted RT schedules (standard fractionation or slight hypofractionation) is discouraged since radiosensitive lymphocytes are cleared out from tumor tissues at each fraction delivery, thus preventing their antitumor function,      (68).With the latter approach, complete responses have been documented earlier than with classic homogeneouslydelivered stereotactic RT doses and before ICI administration, likely implying rapid immune intervention enhanced and maintained by the addition of IT Ferini et al. (69).
All of the above considerations require the modeling of tumor response to RT and IT to help predict the best combination strategy, also given the inadequacy of current radiobiological mathematical models to comprehensively explain the results deriving from this association Ferini et al. (69).
In the proposed mathematical and computational approach, the fractionization effects are taken into account by the functions D n f and W(t n f , t 0 ) and a spatially non-homogeneous behavior can be easily implemented.However, a complete discussion requires a forthcoming devoted analysis.

Conclusions
A comprehensive scope of the combined impact of radiotherapy and immunotherapy is vital for clinical decision-making.Yet, numerous mathematical models in existing literature prove overly complex, characterized by an abundance of differential equations, parameters, and initial conditions, rendering their practical implementation quite challenging.
In this study, we use a macroscopic mathematical approach that does not rely on the underlying microscopic dynamics.Instead, we propose a simplified model of the tumor progression using a Gompertz law, which involves just two parameters.Moreover, utilizing numerical solvers of the equations provided in the appendices is a straightforward process.
Examining how control influences system dynamics under radiation and/or drug therapy sheds light on the disease's temporal progression.This approach holds promise in assisting clinicians to make informed decisions by providing a clearer understanding of treatment outcomes, namely, assessing whether the therapy administered results in full recovery.Progression on the metastatic volume, triggered by RT on the primary, for different values of the coupling between the immune system and RT (see Mathematical Formalism).
Besides its clinical relevance, our approach shows potential for additional experimental validation of synergistic impact of the abscopal phenomenon on treatment outcomes.
We recognize the importance of conducting rigorous investigations to solidify the theoretical basis of our approach.By performing targeted studies and gathering more empirical data, we aim to validate the approach in different conditions and the significance of the abscopal effect in influencing therapy outcomes.This experimental validation will provide a deeper understanding of the interplay between the administered treatment, tumor response, and the abscopal effect.
Moreover, our future research endeavors will focus on elucidating the mechanisms underlying the abscopal effect and quantifying its impact on treatment efficacy.By combining computational modeling with comprehensive experimental studies, we strive to enhance our understanding of this phenomenon and optimize therapeutic strategies accordingly.Ultimately, we aim to improve patient-oriented outcomes by harnessing the full potential of the abscopal effect in cancer therapy.
of ref.Nesseler et al. (58)].The limited role of anti-PD-1 on the untreated primary tumor has been observed [Figure 2B of ref.Nesseler et al. (58)] and it has been checked that the RT on the primary has no direct effect on the implanted secondary.Finally, the synergy between IT and RT has been verified by the regression of the metastasis.Let us first consider the data on the primary tumor, treated by RT only, with three treatments of 8 Gy on days 9,10 and 11 after the injection of the cancer cells.The qualitative analysis [Figure 2A of ref.Nesseler et al. (58)] clearly indicates that the LQM with instantaneous cell killing effect is not able to reproduce the observed effect, due to a delay between the treatments (on days 9,10 and 11) and the regression behavior (i.e., a negative specific rate), starting on day 15, Lim et al. (59); McMahon (60).Therefore, more complex dynamics are in place, which, however, can still be described by a GL pattern.As discussed, the effective parameters can be time-dependent (see Equation S23-S30 in Mathematical Formalism) due to the therapy and the results are reported in Figures 2, 3 respectively for the lower and upper data sets of Figure 2A of ref.Nesseler et al. (58).The corresponding fitted effective parameter values are given in Tables 2, 3.The ( * ) indicates that the fitted parameter a eff has a linear time dependence a eff → a eff (t − t 0 ).This time dependence and the sign change of k eff , compared

Figure 1
(see Equation S23 in Supplementary Material.

FIGURE 2 Effective
FIGURE 2 Effective GL fit of literature data with effective parameters of the lower data set of Figure 2A of ref.Nesseler et al. (58).Black point and the dashed line represent respectively real data and Gompertzian fitting.Blue squares and curve are data and GL fit after the first 8Gy RT treatment.The red rhombus and curve are data and GL after three treatments of 8Gy each.Parameters reported in Table2.

FIGURE 3 GL
FIGURE 3 GL fit of literature data with effective parameters of the upper data set of Figure 2A of ref.Nesseler et al. (58), where the untreated data are in black, blue for one shot of 8Gy of irradiation and red for three 8Gy shots.Parameter reported in Table3.

FIGURE 4
FIGURE 4 Specific rate of the lower data set of Figure 2A of ref.Nesseler et al. (58).The curves represent the specific growth rate day by day, the response to therapy is highlighted by a negative trend.

FIGURE 5
FIGURE 5 Specific rate of the upper data set of Figure 2A of ref.Nesseler et al.(58).The curves represent the specific growth rate day by day, the response to therapy is highlighted by a negative trend.

TABLE 2
GL effective parameters of the comparison in Figure2(see Equation S23 in Supplementary Material).

TABLE 3
GL effective parameters of the comparison inFigure 3 (see Equation S23 in Supplementary Material.
Parameters in day −1 .Statistical error based on c 2 per degree of freedom.