Monte Carlo simulations

The HIMM has four treatment terminals with active and passive beam delivery methods. The treatment terminals were used for irradiation from different directions: horizontal, vertical, horizontal + vertical, and 45° inclination. The dose rate at the isocenter was 0.001-1 $\text{Gy}/\text{s}$ ^[2,20]. The maximum treatment depth in human tissue is 27 cm, which can be used for whole-body irradiation. The main hardware components in the nozzle include a primary collimator, ridge filter unit, range shifter unit, MLC, compensator, and two monitors (dose and beam position detectors). The beam was scanned by the scanning magnets, expanded into a 15 $\times$ 15 ${\text{cm}}^2$ irradiation field in the lateral direction, and then collimated by the primary collimator after passing through the vacuum window. When passing through the ridge filter, the pristine Bragg peak spread in the longitudinal direction, forming an SOBP. After the beam passed through the range shifter to modulate the range, the MLC was used for the lateral conformation. The final three-dimensional distribution of the high-dose zone was consistent with the shape of the target region. The use of the compensator made the dose distribution on the trailing edge of the target area steeper, which helped to protect endangered organs. The shape and position of the beam were measured using a beam position detector, and the intensity of the beam was monitored using a dose detector.

The beamline setting of the HIMM Room 2 in the Monte Carlo simulation is shown in Fig. 1. Prior research has shown that $H/D$ values are strongly influenced by the presence of the patient or phantom, phantom size, anatomical location of the treatment site, and orientation of the patient or phantom ^[16,21]. To eliminate the effects of internal neutrons and focus on machine-only behaviors, the free-in-the-air (without phantom) option was used in most of the simulations. The range shifters and compensators were simplified owing to their patient specificity. In addition, the monitors were not included in the simulation because their effects on secondary neutrons are negligible owing to their thinness compared to other components ^[17]. The materials and corresponding distances of all the hardware components to the isocenter included in the calculational geometry are listed in Table 1.

Fig. 1Full-screen View

Model of the passive beam delivery system

The neutron spectral fluence was tallied into a set of 12-cm-diameter spherical receptors. A total of 16 receptors were located at the isocenter at distances of 25, 50, 100, 150, and 200 cm from the isocenter and angles of 0, 45, and 90° for the incident beam direction (see Fig. 1) ^[15-16]. The neutron spectral fluence $\text{Φ}\left(E\right)$, was tallied into 480 logarithmically spaced bins from $1$ to $1\text{ GeV}$. In most simulation settings, a closed MLC was performed to obtain a conservative $H/D$ value. The SOBP width was set to 6 cm, which is the standard setting for particle therapy ^[22].

The secondary neutron fluence spectrum for the HIMM passive beam delivery system was calculated using the GATE v8.2/Geant4-10-05. Geant4 (geometry and tracking) is a set of Monte Carlo toolkits developed by CERN, using object-oriented technology based on C++, and its purpose is to simulate the physical process of particle transportation in matter ^[23]. Owing to its flexibility, Geant4 is also widely used in high-energy physics, nuclear physics, space and astrophysics, and radiation protection and detection, and it can be applied to the field of medical physics with the addition of expansion packages. GATE was designed by the OpenGATE team, based on the Geant4 kernel. This simulation platform is mainly used for Monte Carlo simulations of computed tomography and can also simulate particle radiotherapy.

The electromagnetic and hadronic physical processes were implemented using the physics list QGSP_BERT_HP_EMV, according to the recommendation from the OpenGATE collaboration. The cut (production threshold) for particle transport was 1 mm. $1\times {10}^7$ carbon-ion histories were tracked for each simulation. Statistical errors were less than 5 % for most receptors. All calculations were conducted on a CPU server (Intel® Xeon® CPU E5-2690 V3 @2.60 GHz $\times$ 48, OS: Scientific Linux release 7.6) and parallel [48] was used to improve computing efficiency.

Calculation of ambient dose equivalent and $H/D$ values

To assess neutron exposure at various locations, the focus was first on the ambient dose equivalent obtained using the fluence spectrum and neutron fluence-to-ambient dose-equivalent conversion coefficients. $H^{*}\left(10\right)={\displaystyle \sum _{i=1}^{\max }{\text{Φ}}_i\cdot h_i},$ (1) where $i$ is the $i$th neutron energy bin in the neutron spectral fluence, ${\text{Φ}}_i$ the fluence in the $i$th energy bin, and $h_i$ is the conversion coefficient of the $i$th neutron energy. Figure 2 shows the response function of WENDI-II, a typical wide-energy neutron detection instrument, together with the ICRP74 fluence-to-ambient dose conversion coefficients and NCRP38 fluence-to-dose equivalent conversion coefficient ^[24-27]. In this study, we chose the coefficients recommended by the ICRP74 to convert the dose. However, its maximum energy is 201 $\text{MeV}$; hence, the calculated results for the neutron region between 201 and 1000 $\text{MeV}$ are shown in Fig. 2. These were generated using the Monte Carlo code Fluka by Ferrari et al.

Fig. 2Full-screen View

(Color online) Response function of WENDI-II with fluence-to-ambient dose conversion coefficients by ICRP74, Ferrari et al., and the NCRP38 fluence-to-dose equivalent conversion coefficient

The $H/D$ value is defined as the ambient dose equivalent per therapeutic dose and is calculated using Eq. (2). It is an important quantitative index for evaluating secondary radiation levels around the patient. $\frac{H}{D}=\frac{H}{p}÷\frac{D}{p},$ (2) where $p$ is the number of carbon ions, and $D$ is the therapeutic dose. The biologically effective dose was obtained by relative biological effectiveness (RBE), weighting the physically absorbed dose on the beam central axis at the depth corresponding to the center of the SOBP. Table 2 lists the RBE values corresponding to the depth at the center of the SOBP under different beam energies. The absorbed dose was determined from separate simulations in which the MLC aperture size was fixed at 10 $\times$ 10 ${\text{cm}}^2$ and a water phantom was introduced. The water phantom (size: 30 $\times$ 30 $\times$ 35 ${\text{cm}}^3$) was located upstream at the isocenter, normal to the incident beam.

$H/D$ value analytical model from Newhauser and co-workers

The spatial distribution law of $H/D$ values in the treatment room has been reported in previous studies ^[15], as follows: ${\left(\frac{H}{D}\right)}_d\propto d^{-\alpha },$ (3) where $d$ is the distance between the position of the detector and the effective neutron point source, and $\alpha$ is the power exponent in which the $H/D$ value decreases as the distance increases. The physical meaning of an effective neutron point source is a point in the MLC where most of the carbon ions have been blocked.

The effective neutron point source distance d is defined as: $d=\sqrt{（d_{\text{iso}}+r\cdot \text{cos} \theta {）}^2+r^2\cdot {\text{sin}}^2\theta },$ (4) where $d_{\text{iso}}$ is the distance from the effective neutron point source to the isocenter of the treatment room, and $r$ and $\theta$ are as shown in Fig. 1.

Inspired by this conclusion, Zheng et al. proposed an ambient dose-equivalent analytical model for a treatment room ^[16]. Using Eq. (5), the $H/D$ values at different spatial positions at distance $d$ from the effective neutron point source can be estimated. ${\left(\frac{H}{D}\right)}_{\text{d},\text{a},\text{m}}={\left(\frac{H}{D}\right)}_{d_{\text{iso}},a_{\text{ref}},m_{\text{ref}}}{\left(\frac{d}{d_{\text{iso}}}\right)}^{-\alpha }{C}_{\text{a}}{C}_{\text{m}},$ (5) where $a_{\text{ref}}$ and $m_{\text{ref}}$ are the aperture sizes of the MLC and SOBP, respectively. The reference values of $a_{\text{ref}}$ and $m_{\text{ref}}$ were 0 ${\text{cm}}^2$ and 6 $\text{cm}$, respectively. The coefficients $C_{\text{a}}$ and $C_{\text{m}}$ were introduced to indicate the effects of the field size and SOBP width on the $H/D$ value.

Characteristics of secondary neutron fluence spectrum and ambient dose equivalent spectrum

The simulated neutron fluence spectrum and its corresponding ambient dose-equivalent spectrum for beam energies of 165, 200, 250, 300, 350, and 400 $\text{MeV}/\text{u}$ at the isocenter are shown in Fig. 3. The simulation results are shown in a linear scale versus the logarithm of neutron energy, which ensures that the relative contribution from each bin is proportional to its area. As expected, the neutron fluence increased with the beam energy. The fluence spectrum contains two pronounced peaks: one low-energy evaporation peak with a position of approximately 1 $\text{MeV}$ and one high-energy direct neutron peak that ranges from 15 $\text{ MeV}$ to twice the energy of the carbon ions. The distribution law of the dose-equivalent spectrum was analogous to that of the fluence spectrum. However, the contribution of low-energy neutrons decreases, and that of high-energy neutrons increases. This is explained by the fact that the ambient dose-equivalent conversion coefficient of neutrons with energy less than 100 $\text{keV}$ is very low. Up to 23 % of the ambient dose equivalent was from neutrons with energy ≥ 15 $\text{MeV}$.

Fig. 3Full-screen View

(Color online) (a) neutron fluence spectrum and (b) neutron ambient dose equivalent spectrum at the isocenter for different carbon ion beam energy.(a)

(SOBP width is 6 cm, MLC is closed)

These results are consistent with previously published results ^[5,28], except for the absence of a thermal neutron peak at approximately ${10}^{-7}\text{ MeV}$. This difference might be related to the simulated parameter settings; the peak can be attributed to scattered neutrons (due to room return), whereas in this study, the shielding wall was simplified. Owing to the low fluence and small corresponding conversion factors of this peak, the effect of its absence on the ambient-dose equivalent calculations is negligible.

Effect of spatial distance on the $H/D$ value

The $H/D$ values for different energies for the 6 $-\text{cm}$ SOBP width, and closed MLC options were calculated and fitted according to the power-law formula. It can be found that the distribution of the $H/D$ value at different spatial locations conforms to the power-law distribution, and the $H/D$ value decreases rapidly with the increase in distance (Fig. 4a). Some data points deviated from the fitting curve; therefore, the data in different directions were normalized and fitted separately. The fitting parameters and goodness-of-fit values are listed in Table 3.

Fig. 4Full-screen View

(Color online) Monte Carlo simulation $H/D$ values and fitting curves (a) at all locations, (b) in the 0° direction, (c) the 45° direction, (d) the 90° direction, and (e) the 90° direction with a linear scale (SOBP width: 6 $\text{cm}$, MLC: closed)

As shown in Fig. 4, the $H/D$ values in all directions fit well with the power-law distribution. However, the $H/D$ values of the locations in the 0° direction exhibited the smallest average error using the fitting formula (Fig. 4b). The $H/D$ values in the 45 and 90° directions were slightly less consistent. This can be explained as follows: at 45 and 90°, the distance between the receptor and the isocenter increased; however, the increase in the distance between the receptor and MLC, which is the main producer of neutrons, was not synchronized with the former. This was particularly noticeable in the 90° direction. Therefore, when the subsequent analysis model predicted the $H/D$ value at different spatial positions, the corresponding fitting parameters were used in different directions to improve accuracy.

Effect of SOBP width on $H/D$ value

The 2-, 4-, 6-, 8-, 10-, and 12-cm SOBP widths and standardized simulation data based on the 6 cm result were selected. Figure 5 shows the variations in the relative $H/D$ values for different SOBP widths. As the width of the SOBP gradually increased, the relative $H/D$ value intensity at each location decreased almost linearly, and the drop range was not significantly different (approximately 30%). This phenomenon demonstrates that the ambient dose equivalent around the isocenter decreases with the thickening of the ridge filter when the MLC is fully closed. Although this finding is unexpected, it is understandable. The multileaf collimator is the main source of secondary neutrons downstream of the beamline. An increase in the thickness of the ridge filter reduced the energy of the beam, thereby reducing the neutrons produced by the interaction of the beam with the MLC.

Fig. 5Full-screen View

$H/D$ values as a function of SOBP width at isocenter θ = 0, 45, or 90° and r = 100 or 200 cm. The black line represented the binomial fit curve of the average $H/D$ values for all simulation locations (the beam energy is 400 $\text{MeV}/\text{u}$, MLC is closed)

The equation in this figure represents the coefficient $C_m$ proposed to indicate the influence of the SOBP width $C_{\text{m}}=0.000112{\text{m}}^2-0.0206487\text{m}+1.12662,\frac{{\chi }^2}{ndf}=0.000285/3.$

Effect of aperture size on $H/D$ value

Based on simulations using aperture sizes α from 0 $\times$ 0 ${\text{cm}}^2$ (closed) to 15 $\times$ 15 ${\text{cm}}^2$ in increments of 3 $\text{cm}$ per side, the effect of the open size of the aperture on the relative $H/D$ value is shown in Fig. 6.

Fig. 6Full-screen View

Ratio of neutron dose equivalent per therapeutic dose for an arbitrary aperture size (f) to that for a closed aperture, ${\left(\frac{H}{D}\right)}_a/{\left(\frac{H}{D}\right)}_0$, as a function of aperture size: (a) $H/D$ value in 0°, (b) $H/D$ value in 45°, and (c) $H/D$ value in 90° (beam energy of 400 $\text{MeV}/\text{u}$, SOBP width of 6 $\text{cm}$)

The secondary neutron dose equivalent is inversely proportional to the opening size of the MLC. The relative $H/D$ values decreased with increasing open aperture, and the minimum was below 5%. Neutrons generated by the beam-limiting device are the main source of secondary neutrons in radiotherapy. The trend of the simulated data in Fig. 6 agrees with the expectations, which also indicates that the results obtained using a closed MLC are conservative. When the opening size was expanded from fully closed to 15 $\times$ 15 ${\text{cm}}^2$, the neutron dose equivalent was reduced by approximately 90% in the 45 and 90° directions, respectively. At r = 200 cm and θ = 0°, the relative $H/D$ value drop was only approximately 60%, which was significantly smaller than that at other locations. Based on the observations in Sect. 3.1, secondary neutrons consist of two parts: high-energy direct neutrons and low-energy evaporation neutrons. The adopted cut-off energy of 15 $\text{MeV}$ subdivided the secondary neutron dose equivalent into a low-energy neutron dose equivalent, $H_1$, and a high-energy neutron dose equivalent, $H_2$.

As shown in Fig. 7, the decreasing laws of the $H/D$ values at different locations are distinct and are related to the neutron energy. At remote locations in the 0° direction, such as at 150 and 200 $cm$, the relative $H_1/D$ intensity first increased and then decreased. When the aperture size was at its maximum, the relative $H_1/D$ intensity was approximately 1 or greater than 1 (Fig. 7a, b). This phenomenon can be understood from the initial generation and eventual fate of the neutrons. As the aperture size increased, more carbon ions passed through the nozzle without creating neutrons in the MLC. In addition, a few secondary neutrons escaped from the upstream beamline. This difference can be attributed to two factors: a reduction in neutron production and an increase in neutron leakage. Because the leaking neutrons (mainly high-energy cascade neutrons) are forward-peaked, the high-energy neutron component at positions in the 0° direction increases. For locations further from the MLC, the change in relative value was more significant.

Fig. 7Full-screen View

(Color online) Ratio of neutron dose equivalent per therapeutic dose for an arbitrary aperture size (f) to that for a closed aperture, ${\left(\frac{H}{D}\right)}_a/{\left(\frac{H}{D}\right)}_0$, as a function of aperture size: (a) $H_1/D$ value in 0°, (b) $H_2/D$ value in 0°, (c) $H_1/D$ value in 45°, (d) $H_2/D$ value in 45°, (e) $H_1/D$ value in 90°, and (f) $H_2/D$ value in 90° (beam energy of 400 $\text{MeV}/\text{u}$, SOBP width of 6 $\text{cm}$).

Because there is no interference from escaped high-energy neutrons in the 45° direction, the falling laws of the $H_1/D$ and $H_2/D$ at different locations are similar (Fig. 7c, d). In the 90° direction, both the relative intensities of $H_1/D$ and $H_2/D$ significantly decreased with an increase in the aperture size. However, the falling laws of $H_1/D$ and $H_2/D$ at the isocenter and 25 and 50 cm from the isocenter were similar to those in the 45° direction. At 100, 150, and 200 cm from the isocenter, the $H_1/D$ and $H_2/D$ gradually decreased with increasing aperture size (Fig. 7e, f). This phenomenon is considered to be related to the relative positions of the receptors and MLC. The locations of the isocenter and 25 and 50 cm from the isocenter in the 90° direction make a small angle with the MLC. Therefore, the relative intensities of the $H_1/D$ and $H_2/D$ falling laws are the same as those of the small-angle area. The locations 100, 150, and 200 cm from the isocenter relative to the MLC were large angles (above 45°), and the falling laws followed a new pattern. Most neutrons arriving at these locations are side-scattered from the MLC. The thicknesses of the shields for the side-scattering neutrons were different. The closer the secondary neutrons are generated to the side surface of the MLC, the thinner the shield required to penetrate and the greater the dose contribution to the receptor. This also explains why the slope of the relative H/D intensity curve increases with aperture size.

Coefficient $C_a$ was proposed to indicate the influence of the MLC aperture size. Considering that the falling laws of $H/D$ in different directions were not completely consistent, we fitted the change in H/D values in different directions separately (Fig. 8). The fitting parameters and goodness-of-fit values are presented in Table 4.

Fig. 8Full-screen View

(a) 0° (beam direction), (b) 45°, (c) and 90° (beam transverse direction) relative $H/D$ value fitting curves. The linear formula was used to fit the change of the $H/D$ value in each direction

Validation of the analytical model

The predicted values ${\left(\frac{H}{D}\right)}_{\text{pred}}$ at various locations under different treatment settings were calculated according to the above equation and compared with the simulated values ${\left(\frac{H}{D}\right)}_{\text{mc}}$ (Fig. 8). The difference between the predicted and simulated values was very small in the 0° direction, and the ratio of ${\left(\frac{H}{D}\right)}_{\text{pred}}$ to ${\left(\frac{H}{D}\right)}_{\text{mc}}$ was approximately 1 (Fig. 8a). However, the ratio between the 45° and 90° directions gradually increased with distance. The reason for this was the systematic error caused by defects in the analytical model at angles other than 0°, as described in Sect. 3.2. A comparison of Fig. 8a, b show that the influence of the SOBP width on the prediction accuracy was insignificant, whereas the error caused by the aperture size was relatively large. Figure 8c shows that the error of the analysis model was approximately 20 % within the classical field size (10 $\times$ 10 ${\text{cm}}^2$) in the 0° direction, which was within the typical neutron measurement uncertainty. The results illustrate the feasibility of the proposed analytical model. If the analytical model can be optimized or the predicted data corrected, the adaptability and accuracy of the evaluation model can be further improved.

Out-of-field neutron dose performance evaluation

The $H/D$ values (perpendicular to the beam direction) from the field were compared with those reported in the literature (Fig. 9). The 300 $\text{MeV}/\text{u}$ beam energy and 6 cm SOBP widths were selected to be comparable to those found in published articles. Table 5 lists the beamline parameters used in previous studies. The $H/D$ value is related to the specific model and design of the equipment and is affected by many factors, such as the beam energy, aperture size, air gap, phantom, location, and method of measurement or calculation. Therefore, these differences must be considered when comparing $H/D$ values of various facilities (Fig. 10).

Fig. 9Full-screen View

(Color online) Ratio of the predicted values, ${\left(\frac{H}{D}\right)}_{\text{pred}}$, to the Monte Carlo simulated values, ${\left(\frac{H}{D}\right)}_{\text{mc}}$. (a) different SOBP width in 0, 45, and 90° directions, (b) different MLC aperture sizes in 0, 45, and 90° directions, (c) different MLC aperture sizes in the 0° direction.

Fig. 10Full-screen View

(Color online) Relationship between the neutron dose equivalent per therapeutic dose and the distance from the edge of the field using different radiotherapy equipment (The upper and lower bounds describe the empirical range of doses associated with different intensity modulated radiation therapy techniques ^[30])

The neutron ambient dose equivalents in carbon and proton radiotherapy are comparable in the horizontal direction of the out-of-field area. A popular explanation for this phenomenon is that the cross section for nuclear interactions is higher for carbon-ion beams than for proton beams, but fewer carbon ions are required to deliver the same therapeutic dose: approximately ${10}^{10}$ carbon ions per treatment, which is a factor of 100 fewer ions than proton therapy ^[31-32].

Compared to the experimental measurement data of the carbon-ion radiotherapy device of the Hyogo Ion Beam Medical Center (HIBMC) by Yonai et al., the simulation results in this study were 20 times greater ^[11]. This discrepancy may originate from the difference between the beamline components and the calculation method. (1) The collimator system in this study is composed of a fixed primary collimator and a 6.48-cm thick tungsten copper alloy MLC. In reference 10, an adjustable four-leaf collimator was used for collimation in advance to reduce the beam profile to the required field size, and a 140-mm thick iron MLC was used, which could significantly reduce the neutron dose. (2) Reference 10 used the WENDI-II wide-energy neutron dose equivalent monitor for measurement, whereas this study used the Monte Carlo method for calculation. The Monte Carlo method results in a calculated value that is larger than the measured value ^[33].