Microdosimetry and Dose-Averaged LET Calculations of Protons in Liquid Water: A Novel Geant4-DNA Application

1 Introduction

The radiobiological effects of ionizing radiation can be traced back to the microscopic patterns of energy deposition in irradiated media, which can be described quantitatively with the formalism of microdosimetry [1–3]. In view of the growing awareness of the variations of proton Relative Biological Effectiveness (RBE) with increasing Linear Energy Transfer (LET) [4–6], several studies have turned to the framework of microdosimetry to approach problems such as the characterization of LET distributions for proton therapy beams and the determination of their RBE [7–10]. Indeed, by obtaining the LET distributions it is possible to optimize treatment plans, taking into account the actual biological impact of radiation by working on LET objectives and constraints. Microdosimetry, on one hand, permits the calculation of RBE from a microscopic approach by means of the Microdosimetric-Kinetic Model (MKM) [11–13] for cell killing endpoint and, on the other hand, provides physical concepts and computational tools to calculate macroscopic dose-averaged LET (${\bar{L}}_{\text{D}}$) distributions. These, in fact, being representative both of the amount of energy imparted and of its spatial concentration, are strictly related to the biological damage induced by radiation.

The microscopic analogue of the LET is called lineal energy (y), which is a stochastic quantity that accounts for the energy imparted per unit length to a certain sensitive volume, called site, by an incident particle and all its secondary products [14]. As for the macroscopic distribution of LET, the microscopic distribution of y can be biologically characterised by its dose-weighted average, ${\bar{y}}_{\text{D}}$ [15]. The difference between ${\bar{L}}_{\text{D}}$ and ${\bar{y}}_{\text{D}}$ is given by energy straggling, which can be accounted for by considering the weighted average δ₂ of the energy imparted to the site per individual collision ɛ_c of the incident particle. The mathematical description of this relation was given by Kellerer, and can be used to extract the macroscopic dose-averaged LET of a beam from its microdosimetric distributions [16–18].

The calculation of microdosimetric distributions and their moments can be carried out both analytically or by means of Monte Carlo (MC) simulations. In this work, we present the main features of a MC code developed with Geant4 [19–21] for the computation of microdosimetric distributions generated by proton track segments in liquid water via track-structure simulations done with the package Geant4-DNA [22–25]. Our code provides two sampling methods, called uniform and weighted, for the computation of the quantities of interest in spherical sites. We also highlight the description of the novel scoring method followed to compute of the energy imparted to the site per electronic collision of the primary proton, in order to connect the obtained ${\bar{y}}_{\text{D}}$ values to the macroscopic ${\bar{L}}_{\text{D}}$ ones avoiding resorting to dedicated calculations as discussed in our previous work [26]. Following the description of the geometry and the sampling algorithm implemented, the results obtained for microdosimetric distributions for different site sizes and energies are presented. Finally, the code is tested and benchmarked by comparing the ${\bar{L}}_{\text{D}}$ calculated with the formula proposed by Kellerer and the macroscopic one. As we found significant deviations from its prediction, probably due to the influence of secondary electrons, here we propose the usage of a quantity, named effective mean chord length, ${\bar{l}}_{\text{F}}^{*}$, with which we calculated effective frequency- and dose-averaged lineal energies, ${\bar{y}}_{\text{F}}^{*}$ and ${\bar{y}}_{\text{D}}^{*}$, respectively, so that we could recover the agreement predicted by Kellerer for protons between 10 and 90 MeV.

2 Materials and Methods

2.1 Theoretical Framework

The use of the LET for the description of the biological effect of different radiation qualities at a small scale presents some limitations, the principal ones being that the concept of LET does not consider the random nature of energy loss along a particle track and that it does not take into account the length of the track relative to a finite target structure. The effects of these limitations, and the ranges of proton energies and site dimensions in which they are relevant to energy deposition, were extensively studied by Kellerer and Chmelevsky [27–30]. For large site dimensions (of the order of 1–10 μm) and small proton energies (below 1 MeV), the energy deposition is strongly influenced by the limited range of protons, since the LET may change significantly along the path through the site. On the other hand, as the energy of the proton increases, energy straggling and eventually δ-rays influx and escape start playing an important role, which becomes more relevant as site dimensions get smaller [14]. The interplay of the various factors affecting the energy deposition in cellular or sub-cellular regions and their contribution to the macroscopic dose-averaged LET is discussed in detail in the following paragraphs.

2.1.1 Basic Microdosimetric Concepts

The basic microdosimetric quantity is the energy imparted to a site by a single energy-deposition event, ϵ_s, where an event is defined as an energy deposition due to particles that are statistically correlated. The energy imparted is a random variable, and its value varies from event to event. Therefore, predictions can be only made on the basis of probability distributions.

The single-event distribution $f\left({ϵ}_{\text{s}}\right)$ of the energy imparted is defined as the distribution of energy deposited in the site by exactly one event. The expectation value of ϵ_s following $f\left({ϵ}_{\text{s}}\right)$ is called frequency-averaged imparted energy per event:

For radiobiological considerations, it is also useful to consider the weighted or dose distribution $d\left({ϵ}_{\text{s}}\right)$ of the energy imparted per event, where the contribution of each event is weighted by the energy it deposits:

The expectation value of ϵ_s following $d\left({ϵ}_{\text{s}}\right)$ is called dose-averaged imparted energy per event ${\overline{ϵ}}_{\text{sD}}$ and it is given by:

From Eq. 3 it is possible to relate both averages through:

where ${\sigma }_{{ϵ}_{\text{s}}}^2$ and $V_{{ϵ}_{\text{s}}}$ are the variance and relative variance of ϵ_s, respectively. In general any continuous stochastic quantity x can be characterized trough its frequency and weighted probability density functions f(x) and d(x), the second being obtained by weighting each event by the value of the stochastic quantity. Therefore, the following general relations can be obtained for the frequency- (${\bar{x}}_{\text{F}}$) and weighted-averaged (${\bar{x}}_{\text{D}}$) values of x:

where σ_x and V_x are the variance and the relative variance of x, respectively.

Another fundamental microdosimetric quantity is the lineal energy, y, defined as:

where $\bar{l}$ is the mean chord length of the site. The chord length l is itself a stochastic quantity and represents the length of a particle track (considered a straight line) within a site. It depends on the shape and dimensions of the site and, if the track is considered as a straight line, the determination of $\bar{l}$ is a purely geometrical problem [18, 31].

The frequency-averaged lineal energy ${\bar{y}}_{\text{F}}={\overline{ϵ}}_{\text{sF}}/{\bar{l}}_{\text{F}}$ is the microdosimetric analogue of the frequency- or track-averaged LET, ${\bar{L}}_{\text{T}}$. If the so called “short track segment condition” holds, i.e., the site dimension is considerably smaller than the proton range so that the LET can be considered constant along the track segment across the site, the frequency-averaged lineal energy and the track-averaged LET coincide: ${\bar{y}}_{\text{F}}={\bar{L}}_{\text{T}}$. It must be noted that, in the specific situations considered in this work, the macroscopic LET values are to be considered as restricted, since the site dimension determines a maximum cut-off to the distance that secondary electrons can travel from the place they are set in motion and therefore, to their kinetic energy. This is especially relevant for clinical applications, where the interest is focused on the energy deposited “locally” by a particle and not on its total energy loss.

2.1.2 The Relative Variance Formula for the Derivation of Macroscopic LET From Microdosimetry

As already mentioned at the beginning of this section, if a medium is exposed to an absorbed dose D, the actual energy imparted to a specified site will depend on various factors [30]. Firstly, the number of independent charged particles traversing the region and their different LET values if a composite beam impinges on the site. Secondly, the chord length associated to each proton crossing the site. Finally, even particles with same LET and chord length can result in a different energy imparted, as both the number of collisions along the chord and the energy lost in each collision might vary. The impact of this last factor, called energy straggling, can be expressed as ϵ_s = L × l ± δ, where L and l are respectively the LET and chord length of a particle, and δ is the stochastic variability of the energy imparted with respect to the expected energy ${\overline{ϵ}}_{\text{s}}=L\cdot l$ due to straggling. The influence of the different random factors on the single-event distribution, under the short track segment condition, can be considered in terms of the relative variance of the different quantities [17], i.e.:

where V_LET, V_l and V_δ are the relative variances of the distributions of LET, chord length and energy straggling, respectively. The computation of the straggling distribution is complicated by the fact that the escape and influx of δ-rays from and into the region of interest presents considerable theoretical difficulties, and the problem needs to be studied with Monte Carlo simulations of δ-ray tracks. The relative variance of the straggling distribution, in particular, is determined by the weighted average of the distribution of energy imparted per individual proton collision:

where ${\bar{\epsilon }}_{\text{c}}$ and ${\sigma }_{{\epsilon }_{\text{c}}}$ are the average and the standard deviation of the distribution of energy imparted per collision, respectively. The relative variance of the straggling distribution can be demonstrated to be equal to [16]:

Consequently, from Eq. 7 the following formula can be obtained:

By using the previous relations, and considering that ${\bar{y}}_{\text{D}}={\overline{ϵ}}_{\text{sD}}/{\bar{l}}_{\text{F}}$, the derivation of the subsequent relation for the dose-averaged LET is straightforward:

As emerges from Eq. 11, identified hereafter as Kellerer’s formula, macroscopic and microscopic quantities are related to each other. Therefore, the calculation through MC track structure simulations of the microdosimetric quantities ${\bar{y}}_{\text{D}}$, δ₂ and the chord length distributions can be used to obtain the dose-averaged LET of a given composite beam.

2.2 MC Track Structure Simulations

In order to calculate the various microdosimetric quantities in Eq. 11 through MC track structure simulations, a particle track needs to be simulated and sampled, placing a site in such a way that at least one energy deposit occurs inside it. There are different types of randomness for the interception of a convex body with straight lines; in the following, two types will be distinguished: uniform and weighted [31, 32]. The first type results when a body is exposed to a uniform isotropic fluence of infinite straight lines, while the second occurs when an uniformly distributed random point is chosen in a body and is traversed by straight lines with uniform random direction. From a computational point of view, uniform sampling consist in randomly selecting a position for the center of the site within the region of interest, while weighted sampling consist in selecting an energy transfer point and then placing a site randomly around it. An immediate drawback of the first method is that it is highly inefficient, since the chances of selecting areas with no energy imparted to them are really high. On the other hand, the second method ensures the occurrence of at least one energy transfer within the site, and it is more suitable for MC track structure calculations. The term “weighted” in this context, refers to the fact that the probability of selecting a point centered at a certain elementary volume dV is not uniform but weighted by a factor that compensates for the non-uniform density of transfer points around the particle track [14, 33].

The chord length is a stochastic quantity that depends on the geometry of the site and the straight line associated to each proton in the site. Therefore, its distribution and average values are strictly connected to the sampling method chosen, affecting the shape of Kellerer’s formula [31]. However, with the selection of the appropriate weighting factor, uniform and random sampling methods can be considered as equivalent and, in the case of spherical sites of diameter d, Eq. 11 takes the form [18]:

Straggling, on its side, is an independent random factor from LET and chord length, as even particles with same LET traversing the site with the same chord length might result in a different energy imparted per event ϵ_s. Therefore, it is not affected by the sampling method considered. However, when simulating a beam composed by particles of different LET, attention must be put to not include the variance of the LET distribution into the straggling one for the computation of δ₂, as explained extensively by Bertolet et al. [26].

In the following paragraphs, the main features of the MC application developed for the computation of the various quantities involved in Kellerer’s formula are presented. All simulations were carried out with the Geant4 toolkit (version 10.5), by using our computing cluster hosted at Centro Informático Científico de Andalucía (CICA, Seville, Spain), consisting of 24 computational nodes, each with 2 × 12C AMD Abu Dhabi 6344 (2.6 GHz/6 MB) L3. For every case considered, a number of 250,000 proton tracks was simulated, divided into 250 parallel jobs. The Geant4-DNA package, available within Geant4, was used for this purpose. Compared to Geant4, which provides various physics list options based on a condensed-history approach for the transport of charged particles, Geant4-DNA includes models developed for track structure simulations, which allows the reproduction of each single interaction of neutral and charged particles with liquid water [22–25]. This is especially key to carry out an accurate modeling of electron transport below 1 keV, which affects to energy deposition tallied at sub-micrometric targets (below ∼100 nm) as pointed out by Kyriakou et al. [34]. Geant4-DNA provides various physics list options as well, which mostly differ on the electron transport models used below 255 keV or its algorithm speed [25]. In this work, intended to highlight the code capabilities and illustrate the accuracy with which the dose-averaged LET can be calculated from microdosimetry quantities, we used the “option 2” of the Geant4-DNA constructor, which is an accelerated version of Geant4-DNA default constructor. This package only allows calculations for proton energies below 100 MeV. Therefore, simulations were performed with mono-energetic proton beams in the range 10–90 MeV, choosing the lower limit to ensure the validity of the short track segment condition. Simulations carried out with weighted sampling ranged from 90 min (protons at 10–20 MeV) to 87 h (90 MeV) per 1,000 primary protons; as for simulations using uniform sampling of site position, they took from 170 min (20–40 MeV) to 87 h (90 MeV) including, to improve statistics, a loop of up to 1,000 tries which finished once at least one energy exchange point was found within the site.

For mono-energetic proton beams, the LET distribution approximates a delta function, centered at the LET of the particle for the specific energy considered. In this specific situation, therefore, the relation ${\bar{L}}_{\text{D}}={\bar{L}}_{\text{T}}={\bar{y}}_{\text{F}}$ (being ${\bar{L}}_{\text{D}}$ and ${\bar{L}}_{\text{T}}$ averages on restricted LET values, as stated previously) holds and can be used for the validation of the code [26]. Furthermore, comparisons were also made with unrestricted macroscopic dose-averaged LET values, computed with the formula proposed by Cortés-Giraldo and Carabe [35]:

where n and s are the indexes for an event and a step, respectively, L_sn is the mean energy loss per unit path length in the material, ɛ_sn is the energy deposited by the primary particle along the step and ω_n its statistical weight. Further details of this calculation can be found in [35].

2.2.1 Geometry

Figure 1 shows the geometry of our microdosimetry code. The world volume, i.e., the volume that represents the experimental area and contains all the other components, consist of a water cubic box having half dimensions R_max slightly greater than the maximum range, R_{δ, max}, of δ-rays emitted by the incident protons. This range is calculated with the formula proposed by Tabata et al. [36], and the condition R_max ≥ R_δ, max is necessary to ensure intra-track electronic equilibrium. The scoring volume, i.e., the volume in which the energy transfer points (hits) are detected and stored, is a water slab positioned in the center of the world volume with same transversal dimensions and thickness Z_S = 4d, where d is the diameter of the site. The value of Z_S was set to ensure enough energy transfer points for each track within the scoring volume, being 4d a good compromise between its thickness and calculations carried out comparing the two site positioning methods explained below; more details can be found elsewhere [37]. Simulations were run with spherical sites of 1 and 10 μm diameter. In order to have the same proton kinetic energy distribution at the center of the world volume for both site diameters, the decision to set R_max = R_δ, max + Z_S, max/2 was taken, where Z_S, max is the thickness of the scoring volume corresponding to a site diameter of 10 μm.

FIGURE 1. Scheme of the geometry used for microdosimetry simulations (2D projection). The primary proton (red track) is generated at the surface (on the left) of a cubic water volume (world volume) and travels along the Z axis. The energy transfer points (hits) are then scored in the middle of the volume (shaded area), in a slab whose thickness (Z_S) varies depending on the dimension of the site. The spherical sites are not physically simulated, as they are virtually selected in the scoring volume with an algorithm that depends on the type of randomness considered. Finally, the half-dimension of the world volume (R_max) varies according to the maximum energy of the incident protons, and is slightly greater than the maximum range of the δ-rays emitted.

In our simulations, the beam (red track in Figure 1) originates at the surface of the cubic volume located at z = − R_max and travels along the z-axis. The physical quantities of interest are scored in randomly placed spherical sites, that are virtually identified as regions of the scoring volume thanks to specific algorithms for the selection of their centers. For the site positioning, two sampling methods have been implemented, uniform and weighted random sampling, which are described in the next paragraph.

2.2.2 Scoring Algorithms

The most intuitive way of scoring energy deposition events in micrometric sites along a proton track would be to randomly select a point ${\vec{P}}_{\text{C}}$ in the region of interest with a uniform probability, and to consider all the hit positions ${\vec{P}}_{{\text{h}}_i}$ occurring at a distance $|{\vec{P}}_{{\text{h}}_i}-{\vec{P}}_{\text{C}}|\le d/2$ from it. This sampling method, identified in this work as uniform, is very robust and it is not subject to any bias. The main drawback of this technique is that it is likely to be highly inefficient, as the probability of selecting volumes with no hits is really high and increases for smaller sites and higher energies.

To increase the sampling efficiency, weighted sampling can be used instead. In this case, an energy transfer point ${\vec{P}}_{\text{H}}$ is randomly chosen with uniform probability among the hits of the simulated track. Then, the center of the site ${\vec{P}}_{\text{C}}$ is sampled with uniform probability in a sphere of radius d/2 around ${\vec{P}}_{\text{H}}$, i.e., at a distance $|{\vec{P}}_{\text{C}}-{\vec{P}}_{\text{H}}|\le d/2$. Finally, all the energy transfer points ${\vec{P}}_{{\text{h}}_i}$, within a distance $|{\vec{P}}_{{\text{h}}_i}-{\vec{P}}_{\text{C}}|\le d/2$ are considered for the computation of the microdosimetric quantities in the site. In this way, the presence of at least one hit in the site is always ensured, and a sampling efficiency of 100% can be achieved. The downside of this method is that it is spatially biased towards regions of high density of hits, and a compensation factor must be introduced to “correct” the results and make them equivalent to uniform ones. As suggested by other studies [33, 38], this factor should be equal to the ratio of the number of hits that can be selected to the number of hits located in the randomly placed site.

The working principle of the two methods, and the way they were implemented in our code, is shown in Figure 2. In both cases, a hit selection region within the scoring volume, where ${\vec{P}}_{\text{C}}$ and ${\vec{P}}_{\text{H}}$ are sorted, can be identified. The dimensions of this region are chosen to ensure that all the randomly selected points are always fully included in the scoring volume, and depend on the scoring algorithm considered. In the case of weighted sampling, the compensation factor is set equal to N_sel/N_int (see Figure 3), where N_sel is the number of hits (with energy different from zero) in the hit selection region that can be randomly selected and N_int is the number of hits in the intersection volume of the site and the hit selection region.

FIGURE 2. Graphical representation of the working principle of the two sampling algorithms: (A) uniform sampling and (B) weighted sampling. In both cases, a hit selection region (maroon shaded area) is identified, where the centers of the sites (yellow circles) are sorted. The dimensions of this region depend on the sampling algorithm used and on the diameter of the site, as lateral margins must be left to ensure that the site is always fully included in the scoring volume. For uniform sampling this margin must be set equal to the site radius, while for weighted sampling the margin is equal to the site diameter.

FIGURE 3. Detail of the working principle of the weighted sampling method in Figure 2B. The compensation factor is computed as the ratio between the number of selectable hits (with deposited energy different from zero) in the hit selection region (N_sel) and the number of hits in the intersection volume of the site with the hit selection region (N_int).

Once the position of the site is sorted, three main quantities are computed: the total energy deposited in the site per event, ϵ_s, the total energy deposited for each single proton collision, ɛ_c, and the trajectory length, s, of the proton traversing the site. In particular, for the computation of the energy imparted per individual proton collision, a Sensitive Detector class was specifically coded to identify uniquely the electronic showers generated by each collision. To this end, a distinctive tag is assigned to every primary electron set in motion by a proton collision. Then, taking advantage of the specific structure of the tracking algorithm of Geant4, all secondary electrons created by further collisions of the electron originating the shower are assigned the same tag, which becames an identifier of the shower. In this way, the energy deposited in the site by each δ-ray shower originating in the site or entering the site can be scored independently, as graphically represented in Figure 4.

FIGURE 4. Schematic representation of the scoring method for the evaluation of the energy deposited per single proton collision in the site, ɛ_c. The primary δ-ray generated by the proton collision is identified and assigned one distinctive tag, that is then inherited by all the secondary electrons in the shower. Only the energy deposition events happening within the site are accounted for in the calculation of ɛ_c, including δ-ray influx from outside the site.

3 Results

Figure 5 shows the results obtained for the frequency-averaged energy imparted per event, ${\overline{ϵ}}_{\text{sF}}$, as a function of the energy and the site diameter, for the two sampling methods considered. Taking the uniform sapling method results as a reference, we obtained relative residuals below 5% for all the energies considered. Regarding the weighted-averaged energy deposited per individual proton collision, δ₂, an intermediate step is necessary to extract its value from the simulations. This quantity is an independent random factor from LET and chord length, therefore their variances should not be included in the computation of δ₂ when calculating ${\bar{L}}_{\text{D}}$ with Eq. 12. This means that only the distribution of energy imparted per individual proton collision corresponding exactly to an energy imparted per event ${\overline{ϵ}}_{\text{sF}}$ and a chord length ${\bar{l}}_{\text{F}}$ should be taken into account to calculate δ₂. In our previous work [26], this was done with a dedicated simulation, positioning the site in such a way that the particle traversed it with a path length equal to ${\bar{l}}_{\text{F}}$. In the work presented here, however, a new approach is proposed to extract δ₂ directly from the main MC simulation, by recording the two dimensional plot of the distribution ${\epsilon }_{\text{c}}^{2}f({\epsilon }_{\text{c}})$ as a function of the primary proton trajectory length s within the site. By doing so, projections could be taken for intervals of s with size 0.1d, to build a graph of δ₂ as a function of s. Examples of this graphs, for the weighted random sampling, are shown in Figure 6. For values of s < d/2, the main contribution to δ₂ is given by the δ-ray influx, whose influence decreases rapidly as s increases. For s > d/2, on the other hand, the value of δ₂ can be considered as constant. Therefore, since ${\bar{l}}_{\text{F}}=2d/3>d/2$ for a spherical site, the final values of δ₂ were extracted from the average of the points in the range s > d/2, weighted by their uncertainty. In Figure 7 the values of δ₂ for the different sampling methods and site diameters considered are represented together as a function of the incident proton kinetic energy: again, taking the uniform sampling results as a reference, we obtained relative residuals below 6%.

FIGURE 5. Frequency-averaged energy imparted per event ${\overline{ϵ}}_{\text{sF}}$ as a function of the proton energy. Uniform (open markers) and weighted (full markers) random sampling results are plot together and compared for sites of 1 (blue triangles and diamonds) and 10 μm (red circles and stars) diameter.

FIGURE 6. Variation of δ₂ as a function of the proton trajectory length s within the site, obtained with weighted random sampling for the different energies. (A) sites of 1 μm diameter, (B) sites of 10 μm diameter. For the lower δ₂ values, the main contribution to δ₂ is given by the δ-ray influx, whose influence decreases quickly as s increases. For s > d/2, the value of δ₂ can be considered as constant. Since ${\bar{l}}_{\text{F}}>d/2$, the final values of δ₂ were extracted by averaging the points (weighted by their uncertainties) in the range s > d/2 (solid lines).

FIGURE 7. Weighted average of the distribution of energy imparted per individual proton collision δ₂ as a function of the proton kinetic energy. Uniform (open markers) and weighted (full markers) random sampling results are compared for sites of 1 (blue triangles and diamonds) and 10 μm (red circles and stars) diameter sites.

The results obtained for the frequency-averaged lineal energy ${\bar{y}}_{\text{F}}$ with both sampling methods are presented in Figure 8. In this figure, ${\bar{y}}_{\text{F}}$ values and ${\bar{L}}_{\text{D}}$ values obtained with Eq. 12 are compared, to verify the relation ${\bar{y}}_{\text{F}}={\bar{L}}_{\text{D}}$, valid for mono-energetic beams. Constant LET condition was fulfilled for all scenarios, with the worst one being 10 MeV protons on 10-μm in-diameter sites, in which a change below 2% was observed. Taking the ${\bar{y}}_{\text{F}}$ values as a reference, the relative residuals of ${\bar{L}}_{\text{D}}$ for the different site diameters and sampling methods considered are represented in Figure 9. While for the 10 μm diameter site we find residuals approximately constant with energy and below 6%, for the 1 μm diameter site deviations increase with energy, reaching values of up to 30%. This effect is more evident looking at the weighted random sampling points, which are not subject to the lack of statistics that affects the uniform random sampling method, especially at higher energies. Moreover, this increasing difference between ${\bar{L}}_{\text{D}}$ and ${\bar{y}}_{\text{F}}$ for the 1 μm site, results in an anomalous behaviour of ${\bar{L}}_{\text{D}}$ if compared to the macroscopic dose-averaged LET obtained from Eq. 13, represented as a black shaded area in Figure 8. Indeed, this last value, being unrestricted, represents the maximum ${\bar{L}}_{\text{D}}$ achievable with Kellerer’s formula and should never be crossed, a statement that is even stronger in the case of the 1 μm site against what happens in the figure.

FIGURE 8. Frequency-averaged lineal energy ${\bar{y}}_{\text{F}}$ as a function of the incident proton kinetic energy. Values obtained for uniform (open markers) and weighted (full markers) random sampling methods are compared for sites of 1 (blue triangles and diamonds) and 10 μm (red circles and stars) diameter. For mono-energetic proton beams the relation ${\bar{L}}_{\text{D}}={\bar{y}}_{\text{F}}$ should hold, therefore ${\bar{y}}_{\text{F}}$ values are here superimposed to the respective ${\bar{L}}_{\text{D}}$ values derived from Eq. 12, represented as shaded areas whose thickness is given by the error associated to ${\bar{L}}_{\text{D}}$ for each diameter site. Finally, the macroscopic unrestricted ${\bar{L}}_{\text{D}}$ computed with Eq. 13 is also reported as a reference.

FIGURE 9. Relative residuals of ${\bar{L}}_{\text{D}}$ from ${\bar{y}}_{\text{F}}$ for the different site diameters and sampling methods considered.

To solve this problem, we proposed the use of new quantities named effective mean chord length, ${\bar{l}}_{\text{F}}^{*}$, and effective frequency- and dose-averaged lineal energy, ${\bar{y}}_{\text{F,D}}^{*}$, respectively. As done for the energy imparted per individual proton collision, a two dimensional plot of ϵ_sf(ϵ_s) as a function of s was produced from the simulation, and projections were extracted for intervals of s equal to 0.1d. In this way, plots of the mean energy imparted per event ϵ_s as a function of the proton trajectory length could be built, as the ones shown in Figure 10 for weighted random sampling. When proton tracks do not cross the site or pass tangent to it, i.e. s ∼ 0, the only contribution to ϵ_s is given by δ-ray influx. Then, for higher s values, ϵ_s grows linearly with the trajectory length. Intersecting the frequency-averaged energies imparted per event obtained considering the full data set and the linear fits of ϵ_s as a function of s shown in Figure 10, the value of s corresponding exactly to an energy imparted per event equal to ${\overline{ϵ}}_{\text{sF}}$ should be obtained. In principle, this value should correspond to the “geometrical” mean chord length ${\bar{l}}_{\text{F}}=2d/3$, however, we found out that this was not the case as the mean chord lengths obtained, which we named effective mean chord lengths ${\bar{l}}_{\text{F}}^{*}$, were slightly different and varied with proton energy (see Figure 11).

FIGURE 10. Plot of the variation of the mean energy imparted per event, ϵ_s, as a function of the proton trajectory length s, described in the site, obtained with weighted random sampling for the different energies and sites: (A) 1 μm diameter; (B) 10 μm diameter. For s ∼ 0, when proton tracks do not cross the site or pass tangent to it, the only contribution to ϵ_s is given by δ-ray influx. Then, for higher s, ϵ_s increases linearly with the trajectory length. The frequency-averaged energies imparted per event ${\overline{ϵ}}_{\text{sF}}$ obtained considering the full data set are represented in this plot as horizontal dashed lines, each colour corresponding to one proton energy. The intersection point of the linear fit of ϵ_s as a function of s with ${\overline{ϵ}}_{\text{sF}}$ corresponds to the effective mean chord length ${\bar{l}}_{\text{F}}^{*}$.

FIGURE 11. Effective mean chord length ${\bar{l}}_{\text{F}}^{*}$ as a function of the proton energy obtained with uniform (open markers) and weighted (full markers) random sampling for sites of 1 μm (blue triangles and diamonds) and 10 μm (red circles and stars) diameter. The “geometrical” mean chord lengths ${\bar{l}}_{\text{F}}=2d/3$ are reported for comparison as yellow horizontal dashed lines.

Using the effective chord lengths for the computation of frequency and dose-averaged effective lineal energies, the results shown in Figure 12 could be obtained, where all the values of ${\bar{L}}_{\text{D}}$ obtained with Eq. 12 lay below the macroscopic unrestricted dose-averaged LET, as expected. In this figure, a general better agreement between uniform and weighted random sampling methods emerges, if compared to Figure 8. Relative residuals of our derived ${\bar{L}}_{\text{D}}$ values from ${\bar{y}}_{\text{F}}^{*}$ were found to be below 6% for both site dimensions and for all energies (see Figure 13) with the only exception of the points obtained for the two higher energies and the smallest site with uniform random sampling, which are strongly affected by the lack of statistics.

FIGURE 12. Effective frequency-averaged lineal energy ${\bar{y}}_{\text{F}}^{*}$ as a function of the incident proton kinetic energy. Values obtained for uniform (open markers) and weighted (full markers) random sampling methods are compared for sites of 1 (blue triangles and diamonds) and 10 μm (red circles and stars) diameter. Values of ${\bar{y}}_{\text{F}}^{*}$ are superimposed to the respective ${\bar{L}}_{\text{D}}$ ones derived considering ${\bar{y}}_{\text{D}}^{*}$ instead of ${\bar{y}}_{\text{D}}$ in Eq. 12, represented as shaded areas whose thickness is given by the error associated to ${\bar{L}}_{\text{D}}$ for each diameter site. Finally, the macroscopic unrestricted ${\bar{L}}_{\text{D}}$ computed with Eq. 13 is also reported as a reference.

FIGURE 13. Relative residuals of ${\bar{L}}_{\text{D}}$ from ${\bar{y}}_{\text{F}}^{*}$ for the different site diameters and sampling methods considered.

4 Discussion

In this work, a novel MC code for the computation of microdosimetric distributions generated by proton track segments in liquid water has been presented. Two different algorithms were implemented for the scoring of the quantities of interest: uniform and weighted random sampling for site positioning. While the first is ideally better as it is not subject to any bias, it is also very inefficient due to the high chances of selecting sites with no energy depositions within. On the contrary, weighted random sampling has a 100% efficiency, but might be subject to biases related to the choice of the weighting factor. One of the main elements of novelty of our code is the inclusion of an algorithm for the computation of the distribution of energy imparted per individual proton collision. This distribution, which is impossible to measure experimentally, reflects the action of energy straggling and δ-ray influx and escape on the energy deposited in the site, and it is fundamental for the computation of dose-averaged LET from microdosimetric quantities with Kellerer’s formula (Eq. 11).

In this study, our weighted random sampling algorithm was optimized so that the results obtained could be considered equivalent to uniform ones, with relative differences generally lower than 6%. Then, the code was used to run various simulations with mono-energetic proton beams, with energies in the range 10–90 MeV impinging on spherical sites of 1 and 10 μm diameter. Since for mono-energetic proton beams the relation ${\bar{L}}_{\text{D}}={\bar{y}}_{\text{F}}$ holds, frequency-averaged lineal energy values obtained with both sampling methods were used as a reference for the validation of the code and the comparison with dose-averaged LET values obtained with Eq. 12. Furthermore, since the site dimensions considered were generally smaller than the maximum range of δ-rays emitted by the incident proton, the ${\bar{L}}_{\text{D}}$ values obtained in this way were expected to be restricted and always lower than the macroscopic unrestricted dose-averaged LET values, which were computed with Eq. 13.

While generally no evident problems emerged in the results obtained for the 10 μm diameter site, the 1 μm site arose questions about the validity of Kellerer’s formula, as the ${\bar{L}}_{\text{D}}$ curve intersected the unrestricted one. This apparent contradiction was solved with our proposal of introducing the concept of effective mean chord length ${\bar{l}}_{\text{F}}^{*}$, defined as the trajectory length corresponding exactly to an energy imparted per event to the site equal to the ${\overline{ϵ}}_{\text{sF}}$ obtained from the MC. Following this definition, ${\bar{l}}_{\text{F}}^{*}$ should be naturally equal to the “geometrical” mean chord length ${\bar{l}}_{\text{F}}=2d/3$. However, this was proved to be wrong, as the effective value of the mean chord length not only differed slightly from the expected one, but varied also with the energy of the incident proton. This behaviour is probably related to the action of δ-ray influx, which affects the distribution of energy imparted to the site per event, increasing the actual value of ${\overline{ϵ}}_{\text{sF}}$ with respect to what would be obtained by only considering direct traversals of proton tracks through the site, especially for the highest energies. By using the ${\bar{l}}_{\text{F}}^{*}$ values for the computation of frequency- and dose-averaged effective lineal energies, ${\bar{y}}_{\text{F}}^{*}$ and ${\bar{l}}_{\text{D}}^{*}$, respectively, a better agreement was observed both between uniform and weighted random sampling results and between ${\bar{L}}_{\text{D}}$ and ${\bar{y}}_{\text{F}}^{*}$ curves. Furthermore, the nature of restricted LET emerged clearly in this case, being all the ${\bar{L}}_{\text{D}}$ values computed with Kellerer’s formula below the macroscopic ones, and generally lower for the 1 μm diameter site. Indeed, the difference between restricted and unrestricted dose-averaged LET values decreases as the energy increases, especially expected for the 10 μm site, since the maximum range of δ-rays starts to be comparable to the site dimensions for proton energies below 20 MeV. Although the agreement found with our alternative approach looks promising, we must point out that this work only involved protons from 10 to 90 MeV (currently, Geant4-DNA simulates protons only up to 100 MeV). Thus, a dedicated analysis (beyond the scope of this work) extending the energy range and/or involving different ions and site shapes seems to be of interest to confirm our proposal, as well as cross-checking with other track-structure codes and experimental measurements. For the latter, it must be remarked that the beam energy spread and detector response introduce further variances in the distributions of microdosimetry quantities which must be carefully taken into account in such an analysis.

Concerning the δ₂-values, we observed deviations from previous works when comparing calculations with sites of 10 μm diameter. While the values obtained in this work were between 2.1 and 2.4 keV for protons from 10 to 90 MeV, lower values were reported previously using other approaches. With dedicated simulations, aiming at isolating its calculation from contributions due to the variance of single-event energy imparted variance as discussed in [26], we reported values of about 1.8 keV; previously, with a more approximated calculation with G4EmPenelopePhysics physics list [35], values of about 1.5 keV were reported for the same energy range. These deviations suggest that the methodology followed to calculate δ₂ must be explained in detail, as it may impact the values of ${\bar{L}}_{\text{D}}$ calculated from microdosimetry quantities by means of Eq. 12. Moreover, different versions of the Geant4 toolkit were used across these works, thus the choice of transport models varied between them as well. This issue shows that a dedicated study to quantify the impact of different transport models on the calculated δ₂-value is of interest, but out of the scope of this article as Geant4 itself offers a wide variety of transport models combinations or physics list.

Using a weighted sampling approach in MC simulations to obtain microdosimetric quantities is undoubtedly advantageous in terms of computation time and efficiency. Indeed, the lower efficiency of uniform random sampling has emerged in our results numerous times, being its effects especially evident for sites of 1 μm diameter and proton kinetic energies above 70 MeV. For the same number of primary protons generated, a loss in statistic going from 60% to more than 90% in the worst case was observed, when comparing uniform sampling with weighted one. Therefore, the observed differences in the computational efficiencies clearly justifies the adoption of a weighted random sampling method for the computation of the microdosimetric quantities of interest, since equivalent results for ${\bar{L}}_{\text{D}}$ could be obtained by using the approach presented in this work. This method brings the advantage of permitting the calculation of RBE by means of models either based on microdosimetry, such as the MKM [11–13] or on LET [39–41] by means of calculating microdosimetry quantities only which, contrary to LET, can be measured. Thus, with the microdosimetry data provided by our code, one can derive analytical parameterisations to calculate the RBE by each approach, without having to derive a specific parameterisation for LET distributions [42–44]. Further, we expect the impact on the uncertainties found on our LET values retrieved from microdosimetry will be relatively small, and will not add up significantly to other more important sources of uncertainty, such as the (α/β) of the irradiated tissue.

The current version of the code, available upon request, is intended to be incorporated within the official examples of the Geant4 toolkit dealing with the extension Geant4-DNA, as it shows the main aspects of the calculation of microdosimetric quantities with track segments of protons and ions (such as He, Li, C, O or Si) using spherical sites, and illustrates how ${\bar{L}}_{\text{D}}$ can be retrieved from ${\bar{y}}_{\text{D}}$ calculated with protons. From a modelling perspective, spherical sites are good enough since sphere-based microdosimetry quantities correlate well with the yield of DSBs, as shown recently [45]. Thus, other shapes such as cylinders were not included in our study to avoid further complications in our analysis of the relations between macro- and microdosimetry distributions. However, they could be useful to benchmark our predictions against experimental results from microdosimetry with other shapes, as done previously in [46].

5 Conclusion

The use of a weighted random sampling approach in MC simulations for the calculation of microdosimetric quantities is widely supported in literature for its undeniable computational efficiency if compared to uniform random sampling. However, especially when considering segments of proton tracks traversing a sensitive volume, particular attention must be put in the computation of the weighting factors, to obtain unbiased results equivalent to uniform ones. In this work, we addressed this problem by developing a MC application for the computation of microdosimetric distributions generated by segments of proton tracks. A remarkable feature of our code, which would constitute a novelty as an official example of Geant4 toolkit, is the possibility of computing the distribution of energy imparted per single proton collision. This distribution, indeed, is strongly connected to the effect of energy straggling and δ-ray influx and escape on the energy deposited in a site per event, and represents the liaison between the microscopic dose-averaged lineal energy and the macroscopic dose-averaged LET, as theorised by Kellerer. To test and benchmark our code, we ran various simulations with mono-energetic proton beams, comparing frequency-averaged lineal energy results with dose-averaged LET ones obtained with Kellerer’s formula. As predicted theoretically for proton beams of this kind in the short track segment condition, we obtained a very good agreement between the two quantities, with percentage differences always lower than 6%. In order to obtain this result, we resorted to a new quantity proposed in this work, the effective mean chord length ${\bar{l}}_{\text{F}}^{*}$, defined from the value of the primary proton trajectory length corresponding to the calculated frequency-averaged energy imparted to the site per single event. In fact, we observed that the action of δ-ray influx in the site strongly affected the distribution of energy imparted to the site per event, displacing its mean value with respect to what would be obtained by only considering direct traversals of proton tracks through the site. For this reason, further calculations at other proton energies and with other ions are of interest to assess the validity of our proposal beyond the limits of this work (protons from 10 to 90 MeV). Such a further study may suggest the introduction of minor corrections of the fitting parameters of our analytical microdosimetric models, described in [42, 43].

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

Conceptualization, AB-R, AB, AC, and MC-G; Methodology, AB-R, AB, AC and MC-G; Simulations, AB-R, AB, GB; Formal Analysis, AB-R, AB and MC-G; Resources, AC and MC-G; Writing-Original Draft Preparation, AB-R. and MC-G; Writing-Review and Editing, AB-R, AB, GB, AC, and MC-G; Project Administration, AC and MC-G; Funding Acquisition, MC-G.

Funding

This work was funded by Grant RTI2018-098117-B-C21 funded by MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”. The work also received funds from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie Grant No. 675265, OMA--Optimization of Medical Accelerators, and from the Council for Economy, Industry, Knowledge and Universities of the Andalusian Government under Grant No. P18-RT-1900.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

All the Monte Carlo simulations presented in this work were carried out at the FIS-ATOM computing cluster hosted at CICA (Seville, Spain).

References