Abstract
Background
The Monte Carlo code GEANT4 was used to implement first steps towards a treatment planning program for fastneutron 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).
Results
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.
Conclusions
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.
Background
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 shutdown in July 2000. Since June 2007, neutron irradiation treatments are performed at the new Medical Application facility (MEDAPP) at the Research Neutron Source Heinz MaierLeibnitz (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 lowenergy neutrons are converted by two uranium plates to highenergy neutrons for medical (and other) applications at beamline SR 10. Filters are used to match the beam characteristics (i.e., neutron spectrum and the neutrontophoton 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 3fold increase in total neutron fluence rate and up to a 6fold increase in field size [3]. 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 [46].
The next step of improvement would be to introduce a computerbased treatment planning system rather than to use waterphantom depthdose 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 [9] 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 [1016], protons [1721], neutrons [22,23] and dose from secondary neutrons during proton and ^{12}C irradiation [2427]. GEANT4 also offers the opportunity to calculate the 3D dose distribution of various particles including neutrons, photons and all secondary particles [28]. In addition it allows to include any weighting function for calculation of biologicallyweighted doses. This was shown e.g. by [29] to be important for future dose assessment of high LETirradiation.
This work presents the results of a detailed simulation of water phantom depthdose curves with GEANT4, based on the measured neutron energy spectrum at the patient treatment position, and its comparison with measurements [30]. First steps towards a patient dose calculation were performed inside a voxel phantom [31] 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 3Ddistributions of neutron and photon absorbed dose.
Material and methods
The Monte Carlo code GEANT4
GEANT4 version 4.8.2 ran on a SuSelinux system. The neutrondata library G4NDL 3.10 was installed which is based on the ENDF/BVI cross section evaluation [32]. The required GEANT4physics 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 [22].
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 energybinned 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 [34].
Geometry of the water phantom
In earlier measurements [1,30] the neutron and photon depthdose curves in a water phantom were determined by a pair of thimble chambers with different neutron sensitivity that were irradiated at the patient position (EXTRADINchamber one out of tissueequivalent (TE) material A150 with TE gas; sensitivity to neutrons/gamma = 48.5/51.5%; chamber two out of Mg/Ar, neutron sensitivity 2%, chamber volume: 0.5 cm^{3}). The phantom consisted of a box of Perspex filled with about 1701 of water with an entrance window for the horizontal beam and a device to position the ionization chambers. This water phantom was simulated here in two different ways (see Figure 1): as a solid box with small measurement chambers inside and in voxelised form consisting of over 10 millions of boxshaped volume elements ("voxels"). Both phantoms consisted of water and were put into a cubic container made of Perspex with dimensions 63.5 × 63.5 × 52 cm^{3 }and 2.0 cm wall thickness (like in the real experiment, material definition see Table1). For the voxelised phantom, the dimensions were slightly enlarged to 64 × 64 × 52 cm^{3 }to simplify the voxelisation algorithm, and voxel sizes of 0.2 × 0.2 × 0.5 cm^{3 }were used in accordance with those used typically in CTimaging for radiation therapy. The Perspex container also included the entrance window of 12 × 34 cm^{2 }size sealed with two aluminum plates of 1.5 mm thickness.
Figure 1. Geometry of the Water Phantom. The two geometry types used for the simulation of the water phantom: left: solid box phantom with the measurement chambers inside (which are also out of water; 30 of the 40 chambers used in the calculations are shown) with one incoming neutron track; right: xyplane of the voxelised water phantom with tracks from 100 primary neutrons of the FRM II beam. For visualization purposes, a 2 cm voxel size was used here, instead of the 2 mm used in the actual GEANT4 simulations.
Table 1. Specification of Perspex in GEANT4
The volume inside the solid water box, where the measurements were performed, was simulated as an array of tubeshaped chambers with a radius of 0.38 cm, a height of 1.21 cm and a volume of 0.54 cm^{3 }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 tubeshaped chambers corresponds well with the ionization chamber used for the real measurements (volume 0.5 cm^{3 }[3]). 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 depthdose 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 10^{6 }primary neutrons on a 3.2 GHz processor pc).
For the solid box geometry, the influence of the surrounding walls, floor and roof of the therapy room was also studied. For this, the patient treatment room was implemented in a schematic way, using G4concrete from the NISTmaterial database as wall material [35]. The walls were simulated to be 83 cm thick (50 cm in case of the beam exit wall) with a beam exit window of 20 × 30 cm^{2}. In Figure 2, the topview of the MEDAPP therapy room as simulated in GEANT4 is shown together with 10 neutron tracks using a neutron spectrum that had been measured by means of a Bonner sphere spectrometer at the patient position [9,33].
Figure 2. Geometry of the MEDAPP Room at the FRM II. Wall material: G4concrete; max. beam window: 20 × 30 cm^{2}; distance to floor/ceiling: 145 cm; thin lines: 10 neutron tracks; thick gray line: neutron therapy beam.
In both geometries, the primary neutron and photon beam cross section was a square of 9 × 9 cm^{2}, which is the same as the one used for the measurements [1,3].
Geometry of the human voxel phantom
The human voxel phantom used here is based on CTscans of a real person. The Hounsfield numbers of the CTslices were translated to 141 organIDs by hand, using anatomical atlases [31]. Each voxel has the size 1.775 × 1.775 × 4.84 mm^{3 }(x, y, zdirection), with the long side of the voxel aligned along the zaxis of the voxel phantom (increasing from toe to head). In total, there are 299 columns, 137 rows and 337 slices (x, y, zdirection) of voxels in the phantom, with vacuumfilled 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 CTscan 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 [36], 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 (xdirection) increase from the right to the left side of the phantom, the row numbers (ydirection) increase from the nose to the back side of the head, and the slice numbers (zdirection) 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 [31] for details on the atomic composition and organtomaterial conversion).
The chosen test case was a salivary gland treatment of the right submandibular gland (lower jaw salivary gland). Because the real field size applied to the patient was almost a rectangle (see left side of Figure 3, a rectangular beam with 6 cm × 7 cm cross section was used in the simulations hitting the patient from the right side with a beam direction along the positive xaxis. In contrast to the real situation, the beam profile was also simulated to be rectangular (see next section ). The position of the beam was simulated in resemblance to that of the real patient case using the skeletal structure as a guideline. It should be emphasized, however, that some approximations had to be made because the headtobody angle of the voxel phantom was somewhat different to that of the real case (see Figure 3). In this way, the calculation algorithm could be tested and a first assessment of the absorbed dose distribution was possible. It should be noted that the field shape can easily be changed in the simulation to adapt it to the shape of the planning target volume (PTV) if necessary.
Figure 3. Beam Position and Size. Comparison of beam position and size for a real patient treatment (left panel; Loeper, MRI/TUM, priv. com.) and the simulated field in the human voxel phantom (gray lines in right panel).
Primary neutron and photon spectrum
The absorbed dose rate (in Gy/s) was determined in GEANT4 by multiplication of the calculated absorbed dose deposited inside the chamber or voxel per primary neutron or photon fluence (resulting in a value with the unit Gycm^{2}) with a total primary neutron fluence rate of 3.2 · 10^{8 }cm^{2}s^{1 }and a total primary photon fluence rate of 2.9 · 10^{8 }cm^{2}s^{1}. The primary neutron fluence rate had independently been determined before by goldfoil activation measurements in a water phantom [37] and by Bonner sphere measurements rescaled to the patient treatment position in the MEDAPP therapy room [9], while the primary photon fluence rate was derived from the comparison of calculated and measured depthdose curves in water (see discussion below). To get the absolute absorbed dose for a specific treatment field inside the voxel phantom, the resulting values were finally multiplied by the actual irradiation time. In Figure 4, the neutron (measured by a BSS spectrometer [9]) and photon (calculated with MCNP [38]) primary spectra are normalized to a total fluence rate of 3.2 · 10^{8 }cm^{2}s^{1 }for neutrons and 2.9 · 10^{8 }cm^{2}s^{1 }for photons, respectively.
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 [37].
Results and Discussion
Calculated depthdose curves in the water phantom
Neutrons
In Figure 5, the calculated absorbed dose rates in the solidbox and the voxelised geometry are shown together with the measured depth dose curve. The results of the solid box and the voxelised geometry are in excellent agreement. The latter were integrated over 4 voxels in lateral and 2 voxels in vertical direction around the central beam axis, resulting in a dose collection volume of 0.16 cm^{3 }from 8 voxels, to get better statistics. Figure 5 demonstrates that the voxelisation algorithm produces the same results as the solidbox calculation.
Figure 5. Neutron Dose Rate. Neutron dose rate in the water phantom (beam profile of 9 × 9 cm^{2 }): comparison of the voxelised (green triangle, integration over 8 central voxels) and the solid box geometry (blue circle = with walls, pink circle = without walls) calculated with GEANT4 (error bars provide one standard deviation as calculated in GEANT4), together with the measurement of Kampfer [30] (red square).
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 burnup of approximately 2% [37]) 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.
Photons
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 [38], 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 [1]. 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.
As for the neutron calculations, the photon calculations in the voxelised and solidbox geometry produce very similar results. Furthermore and similar to the neutron case, inclusion of the surrounding walls does not influence the calculated dose rates. In Figure 6, both the dose rate from primary and secondary photons which were produced by interaction of primary neutrons with the phantom material, are shown. Because the depthdose curve of the neutrons is consistent with the measured data both in relative and absolute terms (see previous section), the dose rate calculated for secondary photons produced by neutrons is also correct on an absolute scale. On average, the contribution from the secondary photons amounts to about 1015% of the total measured photon dose rate in the first centimeters. The total primary photon fluence rate used in the simulations (i.e., 1.8 · 10^{8 }cm^{2}s^{1}, calculated by integrating the photon spectrum) is not correct on an absolute scale (see discussion above) but had to be increased until the sum of the calculated depth dose curves of the primary and secondary photons matched the measurements in the first 10 cm, and a best estimate of the primary photon fluence rate of 2.9 · 10^{8 }cm^{2}s^{1 }was obtained which agrees well with the expected total photon fluence rate (see above). This means that the incident primary total fluence rates of photons and neutrons are rather similar. Though the primary photon spectra is expected to be softer, the shapes of the depth dose curves do not deviate much in the first 10 cm (see Figure 6). Therefore, it is concluded that the primary photon spectrum derived here can be used together with the renormalized total fluence rate for first calculations in the voxel phantom.
Figure 6. Primary and Secondary Photon Dose Rate. Primary and secondary photon dose rate calculated in the solidbox water phantom (error bars provide one standard deviation as calculated in GEANT4), in comparison to measured [30] data (beam profile of 9 × 9 cm^{2}); □ = total dose rate from primary (Δ) and secondary (○) photons; ■ = total dose rate from the renormalized dose rate of primary photons together with that of the secondary photons. Note that while the absolute uncertainties seem to decrease with depth, the relative uncertainties increase with depth.
Calculated doses in the voxel phantom REGINA
Depthdose curves and lateral distribution
In Figure 7, the calculated depthdose curve along the central beam is compared to that obtained in the water phantom. The decrease of dose with depth is similar for both phantoms, but slightly steeper for the voxel phantom, because its atomic composition is different from water (for example, the hydrogen content in adipose tissue (material number 28, see section about the geometry of the voxel phantom ) is 62.5% compared to 66.7% in water). The change of material with depth in the REGINA phantom is also indicated in the figure. When the material changes, e.g., from 19 (skin) to 28 (adipose tissue) or to 21 (muscle), the atomic composition changes and therefore a change in the dose can also be seen. In Figure 8, the lateral absorbed dose distribution is shown for three different columns inside the REGINA voxel phantom: x = 130, 150 and 170, which is equivalent to 3.55 cm distance in between. The neutrons are significantly scattered inside the REGINA voxel phantom and deposit their energy not only inside the beam profile but also in the surrounding tissue. This can be seen at a row (y)value of 40, where the sharp edges of the simulated rectangular beam can still be seen for small depths (x = 130 is equal to 0.3 cm depth at this yvalue) inside the phantom. These sharp edges are washed out with depth, so that in greater depths (e.g. x = 170 or 7.4 cm), the lateral dose is spread over several centimeters (the relative decrease in absorbed dose between y = 39 and y = 41 compared to the dose at y = 41 is 93% for x = 130, 73% for x = 150 and 64% for x = 170). At a yvalue of about 74, the neutrons have already passed through a certain thickness of tissue (i.e., x = 130 is equal to 2.0 cm depth at this yvalue) and the changes of the lateral spread at the same xvalues is less prominent (relative decrease in absorbed dose between y = 73 and y = 75 is 83%/68%/64% for x = 130/150/170, respectively). Note that the beam spread is expected to be even more pronounced for the real beam, because the beam divergence of 12° as well as the penumbra of the different filters and collimators were not included in the calculations, but would further increase the lateral spread of the beam (see the discussion in the section about the neutron and photon spectra). Furthermore, the influence of the material is visible in Figure 8 as well. At the same xvalue, the dose declines towards higher yvalues, which correspond to the back side of the phantom. This is caused by the phantom's uneven surface. Looking at the relevant slice (z = 44, see inset), it can be seen that the neck starts to bulge in the relevant area. Therefore, the radiation is absorbed before reaching the relevant depth whereas at lower yvalues (4055), less material is present. For the green data points (x = 150) an area of high statistical uncertainty can be seen between yvalues of 60 to 66. This is the area where the trachea is located, where fewer particles interact with the air inside and deposit dose there. Behind this area, more energy is therefore deposited which can indeed be seen in the data of (x = 170). Another influence of material can for example be seen at y = 35, 52, 69, 77 in depth x = 150, where the deposited dose is much smaller than in the surrounding voxels. The material at this position is hard bone (cortical), with a hydrogen content of only 3.5% (fractionmass) compared to more than 8 for all other soft matter (also bone marrow and spongiosa). Because the neutron dose deposition mainly depends on the hydrogen content of the present material, the dose in these voxels is about half of the dose in other voxels.
Figure 7. Absorbed DepthDose Curves. Absorbed depthdose curves for the FRM II neutron beam in the REGINA phantom, at slice number 45 and row number 48 (beam center; ■), assuming a rectangular beam profile of 6 × 7 cm^{2}, and in the voxelised water phantom (○); All curves were calculated for 3 min irradiation with a primary neutron fluence rate of 3.2 · 10^{8}n/cm^{2}s (no primary photons); for comparison, material numbers in the REGINA phantom along the central beam are included on the right yaxis (19: skin, 21: muscle, 28: adipose tissue).
Figure 8. Absorbed Lateral Dose Distribution. Absorbed lateral dose distribution of the 6 × 7 cm^{2 }FRM II neutron beam in different depths inside the REGINA voxel phantom after 3 minutes of irradiation with a total primary neutron fiuence rate of 3.2 · 10^{8}n/cm^{2}s (no primary photons) in slice (z) number 44 assuming a rectangular lateral beam profile.
Dose volume histograms
In patient treatment planning, a tool used for plan quality assessment is the dosevolumehistogram (DVH). In this histogram, the detailed spacial data is condensed into a plot with the fraction of volume of an organ irradiated with a certain dose on the yaxis and the corresponding absorbed dose on the xaxis. In the left panel of Figure 9, the cumulative DVH of the studied salivary gland case is shown for six organs: the healthy salivary glands of the REGINA phantom's left side, the treated right submandibular gland, and the (healthy) right sublingual and parotid gland which are also partly in the beam. It should be emphasized that in a real case, the neutron treatment is given mostly as a boost after a normal photon treatment at a linear accelerator. If the doses from the photon treatment were included, the DVH would change significantly. Furthermore, only the physical absorbed dose was considered here without biological weighting of the neutron dose. If such a biological weighting was included, the DVHs would also change because of the different fraction of dose deposited by photons, depending on the organ's position inside the body.
Figure 9. DoseVolumeHistogram. Cumulative (left) dosevolumehistogram (DVH) of the salivary glands after 3 min irradiation with the 6 × 7 cm^{2 }FRM II neutron beam, a total primary neutron fluence rate of 3.2 · 10^{8 }n/cm^{2}s and a total primary photon fluence rate of 2.9 · 10^{8 }7/cm^{2}s; r = right, 1 = left; subman gl = submandibular gland; subli gl = sublingual gland; par gl = parotid gland. Right histogram: dependence of DVH of right and left submandibular glands on RBE, dose normalized to dose at 50% volume of the right/left submandibular gland respectively.
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 deeperlocated organs.
Isodoses
In the top left panel of Figure 10, one slice of the voxel phantom is shown together with the isodose lines of the irradiation with the FRM II neutron spectrum. The isodose lines were achieved by a spline interpolation of the calculated dosematrix. In hard bone material, lower doses can be observed which may be due to the lower hydrogen content (31.2% element abundance compared to 62.5% in adipose tissue). Inside air cavities such as the trachea, the deposited dose is also somewhat reduced, due to the low density and the lack of hydrogen. This can be seen in the middle part of the slice. In the top right panel of Figure 10, the primary photon absorbed dose with depth is shown, while the lower panel shows the total absorbed dose from primary neutrons and photons (here neutrons are not weighted for their increased biological effectiveness compared to that of photons). It can be seen that the total absorbed dose is about 0.5 Gy higher with the primary photons than without. This is particularly important for the healthy left salivary gland, because the depthdose curve of the photons is very shallow. This effect leads to an increased absorbed dose to healthy tissues behind the tumor. Introducing the real distribution of neutrons and photons (including the beam spread) would not change this conclusion significantly. It should be noted here that the biological effectiveness of neutrons is very important to assess the dose in a real clinical case.
Figure 10. Absorbed Dose inside REGINA Voxel Phantom. Absorbed dose inside REGINA voxel phantom after 3 min irradiation with the 6 × 7 cm^{2 }FRM II beam, a total primary neutron fluence rate of 3.2 · 10^{8}n/cm^{2}s and a total primary photon fluence rate of 2.9 · 10^{8}n/cm^{2}s; all contour lines relative to the prescribed dose of 2 Gy (in %); slice 45; in red/blue, the right/left salivary glands are highlighted; top left: dose from primary neutrons; top right: dose from primary photons; bottom: total dose; No biological weighting for the neutrons was applied in the figure.
Conclusions
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 voxelalgorithm of GEANT4 was validated and influences of the surrounding walls on the dose in the main beam were studied.
The voxel phantom REGINA [31] 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 fastneutron 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. [39] 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. [40] 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 [4244] 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, patientspecific treatment planning system at the fastneutron therapy facility MEDAPP. However, further work needs to be done including modeling of beam divergence and exploring the influence of patientspecific parameters on organ doses, before more accurate patientspecific organ doses can be calculated. It is also planned in the future to apply the present approach to physical humanoid phantoms.
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
SG: Writing code, running the calculations, producing the outcome, writing Manuscript, WR: Discussing outcome at various stages of the project and critical revision of the manuscript, MZ: Constructing and providing voxel phantom, FMW: providing details in photon and neutron FRM II spectrum, helping with FRMII questions and revision of the manuscript, HGP: Development of the concept, discussion of the outcome and critical revision of the manuscript. All authors read and approved the final manuscript.
References

Wagner F, Kneschaurek P, Kastenmüller A, LoeperKabasakal B, Kampfer S, Breitkreutz H, Waschkowski W, Molls M, Petry W: The munich fission neutron therapy facility MEDAPP at the research reactor FRM II.

LoeperKabasakal B, Posch A, Auberger T, Wagner F, Kampfer S, Kneschaurek P, Petry W, Lukas P, Molls M: Fission neutron therapy at FRM II: Indications and first results.

Wagner F, Bücherl T, Kampfer S, Kastenmüller A, Waschkowski W: Thermal Neutron Converter for Irradiations with Fission Neutrons.

Risler R, Popescu A: Dosimetry measurements at the fast neutron therapy facility in Seattle.

Kroc T: Preliminary investigations of Monte Carlo simulations of neutron energy and let spectra for fast neutron therapy facilities.

Swanepoel M: The role of the ^{14N}(n, p)^{14}C reaction in neutron irradiation of soft tissues.

Garny S, Leuthold G, Mares V, Paretzke H, Rühm W: GEANT4 transport calculations for neutrons and photons below 15 MeV.

Garny S, Mares V, Rühm W: Response functions of a Bonner sphere spectrometer calculated with GEANT4.

Garny S, Mares V, Roos H, Rühm W, Wagner F: Measurement of the neutron spectrum and neutron dose at the FRM II therapy beamline with Bonner spheres.

Reynaert N, van der Marck S, Schaart D, der Zee WV, VlietVroegindeweij CV, Tomsej M, Jansen J, Heijmen B, Coghe M, Wagter CD: Monte Carlo treatment planning for photon and electron beams.

Bednarz B, Xu X: A feasibility study to calculate unshielded fetal doses to pregnant patients in 6MV photon treatments using Monte Carlo methods and anatomically realistic phantoms.

Jang SY, Liu HH: Underestimation of LowDose Radiation in Treatment Planning of IntensityModulated Radiotherapy.

Jeraj R, Keall P: Monte Carlobased inverse treatment planning.

Laub W, Alber M, Birkner M, Nüsslin F: Monte Carlo dose computation for IMRT optimization.

Sikora M, Muzik J, Söhn M, Weinmann M, Alber M: Monte Carlo vs. Pencil Beam based optimization of stereotactic lung IMRT.

Taschereau R, Roy R, Pouloit J: A comparison of methods to calculate biological effectiveness (RBE) from Monte Carlo simulations.

Fippel M, Soukup M: A Monte Carlo dose calculation algorithm for proton therapy.

Koch N, Newhauser WD, Titt U, D G, Coombes K, Starkschall G: Monte Carlo calculations and measurements of absorbed dose per monitor unit for the treatment of uveal melanoma with proton therapy.

Newhauser W, Fontenot J, Zheng Y, Taddei P, Mirkovic D, Titt U, Zhu X, Sahoo N, Schaffner B, Langenegger A, Koch N, Zhang X, Mohan R: SUFFT25: A MonteCarlo based dose engine for proton radiotherapy treatment planning.

Jiang H, Paganetti H: Adaptation of GEANT4 to Monte Carlo dose calculations based on CT data.

Paganetti H, Jiang H, Lee SY, Kooy HM: Accurate Monte Carlo simulations for nozzle design, commissioningand quality assurance for a proton radiation therapy facility.

Enger SA, af Rosenschöld PM, Rezaei A, Lundqvist H: Monte Carlo calculations of thermal neutron capture in gadolinium: A comparison of GEANT4 and MCNP with measurements.

Pignol J, Slabbert J, Binns P: Monte Carlo simulation of fast neutron spectra: Mean lineal energy estimation with an effectiveness function and correlation to RBE.

Jiang H, Wang B, Xu XG, Suit HD, Paganetti H: Simulation of organspecific patient effective dose due to secondary neutrons in proton radiation treatment.

Taddei P, Mirkovic D, Fontenot J, Giebeler A, Zheng Y, Kornguth D, Mohan R, Newhauser W: Stray radiation dose and second cancer risk for a pediatric patient receiving craniospinal irradiation with proton beams.

Yonai S, Matsufuji N, Kanai T: Monte Carlo study on secondary neutrons in passive carbonion radiotherapy: Identification of the main source and reduction in the secondary neutron dose.

Zheng Y, Fontenot J, Taddei P, Mirkovic D, Newhauser W: Monte Carlo simulations of neutron spectral uence, radiation weighting factor and ambient dose equivalent for a passively scattered proton therapy unit.

Carrier JF, Archambault L, Beaulieu L: Validation of GEANT4, an objectoriented Monte Carlo toolkit, for simulations in medical physics.

Karger CP, Jäkel O: Current Status and New Developments in Ion Therapy.

Kampfer S, Wagner F, Loeper B, Kneschaurek P: Erste dosimetrische Ergebnisse an der neuen Neutronentherapieanlage am FRM IIb. In 37. Jahrestagung der Deutschen Gesellschaft für Medizinische Physik. L. Bogner, B. Dobler; 2006:318319.

Schlattl H, Zankl M, PetoussiHenßN : Organ dose conversion coefficients for voxel models of the reference male and female from idealized photon exposures.

GEANT4: [Http://geant4.web.cern.ch/geant4/support/userdocuments.shtml] webcite

Garny S: Development of a biophysical treatment planning system for the FRM II neutron therapy beamline. In Phdthesis. Technische Universitat München; 2009.

GEANT4: [Http://geant4.web.cern.ch/geant4/support/userdocuments.shtml] webcite

Standards NI, Technology: NIST Database. [http://physics.nist.gov/cgibin/Star/compos.pl] webcite

ICRP: Basic anatomical and physiological data for use in radiological protection: reference values. Report 89, International Commission on Radiological Protection; 2002.

Breitkreutz H, Wagner F, Röhrmoser A, Petry W: Characterisation of the fast reactor neutron beam MEDAPP at FRM II.

Breitkreutz H: Spektrale Charakterisierung des Therapiestrahls am FRM II. In Diplomarbeit. Technische Universitat München; 2007.

Schmid E, Schraube H, Bauchinger M: Chromosome aberration frequencies in human lymphocytes irradiated in a phantom by mixed beam of fission neutrons and γrays.

Magaddino V, Wagner F, Kummermehr J: Preclinical screening of the biological effectiveness of a therapeutic neutron beam at SR10, FRM II.
Experimetelle Strahlentherapie und Klinische Strahlenbiologie 2007, 16:129131.

Garny S, Rühm W, Wagner FM, Paretzke HG: Neutron Therapy at the FRMII  Calculation of Dose inside a Voxel Phantom.
11th World Congress on Medical Physics and Biomedical Engineering 2009.

Oya N, Zölzer F, Werner F, Streffer C: Similar Extent of Apoptosis Induction at Doses of XRays and Neutrons Isoeffective for Cell Inactivation.

Slabbert J, August L, Vral A, Symons J: The relative biological effectiveness of a high energy neutron beam for micronuclei induction in Tlymphocytes of different individuals.

Vandersickel V, Mancini M, Slabbert J, Marras E, Thierens H, Perletti G, Vral A: The radiosensitizing effect of Ku70/80 knockdown Research in MCF10A cells irradiated with Xrays and p(66)+Be(40) neutrons.