1 Introduction

In 2018, the first experimental evidence of proton boron capture therapy was published [2]. Here, a biological effect was observed when DU-145 prostate cancer cells had been exposed to sodium borocaptate \({\hbox {Na}_{2}\hbox {B}_{12}\hbox {H}_{11}\hbox {SH}}\) (BSH) during irradiation. The BSH used in this experiment was consisted of natural boron, i.e., roughly 19 % \({{}^{11}\hbox {B}}\) and 81 % \({{}^{10}\hbox {B}}\). The observed effect was a decrease in the cell survival fraction of the BSH-exposed DU-145 cells. The effect was confined to the spread out Bragg peak and not observed in the entry region of the proton beam. It was originally attributed to the production of additional high-LET alpha particles from the \({\hbox {p} + {}^{11}\hbox {B} \longrightarrow {}^{4}\hbox {He} + {}^{4}\hbox {He} + {}^{4}\hbox {He}}\) (boron proton capture, BPC) reaction. As has been mentioned in the original publication as well as others (e.g., [4]), the yield from this reaction in a proton beam is likely very small, since the maximum total cross section is \(\approx 1\) barn at \(600 \hbox { keV}\) [5]. The question of how many \({{}^{4}\hbox {He}}\)-particles are needed to see a biological effect is still unanswered. Independently, another boron reaction has been proposed as a source of high-LET particles in proton therapy, namely, the \({\hbox {n}_{\textrm{th}} + {}^{10}\hbox {B} \longrightarrow {}^{7}\hbox {Li} + {}^{4}\hbox {He}}\) (Boron neutron capture, BNC) reaction, driven by the thermal secondary neutrons, ubiquitous in targets exposed to clinical proton beams [1]. The previous work comparing the BPC and BNC reaction mechanisms in boron-enriched proton therapy is found in ([6, 7]). They both conclude the dose contribution from the BPC reaction to be higher than the BNC reaction. These studies, however, only investigated small pristine proton beams on small phantoms or partially incomplete physics packages. In this study, we will discuss the use of physics packages and also quantify the neutron fluence dependency on field size, energy, and phantom size. When simulated with clinically more realistic parameters, we find a factor of \(4000 \pm 700\) high-LET particle production from NCEPT than PBCT.

2 Methodology

We used the Monte Carlo particle transport code TOPAS 3.7 [3] for our simulations, which were additionally replicated in SHIELD-HIT12A [8] for cross-checking. The field shape used in these simulations is box-shaped fields. The physics modules used in the TOPAS simulations are as follows:

  • g4em-extra

  • g4em-standardopt4

  • g4h-elasticHP

  • g4h-QGSPBICAllHP

  • g4decay

  • g4stopping

We here highlight the importance of g4h-elasticHP. The thermalization of neutrons occurs through elastic scattering; therefore, if this is not included, only high-energy neutrons from the g4h-QGSPBICAllHP module are generated. By scoring the reactions directly in TOPAS, we are able to examine how the g4h-QGSPBICAllHP package handles the BPC and BNC reactions. This is done by scoring the projectile particle, target nucleus, and produced secondaries. We simulate a thin target of pure \({{}^{10}\hbox {B}}/{{}^{11}\hbox {B}}\) to acquire the cross section directly from TOPAS. The acquired cross sections are shown in Figs. 1 and 2 together with experimental cross section from the EXFOR database. It should be noted that Fig. 2 includes two decay channels to the finale triple \(\alpha\) state, which can be seen as the two groupings of experimental data. Namely, \({\hbox {p} + {}^{11}\hbox {B} \longrightarrow 3\,\, {}^{4}\hbox {He}}\) and \({\hbox {p} + {}^{11}\hbox {B} \longrightarrow {}^{12}\hbox {C} \longrightarrow {}^{8}\hbox {Be} + {}^{4}\hbox {He} \longrightarrow 3\,\,{}^{4}\hbox {He}}\), where the prior is the more dominant. The geometry of our setup consists of a 20 \(\times\) 20 \(\times\) 20 \({\hbox {cm}^{3}}\) phantom, with a cubic-shaped field of 3 \(\times\) 3 \({\hbox {cm}^{2}}\) wide with the SOBP positioned at a depth of 7.7–13.3 cm in the beam direction. Our scoring volume has a width of 1 cm \(\times\) 1 cm in the lateral directions, i.e., lateral equilibrium is retained. Scoring along depth is in 300 steps along the entire phantom (20 cm).

After determining the relevant physics modules need for the simulations, we made a parameter study where we varied the phantom and field size while recording the fluence of low-energy protons and epithermal neutrons. This was done by scoring the averaged fluence inside the target volume covered by the SOBP, that is, to say the average fluence between 7.7 and 13.3 cm in Fig. 5. We did simulations with increasing phantom size of 1 cm in each lateral direction pr step. A similar setup can be made with increasing field sizes. Here, the phantom size is fixed at 100 cm \(\times\) 100 cm \(\times\) 100 cm, and the field is increased in the lateral directions 1 cm for each step. For increasing field/phantom size, we normalize the average dose in the SOBP to 1 Gy. We then investigate the fluence of relevant energies of proton and neutrons needed to drive the BPC and BNC reactions.

Next, we investigate the produced secondaries. We note that other work tend to limit the scope to \(\alpha\) particle production but we will here also include \({}^{7}\text {Li}\) from NBC, since its mean free path range is also on the order of a cell diameter (2–6 \({\upmu \hbox {m}}\)). Including a fragmentation model for the PBC reaction in TOPAS is not trivial. Experimentally obtained cross sections disagree with a factor of two in the relevant energy range \(< 2 \hbox { MeV}\) [5, 9, 10], and a correct quantum mechanical implementation of the BPC reaction with break-up amplitudes described in [11] and [12] is out of scope for this work. Instead, we implement a custom scorer in TOPAS for scoring the yield of a reaction provided a given cross-section table. Despite the experimental cross-section disparity of the BPC reaction, the energy range where the maximum cross section is located is consistent. From the EXFOR database, the largest cross section is below \(1 \hbox { MeV}\), peaking at \(\approx 600 \hbox { keV}\) [5]. For energies larger than \(1 \hbox { MeV}\), the cross section decreases with \(\approx 90\%\). We, therefore, chose to score the fluence of protons \(< 2 \hbox { MeV}\) to quantify the protons available for the PBC reaction. For neutrons, the BNC reaction still has a cross section > 10 barns at \(1 \hbox { keV}\) [13], so scoring fluence of neutrons \(< 1 \hbox { keV}\) is representative of the BPC reaction yield. The maximum cross section for BCN is well established, found at \(3839 \pm 9\) barns for \(0.253 \hbox { eV}\) [14]. Additionally, the yield of the two reactions was estimated. The proton and epithermal fluence were also obtained independently using SHIELD-HIT12A [8].

3 Results

We here present the findings of work described in the methods section. The results will be divided into three sections.

3.1 Geant4 physics package verification

When investigating the Geant4 physics package g4h-QGSPBICAllHP, we start by plotting the cross section for the BNC and BPC reactions. These are shown in Figs. 1 and 2. Figure 1 shows good agreement with the experimental results of the BNC reaction for protons \(< 1 \hbox { MeV}\). Above \(1 \hbox { MeV}\) was not simulated. For the BPC reaction, there is no cross section for protons below \(1 \hbox { MeV}\) available in the g4h-QGSPBICAllHP package (since it is based on TENDL version 1.3.2). Furthermore, the reaction products simulated with this package are inconsistent: From a given nuclear interaction between protons and \({}^{11}\text {B}\), particles different from alphas are frequently generated, violating baryon number conservation. This observation has only been seen in the g4h-QGSPBICAllHP physics package for the BPC reaction. Other reactions from the same and different physics packages have not shown this issue. This makes the g4h-QGSPBICAllHP physics package unfit for simulating BPC reactions (for TOPAS 3.7), since both the reaction yield and dose contribution are likely to be incorrect. We, therefore, provide a custom extension for TOPAS to calculate the yield of a given reaction provided a cross section as an input file. This extension does not generate the reaction products, and therefore, the particle dose contribution from this reaction is not established in this work. The yield can, however, provide a relative comparison between the BPC and BNC reactions.

Fig. 1
figure 1

Total cross section at different energies for boron neutron capture. The scored cross section from TOPAS with the g4h-QGSPBICAllHP physics package is added as a blue line. The experimental data are from the EXFOR database

Full size image
Fig. 2
figure 2

Total cross section at different energies for boron proton capture from the EXFOR database. The blue line is the obtained cross section from TOPAS with the g4h-QGSPBICAllHP physics package. The experimental data are from the EXFOR database

Full size image

3.2 Geometrical parameter study of low-energy neutron and proton fluence

We now turn to the geometrical effects of neutron thermalization in our setup. We plot the average fluence of low-energy neutrons and protons in the SOBP inside the proton beam vs lateral phantom size, which is shown in Fig. 3.

From Fig. 3, it is observed that increasing the phantom size in the lateral directions can increase the average fluence of low-energy neutrons in the SOBP with up to an order of magnitude. This can explain why articles like [6] do not find the BNC reaction to produce more high-LET particles than the BPC reaction, since there simulated phantom size was 4 \(\times\) 4 \(\times\) \(10 {\hbox { cm}^{3}}\). Which from figure 3 means that they are producing close to an order of magnitude less low-energy neutrons, than in a clinically more realistic phantom size. For this specific beam profile, it is shown that you need at least a 10-cm wide and tall phantom to reach lateral equilibrium for neutron thermalization. This effect should be understood in terms of the thermal neutron “soup” as described in [1]. Thermal neutrons are generated from scattering of higher ( \(> 10 \hbox { MeV}\)) energy neutrons. The neutrons are generated mainly from different inelastic scattering on oxygen directly or decay of \({{}^{5}\hbox {He}}\) secondaries (based on simulations). The thermal neutrons move in the phantom without a common direction, thereby filling the phantom as a kind of “soup”. Equilibrium is then reached when there is as many neutron moving out of the SOBP area as moving in. The equilibrium point is specific for each beam size and energy. Meaning you could also have the phantom size constant and increase the field size and have a similar result. This is shown in Fig. 4. In [7], they used a field diameter of 2 mm, which compared to Fig. 4 is not representative of the potential neutron fluence obtained with a more clinically relevant field size. This could explain why they find the BNC reaction to have a lower impact than the BPC reaction.

Fig. 3
figure 3

Fluence vs lateral phantom size. By increasing the phantom size in the lateral directions to the beam axis, we obtain the fluence of neutrons \(< 1 \hbox { keV}\) and protons \(2 < \hbox { MeV}\) in the SOBP. Here, we observe an increase in one order of magnitude for the neutron fluence, when the lateral phantom size is increased from 2.5 to 7.5 cm. The data are normalized to 1-Gy average in the SOBP

Full size image
Fig. 4
figure 4

Average SOBP fluence vs lateral field size. Similar plot as Fig. 3, but for a fixed large phantom size of 100 cm \(\times\) 100 cm \(\times\) 100 cm, and increasing lateral field size. The neutron fluence further increases by at least one order of magnitude. The data are normalized to 1-Gy average in the SOBP

Full size image

We conclude that increasing the field size can gain another order of magnitude increase in neutron fluence in the SOBP, for larger phantoms. The field and phantom size chosen for experiments/simulations can thereby change the produced low-energy neutrons with up to two orders of magnitude.

3.3 \({}^{11}\)B and \({}^{10}\)B reaction yield comparison

We investigate the low-energy neutron fluence inside the proton field only, where the BPC reaction competes with the BNC reaction. A convenient way to visualize this is to plot the fluence of the particles needed for the two reactions, that is, neutrons \(< 1 \hbox { keV}\) and protons \(< 2 \hbox { MeV}\). This can be plotted along the beam axis within the proton field, as shown in Fig. 5.

Here we have a close to two orders of magnitude difference between neutrons \(< 1 \hbox { keV}\) and protons \(< 2 \hbox { MeV}\). The yield for the two reactions can be calculated as follows:

$$\begin{aligned} Y_x=\Phi \sigma _x \rho _A \end{aligned}$$
(1)

Where \(Y_x\) is the yield for reaction x, \(\Phi\) is the fluence of the participating projectile, \(\sigma _x\) is the microscopic cross section for reaction x, and \(\rho\) is the density of atoms in the target. We can calculate the ratio between the two reaction yields as follows:

$$\begin{aligned} R = \frac{Y_{\textrm{BNC}}}{Y_{\textrm{BPC}}} = \frac{\Phi _n}{\Phi _p}\cdot \frac{\sigma _{\textrm{BNC}}}{\sigma _{\textrm{BPC}}}\cdot \frac{\rho _{\textrm{B10}}}{\rho _{\textrm{B11}}} \end{aligned}$$
(2)

We estimate the expected R value based on Fig. 5 and determine each fraction in (2) separately. \(\Phi _n/\Phi _p = 10^{2}\) on average in the SOBP. The \(\sigma _{\textrm{BNC}}/\sigma _{\textrm{BPC}}\) = 10, from choosing the lowest possible \(\sigma _{\textrm{BNC}}=10\) barn and the highest possible \(\sigma _{\textrm{BPC}}=1\) barn. \(\rho _{\textrm{B10}}/\rho _{\textrm{B11}} = 1/4\) since we work with natural boron (\(80\% {{}^{11}\hbox {B}}\) and \(20\% {{}^{10}\hbox {B}}\)). Since BPC produces three high-LET particles and BNC produces two, we multiply R with 2/3 to get the ratio of produced particles and not just ratio of number of reaction. Here, it should be noted that we ignore the increased range of the reactions products from BPC reaction compared to BNC. The result is \(R =17.5\), which means that we produce a factor of 17.5 more high-LET particles from the BNC reaction than from the BPC reaction. If we chose the highest \(\sigma _{\textrm{BNC}}\) and lowest \(\sigma _{\textrm{BNC}}\), we get \(R = 46.67 \cdot 10^{3}\). From the preceding discussion, using a 2-mm thick beam as in [7], the neutron and proton fluences are approximately \(14 \times 10^{4}/\text {cm}^{2}\) and \(6 \times 10^{4}/\text {cm}^{2}\), respectively (derived by fitting a second-order polynomial to Fig. 4). This signifies that compared to an 8-cm field, we see from Eq. 2 that the neutron contribution is underestimated by nearly two orders of magnitude.

Fig. 5
figure 5

Fluence of neutrons and protons for different energy ranges. The simulations are performed in TOPAS with the geometrical setup described in the text and simulation software described in the methods section

Full size image

This suggests that the main reason for the observed enhanced RBE in [2, 15] arrives from the BNC reaction. Note also that the thermal neutron fluence is reduced by one order of magnitude in 1–2 cm from the phantom surface, such that the dominating BNC reaction is suppressed near phantom surfaces (and therefore in particular small phantoms). While the BPC reaction undoubtedly will add to the high-LET particle yield, it is unlikely to be the main contributor by at least an order of magnitude.

To obtain a more precise estimation of the R value, we developed a custom code for calculating the yield of a specified reaction using TOPAS. For the PBC reaction, we utilized the cross-section data from Fig. 2 from the EXFOR database, while for the BNC reaction, we relied on the g4h-QGSPBICAllHP physics package that was verified in Fig. 2. Our boron target was a cell flask placed in the SOBP (at 9-cm depth, see Fig. 5) with dimensions of 92 \(\times\) 51 \(\times\) 29 mm, which corresponds to the size of an actual culture flask used at the Danish Centre of Particle Therapy. We selected a boron medium concentration of 17 mg/mL, which is 100 times higher than the concentration used in the original experiments [2, 15]. This concentration, irradiated with a 10 cm \(\times\) 10 cm field (SOBP positioned at a depth of 7.7–13.3 cm as shown in Fig. 5), produced approximately 1 PBC reaction for every 10.000.000 simulated protons. We simulated 200.000.000 protons, resulting in an R value of \(4000 \pm 700\), with the uncertainty calculated from counting statistics. On average, we generate 4.65 alpha particles per 10.000.000 protons for the PBC reaction, compared to 12,440.6 alphas and \({}^{7}\text {Li}\) per 10,000,000 protons for the BNC reaction. Using the highest possible experimental cross sections from Fig. 2 (approximately 1 barn) from 0.14 to 45 MeV yields an R value of \(22.2 \pm 1\).

4 Conclusion

We show that the fluence of low-energy neutrons (\(< 1 \hbox { keV}\)) strongly depends on the lateral phantom size as well as the field size, which was not taken into account in the previous studies of PBCT. Changing either may affect the low-energy neutron fluence with up to two orders of magnitude. The BNC reaction produced \(4000 \pm 700\) times more high-LET particles than the BPC reaction in the SOBP, when considering clinical dimensions for phantom and field size. Whether the BNC (and BPC) are biologically relevant is out of scope for this work, but remain to be investigated.