Fast and accurate Monte Carlo modeling of a kilovoltage X-ray therapy unit using a photon-source approximation for treatment planning in complex media.

To accurately recompute dose distributions in chest-wall radiotherapy with 120 kVp kilovoltage X-rays, an MCNP4C Monte Carlo model is presented using a fast method that obviates the need to fully model the tube components. To validate the model, half-value layer (HVL), percentage depth doses (PDDs) and beam profiles were measured. Dose measurements were performed for a more complex situation using thermoluminescence dosimeters (TLDs) placed within a Rando phantom. The measured and computed first and second HVLs were 3.8, 10.3 mm Al and 3.8, 10.6 mm Al, respectively. The differences between measured and calculated PDDs and beam profiles in water were within 2 mm/2% for all data points. In the Rando phantom, differences for majority of data points were within 2%. The proposed model offered an approximately 9500-fold reduced run time compared to the conventional full simulation. The acceptable agreement, based on international criteria, between the simulations and the measurements validates the accuracy of the model for its use in treatment planning and radiobiological modeling studies of superficial therapies including chest-wall irradiation using kilovoltage beam.

Introduction

Standard local control rates have been observed in breast cancer patients undergoing post- mastectomy chest-wall irradiation by kilovoltage (kV) X-rays.[1] The technique, used more widely at our center in the past due to a shortage in the availability of linear accelerators, employs a single 120 kVp X-ray beam angled medially by 15 degrees and positioned usually at a 50 cm focus-to-skin distance (FSD). The origins of this treatment technique date back a few decades and its dosimetric basis is unclear. Therefore, it needs to be revisited by accurate calculation of dose distributions. In particular, accurate dose distributions are needed to perform valid radiobiological modeling analyses of this technique.[2]

The complexity in this technique from the dose calculation point of view arises from the presence of beam obliquity, patient contour irregularity and radiation field splash-over, as well as lung and bone tissue heterogeneities. Therefore, the significance of 3D treatment planning to ensure target dose coverage and reduce organ-at-risk dose is undeniable. Most treatment planning systems are, however, designed for megavoltage beams. Similarly, the Monte Carlo (MC) method has been extensively used to perform dosimetric investigations in the megavoltage region and only a limited number of studies have been reported regarding MC modeling of kV X- ray therapy.[3] A few previously published studies have reported on measured X-ray spectra[4567], and various programs have been developed for the calculation of kV X-ray spectra and half-value layer (HVL) values based on calculated spectra.[89] For example, a miniature radiosurgery kV tube was simulated by Yanch and Harte using ITS Version 3.0 ‘p’ codes.[10] Some useful studies were also published by Verhaegen et al. They used MCNP4B and EGS4/BEAM to simulate kV X-ray units to obtain accurate information for absorbed dose calculations in water for kV X-rays. In particular, the photon spectrum, planar fluence and the angular distribution of the photons at the collimator exit of the X-ray unit were studied.[311] Another application in the kV energy region was published by Hill et al., which simulated a kV X-ray beam and water phantom to calculate deep and superficial doses using the EGSnrc code. The X-ray beams were 75 − 135 kVp with field sizes of 2, 5 and 8 cm diameter.[12] Knoos et al., also simulated an orthovoltage unit using EGSnrc. The BEAMnrc code was used to transport electrons, produce X-ray photons in the target and transport them through the treatment machine down to the exit level of the applicator. Further transport in water or CT-based phantoms was simulated by using the DOSXYZnrc code.[13] Most of the studies mentioned above used EGS4, EGSnrc (or BEAMnrc) or MCNP4B. BEAMnrc is faster in calculations and it is easier to use in defining treatment head geometries but MCNP4C offers more flexibility in defining complex geometries.[14] In all previous studies, all the relevant components of the X-ray tube were simulated.

The purpose of this study was to develop a fast reliable dose calculation engine, which can be used in chest-wall treatment planning. In contrast, the main purposes of the previous studies were assessment of kV X-ray beam features, which were mostly time-consuming. This study adds evidence that it is possible to simulate a kV beam without simulating the components and processes leading to production of X-rays with sufficiently high accuracy even in realistic inhomogeneous media. This simplification has benefits in terms of time needed to design the MC model and the simulation run time. This approach is also useful when some details regarding the tube are inaccessible or not known.

In this study, the beam properties in an inhomogeneous phantom together with beam obliquity and contour irregularity were also assessed, which had not been performed in most of the previous studies with kV therapy beams. This study also adds further evidence to validate the accuracy of MCNP4C in inhomogeneous anthropomorphic phantoms for 120 kVp X-rays.

The main purpose of this study was the development and validation of a fast MC model of a superficial X-ray beam to be used as a treatment planning tool for accurate calculation of dose distributions delivered by kV techniques and its optimization in the shortest possible time. In particular, validation of the developed model for the relatively complex irradiation of post- mastectomy chest wall using the above-mentioned technique is of utmost interest. Given its strengths in prediction of dose distributions in heterogeneous media and other complex situations, MC modeling is a desirable tool for treatment planning for challenging geometries.[1516171819]

A routine method to accurately model X-ray tubes is simulating their almost every detail such as target, filter, mirror, etc. The secondary purpose of this study was to check the feasibility of implementing a method that allows us to model the X-ray beam without simulating all the details in the X-ray tube (namely, a photon-source approximation; using a photon source instead of an electron source in the simulations). The reduction in the number of modeled components would be beneficial in terms of (i) Future designing of MC X-ray tube models and (ii) their simulation times.

Materials and Methods

Our kV X-ray therapy machine is the Stabilipan unit (Siemens, Germany) with a tungsten target angled 22 degrees and a total filtration equivalent to 2 mm Al. Collimators (shutters), instead of applicators (cones), were used on this unit to adjust the field size.

Measurements

In order to be usable as a treatment planning tool, the accuracy of the MCNP beam model had to be checked by comparison of its main dosimetric properties such as HVL, percentage depth dose (PDD) and off-axis beam profiles with the actual beam used for therapy. Therefore, the following measurements were performed at 120 kVp:

HVL measurements were carried out using high-purity Al foils and an ion chamber. This was performed as an initial check of the beam spectrum. The first and second HVLs were measured using high-purity Al foils placed at 50 cm FSD, midway between the tube and a Farmer-type ion chamber dosimeter (IBA Medical AB, Sweden). The field size at 50 cm was 2 × 5 cm2 to create as small a field as possible (narrow-beam geometry) while allowing sufficient margin around the chamber for electronic equilibrium

PDDs and beam profiles were measured for two field sizes (8 × 8 cm2 and 14 × 14 cm2) and FSD as 50 cm for normal incidence on a 60 × 60 × 60 cm3 scanning water phantom (IBA Medical AB, Sweden). Electron dosimeters such as diodes and NACP parallel-plate ion chamber are suitable detectors for measuring PDDs and profiles in water medium at 120 kVp.[20] Therefore, PDDs and beam profiles were measured using an electron diode (Type EID; IBA Medical AB, Sweden). The PDD measurements were repeated with an NACP chamber (IBA Medical AB, Sweden) in order to provide complementary measured data. Beam profiles were measured at two depths (2 and 5 cm) and in two cardinal directions (inline and crossline; perpendicular and parallel to the anode-cathode axis, respectively). All data points were normalized on axis at the depth of 2 cm, to avoid dose measurements within the less-reliable depths shallower than the radius of the detector[3]

Further validation experiments were carried out by dose measurements for a more realistic and complex situation of oblique beam incidence on the chest wall using 10 thermoluminescence dosimeters (TLDs) placed within the chest wall, heart and lung regions of an inhomogeneous anthropomorphic Rando phantom. The TLDs were calibrated in the same beam quality by intercomparison with the above-mentioned calibrated Farmer-type ion chamber. Meticulous TLD methodology using a previously described technique was employed to keep the measurement uncertainty within 2%.[21] The treatment field was defined on the Rando phantom exactly like a real patient and a routine treatment set up was performed with 50 cm FSD and 15-degree tube angle medially.

Monte Carlo simulations

The MCNP4C code[22] was used for simulation. MCNP4C has some new features compared to MCNP4B, one of them being electron physics enhancements, which make this code more useful for this study. Improvements made include density effect calculation for stopping power, radiative stopping power, Bremsstrahlung production (spectra intensity and angular distribution) and impact ionization, a new electron library (EL03)], etc.[22] The latest version of MCNP available to us was version 4C. There are more recent versions of MCNP (versions 5 and X). However, the small differences between the results of version 4C and these later versions, together with the reasonably good agreement with experimental measurements presented in this paper, justify our use of MCNP4C.[2324] An updated photon cross-section library, ENDEF/BVI release 8, was used here. For low-energy photon dosimetry, this library offers an improvement over the previously released MCNP4 DLC-200.[24]

In order to reduce the run time significantly, these simulations were simplified by using a photon source instead of an electron source. Instead of modeling the machine components and an electron source to create the X-ray spectrum needed for simulations, a spectrum processor[25] was used. By entering information such as target material, tube voltage, anode angle, and filter material and thickness, the processor could calculate the X-ray spectrum in 0.5 keV intervals.

Using a pinhole technique, the tube's focal spot dimensions measured by film were 7 × 7 mm2. In order to be able to design a photon source that could model the heel effect correctly, a cylinder subroutine source was defined. The photon source was simulated as a 7 mm long and 7 mm diameter cylinder with its long axis along the anode-cathode direction, bisected in length by the beam's central axis, and divided in length into 7 slices by planes at 1 mm intervals. This allowed the poly-energetic source to be biased along its length, i.e., to have a gradual decrease in weight of X-ray production toward its anode-side end, to mimic the loss of X-ray production observed at the anode side in heel effect. This replaced the more routine method of modeling an electron source, target and simulating both electron and photon transport within the target to produce the observed heel effect.

All of the experimental measurements were simulated in MCNP4C using the same geometries:

The first and second HVLs were calculated from simulation results with the same geometry as the experimental measurement. The target-Al-foils and detector-Al-foils distances were both 50 cm. The detector cell was a cylinder of the same dimensions as the sensitive area of the ion chamber (6.25 mm inner diameter and 24 mm length)

PDDs and beam profiles were calculated by simulating a 60 × 60 × 60 cm3 water phantom placed 50 cm from the target. Data cells were in the form of a fine lattice at the central axis for PDD computations, and 4 horizontal lattices in inline and crossline directions at the depths of 2 and 5 cm for profile calculations. In order to have a high spatial resolution in the calculated PDDs, the lattice-cell dimension was set at 1 mm in the direction of each measurement

The exact internal and external contours of the Rando phantom were obtained and measured from CT images and the same geometry was created by carefully defining the planes in the MCNP input file. Data cells were defined exactly in the locations where TLDs were placed during experimental measurements. In order to reduce run time to a practical level while maintaining a low statistical uncertainty, variance reduction techniques were utilized. Geometry splitting was implemented in which the geometry between the tube and data cells was split by concentric cylinders, the importance values of which increased gradually (in relatively small steps: 1, 6, 36, 216) toward the center (to increase the probability of interactions nearer the data cells without biasing the end results).[26] Also, energy cut-offs were used; electrons were transported in the phantom down to the energy of 10 keV and the photon transport cut-off energy was 1 keV.[3] The results of the experimental measurements and MCNP simulations were compared to evaluate the accuracy of the model. To the best of our knowledge, there are no specific criteria of acceptability for treatment planning calculations for kV beams. The criteria of acceptability used in these comparisons were those set out in the International Atomic Energy Agency (IAEA) TECDOC-1583[27] and the American Association of Medical Physicists (AAPM) TG-53 Report 62,[28] the scopes of which do not exclusively cover kV beams.

Results and Discussion

The first and second HVLs were 3.8 and 10.3 mm Al, respectively [Figure 1], which were in very close agreement with the values obtained by MCNP (3.8 and 10.6 mm Al). This constituted the first test of the X-ray spectrum modeled.

Figure 1

Transmission curves resulting from experimental measurement and Monte Carlo simulation to obtain the HVLs of the beam

The comparison between the PDDs measured by the electron diode and the NACP chamber showed that their measured PDDs agreed within 2% at all data points, the agreement for majority of the points were lying within 0.4 − 0.8% [Figure 2]. This close agreement is interesting, given the spectral changes with depth in water and the well-known higher-energy dependence of diodes compared to ionization chambers. However, our results indicate that this effect is negligible with the combination of the beam spectrum, electron diode and depth range (up to 15 cm) used. Therefore, comparisons were made with the diode for the rest of the measurements presented here.

Figure 2

Comparison between PDDs resulting from diode and NACP dosimeters in 8 × 8 cm2 and 14 × 14 cm2 fields

According to the IAEA[27] and the AAPM,[28] acceptable differences between measurements and calculations on the central axis are 2% and 1%, respectively. The differences in both off-axis and outer beam regions in a homogenous phantom should be within 3% according to the IAEA and within 1.5% and 2%, respectively according to the AAPM. In the penumbra region, the differences should be within 2 mm according to both criteria.[2728] The results are presented in Figures 3 to 7. The differences between the measured and calculated PDDs and beam profiles in water were less than 2 mm/2% for all data points, which constitutes generally acceptable agreement between calculations and measurements by at least one of the above-mentioned criteria. The agreement in PDDs is not that good for the first few millimeters, which may be due to increased uncertainty in the measurements at very shallow depths. In oblique chest wall irradiation, the PDD at depths up to 15 cm are of interest. These differences in PDD and profile results are similar to those from a study by Knoos et al.[13]

Figure 3

Comparison of PDD curves resulting from experimental measurements and Monte Carlo simulations in 8 × 8 cm2 and 14 × 14 cm2 fields

Figure 7

Comparison of beam dose profiles resulting from experimental measurements and Monte Carlo simulations in an 8 × 8 cm2 field at depths of 2 and 5 cm in the inline direction

Comparison of beam dose profiles resulting from experimental measurements and Monte Carlo simulations in a 14 × 14 cm2 field at depths of 2 and 5 cm in the inline direction

Comparison of beam dose profiles resulting from experimental measurements and Monte Carlo simulations in a 14×14 cm2 field at depths of 2 and 5 cm in the crossline (anode-cathode) direction. Negative positions refer to the anode side

Comparison of beam dose profiles resulting from experimental measurements and Monte Carlo simulations in an 8 × 8 cm2 field at depths of 2 and 5 cm in the crossline (anode-cathode) direction. Negative positions refer to the anode side

The measured and calculated data for the inhomogeneous Rando phantom were also compared. Table 1 shows the results of dosimetry in a real and simulated Rando phantom. According to the criteria of acceptability recommended by the IAEA and the AAPM,[2728] acceptable differences between measurements and calculations on the central axis for an inhomogeneous 3D phantom are 3% and 5%, respectively. In off-axis regions, the acceptable differences are 3% and 7%, respectively. MCNP-calculated relative doses in the Rando phantom agreed with TLD measurements within 4.5% of the prescription dose; the differences for majority of the data points were within 2%. Again, this shows generally acceptable agreement between calculations and measurements by at least one of the above-mentioned criteria.

Table 1

Mean relative doses resulting from experimental measurements using TLD and Monte Carlo modeling using MCNP

Our findings disagree somewhat with those of Verhaegen et al.[3] who concluded that complete set up (focal spot size, target, inherent filtration and collimators) must be modeled accurately in order to be able to reproduce the measured dose distribution.

Run time and time saving

The simulation of PDD using the MCNP code took about six hours to perform 109 histories with an estimated relative error of about 0.25 to 0.3% in the high-dose regions, when it was running on a Pentium 4, Intel Core 2 Duo CPU P8700 2.53 GHz CPU. To compare the run time when a photon source was designed instead of an electron one, the complete set up of X-ray tube including an electron source was simulated (conventional method) and beam profile computation was repeated using this simulation. An example of the number of histories (NPS) and the computational time (CTME) required to reach a relative error (R) value of 0.1 in MCNP4C when simulating an electron source (conventional method) and a biased photon source (proposed method) are presented in Table 2. From this Table, it can be seen that the run time of the biased photon source is more than 9500 times less than that of the conventional method.

Table 2

Example of the number of histories (NPS) and the computational time (CTME) required to reach a relative error (R) value of 0.1 in MCNP4C when simulating an electron source (conventional method) and a biased photon source (proposed method)

There is no consensus on whether using energy cut-offs should be called a variance reduction technique in MC simulation. Using cut-off energies is certainly not as elaborate as geometry splitting, for instance. However, it can reduce the runtime as it stops the simulation of unnecessarily small steps occurring at very low energies, thereby reducing the variance of the result obtained for a given number of histories. The MCNP4C manual includes it as a variance reduction technique.[22]

Our proposed solution is an alternative to the widely used method of saving phase-space files of the photons emerging from the target (following full simulation of electron and photon interactions in the target). The aim and scope of this paper does not include an investigation of the relative merits of photon-source approximation compared to saving phase-space files. However, obviation of the need to manage huge phase-space files and a reduced requirement for having and implementing detailed geometry information can be mentioned as advantages of the photon-source approximation method.

The superficial X-ray therapy technique described here is one of several possible alternatives for chest-wall radiotherapy.[29] Discussing the advantages and disadvantages of this technique compared to the others is, however, outside the scope of this paper.

Conclusions

The proposed, simplified method of simulating a kV X-ray therapy unit offers huge time savings in terms of modeling and run time, which is necessary for implementation of MC treatment planning. The acceptable agreement between the results of the simulations and the measurements validates the accuracy of the MCNP model despite the exclusion of the components and processes leading to production of X-rays, for use in treatment planning and radiobiological modeling studies of chest-wall irradiation using kV beams. In addition, this model can be used for other (often simpler geometry) superficial treatments.