The energy transfer induced by ionizing radiation in a water medium occurs rapidly (on a scale of femtoseconds (fs)) during the physical stage of water radiolysis and is followed by the formation of radiolytic species. These species are created in a very short time (from femtoseconds (fs) to picoseconds (ps)), mainly through electronic events during the physico-chemical stage [1,2,3,4]. These events—such as thermalization, solvation of sub-excitation electrons, electronic hole migration, and fast electronic recombination—can lead to chemical bond breaks and produce species. After this stage, the chemical stage starts with heterogeneous distributions of radiation-induced spurs along the axis of the tracks from a few picoseconds (ps) to greater than ~100 nanoseconds (ns) in water at room temperature [5]. The spur evolution largely depends on the spatial distribution of these high concentration regions. Methods using Brownian dynamics and Smoluchowski theory have been applied to describe the spur evolution, in which detailed individual trajectories of spherical particles are simulated in chemical reaction-diffusion systems (called "systems" in the following text) [5,6,7,8,9]. In this particle-based representation, only the molecules of interest are explicitly simulated, and the solvent (water) is considered a continuum. Under these conditions, the diffusion of the species occurs until two individual reactants encounter each other (chemical interactions happen when the distance between reactants is less than a pre-defined distance). Over time, the species and their products disperse and are distributed more homogeneously. Their reactions continuously occur until they reach steady-states, where the species concentrations do not change anymore (from 100 microseconds (µs) to several seconds (s)).

Although the detail level is limited when using simplifications (spherical particles, continuum representation), the calculation time remains the main drawback of the particle-based representation when the simulations deal with a large number of species. This makes it difficult to simulate large, absorbed dose radiation systems, where the concentration of produced species may significantly impact the irradiated medium composition. Moreover, the majority of radiation chemistry experiments extend up to the homogeneous regime where radiation-induced species evolve and reach steady-states over long-term scales [5]. Because of its computational cost, the particle-based representation cannot simulate these types of conditions. Instead, compartment-based models can be used because the level of detail is more limited [10].

The compartment-based representation describes species evolution by gathering molecules of the same species as one group placed in a confined volume (or “compartment”). Only molecules in the compartment can react with each other. A basic requirement of compartment-based models is that the distribution of the species is a homogeneous mixture or “well-mixed” in their compartment. Therefore, no information about species positions is described during the evolution. Under this condition, the chemical master equation (CME) [11] is used to describe the time–evolution probability of all the different possible states of the system. In general, the CME is mathematically troublesome. To address this, Gillespie developed a stochastic simulation algorithm (SSA) that equivalently reproduces the CME using a derived Monte Carlo technique [11].

In practice, this “well-mixed” condition is usually not fulfilled since most systems of interest do not begin in homogeneous states. Many of them are first heterogeneously distributed and then tend toward homogeneity as the systems evolve over time, before reaching the fully well-mixed state. For example, in water radiolysis, the transition process between the heterogeneous and homogeneous states may occur continuously from ns to µs or longer, depending on the volume and initial species distribution. Therefore, to be able to apply this well-mixed condition to heterogeneous systems, the compartment-based model proposes dividing the simulation volume into smaller regions (“sub-volume” or “voxels”). In these voxels, the species can be considered as well-mixed (or locally homogeneous). The species can react with each other in the same voxel, and diffusion is modeled by jumps between adjacent voxels. An extension of the CME is the reaction-diffusion master equation (RDME) [12], which is used to account for spatial distributions of species into the states of the system.

2. Simulation of Diffusion

Figure 1. Comparison of spatial (radial) distributions of solvated electrons diffusing in water as computed with the RDME method (blue histogram with unbold error bars) and the SBS method implemented in Geant4-DNA (solid line with bold error bars). The distributions are computed at 10 ns, 20 ns, 50 ns, and 100 ns.

3. Simulation of a Simple Reaction-Diffusion System

For the second benchmark, we investigated a simplified reaction-diffusion system. We considered the same simulation setup with the exception that the solvated electrons did not only diffuse in water but could also interact with each other (Reaction 9 in Table 3). The resolution of the initial mesh in the compartment-based model should be set small enough to ensure that the molecules are well-mixed in each voxel at the start of the simulation. This condition is fulfilled if the voxel size h is smaller than the mean inter-particle distance d. The distance d can be approximated by the expression (V/N)^1/3, with N the number of primary molecules and V the volume over which they are generated. In our simulation case, N = 200 and V=(0.2/21)³ μm³,and we get d = 1.63 nm. We, therefore, selected an initial mesh of 128³ voxels, giving a voxel size of h = 1.5625 nm, which is smaller than the distance d. We then applied the hRDME approach to adapt the size of the mesh during the simulation. It started from the initial mesh and changed over time to coarser meshes of 64³, 32³, 16³, 8³, 4³, 2², and one voxel.

Figure 2 shows the comparison of the time-dependent yield of species (number of molecules) as computed with the particle-based SBS model and the compartment-based model. A very good agreement was observed between the results obtained with both models. Table 2 shows the speed-up factor of the compartment-based model compared to the particle-based SBS model for different time scales covered by the simulations. The speed-up factor is defined as the ratio between the particle-based SBS and the compartment-based model simulation time. While simulations with end times less than 100 ns show a speed-up factor of one order of magnitude, we observed larger factors for longer timescales. The explanation is that, for timescales from 1 ns to 100 ns, the compartment-based model is performed in fine resolution meshes of 128³, 64³, 32³, and 16³ voxels, which makes the compartment-based model less computationally efficient than at longer timescales (from 100 ns to 100 µs) where the system is transferred to coarser meshes. Moreover, at longer timescales, the particle-based SBS is much less efficient than the compartment-based model because it requires using small time-steps when the species reach the volume borders (boundary interactions).

Figure 2. Compartment-based model for the simplified diffusion-reaction system compared with the particle-based SBS simulation.

Table 2. Speed-up factor of the compartment-based model compared to the particle-based SBS model for different simulation end times from 1 ns to 100 µs.

End time	1 ns	10 ns	100 ns	1 µs	10 µs	100 µs
Speedup factor	13	15.7	16.5	2.5 × 102	2.3 × 103	2.3 × 104

4. Full Water Radiolysis Simulation

For the third benchmark, we tested our full SBS-RDME model with a water radiolysis simulation. A water box of 0.2 × 0.2 × 0.2 µm was randomly irradiated on one face by a normal incident beam of 1 MeV electrons. The incident electrons were shot until the sum of all energy deposits in the volume led to an absorbed dose of 100 Gy (about 5 keV for the considered water volume). In this simulation, we used the physics list “G4EmDNAPhysics_option2” for the modeling of the physical stage. The microscopic sub-stage (see Figure 6) extended until 1 ns, at which point, we initiated the compartment-based model with a mesh of 32³ voxels. We then used the hRDME approach to successively transfer species from the finer mesh to coarser meshes of 16³, 8³, 4³, 2³, and finally one voxel where the simulation volume is fully homogeneous.

The list of reactions and reaction rates considered in this work are shown in Table 3. It is worth noting that, in this work, these reactions are considered to be diffusion-controlled reactions where the products are created as soon as the reactants encounter each other.

Table 3. The reactions and reaction rates [15] used in this work that apply at ambient temperature (25 °C).

	Reaction	Reaction Rate (1010 M−1s−1)
1	H• + e−aq + H2O → OH− + H2	2.5
2	H• + •OH → H2O	1.55
3	H• + H• → H2	0.503
4	H2O2 + e−aq → OH− + •OH	1.1
5	H3O+ + e−aq → H• + H2O	2.11
6	H3O+ + OH− → 2H2O	11.3
7	•OH + e−aq → OH−	2.95
8	•OH + •OH → H2O2	0.55
9	e−aq + e−aq + 2H2O → 2OH− + H2	0.636

Figure 3 shows the comparison of time-dependent G-values as computed by the particle-based SBS model of Geant4-DNA and our SBS-RDME model. Note that the comparison is considered from 1 ns when the compartment-based model is activated in the SBS-RDME model. A good agreement was observed for the entire time range between both models, thus, validating our SBS-RDME model implementation. Both SBS-RDME and particle-based SBS models showed that the system reached a steady-state at about 100 µs.

Figure 3. Comparison of time-dependent G-values as computed with the particle-based SBS model and the SBS-RDME model (this work) from 1 ns until 100 µs.

Figure 3 also shows a slight imbalance of neutrality of about 0.18 molecules per 100 eV observed in the G-value of H₃O⁺. The explanation is that during the physical stage, some secondary electrons left the simulated volume as a result of ionizations that took place near the volume borders. These escaping electrons would not produce solvated electrons in the considered volume, while the corresponding ionized water molecules would dissociate in the volume and, therefore, create the slight imbalance of neutrality.

5. Materials and Methods

The main characteristic of our model is the combination of the SBS Brownian dynamics model already available in Geant4-DNA with the compartment-based model using the RDME (so-called “SBS-RDME model”) in water radiolysis simulations. Figure 6 illustrates the simulation scheme of this combination.

Figure 6. The simulation scheme of the combination of the particle-based SBS model with the compartment-based model that uses the RDME (SBS-RDME model).

We began by using Geant4-DNA [19,20,21,22,23,24,25] to model the physical and physico-chemical stage up to 1 ps. Then, the chemical stage was divided into three sub-stages: the microscopic, mesoscopic, and homogeneous sub-stages. The microscopic sub-stage started from the formation of species at the beginning of the chemical stage. As the spatial distribution of the early species (e⁻_aq, H^•, ^•OH, H₃O⁺,…) was concentrated in very small volumes around the tracks, their evolutions largely depended on the species position. Therefore, at this sub-stage, we used the particle-based SBS method of Geant4-DNA to simulate the detailed trajectories of the individual species. Once the simulation reached a time t₁, we stopped the particle-based SBS simulation and initiated a uniform 3D Cartesian mesh for the compartment-based model. This was the beginning of the mesoscopic sub-stage. The initial mesh resolution should be small and selected according to the spatial distribution of species at the end of the microscopic sub-stage. We then used the hRDME approach (see details in Section 4.2.2) to adapt the size of the voxels during the evolution of the system. The system used increasingly coarser meshes over time, to finally merge into one single voxel that covered the full simulation volume at time t_n, when the homogeneous sub-stage started. During the homogeneous sub-stage, we used the CME stochastic process to sample only reactive events. The detail our implementation of the compartment-based model can be found in the original article.

6. Conclusions

In this work, we implemented, in Geant4-DNA, the compartment-based model combined with the SBS Brownian dynamics model already available in Geant4-DNA. We showed that the compartment-based model can reproduce the same species yield obtained by the particle-based SBS of Geant4-DNA but with 100 to 1000 times less computing time. Moreover, our new model can also extend the simulation timescale beyond the microsecond when the system has reached a steady-state. These advantages will allow us to study the production and evolution of reactive oxygen species generated under irradiation with different dose rate conditions, such as in FLASH and conventional radiotherapy.

In this work, we only considered a radiation-induced species system in a cubic water volume. The approach could be applied for more complicated systems, such as complex geometries including biological materials. This may be technically more difficult since we need to handle the two different representations (compartment-based and particle-based) in one simulation. We plan to work on such an extension in the future.