Calculation of biological effectiveness of SOBP proton beams: a TOPAS Monte Carlo study

Objective. This study aims to investigate the biological effectiveness of Spread-Out Bragg-Peak (SOBP) proton beams with initial kinetic energies 50–250 MeV at different depths in water using TOPAS Monte Carlo code. Approach. The study modelled SOBP proton beams using TOPAS time feature. Various LET-based models and Repair-Misrepair-Fixation model were employed to calculate Relative Biological Effectiveness (RBE) for V79 cell lines at different on-axis depths based on TOPAS. Microdosimetric Kinetic Model and biological weighting function-based models, which utilize microdosimetric distributions, were also used to estimate the RBE. A phase-space-based method was adopted for calculating microdosimetric distributions. Main results. The trend of variation of RBE with depth is similar in all the RBE models, but the absolute RBE values vary based on the calculation models. RBE sharply increases at the distal edge of SOBP proton beams. In the entrance region of all the proton beams, RBE values at 4 Gy i.e. RBE(4 Gy) resulting from different models are in the range of 1.04–1.07, comparable to clinically used generic RBE of 1.1. Moving from the proximal to distal end of the SOBP, RBE(4 Gy) is in the range of 1.15–1.33, 1.13–1.21, 1.11–1.17, 1.13–1.18 and 1.17–1.21, respectively for 50, 100, 150, 200 and 250 MeV SOBP beams, whereas at the distal dose fall-off region, these values are 1.68, 1.53, 1.44, 1.42 and 1.40, respectively. Significance. The study emphasises application of depth-, dose- and energy- dependent RBE values in clinical application of proton beams.


Introduction
Proton therapy is one of the fast-growing treatment modalities in radiotherapy.One of the important features of proton therapy is its inherent tissue-sparing capabilities which reduces the burden of treatmentrelated complications on patients and the healthcare system (Newhauser and Zhang 2015).The potential advantage of using proton beams for clinical purpose is related to their distinct energy deposition characteristics.The Bragg curve, generated by monoenergetic protons, exhibits a well-defined, localized peak at the end of the proton track (Raju et al 1978).This results in low dose deposition to the healthy tissue which leads to a reduction in normal tissue toxicity (Lühr et al 2018).Moreover, proton therapy presents the additional benefit of lower integral dose deposition when compared to modern photon-based therapies such as Intensity Modulated Radiation Therapy or Volumetric Modulated Arc Therapy.This characteristic underscore the potential of proton therapy to reduce the risk of secondary neoplasms to develop (Lühr et al 2018).Adhering to clinical requisites, the Spread-out Bragg peak (SOBP) can be generated by modulating the proton energy, ensuring the uniform delivery of the prescribed dose across the targeted volume (Conte et al 2019).
The biological effectiveness of ionizing radiation depends on the Linear Energy Transfer (LET) of the interacting particles (Conte et al 2019).In contrast to low-LET particles, high-LET particles produce clusters of interactions in few of the irradiated cells (Conte et al 2019).These clusters give rise to highly complex patterns of damage.In the context of particles other than electrons and photons, the absorbed dose needs to be weighted by the appropriate Relative Biological Effectiveness (RBE) (Conte et al 2019).
One of the major concerns in proton therapy treatment planning, especially for sites closely located to organs at risk, is the uncertainties associated with the RBE value of clinical proton beams (Koh et al 2022).Upon interaction with matter, high-energy protons give rise to a composite radiation field comprising diverse particles, each possessing varying LET (Conte et al 2019).At any given depth, both primary and secondary particles produced from non-elastic nuclear interactions of protons contribute to the biological damage of the living systems (Paganetti 2002).Consequently, the resulting RBE of this composite radiation field reflects the amalgamation of RBE values attributed to multiple constituents of the resultant mixed radiation spectrum (Conte et al 2019).In current clinical practice, proton therapy treatments are planned assuming a constant RBE value of 1.1 irrespective of beam energy, treatment depth and biological endpoints (Paganetti 2014, Qutub et al 2016, Newpower et al 2019).However, the dependency of RBE on factors like LET, cell or tissue type, and dose per fraction leads to fluctuations in RBE along the depth-dose profile (Giovannini et al 2016, Qutub et al 2016).Lühr et al (2018) extensively explored the existing understanding of the upcoming challenges pertaining to the RBE in the context of proton beam therapy.Qutub et al (2016) introduced a straightforward algorithm for risk assessment in RBE-weighted treatment planning across a two-dimensional plane, and they compared the outcomes with those derived from the conventional RBE assumption of 1.1.
Several radiobiological experimental studies reported an increasing trend of RBE along the SOBP and the RBE values vary significantly from the constant value of 1.1 (Paganetti et al 2002, Chaudhary et al 2014, Paganetti 2014, Conte et al 2019).Paganetti et al (2002) experimentally showed that the RBE of SOBP proton beams (65-250 MeV) varies from 0.9-2.1 for in-vitro scenarios, whereas in-vivo RBE values vary from 0.7-1.6 (Paganetti et al 2002).Chaudhary et al (2014) measured variation in cell killing RBE for both Human fibroblasts (AG01522) and glioma (U87) cells across varying depths, employing 62 MeV monoenergetic and SOBP proton beams.The authors then compared their experimental findings against that calculated using Geant4 Monte Carlo toolkit coupled with Local Effect Model (LEM).The outcomes of the study revealed a significant variance between the predicted biological dose, which was calculated using the variable RBE extracted from the experimental data, and the conventionally utilized constant RBE value of 1.1 in clinical practice (Chaudhary et al 2014).The RBE value derived from radio-biological experiments is a precisely defined parameter and its applicability remains confined to the specific experimental circumstances.Directly extrapolating this value to patient scenarios is generally not possible (Fossati et al This present study endeavours to compute the RBE of SOBP proton beams as a dynamic function of depths in water, utilizing microdosimetry-driven RBE models.These computed RBE values are then compared with the RBE values derived from various LETbased models and the Repair-Misrepair-Fixation (RMF) model.To realize this goal, microdosimetric distributions at 1 μm site size were calculated using TOPAS code.This study also included calculation of LET d profiles with respect to depths in water for SOBP proton beams.Present study considered 50, 100, 150, 200 and 250 MeV SOBP proton beams which covers entire range of clinical proton beams.

TOPAS Monte Carlo code
Present study used TOPAS Monte Carlo code (version: 3.6.1)along with its 'microdosimetric' and 'proto-nRBE' extensions for calculating microdosimetric distributions and RBE of clinical proton beams (Perl et al 2012).TOPAS is a Monte Carlo platform layered on the top of Geant4 toolkit (Agostinelli et al 2003, Allison et al 2006).The lineal energy scorer used to score microdosimetric distributions in the spherical TEPC is developed by Underwood et al (2017. Polster et al (2015) discussed about the implementation and validation of 'protonRBE' extension of TOPAS.All the TOPAS-based simulations used default TOPAS physics list which activates g4em-standard_opt3, g4h-phy_QGSP_BIC_HP, g4decay, g4ion-binarycascade, g4h-elastic_HP, g4q-stopping, and g4radioactivedecay.The electron range cut value was set at 100 nm.

Modelling of SOBP proton beams
For irradiating an extended volume uniformly to a certain dose, SOBP is created by the superposition of several Bragg peaks with different range and thus, initial energy of the protons used to form the SOBP.The weighting of each peak or beamlet takes the proximal dose tail of its next higher energy Bragg curves into account (Velten and Tomé 2020).The SOBP proton beams were modelled using the algorithm developed by Velten and Tomé (Velten and Tomé 2020).The key points of the algorithm are discussed here.The relation between range and initial energy of a proton in a medium can be approximated using a power law: 0 where a and p 0 depend on the medium and to a lesser extent on proton energy.For water medium, a = ´--cmMeV 2.2 10 3 1 and » p 1.77 0 (Bortfeld 1997).The SOBP extends from where c is the maximum R .0 For the k-th beamlet of SOBP, range (r k ), energy (e k ), and weight (w k ) were calculated using the modified formulas of Jette and Chen (Jette and Chen 2011) as given below: In TOPAS, w k corresponds to the amount of the time a given beamlet energy is being used for the simulation.Thus, beam on time for the i-th beamlet can be calculated as per the equation and total beam on time, T = 1000 ms is considered.For this purpose, time features of TOPAS code is used.The number of beamlets (n) varies based on SOBP width using a uniform range spacing of 0.5 mm.Velten and Tome used brute-force approach for optimization of the power law parameter (p) yielding flat SOBPs of protons (E 0 = 50-250 MeV) in water (Velten and Tomé 2020).The value of p varies depends on E 0 and c.The optimized values of p are reported by Velten and Tomé (2020).This algorithm is published as a python script as an extension of TOPAS code (https://github.com/cvelten/TOPAS-SOBP).The initial kinetic energies (E 0 ) and SOBP width (W SOBP ) of the SOBP proton beams considered in the present study are presented in table 1. W SOBP is defined as the fraction of the maximum range of the protons.Note that there are several definitions for SOBP width (Levin et al 2005).However, as per ICRU report 78, SOBP width is defined as the distance between the depth of 94% of plateau dose (flat region of the SOBP) at the proximal and distal side (Newhauser 2009).

Calculation of depth-dose profile
In this study, a parallel proton beam with a beam radius of 5 cm was incident on a 40 × 40 × 40 cm 3 unit density liquid water phantom.The proton beams entered from vacuum and incident on the surface of the water phantom.Note that most of the Monte Carlo-based published studies used parallel SOBP proton beams (Böhlen et al 2011, Paganetti 2017, Vassiliev et al 2019, Velten and Tomé 2020) and for this purpose, SOBP proton beams as mentioned in para 2.3 were considered.Depth along the central axis was voxelized into cylindrical voxels of dimensions 1 cm × 2 mm (radius × height).The absorbed dose to water was scored in each voxel using the TOPAS Monte Carlo code.The absorbed dose at each depth was normalized with respect to the dose at the depth of dose maximum.This normalized percentage dose was plotted against depth in water.
2.4.Calculation of phase space in water phantom TOPAS-based microdosimetric calculation of SOBP proton beams as a function of depth in water phantom is a two-step process: (a) calculation of phase-space at the depth of interest in water, and (b) calculation of microdosimetric distributions in the LET-1/2 TEPC using the phase-space as a source.In the calculations, parallel SOBP proton beams of radius 5 cm were incident on a 40 cm × 40 cm × 40 cm water phantom.Thin boxes having dimension 2 cm × 2 cm × 10 μm filled with water were placed at the depths of interest along the central axis in water phantom.The schematic simulation geometry is shown in figure 1.For a given SOBP beam, present study calculated phasespaces at five depths in water.Figure 2 shows the positions of these calculation depths on a typical SOBP beam.These depths represent entrance region (P), dose-gradient region (Q), proximal (R), central (S) and distal end (T) of the SOBP region.The regions T -U of figure 2 represents distal dose fall-off region.In addition, phase-spaces were also scored for 50 and 250 MeV monoenergetic proton beams (beam radius = 5 cm) at selective depths in water.The phasespace was scored across the 'ZplusSurface' of the thin box.Protons (primary and secondary) and secondary particles such as neutrons, alphas, electrons, and gammas were stored in the phase space.The phasespace file stored the position, particle type, energy and momenta of all the particles.The original number of simulated primaries for each energy component of the SOBP beam is 4 × 10 8 .The phase-space data were saved in a binary format.Note that the co-ordinates of the phase-space are always with respect to the WORLD volume.The electron range cut was set at 100 nm in water phantom.In clinical proton energy range, water closely mimics the properties of human tissues in terms of energy loss, Multiple Coulomb Scattering and nuclear interactions (Newhauser and Zhang 2015).As per IAEA andICRU reports (Kacperek 2000, Radiology andESTRO. 2000), recommended phantom material for dose and range measurements is water.
In order to get confidence on the simulation parameters and physics models used in TOPAS-based calculations, auxiliary simulations have been carried out to compare with the measured microdosimetric distributions of 62 MeV SOBP proton beams by (Conte et al 2019) at the CATANA proton therapy centre of the Southern National Laboratories of INFN in Catania.At this facility, the initial energy of the protons was 62 MeV, which underwent energy degradation as they traversed different elements of the beam line.As a result, the effective energy of the 62 MeV protons at the exit point of the CATANA proton beam line was approximately 53.1 MeV which was directed onto the rotating modulator wheel to generate the SOBP (Chiriotti Alvarez 2015).Thus, SOBP proton beams with initial kinetic energy of 53.1 MeV was modelled in TOPAS using the algorithm described in section 2.3.and phase spaces were calculated at d = 1.4,24.6 and 29 mm depth in water.

Microdosimetric distributions
The microscopic distribution of the energy deposition by a single particle is important for determining the biological effect (Conte et al 2019).Microdosimetry evaluates stochastic distribution of energy depositions at cellular and sub-cellular targets (Rossi 1979, Griffiths 1985, Zaider et al 1996).Mcrodosmetric Kinetic Model (MKM) and biological weighting functionbased biophysical model used microdosimetric distributions as input parameters for determining RBE of ionizing radiations (Kellerer and Rossi 1974, Rossi 1979, Griffiths 1985, Zaider et al 1996).The lineal energy (y) is an important quantity in microdosimetry and is defined as l, where e is the imparted energy by a single energy deposition event within the microscopic volume and l is the mean chord length (Griffiths where d is site The on-axis microdosimetric distributions at 1 μm site size for SOBP proton beams as a function of depth along the central axis (d) were calculated using pre-calculated phase-space as source and TEPC in vacuum.In TOPAS, the co-ordinates of the phasespace are defined with respect to the 'WORLD' volume.In the simulation set up, the centre of the LET-1/2 TEPC should be at least 0.762 cm (external radius of the TEPC) away from the centre of the phasespace to avoid overlap.Note that present study included simulation of 62 MeV SOBP proton beams for comparing TOPAS calculated ̅ y F and ̅ y D at selective depths to the corresponding measured values by Conte et al (2019).Additionally, auxiliary simulations were carried out for calculating microdosimetric distributions of 100 MeV SOBP protons at d = 5 and 7 cm using full simulation set up and compared with the corresponding distributions calculated using the phase-space as described above Full simulation setup consists of parallel beam of protons (beam radius = 5 cm) and TEPC placed at d = 5 and 7 cm in 40 cm × 40 cm × 40 cm water phantom.This simulation is considered to ensure that there is nothing wrong with the phase-space model.r d is the domain radius and a 0 is the slope in the limit of LET = 0 The saturation-corrected dose-mean lineal energy, * y , is defined by the equation: y 0 is the saturation parameter, which is used to correct the overkilling effect for high LET radiation.The y 0 can be calculated using the following formula (Kase et al 2006) Where, R n is the nuclear radius.The parameters of the MKM for V79 cells used in the present study are: b= b x = 0.02 Gy −2 , r d = 0.26 μm, R n = 4.1 μm, a 0 = 0.184 Gy −1 , r t = 1 g cm RBE of the proton beams is defined as the ratio of isoeffect doses of a reference low LET radiation D x to the dose of the proton beams D p .As per LQ model (based on the parameters α and β), RBE is given as

Repair-misrepair-fixation
The details formalism of the repair-misrepair-fixation (RMF) model is described by Carlson et al (2008) and Frese et al (2012).In this model, the effects of particle type and kinetic energy on α and β are explicitly linked to the initial formation of DNA double strand breaks (DSBs), i.e., in the limit when the dose is small compared to α/β.As per this model, Where, S: number of DSB in Gy −1 Gbp −1 (or per cell); ̅ z : F frequency-mean specific energy in cell nucleus of diameter, d. q and k are two cell-specific adjustable parameters that relate to the biological processing of DSBs into more lethal forms of damage, such as chromosome aberrations.As per first order approximation, @ ̅ z 0.204 In TOPAS, α and β parameters of equation (v) are estimated using the following equations, where RBE DSB is defined as the ratio of S for the protons to that for the 60 Co reference radiation.
The parameter S required in the RMF model were determined in Monte Carlo Damage Simulations for a nuclear diameter of 5 μm.The present study considered 6.257 Gbp in the DNA of a V79 cells.S = 8.32 Gy −1 Gbp −1 is considered for reference radiations which represents normoxic V79 cell (Carlson et al 2008).A look-up table of the DSBs depending on the kinetic energy of the protons is included in the cell line file (Polster et al 2015).
In the present study, 40 × 40 × 40 cm 3 water phantom was voxelized along the central axis depth and the dimension of each voxel was 4 × 4 × 2 mm 3 .
Proton beams (50-250 MeV SOBP and, 50 and 250 MeV monoenergetic) were incident on the phantom surface.At first, the RBE scorers calculate the dose and LET d separately, and after completion of the simulation, the function 'CombineSubScorers' are called to combine these separate quantities on a voxel-by-voxel basis to calculate RBE or survival fraction as a function of depth.The functions listed in models 1-5 are applied to convert the calculated dose and LET d maps into RBE.The parameter 'Sc/SimultaneousExposure = 'False' was set which appropriately represents the cell experiments, where each scorer bin calculates RBE for a separate irradiation.That is, the simulation of the cell experiments is repeated at each depth using same prescribed dose and TOPAS has capability to simulate all these experiments in a single run, by normalizing the dose at each depth to the prescribed dose appropriately.In TOPAS, the LET d scorer calculates the LET d in the unit of The present study considered D p = 2 and 4 Gy for calculating RBE because the typical prescribed dose in proton therapy for deep sited lesions is in the range of 1.8-4 CGE/fraction where 'CGE' represents 'Cobalt Gray Equivalent' (Paganetti et al 2002).RBE at D p = 4 Gy is calculated using all the models.However, RBE at D p = 2 Gy is calculated using selective models such as MKM and model by Wedenberg et al (2013).Hereafter, the RBE at D p = 4 and 2 Gy will be referred as RBE(4 Gy) and RBE(2 Gy), respectively.

Depth-dose profile
The SOBP dose distribution is modelled by summing up the contributions from individually modulated pristine Bragg peaks (Levin et al 2005).Figure 3 presents central axis depth-dose profiles of SOBP proton beams in water.The width of SOBP (as derived from the depth-dose profile) for 50, 100, 150, 200 and 250 MeV SOBP proton beams are 0.9, 3, 5.4, 5.2, 5.7 cm, respectively.Beyond the SOBP region, the protons will stop due to which the dose will sharply drop to zero (Jäkel 2020) and this region is referred as distal dose fall-off region.The depths at which the microdosimetric distribution (positions P, Q, R, S and T as shown in figure 2) to be calculated are determined from the depth-dose profile of each SOBP beam.The values of P, Q, R, S and T in water phantom along the central axis for the considered SOBP proton beams are presented in table 2.

LET d versus depth
Figure 4 presents on-axis TOPAS-calculated LET d values as a function of depth in water for SOBP proton beams.At a given depth, LET d is higher for smaller proton energy.For a given proton beam, LET d gradually increases with depth, reaches a maximum value and then falls-off rapidly.LET d reaches maximum beyond the peak position i.e. at the distal dose fall-off region.For example, for 50 and 250 MeV SOBP proton beams, the maximum values of LET d are 27.4 keV/μm at d = 2.7 cm and 12.5 keV/μm at d = 39.1 cm, respectively.The LET increases with the decrease of proton energy.At the distal dose fall-off region, energy of proton is smaller and hence higher LET d .Within the SOBP region, the LET d values increase gradually with depth.For a given beam, within the flat region of the SOBP beam (see 'R-T' of figure 2), LET d value varies significantly.As compared to the proximal end, the LET d value at the distal end of the SOBP region is higher by a factor of 2.6, 3.6, 2.6, 2.3 and 2, respectively for 50, 100, 150, 200 and 250 MeV SOBP beams.Except at distal dose fall-off region, the variation of LET d with depth shows similar trend as observed by Wilkens and Oelfke (2003).Note that the analytical model by Wilkens and Oelfke (2003) is based on approximation of stopping power tables by power-law which neglects nuclear reaction.On the contrary, TOPAS includes nuclear reaction.Due to the nuclear reaction, the primary protons can generate neutrons which can travel much further than primary beam and these neutrons in turn generate secondary protons at far behind the peak of depth-dose curve.As the 'ProtonLET' scorer of TOPAS considers both primary and secondary protons as well as secondary electrons (Perl et al 2012), the LET can be calculated even at the distal dose fall-off region where only few secondary protons are available which have negligible contribution to the total dose.Hence, TOPAS-calculated LET d falls-off of after reaching maximum value.

Microdosimetric distributions
The microdosimetric distribution is the plot of yd(y) on a linear scale versus y on a log-scale.The distribution is normalized such that ò = ( ) d y dy 1.In this type of plot, equal areas under the curve represent equal doses delivered at the intervals of y-values.The present study used 20 lineal-energy bins per decade for plotting the distributions.

Comparison with published data
Table 3 presents a comparison between TOPAScalculated ̅ y F and ̅ y D values with the corresponding

Microdosimetric distributions of SOBP beams
Figure 5 presents on-axis microdosimetric distribution of SOBP proton beams at different depths in water.For a given beam, the proton energy decreases with depth and hence the distribution shifts towards the higher y value.The microdosmetric distribution at the entrance region (d = 0.5 cm) of 50 MeV SOBP protons shows that the majority of absorbed dose is delivered by events of size about 2.4 keV/ μm (peak position) which corresponds to the stopping power of 22 MeV protons in water.However, at the entrance region (position 'P' in figure 2), the microdosimetric distributions of 100-250 MeV SOBP proton beams show double maxima and the positions of these maxima vary depending on the beam energy.For example, positions of these maxima at d = 2 cm are: (a) y = 0.8 and 2 keV/μm for 100 MeV SOBP protons, and (b) y = 0.34 and 1.7 keV/μm for 250 MeV SOBP protons.To investigate the origin of this double maxima in the microdosimetric distribution, particlewise relative contributions to the total microdosimetric distributions are calculated for 100 and 250 MeV SOBP protons at the entrance region (d = 2 cm) and shown in figure 6.This figure shows  that both the peaks are due to protons.As y = 0.34 keV/μm corresponds to the stopping power values of 250 MeV protons in water, this peak occurs due to the primary 250 MeV protons.The second peak is at y = 1.7 keV/μm that corresponds to the stopping power values of 35 MeV protons which is possibly produced from inelastic scattering and non-elastic nuclear reaction with water.Similar analysis for 100 MeV SOBP protons shows that the two peaks correspond to the 87 and 27.5 MeV protons.Moving along the depth-dose profile toward the distal end, the protons slow down and the distributions shift correspondingly to higher y values, as it can be seen in figure 5.The microdosimetric distribution varies significantly even within the SOBP region.At the centre of the SOBP, for 50 and 250 MeV SOBP proton beams, the majority of absorbed dose is delivered by events of size about 3.8 and 2.4 keV/μm, respectively, which correspond to stopping power of 16 and 12.5 MeV protons in water, respectively.Downstream this depth (towards the distal end), the maxima of the distributions move further to higher y values, corresponding to slowed-down protons of increased LET.For a given SOBP beam, the peak height of the distribution is maximum at the distal end of the SOBP region (position 'T' in figure 2).In the distal end, the proton  edge, corresponding to the maximum energy that protons can release in the simulated water-equivalent site of 1 μm diameter, is clearly visible at about y = 75 keV/μm which is due to low-energy protons of about 50 keV of energy.
All the microdosimetric distributions as shown in figure 5 show a long tail portion extended up to y = 1000 keV/μm which is due to the secondary particles produced in non-elastic nuclear reactions by high energy protons.These events have a low probability of occurrence.However, their contribution to the dose is relatively higher because of their high LET.Note that the majority of events at the tail portion is from the secondary alpha particles (see figure 6).
Table 4 presents ̅ y F and ̅ y D values for SOBP proton beams as a function of depth in water.For a given SOBP beam, ̅ y F value gradually increases with depth.However, variation in ̅ y D value with depth depends on the energy of the SOBP beam.This variation is significant at smaller energy (50 MeV SOBP protons) and the variation reduces with increasing energy.For example, in the case of 50 MeV SOBP beam, ̅ y D value at distal end of SOBP region (position 'T' in figure 2) increases by a factor of 2.5 as compared to that at the entrance region.Whereas, variation of ̅ y D value between entrance region and distal end of SOBP region is about 31% and 6%, respectively for 100 and 250 MeV SOBP beams.The increase in ̅ y D value at the distal end of the SOBP region as compared to that at the centre of SOBP (positions 'S' and 'T' in figure 2) are 90%, 8%, 14%, 11% and 22% respectively, for 50, 100, 150, 200 and 250 MeV SOBP proton beams.

RBE of SOBP protons
Figures 7(a)-(e) present on-axis RBE values calculated based on different models for 50-250 MeV SOBP beams as a function of depth in water, respectively.Each figure presents (i) RBE MKM for D p = 2 and 4 Gy, (ii) RBE r , (iii) RBE (4 Gy) according to the models described in sections 2.7.2 and 2.7.3 and (iv) RBE w (2 Gy).Comparison of RBE(4 Gy) and RBE(2 Gy) demonstrates the influence of dose/fraction on the RBE of SOBP proton beams.The figures show that the microdosimetry-based RBE values i.e.RBE MKM and RBE r : (a) are insensitive to depth in the entrance region and (b) gradually increase with depth within the SOBP region.For a given beam, the value of RBE MKM (4 Gy) at a given depth is comparable to the corresponding values of RBE r with a maximum difference of about 6% at the distal end of the SOBP region.Values of RBE MKM (2 Gy) and RBE MKM (4 Gy) are comparable for 100-250 MeV SOBP beams (difference is 3%-4%).Whereas RBE MKM (2 Gy) is higher than RBE MKM (4 Gy) by about 10% at the distal end of 50 MeV SOBP beam.The variation of RBE MKM or RBE r in the SOBP region is sensitive to beam energy and D p .Table 5 presents percentage increase in RBE MKM (4 Gy), RBE MKM (2 Gy) and RBE r at distal end as compared to that at the proximal end of SOBP region (increase at position 'T' as compared to 'R' in figure 2  2.30 ± 0.04 8.84 ± 0.17 6.09 ± 0.15 using any of the models) at the distal end of the SOBP region (position T in figure 2) are higher by: (a) about 20% for 50 and 100 MeV beams and (b) about 11% for 150-250 MeV beams.Similar analysis for RBE RMF shows an increase of about 15% for 50 and 100 MeV beams, and about 8% for 150-250 MeV SOBP beams.
The values of RBE(4 Gy) is sensitive to the calculation models (see figure 7).Despite this sensitivity, all models exhibit a similar trend in the variation of RBE with depths in water.Interestingly, the RBE w (4 Gy) value for a given depth and beam energy is comparable with the average RBE(4 Gy) value obtained from the corresponding RBE(4 Gy) values (same depth and beam energy) using the LET-based and RMF models.
To comprehensively understand the average trend of RBE(4 Gy) values as a function of depths in water, table 6 presents RBE w (4 Gy) values at different regions of the depth-dose profile for the investigated SOBP beams.The table highlights that RBE w (4 Gy) sig- nificantly deviates from a constant value of 1.1, particularly in the distal half of the SOBP region and the distal dose fall-off region.
The trend of variation of RBE w (2 Gy) with depth is similar to that of RBE w (4 Gy) (see figure 7).For a given beam RBE w (2 Gy) is comparable to RBE w (4 Gy) at the entrance region (within 2%-4%).In the SOBP region, the variation is sensitive to depth.Variation is less in the proximal end and more in the distal end of the SOBP region.For example, when depth varies from the proximal to the distal end of the SOBP region, the ratio of RBE w (2 Gy) to RBE w (4 Gy) is in the range of 1.1-1.18for 50 MeV beam and about 1.06-1.15for 100-250 MeV beams.It is important to highlight that figure 7 presents the models RBE MKM (2 Gy) and RBE w (2 Gy), representing microdosimetry-based and LET-based approaches, respectively.Notably, another microdosimetric modelbased RBE i.e.RBE r , is independent of dose/fraction.Among the various LET-based models, RBE w is chosen because, similar to RBE w (4 Gy), RBE w (2 Gy) at a given depth and beam energy is comparable with the average RBE(2 Gy) value derived from corresponding RBE (2 Gy) values (at the same depth and beam energy) using all the LET-based and RMF models.

Microdosimetry-based RBE versus RMF and LET-based RBE
For all the SOBP beams (except for 50 MeV at d = 2 cm), at a given depth RBE r is comparable to the corresponding RBE RMF (see figure 7).RBE r shows a maximum deviation with RBE Ch (4 Gy) in the SOBP region (about 13% for 50-200 MeV beams and 19% for 250 MeV beam).Whereas in the entrance region, the deviation is 7%-8%.RBE w (4 Gy) compares reasonably well with the corresponding RBE MKM (4 Gy) and RBE r .The maximum deviation between RBE MKM (4 Gy) and RBE w (4 Gy) is about 8% in the SOBP region of 250 MeV beam and for other SOBP beams, deviation is less than 4%.Whereas the maximum deviation between RBE MKM (2 Gy) and RBE w (2 Gy) in the SOBP region is about (a) 15% for 250 MeV, (b) 10% for 200 MeV and (c) 7% for 150-50 MeV beams.

RBE of SOBP and monoenergetic proton beams
Table 7 presents a comparison of TOPAS-calculated RBE values for 50 and 250 MeV monoenergetic protons with those of SOBP protons having the same nominal energy.Due to the difference in shapes of the depth-dose profiles for SOBP and monoenergetic protons, it becomes challenging to designate the same depth for both types of beams for comparison of RBE values.Therefore, the comparison is conducted in three specific regions: (a) the entrance region, (b) the peak position of monoenergetic beams versus the center of the SOBP region for SOBP beams, and (c) the distal end of the depth-dose profile.The findings in the table indicate that the RBE values for monoenergetic proton beams are comparable to corresponding RBE values for SOBP beams with the same nominal energy in the entrance region.When comparing RBE values at the peak position of 50 MeV monoenergetic beams with those at the center of the SOBP region for 50 MeV SOBP beams, the maximum difference (about 10%) is found for RBE MKM (4 Gy).A similar analysis for 250 MeV monoenergetic and 250 MeV SOBP beams reveals a maximum difference of about 8% in the case of RBE w (4 Gy).At the distal end, RBE values for monoenergetic protons, excluding RBE MKM (4 Gy) for 50 MeV, significantly differ from corresponding RBE values for SOBP proton beams.

Biological dose profile
Figures 8(a)-(e) compare biological dose calculated using depth-specific RBE w (2 Gy) and RBE w (4 Gy) values for 50-250 MeV SOBP proton beams.The figure also includes biological dose based on a constant RBE value of 1.1 as used in current clinical practice.Note that depth-specific RBE-weighted biological dose is defined as BD w (D p ) = relative dose (%) × RBE w (D p ) whereas constant RBE of 1.1 weighted biological dose is defined as BD 1.1 = relative dose (%) × 1.1.BD w (2 Gy), BD w (4 Gy) and BD 1.1 are comparable at the entrance region whereas they differ significantly in the SOBP region.BD w (2 Gy), BD w (4 Gy) exhibit a peak close to the distal dose fall-off region of the SOBP.In the SOBP region, BD w (2 Gy) is systematically higher than BD w (4 Gy).The present study calculated depth-specific RBE values using various LET-based models.However, variation of BD(D p ) with depths is shown using depthspecific RBE w values.As discussed in section 3.4, RBE w at a given depth and beam energy is comparable to the corresponding average RBE values estimated from all the LET-based and RMF models.Consequently, variation of BD w (D p ) with depths provides the average trend of the variation of BD(D p ) with depth.Although there is a reasonably good agreement between the two sets of values a minor deviation is observed which may be attributed to the methodology used in the radiobiological measurements which involved irradiating V79 cell lines under the beam line of a proton therapy facility.In the experimental studies (referenced in table 5), the initial kinetic energy of the primary proton beam was subjected to degradation while traversing various components of the beam line.Consequently, the resulting SOBP was composed of effectively lower energy protons compared to the initial kinetic energy.In contrast, TOPAS simulations did not account for this particular form of energy degradation during beam line traversal.This may contribute to the slight variation between the measured and calculated RBE values.
Present study shows that for a given proton beam, RBE w (2 Gy) are higher than the corresponding RBE w (4 Gy).As compared to the proximal end, the deviation of RBE w (2 Gy) from RBE w (4 Gy) is more at the distal end of the SOBP.This trend is consistent with several in vitro studies which showed an increase in RBE at low doses for late responding tissue Based on the recommendations of a Joint IAEA/ ICRU Working Group, at present a generic RBE of 1.1 is employed for proton therapy regardless of the beam energy, dose, depth, tissue type, overall treatment time, prior treatment and other factors (Fowler 1989, Goodhead et al 1992, Wedenberg et al 2013, Wouters et al 2015).This RBE value was derived from laboratory experiments that provided a mean RBE for animal tissue extrapolated to humans (Wouters et al 2015).The use of this generic RBE value of 1.1 produces uniform biological effect in the SOBP region as shown in figure 8.However, when depth-specific RBE value for a given beam is used, a uniform physical dose (absorbed dose) in the SOBP region does not deliver a uniform biological effect (see figure 8) and the biological dose shows sharp increase at the distal end of the SOBP (Fowler 1989, Goodhead et al 1992, Wedenberg et al 2013).Thus, the variation in biological dose across different SOBP energies and regions underline the importance of considering RBE variation in treatment planning, especially for tumours situated within or near critical structures.
In the present study, the calculated RBE values are for V79 cell lines and are not necessarily applicable to clinically relevant tissues (Polster et al 2015).However, as the / a b value for V79 cells is comparable to some late-responding tissues such as spinal cord, optic nerve, kidney, lung and liver, and hence it is possible that V79 cell lines may be predictive for the biological response of such tissues (Wouters et al 2015).If the cell-or tissue-specific parameters (a b) and ) are derived for reference radiation based on radio-biological experiments, then these models as described in section 2.7 can provide a good estimation of RBE of proton beams for any organ or tissue.

Conclusions
The present study investigated biological effectiveness of SOBP proton beams as a function of depth in water based on different models such as microdosimetrybased models, LET based models and RMF model using TOPAS Monte Carlo codes.The study also investigated the influence of dose per fraction on the calculated RBE.Microdosimetric distributions calculated at 1 μm site size were utilized for estimating RBE MKM and RBE r .The study shows that the shapes of the microdosimetric distributions including peak-height and peak position vary significantly depending on the energy of the SOBP beam and depth in water and this difference is reflected in the calculated values of RBE MKM and RBE r .For a given beam, while the absolute values of RBE may differ based on the model, the pattern of how RBE changes as a function of depth exhibits a similar trend across all these models.This suggests that the models are collectively capturing the fundamental behaviour of RBE with varying depths, even if they may diverge in their specific numerical predictions.
The study reveals that the RBE of a given SOBP beam exhibits minimal sensitivity to depth at the entrance region, suggesting that RBE ≈ 1.1 serves as a suitable approximation in this region.However, within the SOBP region, the RBE gradually increases and reaches its maximum value at the distal dose falloff region.It is noteworthy that in the SOBP region, the RBE exceeds 1.1, and this value varies with beam energy, depth, and dose per fraction.The study highlights the complex and dynamic nature of RBE across different regions of SOBP beams.
The calculated biological dose, which is obtained by weighting the physical dose with depth and dose-specific RBE values, is not uniform across the SOBP region and sharply increases towards the distal end of the SOBP.The study highlights that adopting a constant RBE value of 1.1 results in a systematic underestimation of the biological dose within the SOBP region.The extent of this underestimation depends on both the depth in water and the dose per fraction used.
In essence, this study emphasizes the need for precise consideration of RBE values in proton therapy treatment planning, accounting for variations in proton energy and different regions within the SOBP.It also highlights the potential significance of biological effects at the distal dose fall-off region, which could have implications for tumour control and healthy tissue sparing.Thus, in a clinical context, it is advisable to use variable RBE models over the simplistic constant RBE model of 1.1.However, the challenge emerges in selecting the precise RBE model to employ, particularly when various models yield different RBE values.This underscores the necessity for additional research to identify the most accurate RBE model for proton therapy treatments.
2018).Furthermore, radiobiological experiment for measuring RBE is expensive and very time consuming.Several phenomenological and biophysical models based on microdosimetric distributions and dose-average LET (LET d ) are proposed for predicting RBE more conveniently (Griffiths 1985, Paganetti et al 2002, Paganetti 2014, Giovannini et al 2016, Qutub et al 2016, Lühr et al 2018, Takada et al 2018, Conte et al 2019, Newpower et al 2019, Colautti et al 2020).The LET d is not a directly measurable quantity and it can be calculated using Monte Carlo methods as they are capable of transporting all the primary and secondary particles (Conte et al 2019).The precision and accuracy of simulated LET d values can be influenced by uncertainties in cross sections and the sensitivity of LET distributions to frequently adjustable simulation parameters (Conte et al 2019).Microdosimetric technique is an efficient approach to model RBE of ionizing radiations (Griffiths 1985, Anderson et al 2017, Conte et al 2019, Vassiliev et al 2019).Microdosimetry measures not only the average dose but the pattern of dose deposition at the micrometric scale (Conte et al 2019).Microdosimetric distributions can be measured and it can also be calculated using Monte Carlo codes.Several published studies are available on estimation of RBE of protons based on microdosimetric distributions and LET d(Takada et al 2018, Conte et al 2019, Newpower et al 2019, Colautti et al 2020).The majority of the microdosimetry-based studies are based on measured microdosimetric distributions and for limited SOBP proton beams such as 62 and 155 MeV SOBP beams (Kase et al 2013, Takada et al 2018, Conte et al 2019, Newpower et al 2019, Colautti et al 2020).It's worth highlighting that the instrument dependent threshold value for y in the range of 0.3-0.5 keV/μm for reducing the electronic noise may affect the bin-wise yield of measured microdosimetric distributions, particularly in the low y region (Rollet et al 2004, Chattaraj and Selvam 2020).Additionally, the measurements may suffer from inherent uncertainties such as detector position accuracy, experimental setup variations, variation in gas pressure of the filling gas of the Tissue Equivalent Proportional Counter (TEPC) etc. Notably, these limitations are absent in Monte Carlo calculations (Chattaraj et al 2019).Numerous published investigations have explored the assessment of RBE for clinical proton beams (Paganetti et al 2002, Chaudhary et al 2014, Paganetti 2014, Conte et al 2019).Predominantly, these studies are based on radiobiological experiments.A more limited subset of published works has focused on measured microdosimetric distributions specifically for 62 and 155 MeV SOBP proton beams (Kase et al 2013, Conte et al 2019).Yet, a significant gap remains in the literature, as there lacks a report on Monte Carlo-based microdosimetric distributions for SOBP proton beams aimed at RBE estimation as a function of proton energy and depths in water.

Figure 1 .
Figure 1.Schematic simulation geometry for calculation of phase-space at various depths.Green colour box represents 40 cm × 40 cm × 40 cm water phantom; Magenta colour boxes represents 2 cm × 2 cm × 10 μm water filled thin box for scoring phase-spaces.P, Q, R, S and T represent entrance region, dose-gradient region, proximal, central and distal end of the SOBP region, respectively.

2. 6 .
Calculation of RBE 2.6.1.Microdosimetry-based RBE calculation RBE can be calculated by convolving microdosimetric distributions with biological weighting function, r (y) (Pihet et al 1990, Coutrakon et al 1997, Paganetti et al 1997) or using MKM (Coutrakon et al 1997, Hawkins 2003, Kase et al 2006) or Theory of Dual Radiation Action (TDRA) (Kellerer and Rossi 1974).The r (y)-based approach is a biophysical model for predicting RBE (Coutrakon et al 1997, Paganetti et al 1997, Zhu et al 2019).The r(y)-based RBE is given by: was derived experimentally using neutron therapy beams with different radiation qualities and the acute effects in the intestinal crypt cells of mice was considered as the biological end point (Pihet et al 1990, Coutrakon et al 1997).Coutrakon et al (1997) demonstrated that RBE r shows reasonably good agreement with the experimentally obtained RBE at the 10% survival of V79 hamster cells in the absorbed dose range of 2-12 Gy of protons and 60 Co (reference radiation).The MKM is a phenomenological model introduced by Hawkins (Hawkins 2003) and then modified by Kase et al (2006) (Kase et al 2006, Zhu et al 2019).MKM-based RBE as a function of dose is calculated using the following formula:

Figure 2 .2
Figure 2. A typical SOBP proton beam representing various position of calculation depths.P, Q, R, S and T represent entrance region, dose-gradient region, proximal, central and distal end of the SOBP region, respectively.Portion TU represents distal dose fall-off region.

Model 1 :
Model by Wedenberg et al (2013): This model assumed β to be constant (β = β x ) and a a

d2
(1) gives the RBE W using the above conditions.Model byCarabe-Fernandez et al (2007) andCarabe et al (2012):This model used concept of maximum RBE (RBE max ) and minimum RBE (RBE min ).RBE max and RBE min correspond to the asymptotic values of RBE at Dose = 0 and ∞, respectively.The model explicitly assumed a dependency of β on LET d .As per this model, RBE as a function of D p is given by The constants of equations (4) and (5) is based on the experimental data of V79 Chinese hamster cells.Model 3: The model by Chen and Ahmad (2012) assumed a constant value of β (β = β x ) and a more complex dependency of α on LET d .The equation (i) is based on experimental data of V79 Chinese hamster cells.This model calculates RBE as a function D p according to the following equation: The model by Wilkens and Oelfke (2004): This model assumed β = β x and aThe relationship between a and LET d is based on the survival data of V79 Chinese hamster cells According to this model, RBE as a function of D p is given by All the cell line (or end point) dependent parameters required for the above-mentioned biophysical models (models 1-4) are defined in a parameter text file for V79 Hamster lung fibroblasts cells.The radiosensitivity parameters for reference radiation a = this parameter file were obtained from the experimental cell survival curve of V79 cells after 60 Co irradiation using non-linear least squares regression (Wouters et al 2015).
et al 2015).Note that q and k are independent of proton LET up to at least 75-100 keV/μm, which covers entire clinical energy range.RBE at iso-survival of reference radiation, according to this model, is calculated using the following equation (Frese et al 2012): following equitation (Grassberger and Paganetti 2011, Polster et al 2015): water medium, the above-mentioned unit of LET d is numerically equivalent to keV/μm.Where, dE is the energy loss by an event, r dE dx is the mass-collision stopping power.All primary and secondary protons and secondary electrons are taken into account for calculating LET d .Note that in the LET d scorer, the energy of secondary electrons is deposited at the position of their creation.

Figure 4 .
Figure 4. TOPAS-calculated depth versus LET d profile for SOBP proton beams with initial kinetic energy 50-250 MeV.

Figure 6 .
Figure 6.Contributions from different particles in the microdosimetric distributions of SOBP proton beams at 2 cm depth in water: (a) 250 MeV (b) 100 MeV.
).A comparison of the RBE(4 Gy) values calculated using various LET-based models and RMF model (see figures 7 (a)-(e) for SOBP beams of 50-250 MeV, respectively) shows that for a given beam energy and depth: (a) except in the distal dose fall-off region of 50 MeV SOBP beam, ( ) RBE 4 Gy Ch is always higher than the value produced by other models, and (b) as compared to the value produced by other LET-based models.In a nutshell, the comparison of the estimated RBE(4 Gy) values from all the LET-based RBE models indicate that ( values of RBE(4 Gy), respectively for a given SOBP beam at a specific depth, and the RBE(4 Gy) values resulted from other LET-based models lie within this two values.The difference between ( approaches towards the distal end of the SOBP region.For a given SOBP beam, variation in RBE (calculated using LET-based models and RMF model) with depth follows the similar trend as that of RBE MKM and RBE r .In the entrance region, depending on the calculation model and beam energy, LET-based RBE(4 Gy) value lies in the range of about 1.03-1.15( ( ) RBE 4 Gy Ch = 1.15 for 50 MeV SOBP protons) whereas ( ) RBE 4 Gy RMF is about 1.04, independent of beam energy.As compared to the proximal end (position R in figure 2), the LET-based RBE values (calculated

Figure 7 .
Figure 7. TOPAS-calculated RBE values at a prescribed dose of 4 Gy according to the models by Wedenberg et al, Wilkens and Oelfke, Carabe et al, Chen and Ahmad, RMF model and MKM for SOBP beams with initial kinetic energy: (a) 50 MeV, (b) 100 MeV, (c) 150 MeV, (d) 200 MeV and (e) 250 MeV.The figures also include RBE r values and RBE for a prescribed dose of 2 Gy according to the models by Wedenberg et al and MKM model.
RBE r = 1.14 RBE r = 1.49RBE r = 1.35 RBE r = 1.10 deviation for the majority of the measurement points (Polster et al 2015).With the parameters implemented for V79 cell line in TOPAS, the model by Wedenberg et al (2013) reproduces the experimental data best (Polster et al 2015).This comparison provide confidence to use these models for estimating RBE values at another clinical proton beam.Table 8 compares measured RBE values of V79 cell lines obtained from various radiobiological experiments (Robertson et al 1994, Paganetti 2014, Polster et al 2015, Wouters et al 2015) and the corresponding values calculated in the present study using TOPAS for different SOBP proton beams.
(Hall et al 1978, Blomquist et al 1993, Gueulette et al 1996, Wouters et al 1996, Bettega et al 2000, Wouters et al 2015).Note that the representative values of / a b for late and early-responding tissue are ~3 and 10, respectively (Goodhead et al 1992, Wedenberg et al 2013).The / a b value of V79 cells for 60 Co gamma rays is 1.44 ± 0.06 Gy.Thus, it is expected that V79 cells will exhibit increase in RBE at lower doses.RBE also shows dose dependency for high LET radiations due to the increased contribution of single-hit inactivation by densely ionizing particles (Wouters et al 2015).Note that high energy proton tracks have lower ionization densities whereas at the distal end of SOBP, the protons are slowed down (high LET) and hence ionization density is high (Wouters et al 2015).

Figure 9 .
Figure 9. Influence of upper cut-off value of y on the ȳ , F ȳD and y * values at the entrance region of (a) 50 MeV SOBP beam and (b) 250 MeV SOBP beams.

Table 1 .
Initial kinetic energies and SOBP width of the SOBP proton beams considered in the present study.
Carabe et al (2012), Chen and Ahmad (2012), and Wilkens and Oelfke (2004) which depend on dose, LET and α/β ratio.Note that α and β are specific parameters of the linear-quadratic (LQ) model.The details of these models are described byPolster et al  (2015).

Table 5 .
Percentage increase in RBE MKM (4 Gy), RBE MKM (2 Gy) and RBE r at distal end as compared to the proximal end of SOBP region as a function of proton beam energy.

Table 6 .
TOPAS-calculated RBE w (4 Gy) at different region of the depth-dose profile for different SOBP proton beams.