MagicDose program development & framework

Dose calculation is the core part of the whole BNCT treatment planning system, which is based on the MC particle transport method. MagicDose, a dose calculation program dedicated to the BNCT treatment planning system, is developed in C++, which is convenient for future development and maintenance of the BNCT treatment planning system. The MagicDose program consists of seven functional modules, including the Geometry module, Source module, Database module, Particle Transport module, Tally module, Output module, and Auxiliary module. Figure 3 illustrates the MagicDose development framework.

Fig. 3Full-screen View

(Color online) Development framework of MagicDose

Fig. 4Full-screen View

(Color online) (a) Lattice in input card; (b) Model

Figure 5(a) shows the schematic of an in-hospital neutron irradiator type I reactor (IHNI-I), which is the first prototype in the world that meets the requirements of the IAEA-specialized neutron source device BNCT [20]. The site selection is located approximately 2 kilometers from the Beijing Nuclear Industry Hospital rather than being directly installed within the hospital premises. This decision stems from the consideration of various uncertainties, such as the distribution and control of radiation dose levels, with the aim of accumulating experience and establishing a technical foundation for the future construction and application of BNCT neutron sources [21]. The IHNI-I has a “tank-pool” structure with a rated power of 30 kW and two beam holes for thermal and epithermal neutrons. To avoid the significant computational time required to track particles from the reactor core for simulating the dose distribution in the human body, MagicDose establishes equivalent surface sources at the exits of the two beam flow orifices of IHNI-I to facilitate rapid dose calculations. The equivalent surface source, denoted as φ(R, E, μ), is a function of the spatial position, energy, and direction, where R represents the radial distance (cm), E denotes the energy (MeV), and μ is the cosine of the scattering angle (ranging from -1 to 1). Given that the collimator aperture diameter of IHNI-I is 12 cm, to more accurately capture the particle information emitted from the beam ports in constructing the equivalent surface source, this study divides the radial direction R into three intervals: 0<R<6, 6<R<10, and 10<R<15, as shown in Fig. 5(b) [22, 23]. The energy spectra shown in Fig. 5(c) originate from an equivalent surface source at the thermal neutron beam port of IHNI-I, whereas the energy spectrum in Fig. 5(d) originates from the equivalent surface source at the epithermal neutron beam port of IHNI-I [21].

Fig. 5Full-screen View

(Color online) (a) Schematic of in-hospital neutron irradiator type I reactor (IHNI-I); (b) Schematic of the IHNI-I equivalent surface source; (c) Energy spectrum of the equivalent surface source at the thermal neutron beam port of IHNI-I; (d) Energy spectrum of the equivalent surface source at the epithermal neutron beam port of IHNI-I

Fig. 6Full-screen View

(Color online) Kerma value for each dose component

Fig. 7Full-screen View

(Color online) Particle transport flow diagram

The MagicDose program adopts the particle transport module as the core, interconnects with the upstream geometry, source, and database modules, and simultaneously supports the downstream tally and output modules. Each module is independent of the others, and module functions are realized by messages passing through the interface, which is conducive to the maintenance and expansion of the program.

Calculation conditions and models

To verify the accuracy and efficiency of MagicDose, three models are selected for testing. First, because the MCNP program can be used as a standard for all types of MC calculations and is often used to calibrate the correctness tests of other MC programs, the results of MagicDose calculations are compared with MCNP. The two programs select the voxel phantom of a modified Snyder head with 16 and 8 mm resolution, which is commonly used internationally. Under the same irradiation conditions, we compare the calculated spatial neutron and photon dose rates to see if they are consistent and count the computational efficiency of both. Subsequently, to explore the influence of voxel phantom on the calculation results under different spatial resolutions, the MCNP computational analytical model is used as the reference result. The voxel phantom of the modified Snyder head with a tumor having three different spatial resolutions of 16, 8, and 1 mm was constructed using MagicDose combined with the material-located method using the central point algorithm. The depth-dose-rate curves and spatial value distribution plots are used to analyze the effect of model voxel size on the results and to statistically calculate the computational efficiency of MagicDose with different resolutions and the number of threads. Finally, based on the CT head model of a patient given by DICOM data, we initially verify whether MagicDose can be applied to the clinical case of BNCT.

All the calculations in this study are performed on an Intel (R) Core(TM) i7-10700 K CPU @ 3.8 GHz with 32 GB RAM and 16 threads. To minimize the impact of the results from the deep penetration of the MC program, the number of particles is set to 1 × 10⁸ for both MagicDose and MCNP to ensure that the number of samples is adequate [32, 33]. In the actual BNCT treatment, some photons are produced by the accompanying neutron source, which cannot be avoided. We assume the absence of accompanying photon contamination from the neutron source and only discuss the induced photons produced by the reaction of neutrons with nuclides in the voxel.

(1) Voxel phantom of the modified Snyder head with 16 and 8 mm

Fig. 8Full-screen View

(Color online) Schematic of modified Snyder head

The surface analytical equations for the modified Snyder head phantom with tumor at each interface are expressed as follows:

Boundary surfaces of the scalp and air(5)

Boundary surfaces of the skull and scalp(6)

Boundary surface between brain and skull(7)

In clinical treatment, owing to the different shapes, locations, and sizes of patients’ tumors, it is not possible to use analytical equations to describe the phantom; therefore, it is common to use a continuous and uniform cubic mesh to approximate the phantom structure in international clinics, which are generally called voxel phantoms. Goorley et al. [3] constructed 16 and 8 mm voxel phantoms as an international benchmark problem based on the modified Snyder’s analytical model, which is mixed with four basic materials (brain, skull, scalp, and air); the volume share of each basic material in the mix is a multiple of 10% to produce the mixed materials and corresponding mixing densities, and a total of 286 hybrid materials are derived. Therefore, to verify the accuracy of MagicDose, identical 16 and 8 mm voxel phantoms are constructed based on MagicDose and MCNP respectively (Fig. 9).

Fig. 9Full-screen View

(Color online) Voxel phantom of modified Snyder head with a 16 mm, b 8 mm

For the irradiation condition, in this paper, the international wide-spectrum epithermal mixed neutron beam proposed for the modified Snyder head phantom is adopted, in which the thermal neutron with energy less than 0.5 eV accounts for 10%, the epithermal neutron with energy between 0.5 eV and 10 keV accounts for 89%, and the fast neutron with energy between 10 keV and 2 MeV accounts for 1%. The overall neutron source spectra obeyed the distribution according to the energy-dependent probability function 1/E, uniformly distributed on a disk surface with a radius of 5 cm and a source strength of 10¹⁰ n/cm² ·s [3]. Figure 10 presents the energy spectrum of the neutron source.

Fig. 10Full-screen View

(Color online) Energy spectra of the wide-spectrum epithermal mixed neutron beam

(2) Voxel phantom of the modified Snyder head with tumor based on the material-location method with a central point algorithm

To study the influence of the voxel phantom on the calculation results under different spatial resolutions and to simulate the situation of glial tumors in the brain in the clinic. Based on the model illustrated in Fig. 8, a sphere with a radius of 1.5 cm is added to the brain tissue as the tumor boundary, and the coordinate origin is at the center of the entire analytical model. Fig. 11 shows the modified Snyder head phantom with tumor [34, 35], in which the equation for the boundary surface of the tumor and brain is given by Eq. (8).(8)

Fig. 11Full-screen View

(Color online) Schematic transverse of modified Snyder head with tumor

For the analytical model, the materials can be accurately described, but the voxel phantom will be difficult to deal with because the voxel mesh may be at the interface of the two substances. The traditional way to deal with this is to fill by obtaining the exact volume ratio of different materials inside each cubic mesh (Fig. 9), which is extremely time consuming and results in bias when the number of voxel meshes is large. As described above, the model material was derived 286 from 4 materials as the hybrid materials. In this study, to reduce the complexity of the simulation problem, the material of the center point of the mesh is used in the phantom construction of MagicDose, and the material density of the whole mesh is determined according to the location of the center point of each mesh, which is reduced from the original 286 kinds of material to four types of materials [36], Fig. 12 shows the material-location method with a central point algorithm. In Fig. 12, (X_m, Y_m, Z_m) denotes the center point position of each voxel mesh, and f₁(X_m, Y_m, Z_m), f₂(X_m, Y_m, Z_m), f₃(X_m, Y_m, Z_m) and f₄(X_m, Y_m, Z_m) denote the boundary equations Eq. (5)–(8) as follows: When the position of the center point is substituted in the boundary equations with a value less than 1, it is “True”, the mesh material is the tissue component inside the boundary equation, otherwise it is “False”, the mesh material is the tissue component outside the boundary equation.

Fig. 12Full-screen View

(Color online) The material-located method with central point algorithm

Based on the material-located method with the central point algorithm, MagicDose constructed three voxel phantoms of the modified Snyder head with tumor of different spatial resolutions (16, 8, and 1 mm), and the corresponding voxel mesh numbers of M_x×N_y×L_z are 2744, 21952, and 11239424, respectively [37]. As shown in Fig. 13 show transverse views of different voxel phantoms.

Fig. 13Full-screen View

(Color online) Transverse view of voxel phantom of modified Snyder head with tumor at (a) 16 mm, b 8 mm, c 1 mm spatial resolution

In this study, MCNP is used to construct a modified Snyder head analytical model with a tumor. The results of the MCNP computation are used as the reference result, which is compared and verified with the results of MagicDose under different resolutions, and are used to explore the influence of voxel phantom on the computational results under different spatial resolutions. The tissue compositions of the voxel phantom of the modified Snyder head with tumor and the modified Snyder head analytical model with tumor are shown in Table 3. The tissue compositions of the tumor are consistent with that of normal brain tissues, except for the different boron contents, according to the relevant literature [38]. According to the current clinical trials requirements, the ¹⁰B content (T) of the tumor must be higher than the ¹⁰B content (N) of the normal tissue by a factor of 2.5 or more (T/N > 2.5), and the calculations in this paper set the boron concentration of scalp tissue and normal tissue in the brain at 10 ppm and the boron concentration within the tumor tissue at 30 ppm, with no boron in the skull [39]. MagicDose and MCNP both utilize the wide-spectrum epithermal mixed neutron beam spectrum shown in Fig. 10 for the simulated dose calculations. To make the calculated results comparable to international BNCT dose results, according to Eq. (4). The RBE and CBE factors are selected with reference to the values obtained from the MIT and BNL in BNCT clinical trials [40, 41], as shown in Table 4.

(3) CT head model of a patient based on DICOM data

In a preliminary application, MagicDose was used to construct a voxel phantom based on DICOM data from a patient’s head (Fig. 14 shows coronal, transverse, and sagittal CT images [42]. The dimensions of the phantom in the X-, Y-, and Z-directions are 18, 18 and 16.5 cm, respectively, and the size of the constructed voxel is 1 mm×1 mm×1 mm. An equivalent source of the epithermal neutron beam plane of the IHNI-I reactor with the MagicDose source module is used for the irradiation conditions (Fig. 5(d)), and the irradiation direction was vertical to the back of the head (Fig. 14(g)).

Fig. 14Full-screen View

(Color online) (a) Coronal, (b) Transverse, and (c) Sagittal of CT images; (d) Coronal, (e) Transverse, and (f) Sagittal of MagicDose constructed model; (g) Schematic of irradiation direction

Correctness testing of MagicDose

Figure 15 shows comparison plots of the thermal neutron dose rate, epithermal and fast neutron dose rates, and induced photon dose rate values for all voxels in the voxel phantom of the modified Snyder head at spatial resolutions of 16 and 8 mm for MagicDose and MCNP. Using MCNP as a reference result [3], the slope ratio method is used in this study to distinguish how close the dose rate results calculated by MagicDose are to the MCNP results. For any voxel n in the model, the dose rate result of MCNP is taken as the horizontal coordinate value Xn, and the dose rate result of MagicDose is the vertical coordinate value Y_n, and the points (X_n, Y_n) are labeled sequentially on the coordinate system, and it is observed that all the images formed by scattering the points are scattered near the straight line equation of y=x. As illustrated in Fig. 15, the number of model voxels with a spatial resolution of 16 mm is fewer than that of the 8 mm model voxels. Consequently, the number of coordinate points in the 16 mm model are lower than that in the 8 mm model. However, the points (X_n, Y_n) generated by MagicDose and MCNP are predominantly concentrated near the line y=x. Furthermore, the locally magnified graphs of the respective dose rates reveal that the deviation of the dose rate results from y=x is minimal. This observation indicates that the results computed by MagicDose and MCNP are in strong agreement, thereby validating the accuracy of MagicDose.

Fig. 15Full-screen View

(Color online) (a)(b) Thermal neutron dose rate based on voxel phantom of modified Snyder head; (c)(d) Epithermal and fast neutron dose rate based on voxel phantom of modified Snyder head; (e)(f) Induced photon dose rate based on voxel phantom of modified Snyder head

In the model calculations, both MagicDose and MCNP are performed under the same conditions, as listed in Table 5. At a spatial resolution of 16 mm, the calculation times for MCNP and MagicDose are 1373 and 944 s, respectively, with MagicDose achieving an efficiency improvement of approximately 31.24% compared to MCNP. At a spatial resolution of 8 mm, the calculation times for MCNP and MagicDose are 2042 and 1457 s, respectively, with MagicDose demonstrating an efficiency improvement of approximately 28.65% over MCNP. The numerical results demonstrate that MagicDose outperforms MCNP in terms of computational efficiency, while ensuring the accuracy of the BNCT calculations.

Exploration of MagicDose at different spatial resolutions

The result of each dose rate of the modified Snyder head analytical model with tumor simulated by MCNP calculation is used as the reference in this section. In the comparative validation calculations, using the parameter setup described above (including material and source terms, etc.) as well as the nuclear cross-section database, MagicDose is used to construct voxel phantoms with three different spatial resolutions of voxel sizes of 16, 8, and 1 mm, and to compare the variability in the results of the phantom-depth correlation dose-rate curves along the Z-axis computed by the two MC particle transport programs. Simultaneously, the relevant spatial dose rates of the head phantom at different voxel resolutions are statistically calculated to explore the effect of the voxel size of the voxel phantom on the results and to provide reference opinions for the optimization of the voxel phantom design for clinical treatment.

Figure 16, 17, 18 show the boron dose rate, thermal neutron dose rate, epithermal and fast neutron dose rates, induced photon dose rate, total relative biological effect dose rate, relative deviation with depth, and spatial dose rate distributions at different voxel resolutions. The relative deviation is the percentage of deviation obtained by comparing the MagicDose voxel phantom results with the MCNP analytical model results as the reference result, and the ICRU24 report suggests that the relative deviation between the actual treatment and the planned dose of BNCT should not exceed 5%; otherwise, the risk of tumor recurrence persists, and radiological complications increase in normal tissues [43].

Fig. 16Full-screen View

(Color online) (a)(b) Boron dose rate and relative deviation vs. depth; (c)(d)(e) Boron dose rate at 16, 8, and 1 mm; (f)(g) Thermal neutron dose rate and relative deviation vs. depth; (h)(i)(j) Thermal neutron dose rate at 16, 8, and 1 mm

Fig. 17Full-screen View

(Color online) (a)(b) Epithermal and fast neutron dose rate and relative deviation vs. depth; (c)(d)(e) Epithermal and fast neutron dose rate at 16 mm, 8 mm and 1 mm; (f)(g) Induced photon dose rate and relative deviation vs. depth; (h)(i)(j) Induced photon dose rate at 16 mm, 8 mm and 1 mm

Fig. 18Full-screen View

(Color online) (a)(b) Total relative biological effect dose rate and relative deviation vs. depth; (c)(d)(e) Total relative biological effect dose rate at 16 mm, 8 mm and 1 mm

As shown in Fig. 16(a), the boron dose rate calculated by MagicDose at different resolutions of the voxel is compared with the reference result, which is consistent with the results at 1 mm. However, the variability of the results is more obvious with an increase in the resolution of the voxel. Fig. 16(b) shows that the relative deviations of 1 mm are within ±5% of the clinical allowable for BNCT, whereas 8 and 16 mm have greater relative deviations between depths of 0.5–1.5 cm, 5.4–8.2 cm, and 14.2–18.2 cm. By integrating the phantoms depicted in Fig. 13 with the spatial dose-rate distribution derived from Fig. 16(c)–(e), it becomes evident that significant relative deviations can be ascribed to the depth ranges of 0.5–1.5 cm, 14.2–18.2 cm, which corresponds to the skull, and 5.4–8.2 cm, which pertains to the tumor. The phantoms with dimensions of 8 mm and 16 mm lack high-resolution representation of tissues in these specific regions, ultimately resulting in inaccuracies in the dose rate calculations. Figure 16(c)–(e) demonstrate that the results derived from the voxel phantom at high resolution exhibit greater accuracy. Analysis of the spatial dose rate values clearly indicates that the skull does not contribute to the boron dose rate, as it lacks boron content. Conversely, in the tumor region, the elevated boron content leads to a substantial increase in the boron dose rate, thereby facilitating more precise tumor localization in clinical settings.

From Fig. 16(f)–(j) and Fig. 17(f)–(j), we observe that the thermal dose and induced photon dose rates for 1 and 8 mm are consistent with the reference result, and the relative deviations are limited to ±5%. However, the calculated results are biased owing to the large resolution of the 16 mm voxel. The distribution curves of the thermal dose rate and induced photon dose rate exhibit the same trend, because the induced photon comes from the neutron capture reaction and the γ-ray is generated by the capture reaction between thermal neutrons and H-atoms in the body; the photon is generated in place of the higher neutron dose, and the greatest damage is produced in the superficial tissues of the head.

This is shown in Fig. 17(a), (c), (d), and (e), because 89% of the neutron source radiation is epithermal neutrons, neutrons slowing down to thermal neutrons occur initially and the dose rate gradually decreases as the depth increases. From Fig. 17(b), it is understood that when the depth is less than 12 cm, the relative deviations of both 1 mm and 8 mm are basically kept within 5%; when the depth is greater than 12 cm, it will cause the problem of deep penetration, at which time the MC program calculates the number of particles arriving at the depth to be less, resulting in the statistical results of the relative deviation to be larger.

The dose rates are shown in Fig. 16-17 are combined according to the corresponding factors of CBE or RBE to obtain the total dose rate of the relative biological effect, as shown in Fig. 18. As shown in Fig. 18(a)(b), the total dose rate at 1 mm is consistent with the reference result, and the relative deviations are kept within 5%; the discrepancy of the results at 8 mm and 16 mm is due to the boron dose as the main dose component, and the effect of the resolution of the voxel phantom on the boron dose, which ultimately leads to the deviation of the total dose rate. Overall, Fig. 16, 17, 18 show that the dose rates and relative deviations of the high-resolution voxel phantoms calculated from MagicDose are consistent with the reference result and remain within ±5%, indicating that as the size of the constructed voxel mesh gets smaller, the better it converges to the analytical model, the better the computed results match, highlighting the effect of the voxel size of the voxel phantoms on the results

Because the boron concentration of BNCT changes rapidly with the body’s metabolism after boron injection, BNCT is highly demanded owing to its speed of dose calculation, which ensures that the dose calculation can be completed in a very short period of time (<3600 s). The number of voxels also affects the efficiency of the calculation. Therefore, to verify the computational efficiency of MgaicDose, simulations are performed using 1, 4, 8, and 16 threads for the voxel phantom of the modified Snyder head with tumor at 16, 8, and 1 mm, respectively, to test the computational efficiency of the program at different spatial resolutions and with different numbers of threads. Figure 19 compares the computation time (s) for different mesh scales and different numbers of thread 3D histograms. Evidently, under the same computational conditions, the smaller the voxel size, the longer the simulation computation time. However, with an increase in the number of parallel threads, the time used for the simulation becomes increasingly shorter, with increasingly obvious acceleration effects.

Fig. 19Full-screen View

(Color online) Comparison of computation time (s) for different mesh scales and different number of threads 3D histograms

Initial application of MagicDose based on DICOM data

MagicDose constructes a voxel phantom of the head using DICOM data, and Fig. 14 shows that the schematics of MagicDose in the coronal, transverse, and sagittal planes are consistent with the CT images of the corresponding locations. Based on the irradiation from an equivalent surface source at the epithermal neutron beam port of IHNI-I, as shown in Fig. 15(d), Fig. 20(a), (b), (c), and (d) show the boron dose rate, thermal neutron dose rate, epithermal and fast neutron dose rates, and induced photon dose rate, respectively. The total relative biological effects of the dose rate (Fig. 20(d)) is derived by integrating Table 4 with the aforementioned dose-rate values. The spatial distributions of the boron and thermal neutron dose rates are similar to those shown in Fig. 20(a), and (b) because both doses are related to thermal neutrons, and the boron dose rate is significantly higher than the contribution of the epithermal neutron dose rates presented in Fig. 20(c) is mainly located in the shallow region of the posterior cerebral epidermis. Figure 20(d) shows that the induced photon dose rate is distributed throughout the space and presents the complete human head profile; combining the CBE and RBE factors, the total relative bioeffect dose rate shown in Fig. 20(e) is obtained, and the highest dose region appears in the irradiated region of the back of the head.

Fig. 20Full-screen View

(Color online) (a) Boron dose rate distribution based on head CT; (b) Thermal neutron dose rate distribution based on head CT; (c) Epithermal and fast neutron dose rate distribution based on head CT; (d) Induced photon dose rate distribution based on head CT; (e) Total relative biological effect dose rate based on head CT