FLUKA Monte Carlo simulation for the Leksell Gamma Knife Perfexion radiosurgery system: Homogeneous media

Gamma Knife radiosurgery is worldwide recognized as one of the most relevant techniques in order to perform radiation treatments for brain tumours, arteriovenous malformations and brain dysfunctions []. The Leksell Gamma Knife Perfexion stereotactic radiosurgery system, introduced in 2006 by Elekta Instruments AB Stockholm Sweden, offers a non-invasive alternative for many patients for whom traditional brain surgery is not an option. The LGK-PFX is a ⁶⁰Co based medical device, able to deliver targeted radiation on a specific area or areas of the brain. The Perfexion model has a completely redesigned radiation unit geometry with respect to its predecessors. In previous Gamma Knife models the beam collimation system consisted in two separate parts: internal collimators which are arranged in a hemispherical pattern and embedded within the shield, and interchangeable collimator elements located in hemispherical external helmets. When the helmet is in the treatment position the entire collimator system forms 201 conical channels with a circular cross-section. The size and the shape of the dose deposition is determined by the collimator elements located in the helmet, which are available with four different hole diameters of 4, 8, 14 and 18 mm respectively. Further details of previous Gamma Knife geometry models can be found in references [, , ]. The most striking difference between Perfexion and its predecessors is the absence of the external part of the collimation system. The LGK-PFX is equipped with three different collimator sizes, able to focus the gamma rays on a common spot, the isocenter of the field, with a radial dimension of 4 mm, 8 mm and 16 mm respectively. Gamma radiation is provided by 192 ⁶⁰Co sources arranged on 8 independent sectors of 24 elements. In order to change the beam diameters, each sector can be placed in correspondence of the selected collimator set. This action is performed by servo-controlled motors located at the rear of the radiation unit []. Each Gamma Knife is accompanied by a computer based treatment planning system: the Leksell GammaPlan, which is specifically designed to simulate and plan the dose delivered by the device. The entire treatment procedure, including the surgical patient head fixation in the stereotactic frame, the image acquisition, the planning and the delivery of the radiation, is typically performed in a single session. In order to obtain the calculated dose distributions with the required accuracy, in clinical time, LGP employs simplified physical models. To perform the calculation, LGP version 8.2 uses a known beam profile and a pre-measured dose at a reference point and applies inverse square law and linear attenuation exponential formula. Furthermore, it takes into account photon scattering simply by including a measured off-axis ratio and considers the crossed media as a unit density homogeneous material. Although the healthy brain tissue can be assumed to be relatively homogeneous, when the radiation beam crosses inhomogeneous zone interfaces, the dose distribution can be affected by electronic disequilibrium. The verifications of the TPS output, as required by the quality assurance program, and the intrinsic difficulties of dose measurements, especially in small fields and/or in the presence of heterogeneous zones, make the Monte Carlo method a suitable and powerful tool in order to perform accurate comparisons of dosimetric parameters as well as dose distribution calculations.

In the literature several studies concerning the implementation of Monte Carlo simulation with the purpose of LGP comparison can be found, both in homogeneous and inhomogeneous conditions. Since the LGK-PFX is provided without detailed geometry drawings, and since the collimation system geometry measurements are substantially impracticable, all these studies are about previous Gamma Knife models (model C and model 4C). Many studies have been performed with EGS4 Monte Carlo code [, , ], in which the properties of Gamma Knife radiation fields in different conditions have been considered. Other studies have been performed with PENELOPE Monte Carlo code, using homogeneous as well as inhomogeneous phantoms [, , , ]. For the homogeneous problem, the authors observed correspondence within 1.5% between their calculation for the ROF and the data calculated by EGS4 studies as well as results provided by Elekta. For the heterogeneous problem, the authors observed that LGP underestimates the dose distribution near the tissue interface behind an air cavity []. In the study of Romano and Cuttone, the authors focused their attention on the simulation of relative dose distributions for the Gamma Knife model 4C, using Geant4 Monte Carlo code []. For homogeneous phantoms, they noticed that LGP performs a correct calculation of dose distributions, while, for inhomogeneous phantoms, the TPS underestimates the deposited dose near air cavities. MCNP MC code has been used, for homogeneous media, in order to perform the ROF and relative dose profile calculations, in the work of Trnka, Novotny and Kluson []. FLUKA is a multi-purpose Monte Carlo code which can transport a variety of particles over a wide energy range in complex geometries. The code is a joint project of INFN (National Institute for Nuclear Physics of Italy) and CERN [, ]. FLUKA has different applications in many fields of physics such as high energy experimental physics, engineering shielding, detector and telescope design, cosmic ray studies, dosimetry, medical physics and radiation therapy. We have chosen to test the FLUKA MC code as an alternative to EGS4, MCNP and PENELOPE, which have been widely used successfully for the Monte Carlo calculation of Gamma Knife parameters. It might be interesting therefore to have a qualitative comparison among these codes for photon beams such as from ⁶⁰Co, for which the code is perfectly suited []. Therefore the LGK-PFX dosimetry problem can be a good candidate in order to investigate the code capability on this energy range.

On the basis of the Perfexion design parameters, made available by Elekta Instruments AB Stockholm Sweden under a confidentiality agreement, an accurate geometrical model of the radiation unit has been developed.

The MC implemented model describes the entire radiation beam and all the 192 sources at the same time. We have taken into account the technical construction details of the sources, the collimation channels, and all the involved materials. Figure 1 shows the FLUKA geometrical model implemented in the simulation compared with the LGK-PFX radiation unit. Each radioactive source consists in ⁶⁰Co cylindrical pellets 1 mm in diameter and 1 mm in length encapsulated in metallic bushing. The capsules contain an appropriate number of pellets in order to ensure an average activity of about 1 TBq to each source. The ⁶⁰Co nuclei transforms itself through β⁻ decay to an excited state of ⁶⁰Ni, which in turn decays to its ground state emitting two photons with energy of 1.17 MeV and 1.33 MeV respectively. In our simulation, we have modelled each source body as a cylinder with a diameter of 1 mm, and a length of 18 mm, made of homogeneous cobalt material. The cylindrical capsules, in which the ⁶⁰Co pellets are housed, have also been simulated as made of homogeneous metallic material. In order to take into account the high number of sources involved in the simulation process, the initial kinematics in FLUKA has been programmed codifying an external user routine. Based on random number generators, the routine first selects a source among the 192 available, and then performs a uniform sampling of the beam starting point within the source volume. This process is reiterated for each simulated ⁶⁰Co nuclear decay. The β⁻ electron (average energy of about 315 keV) is supposed to be absorbed by the metallic source cover and/or by the source itself, therefore, for each sampled starting point, only the two photons have been produced and transported in sequence. With respect to the whole solid angle, only a small fraction of the photons emitted by the radioactive decay can pass through the collimator holes without being absorbed by the outer shield. Regarding the simulation process, in principle, the solid angle value doesn't affect the time consumption directly, but a statistics with errors ≤1% is reached with a number of primary particles which largely increases with the solid angle value. Therefore, implementing the isotropic sampling of the beam direction, which involves the loss of the large amount of the relevant events, would lead to a large increase of the primary photons and, indirectly, of the calculation time. In order to have a good balance between the quality of the results and the CPU time consumption (expressed through the number of primary photons) we decided to sample the beam direction in a restricted solid angle with a given semi-aperture of about 5°. This value has been chosen after a series of tests. We made several simulations with different solid angles, starting from Ω = 4π, comparing the calculated ROF values. The dose deposition for each collimator has been scored in a cylindrical volume element, length 1 mm and radius 1 mm, placed in the geometrical centre of the dosimetry phantom. For each implemented solid angle the calculated ROF values were within 1% of variation. The semi-aperture of 5° represents also the inferior limit for the acceptable solid angle values, in order to reproduce correctly the relative linear dose profiles. We noticed that, for lower values, the profiles are narrower than the Elekta reference ones. For values equal to or higher than 5° the contribution of photons scattered in the collimation channels to the total dose is properly taken into account. The transport of secondary charged particles (electrons and positrons) in the non-target regions, i.e. the collimation system, the source volumes and the capsules containing the sources, requires a large amount of calculation time. The FLUKA energy cut-offs have been differently set depending on the regions in which the gamma rays and the secondary particles interact. For the phantom, in which the dose calculation is performed, the cut-off energy thresholds for transport and production of photons and secondary electrons have been set to 10 keV. Instead, in the non-target materials, the cut-off threshold values have been set to 2 MeV. This value has been chosen in order to neglect the production and the transport of the secondary charged particles. This choice has been tested in order to ensure that no significant bias was introduced in the calculated dose deposition scored in the centre of the dosimetry phantom. This test has been performed as a comparison between scored dose values in a cylindrical volume element, length 1 mm and radius 1 mm, placed in the geometrical centre of the dosimetry phantom. We implemented several simulations, varying the thresholds for transport and production of secondary charged particles in the non-target materials. We noticed that, for this particular choice of scoring volume, the secondary electrons and positrons, generated by photon interactions in the collimation channels, do not have any relevant effect on the scored dose deposition, therefore we decided to neglect their production and transport. Eventually, a total of three different geometry set-up have been prepared, one for each collimator type. The MC calculation has been performed in the homogeneous condition of the phantom target, simulated as a sphere of 160 mm in diameter made of water equivalent PMMA Gammex-457 Solid Water material. The phantom has been exposed to the field generated by all the 192 sources on, with geometrical centre superimposed to the isocenter of the device. For each collimator size, we have investigated the relative dose distribution along the three coordinate axes passing through the isocenter, and the relative output factors. The FLUKA relative dose distributions have been obtained by binning the dosimetry phantom by means of VOXELs. In order to average the statistics fluctuations in the Monte Carlo calculations, the output data distributions of FLUKA calculation have been smoothed by means of a spline, importing the data in MATLAB. The FLUKA profiles have been then normalized with respect to the dose value scored in the VOXEL placed in the isocentre of the field, which has coordinate (0,0,0) and corresponds to the origin of the frame of reference.

The validation of the FLUKA model has been performed by comparing simulated results with experimental measurements for each collimator size, superimposing one-dimensional relative dose profiles. The experimental data have been acquired at the Gamma Knife centre of Niguarda Ca'Granda Milan Hospital. The dosimetric measurements have been performed using radiochromic EBT films, made by ISP (International Specialty Products). GafChromic EBT film consists of two active layers 17 μm thick separated by a surface layer 6 μm thick protected by two polyester layers each one 97 μm thick. The active layers are composed by water equivalent material with Z_eff = 6.98. EBT GafChromic films have been placed in the Elekta dosimetry phantom and exposed to the field generated by the 192 collimators open. For each film measurement, the exposure time has been set to 100.8 s, corresponding to an absorption dose of 5 Gy in the configuration with all the 16 mm collimators open. The proportionality between absorbed dose and exposition time has been obtained by means of a free-air ion chamber with a small sensible volume (31016-PTW Freiburg), inserted in the dosimetry phantom. For each collimator set two orientation for the film exposure have been selected, film lying in coronal plane (corresponding to the ZX-plane the in Elekta frame of reference) and film lying in axial plane (corresponding to the XY-plane the in Elekta frame of reference). A set of three measurements have been performed for each film orientation. When exposed to gamma rays the active layer responds to radiations forming a dark polymeric dye. Since the grey level is not linearly proportional to the value of the delivered dose, a calibration procedure is also required. For this reason, a total of 16 films have been exposed to a rising radiation dose from 0 to 10 Gy in order to obtain a calibration curve, necessary to correct this non linear response. After 48 h from exposure, EBT films were scanned using an Epson V700 scanner. GafChromic films have been digitalized with a 48 bit pixel depth, and a resolution of 400 dpi, as not compressed TIFF images. After the scanning process the images were analysed by means of the Image processing and analysis in Java (ImageJ) software. We found that a third degree polynomial calibration curve is suitable to fit the data. The experimental error is affected by many factors such as the difference of response from a film to another, the error in the determination of the net optical density, the time for development after exposure, the scanner parameters and the errors of the calibration curve. The experimental error was determined using the Devic method [], that, for a given dose over 0.3 Gy, estimates the uncertainty smaller than 5%. The validation of the simulation has been performed comparing the relative linear dose profiles for each available collimator size and for each coordinate axis. The acceptance criterion between measurements and FLUKA calculations is based on the gamma index (GI) method []. This index takes into account the relative shift both in terms of intensity and in terms of position, combining together the dose difference (DD) and the distance to agreement (DTA) methods. The DD is the difference between the dose value in a measured data point and the dose value in a point of the calculated distribution which has the same coordinates. The DTA is the distance between a measured data point and the nearest point in the calculated dose distribution that has the same dose value. By definition, when the gamma index results ≥1 the dose distributions are considered to be not in good agreement. GI test is instead considered passed if the calculated value results <1. In this work the acceptance threshold values for DD and DTA are derived from the Report of the American Association of Physicists in Medicine for stereotactic radiosurgery, respectively as DD = 3% and DTA = 2.4 mm []. After the validation process, the FLUKA relative dose profiles have been compared with LGP data, provided by Elekta Instruments AB and calculated by TPS in the same physical set-up of the phantom target. Also in this case the quality index is based on the gamma index method.

For all Gamma Knife models, the absolute output is calibrated with reference to the dose rate of the largest open collimator, i.e. the 16 mm in the case of the Perfexion model. The Relative Output Factors for the LGK-PFX are defined as the ratio between the dose given by a set of collimators and the dose given by the 16 mm collimators, which are taken as reference size, with ROF value set equal to 1.000. The FLUKA calculation for the ROF has been performed considering the ratio of the mean energy values, normalized to one source particle, deposited in a proper volume of the dosimetry phantom. For the purpose of ROF determination, in order to perform the correct calculation, the volume must be the same for each collimator set and must be uniformly covered by the radiation beam. For these reasons, and in order to take carefully into account the focal spot generated by the 4 mm collimators, the FLUKA scoring volume has been chosen as a cylinder of 1 mm in radius an 1 mm length placed in the isocenter of the field with longitudinal axis lying along the Z axis. Monte Carlo results have been compared with the data provided by Elekta Instruments AB, calculated by LGP in the homogeneous condition of the phantom target. The FLUKA calculations have been performed in the cluster made available at Medical Physics Department of Niguarda Ca'Granda Milan Hospital, made up of twelve quadric-processor nodes, using the version 2011 of the code. For each simulation a total number of 4 × 10⁹ primary histories have been performed for a total calculation time of about 20 h.

In this work, using the FLUKA Monte Carlo code, a complete simulation for the Leksell Gamma Knife Perfexion radiosurgery system has been successfully implemented. MC model is fully validated by means of comparisons with experimental measurements performed with EBT GafChromic films. FLUKA simulations and experimental measurements have been both performed in the simplest situation of 192 collimators opened and for a homogeneous water equivalent phantom. The quality index for the comparison is based on the gamma index method, which calculation results <1 for each analysed case. The FLUKA Monte Carlo model has then been used to calculate the LGK-PFX linear dose profiles, along coordinate axes, and the Relative Output Factors. The comparisons between FLUKA and LGP calculated linear profiles, for the 16 mm, 8 mm and 4 mm collimator size, are shown in Figure 2, Figure 3, Figure 4 respectively. The agreement between the two calculations is expressed by means of the gamma index method, which results <1 for each analysed case, according to the acceptance thresholds. The good agreement between FLUKA calculation for the ROF and the values provided by Elekta is reported in Table 1.

These results prove that the FLUKA code is capable to simulate successfully the LGK-PFX electromagnetic interactions. The accuracy of FLUKA results, especially in the case of the ROF determination, for which the difference between calculated and reference is under 1%, is fully comparable with the accuracy achieved by PENELOPE, EGS4 and MCNP calculation. For the latter codes the percentage differences between Elekta provided values and MC calculations are within 1.5% [, , ]. These results, however, cannot be directly compared to the previously reported in literature [, , , , , ], since they refer to different Gamma Knife models which have different radiation unit geometry.

The next step will be the study of the inhomogeneous condition of the target, performing new simulations after introducing heavy media or air cavities in the phantom. If the model will prove its goodness also in non-standard situation, it will be employed as a proper computational testing tool to be compared with LGP in real clinical situations. This work is in progress. In order to speed up the simulation time and increase the quality of the results, acquisitions of several phase space files of the current simulations are planned at the Medical Physics Department of Niguarda Ca' Granda Hospital. Eventually, importing a CT scan in the FLUKA simulation, it will be possible to analyse and perform various comparisons with the LGP software calculations for real clinical cases.

The authors would like to express their gratitude to Elekta Instruments AB Stockholm Sweden and, in particular to Dr. Jonas Johansson, for their help during the preliminary part of the work. We would like to express our thanks to Dr. Alessandro La Camera, person in charge at the Milan Niguarda Ca' Granda Hospital Gamma Knife centre, and his staff. We are grateful to the FLUKA collaboration for the support.