Development and validation of an open source Monte Carlo dosimetry model for wide‐beam CT scanners using Fluka

Abstract Purpose Development and validation of an open source Fluka‐based Monte Carlo source model for diagnostic patient dose calculations. Methods A framework to simulate a computed tomography (CT) scanner using Fluka Monte Carlo particle transport code was developed. The General Electric (GE) Revolution scanner with the large body filter and 120 kV tube potential was characterized using measurements. The model was validated on benchmark CT test problems and on dose measurements in computed tomography dose index (CTDI) and anthropomorphic phantoms. Axial and helical operation modes with provision for tube current modulation (TCM) were implemented. The particle simulations in Fluka were accelerated by executing them on a high‐performance computing cluster. Results The simulation results agreed to better than an average of 4% of the reference simulation results from the AAPM Report 195 test scenarios, namely: better than 2% for both test problems in case 4 using the PMMA phantom, and better than 5% of the reference result for 14 of 17 organs in case 5, and within 10% for the three remaining organs. The Fluka simulation results agreed to better than 2% of the air kerma measured in‐air at isocenter of the GE Revolution scanner. The simulated air kerma in the center of the CTDI phantom overestimated the measurement by 7.5% and a correction factor was introduced to account for this. The simulated mean absorbed doses for a chest scan of the pediatric anthropomorphic phantom was completed in ~9 min and agreed to within the 95% CI for bone, soft tissue, and lung measurements made using MOSFET detectors for fixed current axial and helical scans as well as helical scan with TCM. Conclusion A Fluka‐based Monte Carlo simulation model of axial and helical acquisition techniques using a wide‐beam collimation CT scanner demonstrated good agreement between measured and simulated results for both fixed current and TCM in complex and simple geometries. Code and dataset will be made available at https://github.com/chezhia/FLUKA_CT.


| INTRODUCTION
The use of radiation based imaging systems has become indispensable for medical diagnosis leading to a rise in diagnostic radiation.
However, long term health effects of cumulative diagnostic radiation exposure are not yet completely understood. The evaluation of possible long-term effects from diagnostic radiation will only be possible once fast and accurate patient-specific dose quantification is available. Current practice simply archives machine-dependent, radiation output reports containing examination-specific and phantom-based dosimetry information such as volume computed tomography dose index (CTDI vol ) and dose length product (DLP). Calculations exist to convert these phantom-based radiation output metrics to dose estimates more representative of patient dose, such as size specific dose estimates (SSDE) 1,2 which have been used to estimate organ specific doses 3,4 and population-specific effective dose from SSDE and DLP. 5 These techniques, however, are still based on approximations and are not patient-specific. Since patient size, anatomy, geometry, and tissue composition can vary greatly among patient populations, even within similar ages and weight groups, quantifying patient-specific radiation dose is highly preferable.
For radiation transport problems with complex geometry and varying material properties, Monte Carlo simulation-based dose quantification techniques remain the gold standard despite long computational times. With the recent advent of lower-cost high performance computing resources, using Monte Carlo simulations for patient-specific radiation dose reporting within the clinical environment is more promising than ever before. Additionally, recent studies have implemented Monte Carlo based dose simulations in parallel computing environments, such as Graphical Processing Units (GPUs). 6,7 The implementation of high-performance, parallel computing may lead to patient-specific radiation dosimetry calculation times that become feasible even during demanding clinical schedules.
The purpose of this study is to develop a fast, multi-mode (i.e., helical and axial) CT scanner simulation model of a wide-beam CT scanner including tube current modulation (TCM). The radiation transport model for the CT scanner is developed using the open source Monte Carlo code Fluka. Although previous studies have utilized Fluka in applications such as radiotherapy 8 and nuclear medicine, 9 Fluka has not been validated for diagnostic dose calculations in CT. As the same voxel geometry and organ segmentation generated for radiotherapy dose simulations can be used for diagnostic dose calculations, the model developed in this work can be extended to calculate total radiation dose for radiotherapy patients using Fluka is a general purpose particle transport calculation tool that can track a wide range of charged and neutral particles including photons and electrons with energies ranging from 100 eV to thousands of TeV. Fluka also includes stable physical models to handle various particle interaction mechanisms across a variety of material types.
Fluka (Version 2011.2c-5) supports complex geometries through a combinatorial geometry package and has options for user written routines that were used to customize the simulation to suit specific problems in this study. The source routine in Fluka was highly modified to support various operation modes of the CT scanner, for example, axial and helical. The input simulation files for benchmarking Fluka accuracy were generated using an advanced user interface called Flair. Flair has an inbuilt routine that allows voxelization of CT images into a format that is understood by Fluka.

2.A.2 | Source spectrum generation
The energy spectrum of the x-ray tube for the GE Revolution CT scanner was not readily available from the vendor as it was proprietary information. To generate the source x-ray photons for the Monte Carlo simulation, the tube x-ray spectrum was calculated using SPEKTR 3.0 10 and optimized to match the measured spectrum of the GE Revolution scanner. The SPEKTR toolkit calculated x-ray spectra based on the Tungsten Anode Spectral Model using Interpolating Cubic Splines (TASMICS) algorithm. It generated the x-ray spectra in 1 keV energy bins over beam energies from 20 to 150 keV. The spectrum calculated from SPEKTR was given in photons/mm 2 /mAs at 100 cm from the source. The spectrum at the isocenter of the scanner (62.56 cm from the tube for the GE Revolution) was then calculated by accounting for geometric attenuation. The TASMICS spectrum from SPEKTR represented the unfiltered spectrum from the x-ray tube. In modern CT scanners, the x-ray spectrum undergoes some inherent filtration after exiting the x-ray tube housing and to determine this filtered spectrum, SPEKTR included an optimization function that matched the exposure of the final spectrum to the measured x-ray tube exposure (mGy/mAs). In this study, the SPEKTR generated spectrum was optimized to match both exposure and Half Value Layer (HVL) of the measured spectrum using MATLAB's [2016a] optimization toolbox. The optimization involved generating an equivalent filter that provides a filtered spectrum that matches the measured spectrum. 11 The x-ray spectrum at 120 kV tube potential measured at the iso-center of the scanner was optimized with an equivalent filter made of carbon and aluminum with 3.6 and 9.3 mm thickness, respectively. The electron and photon cut-off SOMASUNDARAM ET AL. | 133 energies for the Fluka simulation were set to be 1 and 0.5 KeV for all material regions.

2.A.3 | Bowtie filter characterization
In addition to the inherent tube filtration, a bowtie shaped filter is used to absorb the lower energy photons of the spectrum such that the periphery of the body is not over exposed to radiation and to maintain uniform image quality across the scan volume. The shape of the bowtie filter was characterized by measuring the beam intensity along the radial axis of the scanner gantry. 11 The tube was parked at 90°using the service mode of the scanner and a 0.6 cc ion chamber (10x6-0.6, Radcal, Monrovia, CA) was mounted on a jig placed on the patient table such that the center of the chamber aligned with the isocenter of the scanner. Measurements of the in-air exposure were then made at intervals of 1 cm from isocenter by adjusting the table height upwards until the ion chamber was at the edge of the bowtie filter towards the top of the CT gantry. The measured exposure values were normalized and used to optimize the thickness profile of the bowtie filter such that the equivalent spectrum matched that measured spectrum along the radial axis. Once the thickness profile of the bowtie was determined, the filtered spectrum exiting the bowtie at each of the 1 cm intervals was calculated. Based on the total number of photons emitted in each bowtie interval, a discrete probability distribution function of the bowtie profile was created. During the Monte Carlo simulation in Fluka, the initial direction of the x-ray photons was sampled using the bowtie distribution function. The bowtie interval, through which the particle would exit the filter, was sampled from the probability distribution function, and uniform random sampling was used within the interval to select the exact coordinates (XY) that would be traversed by the exiting particle. The initial polar angle of the particle was then calculated from the XY coordinates.

2.A.4 | Beam profile characterization
The heel effect, that is, the variation in beam intensity (or beam profile) along the long axis (z-axis) perpendicular to the radial axis of the scanner, was more pronounced with the GE Revolution scanner when operated using wide collimation (up to 16 cm) compared to previous publications covering only 4 cm along the z-axis. 12 The beam profile was measured by exposing a radiochromic film at the desired collimation with 375 mA, 120 kV and 1 s tube rotation. Exposure values in the range between 350 to 375 mAs were optimal for beam width measurements using radiochromic film and were verified using ionchamber measurements. 13 The exposed film was digitized using an Epson flatbed scanner (10000XL) and processed in a custom-built MATLAB [2016a] script to determine the FWHM of the beam profile.
The beam profile FWHM was considered the actual collimation width for the simulation, and was discretized into 1 cm intervals and converted into a discrete probability distribution of the x-ray intensity.
The probability function was used to sample the z-coordinate, and the initial azimuthal angle, with which the source photons would reach the patient volume. Combining the initial azimuthal angle along with the polar angle, sampled from the bowtie profile distribution, and yielded the initial particle direction of the source photons.

2.A.5 | Tube rotation and table movement
The actual tube rotation in a scanner can usually be approximated as small fixed steps around the circumference of the gantry, but the exact number of steps to complete a single rotation was unavailable.
Therefore, in this study, the source photons were sampled uniformly along the circumference of the rotation path with the tube-to-isocenter distance as the radius.
For helical scans, the table movement was simulated by translating the x-ray tube along the z-axis with the table and the patient held at a fixed position. This was achieved by uniformly sampling the z-location for the source particle along the range of table movement with respect to the x-ray tube. The range of the table movement T z was given by: where, P was the pitch of the helical scan, Col was the x-ray tube collimation, Exp time was the total exposure time for the scan and Rot time was the time taken by the tube to complete one full rotation.
The z-location for the source particle, z p was sampled as: where, ε was a random number between 0 and 1, z i was the starting location of the tube with respect to the patient/table. The ± sign was used to indicate that the table movement could be in the positive or negative direction depending on a head first or feet first scan.
Once z p was sampled, the radial location (x p , y p ) of the particle was determined as: Where, R iso was the tube to isocenter distance for the scanner and β p was calculated as: β i was the initial angle (in radians) made by the x-ray tube with respect to the y-axis of the scanner.

2.A.6 | Tube current modulation (TCM)
The parametric details to simulate TCM were not available in the DICOM header metadata. To model TCM in this study, the average beam current (mA) recorded for each reconstructed slice was used to construct a probability distribution along the z-axis of the scan and the source photons were sampled from this distribution to simulate the TCM profile. The TCM algorithm implemented in the Fluka model did not simulate the rotational beam current modulation within each reconstructed slice (i.e., modulation in the x-y plane).

2.B | Primary fluence calculation
The primary x-ray fluence of the x-ray tube derived from Fluka was tallied as dose per primary particle and this quantity was scaled by the total number of primary photons emitted based on kV, mA, scan duration, filter settings, collimation, and heel effect. The total number of primary photons was used to estimate the total simulated dose/exposure for comparison with the measured dose/exposure values to validate the simulation model. The following formula was used to calculate the total source strength (Q): where, N B was the total number of bins used to characterize the bowtie filter, N H was the total number of bins used to characterize the beam profile, H j was the normalized beam intensity in bin j, N K was the total number of energy groups handled by SPEKTR, A ij was the rectangular area mm 2 formed by the length of the bowtie filter bin i and length of the beam profile bin j, D ij was the distance from the x-ray tube to the center of the rectangular area A ij , SDD was the distance from the x-ray tube to isocenter of the scanner, and finally,  The interpolation functions were used to determine the conversion maps to calculate densities and elemental weights for the entire range of Hounsfield values (HU) in the DICOM images. Figure 1 shows the density calibration curve calculated for the GE revolution with phantom measurements made at a tube voltage of 120 kVp, while Table 1 shows the material composition in each HU range.
The exact geometry of the patient respectively of AAPM Report 195 and the corresponding simulation parameters for tests 1 and 2 are given in Table 2. The simulation parameters are shown in Table 3, and the geometry setup in fig. 17   from NIST XCOM photon cross section databases. 16 The method of tallying the energy fluence instead of air kerma (exposure) directly helped speed-up the simulation by allowing the use of a smaller particle history size of 10 7 compared to~10 9 while maintaining the statistical precision to better than 1% for all the test cases.

2.D.3 | CTDI phantom measurements and simulation
To test the performance of the proposed Fluka source model in simulating a CTDI measurement experiment, a 15 cm long CTDI body phantom (32 cm diameter) was exposed using the same scanner settings in Table 4. The exposure measurements were made using the  Figure 4 demonstrates the phantom fully loaded with MOSFET dosimeters, and the respective prescribed CT scan range.  80 mm using both TCM as well as fixed current modes. CT scan acquisition parameters for both modes were shown in Table 5; for the helical scan with TCM, the mA values varied between 150 and 450 mA and the Noise Index of the image was 4.
The MOSFET detectors were cross-calibrated with an x-ray machine using the methods previously described. 17 The HVL of the x-ray machine was matched to that of the CT machine used by placing additional aluminum plates at the exit port of the x-ray machine to harden the x-ray beam for a HVL of 7.62 mm Al at 120 kV. The average calibration error was 10% (range 4%-16%) for the 25 MOS- the tally size was chosen such that they were large enough to achieve a statistical precision to better than 7%, and small enough such that they did not overlap with other neighboring tissues; since the MOSFET measurement errors had a standard deviation of~10%, a more accurate simulation was not required.

2.E | Parallel execution of Fluka
The

3.B | AAPM Report 195 benchmark tests
Fluka results for the four contiguous cylindrical detectors of case 4 test 1 agreed to better than 1.5% of the benchmark result for both monoenergetic and spectral source at 10 and 80 mm collimations.
Fluka results for case 4 test 2 agreed to better than 2% of the benchmark result for the detector in the periphery and better than 1% for the center detector at both 10 and 80 mm collimations. For case 5, the static projection of the monoenergetic source at 0°resulted in a better than 3% agreement between Fluka results and the benchmark for all organs except adrenal gland and breasts for which the results agreed to better than 4%; and for the spectral source, Fluka results agreed to better than 3% for all the organs except the thyroid and adrenal gland which differed by 5% and 6%, respectively. Fluka results for case 5 with tube rotation and monoenergetic source agreed to better than 5% for all the organs except thyroid (9.8%) and breasts (5.4%) in comparison to the benchmark results. For the spectral source, Fluka results agreed to within 5% of the benchmark for 14 of the 17 organs. The results for the adrenal gland, the thyroid and the breasts differed by 8.7%, 7.1%, and 5.9%, respectively. Table 7 shows the relative deviation of Fluka results to the benchmark for all the source types.
Overall, Fluka results agreed to better than 4% to the reference results for the cylindrical phantom simulation and for 14 of the 17 organs in the XCAT phantom simulation.  Table 9 shows the comparison of the simulated exposure values in the CTDI phantom to the measured exposure using the pencil chamber placed at the center hole of the phantom along with the correction for the patient table.

3.D | CTDI phantom measurements and simulation
The simulation results exceeded the measured exposure by approxi-mately~7.5% for both tube collimations (7.1% for 8 cm collimation and 7.6% for 16 cm collimation). To correct for this effect, a correction factor (CF CTDI ) based on the ratio of the measured and simulated exposure values at the center of the CTDI phantom was calculated: Since the dose in the anthropomorphic phantoms is also measured in a similar scattering environment to a CTDI phantom, the simulation results for the anthropomorphic phantoms were multiplied by the CTDI correction factor (CF CTDI ).
3.E | Physical anthropomorphic phantom measurements and simulation Where, r cf is the relative error in the calibration factor with respect to the measured dose.

| DISCUSSION AND CONCLUSION
In this study, the characterization of the GE Revolution CT scanner and development of a Fluka Monte Carlo package were created for simulating tissue absorbed doses in both a wide-beam axial mode (i.e., 160 mm) and 80 mm helical mode with and without TCM.
Experimental measurements and computational tools were utilized to quantify the x-ray spectrum, bowtie filter, beam profile, and primary photon fluence of the scanner using the large body filter and a tube voltage of 120 kV.  and 9.8%, respectively) from the reference results is likely due to the differences in particle interaction physics between Fluka and other codes. This was further confirmed by a high accuracy simulation (10 9 photons) of the test problem with spectral source and random tube rotation. When compared to the simulation with 10 8 photons, the statistical errors for the 10 9 simulation for the adrenal gland, thyroid, and breasts dropped from 4.8%, 1.2%, and 0.4% to 1.5%, 0.5%, and 0.2%, respectively. The results for the 10 9 simulation in comparison to the reference means were within 5.5%, 6.5%, and 5.6% for the adrenal gland, thyroid, and breasts. Except for the adrenal gland, the high accuracy simulation did not cause a significant improvement (<1%) to Fluka results in comparison to the reference results. For the adrenal gland there was a 3.1% (5.5% for 10 9 and 8.6% for 10 8 histories) improvement.
The discrepancy in the simulated exposure at the CTDI phantom center (leading to the introduction of the CF CTDI factor) is due to systematic errors in the estimated x-ray spectrum. These errors could be due to subtle differences in beam hardening caused by the actual filter materials compared to assumptions made to generate an equivalent spectrum. Air kerma measurements at isocenter agreed well with simulated exposures indicating that the ion chamber does not capture the difference in spectrum when the significant dose contribution is from the primary x-ray fluence emitted directly from the tube. But for air kerma measurements in the center of the CTDI phantom, a significant dose contribution comes from the secondary radiation caused by scattering interactions in the phantom material.
Since the spectral differences are likely to affect the amount of scattered radiation produced, the measurement and simulation results at the center of the phantom could lead to a larger discrepancy.
When comparing the MOSFET measurements in the anthropo-

This work was supported by American Lebanese Syrian Associated
Charities (ALSAC).

CONFLI CT OF INTEREST
The authors declare no conflict of interest.