Investigating the feasibility of TOPAS-nBio for Monte Carlo track structure simulations by adapting GEANT4-DNA examples application

At the beginning of this section, we will briefly introduce the examples that were re-simulated during this study and previously carried out with GEANT4-DNA by Incerti2018. All simulations were performed in liquid water. Details such as simulation parameters will be further discussed in the following subsections.

2.1.1. Dnaphysics example

2.1.2. Range example

The energy-dependent range r of one particle is determined as the sum of all step lengths between two interaction points. The step length is calculated using the three-dimensional Pythagorean theorem and the coordinates x, y, z of each trajectory point s:

$$\begin{eqnarray}&&r=\sum _{s=1}^{{N}_{s}-1}\sqrt{{({x}_{s}-{x}_{s+1})}^{2}+{({y}_{s}-{y}_{s+1})}^{2}+{({z}_{s}-{z}_{s+1})}^{2}},\end{eqnarray} \tag{ 1 }$$

where N_s is the maximum step number.

2.1.3. Spower example

The energy-dependent stopping power s of the particles is studied. For this purpose, the sum of the deposited energy E_s of each step s per primary particle is divided by the corresponding range r (see above):

$$\begin{eqnarray}&&s=\displaystyle \frac{{\sum }_{s=1}^{{N}_{s}}{E}_{s}}{r}.\end{eqnarray} \tag{ 2 }$$

2.1.4. Mean free path (Mfp) example

The mfp is the average length of the first step of a particle with a certain kinetic energy. The inelastic mean free path (imfp) is the average length of the first step of a particle considering only inelastic interactions, in which case elastic interactions are not simulated. Since the coordinates of the first trajectory point are located in the origin, this results in the following expression:

$$\begin{eqnarray}&&\mathrm{mfp}=\displaystyle \frac{{\sum }_{e=0}^{N}\sqrt{{x}_{1,e}^{2}+{y}_{1,e}^{2}+{z}_{1,e}^{2}}}{N}.\end{eqnarray} \tag{ 3 }$$

Here N is the number of events e. This equation is also valid for the calculation of the imfp.

2.1.5. W value example

The W-value w is determined by the ratio of the initial particle energy E and the mean number of ionizations $\overline{{N}_{{ion}}}$ per event:

$$\begin{eqnarray}&&w=\displaystyle \frac{E}{\overline{{N}_{{ion}}}}.\end{eqnarray} \tag{ 4 }$$

In the simulations performed with TOPAS (Perl et al 2012), we used TOPAS version 3.2.p2, which is based on the GEANT4version 10.05.p01. The features of TOPAS-nBio are made accessible by several extensions provided by the TOPAS collaboration. We used TOPAS-nBio v1.0-beta.1. Even though the latest GEANT4 version (GEANT4 10.07) is already released, there is no disadvantage of running the simulations based on GEANT4 version 10.05, as there are no significant differences when comparing the example results from those two different GEANT4 versions as shown in the appendix.

Most of the necessary simulation parameters, such as those concerning geometry and radiation source, are described in Incerti2018 and could thus be easily adapted in the TOPAS simulations. All remaining parameters were extracted from the simulation files provided on GitHub ⁵ by the GEANT4 collaboration. These are summarized in table 1. In each simulation, except of the dnaphysics example, an isotropic particle source was placed at the center of a sphere with a radius of 1 m. This was filled with liquid water (G4_WATER) since cross sections data in GEANT4-DNA are only available for liquid water and DNA-related materials. Using G4_WATER, water has a density of 1.0 g cm⁻³ and a mean ionization potential of I_W  = 78 eV. In the dnaphysics example application the particle source was placed in the middle of a box with an edge length of 100 μm. The physics lists G4EmDNAPhysics_option2, G4EmDNAPhysics_option4 and G4EmDNAPhysics_option6 were used in separate runs for all simulations concerning electrons. Simulations for protons and alpha particles in the range example were performed using G4EmDNAPhysics_option2. For a correct calculation of the stopping power of electrons, G4EmDNAPhysics_stationary_opt2, G4EmDNAPhysics_stationary_opt4 and G4EmDNAPhysics_stationary_opt6 needed to be used, since here the kinetic energy of the particle is constant in each step as described by Incerti2018. In most simulations the default tracking cut of the particles was applied; details are listed in table 1. The number of primary particles was set between 10² and 10⁴ depending on the particle type and example application similar to Incerti2018. In the mfp example, the mfp was once investigated including all electron processes and a second time considering only inelastic interactions (imfp) as done by Incerti2018. For the calculation of the imfp, we used a physics list adopted from TsEmDNAPhysics, so that elastic interactions could be inactivated. All other physics models and processes were implemented as it is the case in G4EmDNAPhysics_option2, G4EmDNAPhysics_option4, and G4EmDNAPhysics_option6 which is specified in table 2.

Contrary to GEANT4, the simulations are not pre-programmed in TOPAS, and the results are not directly obtained after running the simulation as it is the case in GEANT4. In order to still be able to calculate these corresponding quantities, we used the tuple scorer in our simulations. This scorer stores information related to a particle track such as the particle type and the x-, y- and z-position of the trajectory step point in one file. In addition, the user can decide whether additional tracking information should be included in the output file by selecting items from the following list:

• EventID (number of the primary particle)
• TrackID (identification number to differentiate between tracks of primary and secondary particles)
• Step number (number of the step of the corresponding track)
• Particle name
• Process name
• Volume name
• Volume copy number (copy number of a volume, if its geometry is copied in the simulation)
• ParentID (track identification numbers of the parent tracks of chemical reactions)
• Vertex position x/y/z (coordinates of the step points of secondary particle tracks in relation to their origin)
• Global time (time of flight)
• Energy deposited
• Kinetic energy (before interaction)

The respective quantity of each example could then be calculated from the corresponding tracking information. For some simulations, it was also useful to add filters to the scorers, e.g. when focusing exclusively on primary particles or particles of a given type. This can be a helpful feature, for example, to keep the output files small and to simplify and speed up the analysis. An overview of possible scoring filters is provided on the TOPAS documentation website ⁶ by the TOPAS collaboration. The tuple scorer provided by TOPAS does only account for interactions with an energy transfer larger than 0 eV. In order to score elastic scatterings also, the tuple scorer was adapted to score events with an energy transfer of 0 eV. All data were stored in a ROOT file.

To validate our method, simulation results using TOPAS-nBio based on GEANT4 version 10.05 were benchmarked with those of GEANT4-DNA using the same software version. The simulations with GEANT4-DNA were performed using the pre-programmed example files included in GEANT4. When running the simulation files, the user directly receives a file in which the results are listed. Thus, there is no need for additional calculations.

The occurrence of physical processes stated in table 3 for electrons of energies between 10 and 60 eV was studied using the physics list G4EmDNAPhysics_opt2. An isotropic source emitting 100 electrons was placed at the center of a cube with an edge length of 100 μm filled with liquid water. The tracking of the particles was terminated by applying the default tracking cut of 7.4 eV. The tuple scorer was applied to investigate detailed track structure information. Post-processing those data, the following investigations were made:

• (i)  
  The frequency of electron processes was investigated by counting each process in total of all 100 primary particles and their secondaries.
• (ii)  
  Since the deposited energy is of particular interest when considering the effect of radiation on DNA scales, this parameter was studied in more detail. The mean deposited energy of a single interaction was calculated per process.
• (iii)  
  Additionally, the sum of the deposited energy of each process per primary particle was calculated.
• (iv)  
  The maximum energy deposited for each process was investigated.
• (v)  
  The chronological order in which the processes of different events and for different initial energies occur along one track was investigated.

In order to evaluate the output of the simulations from TOPAS, we used ROOT version 6.20.04, an open-source analysis framework developed and provided by CERN (Brun and Rademakers 1997), and Python 2.7.17. ROOT offers the user the possibility to choose Python for an evaluation of the ROOT files by importing the module root_numpy into Python. More precisely, we used the function root2array for converting a tree, a special structure in a ROOT file containing the stored data, into a numpy structured array. Then, all tracking data could be commonly accessed and processed with Python.

In this study, the corresponding properties were determined in each event individually and then averaged over the complete set of events to determine a mean value and the associated standard deviation. Since no statistical values can be estimated by simulating just one run of the dnaphysics example, the standard deviation was calculated using the results of 20 independent runs performed with different random seeds.

3.1.1. Dnaphysics

Figure 1 shows the frequency of individual processes simulated with TOPAS-nBio and GEANT4-DNA using different physics lists. In the bottom panels, the relative deviation from TOPAS-nBio to GEANT4-DNA is visualized. For each process the deviations are smaller than 1.2%. Comparing the counts of each process using G4EmDNAPhysics_opt2, the greatest deviation is 1.2% for electron attachment and the average deviation of all processes using is physics list option 0.2%. Using G4EmDNAPhysics_opt4 and G4EmDNAPhysics_opt6, a maximum deviation of 1.2% and 1.1% respectively is observed for hydrogen excitation. On average the values of TOPAS-nBio and GEANT4-DNA deviate by 0.1% using these two physics lists. However, all values between TOPAS-nBio and GEANT4-DNA coincide within one standard deviation.

Zoom In Zoom Out Reset image size

In addition, figure 2 shows the comparison of the counts per process simulated with TOPAS-nBio using these three physics list options by calculating the percentage deviation between two options. The physics list that is listened first in the legend was always selected as the 'original data'. For example, the calculation of the data series 2 versus 4 was performed as follows:

$$\begin{eqnarray}&&{d}_{p}=\displaystyle \frac{{x}_{p,{opt}2}-{x}_{p,{opt}4}}{{x}_{p,{opt}2}}\ast 100.\end{eqnarray} \tag{ 5 }$$

Here, d_p is the percentage difference of process p, ${x}_{p,{opt}2}$ the counts of process p using physics lists option 2 and ${x}_{p,{opt}4}$ the counts of process p using physics lists option 4.

Zoom In Zoom Out Reset image size

The differences between the individual physics lists are quite low for proton and hydrogen processes since the models of these processes are the same for each physics list option. The deviations are of statistical nature. For electron processes, the absolute differences are between 3.6% and 255%. On average, electron processes deviate by 46% for the different physic list options, thus, the correct choice for simulations in TOPAS-nBio is of great importance since simulation results can differ significantly depending on the applied physics list.

3.1.2. Range

In figure 3(a), the range of the investigated particles is shown as a function of initial energy calculated with TOPAS-nBio and GEANT4-DNA using different physics lists. In TOPAS-nBio, more initial energies were investigated compared to the GEANT4-DNA simulations. Hence, the individual data points produced with TOPAS-nBio are connected by a line to guide the eye. Additionally, the ICRU90 (ICRU 2016) ranges of the particles are illustrated. In figure 3(b), the deviations between the TOPAS-nBio simulations and GEANT4-DNA are portrayed. Considering electrons, large relative deviation can only be observed at small, absolute ranges. For initial electron energies above 500 eV, deviations are smaller than 2%. Regarding option 4, the maximum deviation is slightly above 20%; considering option 2, the deviations range from −7.7% to 14.2%, and using option 6, a maximum deviation of 11.9% at an energy of 40 eV is observed. Nevertheless, the mean deviation for electrons using all three options is statistically disturbed around zero and all deviations are in order of one standard deviation. Here, no systematic deviations can be observed. In figure 3(c), the ranges for electrons of energies from 10 to 100 eV are shown in a zoom-in view. As described by Francis et al (2011), this range profile results from cross sections of inelastic scattering processes implemented in the physics list, which are significantly changing for small energies. This is a consequence of the tresholds for inelastic interactions located in this energy range. By undergoing an inelastic interaction, the particle looses kinetic energy reducing the residual range. Thus, these peaks at lower energies occur as the initial energy exceeds new tresholds for inelastic interactions. Since the tresholds and cross sections of each process differ in the different physics list G4EmDNAPhysics_opt2, G4EmDNAPhysics_opt4 and G4EmDNAPhysics_opt6, the profile of the range at low energies is also different for each simulated physics list.

Zoom In Zoom Out Reset image size

When comparing the proton ranges in figure 3(d), it is interesting to note that the deviations between TOPAS-nBio and GEANT4-DNA are largest at an initial energy of 1 keV (30%) and then decrease systematically down to a minimum value of 0.6%. These deviations result from the fact that in the GEANT4-DNA pre-programmed examples the tracking cut was not adjusted as described in the study by Incerti2018. In the TOPAS-nBio simulations the tracking cut was set higher resulting in a lower range. Comparing the ranges calculated with TOPAS-nBio and GEANT4-DNA with ICRU90 data, it can be seen that the agreement between TOPAS-nBio and ICRU values is more accurately compared to GEANT4-DNA. For a more accurate comparison between GEANT4-DNA and TOPAS-nBio, the simulations were repeated in TOPAS-nBio, this time using the same tracking cut (100 eV) as in GEANT4-DNA. With this adjustment, the values agree well with those of GEANT4-DNA (see figure 3(e)) and a maximum deviation of 0.6% could be obtained.

The ranges of alpha particles (see figure 3(a)) comply well with both the GEANT4-DNA values and the ICRU values, while deviations are smaller than 3% for all considered energies and coincide within one standard deviation.

3.1.3. Spower

Figure 4(a) shows the mean value of the stopping powers of all three physics list options as a function of the initial particle energy for electrons. Due to practical reasons, only four values were simulated in TOPAS-nBio per physics list. Therefore, the x-axis in 4(a) is discrete and the values of TOPAS-nBio and GEANT4-DNA are plotted side by side for a more concise illustration of the results. While using G4EmDNAPhysics_stationary_opt2 and G4EmDNAPhysics_stationary_opt6 a percentage deviation of 0.12% at maximum of the electron stopping powers is obtained, the deviation for G4EmDNAPhysics_stationary_opt4 is ten times larger with 1.3%. It is noticeable that for G4EmDNAPhysics_stationary_opt2 and G4EmDNAPhysics_stationary_opt6, the stopping powers obtained with TOPAS-nBio are always larger than those from GEANT4-DNA, and contrarily for G4EmDNAPhysics_stationary_opt4, all values from TOPAS-nBio are smaller than those from GEANT4-DNA. However, all values between TOPAS-nBio and GEANT4-DNA agree within one standard uncertainty. For transparency, the uncertainties of the spower example are shown in figure 4(a).

Zoom In Zoom Out Reset image size

3.1.4. Mfp

3.1.5. W value

In figure 5(a) the frequency of each process in relation to the initial particle energy is illustrated. The most common process for all energies is elastic scattering followed by vibrational excitation. In general, the number of those processes increases with the initial energy, due to the increasing range of the particles noting that in figure 5(a) the y-axis is plotted logarithmically. Both, elastic scattering and vibrational excitation, show a minimal occurrence at an energy of ∼18 eV. Above an energy of 11 eV, the probability of ionization increases since more ionization energy tresholds are exceeded and the related cross sections also increase (Incerti et al 2010b, Bordage et al 2016). The energy tresholds required for an ionization as well as for an electronic or vibrational excitation are listed in table 4. These processes reduce the range and hence the amount of elastic interactions significantly due to a noticeable energy loss. The number of elastic interactions fluctuate further with increasing initial energy because on the one hand, starting with a higher initial energy more interactions, i.e. also usually a greater range, are possible which results in a higher number of elastic interactions. On the other hand, the kinetic energy of the particle and thus the range is reduced since with increasing initial energy more tresholds for inelastic interactions are exceeded and a higher number of inelastic interactions might occur in one event. The electron solvation process operates as a tracking cut because no chemical interactions are taken into account in any simulations. Since secondary particles were also considered in this simulation, more than 100 electron solvations can take place in each simulation with 100 primary particles as with increasing initial energy more secondary particles are being produced. The tracking of those produced electrons is again terminated by this process.

Zoom In Zoom Out Reset image size

Table 4. Deposited energies in eV by vibrational excitation, electronic excitation and ionization found in our simulations. These energies correspond to the energy tresholds of the associated processes (E ${}_{T}^{vib}$ for vibrational excitation, ${E}_{{T}}^{{e}{x}}$ for excitation, ${E}_{{T}}^{{i}{o}{n}}$ for ionization) implemented in the different models. The vibrational modes are included in the Sanche Excitation Model described in the study of Michaud and Sanche (1987). The energy tresholds for electronic excitation are taken from the Emfietzoglou Dielectric Model (Emfietzoglou and Moscovitch 2002) and ionization levels are described by Dingfelder et al (1998).

Vib. mode	${E}_{{T}}^{{v}{i}{b}}$	Electr. excitation level	${E}_{{T}}^{{e}{x}}$	Molecular orbital	${E}_{{T}}^{{i}{o}{n}}$
2(${\nu }_{\mathrm{1,3}}$)	0.835	A1B1	8.22	1b1	10.79
${\nu }_{T}^{{\prime\prime} }$	0.024	b1a1	10.00	3a1	13.39
${\nu }_{L}^{{\prime} }$	0.061	Ryd A + B	11.24	1b2	16.05
${\nu }_{L}^{{\prime\prime} }$	0.092	Ryd C + D	12.61	2a1	32.29
${\nu }_{2}$	0.204	Diffuse bands	13.77
${\nu }_{\mathrm{1,3}}$	0.417
${\nu }_{3}$	0.46
${\nu }_{\mathrm{1,3}}$+${\nu }_{L}$	0.5
${\nu }_{T}^{{\prime} }$	0.01

In figure 5(b) the average deposited energy in one single inelastic interaction is shown as a function of the initial energy separately for each process. For most of the processes, the mean deposited energy remains almost constant over the whole energy range as a consequence of the corresponding energy tresholds (see table 4). The latter is also observed for ionization, which has the highest mean deposited energy in this energy range. Only for electron solvation an energy-dependent mean deposited energy can be observed. The amount of the energy deposition by this process depends on the kinetic energy, which remains after an inelastic interaction if the treshold of 7.4 eV is reached. Since at low energies, only a few eV are left after an ionization or excitation, the mean deposited energy by electron solvation is very small. In figure 5(c), the mean energy deposited in total per primary particle is shown as a function of the initial energy for each process. Those profiles result from the combination of the frequency of each electron process (figure 5(a)) and the corresponding mean deposited energy (figure 5(b)). Altogether, it can be seen that ionization constitutes the largest amount of deposited energy in this energy range.

Figure 5(d) shows the maximum deposited energies for each process as a function of the initial energy. The steps in the line of ionization and excitation visualize the initial energy needed to reach the different tresholds of these processes. The deposited energy always corresponds exactly to the energy tresholds listed in table 4, since the deposited energy in TOPAS-nBio refers to the locally deposited energy, not the kinetic energy loss of the particles. This becomes especially apparent for ionization because the energy tresholds in this case vary the most. Electron solvation results in a maximum deposited energy of 7.4 eV, since this is the tracking cut applied in G4EmDNAPhysics_opt2.

To further examine the interactions of low-energy electrons, the frequency and sequence of different processes were investigated for electrons with an initial energy of 24.2 eV. In figure 6, the number of steps is illustrated as bars with different colors for each process for electrons with an initial energy of 24.2 eV for 100 single events. On a second y-axis, the corresponding range of the electrons is plotted. An electron energy of 24.2 eV was chosen because, on the one hand, particles with this initial energy have sufficient kinetic energy so that all processes occur, and on the other hand, not too many inelastic interactions take place. In figure 6, it can be seen that the more steps occur per event, the larger the range and the more elastic scattering and vibrational excitation takes place. Furthermore, it is noticeable that when two ionizations or excitations appear in one event, the number of steps and the range are small, as it can be observed exemplarily in the zoom-in view. In a further step, the sequence of the processes and the deposition of the energy of individual events were examined. In figure 7, the order of the occurring processes of three different events, which are marked with a red arrow in figure 6, is shown. In this representation, the processes are illustrated as consecutive bars, filled with a different color for each process, and plotted against their corresponding step number on the x-axis. The width of the bars has no further meaning, it only gets smaller the more steps occur in the track. Additionally, the loss of kinetic energy along the steps is shown in each graph.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Beginning with the track of Event ID 3, the track starts with an ionization, for which reason the particle already loses quite a lot of its kinetic energy. After that, a kinetic energy of only roughly 12 eV is left. This process is followed by a number of elastic interactions alternating with vibrational excitations. The vibrational excitations cause the particle to lose a tiny fraction of its kinetic energy until it reaches the tracking cut of 7.4 eV. The track is then terminated with electron solvation, whereby the remaining kinetic energy is deposited locally. In event 32, the number of steps is significantly smaller and fewer elastic interactions and vibrational excitations are present. However, two additional inelastic processes, ionization and excitation, occur for this purpose through which the particle loses a lot of kinetic energy simultaneously. It is also typical for particles in this energy regime that if there are two inelastic interactions (except of vibrational excitations) in one event, the second one occurs at the end of the track. In this case, the remaining kinetic energy is smaller than the tracking cut, which means that the energy is deposited locally by electron solvation as a last step of the track.

In addition to these two characteristic histories of a particle with an initial energy of 24.2 eV, in Event ID 59, the particle reaches the tracking cut already after having conducted only one step. Such a scenario is only possible if ionization occurs as the first process and such a large amount of energy is transferred to the secondary particle that the residual energy of the primary particle is smaller than the tracking cut. Since in most cases only small fractions of kinetic energy are transferred to secondary particles, this kind of track is rather unlikely.

In conclusion, the simulations regarding electrons of energies between 10 and 60 eV showed that the most frequent process is elastic scattering, while the largest amount of energy is deposited by ionization. Thus, the selection of the correct model for ionization is very important for simulations on molecular levels and for DNA damage studies, and should therefore be well considered. The models vary with the different physics lists and can also be selected independently in TOPAS using TsEmDNAPhysics. The model of the processes could not be validated yet due to the absence of experimental data below 1 keV. However, this investigation should give an insight of how damage on molecular levels and DNA is produced. Especially, when further simulations are made with these models to determine, for example, the distribution and amount of SSB and DSB in relation to the effect of radiation at biological levels.