First steps towards a fast-neutron therapy planning program
© Garny et al; licensee BioMed Central Ltd. 2011
Received: 9 June 2011
Accepted: 25 November 2011
Published: 25 November 2011
Skip to main content
© Garny et al; licensee BioMed Central Ltd. 2011
Received: 9 June 2011
Accepted: 25 November 2011
Published: 25 November 2011
The Monte Carlo code GEANT4 was used to implement first steps towards a treatment planning program for fast-neutron therapy at the FRM II research reactor in Garching, Germany. Depth dose curves were calculated inside a water phantom using measured primary neutron and simulated primary photon spectra and compared with depth dose curves measured earlier. The calculations were performed with GEANT4 in two different ways, simulating a simple box geometry and splitting this box into millions of small voxels (this was done to validate the voxelisation procedure that was also used to voxelise the human body).
In both cases, the dose distributions were very similar to those measured in the water phantom, up to a depth of 30 cm. In order to model the situation of patients treated at the FRM II MEDAPP therapy beamline for salivary gland tumors, a human voxel phantom was implemented in GEANT4 and irradiated with the implemented MEDAPP neutron and photon spectra. The 3D dose distribution calculated inside the head of the phantom was similar to the depth dose curves in the water phantom, with some differences that are explained by differences in elementary composition. The lateral dose distribution was studied at various depths. The calculated cumulative dose volume histograms for the voxel phantom show the exposure of organs at risk surrounding the tumor.
In order to minimize the dose to healthy tissue, a conformal treatment is necessary. This can only be accomplished with the help of an advanced treatment planning system like the one developed here. Although all calculations were done for absorbed dose only, any biological dose weighting can be implemented easily, to take into account the increased radiobiological effectiveness of neutrons compared to photons.
At the former research reactor FRM at Garching, Germany, neutrons have been used for radiation therapy since 1985, and 715 patients with different types of tumors were treated until its shut-down in July 2000. Since June 2007, neutron irradiation treatments are performed at the new Medical Application facility (MEDAPP) at the Research Neutron Source Heinz Maier-Leibnitz (FRM II) [1, 2]. This new reactor was designed to provide a very high neutron fluence rate, primarily at thermal energies. In the moderator tank of the reactor, the low-energy neutrons are converted by two uranium plates to high-energy neutrons for medical (and other) applications at beamline SR 10. Filters are used to match the beam characteristics (i.e., neutron spectrum and the neutron-to-photon ratio) to those of the first facility at the former FRM. For this reason, medical knowledge acquired during the operation of the first reactor can be used for the treatment at the FRM II. As main improvements to the earlier facility, the new facility includes a 3-fold increase in total neutron fluence rate and up to a 6-fold increase in field size . This leads to a decrease in treatment time and to an improvement of treatment conditions and quality, which are now comparable to those at other neutron therapy facilities like for example at Seattle, Fermilab or iThemba [4–6].
The next step of improvement would be to introduce a computer-based treatment planning system rather than to use waterphantom depth-dose curves for dose estimation which are used at the moment in a certificated (CE) procedure. The physical basics for such a planning system have already been thoroughly explored both by using the Monte Carlo Code GEANT4 for simulations of neutron and photon transport [7, 8], and by using a Bonner sphere spectrometer  to measure the neutron energy spectrum at the patient position.
It is noted that in the past, different Monte Carlo algorithms have been applied for dose assessment in radiation therapy. For example, several different programs (including MCNP, MCNPX, EGS4, BEAMnrc VMC, PHITS, GEANT4) were used to simulate photons and electrons [10–16], protons [17–21], neutrons [22, 23] and dose from secondary neutrons during proton and 12C irradiation [24–27]. GEANT4 also offers the opportunity to calculate the 3D dose distribution of various particles including neutrons, photons and all secondary particles . In addition it allows to include any weighting function for calculation of biologically-weighted doses. This was shown e.g. by  to be important for future dose assessment of high LET-irradiation.
This work presents the results of a detailed simulation of water phantom depth-dose curves with GEANT4, based on the measured neutron energy spectrum at the patient treatment position, and its comparison with measurements . First steps towards a patient dose calculation were performed inside a voxel phantom  recommended by the International Commission on Radiation Protection (ICRP) for use in radiation protection (ICRP 103/Match 2007). It is shown that the implemented program allows calculation of 3D-distributions of neutron and photon absorbed dose.
GEANT4 version 4.8.2 ran on a SuSe-linux system. The neutron-data library G4NDL 3.10 was installed which is based on the ENDF/B-VI cross section evaluation . The required GEANT4-physics list was built including low energy processes. In particular, the S(α, β)-matrix was implemented to simulate thermal elastic scattering. A detailed description and testing of this physics list is given in [7, 33]. The importance of the S(α, β)-reaction was also discussed by Enger et al .
The deposited energy was calculated with a scorer (G4DoseDeposit), which was modified to determine also the statistical error of the dose. The G4SDParticleFilter was applied to score the absorbed dose deposited by photons, electrons, protons and neutrons separately. In order to calculate the neutron fluence rate, an energy-binned G4CellFlux scorer was used, which was also modified to calculate the statistical error. The binning was done with G4SDParticleWithEnergyFilter. Details on the usage of GEANT4 can be found in .
Specification of Perspex in GEANT4
atomic composition *
ratio of elements:
temperature = 300 K
state of aggregation = solid (kStateSolid)
density = 1.19 g/cm3
for energy loss purposes:
ChemicalFormula = "(C_2H_4)_N_Polyethylene"
The volume inside the solid water box, where the measurements were performed, was simulated as an array of tube-shaped chambers with a radius of 0.38 cm, a height of 1.21 cm and a volume of 0.54 cm3 that also consisted of water, in order to simulate an undisturbed measurement. Note that after the measurement, the doses measured inside the water phantom by means of the ionization chambers were corrected for the influence of the chamber on the measurement to provide the doses in an undisturbed phantom [3, 30]. In the simulation, the chambers were parameterized in such a way that the first chamber was placed at 0.5 cm depth in water, followed every centimeter by another detector chamber. This array of tube-shaped chambers corresponds well with the ionization chamber used for the real measurements (volume 0.5 cm3 ). In total there were 40 measurement volumes, which corresponds to a maximum depth of 39.5 cm. In the case of the voxelised phantom, the parameterized chambers were removed and the depth-dose scoring was done directly inside the voxels.
In the voxelised geometry, the total absorbed dose and the absorbed dose from primary and secondary photons and electrons were scored separately, simulating the real measurement, where both the neutron and photon doses were determined. In the case of the solid box geometry it was possible to collect more information, because less scoring volumes were present and therefore less data storage space was needed for each quantity. The total, neutron, proton, electron/positron and photon dose were all scored separately as well as the local moderated neutron fiuence rates in 58 energy bins. In spite of this, the calculation time for the solid box is about 3.5 times faster (4h compared to 14h for 106 primary neutrons on a 3.2 GHz processor pc).
The human voxel phantom used here is based on CT-scans of a real person. The Hounsfield numbers of the CT-slices were translated to 141 organ-IDs by hand, using anatomical atlases . Each voxel has the size 1.775 × 1.775 × 4.84 mm3 (x, y, z-direction), with the long side of the voxel aligned along the z-axis of the voxel phantom (increasing from toe to head). In total, there are 299 columns, 137 rows and 337 slices (x, y, z-direction) of voxels in the phantom, with vacuum-filled voxels surrounding the human body. The phantom represents an idealized female of 163 cm height and 60 kg body mass, based on the voxel phantom segmented from the CT-scan of a woman (Laura: 167 cm height, 59 kg). Compared to the original phantom Laura, the voxel size was changed, and the masses of individual organs were adjusted to reference values, to fulfill the requirements of ICRP 89 , using anatomical books for guidance. In this way, the Reference Female REGINA that was used in the present task was constructed.
Because the application of the FRM II neutron beam is most promising for the head and neck region, only the upper quarter of the voxel phantom (including 87 slices) was used for the calculations. The orientation of the voxel phantom in GEANT4 is the following: the column numbers (x-direction) increase from the right to the left side of the phantom, the row numbers (y-direction) increase from the nose to the back side of the head, and the slice numbers (z-direction) increase from the breast to the head. The voxel phantom was implemented in GEANT4 using the fast phantom parameterization (imported from version 4.9.0 into the used version 4.8.2), which includes all voxels (also the surrounding vacuum voxels) and identifies a specific voxel by its position in the grid. This algorithm is very fast but requires a lot of memory because the scored information is saved for all voxels that are placed.
The material of the voxel at (x, y, z)-position is set according to the number given in the phantom's datafile: the 141 organ IDs are projected onto 30 different materials (with arbitrary numbers) which are then used in the calculations (see  for details on the atomic composition and organ-to-material conversion).
In order to sample the neutrons and photons in the calculation according to the measured/calculated spectra, these spectra were integrated and normalized. From the resulting probability function, the primary neutrons and photons were sampled using random numbers.
The primary neutrons and photons were simulated to be homogeneously distributed over the whole beam. Furthermore, the beam profile was approximated to be rectangular with no decrease towards the edges. This is a simplification as the real beam is spread because of neutron scattering inside the beam tube and the large lateral dimensions of the source. Therefore, the real beam has a shallow decline towards the edges .
Figure 5 also demonstrates the excellent agreement of the calculations with the measurements. For example, the decrease of the dose rate with depth is very similar for both cases. Small differences between the measured and the calculated curves in terms of their absolute values could be - among other reasons - due to the different date of the measurement (the Bonner sphere measurements were performed approximately 1 year after the depth dose measurements, resulting in a small converter plate burn-up of approximately 2% ) and, for greater depths, due to calculation statistics. Furthermore, changes in the humidity of the air in the treatment room and the beam tube as well as the exact position of the control rods have an influence on the neutron spectrum and therefore can also change the deposited dose.
Finally, the comparison between the depth dose curves calculated with and without the concrete walls of the treatment room (Figure 5) demonstrate the influence of these walls to be of minor importance. For example, for small and medium depths (up to 10 cm) it is less than 2%. Thus, it is not necessary to consider the treatment room when calculating the dose inside the beam in the voxel phantom.
The primary FRM II treatment beam also includes a photon component, which is caused by prompt gammas during fission, by delayed gammas from fission products in the converter plate, by gammas from activated structures and by Compton scattered gammas from the reactor core. This primary photon spectrum is not yet characterized experimentally. However, an MCNP calculation has already been performed , transporting the photons through the beam tube up to the patient treatment position. With this program, it was not possible to simulate exactly the production of photons from the decay of activation and delayed fission products in the converter plates: The simulation resulted in an overall photon fluence rate which was 40% too low, and in a poor reliability of the photon spectrum. Nevertheless, this spectrum was used as a first estimate, to obtain information on the photon depth dose distribution inside the phantom.
For the calculations inside the water phantom, both the solid box and voxelised geometry were used (with and without surrounding walls), and the results compared to the measured data . It should be noted that the measured photon dose rate also includes contributions from secondary photons produced by the neutrons inside the water phantom, as the twin method applied in the measurements does not distinguish between primary and secondary photons.
The right panel of Figure 9 shows the dependence of the DVH of the two submandibular glands on the inclusion of the RBE. For this qualitative analysis, a clinical RBE of three (derived from old clinical data (Wagner, priv com.)) for the neutrons was included after the Monte Carlo calculation, to estimate the influence of a biological weighting for neutrons and their secondary particles. However, at the current state of investigation a proper RBE as a function of neutron energy is not available. The fraction of dose deposited by photons is larger for the left than for the right submandibular gland. Therefore, the ratio of the weighted dose compared to the absorbed dose is larger for the right submandibular gland than for the left one. This reduces the effective dose in the healthy submandibular gland and similarly in other deeper-located organs.
The absorbed depth dose distribution of the primary neutrons and photons of the FRM II beam was calculated using GEANT4 in a voxelised water phantom, and a very good agreement with measured data was found. The voxel-algorithm of GEANT4 was validated and influences of the surrounding walls on the dose in the main beam were studied.
The voxel phantom REGINA  was used as a first test and a 3D dose distribution of a salivary gland was calculated. The depth dose curve obtained was compared to that in the water phantom and the lateral distribution in different depths was discussed in detail. Furthermore, the dose volume histogram of this test case was calculated and the isodose representations in one slice given as an example. It is emphasized that the calculations described in the present work were done using a parallel neutron beam impinging on the human voxel phantom. Any allowance for beam divergence is expected to modify the calculated isodose lines and DVHs somewhat. Thus, before any more realistic 3D distributions can be given for an individual patient, the beam characteristics and its influence on the 3D dose distribution within a patient should be investigated in detail. In this respect, the given dose distributions (see Figures 7, 8, 9 and 10) should be interpreted still with care. Nevertheless, it is noted that the given dose distributions represent a major step forward in planning of fast-neutron therapy at FRM II, as they are based on quite realistic neutron and photon input spectra, and a realistic human morphology. Note that, to the best of our knowledge, GEANT4 has never been used for this application, and so far only simple water phantoms are being used in the certified dose estimate procedure applied at the MEDAPP facility.
GEANT4 allows to assess a biological dose using radiation weighting factors. Such factors can include particle type, energy and energy loss as well as dose deposition or other quantities like oxygen content or radiation sensitivity of the tissue involved. For the case of chromosomal aberrations Schmid et al.  have shown that the irradiation is more effective in the first few centimeters of tissue compared to deeper regions. Similar results were obtained by Magaddino et al.  for the tumor control probability and the RBE for permanent colony control. The effect was qualitatively shown in this paper by applying a clinical RBE and could also be seen when artificial RBEs were implemented in the Monte Carlo calculation [33, 41]. This helps to lessen the problem of high doses in healthy organs on the patient's side which is opposite to the beam. In the future, cell experiments like [42–44] will be used to improve the biological weighting of the dose.
It is concluded that the results presented here represent a significant step further towards a reliable, patient-specific treatment planning system at the fast-neutron therapy facility MEDAPP. However, further work needs to be done including modeling of beam divergence and exploring the influence of patient-specific parameters on organ doses, before more accurate patient-specific organ doses can be calculated. It is also planned in the future to apply the present approach to physical humanoid phantoms.
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.