In this model, particle tracks are assumed to be perfectly straight and homogeneous along the beam axis, and can thus be described in the 2-dimensional plane orthogonal to the beam axis with the position vector r ⃗ = ( x , y ) . We consider a system of N particle tracks distributed across a target. The spatio-temporal distribution of the RS of an arbitrary track at position r i ⃗ and arrival time t _i is defined by a probability density function (PDF) c i ( r ⃗ , t ) . We do not differentiate between different RS and instead describe the net dynamics of all RS of the track together. The track diffuses, interacts with its environment, and potentially interacts with other tracks. Thus, the evolution of a track is governed by a reaction-diffusion equation comprising a radially symmetric diffusion term with diffusion coefficient α, an environmental scavenging term with rate constant k _e , and a track-track interaction term with rate constant k _r : ∂ ∂ t c i r ⃗ , t = α ∇ 2 c i r ⃗ , t − k e c i r ⃗ , t − k r ∑ i , j N c i r ⃗ , t c j r ⃗ , t . (1) Track-track interactions are assumed to have a net reductive effect on the total number of RS, since they mostly comprise primary radical combination into secondary molecules. The track-track interaction term can be decomposed into an intratrack interaction term (where i = j) and an intertrack interaction term (where i ≠ j). Since the model follows track dynamics up till the point at which intertrack interactions become significant, it is assumed that intertrack interactions negligibly affect the evolution dynamics of individual tracks up until that point. Thus, the intertrack interaction term is neglected, and intertrack interaction is quantified separately in Section 2.1.3.

The intratrack interaction term and the environmental scavenging term, both of which contribute to a decline in the number of RS in the track, are combined into one decay term with effective decay constant k _s . This is done in order to keep Eq. 1 linear and thus allow for an analytical solution. Thus, the simplified diffusion-decay reaction defining individual tracks in this model is ∂ ∂ t c i r ⃗ , t = α ∇ 2 c i r ⃗ , t − k s c i r ⃗ , t . (2)

Setting c ₀ to be the initial total number of RS in the track, the solution to this partial differential equation (see Supplementary Appendix SA1) yields the PDF of an arbitrary track centered at r i ⃗ for t > t _i , a normal distribution about the arrival point of the track that broadens and decays with time: c i r ⃗ , t = c 0 ⋅ e − | r ⃗ − r ⃗ i | 2 4 α t − t i + τ 0 − k s t − t i 4 π α t − t i + τ 0 . (3) Here, we applied a temporal shift parameter τ ₀, such that the track has a finite width upon its arrival t = t _i , which represents the spatial variance along the primary particle’s path.

The beam is assumed to have a constant dose rate over the exposure time T, and thus the PDF of the arrival time t _i of the ith track at time t is P t t i = 1 min t , T for 0 ≤ t i ≤ min t , T 0 elsewhere . (4) The expected total number of tracks on the target over time N(t) can be expressed with the total energy deposited in the target E and the average energy deposited per track E _t for a target of density ρ, depth z (measured along the beam axis), and cross-sectional area A. N t = E E t ⋅ min t , T T (5) = D ⋅ A ⋅ ρ ⋅ z E t ⋅ min t , T T (6) = D ⋅ A ⋅ ρ L ⋅ min t , T T . (7) The LET (L = E _t /z) is introduced in the last step and will depend on the radiation quality.

Assuming a spatially homogeneous beam across the target area A, the expected number of tracks within any circular area of radius R can be expressed as a fraction of the total number of tracks on the target: N R t = N t ⋅ π R 2 A (8) The PDF of the displacement s between one arbitrary track and any other track within a radius R is given by the ratio between the area element 2πs and the total area πR ²: P s s = 2 s R 2 , s ≤ R . (9)

The quantification of track-track interaction is a critical component of this model. A measure is needed that accounts for the spatial variations of different tracks, and reflects a physical quantity relevant to the radiation chemistry at hand. To this end, first the interaction rate ω is first defined, which yields the rate of change of the quantity of species at time t due to the interaction of two tracks assuming second order reactions with reaction rate k _r . ω 1,2 t = ∫ k r ⋅ c 1 r ⃗ , t ⋅ c 2 r ⃗ , t d V (10) = k r ⋅ c 0 2 ⋅ e − | r 1 ⃗ − r 2 ⃗ | 2 8 α 2 t − t 1 − t 2 + 2 τ 0 − k s 2 t − t 1 − t 2 8 π α 2 t − t 1 − t 2 + 2 τ 0 . (11) Different reaction rates will be specified for intertrack interactions k _r,inter and intratrack interaction k _r,intra .

Integrating the interaction rate in Eq. 11 over time yields the total track-track conversion I, which represents the net change in the quantity of species due to the interaction over the relevant time period: I 1,2 t = ∫ 0 t ω 1,2 t ′ d t ′ . (12) I is the measure of interest in this model, as it represents an empirically measurable quantity.

The expected interaction rate of an arbitrary track a with all neighboring tracks within an interaction volume of radius R can be found by multiplying the number of tracks inside that interaction volume N _R by the expected interaction rate between one arbitrary track and another within that volume: ∑ i ω a , i t = N R t ⋅ ⟨ ω 1,2 t ⟩ , (13)

The total intertrack conversion of all tracks in the target volume is then I inter t = N t 2 ⋅ ∫ 0 t lim R → ∞ N R t ′ ⋅ ⟨ ω 1,2 t ′ ⟩ d t ′ , (14) where R approaches infinity under the assumption that each track is effectively within an infinite isotropic volume. The factor 1 2 ensures that each intertrack interaction is not double-counted.

Evaluating this expression (see Supplementary Appendix SA2) yields I inter t = lim ϵ → 0 B 2 ⋅ t ⋅ ln t ϵ − Γ 0,2 k s t + 2 Γ 0 , k s t + Γ 0,2 k s ϵ − 2 Γ 0 , k s ϵ for t ≤ T I inter T + B 2 ⋅ e k s T − 1 2 2 k s ⋅ e − 2 k s T − e − 2 k s t for t > T , (15) where B = N t 2 A ⋅ k r , inter k s 2 ⋅ c 0 2 ⋅ 1 min t , T 2 , (16) and Γ(s, x) is the incomplete upper gamma function given by Γ s , x = ∫ x ∞ t s − 1 e − t d t . (17)

The total intratrack conversion across the entire target volume is simply a sum of each track’s self-interaction I intra t = ∑ i N t ∫ 0 t ω i , i t ′ d t ′ . (18) This expression can be simplified with a suitable approximation, removing the need for tedious summation over all tracks. A track’s outward diffusion means that the amount of intratrack conversion of a single track asymptotically converges over time to a maximum value. If the time it takes for the amount of conversion to almost fully converge is negligible relative to the total exposure time, the total intratrack conversion can be approximated as a linear increase over the exposure time towards the total maximum amount of intratrack conversion of all tracks on the target. I intra t ≈ N t ⋅ ∫ 0 ∞ ω i , i t d t (19) = N t ⋅ ∫ 0 ∞ k r , intra ⋅ c 0 2 ⋅ e − 2 k s t 8 π α t + τ 0 d t (20) = N t ⋅ k r , intra ⋅ c 0 2 ⋅ e 2 k s τ 0 8 π α ⋅ Γ 0 , 2 k s τ 0 . (21) If the convergence time is long relative to the total exposure time, this approximation will underestimate the amount of intratrack interaction at short time scales, but will converge on the correct value at longer time scales.

Although I _inter gives the absolute effect of intertrack interaction, we were interested in quantifying its relative effect on the total ongoing radiation chemistry, which normally consists of only intratrack interactions. We therefore calculated the relative change that intertrack interactions cause to the total radiation chemistry. Φ t = I inter t + I intra t I intra t − 1 (22) = I inter t I intra t . (23) Φ(t) thus remains 0 until intertrack interactions begin to affect the radiation chemistry; once Φ(t) has surpassed some critical value Φ _c , we may consider intertrack interactions to be significant. Φ(t) is the main measure of interest in this model.

Evaluating this measure with Eqs 7, 14, 18 yields Φ t = K ⋅ 2 π α ⋅ D ⋅ ρ L ⋅ e − 2 k s τ 0 k s 2 ⋅ Γ 0,2 k s τ 0 ⋅ 1 T ⋅ f 1 t , k s for t ≤ T f 1 T , k s + e k s T − 1 2 k s T e − 2 k s T − e − 2 k s t for t > T (24) where f 1 t ′ , k s = 2 ⋅ lim ϵ → 0 ln t ′ ϵ − Γ 0,2 k s t ′ + 2 Γ 0 , k s t ′ + Γ 0,2 k s ϵ − 2 Γ 0 , k s ϵ , (25) K = k _r,inter/k _r,intra is the ratio of reaction rates between intertrack and intratrack interactions, and Γ(s, x) is defined as in Eq. 17. Table 1 summarizes all variables and their definitions.

The primary particle fluence can be expressed using Eq. 7 as. F = N A (26) = D ⋅ ρ L . (27)

The value of Φ(t) for t ≫ T is of interest for analyzing the final ratio of intertrack-to-intratrack interactions after all tracks have decayed and no more interactions can occur. This value, using Eqs 24, 27, is given by. Φ t ≫ T = lim t → ∞ Φ t (28) = K ⋅ α ⋅ F ⋅ f 2 T , k s , τ 0 , (29) where f 2 T , k s , τ 0 = e − 2 k s τ 0 k s 2 ⋅ Γ 0,2 k s τ 0 ⋅ 2 π T ⋅ f 1 T , k s + e − k s T − 1 2 k s T . (30)

To analyze the spatio-temporal evolution of a single track, modulated by the model parameters α and τ ₀, the radial variance of the first 1 μm of a simulated track in a 10 μm diameter water volume is compared with its variance over time as predicted by the model: 2α(t + τ ₀). Only the first 1 μm of the simulated track is analyzed to avoid the effects of beam straggling deeper in the target, which would erroneously affect the analysis of the track’s variance due to RS diffusion.

To analyze the decay in number of RS caused by the combination of intratrack interactions and scavenging, modulated by the model parameter k _s , the number of RS of the first 20 μs of the simulated track in a 10 μm diameter water volume is compared with the model’s simplified decay term, e − k s t , where k _s is the effective track decay constant.

MC simulations of multiple particles and their interacting tracks were performed to compare with the model’s estimation of intertrack interaction. The impinging particles were distributed uniformly randomly across the face of the 1 μs diameter, 1 μs deep cylindrical target volume. Other than the stated differences, all other simulation parameters remained equal to those in the aforementioned simulations (Section 2.3.1). Two different doses, 12 Gy and 58 Gy, were simulated by adjusting the number of primary electrons to 393 and 1963, respectively. The CPU times for these simulations are reported in Table 3. The exposure time did not vary because within TOPAS n-Bio, all tracks arrive in the target simultaneously, yielding an effectively infinite dose rate. These simulations were repeated for all three scavenging modes.

Two otherwise identical simulations were performed: one in which tracks were simulated independently of one another and thus unable to interact, and one in which all tracks were simulated together, such that tracks were able to chemically interact. These simulation modes are referred to as ‘non-interacting’ and ‘interacting’, respectively. The ‘non-interacting’ simulations represent the radiation chemistry at the low-dose-rate limit, where interactions consist only of intratrack and scavenging interactions. The ‘interacting’ simulations add on the effects of intertrack interactions. The number of RS in each simulation type was recorded over time and compared.

MC simulations of many interacting tracks were used for comparison with the model’s estimation of influence of intertrack interaction. This was repeated for two doses, 12 Gy and 58 Gy, and all three scavenging modes. Figures 5A, B show the numbers of the different RS of interest during a simulation where tracks were not allowed to interact (solid lines) or allowed to interact (dotted lines). Results are shown here only for the 58 Gy simulations and the physoxic and cellular scavenging modes. The normoxic scavenging mode results tended to lie between the physoxic and cellular results (Supplementary Appendix SA3, Supplementary Figures S7A, B), whereas the 12 Gy simulation showed far smaller differences between the interacting and non-interacting simulations for all scavenging modes.

For the physoxic scavenging mode, intertrack interactions resulted in a decrease in the yields of ^•OH and e a q − , and a slight increase in the yield of H₂O₂, beginning at ≈ 10⁻⁷s (Figure 5A). Figure 5C reveals the individual reactions most responsible for these differential yields. Intertrack depletion of ^•OH, and e a q − is almost solely caused by their reaction with each other, whereas the self-recombination of ^•OH, and the consumption of e a q − by H₃O⁺, contribute only slightly. For the cellular scavenging mode, the quick depletion of ^•OH and e a q − due to high background scavenging is evident in Figure 5B. Here, intertrack interactions consist mostly of the reaction of H₃O⁺ with OH⁻ (Figure 5D) and are negligible with respect to the yields of the RS of interest.

To compare the simulation results with those of the model, the ratio of the sums of the RS of the two simulation types shown in Figures 5A, B are taken and depicted by the solid lines with shaded in errors in Figures 5E, F. This ratio is a measure of the influence of intertrack interactions; it yields the relative difference of the RS counts caused solely by intertrack interactions. At very short time scales, this ratio remains 1 as there is no difference between the simulation types. At longer time scales, once intertrack interactions begin to affect RS yields, this ratio deviates from 1. Unsurprisingly, the higher the dose and the lower the amount of scavenging, the larger the effect of intertrack interactions. Notably, the slight deviation from 1 in the cellular scavenging results is almost entirely due to the change in H₃O⁺ yields, as shown in Figure 5D. Results for the normoxic scavenging mode are not shown for brevity; as expected, they tended to fall between the results of the cellular and physoxic scavenging modes, remaining closer to the physoxic results (Supplementary Appendix SA3, Supplementary Figure S7C).

This ratio is directly compared with an analogous measure from the model, that is, 1 − Φ(t). The ratio aligns quite well with what 1 − Φ of the model actually represents based on Eq. 23. Φ is subtracted from 1 because intertrack interactions resulted in a net decrease in the number of the selected RS. The final free parameter of the model, the reaction rate ratio K, was adjusted to fit 1 − Φ to the simulation results. The variables used to calculate Φ(t) (see Eq. 24) were: T = 10⁻¹⁵s, L = 0.14 keV/μm, ρ = 997 kg/m³, D was adjusted to match the simulation dose, and K, α, τ ₀, and k _s were adjusted depending on the scavenging mode to the values in Table 4. The value of T was chosen to simulate an effectively infinite dose rate since all particles arrived simultaneously in the simulation; decreasing T to even lower values had negligible effects on the resulting values of Φ(t).

The model output was shown to align well with MC simulation data when fitted with the variable K = 0.001. We also performed a similar comparison, now with experimental results of H₂O₂ yields after electron irradiation of physoxic water. Kacem et al irradiated 4% O₂ water with a pulsed electron beam in 3 different setups, which we labeled A-C, and observed a decrease in the resultant G-value of H₂O₂ ( G H 2 O 2 ) with increasing dose rate, as shown in Table 5 [35]. By following the argument in Section 2.2, we assume that there are no interactions between the pulses, which were delivered at 100 Hz. Thus, the mean dose rate and number of pulses are irrelevant, and the only parameters that matter are the dose and dose rate within the pulse.

Φ was calculated for these pulse parameters by using the values for model parameters α, τ ₀, k _s , and K derived in the above sections for physoxic scavenging conditions, and its value is compared with the measured change in G H 2 O 2 relative to Setup A. A 17% decrease in G H 2 O 2 was observed in Setup B relative to Setup A, although the model predicts negligibly low intertrack interaction (Φ = 4.1 × 10⁻³ ≪ 0.17) at Setup B.

In the case the derived values for the model parameters are off, we examined what changes to these parameters could result in a fit to this experimental data; i.e., Φ = 0.17. We saw that either increasing K, increasing τ ₀, or decreasing k _s could yield agreement between the value of Φ and the experimentally observed changes. We are not as interested in matching the values at Setup C because we expect Φ to be overestimated as intertrack interaction increases, due to the weakly-interacting tracks assumption made in Section 2.1.1.

The pulse parameters of this study can also be visualized on the heatmap (Figure 6) introduced in the following section, which displays Φ over the pulse parameter space for cellular scavenging conditions. The heatmap of Φ for physoxic scavenging conditions looks similar, but the isovalue curve is shifted slightly down towards lower doses and to the right towards higher exposure times. The ‘corner’ of the white isovalue curve occurs in this case at 13 Gy and 2 μs for K = 0.003, and 0.3 Gy and 2 μs for K = 0.13.

To compare the predictions of the model with experimental in vivo data from the literature, we used a heatmap of the value of Φ over the dose and exposure time parameter space (Figure 6). Again, by following the argument in Section 2.2, we can consider this to be the pulse dose and pulse width of each pulse in a pulsed electron delivery. Results are shown using the values for model parameters α, τ ₀, and k _s derived in the above sections for cellular scavenging conditions, with the notable exception of K = 0.33, which was chosen as it best fits these experimental data. The white and blue regions of the heatmap indicate where the model predicts intertrack interaction to have a ≥ 5 % effect on the yields of reactive species. The pulse parameters from in vivo electron FLASH data are superimposed on this heatmap.