Experimental comparison of clinically used ion beams for imaging applications using a range telescope

The detector consists of 61 parallel plate ionization chambers interleaved with 3 mm thick absorber plates of polymethylmethacrylat (PMMA). The combined RSP of one absorber plate and the corresponding ionization chamber is 1.192, which gives the detector the possibility to register WET differences of up to 21.5 cm. Along with its low cost and maintenance requirements, a key feature of the detector is its large active area of 30 cm × 30 cm, which is suited for clinical usage as most particle facilities have maximum delivery fields of around 20 cm × 20 cm (Haberer et al 2004, Pedroni et al 2004, Rossi 2011). The range telescope is a further development of a system initially designed to verify treatment plans (Brusasco 1999), similar to the commercial multilayer ionization chamber (MLIC) (Giraffe, IBA Dosimetry GmbH, Schwarzenbruck, Germany), and has been previously used for iRAD acquisition with carbon ions at HIT (Rinaldi et al 2014). Since then, the system has been upgraded with new readout electronics and different approaches for noise reduction to enable investigations of low-dose carbon ion iRADs (Magallanes et al 2020). To read out all 61 parallel plate ionization chambers simultaneously, an I-128 electrometer with 128 channels is used (Pyramid Technical Consultants Inc. Lexington, Ma, USA). A change in raster point position is detected by the electronics through a specially designed trigger signal that is produced from the beam delivery system. As the detector is not position-sensitive, the raster point position at the isocenter plane of the beam delivery system is used to locate the pencil beam at the detector.

For this study, two phantoms in different configurations have been irradiated with protons, helium and carbon ions under comparable imaging dose and beam settings. Beam energies of similar range in water were chosen in order to fully traverse the phantoms. The narrowest available beam spot size (Parodi et al 2012) was used to minimize range mixing effects. Additionally, the particle number per raster point was adjusted such that the physical dose delivered to the phantom was comparable for all ion irradiations (about 10 mGy per iRAD of the slab phantom). Moreover, recent works observed an average 18.9% higher relative biological effectiveness in the entrance channel for carbon ion CTs compared to proton CTs (Meyer et al 2019), suggesting a reduced effective dose difference between the here delivered ion beams. The most important irradiation parameters are summarized in table 1. All investigated ions are currently being used (protons, carbon ions) or are considered (helium ions) for clinical usage (Krämer et al 2016, Tessonnier et al 2017c). Moreover, the irradiation settings are readily available at HIT without substantial adjustment of the beam delivery or other settings.

For the iRAD acquisitions, the phantoms were manually positioned at the isocenter with the laser alignment system at HIT. Subsequently, a pRAD, heRAD and cRAD was acquired in sequence, without modifying the experimental setup. For a fair comparison between different iRADs, the same scanning irradiation pattern was used. For all acquisitions, a 1 mm image resolution, corresponding to the irradiation scanning step size, has been used. For the slab and the stepped wedge phantom, a 10 × 10 cm² and 20 × 15 cm² scanning field has been chosen, respectively.

The so-called slab phantom consists of six 1 cm thick slabs of five different tissue equivalent materials (Gammex, Middleton, WI, USA) that are combined to form a 6 × 10 × 10 cm³ phantom. Table 2 shows the experimentally validated material parameters of every slab (Hudobivnik et al 2016, Meyer et al 2017). The stepped wedge phantom consists of PMMA with an experimentally determined RSP of 1.165 (Magallanes Hernández 2017). It is carved by 12 steps of 7 × 20 mm² and one 7 × 10 mm² and has been irradiated from two orthogonal sides yielding two distinguished iRADs: the first one resulting in a simple geometrical object (square) with 13 different WET values (configuration 1) and the second one resulting in an homogeneous WET with sharp stepped edges (configuration 2).

The two phantoms were chosen because they present different challenges to our imaging method. The slab phantom consists of several sharp tissue interfaces that have to be accurately resolved and features five different materials with varying RSP values to be retrieved. On the other hand, the stepped wedge phantom has a smoother transition in WET in configuration 1 and a complicated geometry with exposed steps in configuration 2. The investigated phantoms allow us not only to (qualitatively) evaluate the spatial resolution of the different ions, but also the WET accuracy for the different steps and materials. Figure 1 shows schematics with dimensions for both phantoms and the beam direction from which the iRADs have been taken. For visualization purposes, the images have been cropped to reduce the amount of surrounding air.

Zoom In Zoom Out Reset image size

To benchmark the acquired images, MC simulations have been performed using a computational framework established at HIT on the basis of the FLUKA code (Parodi et al 2012, Tessonnier et al 2016, 2017a, Ferrari et al 2005, Böhlen et al 2014). This framework has already been used for cRAD simulations of the considered detector in (Krah et al 2015, Meyer et al 2017). For the reference simulations, 50% of the experimentally delivered particles have been used to create iRADs of comparable quality for both phantoms using FLUKA v2011.2x.0. All phantoms have been reproduced in-silico with RSP and density values that reflect the experimentally acquired values shown in table 2 (Meyer et al 2017). If necessary, the mass density and more importantly the ionization potential of the individual materials was adjusted in a trial-and-error process to match the corresponding RSP values with the experimental ones (Meyer et al 2017). To best reproduce the experimental data, the energy deposition inside the ionization chambers is scored for all 61 detector channels, as described in (Meyer et al 2017).

As the range telescope can only measure the residual range of the particle beams, the range information must be converted to WET. Since it provides better results than the previously used approach attributing the WET simply to the location of the signal maximum, the Bragg peak decomposition method (BPD) was used to identify the dominant WET value enclosed in the detector signal and to determine the pixel values in the iRADs (Meyer et al 2017) in a two-step process:

• (a)  
  The detector signal for each raster point is decomposed into its individual contributions, or channel weights, by means of a MC-based look-up table.
• (b)  
  The channel provided with the maximum weight is identified and converted to WET by an ion-specific channel-to-WET conversion, based on the previous work (Meyer et al 2017) for carbon ions and extended to helium ions and protons in this work.

For the heRAD and pRAD channel-to-WET conversion, the readout channel with maximum signal for energies in the range of 48 MeV u⁻¹ to 180 MeV u⁻¹ has been simulated in-silico and fitted with a second order polynomial function to describe the relationship between initial beam energy and maximum channel readout. The WET (in mm) as a function of initial energy and maximum readout channel for the calibration can be described by:

$$\begin{equation*}WET\left( {E,ch} \right) = rWE{T_{PMMA}} \times {t_{det}} \times \left( {\left( {{p_1} \times {E^2} + {p_2} \times E + {p_3}} \right) - ch\left( E \right)} \right),\end{equation*}$$

where $rWE{T_{PMMA}} = 1.192$ is the combined RSP of one absorber slab, ${t_{det}}$ is the physical thickness of one detector absorber slab (i.e. 3 mm), ${p_1},{p_2}$ and ${p_3}$ (table 3) are ion-dependent fitting parameters and $ch\left( E \right)$ is the channel with the highest Bragg Peak weight for an iRAD acquisition with particle beam energy $E$ [MeV u⁻¹]. The fitting values for cRADs are reported in the previous works of (Meyer et al 2017, Magallanes Hernández 2017).

Finally, the WET value is assigned to the pixel in the iRAD at the given raster point position. The above described image generation process has been applied to all reported iRADs in this paper.

The image quality of iRADs was determined by calculating the normalized root mean square error ($NRMSE$)

$$\begin{equation*}NRMSE = \frac{{\sqrt {\mathop \sum \nolimits_{x = 1,y = 1}^{{N_x},{N_y}} \frac{{{{\left( {{W_{xy}} - {\omega _{xy}}} \right)}^2}}}{{{N_x} \times {N_y}}}} }}{{{W_{max}} - {W_{min}}}},\end{equation*}$$

where ${W_{xy}}$ and ${\omega _{xy}}$ are the WET values at position (x, y) for the experimental and reference image, respectively, ${N_x}$ and ${N_y}$ are the dimension of the image in x and y and ${W_{max}}$ and ${W_{min}}$ are the maximum and minimum WET value of the reference image, respectively. The reference image was obtained from the phantom geometry and RSP estimation. For analysis, it has been manually mapped to the experimentally acquired and in-silico iRADs for each phantom. For all cases, the shift was ∼1 mm. By using the NRMSE, which is normalized with respect to the range of observed WET values, one is able to directly compare the imaging qualities of different phantoms, exhibiting different WET values.

For the slab phantom, the WET mean and standard deviation of each pixel column (i.e. vertical direction) in the iRAD have been calculated. To countermeasure effects that arise from vertical positioning uncertainties that could affect the standard deviation and mean, the iRADs were additionally cropped along the vertical axis. The mean and standard deviation are then compared to the calculated reference image. The WET accuracy of the individual slabs has been determined by the relative error (ℜ)

$$\begin{equation*}\Re = \frac{{{q_2} - {W_{Slab}}}}{{{W_{Slab}}}},\end{equation*}$$

where ${q_2}$ is the median value of the WET distribution in one slab of an iRAD and ${W_{Slab}}$ is the corresponding reference WET value. Moreover, ℜ was calculated for every slab within a region of interest defined laterally by the reference image and vertically by a margin of 10 mm added to the image border. This approach ensures that only effects of tissue interfaces and not positioning uncertainties are accounted for in the ℜ calculation. The same approach has been used for the individual steps of the stepped wedge phantom in configuration 1 without the vertical margin to the image border, as the transition between steps was smoother and allowed for an easier alignment.

Figure 2 shows the experimental iRADs (left column) for the slab phantom. Compared to the reference image (top row), the cRAD yields the highest visual agreement of all three iRADs in terms of WET (colour-coded) and geometry, followed by the heRAD and pRAD. Especially for the pRAD, the sharp transitions between the individual slabs are fading, whereas for the heRAD and the cRAD this effect is less pronounced. Overall, the qualitative difference between the heRAD and cRAD is hardly appreciable in comparison to the pRAD. This observation is also reflected in the $NRMSE$ of the iRADs being 7.8%, 10.5%, and 16.3% for the experimentally acquired cRAD, heRAD, and pRAD, respectively. The simulated iRADs (right column) are visually in good agreement with the experimental results. The main difference between experimental data and simulation is that the iRADs are less noisy at the boundaries of the individual slabs. The $NRMSE$ for the simulated iRADs is 10.7%, 12.7%, and 18.2% for the cRad, heRAD, and pRAD, respectively. The NRMSE of the experimental iRADs is within 3% to the simulated images, demonstrating a reasonable agreement between simulation and experimental data. Additional checks indicate that this difference is due to the extreme interface of the slab phantom, challenging the BPD method at the boundary to air.

Zoom In Zoom Out Reset image size

Line profiles along the horizontal dimension for the slab phantom are presented in figure 3. WET mean and standard deviation are determined by averaging along the vertical x-axis (cf Figure 2) of the iRADs. For all iRADs, the error increases at the transition between the individual slabs. Overall, the pRAD has the highest errors, followed by the heRAD and the cRAD.

Zoom In Zoom Out Reset image size

The ℜ of the individual materials of slab phantom is listed in table 4. Except for the lung slab, ℜ is always within ± 3%, demonstrating satisfactory WET resolution. The mean absolute ℜ of slabs 1–5 (muscle slab is averaged for both slabs) for the experimental iRADs is 1.8% for cRAD, 2.2% for heRAD and 2.0% for pRAD, whereas for the simulated iRADs they are 1.9%, 2.5% and 2.3%.

Configuration 1

The experimentally acquired and simulated iRADs for the stepped wedge phantom in configuration 1 are displayed in figure 4. As for the slab phantom, the cRAD yields the highest visual agreement with respect to the reference image, followed by the heRAD and pRAD. For the pRAD, the straight sides of the phantom form an almost wave-like pattern. This effect is less pronounced for the heRAD and negligible for the cRAD. The experimental $NRMSE$ is 5.2% for the cRAD, 6.0% for the heRAD and 8.6% for the pRAD. The simulated data exhibits improved results compared to the experimental data with a corresponding $NRMSE$ of 3.3%, 5.7% and 7.5 %. For both, experimental and simulated iRADs, the individual steps of the phantom and the transitions are sharply resolved. Moreover, the average absolute ℜ of the individual steps was 2.3%, 2.4% and 2.6% for experimentally acquired cRAD, heRAD and pRAD, respectively, while for in-silico iRADs the ℜ was 2.4%, 2.2% and 2.1 %. Analysing the average ℜ of the simulated iRADs, the ℜ of the step with the smallest WET is approximately 3.5 times bigger than the one of the largest WET.

Zoom In Zoom Out Reset image size

Configuration 2

Figure 5 shows the experimentally acquired and simulated iRADs for the stepped wedge phantom in configuration 2. In both cases the visual agreement with respect to the reference is worst for the pRAD. The step-like structure appears blurred, predominantly in the pRADs. The $NRMSE$ for the experimental iRADs is 8.3% for the cRAD, 8.6% for the heRAD and 13.5% for the pRAD and the corresponding errors for the simulated iRADs are 5.6%, 6.1% and 8.6%. When analysing the WET difference within the phantom, a ℜ of 0.1 %, −0.8% and −0.9% was found for experimentally acquired cRAD, heRAD and pRAD, respectively, while for in-silico generated iRADs ℜ was 0.1%, −0.8% and −0.9%.

Zoom In Zoom Out Reset image size

As all iRAD have been taken under comparable dosimetric and beam settings, as well as in a clinical quality, a fair comparison between the imaging capabilities of carbon ion, helium ion and proton beams cost-effective range telescope (very similar to nowadays commercially available dosimetric MLIC systems) is ensured. In patients, there are different types of tissues, e.g. bone and muscle, which have to be resolved by every employed imaging modality. In this regard, the slab phantom can be considered as a worst-case scenario for boundary transitions. Our analysis has shown that carbon ions achieve the best results at resolving a geometry with extreme, well-defined boundaries using the range telescope. More specifically, the visual comparison and quantitative $NRMSE$ and ℜ values of the cRAD were all closest to the reference. It is worth noting though that helium ions showed almost equal imaging capabilities compared to carbon ions. Only pRAD exhibited considerable inaccuracies at slab transitions due to their increased beam size and MCS, indicating that they might be less suited than carbon or helium ions due their increased lateral scattering, unless special algorithms are employed to counteract these effects such as (Krah et al 2015) or (Gianoli et al 2016). Nevertheless, a smoother transition between tissue interfaces as present in the slab phantom can be expected in patients. The resulting inaccuracies should therefore be less pronounced in clinical practice.

For some experimental iRADs, minor image artifacts were observed. They are caused by so-called pickup noise that originates from resonance frequencies of single ionization chambers. Although the BPD is able to eliminate some of the pick-up noise (Magallanes Hernández 2017), with our current detector setup and readout electronics the issue can not be completely avoided. Switching the here used I-128 readout electronic to a model with continuous charge integration could potentially decrease imaging artifacts (Magallanes et al 2020). In our analysis, a substantially increased ℜ was found for lung-like material. The increased ℜ is partly due to the detector's intrinsic WET resolution, i.e. due to its discrete energy measurement capabilities of up to half of the absorber plate thickness. Consequently, when calculating the ℜ normalized to the reference WET value (in the case of the lung slab only 2.84 cm), the limited WET granularity affects ℜ for all ions. A similar trend was seen when looking at the WET ℜ of the individual steps of the stepped wedge phantom in configuration 1, where the largest average ℜ was found for the step with the smallest WET. Moreover, for some inserts (e.g. adipose) cRADs yielded a higher ℜ than heRADs or pRADs, going against the otherwise overall trend of improved analysis metrics for cRADs. In these cases, we think that a one-channel detector readout shift could be the cause, as remaining WET differences between expected value and measured are smaller than 3.5 mm WET (the current WET granularity of the detector). Over one central pixel row averaged detector readouts are provided in the supplementary information for every slab and both simulation and experimental data. For visualization purposes, the averaged detector signal is normalized by its maximum.

Although the WET retrieval of the detector should not exhibit a systematic over- or underestimation, most retrieved ℜ values were smaller than 0, indicating an underestimation of the WET. It would be worth investigating, whether this effect reduces when the granularity of the detector is improved, e.g. by changing the absorber plate thickness from 3 mm to 1 mm as previously proposed (Meyer et al 2017).

The line profiles (figure 3) revealed that the overall WET resolution for all ions largely varies at the boundaries between the slabs due to the finite beam spot size and range mixing effects at tissue interfaces, while there is little error in the homogenous regions of the phantom. As already mentioned, the main source of error in the homogenous regions is the intrinsic WET resolution of the detector. This effect dominates over inaccuracies due to range straggling of the investigated ions. By decreasing the absorber plate thickness of the detector, as already investigated in-silico by (Meyer et al 2017), the overall WET error could be further reduced.

The complex geometry of the stepped wedge phantom challenges our method to resolve sharp edges in configuration 2 and different WET values in configuration 1, where the latter are due to different physical thickness rather than different RSP. Even though the observed NRSME was smaller for the stepped wedge phantom in configuration 1 than for the slab phantom (NRMSE is naturally lower at tissue interfaces with similar RSP) configuration 2 was especially challenging for protons. The individual steps of the phantoms are hardly resolved in the pRADs, resulting in the apparent triangular shape of the phantoms body. As the BPD is sensitive to larger beam shapes, the combination of an increased beam spot size and particle scattering with a finite image resolution leads to inaccuracies in the correct reproduction of sharp edges. With integration mode detectors relying on the identification of the WET component with maximum weight, the blurring of the signal translates into inaccuracy of the image. As a consequence, the iRAD does not result blurred, rather inaccurate (Gianoli et al 2019). For helium and carbon ions, this effect was less pronounced, due to their reduced beam spot size (table 1) and MCS (Tessonnier et al 2017b, Gehrke et al 2018a) compared to lighter ions. This resulted in a reduced $NRMSE$, confirming almost equal imaging capabilities for the experimental carbon and helium iRADs. To overcome the limitations due to the increased beam spot size and particle scattering with the considered integration mode detector, deconvolution kernels taking the initial beam shape into account could be investigated to potentially increase the achievable image resolution.

Except for the slab phantom, the simulated iRADs always yielded a decreased $NRMSE$ compared to experimental acquisitions since charge collection efficiency and fluctuations in the detector signals were not simulated. The simulations thereby function as a best-case approximation of the experiment. Nevertheless, overall good agreement for simulation and experimental iRADs both visually and in terms of $NRMSE$ have been found in all cases. This is attributed to the general good agreement between MC simulations and experiments at HIT, as the raster point scanning as well as the characterization of materials and phantoms have been experimentally validated in previous works (Meyer et al 2017, Magallanes Hernández 2017). Overall, the extension of the MC ion-imaging framework in combination with the BPD, originally developed in the context of carbon ion imaging, to protons and helium ions has been successfully validated. With this framework, extensive MC studies on new detector designs, including also the biological implications of ion-imaging on the human body, along with improvements of data processing approaches can be investigated in-silico prior to further experimental work (Meyer et al 2019).