Mathematical Analysis of FLASH Effect Models Based on Theoretical Hypotheses

where $p(t)$ is the concentration of oxygen at the time point $t$; $Recovery(t,p)$ is the term related to oxygen recovery in the tissue; $ROC(t,p)$ is the radiolytic oxygen consumption rate. The differential equation has an initial condition as listed in Equation (4).

where $p_{0}$ represents the initial oxygen concentration of the tissue before irradiation, which is influenced by factors such as the oxygen concentration in the capillaries, the density of capillaries, the rate of oxygen diffusion, the rate of metabolic oxygen consumption, and the distance from the tissue point to the nearest capillary.

Previous studies indicate that the process of oxygen recovery is mainly attributed to oxygen diffusion. The temporal profile of oxygen concentration approximately follows an exponential decay pattern (Hu et al.; 2022; Zhu et al.; 2021; Pratx and Kapp; 2019). The term of oxygen recovery in Equation (3) can be described by Equation (5).

where $\lambda$ is the characteristic constant for oxygen diffusion, determined by factors such as the density of microvessels in the tissue, the diffusivity of the medium and the tortuosity of cells. The parameter $\lambda$ encapsulates the spatial aspects of oxygen concentration variations in the model in an implicit manner.

The radiolytic oxygen consumption is mainly attributed to the reactions between radiation-induced radicals and oxygen (Wardman; 2020, 2021, 2022). In many related studies, the radiolytic oxygen consumption rate ($\mu$M/Gy) is often treated as a constant (Pratx and Kapp; 2019; Zhu et al.; 2021). However, this assumption is invalid in the cells where the oxygen concentration is low. Some studies set the oxygen consumption rate proportional to the oxygen consumption (Petersson et al.; 2020), which also cannot reflect the real condition when the oxygen concentration is high. Experiments have observed different radiolytic oxygen consumption rates under different initial oxygen concentrations(El Khatib et al.; 2023; Jansen et al.; 2021; Slyke et al.; 2022). The difference can be attributed to two effects: reactions between radiation-induced radicals and competition between oxygen and other cellular compounds (for example, antioxidants) reacting with radicals(Wardman; 2016). The short lifetime of radicals (e.g., ·OH, ·H, $e_{aq}$) resulting from water radiolysis, typically in the nanosecond range (Roots and Okada; 1975), does not align with the microsecond timescale of radiolytic oxygen consumption. Reactions among these radicals cannot play the primary role in the process. Thus, it is derived that the competition between oxygen and cellular compounds is the main factor influencing radiolytic oxygen consumption. The amount of radicals is proportional to the radiation dose.

Based on the above analysis, we use a fractional formula to describe the competition between oxygen and other cellular compounds as the formula used by Zou et al.(Zou et al.; 2022), which is the solution to the equation describing the competition between oxygen and antioxidants to react with radiation-induced radicals. This formula also corresponds to the results of radiolytic oxygen consumption vs. dose (Jansen et al.; 2021; Slyke et al.; 2022). The term of radiolytic oxygen consumption can be described by Equation (6).

where $g_{1}$ is the yield of radiolytic oxygen consumption, which is related to the radiation-induced radicals; $\dot{D}(t)$ is the dose rate at the time point $t$; $A$ represents the equivalent concentration of other cellular compounds competing with oxygen to react with radiation-induced radicals.

Previous studies have shown that the average dose rate during irradiation plays a crucial role in the FLASH effect (Hu et al.; 2022; Karsch et al.; 2022). In this study, we simplify the equation by assuming a constant dose rate throughout the irradiation, represented by Equation (7). This simplification proved valid for beams with much shorter pulse intervals than characteristic time of oxygen recovery ranging from hundreds of milliseconds to several seconds (Hu et al.; 2022). The requirement of the irradiation time structure well applies to the two types of beams commonly used in the FLASH experiments currently, i.e., linac electron beams with the pulse interval of $\sim$5 ms, and cyclotron protons with continuous wave mode. On the other hand, if the pulse interval is long enough, the whole therapy can be treated as the combination of several independent FLASH irradiations. It should be noted that if the pulse interval is neither too long nor too short, or the pulse interval is short but the pulse intensity is also changing, the analytical method is not suitable and numerical solutions are preferable.

where $D$ is the total dose and $T$ is the total irradiation time.

With the above assumptions and simplifications, Equation (3) is converted into Equation (8).

The relative amount of damage induced by the FLASH irradiation is calculated by the dose-averaged $OER$ over the whole irradiation period. For a constant dose rate scenario, damage can be calculated by Equation (10).

As for the CONV condition, because of the low dose rate delivery, the oxygen concentration remains almost unchanged during the irradiation. Thus, the biological effect can be represented as Equation (11),

To estimate the maximal FLASH effect, we calculate the limit of the function defined by Equation (10) by setting $T\to 0$.

Within the extremely short irradiation duration, the oxygen recovery can be ignored. The relationship between oxygen concentration at the end of irradiation and the total dose can be derived by solving Equation (12) with initial condition Equation (4).

By setting $t=T$ in the solution of Equation (12), we can calculate the oxygen concentration at the end of the irradiation, $p_{T}$, to estimate the maximal change of oxygen concentration during FLASH irradiation.

According to the same reason described in the model of oxygen depletion, we consider the radiation delivered at a constant dose rate. Thus, Equation (13) can be divided into two stages: Equation (14) during irradiation and Equation (15) after irradiation.

where $R_{T}$ is the concentration of peroxyl radicals at the end of irradiation, which can be calculated by the solution of Equation (14).

For CONV irradiation, we calculate the limit of the damage function defined by Equation (17) to estimate the radiobiological effect of CONV irradiation.

Then we calculated the integral in Equation (10) by changing the integration variable from $t$ to $p$ as shown in Equation (20u).

The result of the integral is shown in Equation (20v).

In the calculation with computer programs like MATLAB, we find that for the conditions of low total dose or long irradiation time, effective results cannot be obtained through the forms of Equation (19) and (20v) due to the limited numerical precision of the computer. Thus, we use some tricks to conduct the calculation as follows.

For low total dose, we use the first-order Taylor expansion of Equation (19) to calculate $p(t)$, which is shown in Equation (20w).

With a similar mathematical trick with Equation (20w), the biological effect of the above-mentioned condition ($p_{T}\to p_{0}$) can be calculated by Equation (20x).

For long irradiation time, the oxygen recovery and radiolytic oxygen consumption reach equilibration, i.e., $\lambda(p_{T}-p_{0})\approx g_{1}Dp_{T}/[T(A+p_{T})]$. The term, $p_{T}$, can be calculated by Equation (20y).

The damage can be calculated by Equation (20z) according to Equation (20y), Equation (20v) and Equation (19).

Radiolytic oxygen consumption is one of the key factors in the oxygen depletion hypothesis. By setting $t=T$ in the Equation (19), we can calculate the radiolytic oxygen consumption (net change of oxygen concentration during the irradiation, which includes the effect of oxygen recovery) induced by a given dose ($D$) and irradiation time ($T$). Here, we define $\Delta p=p_{0}-p_{T}$ to represent the radiolytic oxygen consumption. The maximal radiolytic oxygen consumption, $\Delta p_{max}$, is calculated by taking the limit for $T\to 0$.

First, to show the impact of irradiation time on the radiolytic oxygen consumption, by setting the initial oxygen concentration ($p_{0}$) as 10 $\mu$M (7 mmHg), we calculate a series of radiolytic oxygen consumption ($\Delta p$) for a range of irradiation times spanning milliseconds to several seconds, under total doses ranging from 5 Gy to 60 Gy. The resulting curves of normalized oxygen consumption, i.e., $\Delta p/\Delta p_{max}$, are plotted against irradiation time for different dose deliveries in figure 1.

The results of the mechanism hypothesis-based calculations highly reflect the core principle of the oxygen depletion hypothesis. The curves exhibit an identical inverse “S” shape, indicating that shorter irradiation time results in greater oxygen consumption due to the limited rate of oxygen recovery. Additionally, the curves are nearly overlapping across different doses. This suggests that the temporal characteristic of radiolytic oxygen consumption is largely dose-independent. The oxygen recovery process determines how the irradiation time impacts the radiolytic oxygen consumption, which is almost independent of dose.

Furthermore, we study the influences of initial oxygen concentration ($p_{0}$) and dose ($D$) on the maximal oxygen consumption, $\Delta p_{max}$. The implicit solution Equation (12) is listed in Equation (20aa).

With different initial oxygen concentrations ranging from 1 to 100 $\mu$M (0.7-70 mmHg), the values of $\Delta p_{max}$ against total doses from 0 to 100 Gy are calculated based on the Equation (20aa). As shown in figure 2, the resulting curves indicate that $\Delta p_{max}$ is no longer proportional to the dose magnitude when the dose has already exceeded the amount required to make the cell hypoxic. The saturation effect occurs significantly for low values of initial oxygen concentration. For the simulation parameters used in this study, it can be seen in the curves with the initial oxygen concentration lower than 20 $\mu$M (14 mmHg).

Additionally, to show the influence of the competing reactions for the radiolytic oxygen reaction, we change the equivalent concentration of other cellular compounds ($A$) from 1.0 to 20.0 $\mu$M (0.7-14 mmHg) and calculate the corresponding $\Delta p_{max}$ induced by 20 Gy FLASH irradiation under different initial oxygen concentrations.

As shown in figure 3, the competition between oxygen and other cellular compounds influences radiolytic oxygen consumption greatly when the initial oxygen concentration is low ($<20\mu$M for the parameters in this study). The radiolytic oxygen consumption cannot be regarded as a constant for different conditions of the cellular composition.

The impact of irradiation time ($T$) on the FLASH effect ($DRF$) can be derived as Equation (20ab) by Equation (20v) (damage of FLASH) and Equation (11) (damage of CONV). Data points of $DRF$s vs. $T$ for an example dose setup are shown in figure 4. The total dose is set to 20 Gy, and the initial oxygen concentration is set as $p_{0}=10\mu$M.

The relationship of $DRF$ vs. $T$ can be approximated by a relatively simple formula (Equation (20ac)). The fitting curve of this formula is plotted in figure 4 with fitting parameters listed in table 3.

where $DRF_{min}$ is the minimal $DRF$ for a given dose with extremely high dose rate irradiation, i.e., $T\to 0$; $\xi_{1}$ and $\xi_{2}$ are dose-independent constants, corresponding to the dose-independence of the time characteristic of radiolytic oxygen consumption shown in figure 1.

As for the choice of initial oxygen concentration, $10\mu$M, in this example, though it only represents the situation for hypoxic tissue and tumors (Vaupel et al.; 1991), we present the resulting curve because the change of biological effect is numerically larger and the trend can be observed more obviously than that for normal tissue situation. The average oxygen concentration of normal tissue is $\sim 70\mu$M (50 mmHg) (Vaupel et al.; 1991), under which the change of biological effect is too slight to show in the figure.

Then, the minimal damage of FLASH irradiation by a given dose is estimated by calculating the limit of Equation (20v) when $T\to 0$. The result is shown in Equation (20ad)

The $DRF_{min}$ is given with Equation (20ad) deciding the $Damage$ induced by the same dose of CONV irradiation as shown in Equation (20ae).

The $DRF_{min}$ for different doses and initial oxygen concentrations are shown in figure 5. The curves indicate that the FLASH effect predicted by the oxygen depletion hypothesis is influenced by the initial concentration greatly. For normal oxygenated tissues with the median oxygen concentration of $70\mu$M (50 mmHg) (Vaupel et al.; 1991), the biological effect changes slightly. However, in hypoxic tissues or hypoxic parts of tumors where the typical oxygen concentration ranges from 0 to $10\mu$M (7 mmHg)(Vaupel et al.; 1991; Menon et al.; 2003), the biological effect alters greatly, which contradicts the almost unchanged tumor control shown in FLASH effect experiments. Our results indicate that the oxygen depletion hypothesis cannot feasibly explain the FLASH effect that spares the normal tissue but keeps the equivalent tumor control.

The second part is the solution of the equation which describes the reactions after irradiation ($t\geq T$). The solution is shown in Equation (20ah).

where $C_{2}$ is a constant determined by the concentration of peroxyl radical at the end of irradiation ($R(t)$), which can be calculated using Equation (20af).

Then, we calculate the integral of the solutions of these two parts based on Equation (16). For the first part (integral of Equation (20af)), the integral interval is $[0,T)$, and the result is shown in Equation (20aj).

For the second part (integral of Equation (20ah)), the integral interval is $[T,+\infty)$. The result is shown in Equation (20ak).

The sum of these two parts (AUC[ROO·]) represents the damage related to peroxyl radicals. The magnitude of AUC[ROO·] against the total irradiation time, $T$, is shown in figure 6. In this figure, each data point is an individual simulation, not the cumulative effect of results of different irradiation times.

The result indicates the time characteristic of the FLASH effect. The damage to normal tissue decreases with the decline of irradiation time, which forms an “S” shape curve. The steepest time interval of the curves is around hundreds of milliseconds corresponding to (Montay-Gruel et al.; 2017). The “S” shape curve predicted by the hypothesis model indicates that the radical recombination-antioxidants hypothesis can explain the time characteristic of FLASH to a certain extent. The characteristic time is mainly determined by the reaction between peroxyl radicals and antioxidants, which is represented by the parameter $k_{1}$ in the model. As shown in the three curves with different $k_{1}$, the parameter can also greatly influence the difference in radiobiological effect between FLASH and CONV irradiation with the same dose delivery. The simulation results align with the radical recombination-antioxidants hypothesis, which deems that peroxyl radical-related reactions dominate the FLASH effect. It underscores that antioxidants compete to react with peroxyl radicals, which determines the time characteristic of the FLASH effect and the maximal difference between FLASH and CONV for a specific dose. This is due to the fact that the rates of these reactions are governed by the concentrations of the reacting species. The hypothesis infers that tissue with different antioxidant levels may exhibit different time characteristics of the FLASH effect (Hu et al.; 2022). The influence of $k_{1}$ on the time characteristic and difference between FLASH and CONV reflects this main point of the hypothesis.

However, the formula is too complex to explicitly reflect the relationship. We find a simple form in Equation (20am) as an excellent approximation ($R^{2}>0.99$) of $DRF(T)$ in a wide range of parameters (we tested it by setting dose to 0.1-100 Gy and setting $k_{1}$ to 0.01-100 s^-1).

where $DRF_{min}$ is the minimal value of $DRF$ for a given dose, which is equal to the value calculated by Equation (20anaq); $T_{m}$ is a characteristic time, which is determined by $k_{1}$ and $D$. We also find that the formulas of $T_{m}(k_{1})$ and $T_{m}(D)$ have concise approximate forms listed in Equation (20ana) and Equation (20anb). The formula of $T_{m}(k_{1})$ is obtained by excluding the data points where $k_{1}<1$, because the data points of $k_{1}>1$ can be approximated by a linear relationship on a log-log scale. If the radical recombination-antioxidants were the main mechanism of the FLASH effect, the value of $k_{1}$ should be higher than 1.0 s^-1 because it should explain the time characteristic of the FLASH effect observed by experiments (Montay-Gruel et al.; 2017). These formulas can be used as empirical formulas to predict the FLASH effect in practical application.

where $\varepsilon_{1},\varepsilon_{2},\varepsilon_{3},\vartheta_{1},\vartheta_{2}$ are parameters that can obtained by fitting. The data points used in fitting and fitted curves of $T_{m}(D)$ and $T_{m}(k_{1})$ in figure 7. The results of curve fitting are listed in table 4.

For clinical application, the $T_{m}$ for a given dose can be obtained by fitting the data points of $(\log(T),DRF(D,T))$. The $T_{m}(D)$ can be obtained by fitting the data points from a series of doses.

Based on the methods list in section 2.2.3, we calculated the minimal damage induced by extremely short time ($T\to 0$) irradiation (Equation (20anao)) and CONV (Equation (20anap)) irradiation to obtain the minimal value of $DRF$, i.e., $DRF_{min}$, to evaluate the maximal change of biological effect by a given dose of FLASH irradiation.

Then the $DRF_{min}$ can be calculated and the result is shown in Equation (20anaq).

We obtain a concise formula and it can be utilized to fit the experimental data conveniently. According to this formula, Figure 8 shows the $DRF_{min}$ irradiated by different doses under cells with different levels of antioxidants, $k_{1}$.

The curves indicate that the antioxidants in the cell greatly influence the change of biological effect attributed to the FLASH irradiation. This also corresponds to the hypothesis’ explanation of the non-protective effect in tumors, which posits that the elevated antioxidant levels within tumors inhibit radical recombination, thereby narrowing the difference between the FLASH and CONV modalities. The value of $k_{1}$ for normal tissue is estimated to be 8.5 s^-1 (Hu et al.; 2023) whereas in tumors, this value is several times higher. Moreover, the curves show that the difference between FLASH and CONV irradiation increases with the increase of the total dose when the total dose is low. Then the $DRF_{min}$ changes relatively slightly with the increase of the total dose when the total dose is high enough.

For the oxygen depletion hypothesis, oxygen is the main point and the factors related to oxygen should remain the same. The density of capillaries and the types of cells in the tissue mainly determine the factors related to oxygen consumption and recovery, which are the main considerations in data filtration. According to the oxygen depletion hypothesis, the criteria for the same group experiments are listed below.

1. 1.

   The tissues should have the same density of capillaries and the same tortuosity of cells for the same time characteristic of oxygen recovery, i.e., $\lambda$.

2. 2.

   The tissues should have the same metabolic oxygen consumption rate and the same blood supplement for the same initial oxygen concentration, i.e., $p_{0}$.

3. 3.

   The tissues should have the same level of other cellular compounds, i.e., $A$ in the model, competing to react with radiation-induced radicals alongside oxygen, among which cellular antioxidants are the most significant.

4. 4.

   The OER curves of the tissues should be similar, ensuring the same level of $m$ and $K$.

5. 5.

   The types of irradiation should be similar for the same yields of radicals, i.e., $g_{1}$.

The radical recombination-antioxidants hypothesis mainly focuses on reactions of radicals and the level of antioxidants, not the oxygen in the tissue. According to the radical recombination-antioxidants hypothesis, the criteria for the same group experiments are listed below.

1. 1.

   The tissues should have the same level of antioxidants, i.e., $k_{1}$.

2. 2.

   The responses of cells to irradiation should be similar for the same level of $f_{1}$ and $f_{2}$.

3. 3.

   The types of irradiation should be similar for the same yields of radicals. It should be noted that the tissues should not be extremely hypoxic because low oxygen concentration could greatly compress the generation of peroxyl radicals, i.e., $g_{2}$.

The experimental design or the retrospective study should consider the above criteria to avoid the influences on the conclusion brought by these factors.

To overcome the limitation, we propose a series of experiments using fractionated FLASH dose setups. In these experiments, the definition of fraction differs from the general definition in the field of radiotherapy. The duration between two consecutive fractions is several minutes rather than the common practice of one day. The total dose of FLASH irradiation is aimed to cover the high-dose interval where NTCP changes greatly, which is different from that of CONV irradiation because of the FLASH effect. However, for the limited granularity of fractionation dose, simple fractionated irradiation may not be able to achieve the total dose target. Thus, a hybrid irradiation strategy is proposed to supplement the limitation. The detailed description of the systematic experiments is listed below.

1. 1.

   CONV group. For CONV irradiation, the total dose is set to reach the interval where NTCP changes greatly. Then an endpoint is chosen to obtain the reference dose, $D_{CONV}$.

2. 2.

   Fractionated FLASH dose setup. For FLASH irradiation, set a series of doses, $D_{FLASH,f}$, e.g., 1, 2, 5, 10,… Gy. For each setup of fractionated dose FLASH irradiation, we can define a virtual iso-effective dose, $D_{FLASH,ISO}(D_{FLASH,f})$, which induces the same radiobiological effect with the $D_{CONV}$ dose of CONV irradiation. The virtual iso-effective dose, $D_{FLASH,ISO}$, often cannot be obtained by fractionated FLASH irradiation solely because of the limited granularity of the fractionated dose.

3. 3.

   Hybrid irradiation strategy. For each dose setup of fractionated FLASH irradiation, to start with, $D_{FLASH,f}$ is delivered at FLASH dose rate for $n$ times to the experimental group with the duration of several minutes between two fractions, where $n$ is the maximal integer for $n\cdot D_{FLASH,f}\leq D_{FLASH,ISO}$. Then, the group is irradiated by a residual CONV dose, $D_{CONV,residual}$, to obtain the iso-effect dose for the fractionation dose setup, $D_{hybrid,ISO}(D_{FLASH,f})$, which induces the same radiobiological effect with $D_{CONV}$ CONV irradiation. The $D_{hybrid,ISO}(D_{FLASH,f})$ is the sum of $n\cdot D_{FLASH,f}$ and the residual CONV dose.

   With the systematic experiments, the DRF of the fractionation dose for FLASH effect can be calculated by Equation (20anar).

These experiments enable the acquisition of $DRF$ versus dose and irradiation time curves across a wide range. The established dataset is expected to provide a fundamental base for advancing clinical FLASH radiotherapy.

It should be noted that the two mathematical models have different features, and the implementation pipelines are somehow different.

For the oxygen depletion hypothesis, the model parameters are generally accessible through measurements. They can be obtained by the steps listed below and then used to calculate the $DRF$ value by the analytical Equation (20ab) for an exact dose,$D$, within the irradiation duration, $T$, and also the $DRF_{m}in$ value by Equation (20ae) for the extreme condition as $T\to 0$.

1. 1.

   Obtain the OER curve of the tissue to get the parameters, $m$ and $K$, within the curve.

2. 2.

   Measure the maximal radiolytic oxygen consumption irradiated by a series of given doses to obtain the yield of radiolytic oxygen consumption, $g_{1}$, and the equivalent concentration of cellular compounds, $A$, according to Equation (20aa).

3. 3.

   Measure the oxygen recovery curve to obtain the feature constant, $\lambda$, or fit the data points with the same dose but a variety of irradiation times to obtain $\lambda$, fixing the parameters from the above two steps.

With the parameters measured above, the FLASH effect can be predicted with the following steps in the clinical application.

1. 1.

   Calculate the $DRF_{min}$ vs. dose curve based on Equation (20ae).

2. 2.

   Set a dose and calculate a series of $DRF$ for different irradiation durations according to Equation (20ab).

3. 3.

   With the calculation results,the parameters, $\xi_{1}$ and $\xi_{2}$, in the approximate formula Equation (20ac) can be fitted. It should be noted that these two parameters are independent to doses and can be used for all doses because the variation of net oxygen concentration by the irradiation is independent to doses, as shown in figure 1.

4. 4.

   For a specific dose, $D$, one can use the approximate formula, Equation (20ac), and the $DRF_{min}$ calculated from step 1 to calculate the $DRF$ vs. irradiation time curve and predict the FLASH effect for clinical use.

According to the radical recombination-antioxidants hypothesis, we can predict the FLASH effect solely based on the experimental data points of $DRF$ with different doses and irradiation times, other than the analytical calculations of $DRF$ and $DRF_{min}$ as above, because of the concise form of the approximate formula for $DRF$ (Equation (20am)), and the inaccessibility of the model parameters, especially for the exact values of $f_{1}$ and $f_{2}$. The prediction steps are listed below.

1. 1.

   Divide the experimental data points, $DRFs$, into several groups according to the dose.

2. 2.

   For the data group with the same dose, fit the parameters $T_{m}$ and $DRF_{min}$ according to the form of Equation (20am) using the experimental data points ($T$, $DRF$) classified in the first step.

3. 3.

   Fit the curve of $DRF_{min}$ vs. dose as Equation (20anaq) using the data points ($D$, $DRF_{min}$) from the second step.

4. 4.

   Fit the curve of $T_{m}$ vs. dose to obtain parameters in Equation (20ana) using the data points ($D$, $T_{m}$) from the second step.

5. 5.

   Predict the FLASH effect based on Equation (20anaq), Equation (20ana) and Equation (20am) for a given dose and irradiation time.