A 3D model to calculate water-to-air stopping power ratio in therapeutic carbon ion fields.

INTRODUCTION

Air-filled ionization chambers (ICs) are extensively used in the dosimetry of charged particle radiotherapy [1]. The calibration procedure of ionization chambers for the determination of absorbed dose to water, which is the standard quantity used for dose determination in external radiotherapy [2] is known as N formalism. In this formalism, the readout of the chamber is converted into absorbed dose to water via two factors: the calibration factor of the chamber, and a quality factor that accounts for the specificity of the beam. The water-to-air stopping power ratio, or s, is one of the main components of these quality factors, and, in the case of carbon ion beams, its biggest source of uncertainty [2].

In a previous work by our group [3], an expression was proposed to calculate s for carbon ion beams at different residual ranges, based on a set of Monte Carlo calculations and experimental measurements, namely: where R is expressed in cm and calculated using a practical range at the 50% dose level R(z) = R – z, where z is the depth in water. This expression is based on a 1D analysis of dose and s distributions, which is enough to model the variations in s for homogeneous dose distributions, like the ones mostly used for calibration and quality assurance (QA) purposes. However, this 1D description might be insufficient in some cases. An example of this would be treatment plan verification with a matrix of ionization chambers [4], a protocol often used in scanning-beam facilities where a patient plan is shot into a water phantom and the deposited dose is measured at several points (see Fig. 1). In such a case, the residual range R is not defined at every point, so the application of equation (1) is not possible.

Fig. 1.

PTW matrix of ionization chambers (PTW Freiburg, Germany) used for plan verification and QA measurements of proton and light ion beams, inside a PTW Water Phantom.

This contribution addresses the applicability of the proposed expression in 3D homogeneous and inhomogeneous dose distributions. In particular, we present a procedure to calculate s from R at every point, in realistic treatment plans of arbitrary geometry. The validity of the procedure is cross-checked with detailed Monte-Carlo calculations of s carried out with FLUKA [5, 6] in two 3D geometries, a QA plan consisting of a cube with homogeneous dose, and a real field of a treatment plan from a head-and-neck patient. The resulting s maps are also studied on a voxel-by-voxel basis, in order to quantify how well the simplified procedure matches the values calculated with the full Monte Carlo simulation.

MATERIALS AND METHODS

Extension of the model to 3D

In order to apply expression (1) to calculate stopping power ratios at an arbitrary position of a treatment field, a value of R, and consequently a value of R, have to be determined for every voxel. The process is depicted in Fig. 2. For each voxel (small box in the figure), its transversal coordinates {x,y} are matched to the closest raster point from the treatment plan (shown with stars in the diagram). The raster points are the control points of the scanned beam plan. For each point, a set of pencil beam energies is prescribed. Of these energies, the highest will determine the range associated with that raster point, and therefore, the derived value of R will be assigned to the original voxel at {x,y}, and used to calculate its R This algorithm is designed to be implemented in the treatment planning system (TPS), as it uses information from the treatment planning physical base data.

Fig. 2.

Scheme depicting the determination of the R at an arbitrary point of the field by looking for the closest raster point in the treatment plan.

Fig. 3.

Stopping power ratios calculated with FLUKA (left) and using the proposed model (right), for a QA plan consisting of a 3-cm side cube of homogeneous physical dose of 1 Gy, centered at a depth of 5 cm, on a water phantom. The vertical section at the isocenter plane is displayed. Continuous lines represent isodose levels in the physical dose distribution.

Fig. 4.

Stopping power ratios calculated with FLUKA (left) and using the proposed model (right), for a real treatment plan from a head-and-neck patient, simulated on a water phantom. The vertical section at the isocenter plane is displayed. Continuous lines represent isodose levels in the physical dose distribution.

Calculation of s with FLUKA

A full calculation of s at the point has to take into account the complete energy and fragment spectra of the field [2], weighting the stopping powers, S(E)/ρ, in air and water, with the point fluence in water at any given point, , yielding:

Using equation (2), we applied the implementation described in [2] in the Monte Carlo code FLUKA 2011.2 [5, 6]. The most relevant parameters of the simulation were the HADROTHERAPY set of defaults from FLUKA, including a transport threshold Δ of 25 keV/u for charged hadrons, 10 keV for delta ray production and transport (as discussed in [3]), and mean ionization potentials of Iw = 78 eV for water and Iair = 90.2 eV for air [7, 8].

The proposed algorithm was tested on two carbon ion fields. The first one (Fig. 3) featured a cube with an edge length of 3 cm placed at a depth of 5 cm in water, irradiated with a dose of 1 Gy. The second one (Fig. 4) consisted of a real plan from a head-and-neck patient irradiated at the Heidelberg Ion-Beam Therapy Centre (HIT). For both fields, the water-to-air stopping power ratio was calculated on every voxel, using both the described algorithm (see ‘Extension of the model to 3D’) and the full Monte Carlo simulation (see ‘Calculation of s with FLUKA’). The voxel sizes were set at 4 mm in the transversal directions and one-third of a millimeter in the longitudinal direction, in order to accurately describe the sharp rise in s at the distal edge of the peak.

RESULTS

In order to further quantify the similarity of the two s maps, we analyzed them on a voxel-by-voxel basis. The same analysis was done with the currently recommended constant value of sw,air = 1.13 [2], in order to evaluate the performance of the two approaches. In order to do so, the histograms displayed in Fig. 5 were constructed. They represent the relative difference between the FLUKA-calculated stopping power ratios, and the s values calculated by the described procedure (labeled as ‘Model’) and given in the TRS-398 recommendation [2] (labeled as ‘Fixed value 1.13’). Each field was divided into four areas, which were considered separately, namely: target (T), entrance channel (EC), distal area (D) and lateral area (L). The target area was defined as the set of points situated between the proximal 95% dose limit of the least energetic peak, and the distal 95% dose limit of the most energetic peak, for a given point in the transversal plane. EC and D areas were consequently defined as the areas proximal and distal to the target area. Finally, the lateral areas were identified as those lying outside of the transversal projection of the field. The determination of the field areas has been done automatically with the treatment plan data. Only the voxels with at least 1% of the maximum dose were considered for the histograms.

The results are very similar for the two cases considered. While the fixed value of 1.13 is between 1% and 1.5% lower than the FLUKA-calculated s values, the proposed model can reproduce them within 0.3%.

DISCUSSION AND CONCLUSIONS

It has been shown that the described fast computational method applying expression (1) is in good accordance with the Monte Carlo calculations of s for a complete carbon ion treatment field, using the treatment plan information to determine the residual ranges at every point.

The validity of the results is closely tied to the validity of equation (1). This study does not deal with the correctness of this expression, already discussed in [3], but with its possible application to 3D environments. Moreover, other simplified s expressions are available in the literature [9, 10] and the algorithm described here could be used with them as well.

The uncertainty of the stopping power ratios given by our model is directly the uncertainty of expression (1), estimated at 0.6% for R ≥ 0.3 cm, and 0.8% for R < 0.3 cm. From the results in [3] regarding variability of s among different SOBPs of different widths situated at different ranges, we do not expect our results to be sensitive to a particular choice of treatment plan. However, uncertainties have purposefully been left out of the discussion (and the figures), as the estimated error of s in the IAEA formulation [2] of 2% spans all the proposed values.

Most of the differences between the IAEA approach and the suggested model (shown in Fig. 5) do not originate from the described position dependence of s, but rather with the use of different stopping power tables, with different mean ionization potentials for air and water. With a fixed value approach, but using s = 1.145, the discrepancy with the Monte Carlo values would be significantly reduced. However, the model is still superior at reproducing the variations of s within the treatment field, as shown in Fig. 5 by the clear reduction in the separation between the different areas of the field.

Fig. 5.

Histograms showing the deviation between the FLUKA-calculated s, the s calculated with the proposed algorithm and the fixed value of 1.13, for the homogeneous cube (left) and for the head-and-neck field (right). The insets depict the division of the field into entrance channel (EC), target (T), distal (D) and lateral (L) areas, labeled on horizontal sections of the respective dose distributions.

When performing treatment plan verification with ionization chambers, a chamber-specific correction accounting for the dependence of s on the effective point of measurement could be implemented in the TPS and would result in an improvement of 0.2–0.3% of the overall accuracy of the plan verification module.

FUNDING

This work was supported by the European Community's Seventh Framework Program 2007–2013 [Grant number 215840-2, ‘PARTNER’].