{'fraction_non_alphanumeric': 0.03151505072343256, 'fraction_numerical': 0.016339597538666224, 'mean_word_length': 4.567824074074074, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Estimation of relative biological effectiveness for boron neutron capture therapy using the PHITS code coupled with a microdosimetric kinetic model

Abstract The absorbed doses deposited by boron neutron capture therapy (BNCT) can be categorized into four components: α and 7Li particles from the 10B(n, α)7Li reaction, 0.54-MeV protons from the 14N(n, p)14C reaction, the recoiled protons from the 1H(n, n) 1H reaction, and photons from the neutron beam and 1H(n, γ)2H reaction. For evaluating the irradiation effect in tumors and the surrounding normal tissues in BNCT, it is of great importance to estimate the relative biological effectiveness (RBE) for each dose component in the same framework. We have, therefore, established a new method for estimating the RBE of all BNCT dose components on the basis of the microdosimetric kinetic model. This method employs the probability density of lineal energy, y, in a subcellular structure as the index for expressing RBE, which can be calculated using the microdosimetric function implemented in the particle transport simulation code (PHITS). The accuracy of this method was tested by comparing the calculated RBE values with corresponding measured data in a water phantom irradiated with an epithermal neutron beam. The calculation technique developed in this study will be useful for biological dose estimation in treatment planning for BNCT.


INTRODUCTION
Boron neutron capture therapy (BNCT) is a binary treatment modality that was advocated by Locher in 1936 and that can selectively target tumor tissues [1]. The two particles (an α particle and a lithium ( 7 Li) particle) generated in the neutron capture reaction with boron) have high linear energy transfer (LET) and a short range (of the order of one cell diameter). In principle, BNCT allows the preferential destruction of 10 B-loaded tumor cells, sparing normal tissues without 10 B. The radiation field of BNCT produced in normal tissues and tumors consists of a mixture of several particles, such as neutrons, protons, photons, α particles and 7 Li particles. In the estimation of the physical dose, these particles are categorized into four major dose components [2]: α particles and 7 Li particles from the 10 B(n, α) 7 Li reaction, 0.54-MeV protons from the 14 N(n, p) 14 C reaction, the recoiled protons from the 1 H(n, n) 1 H reaction, and photons from the neutron beam and 1 H(n, γ) 2 H reaction; these dose components are called the boron, nitrogen, hydrogen and gamma doses, respectively, in this paper. The boron dose from the neutron capture reaction accounts for the majority of the biological effect in tumor tissues [3]. The other dose components deliver background doses (from a mixture of high-and low-LET radiation components) to both normal and tumor tissues [4].
In order to estimate the treatment effect, the physical dose needs to be converted to the biological dose. Therefore, each of the four major dose components is multiplied by their relative biological effectiveness (RBE), which is determined on the basis of certain endpoints of clinical or experimental irradiation conditions. The RBE of each dose component depends on the variation in the physical dose according to the dose dependence of the cell survival curves. In addition, the RBE of the recoiled proton component depends on the variation in the spectral characteristics of the neutron beam. Therefore, it is necessary to evaluate the biological dose by taking into account the dependence not only on the dose, but also on the neutron spectrum. These dependencies are not currently considered in treatment planning for BNCT. Hence, development of a biological dose estimation model that can take these dependencies into account is needed.
A number of studies based on models that use a microdosimetric cellular scale have been devoted to the evaluation of the biological effectiveness of BNCT. For example, Bond et al. and Gabel et al. proposed the use of a 'hit size effectiveness function' (HSEF) that gives a cell survival fraction as a function of the energy deposition in a cell nucleus [5,6]. The HSEF was implemented in the microdosimetry model developed by Vliet-Vroegindeweij et al., which can calculate the cell survival fraction of all the important high-LET radiation components for BNCT [7]. Santa Cruz et al. developed a model for estimating the biological dose for BNCT, taking into account the synergistic interactions between different radiations [8,9]. However, none of the existing models can express the RBE of each of the four major BNCT dose components in the same framework. We have, therefore, developed a calculation method for estimating the biological dose for each of the four major BNCT dose components. This method is based on the microdosimetric kinetic (MK) model [10], which was derived from the theory of dual radiation action [11]. The MK model can calculate the cell survival fractions of various charged particles from the probability density of specific energy in a conceptual subcellular (submicron) structure of the cell nucleus (called a 'domain'). The MK model is adapted to the treatment planning for heavy charged-particle therapy in the Heavy-ion Medical Accelerator in Chiba [12,13].
In this study, the RBEs of each of the four major dose components were estimated using the MK model. For this purpose, it was necessary to calculate the microscopic probability density of lineal energy, y, in a macroscopic matter. In our method, the probability density of y in the domain was calculated using the microdosimetric function [14] implemented in the macroscopic particle transport simulation code (PHITS) [15]. Detailed procedures based on the MK model for estimating the RBE of each of the four major BNCT dose components are presented herein, together with the results of their verification using a past biological experiment for V79 Chinese hamster cells [16,17].

RBE calculation by the MK model
Using the MK model, the cell survival fractions of various charged particles can be estimated from the probability densities of specific energies in the domain. The procedures based on the MK model for estimating the cell survival fraction and RBE for each of the four major BNCT dose components are provided in this section.
According to the linear-quadratic (LQ) model, the cell survival fraction can be expressed as follows: where S is the survival fraction for the absorbed dose, D. The parameters α and β are the linear coefficient and quadratic coefficient of LQ parameters, respectively. The RBE of various charged particles for the cell survival fraction, S, is calculated as follows: where D γ (S) and D(S) are the absorbed dose for a reference radiation and a radiation of interest, respectively, for the similar value of the survival fraction, S. Note that the γ-rays emitted from the 137 Cs source are regarded as the reference radiation for estimating the RBE in this study. The parameters α γ and β γ are the LQ parameters that can be experimentally obtained from the survival curve of cells irradiated by γ-rays. Thus, the RBEs of any radiation can be estimated by determining their LQ parameters, α and β. The β parameter for each cell can be mathematically derived as a constant in the MK model [10], and its numerical value can be obtained from the survival curve of the cells irradiated by a reference radiation. In the Kase-proposed MK model, the parameter α can be estimated from the obtained β as follows [18,19]: where α 0 is a constant that represents the initial slope of the survival fraction curve in the limit of LET = 0. The parameter z Ã 1D denotes the saturation-corrected dose-mean specific energy defined in ICRU Report 36 [20], which can be calculated as follows: where l, m, ρ and r d are the mean chord length, mass, density and radius of the domain, and d(y, 2r d ) is the absorbed-dose probability density of y for site diameter 2r d . The parameter y* is the saturation-corrected dose-mean lineal energy, and y 0 is a so-called saturation parameter that indicates the lineal energy above which correction for the saturation due to the overkill effect becomes very important. If one can determine ρ, r d , α 0 , y 0 and β, the linear coefficient of α for any radiation field can be estimated by calculating d(y, 2r d ) in the field. Thus, preceding the calculation of the survival fraction using the MK model, these five parameters should be determined. The domain density ρ can be assumed to be 1.0 g/cm 3 . The parameter β can be obtained from the survival curve irradiated by a source of reference radiation, such as γ-rays, as mentioned before. The other three parameters r d , α 0 and y 0 can be determined by substituting the experimentally obtained α parameters for various radiation fields and their d(y, 2r d ) values into Eqs. (3) and (4), respectively.
In this study, the parameter β and the other three parameters r d , α 0 and y 0 for V79 Chinese hamster cells were determined from a past experiment that was performed to measure their cell survival fractions with γ-rays and a neutron beam, respectively. The d(y, 2r d ) was calculated using the microdosimetric function in the PHITS code [14]. The details are described below.

Estimation of cell survival fraction for the free-air experiment
To determine the MK model parameters, we estimated the α parameter for the boron, proton and gamma doses by reevaluating a past experiment measuring the cell survival fraction of V79 Chinese hamster cells incubated with boric acid of 0, 5 and 10 ppm 10 B concentrations irradiated with γ-rays from a 137 Cs source and an epithermal neutron beam in the Japan Research Reactor No. 4 (JRR-4) [17]. Boric acid serves as a reference boron carrier for the boron dose component because it is assumed that this compound distributes boron atoms uniformly throughout the cell, nucleus and external compartment [21,22]. The biological response of the boron dose component cannot be directly determined from the measured data. To estimate the survival fraction of the boron dose component, we compared the irradiation responses of boronated cells with those of non-boronated cells. The biological responses of protons from nitrogen capture reactions and hydrogen recoils cannot be experimentally distinguished, because they are attributed to protons emitted from the natural composition of the cell nucleus (i.e. unlike the attribution from boron, which is artificially added to the cells). This dose component is termed the 'proton dose' in this report. The cell survival fraction of the proton dose component was estimated by comparing cell survival fractions for non-boronated cells irradiated with neutron beams with those irradiated with γ-rays.
The β term in the LQ model can be ignored for the limit of D = 0, when the parameter α is higher than β. Then, the survival fraction can be simply approximated as follows: where s is the survival fraction for only the very low-dose irradiation. Figure 1 shows the survival curves of V79 Chinese hamster cells directly irradiated by the JRR-4 neutron beam; this direct irradiation is called a 'free-air experiment', in contrast to a 'phantom experiment', in which cells are irradiated in a water phantom, as described later.
The survival fraction of the boron dose component is expressed as the ratio of the survival fractions with and without boron. For very low-dose irradiation, the survival fraction of the boron dose component, s B , can be written as follows: where s wtB and s woB are the survival fractions of boronated and non-boronated cells, respectively. α wtB and α woB are α parameters for boronated and non-boronated cells, respectively, and α B is the α parameter for the boron dose only. The absorbed dose for boronated cells, D wtB , is expressed as the sum of the boron dose component, D B , and the dose for non-boronated cells, D woB ; i.e.: From Eqs. (6) and (7), α B is expressed as follows: Therefore, the α parameter of the boron dose component, α B , is calculated using the value of D woB /D B , and α wtB and α woB obtained from the survival curves of V79 Chinese hamster cells. The dose components D B and D woB are expressed as follows: where K B,th is the neutron kerma factor for the boron dose component, for which the numerical value is evaluated to be 7.51 Â 10 −14 Gy·cm 2 /ppm [23], φ th is the thermal neutron fluence measured by gold foils, and C B is the 10 B concentration in the conditioned medium in parts per million. D N , D H and D γ are the nitrogen, hydrogen and gamma dose components, respectively. D woB is expressed as the sum of the nonboron dose components. Each dose component D N and D H is estimated as follows: where K N,th is the neutron kerma factor for the nitrogen dose component, for which the numerical value was evaluated to be 1.50 Â 10 −13 Gy·cm 2 [23]. K H,th is the special kerma factor that converts thermal neutron flux to the hydrogen dose; this value depends on the neutron spectrum, and its numerical value for the JRR-4 epithermal neutron beam was evaluated to be 4.03 Â 10 −13 Gy·cm 2 [23]. The gamma dose was measured using a thermoluminescent dosimeter (TLD) that had beryllium oxide as a thermoluminescence material and was covered with quartz glass (Panasonic 170LS). According to Eq. (1), the survival fraction of the boron dose component for a wide dose range, S B , can be expressed as follows: where α B can be determined by Eq. (8), and β can be obtained from the survival curves of γ-ray irradiations. Similar to Eq. (8), the α parameter of the proton dose, α p , for low-dose irradiation is expressed as follows: where D p is the absorbed dose for the proton dose expressed as the sum of the nitrogen dose component and the hydrogen dose component.

Determination of MK model parameters
Probability densities of the lineal energy in microscopic sites cannot be directly calculated from the conventional particle transport simulation of the Particle and Heavy Ion Transport code System (PHITS), since it is designed to simulate particle motions on a macroscopic scale and employs a continuous-slowing-down approximation for calculating the energy loss of charged particles. Therefore, a special tally function was introduced into the code for calculating the microdosimetric probability densities (using a mathematical function that can instantaneously calculate quantities around trajectories of charged particles). The function was developed on the basis of track structure simulation, taking into account the productions of δ-rays and Auger electrons [14,24]. Using this function, we can calculate the probability densities of y for all types of charged particles and target sizes with a precision equivalent to that of the microscopic track-structure simulation. In our calculation, the probability densities of the deposited energy in water target sites for the boron, nitrogen, hydrogen and gamma doses were calculated separately using the sources defined as follows: as protons for which the energy spectra were the same as those of protons generated in a test tube filled with water and a tissue-equivalent material [25] with a diameter of 4 cm, respectively, irradiated by the epithermal neutron beam in the JRR-4. Here, tissue composition was assumed to be: H (10.1%), C (11.1%), N (2.6%) and O (76.2%) by weight percentage. (iv) The source of the gamma dose was defined as monoenergetic photons with an initial energy of 0.662 MeV.
The parameters necessary in the MK model for estimating the cell survival fraction are β, r d , α 0 and y 0 . The parameter β can be obtained from the survival curve of γ-ray irradiations. The other three parameters, r d , α 0 and y 0 , can be estimated from the least-square fittings of all the experimental data for the free-air condition. These fittings were performed: (i) to calculate the d(y, 2r d ) for the boron, proton and gamma doses from the PHITS code using the temporary r d value; (ii) to estimate the α 0 and y 0 from the least-square fittings of all the experimental data based on Eqs. (3) and (4); (iii) to update the r d value as the minimum condition of the error for the α parameter between the measured and calculated values; and (iv) to go back to procedure (1) until the values of the three parameters were converged.

RBE estimation in a water phantom
To validate the applicability of our method for estimating the biological dose of BNCT, we estimated the RBE of the boronated and non-boronated cells in a water phantom irradiated by the epithermal neutron beam in the JRR-4. A schematic setup for the phantom experiment is illustrated in Fig. 2. The RBE can be estimated by substituting α values calculated using the MK model parameters (obtained from the free-air experiment) into Eq. (2). The estimation of y* at the cell locations in the phantom was calculated by performing the PHITS simulation using the microdosimetric function. The accuracy of this method was tested by comparing the calculated survival fraction and RBE with the corresponding experimental data obtained by Yamamoto et al. [16].

α value for the free-air experiment
The survival curves for the absorbed dose under free-air conditions both with and without boron are also shown in Fig. 1. The LQ parameters for the γ-ray radiations were determined to be 0.115 ± 0.004 Gy −1 and 0.0220 ± 0.0004 Gy −2 for α γ and β, respectively. The numerical value of α woB was estimated to be 0.626 ± 0.018 Gy −1 , and that of α wtB was estimated to be 0.779 ± 0.024 Gy −1 and 0.895 ± 0.054 Gy −1 for boron concentrations of 5 ppm and 10 ppm, respectively.
The value of D woB /D B obtained was 2.86 and 1.43 for 5 ppm and 10 ppm, respectively. By substituting these parameters into Eq. (8), α B was calculated to be 1.214 ± 0.105 Gy −1 and 1.281 ± 0.133 Gy −1 for 5 ppm and 10 ppm, respectively. However, the α B should be independent of the 10 B concentration, and thus we adopted the average value, 1.247 Gy −1 , in the following analysis. The parameters α woB and α γ were estimated to be 0.626 ± 0.018 Gy −1 and 0.115 ± 0.004 Gy −1 , respectively. The value of D γ /D p obtained was 0.94. By substituting these parameters into Eq. (14), α p was calculated to be 1.106 ± 0.035 Gy −1 . Figure 3 illustrates the comparison of the cell survival fraction for the boron and proton dose components estimated from the α value obtained by re-evaluation of the past experimental results using the method described above with our simulation using LQ parameters obtained from the MK model. Figure 4 shows the calculated result of the absorbed dose probability densities, d(y, 2r d ), for the boron, nitrogen, hydrogen and proton doses. The numerical value of r d was set to 0.222 μm. It was found from the graph that the direct and δ-ray contributions could be clearly observed in the higher and lower y values, respectively, for the boron dose. The d(y, 2r d ) spectrum for the hydrogen dose shifted to lower y events than that for the nitrogen dose because of the higher energy deposition by recoil protons.

Determination of MK model parameters
The parameters necessary in the MK model for estimating the cell survival fraction are β, r d , α 0 and y 0 . The evaluated values for β, r d , α 0 and y 0 were 0.022 Gy −1 , 0.222 μm, 0.052 Gy −1 and 120 keV/ μm, respectively, in this study.    Figure 6a illustrates the calculated and measured [26] thermal neutron fluxes for the incidence of the epithermal neutron beam to the phantom. The statistical uncertainties of the calculated data were <1%. The calculated data were normalized to reproduce the maximum value of the measured thermal neutron flux. These distributions agreed with one another within 5%, indicating the accuracy of dose calculation by the PHITS code. Figure 6b illustrates the calculated absorbed dose for each dose component. The boron doses were estimated by setting the boron concentration at 5 ppm. The calculated data were normalized by the same scaling factor used in panel (a). The gamma dose made the dominant contribution to the total dose in this calculated condition, because the boron concentration (5 ppm) was lower than that used in the practical situation of BNCT.

Cell survival fractions for the free-air experiment
The survival curves of boron and proton doses estimated from the MK model using the obtained parameters are shown in Fig. 3. It is evident from the graphs that the calculated data for both the boron and proton doses agree fairly well with the corresponding measured data. These agreements indicate that the cell survival fraction for the free-air experiment of each dose component can be simply estimated from y* calculated using the MK model.

Boron carrier
Boric acid served as a reference boron carrier for the boron dose component because it is assumed that this compound distributes boron atoms uniformly throughout the cell, nucleus and external compartment [21,22]. In this study, the RBE of the boron dose component is estimated on the basis of this assumption.
The biological effects strongly depend on the localization of the boron compound because of the short range of the α and Li particles generated in the neutron capture reaction in boron [6]. Therefore, more accurate characterization of the microscopic distribution of the boron compounds and the ability to consider heterogeneous distributions are required. In our model, the approach would be applicable only for compounds that are homogenously distributed throughout the irradiated volume. Therefore, further studies are required to estimate the biological effect, taking into account the microscopic distribution of the boron compounds (e.g. p-boronophenylalanine or sodium borocaptate).

RBE estimation in a water phantom
As can be seen in Figs, 7 and 8, survival fractions and RBEs for various depths in the phantom can be simply estimated from the PHITS simulation coupled with the MK model, whose parameters are determined from the free-air experiment. It is evident from the graphs that the calculated data for both the boronated and non-boronated cells agree fairly well   Fig. 7). Panel (a) represents the data for cells with 5 ppm 10 B concentration, whereas Panel (b) represents the data for the non-boronated cells.
with the corresponding measured data. However, the discrepancies found at positions deeper than 40 mm in Fig. 7a and at both less than 15 mm and deeper than 40 mm in Fig. 8a are probably attributable to experimental uncertainty and the difference between measured and calculated conditions, such as the neutron spectrum in the phantom. The depth dependences of these data are predominantly attributable to the depth-dose characteristics for the each dose component in the water phantom. In addition, the cell survival fraction and RBE for the hydrogen dose component depends not only on the dose but also on the spectral characteristics of recoiled protons (attributable to variation in the neutron spectrum with depth in the phantom).
We have established a new method based on the MK model for estimating the cell survival fraction or RBE for all BNCT dose components simply from y* calculated using the PHITS code. The RBE values can be determined by taking into account the dependence not only on the dose but also on the neutron spectrum. Using this method, the biological dose can be estimated from the RBE multiplied by the physical dose, which is also calculated using the PHITS code. An advantage of applying this method to treatment planning for BNCT is that it can calculate the required physical dose and RBE for the BNCT dose components in the same framework. Before using the model practically in treatment planning for BNCT, further studies are required to estimate the biological effect while taking into account the localization of the boron compound in the target cells [6,8,27]. It is also necessary to consider the stochastic nature of absorbed doses, not only for the domain but also for the cell nucleus, which is an important concept in high-LET and hypofractionated radiation therapy [28].