A New Approach to Modeling the Microdosimetry of Proton Therapy Beams
Ramin Abolfath Yusuf Helo David J. Carlson Robert Stewart David Grosshans Radhe Mohan
physics.med-ph physics.bio-ph






^1Department of Radiation Physics and Oncology, University of Texas MD Anderson Cancer Center, Houston, TX, 75031, USA 

^2Department of Radiation Oncology, University of Pennsylvania, Philadelphia, PA 19104, USA 

^3Imanova Centre for Imaging Sciences, Invicro LLC, London, UK 

^4Department of Medical Physics and Bioengineering, University College London, UK 

^5School of Medicine, Department of Radiation Oncology, University of Washington, Seattle, WA 98195, USA






Introduction:
To revisit the formulation of the mean chord length in microdosimetry and replace it by the particle mean free path appropriate for modelings in radio-biology.

Methods:
We perform a collision-by-collision following by event-by-event Geant4 Monte Carlo simulation and calculate double-averaged stepping-length, ⟨⟨ l ⟩⟩, for a range of target sizes from mm down to μm and depth in water.
We consider ⟨⟨ l ⟩⟩ to represent the particle mean free path.

Results:
We show that ⟨⟨ l ⟩⟩ continuously drops as a function of depth and asymptotically saturates to a minimum value in low energies, where it exhibits a universal scaling behavior, independent of particle nominal beam energy. We correlate ⟨⟨ l ⟩⟩ to linear density of DNA damage, complexities of initial lethal lesions and illustrate a relative difference between predictive RBEs in model calculations using mean chord length vs. the proposed mean free path.
We demonstrate consistency between rapid increase in RBE within and beyond the Bragg peak and ⟨⟨ l ⟩⟩, a decreasing function of depth.

Discussion and conclusion:
An interplay between localities in imparted energy at nano-meter scale and subsequent physio-chemical processes, causalities and pathways in DNA damage requires substitution of geometrical chord length of cell nuclei by mean-free path of proton and charged particles to account for a mean distance among sequential collisions in DNA materials.
To this ends, the event averaging over cell volume in the current microdosimetry formalism must be superseded by the collision averaging scored within the volume. The former, is fundamentally global attribute of the cell nuclei surfaces and boundaries and is characterized by their membrane diameters, hence such global indices are not appropriate to quantitatively represent the radio-biological strength of the particles and their RBE variabilities that is associated with the sensitivities to local structure of the collisions and their spatio-temporal collective patterns in DNA materials.




A New Approach to Modeling the Microdosimetry of Proton Therapy Beams

    Ramin Abolfath^1,2†, Yusuf Helo^3,4, David J. Carlson^2, Robert Stewart^5, David Grosshans^1, Radhe Mohan^1,*
    December 30, 2023
=================================================================================================================




§ INTRODUCTION

ICRU Report 36 on microdosimetry [ICRU36] provides general formulations for the mean chord length in terms of the geometrical dimension of the sensitive volumes of interest (SVOI).
Accordingly, ⟨ l ⟩≡l is the mean length of randomly oriented chords in that volume and can be calculated based on the geometrical structure of the target of interest. 
The mean chord length that is an intrinsic geometrical character of SVOI, results from the random interception of the SVOI by a straight geometrical line, superseded for a (single) physical track of a charged particle.
The chord length calculated as such is the mean of shortest Euclidean distances between pairs of end points, crossing closed surface of SVOI.

As originally introduced by Rossi, Kellerer and colleagues (e.g., see for example Refs. [Kellerer1971:RR,Kellerer_Rossi1972:CTRR,Kellerer1984:RR,Kellerer1985:Book,Kellerer1975:REB]), the lineal energy, y = ε / ⟨ l ⟩,
is the energy imparted in a single event, ε, divided by the mean chord length, ⟨ l ⟩.
y is a random variable analogous to particle linear energy transfer (LET). Accordingly, the randomness stems from ε as the chord-length has been averaged out and it is a predetermined number.
For a convex volume, V, such as spherical geometries with area A, a pure geometrical calculation, first performed by Cauchy [Cauchy1908], yields ⟨ l ⟩ = 4V/A.
For a rectangular slab with the transversal thickness d perpendicular to the unidirectional beam central axis, ⟨ l ⟩ = d [Magrin2018:PMB].

The mean chord length, as such, which is based on straight tracks in a micro-meter size in a SVOI, has been vastly used in microdosimetry literature as well as in modelings in the radio-biological effects of any field of radiation.

One of the first studies based on MC technique to calculate chord-length distributions was reported by Birkhoff et al. [Birkhoff1970:HP].
In this type of MC calculations, particles are not scored via a collision-by-collision approach using coarse-grained quantum mechanical scattering cross-sections, as we perform in now-a-days MC models such as Geant4 [Agostinelli2003:NIMA,Allison2006:TNS,Allison2016:NIMPRA] and Geant4-DNA [Incerti2010:IJMSSC] where the tracks, including primary and secondary particles, exhibit the inchoate distribution of energy transfers [Conte2012:NJP,Abolfath2011:JPC,Abolfath2013:PMB,Abolfath2016:MP,Abolfath2017:SR,Abolfath2019:EPJD].
Instead, Birkhoff et al. [Birkhoff1970:HP] considered straight geometrical lines intercepting a sensitive volume of energy-proportional devices such as spherical / cylindrical gas proportional counters.
The rationale behind these simplifications is based on uniform and continuous energy transfer to matter along the tracks of charged particles.
Moreover, the tracks considered in simulations performed by Birkhoff et al., were assumed to be part of a uniform source beam of radiation.

The MC techniques as such,



are therefore analogous to the classical Metropolis MC algorithm proposed for calculation of the numerical value of π, or numerical multi-variable integrations [NumericalRecipes].
By taking a large random sample, the resulting chord-length distribution approaches the correct geometrical distribution in a statistical manner.





This is seemingly a plausible approach, applicable to cavities of the instruments designed for measuring pulse height distributions in a field of radiation such as energy-proportional devices, as details in localities of the collisions, their spatial distribution and compactness within the cavity sensitive volume are irrelevant to overall observable electrical responses. The mathematical model including the calculation of the chord length,
as presented in Birkhoff et al. [Birkhoff1970:HP] and recent modifications and refinements in Refs. [Bolst2017:PMB,Anderson2017:MP,Bolst2018:PMB], attempted to introduce an alternative approach to the Cauchy path length, is consistent with the occurrence of the physical processes in micro-dosimetry devices, considering Rossi's formulation of lineal energy y.

The responses from biological systems, however, seems slightly different from the physical responses we expect to observe from the design of micro-dosimetry proportional counters.
For example, the nano-meter scale linear density of DNA double strand breaks (DSBs) induced by traversing of a track of charged particle in a cell nucleus that determines complexities in lethal pathways of the lesions, cell death and tissue late effects, is expected to be more sensitive to the compactness and linear distribution of the individual collisions, rather than the geometrical dimensions of the cell nuclei [Abolfath2019:EPJD].
The former roots in nano-meter scale localities in physio-chemical phenomenon and effects occurring in DNA material in cell nuclei whereas the latter is global characters associated with the outer structure of the cell nuclei boundaries and their membranes, prone to macroscopic-scale volume effects.
More precisely, such model calculations, as also pointed out, e.g., in Ref. [Carlson2008:RR], only provide an estimate of overall imparted energy within cell nucleus volume without offering an appropriate analysis in resolving nano-scale biologically relevant events, as the geometric chord length of a cell nucleus is in μ m-scale, three orders of magnitude larger than nano-meter scale, the microscopic resolution of DNA-damage.

Moreover, as particles penetrate in tissue, they lose their kinetic energy, hence the compactness of collisions and the complexities in DNA damage are expected to rise. The geometrical chord-length as calculated by Cauchy and implemented in microdosimetry, is a fixed number, e.g., a factor proportional to the dimension of cell nucleus and does not vary as a function of depth.


Let us demonstrate the subtleties and lack of sufficiency in the standard microdosimetry formalism to fully describe and capturing the microscopic / nano-meter responses and critical variabilities in biological complexities, 
by considering passage of a single charged particle in a μm-scale SVOI, a typical representation of a cell nucleus, as schematically shown in Fig. <ref>.
Hence occurrence of various damages in DNA-materials and chromosomes such as double strand breaks (DSBs) are expected to be seen in these volumes [Abolfath2011:JPC,Abolfath2013:PMB,Bianco2015:RRD].

In Fig. <ref>, the black cross-lines (×) represent locations of the site of damages in DNA materials connected diagrammatically to the charged particle by the wiggly lines. The lines are Feynman diagram representation of the photon field propagators in quantum electrodynamics (QED) that describe interaction of charged particles in scattering processes with random interaction sites on DNA and/or the environment of DNA in a cell nucleus.
A process that describes release of OH-free radicals and/or reactive oxygen species (ROS) in indirect DNA damage (the latter), or shell electrons localized initially in DNA in direct damage processes (the former).
For further details in simulating such microscopic events at nano-scales, we refer the readers to our previous studies and publications, for example Ref. [Abolfath2013:PMB].

For the present discussion, it is critical to recognize that the locations of the sites of interaction are random, hence the path lengths and distances between two sequential interaction sites along the beam central axis, l, are random variable. This is in addition to randomness in energy imparted at the location of interaction, ε.

Although it should be obvious that both l and ε are two random independent variables, following two independent distribution functions, in the standard microdosimetry formalism, only ε was treated directly a random variable associated with the beam quality.
The sub-micrometer randomness in steps in l was neglected, because such information was not necessary to be recorded from the measurement theory standpoints.
Instead a constant chord-length that is a characteristic / geometrical length of microdosimetric sites was superseded.

More rigourously, the statistical fluctuations in l with regards to all collisions scored at nano-scale, crucial to DNA-damage statistical analysis, cannot be captured by micrometer-scale chord-length calculated by averaging l over the track-distribution functions.
In fact, in the microdosimetry distribution functions, the location of the collisions over a single track were already traced out / integrated over thus any information regarding to individual collision was wiped out.
In our approach, presented in this work, we substituted the standard microdosimetry distribution functions by collision distribution functions.
Because a single track consists of several collisions in a microdosimetry site, the collision-based distribution functions of ε and l are more informative and appropriate for modeling biological responses.


In addition to aforementioned shortcomings, substantial number of studies and publications, including the microdosimetric kinetic model (MKM) [Hawkins1998:MP,Hawkins2003:RR,Kase2006], implemented in the Particle and Heavy Ion Transport Code System
(PHITS) [Takada2018:JRR], formulated based on an alternative approach by introducing a concept of virtual spherical domains in cell nuclei with a size that varies within nano- to μ-meter.
In these models, it is not clear the rational for a choice of domain size but considering it spuriously an additional phenomenological parameter.
Therefore the application of current micro-dosimetry formulation in radio-biology does not seem appropriate and is prone to mixing the above discernible aspects.




We note that the Cauchy formula for the mean chord length is valid only for isotropic radiation (for non-spherical detectors) [Magrin2018:PMB].
Hadron therapy beams are, at the contrary, non-isotropic and the primary particles have essentially a unique direction.
Thus several microdosimetry publications were developed based on non-Cauchy algorithms for the mean chord length.
However, the extension of Cauchy to non-Cauchy calculation of chord length does not seem to overcome the shortcomings of application of chord length in biological models, because the chord-length, by definition, is a character of a surface boundary of microdosimetry sites, whether it is a detector, cell nucleus or MKM domains.
The chord-length as such does not contain detailed information on occurrence of the biological events and complexities inside the bulk of the site.

Thus it is appealing to revisit and refine these well established micro-dosimetry models and generalize them for modelings in biological response theories, radio-biological applications, and nano-dosimetry.
This is the main goal of the present study to propose necessary modifications essential in tweaking micro-dosimetry and suit it for nano-dosimetry formulations.

Along this line of thoughts, we propose substitution of the cellular geometrical chord-length (regardless if it was computed using Cauchy or non-Cauchy techniques) by the tracks diffusive mean free path, ⟨⟨ l ⟩⟩.



































We demonstrate that to simulate biologically relevant passage of a charge particle, it is necessary to collect the events, collision-by-collision, and combine them together to shape a track-structure, a bottom-up approach.
This led us to perform a double averaging over collisions and tracks to calculate the track mean-free path, denoted by ⟨⟨ l ⟩⟩.
Similar to a gas phase of matter, the molecules collide with one another, hence the diffusive mean free path is the average distance a particle travels between collisions. The larger the particles or the denser the gas, the more frequent the collisions are and the shorter the mean free path.
Turning back to original problem, the shorter the mean free path of the charged particle tracks manifests in higher complexities in DNA damage, hence the higher biological impact.
We thus supersede the track-averaged chord length, ⟨ l ⟩, by the mean-free path, ⟨⟨ l ⟩⟩.
In calculation of lineal energy and LET, to be consistent with changes in calculating ⟨⟨ l ⟩⟩, we also supersede the single event energy imparted, ε, by energy imparted in a single collision in the cellular SVOI. We then average over individual collisions and tracks in SVOI, wherever is necessary. 

We substantiate our proposal by presenting numerical illustrations on recently reported experimental data based on in-vitro clonogenic cell survival assay of non-small cell lung cancer (NSCLC) cells, i.e., an observation which was obtained by performing a high-throughput and high accuracy clonogenic cell-survival data acquired under exposure of the therapeutic scanned proton beams [Guan2015:SR].
We compare two RBE models using (1) ICRU 36 geometrical chord length of typical spherical cell nuclei with radius of ⟨ l ⟩ = 5 μ m and (2) charge particle diffusive mean-free path, ⟨⟨ l ⟩⟩ that is a variable and a function of depth.
We show that a monotonic decrease in ⟨⟨ l ⟩⟩ as a function of depth leads to a monotonic increase up to a factor of 4 in RBE.






§ METHOD AND MATERIALS



 §.§ Mathematical formalism

In our recent study [Abolfath2019:EPJD], we demonstrated that the linear density of DNA DSBs, Δ_l,
a fundamental biological response function in a radiation field, is a variable of linear density of the collisions.
To this ends, we recall the following relation derived in our first principle multi-scale study [Abolfath2019:EPJD]

    Δ_l = μ/m⟨⟨ε^2 ⟩⟩/⟨⟨ε⟩⟩1/⟨⟨ l ⟩⟩,

where ε and l are single collision energy imparted and stepping-length, calculated by Geant4 MC toolkit [Agostinelli2003:NIMA,Allison2006:TNS,Allison2016:NIMPRA].
Here μ and m are the average number of DSBs per deposition of 1 Gy of ionizing dose and cellular / DNA mass, respectively.


Hence, the relevant quantity in unit of lineal-energy (energy per length), consistent with Eq.(<ref>), was defined

    y_1D = ⟨⟨ε^2 ⟩⟩/⟨⟨ε⟩⟩1/⟨⟨ l ⟩⟩,

such that

    Δ_l = μ/m y_1D.


In terms of these quantities, calculation of cell survival was performed

    -ln(SF) = ∑_i=1^N ∑_j=1^i b_i,jΔ_l^i-j D^j.

Expanding Eq.(<ref>) and keeping the series up to quadratic term in dose, D, we found radio-biological α, β indices, as following

    α = ∑_i=1^N b_i,1Δ_l^i-1,

and

    β = ∑_i=1^N-1 b_i+1,2Δ_l^i-1,

The dependence of α and β on y_1D can be derived from the dependence of Δ_l on y_1D.
The coefficients in Eq.(<ref>) are identified by fitting to the experimental cell-survival data.



 §.§ Geant4 Monte Carlo simulations

In Geant4 each pencil proton beams were simulated by irradiating a cylinder water phantom with 20 cm radius and 20 cm length. The mean deposited energies ⟨⟨ε⟩⟩ and ⟨⟨ε^2⟩⟩, the mean track length ⟨⟨ l ⟩⟩ were scored within a linear array of voxels with 0.5 mm down to 5 μ m thickness.
In our notations, averaging over all single point collisions, following by averaging over all single primary and secondary particle tracks are denoted by double averaging, ⟨⟨⋯⟩⟩.
In MC we scored y_1D, y_D and LET_d from collision-by-collision following by event-by-event averaging and energy deposition ε_j and stepping length l_j, using the following identities
y_1D = (∑_j ε_j^2/∑_j ε_j)/(∑_j l_j/∑_j 1_j), y_D = ∑_j (ε_j/l_j)^2/∑_j (ε_j/l_j), and
LET_d = ∑_j (ε_j^2/l_j) / ∑_j ε_j
where sum over j includes all energy deposition events from primary and secondary processes in all steps and tracks in a specific voxel, hence ∑_j 1_j represents total number of single-collisions, scored in all energy imparted events.

Therefore y_1D=(⟨⟨ε^2⟩⟩/⟨⟨ε⟩⟩)/⟨⟨ l ⟩⟩ as well as other types of LET's were calculated.
The number of primary protons and the number of interactions per track were saved, in the same volume.  Then, the energy deposition, the track length, the number of primary proton and the number of interactions were accumulated, in each cell. All simulations used 10^6 protons with series of energy cut-off, corresponding to particle range that vary within 1 mm and 1 μ m.

In Geant4, any particle with energy below the cut-off value is assumed to not produce secondary particles (i.e. production threshold).
Below these cuts, the particle is transported further according to the CSDA approximation which will still imply a varying energy loss, i.e., no tracking cuts  [Agostinelli2003:NIMA].
All simulation results presented used the QGSP-BIC-EMY physics list.
We used Gaussian proton energy spectrums with very small FWHM (0.18 MeV). Because of small divergence the simulated beam is mono-energetic.







In Geant4, when primary particle (in our case proton) collide with other particle, secondary particles are generated.
This includes photon, electron, proton, neutron, He and heavier ions. The primary particle will not be processed until all secondaries are dealt with.
The importance of the energy cut-off, requires the program to stop producing more secondaries when the energy of the secondaries become lower than energy cut-off, defined by the user as a cut-off in length-scale.
Otherwise the simulation will not be practical as will take considerably long time to simulate the events.
However, the very low energy secondaries will not be killed, but it will follow the continuous slowing down approximation, CSDA.
The CSDA works fine with low energy particles as they tend to not travel far in the material.
In other words, the cut-off represents the accuracy of the stopping position, and any particle will always be tracked down to zero kinetic energy.
In our previous studies, we did validate the Geant4 simulation using the cut-off values (as we used in current study) and experimentally measured Bragg peak via production of Cerenkov light [Helo2014:PMB].
The difference was found to be less than 1% for cut-off of 0.01 mm.

In the current manuscript, we used the terms cut-off, and cut-all interchangeably.
The term cut-all used in the legend of figures to emphasize a specific cut-off value applied uniformly for “all" primary and secondary particles.
Note that the corresponding cut-off values in energy depends on the type of secondary particles.
For example the length-scale cut-off of 0.01 mm, approximately corresponds to 0.025 MeV cut-off for electrons in water.
Thus this value of energy for proton and the rest of particles is different than the corresponding value of cut-off in energy for electron.

Finally we note that the mean free path calculated by Geant4 methodology, as described above, should not be confused with the distance between individual ionizations, as this would be a true track structure approach. But from the Geant4 MC simulations, it must be clear what we calculated is a condensed history simulation where a large number of ionizations are grouped in a single energy loss step so it is rather the mean free path between individual energy loss events above a certain production cut.












§ RESULTS




In Fig. <ref>, we illustrate the effect of particle cutoff length on ⟨⟨ l ⟩⟩.
The calculation was performed for a pencil beam of proton with nominal energy 80 MeV, hitting a cylindrical water with radius and length of 20 cm, and slice thickness of 0.5 mm.
As shown in Fig. <ref>, the depth dependence of ⟨⟨ l ⟩⟩ is strongly influenced by the user defined cutoff choice.
The cutoff larger or equal to size of slice thickness, as shown by the black circles and red triangles respectively, yield the collision averaged chord-length equal to the slice thickness, a result that coincides with the geometrical chord-length of the water slice-thickness.
In this limit of large cut-off's, our numerical results recover the CSDA limits and the geometrical chord length calculated by Cauchy, Kellerer et al. and ICRU 36.
With lowering the cutoff down to 10th of slice thickness (green triangles) or even further (blue squares), the track mean-free path becomes significantly smaller than the geometrical slice thickness.
By lowering cutoff, below 0.01 mm, we observe negligible change in ⟨⟨ l ⟩⟩,
i.e., if we lower the cut-off below 0.01 mm, the generated curve exhibits slight difference relative to the curve corresponding to 0.01 mm cut-off, shown by squares in Fig. <ref>.
Numerically, we reached to a domain that resulted in convergence of the output data with respect to variations in cutoff values.
A continuous decrease in ⟨⟨ l ⟩⟩ as a function of depth reveals an increase in linear compactness of collisions.

Note that, for slice thickness of 0.5 mm, as in Fig. <ref>, and large cut-off values, no secondary particles are generated.
The primary particle, i.e., proton, loss energy according to CSDA and in results the simulated chord length turns out to be identical to the slice thickness, e.g., 0.5 mm as in Fig. <ref>, which fits very well with Rossi-Kellerer's theory of chord-length.
However, if we lower the cut-off value the simulated chord length will become considerably less (0.15 mm as seen in Fig. <ref>).



In Fig. <ref>, we show ⟨⟨ l ⟩⟩ (black lines) and relative energy deposition, ⟨⟨ε⟩⟩, (red lines) calculated for a pencil beam of 10^6 protons with nominal energy 80 MeV, traversing a cylindrical water phantom with slice thickness equal to 5.0 μ m.
To score such fine slices, we divided the water phantom thickness, equal to 20 cm, into 40,000 divisions and scored collision-by-collision energy deposition, ε, and stepping length, l, in each slice, using a cutoff values equal to 1 μ m.
The vertical dash line indicates the position of Bragg peak at approximate 5.2 cm depth corresponding with LET_d ≈ 10 keV/μ m.

Note that because of 5 μm slice thickness and use of 1 μm cut-off, the effect of delta-rays and secondary particles, imparting their energies outside of the slices they were generated, effectively alter the shape of PDD. However, the location of the Bragg peak does not change as the slice thickness reduces from 1 mm to 5 μm.
We should also remark that if we lower the thickness of the SVOI to nano-scales, we must lower the cut-off values to below 1 nm to obtain sensible results.


It is intriguing to note that from beam entrance to the end of the proton range, ⟨⟨ l ⟩⟩ decreases by a factor of eight.
Coincidentally, a factor of 4 has been reported experimentally in increase in proton RBE distal and proximal to Bragg peak (for more discussion see Fig. <ref>, below).
The wiggling lines in the vicinity of proton range stem from statistical noises due to particle energy straggling and lack of enough statistics because of small number of particle fluence.




In Fig. <ref>, y_1D = [⟨⟨ε^2 ⟩⟩ / ⟨⟨ε⟩⟩] / ⟨⟨ l ⟩⟩ vs. LET_d = ⟨⟨ε^2/l ⟩⟩ / ⟨⟨ε⟩⟩ is shown.
The trace shown by black dots obtained from variable ⟨⟨ l ⟩⟩, as shown in Fig. <ref> whereas the trace by the red dots uses an ICRU 36 type of chord-length estimate, i.e., ⟨⟨ l ⟩⟩ = 5 μ m everywhere, regardless of the depth.
Note that ⟨⟨ l ⟩⟩ = 5 μ m is a geometrical chord-length for a slab with thickness equal to 5 μ m.

The green solid line is a CSDA line, implying LET_d = y_1D under assumption ⟨⟨ 1/l ⟩⟩≈ 1 / ⟨⟨ l ⟩⟩.
A similar non-linear dependencies between different types of LET's were previously presented in [Abolfath2019:EPJD].
A non-linear dependence between LET_d and y_1D suggests advantageous in using y_1D over LET_d for RBE studies because of linear dependence between y_1D and Δ_l, as is given by Eq.(<ref>).

We note that the scattered dots in Fig. <ref> are due to the statistical fluctuations in a MC simulation performed under limited number of primary protons. By increasing number of protons, the scattered dots merge together and form a single curve for each value of cutoff as shown in Fig. <ref>.
In this calculation we used 10^6 protons. Because of small thickness of slices and small cut-off value, the MC calculation is very slow. To generate similar curves with lower statistical variations, a simulation based on larger number of protons, e.g., 10^7 protons is required. However, because in Fig. <ref>, we only intend to show a non-linear dependence between LET_d and y_1D, the envelope dependence obtained from scoring a million of protons is adequate for our current qualitative analysis. We note that, the red and black envelop-curves in Fig. <ref> were collected through merging the scattered points. We also note that we did not perform any analytical fitting procedure to derive a specific curve, to investigate the analytical dependencies of LET_d on y_1D.




We now turn to present RBE of in-vitro clonogenic cell survival assay of H460 cell lines, a type of non-small cell lung cancer (NSCLC).
The dependence of SF on dose and LET were measured experimentally in our group [Guan2015:SR]
and fitted by 3D global fitting method by the present authors [Abolfath2013:PMB,Abolfath2016:MP,Abolfath2017:SR,Abolfath2019:EPJD].
Fig. <ref> shows the result of RBE calculated using variable ⟨⟨ l ⟩⟩ (red dots) and constant ⟨⟨ l ⟩⟩ = 5μ m (black dots),
corresponding a biological endpoint with 10% cell-survival fraction, as reported by Guan et al. [Guan2015:SR].
The numerical values of the coefficients in the polynomial expansion of α and β on Δ_l and LET are given by Eqs.(<ref>-<ref>).
For LET ≤ 10 keV/μ m we obtained the optimal fitting to the following polynomials with
N=2, b_1,1 = 0.05237, b_2,1 = 0.0151 and b_2,2 = 0.03265, corresponding to reduced χ^2 = 0.00541, R^2(COD)=0.98233 and adjusted R^2 = 0.98233.
This is a limit that dependence of α on Δ_l and LET is linear whereas β shows no dependence on Δ_l (and LET),
Note that this is a limit frequently used in literature for all range of LET
(see e.g., Eq.(II.28) in Ref. [Hawkins1998:MP] or Eq. (8) in Ref. [Kase2006]).
For LET > 10 keV/μ m, we obtained N=6 with
b_1,1 = 0.12288, b_2,1 = 0.00624,
b_3,1 = 1.322× 10^-6, b_4,1 = 9.449× 10^-7,
b_5,1 = 8.812× 10^-7,
b_6,1 = 7.725× 10^-8, and
b_2,2 = 0.00382, b_3,2 = 0.0007821, b_4,2 = 9.275× 10^-5, b_5,2 = 9.511× 10^-7,
b_6,2 = 1.68× 10^-7
corresponding to
reduced χ^2 = 0.05472, R^2(COD)=0.96777 and adjusted R^2 = 0.93911.



The latter is a method of calculation, recommended by ICRU 36.
A comparison between two methods (a) considering variability in ⟨⟨ l ⟩⟩ as a function of depth by double averaging over collisions and tracks and (b) considering geometrical chord length equal to the cell line thickness, 5μ m, evidently exhibits significant difference in predicting RBE by a factor of 4 in domains distal and proximal to the Bragg peak.

Similar variabilities in RBE as a function of proton range were reported in other experimental works using conventional microdosimetry approach and a methodology based on an empirical or phenomenological “biological waiting functions" in calculation of RBE [DeNardo2004:RPD].
Accordingly, the biological waiting functions fitted to the spectrum of a spread-out Bragg peaks (SOBP) beam of proton collected by a 2.3 mm a microdosimetric prob, using a tissue-equivalent proportional counter (TEPC), resulted in a monotonic increase in RBE up to 2.5, a value close to the RBE reported in this work.
We note that the proposed modifications presented in this work, is free from such phenomenological convolution between the beam spectrum and RBE.



§ DISCUSSION AND CONCLUSION

There are substantial evidence that proton (similar to heavier charged particles) RBE increases as a function of depth.
This is partly due to increase in spatial density of collisions as protons pass through tissue and lose energy at an increasing rate.
To make connection between experimental radiobiological data and mechanistic models, it is customary approach in literature to use microdosimetry models, developed by Rossi, Kellerer and colleagues [Kellerer1971:RR,Kellerer_Rossi1972:CTRR,Kellerer1984:RR,Kellerer1985:Book,Kellerer1975:REB].
Embedded in these models, geometrical chord length of a SVOI plays a crucial role.
We discussed that although this quantity is a good parameter to describe microdosimetry processes and charge collections in tissue proportional chambers, but because it is intrinsically character of the surface and geometrical boundaries of cell nuclei, it cannot directly describe nano-meter scale localities in stochastic microscopic physio-chemical processes in DNA materials and DNA-damage.
Thus we proposed the diffusive mean-free path length of the particle tracks as a new metric that describes appropriately variabilities in collisions compactness as a function of depth in tissue to substitute the geometrical chord-lengths considered in micro-dosimetry modelings.
This variation of micro-dosimetry is potentially more appropriate for radio-biological studies.

Acknowledgement:
The authors would like to acknowledge useful discussion and scientific exchanges with Drs. Alejandro Carabe-Fernandez and Alejandro Bertolet Reina.
The work at the University of Texas, MD Anderson Cancer Center was supported by the NIH / NCI under Grant No. U19 CA021239.

Authors contributions:
RA: wrote the main manuscript, prepared figures, performed mathematical derivations and computational steps including Geant4 and Geant4-DNA Monte Carlo simulations and three dimensional surface fitting to the experimental data.
YH: contributed to Geant4 Monte Carlo simulations and writing the manuscript.
DC, RS, DG and RM: wrote the main manuscript, contributed to scientific problem and co-supervised the project.

Corresponding Authors:

^† ramin1.abolfath@gmail.com / Ramin.Abolfath@pennmedicine.upenn.edu 

^* rmohan@mdanderson.org

























































99

ICRU36
ICRU 1983 Microdosimetry ICRU Report 36 (Bethesda, MD: ICRU).

Kellerer1971:RR
Kellerer A M, Considerations on the Random Traversal of Convex Bodies and Solutions for General Cylinders1, Radiat. Res. 47, 359-376 (1971).

Kellerer_Rossi1972:CTRR
A. M. Kellerer, H. H. Rossi, Curr. Top. Radiat. Res. Q. 8, 85 (1972).

Kellerer1984:RR
Kellerer A M, Chord length distributions and related quantities for spheroids Radiat. Res. 98, 425–37 (1984).

Kellerer1985:Book 
A. M. Kellerer, Fundamentals of microdosimetry The Dosimetry of Ionizing Radiation vol. 1 ed Kase, K. R., et al. (London: Academic) pp 77–161 (1985).

Kellerer1975:REB
A. M. Kellerer, D. Chmelevsky, Rad. and Environm. Biophys. 12, 205 (1975).

Cauchy1908
A. Cauchy, Oeuves Completes, Vol. 2. Gauthier-Villars, Paris (1908).

Magrin2018:PMB
G. Magrin, A method to convert spectra from slab microdosimeters in therapeutic ion-beams to the spectra referring to microdosimeters of different shapes and material, Phys. Med. Biol. 63, 215021 (2018).

Birkhoff1970:HP
R. D. Birkhoff, J. E. Turner, V. E. Anderson, J. M. Feolap and R. N. Hamm,
The determination of let spectra from I. track-length distributions in cavities,
Health Physics Pergamon Press 18, 1-14 (1970).

Agostinelli2003:NIMA
S. Agostinelli et. al., Nucl. Instrum. Meth. A 506, 250 (2003).

Allison2006:TNS
J. Allison et. al., Geant4 developments and applications IEEE Trans. Nucl. Sci. 53, 270–8 (2006).

Allison2016:NIMPRA
J. Allison et. al., Recent developments in Geant4 Nucl. Instrum. Methods Phys. Res. A 835, 186–225 (2016).

Incerti2010:IJMSSC
S. Incerti S. et. al., Int. J. Modelling Simul. Sci. Comput. 1, 157 (2010).

Conte2012:NJP 
V. Conte, P. Colautti, B. Grosswendt, D. Moro, L. De Nardo, Track structure of light ions: experiments and simulations, New J. Phys. 14, 093010 (2012).

Abolfath2011:JPC
Abolfath R M, van Duin A C T and Brabec T, Reactive molecular dynamics study on the first steps of DNA damage by free hydroxyl radicals, J. Phys. Chem. A 115 11045 (2011).
See the real-time simulations and movies at:
http://qmsimulator.wordpress.com/

Abolfath2013:PMB

Abolfath RM, Carlson DJ, Chen ZJ, Nath R. A molecular dynamics simulation of DNA damage induction by ionizing radiation, Phys. Med. Biol. 58, 7143-7157 (2013).

Abolfath2016:MP

Abolfath, R., Bronk, L., Helo, Y., Schuemann, J., Titt, U., Grosshans, D. and Mohan, R., WE-H-BRA-01: BEST IN PHYSICS (THERAPY): Nano-Dosimetric Kinetic Model for Variable Relative Biological Effectiveness of Proton and Ion Beams. Med. Phys. 43, 3842 (2016).

Abolfath2017:SR

Abolfath R, Peeler CR, Newpower M, Bronk L, Grosshans D, and Mohan R,
A model for relative biological effectiveness of therapeutic proton beams based on a global fit of cell survival data,
Sci. Rep. 7, 8340 (2017).

Abolfath2019:EPJD
R. Abolfath, Y. Helo, L. Bronk, A. Carabe, D. Grosshans and R. Mohan Renormalization of radiobiological response functions by energy loss fluctuations and complexities in chromosome aberration induction: deactivation theory for proton therapy from cells to tumor control, European Physical Journal D 73, 64 (2019), DOI: 10.1140/epjd/e2019-90263-5

NumericalRecipes
W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipes in Fortran 77: the Art of Scientific Computing. Second Edition, vol. 1 (1996).

Bolst2017:PMB
Bolst D, Guatelli S, Tran L, Chartier L, Lerch M, Matsufuji N and Rosenfeld A, Correction factors to convert microdosimetry
measurements in silicon to tissue in 12C ion therapy Phys. Med. Biol. 62, 2055 (2017).

Anderson2017:MP
S. E. Anderson, K. M. Furutani, L. T. Tran, L. Chartier, M. Petasecca, M. Lerch, D. A. Prokopovich, M. Reinhard, V. L. Perevertaylo, A. B. Rosenfeld, M. G. Herman, C. Beltran, Microdosimetric measurements of a clinical proton beam with micrometersized solid-state detector, Med. Phys. 44, 6029 (2017).

Bolst2018:PMB
D. Bolst, S. Guatelli , L. T. Tran, A. B. Rosenfeld
Optimisation of the design of SOI microdosimeters for hadron therapy quality assurance
Phys. Med. Biol. 63 215007 (2018).

Carlson2008:RR
D. J. Carlson, R. D. Stewart, V. A. Semenenkoa, G. A. Sandison,
Combined Use of Monte Carlo DNA Damage Simulations and Deterministic Repair Models to Examine Putative Mechanisms of Cell Killing
Radiat. Res. 169, 447–459 (2008).

Hawkins1998:MP
R. B. Hawkins, Med. Phys. 25, 1157 (1998).

Hawkins2003:RR
R. B. Hawkins, Radiat. Res. 160, 61 (2003).

Kase2006
Y. Kase, et al. Radiat. Res. 166, 629 (2006).

Bianco2015:RRD
D. Bianco, C. Villagrasa, M. Dos Santos, Multi-scale analysis of simulated proton and alpha irradiation, Radiat. Protect. Dos. 164, 444–448 (2015).

Takada2018:JRR
K. Takada, T. Sato, H. Kumada, J. Koketsu, H. Takei, H. Sakurai, T. Sakae,
Validation of the physical and RBE-weighted dose estimator based on PHITS coupled with a microdosimetric kinetic model for proton therapy.
J. Radiat. Res. 59, 91-99 (2018). 

Guan2015:SR 
F. Guan, L. Bronk, U. Titt, S. H. Lin, D. Mirkovic, M. D. Kerr, X. R. Zhu, J. Dinh, M. Sobieski, C. Stephan, C. R. Peeler, R. Taleei, R. Mohan, D. R.  Grosshans, Sci. Rep. 5, 1 (2015).

Helo2014:PMB
Y. Helo, A. Kacperek, I. Rosenberg, G. Royle, A. P. Gibson, The physics of Cerenkov light production during proton therapy,
Phys. Med. Biol. 59, 7107-23 (2014).

DeNardo2004:RPD
L. De Nardo, D. Moro, P. Colautti, V. Conte, G. Tornielli, G. Cuttone
Microdosimetric investigation at the therapeutic proton beam facility of CATANA,
Radiat. Protect. Dos., 110, 681-686 (2004).







