Radiotherapy remains one of the most effective local treatments for cancer. In addition to directly inducing cell death, ionizing radiation is now recognized to play a complex modulatory role on the tumor microenvironment [for a recent comprehensive review, see ref (1)]. In particular, radiation-induced DNA damage and the subsequent release of damage-associated molecular patterns (DAMPs) activate proinflammatory pathways and recruit innate and adaptive immune effectors, including macrophages and cytotoxic T lymphocytes. This immunomodulatory effect may have a significant impact not only on the primary tumor but also on distant, non-irradiated lesions, a phenomenon known as the abscopal effect (2, 3).

Experimental and clinical observations have reported abscopal responses under particular conditions, including high-dose irradiation schemes (4–14). In this context, FLASH radiotherapy—a novel technique delivering ultra-high rates—has emerged as a promising modality capable of triggering such systemic responses while sparing normal tissues (15–19).

Notably, preclinical studies on mice have documented a biphasic response following a single FLASH dose of 25 Gy to the primary tumor, characterized by an initial regression, a transient regrowth, and a delayed second phase of tumor reduction (20). For a recent review on preclinical studies with FLASH therapy, see ref (21).

Understanding the mechanisms underlying this temporal pattern is essential to elucidate the potential of FLASH therapy as a systemic treatment. However, direct modeling of the full complexity of the biological cascades involved—including immune cell trafficking, antigen presentation, and tumor–immune interactions—remains a formidable challenge due to the limited availability of detailed in vivo data.

To address this, we adopt a phenomenological modeling strategy grounded in experimental observations. The present work introduces a computational framework that builds upon prior in vivo studies of FLASH therapy and aims to capture both local and systemic tumor responses. The approach relies on two key components: 1) a previously proposed parametrization of FLASH-induced effects based on tumor volume data (22) and 2) a dynamical model incorporating tumor growth and its modulation via radiation-activated immune response. Importantly, the model is formulated to account for both primary and metastatic tumor sites, with attention to the eventually delayed, systemic effects of localized irradiation. The method is completely general and can be applied to any tumor phenotype.

Tumor growth is modeled using the Gompertz law (23–26), which captures the decelerating expansion of tumors due to intrinsic regulatory mechanisms. Radiotherapy effects are introduced as a time-dependent term modifying the growth rate (22), calibrated against FLASH experimental data in which the dose is given in short pulses at ultra-high dose rate on lung fibrosis in mice, based on a linear electron accelerator (20). The phases observed after irradiation are hypothesized to result from immune activation, and their impact on distant metastases is incorporated through a coupling term modulated by parameters governing the strength and timing of the systemic immune response.

By comparing the predicted dynamics of tumor volumes at primary and metastatic sites, the model allows for a quantitative assessment of how FLASH-induced immune mechanisms may lead to the observed abscopal effects. Therefore, this study provides a theoretical framework to interpret experimental data and supports the hypothesis that FLASH therapy, under specific dosing conditions, can induce an effective systemic antitumor response.

More precisely, the progression of the metastasis depends on two independent parameters, which describe respectively the local initial conditions (of the distant site) and the initial size of the metastasis at the onset of the systemic immune effect, induced by the FLASH therapy applied to the primary tumor. A possible correlation between them implies a metastasis size dependence of the immune response. Moreover, a time delay in immune response between the primary tumor and its metastasis suggests a dynamics not reducible to a one-compartment pharmacokinetic model.

The proposed approach relies on two key components: 1) a previous analysis of in vivo data on FLASH therapy, which supports a phenomenological parametrization of its effects (20, 22), and 2) a computational framework that models the immune response activated by radiotherapy. See refs (27–30). for the mathematical model of the synergy between radiotherapy and the immune system.

In this section, we clarify how measurements of the volume changes in the primary tumor and in the metastasis can provide a quantitative indication of the systemic immune response activated by FLASH therapy.

Indeed, the computational method of the previous section permits to evaluate the time evolution of the ratio V(t)/V(t ₁) for the primary tumor and metastasis, from time t ₁.

The parameters of the primary tumor evolution, for 25 Gy of FLASH therapy (20), have been determined in ref (22). (V ^{p (0} ) = 1): ln ( V ∞ p / V p ( 0 ) ) = 2.83, k^p = 0.0421 day⁻¹, c ₀ = 0.056, c ₁ = 0.166, c ₂ = 0.158 day⁻¹, and c_f = −0.00043 day⁻². The regrowth starts at t ₁ = 5 days after FLASH therapy and V^p (t ₁)/V^p (0) = 0.76.

For no delay, the results for the different values of λ for fixed y ₀ = 1, based on Equation 11, are depicted in Figure 3 . Also, a positive λ = 1.4 is reported to show its effect. The reference value is V^M (t ₁) = 1. The abscopal effect strongly depends on λ.

The analogous results as a function of y ₀, for λ = 1, are shown in Figure 4 . As previously discussed, y ₀ > 1 triggers the abscopal effect, enhancing the late time effect.

The two parameters may be correlated. In fact, although y ₀ is time independent, it can depend on the value of the ratio V ∞ M / V M ( t 1 ) , i.e., on λ.

The correlation between the two parameters has significant and testable implications. Indeed, the abscopal effect appears to be dependent on tumor size (37). In models with large distal tumors, treatment of the primary tumor produced a strong abscopal response, significantly inhibiting both primary and distant growth. In contrast, in models with small distal tumors, the abscopal effect was minimal, with a little difference observed in combination with irradiation. This difference may reflect variations in the strength of the immune response required for residual inhibition or size-dependent differences in immune cell infiltration.

It is useful to illustrate them with an example.

Let us suppose that, for the untreated tumor, the ratio V ∞ p / V p ( t 1 ) and k are experimentally determined by fitting the growth curve. If we consider distant metastases located at the same site, it is reasonable to assume that the local carrying capacity,   V ∞ M , is the same for both larger and smaller metastases, i.e., those with larger or smaller volume V M ( t 1 ) . Under these conditions, λ is inversely proportional to V M ( t 1 ) , and therefore, if the two parameters are not correlated, the abscopal effect is more pronounced for the larger metastasis, as shown in Figure 5 for λ = 0.7 and λ = 0.875, corresponding to a volume reduction of 20% at t = t ₁, for two independent values of y ₀. The abscopal effect is stronger for larger metastases (the continuous line is below the dotted line).

This result is reversed (the abscopal effect is greater for smaller metastases) if the two parameters are correlated, as shown in Figure 6 . Therefore, if it is confirmed that in the case of FLASH therapy, the abscopal effect is more pronounced for larger metastases, it would imply that the size of the metastasis is correlated with the effects of the immune system (37).

Notice that Figures 3 , 4 , and 6 also present the variability of the results, highlighting the sensitivity of the model to modifications in y ₀ and λ. Regarding the parameters c ₀ ,c ₁,c ₂, and c_f , their determination from the data in ref (20). is fairly accurate [see Table 2 in ref (22).], and a variation of one standard deviation results in effects of approximately 10% in the previous figures.

A time delay, τ, in the immune response between the immune effect in the primary tumor and the metastatic site implies that Y(t) = 0 for t < t ₁ + τ and Y(t) = y ₀ for t ≥ t ₁ + τ. In this scenario, the primary tumor begins to regrow at t = t ₁, while the growth rate of the metastasis is altered at t = t ₁ + τ due to the delayed onset of the immune response. Therefore, it is more useful to compare the metastasis and primary tumor evolutions after t = t ₁ using the previous formula (Equation 11) with t ₁ → t ₁ + τ for the metastasis.

The effect of the delay is depicted in Figure 7 for different values of λ and correlated y ₀ = 2λ. Notice that the conventional value V_M (t ₁) = 1 for τ = 0 implies that V_M (t ₁) > 1 for τ ≠ 0 .

This work aimed to propose a method to determine to what extent FLASH therapy induces a systemic immune response by monitoring the abscopal effect.

The computational approach for quantitative analysis is based on the following steps, which are informed by experimental data and guided by specific assumptions:

Steps (1) and (2) rely on the availability of experimental data. Steps (3) through (5) are model assumptions. The quantitative analysis requires the determination of only two to three additional parameters, which can be easily extracted using standard software tools applied to the metastatic growth data. The correlation between λ and y ₀ implies a direct connection between the abscopal effect and the size of the metastasis at the onset of the immune response due to FLASH therapy. The delay time gives direct information on the one- or multicompartment pharmacokinetic models.

The proposed computational method is based on the Gompertz growth law (point 1); however, any type of two-parameter sigmoidal curve can be used, such as the logistic curve. The Gompertz law provides a better fit to the available data [see ref (26)]. For completeness, the formulas applicable to the logistic growth curve are reported in Supplementary Material S3 .

Regarding point 4, the immune response to distant metastases may vary over time rather than remaining constant. However, this possibility can be readily incorporated into the computational model if a single value of y ₀ does not adequately fit the metastatic evolution data. This can be achieved by introducing a specific time-dependent form of Y(t) into Equation 6. For instance, an exponential decrease of the immune response over time is presented in Supplementary Material S2 .

The computational model relies on a set of parameters that must be specified for the given tumor phenotype in order to quantitatively assess the abscopal effect induced by radiotherapy. Therefore, the experimental procedure requires that the primary tumor is first inoculated in the test animals, followed by injection of the same tumor cells at a distant site to allow metastatic development. FLASH therapy is then applied to the primary tumor, and the evolution of both primary and metastatic tumor volumes is subsequently monitored. Initially, it is essential to fit the primary tumor growth data from two control groups—one untreated and one treated with FLASH therapy alone, without metastasis. This fitting procedure allows for the determination of the parameters k, V ∞ p , and c ₀, c ₁, c ₂, c ₃, and c_f for the tumor phenotype. Subsequent measurements of metastatic volume changes enable the determination of the two parameters λ, y ₀ (assuming no delay). However, by measuring both the primary and metastatic tumor volumes at the onset of the regrowth phase, the parameter λ becomes fixed—since the carrying capacity for the primary tumor has already been determined—and the fitting procedure then depends only on a single parameter y ₀. The delay time τ is an additional parameter that can be introduced during data fitting to test whether it plays a significant role. Its relevance can be assessed by evaluating the improvement in the quality of the fit (χ ² per degree of freedom) when the delay is included. As discussed, the assumption in point (4) can be relaxed at the cost of introducing additional parameters.

The difference between tumor phenotypes is reflected in the specific values of the parameters in the function F(t) (regarding the primary tumor) and in the values of the additional parameters y ₀ and λ, which characterize the induced immunity response on the metastases.

In conclusion, we have proposed a quantitative method to evaluate the abscopal effect associated with FLASH therapy. In our view, this type of treatment should trigger a systemic immune response. Otherwise, it would be difficult to explain why, following the initial sharp reduction in tumor cells, a regrowth phase occurs with a specific rate significantly lower than in the untreated case—and, more importantly, why a subsequent tumor regression takes place long after the treatment has been administered.

The computational method, while open to improvement in several respects, offers a quantitative estimate of deciphering the activated immune response.