2.1 Principles for detecting three-dimensional boron concentration and dose in BNCT

In this study, a detection method for the three-dimensional distribution of boron concentration is proposed based on the reaction mechanism of neutrons captured by ¹H and ¹⁰B isotopes. The detection method can be expressed through Eq. (1), in which ${\omega }_{{\text{10}}_{\text{B}}}$ (i.e., the ¹⁰B concentration in a voxel of interest) can be calculated as follows:

where $M_{{\text{1}}_{\text{H}}}$ is the molar mass of the ¹H isotope (g/mol), $M_{{\text{10}}_{\text{B}}}$ is the molar mass of the ¹⁰B isotope (g/mol), $N_{0.478}$ is the yield of PG rays with an energy of 478 keV, $N_{2.224}$ is the yield of PG rays with an energy of 2.224 MeV, ${\omega }_{{\text{1}}_{\text{H}}}$ is the ¹H isotope concentration in the voxel of interest, ${\omega }_{{\text{10}}_{\text{B}}}$ is the ¹⁰B isotope concentration in the voxel of interest, and K is the ratio of the cross section of neutrons captured by ¹⁰B and ¹H isotopes. In combination with the average energy released in a boron neutron capture reaction (i.e., E_ave of 2.3388 MeV) and the mass of voxel of interest (m), the boron dose deposited in the voxel of interest ($D_{{\text{10}}_{\text{B}}}$) can be calculated using Eq. (2).

The detailed derivation of the proposed three-dimensional boron concentration detection method used in this study is as follows. As shown in Fig.1, when an epithermal/thermal neutron beam travels through the tissues of interest, neutrons interact with different isotopes in the tissues. Considering the content of different isotopes in the tissues and the neutron capture cross section, with BNCT, incident neutrons mainly interact with ¹H and ¹⁰B isotopes. These reactions emit PG rays with different energies (i.e., 478 keV and 2.224 MeV).

Fig. 1Full-screen View

(Color online) Schematic diagram of PG rays emitted from neutron capture reaction between ¹⁰B and ¹H isotopes

To measure the three-dimensional distribution of the boron concentration, we calculated the boron concentration at the voxel level. For a voxel of interest, the distribution of ¹⁰B and ¹H concentrations can be uniformly distributed. Thus, the production rate of different PG rays in the voxel of interest within the neutron irradiation can be described as follows:

where ${\sigma }_{{\text{1}}_{\text{H}}}(E)$ is the capture cross section of the ¹H isotope (barn), ${\sigma }_{{\text{10}}_{\text{B}}}(E)$ is the capture cross section of the ¹⁰B isotope (barn), $\varphi (E)$ is the fluence rate of neutrons incident into the voxel, E₀ is the minimum value of the neutron energy, E₁ is the maximum energy (MeV), and N_A is Avogadro’s constant.

For a given voxel, the concentrations of different elements in a voxel can be considered as constants. Thus, the ratio of PG ray production in a voxel can be described as

where K(E) is the ratio of the cross section of neutrons captured by the ¹⁰B and ¹H isotopes. The cross-section of the reaction between neutrons and elements is only related to the energy of the neutron; thus, the physical quantity K(E) can be considered to be related only to the neutron energy.

For a BNCT epithermal/thermal neutron beam, as shown in Fig. 2, the ratio of neutron cross sections that produce the PG rays with 478 keV and 2.224 MeV is almost a constant value over the epithermal/thermal neutron energy range (i.e., less than 10 keV). In the clinical treatment of BNCT, the neutron beam used is an epithermal/thermal neutron beam. According to the IAEA regulations on neutron beams, it is necessary to ensure that the ratio of epithermal neutrons to fast neutrons in the neutron beam is greater than 20. The influence of fast neutrons was ignored when considering the negligible production of PG rays from fast neutrons.

Fig. 2Full-screen View

(Color online) Relationship of neutron capture cross section between ¹⁰B and ¹H isotopes. The blue and black lines show the capture cross section of ¹⁰B(n, α)⁷Li and ¹H(n, γ)²H, respectively. The red line shows the ratio of ¹⁰B and ¹H neutron capture reactions of the cross section within the range of a thermal to epithermal neutron [27, 28].

With the above assumption, Eq. (5) can be simplified as

Equation (1) can be obtained by adjusting Equation (7). In Eq. (1), K, $M_{{\text{1}}_{\text{H}}}$, and $M_{{\text{10}}_{\text{B}}}$ are all constants, and $N_{0.478}$, and $N_{2.224}$ are the yields of different gamma rays that can be theoretically detected by the three-dimensional PG ray detection system (e.g., using SPECT or a Compton camera) during BNCT. The last remaining parameter, ${\omega }_{{\text{1}}_{\text{H}}}$, is another important aspect for the feasibility of this boron concentration measurement method. According to previous studies, the hydrogen concentration can be inferred through spectral computed tomography (CT) and material decomposition algorithms [22, 23], allowing the boron concentration to be determined by detecting the PG yield of 478 keV and 2.224 MeV PG rays when using the proposed approach. With the development of radiation detection technology [24–26], the boron concentration measurement method has significant potential for clinical applications.

In practical applications, the yield distribution of different PG rays in all voxels should be analyzed to obtain the three-dimensional distribution of the ¹⁰B concentration in the patient.

2.2 Monte Carlo toolkit and configurations

The Monte Carlo simulation software TOPAS 3.1.2 [29] was used to study the multiparticle transport in this study. TOPAS is a particle therapy research oriented Monte Carlo platform based on Geant4 [30–32]. The physical models used in this research include “G4EMStandardPhysics_option4,” “G4HadronPhysicsQGSP_BIC_ AllHP,” “G4DecayPhysics,” “G4IonBinaryCascadePhysics,” “G4HadronElastic PhysicsHP,” “G4StoppingPhysics,” and “G4EmExtraPhysics.” [33, 34].

The component elements of the tumor and normal tissues set in the simulation in this study were all based on the ICRU report 46 [35]. The neutron energy spectrum simulated in this work was obtained from the Massachusetts Institute of Technology (MIT) reactor neutron source [33, 36]. The number of simulated neutrons used during each simulation was 5 × 10⁸, and the calculations for each data were repeated 10 times to determine the statistics of the results. In the simulations, sensitive detectors and filtering scorers were set to obtain the actual PG ray yield information of each voxel under neutron irradiation. During the simulations, each voxel was set as a sensitive detector to detect the number of PG rays generated in each voxel. In addition, a certain energy divergence is considered when recording the PG rays; that is, the yield of each PG ray recorded is the yield of the PG rays within a certain energy range. For PG rays with an energy of 478 keV, the energy range was set as 477 to 479 keV, and the energy range was set as 2.223 to 2.225 MeV for PG rays with an energy of 2.224 MeV.

In this study, three geometric structures were set during the simulations to analyze the different influencing factors on the accuracy of the measured boron concentration, as shown in Figure 3. When analyzing the influence of different neutron components in the incident neutron beam, the geometric structure set during the simulation is as shown as Figure 3a). The mass concentration of hydrogen for all tumor tissues set in this simulation was 0.107. The radius of the tumor was 0.5 cm, and the distribution of the ¹H and ¹⁰B B isotope concentrations in the tumor voxel was uniform. The boron concentration in the tumor tissue was increased from 10 to 100 ppm at 10 ppm intervals, and only one small tumor tissue was set during each simulation. Figure 3b) shows the structure set used to analyze the influence of the tumor depth and surrounding tissue on the accuracy of the measured boron concentration. Two tumor layers were designed at different depths in the tissue. Each tumor layer was composed of 10 tumor tissues with different boron concentrations. The radius of the tumor tissue was set as 0.5 cm, and the depths of the two layers were 3.5 and 5.5 cm from the surface, respectively. At the same depth from top to bottom, the ¹⁰B concentration in the tumor increased from 10 to 100 ppm at 10 ppm intervals. To consider a heterogeneous situation, a bone tissue with a thickness of 1 cm was placed in front of the first tumor layer and the high boron concentration parts (>50 ppm) of the second tumor layer. Based on the structure of the tumor depth and bone tissue set, the neutron spectrum and fluence rate will be different for each tumor voxel. Figure 3c shows the structure used to analyze the influence of the size and heterogeneity of the tumor voxel on the accuracy of the measured boron concentration. To investigate the effect, a spherical geometry composed of two regions (i.e., inner and outer layers) with different ¹⁰B concentrations was set during the simulations. The density of the voxels with different boron concentrations set during this simulation is considered to be unchanged, i.e., 1.04 g/cm³. By adjusting the size of the inner and outer layers, the average ¹⁰B concentration of different voxels was 50 ppm, and only one small voxel was set during each simulation.

Fig. 3Full-screen View

(Color online) Geometric structure used to analyze the influencing factors of boron concentration estimation for different influence factors: a) The structure used to analyze the influence of the neutron energy on the accuracy of boron concentration measured, b) the structure applied to analyze the influence of the tumor depth and surrounding tissue on the accuracy of boron concentration measured, and c) the structure used to analyze the influence of the size and heterogeneity of the voxel on the accuracy of boron concentration measured.

3.1 Influence of neutron energy on the accuracy of the measured boron concentration

In Eq. (1), it can be seen that $N_{0.478}/N_{2.224}$ should only be related to ${\omega }_{{\text{10}}_{\text{B}}}$ and ${\omega }_{{\text{1}}_{\text{H}}}$; however, it can also be affected by other factors. Therefore, this study focuses on analyzing the influence of different factors on $N_{0.478}/N_{2.224}$, which can also clarify the influence on the accuracy of this method. For epithermal/thermal neutrons, a linear correlation can be observed between $N_{0.478}/N_{2.224}$ and ${\omega }_{{\text{10}}_{\text{B}}}$ when ${\omega }_{{\text{1}}_{\text{H}}}$ in a small voxel of interest is known through a theoretical analysis. This finding indicates that K is a fixed value within this energy range. However, the energy range of neutron beams used in clinical BNCT is wide. In high-energy neutrons, K is no longer a fixed value, and the linear correlation may be affected by the high-energy neutrons in the actual BNCT beam. The PG yields corresponding to different neutron components (i.e., thermal neutrons, E < 1 eV; epithermal neutrons, 1 eV < E < 10 keV; fast neutrons, E > 10 keV; and high-energy neutrons, E > 1 MeV) were analyzed to clarify whether K holds a fixed value and determines the effect on the accuracy of the boron concentration prediction.

Figure 4 shows the relationship between $N_{0.478}/N_{2.224}$ and ${\omega }_{{\text{10}}_{\text{B}}}$, where it can be seen that $N_{0.478}/N_{2.224}$ linearly increases with ${\omega }_{{\text{10}}_{\text{B}}}$. Thus, the actual ${\omega }_{{\text{10}}_{\text{B}}}$ can be calculated using $N_{0.478}/N_{2.224}$ for the actual BNCT irradiation.

Fig. 4.Full-screen View

Relation between $N_{0.478}/N_{2.224}$ and ${\omega }_{{\text{10}}_{\text{B}}}$ under uniform ¹⁰B distribution

Table 1 shows the yields of PG rays corresponding to neutrons with different energies (the values in parentheses are counting errors). Under the simulation setting, the proportion of PG rays generated by fast neutrons was only approximately 0.1% of the total PG yield. The proportion of PG rays generated by high-energy neutrons is less than 0.005% of the total PG yield, and the K corresponding to this energy band is no longer fixed. The results indicate that the PG rays produced by fast neutrons will not affect the three-dimensional boron concentration measurement method proposed in this study, and even within this energy region, K is no longer a fixed value assumed in the theoretical derivation.

Table 2 shows the relationships between $N_{0.478}/N_{2.224}$ and the boron concentrations from different neutron components. These results indicate that the $N_{0.478}/N_{2.224}$ produced by thermal neutrons and epithermal neutrons are almost equal to the average value produced by all incident neutrons, and the relative deviations between them are almost less than 1% at the same boron concentration in the tumor tissue. Overall, the results indicate that the yield of different PG rays is affected by the outcomes of the thermal and epithermal neutrons. The simulation results prove that although K is not a constant assumed in the theoretical derivation for the actual BNCT neutron beams, the proportion of PG rays and capture cross sections of different elements with high-energy neutrons are sufficiently small. Consequently, these factors did not have a significant negative effect on $N_{0.478}/N_{2.224}$.

3.2 Influence of tumor depth and surrounding tissue on the accuracy of the measured boron concentration

In the actual BNCT treatment, the neutrons are scattered when neutrons travel through human tissue, and the energy and flux of neutrons reaching the tumor area will become complicated. The yield of different PG rays and the energy deposited in the tumor voxels by complex neutron beams will be significantly different from the results produced by the ideal neutron beams. Thus, in this section, the influences of the tumor depth and surroundings on the accuracy of the measured boron concentration are analyzed.

Based on the geometric structure shown in Fig. 3b, the effect of tissue structure is studied by comparing the PG ratios produced in different tumor layers, the results of which are shown in Fig. 5. The generation of 2.224 MeV PG is only related to the neutron fluence rate and ¹H concentration; therefore, the neutron fluence rate at the middle position is high, and the two ends are low because of the influence of neutron scattering from the $N_{2.224}$ distribution. By contrast, the $N_{0.478}$ distribution is completely different from the $N_{2.224}$ distribution. With the changes in location and ¹⁰B concentration of the tumor, $N_{0.478}$ gradually changes. Based on the effect of the neutron fluence rate and ¹⁰B concentration, $N_{0.478}$ of the tumor with a ¹⁰B concentration of 80 ppm is the highest. The results in the second layer are similar to those in the first layer; however, the yield of PGs with the same ¹⁰B concentration is significantly reduced because of the overall decrease in the neutron fluence rate. The comparison results of $N_{0.478}/N_{2.224}$ of the two layers (Fig. 5b) indicate that even if tumors are in different neutron fields, $N_{0.478}/N_{2.224}$ of the two layers are nearly the same, and the deviations are nearly less than 1% when the values of ${\omega }_{{}_{}{}^{10}B}$ are the same. This result shows that the changes in the depth of the tumor and the tissue structure will only affect the yield of the PG rays and not the ratio ($N_{0.478}/N_{2.224}$) for a voxel with a uniform ¹⁰B. Therefore, the method of inferring the ¹⁰B concentration based on the ratio of dual-energy PG yields still applies for tumors at different depths in the human body.

Fig. 5Full-screen View

(Color online) Influence of depth and tissue structure on $N_{0.478}/N_{2.224}$: (a) PG yield under different configurations with boron and hydrogen concentration ratios. The red and blue points in the figure represent the PG yield of the first and second tumor layers, respectively. The dot and square points represent $N_{0.478}$ and $N_{2.224}$, respectively. (b) $N_{0.478}/N_{2.224}$ for different boron and hydrogen concentration ratios and the deviation of different tumor layers. The purple point represents $N_{0.478}/N_{2.224}$ from the first tumor layer. The purple line represents $N_{0.478}/N_{2.224}$ from the second tumor layer. The red line represents the relative deviation of these two values.

3.3 Influence of the size and heterogeneity of the voxel on the accuracy of the measured boron concentration

In the above study, the distribution of boron concentration in the tumor tissue is even, and thus the influence of boron concentration distribution was ignored. In this case, the average boron concentration of the tumor tissue was equal to the boron concentration of each tumor voxel. In clinical studies, the distribution of boron drugs in patients shows individual differences, and the absorption capacity of boron drugs is also different inside the patient. These factors lead to a non-uniform distribution of boron drugs in the ROI. In this case, the average boron concentration of tumor tissue is not equal to the boron concentration in each tumor voxel, and the yields of different PG rays are different in different tumor voxels. Thus, the influence of boron distribution in a voxel on the accuracy of boron concentration measurements should be clarified.

Figure 6a shows $N_{0.478}/N_{2.224}$ with different voxel sizes and different compositions of the ¹⁰B concentration. For a uniform distribution of the boron concentration (i.e., where the boron concentration of the inner and outer layers are both 50 ppm), $N_{0.478}/N_{2.224}$ is nearly unchanged as the voxel size increases. Therefore, the voxel size does not affect $N_{0.478}/N_{2.224}$ when the boron concentration is evenly distributed in the voxels. However, when the distribution of boron concentration in a voxel is uneven, the greater the unevenness of the boron concentration distribution is at the same voxel size, the greater the difference in $N_{0.478}/N_{2.224}$. For example, when the voxel size is 2.8 cm, the results of 10/90 ppm and 40/60 ppm are 0.3344 and 0.3089, respectively. In addition, in the case of a non-uniform distribution of boron concentration in the voxel, the voxel size will also affect $N_{0.478}/N_{2.224}$. For example, when the boron concentration in the voxel is composed of 10/90 ppm, $N_{0.478}/N_{2.224}$ is increased from 0.2979 to 0.3344 when the voxel size is increased from 0.4 to 2.8 cm. These changes in $N_{0.478}/N_{2.224}$ indicate that the deviations between the boron concentration measured by this method and the true boron concentration reached more than 10% when the boron distribution in the voxel was non-uniform. Moreover, when the voxel size is large and the degree of non-uniformity of the ¹⁰B concentration distribution is large (i.e., the voxel size is 0.8 cm and is composed of 40/60 ppm), the deviation between the actual ¹⁰B concentration and the average ¹⁰B concentration obtained using the PG yields of the entire voxel is considerable. Therefore, when using this method in clinical applications, selecting a reasonable voxel size is important to ensuring the accuracy of the boron concentration obtained. Figure 6b shows the deviations of $N_{0.478}/N_{2.224}$ of different voxels with different concentrations, and uniform ¹⁰B distribution results were selected as a reference. Considering the calculation error, when the deviation is less than 2%, this result is considered to be approximately equal to the reference. When the voxel size is 0.4 cm, the deviations of the results for the four ¹⁰B concentration compositions are all less than 2%. In some cases, the deviations are greater than 2% when the voxel size increases. For example, when the boron concentration in the voxel is 10/90 ppm, the deviation of the $N_{0.478}/N_{2.224}$ is increased from 0.60% to 3.02% when the voxel size is increased from 0.4 to 0.8 cm. In actual clinical practice, selecting the optimal voxel size based on different treatment plans is the key to obtaining an accurate boron concentration.

Fig. 6Full-screen View

(Color online) Influence of voxel sizes and boron concentration distribution of tumor voxels on $N_{0.478}/N_{2.224}$: (a) $N_{0.478}/N_{2.224}$ under different boron concentration compositions and voxel sizes and (b) deviations between $N_{0.478}/N_{2.224}$ with non-uniform and uniform boron distribution of tumor voxels.