Theoretical tumor edge detection technique using multiple Bragg peak decomposition in carbon ion therapy

Carbon ions beam lose mainly their energy through Coulomb interactions with electrons, described by the Bethe-Bloch theory. The increased energy loss per unit length as the energy decrease leads to a high energy deposition peak at the end of their range, the BP. The range of a carbon ion beam (R) can be computed using the Bragg-Kleeman rule [10]:

$$\begin{eqnarray}&&R=\alpha {E}_{0}^{p},\end{eqnarray} \tag{ 1 }$$

where α and p represent fitting parameters and E₀ is the initial beam energy. The range can also be approximated by the position of its BP. The beam deposits energy in the medium between z = 0 and R. At any point, z < R, the remaining energy, E(z), must be enough to travel the distance R − z. In this region, E(z) can be described by:

$$\begin{eqnarray}&&R-z=\alpha E{\left(z\right)}^{p}\iff E(z)={\left(\displaystyle \frac{R-z}{\alpha }\right)}^{\tfrac{1}{p}}.\end{eqnarray} \tag{ 2 }$$

As Bortfeld showed in 1999 [11], the energy deposited at any point z < R can be computed by the derivative of equation (2):

$$\begin{eqnarray}&&\displaystyle \frac{{dE}}{{dz}}(z)=\displaystyle \frac{{\left(R-z\right)}^{\tfrac{1}{p}-1}}{p{\alpha }^{\tfrac{1}{p}}}.\end{eqnarray} \tag{ 3 }$$

Assuming that λ represents an infinitesimal distance from the BP position, the energy deposited by a single carbon ion at z = R can be approximated to:

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{z\to R}\left(\displaystyle \frac{{dE}}{{dz}}\right)=\displaystyle \frac{{\lambda }^{\tfrac{1}{p}-1}}{p{\alpha }^{\tfrac{1}{p}}}.\end{eqnarray} \tag{ 4 }$$

Equations (3) and (4) allow the computation of the energy deposited/intensity at any point of a carbon ion beam crossing a homogeneous medium.

Let one consider an object composed of two materials. The first material is water and represents the surrounding medium. The second material is a cubic block of bone placed inside the surrounding material, as seen in figure 1(a). Therefore, a high gradient interface is created within the surrounding medium. This interface is irradiated with a carbon ion beam parallel to the edge direction. A range detector records the energy of the carbon ion beams at the distal end of the object (figure 1(b)). Due to the heterogeneous medium, range mixing will distort the initial pristine BP signal (expected in a homogeneous medium) into a two BPs signal. The first one corresponds to carbon ions that crossed the densest material (shallow BP at depth R₁ after crossing the medium M1) and the other one to carbon ions that crossed the least dense material (deeper BP at depth R₂ after crossing the medium M2). Thus, the detector's signal at any point is given by the sum of the deposited energy distributions of carbon ions crossing each interface:

$$\begin{eqnarray}&&\displaystyle \frac{{dE}}{{dz}}(z)={W}_{1}\left(\displaystyle \frac{{{dE}}_{M1}}{{dz}}(z)\right)+{W}_{2}\left(\displaystyle \frac{{{dE}}_{M2}}{{dz}}(z)\right),\end{eqnarray} \tag{ 5 }$$

where W₁ represents the fraction of carbon ions crossing the densest medium (M₁) and W₂ the fraction that crosses the least dense medium (M₂). Using the prior definitions (i.e. equations (3) and (4)), the intensity of the first BP, located at R₁, can be calculated as:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{dE}}{{dz}}({R}_{1}) & = & {W}_{1}\left(\displaystyle \frac{{{dE}}_{M1}}{{dz}}({R}_{1})\right)+{W}_{2}\left(\displaystyle \frac{{{dE}}_{M2}}{{dz}}({R}_{1})\right)\\ & = & {W}_{1}\left(\displaystyle \frac{{\lambda }^{\tfrac{1}{p}-1}}{p{\alpha }^{\tfrac{1}{p}}}\right)+{W}_{2}\left(\displaystyle \frac{{\left({R}_{2}-{R}_{1}\right)}^{\tfrac{1}{p}-1}}{p{\alpha }^{\tfrac{1}{p}}}\right).\end{array}\end{eqnarray} \tag{ 6 }$$

Whereas the intensity of the second BP, located at R₂ is defined as:

$$\begin{eqnarray}&&\displaystyle \frac{{dE}}{{dz}}({R}_{2})={W}_{2}\left(\displaystyle \frac{{{dE}}_{M2}}{{dz}}({R}_{2})\right)={W}_{2}\displaystyle \frac{{\lambda }^{\tfrac{1}{p}-1}}{p{\alpha }^{\tfrac{1}{p}}}.\end{eqnarray} \tag{ 7 }$$

The fraction of carbon ions that crosses each medium is defined by the distance between the edge and the central position of the irradiation beam. Considering a beam whose lateral distribution is a Gaussian distribution, the distance between the beam and the edge lateral position that would produce the measured fraction of carbon ions crossing each medium can be calculated through the integral of the 2D normal distribution (${ \mathcal N }$). The relationship is:

$$\begin{eqnarray}&&{W}_{2}=1-{W}_{1}={\int }_{Y}^{\infty }{ \mathcal N }(0,{\sigma }^{2}){dx},\end{eqnarray} \tag{ 8 }$$

where Y is the lateral distance between the beam propagation axis and the edge position. This equation can be solved for Y with the knowledge W₁ and W₂, which are extracted from equations (6) and (7).

The derived formalism can be extended to multiple interfaces and related BPs by measuring all the BP intensities and determining the fraction of carbon ions crossing each interface. Figure 2 represent a detected signal for a phantom producing N BPs. The intensity of the i^th BP, located at depth R_i, with $i\in [1,N]$ can be calculated by:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{dE}}{{dz}}({R}_{N}) & = & {W}_{N}\left(\displaystyle \frac{{\lambda }^{\tfrac{1}{p}-1}}{p{\alpha }^{\tfrac{1}{p}}}\right)\\ \displaystyle \frac{{dE}}{{dz}}({R}_{N-1}) & = & {W}_{N-1}\left(\displaystyle \frac{{\lambda }^{\tfrac{1}{p}-1}}{p{\alpha }^{\tfrac{1}{p}}}\right)\\ & & +{W}_{N}\left(\displaystyle \frac{{\left({R}_{N}-{R}_{N-1}\right)}^{\tfrac{1}{p}-1}}{p{\alpha }^{\tfrac{1}{p}}}\right)\\ \displaystyle \frac{{dE}}{{dz}}({R}_{i}) & = & {W}_{i}\left(\displaystyle \frac{{\lambda }^{\tfrac{1}{p}-1}}{p{\alpha }^{\tfrac{1}{p}}}\right)\\ & & +\displaystyle \sum _{j=i+1}^{N}{W}_{j}\left(\displaystyle \frac{{\left({R}_{j}-{R}_{i}\right)}^{\tfrac{1}{p}-1}}{p{\alpha }^{\tfrac{1}{p}}}\right).\end{array}\end{eqnarray} \tag{ 9 }$$

Zoom In Zoom Out Reset image size

For every BP, the R_i position and the energy deposited at this position (${dE}/{dz}({R}_{i})$) can be measured through the range detector. If λ, p and α are known, the fraction of carbon ions in the last BP (W_N) can be calculated. The remaining fractions (${W}_{N-1}$,...,W_i) are then iteratively calculated from the second to last ${W}_{N-1}$ BPs, up to the i − th BPs (equation (9)).

The final step of the method is to identify the edge of interested and which BPs correspond to carbon ion paths below and above the edge. To do so, prior-knowledge of the edge whereabouts is required. This information can be obtained through the planning CT which will provide an approximation of the edge position. Then, by scanning the assumed edge position with two irradiations beams at two different perpendicular distances (Y) from the edge, the exact edge position can be measured using the technique presented in figure 3. In this example, four BPs were detected in the exit detector. Thus, the incident pencil beam crossed three high-gradient edges parallel to the beam propagation axis. When moving the beam towards the edge of the tumour, W₄ decreased, while W₃, W₂, W₁ increased. By assessing the variation of the carbon ions fractions, one can determine the distance to the tumor edge Y, computed using equation (8).

Zoom In Zoom Out Reset image size

The parameters α, p and λ of equation (7) depend on the detector system and the beam accelerator and transportation system. Thereby, a prior calibration is required. To retrieve them, two carbon ion beams with different energies were propagated through a homogeneous water medium into the range detector. By applying the fit from equation (1), α and p were obtained. λ was calculated using equation (4) with the intensity of the measured BP for a known dose deposition. MC simulations with two beam energies, 300 MeV/u and 400 MeV/u, were considered to simulate such scenario and fit α, p and λ.

In order to assess and validate the derived model, two parametric phantoms with a single high-gradient edge were considered. The first parametric phantom was a water tank with a cubic bone insert of 2 cm thickness (CB₂–50% CaCO₃ [12]) in the center of the water tank (figure 1). In the second phantom, the cubic insert was replaced by a semi-cylindrical bone insert (radius = 2cm) to produce range dilution effects. Both phantoms were irradiated with carbon ion pencil beams with three different lateral spreads with full width at half maximum (FWHM) of = 4, 8 and 10 mm respectively. To investigate the precision achievable as a function of the distance for each FWHM, the irradiation beam was placed at different lateral distances from the edge (Δy = [−3, −2, −1, 0, 1, 2, 3] mm).

To present a scenario closer to the complexity of a human body, a CT dataset of a patient with a lung tumor from the The Cancer Imaging Archive [13] was used. From the CT, a high-gradient edge was identified between the tumor and the lung tissue. Two directions were selected (figure 5(a)) to crosses the edge. For each direction, five parallel irradiation beams (FWHM = 4mm and spaced 1 mm apart from each other), were propagated through the patient, covering the domain y = [−2, 2] mm around the chosen high-gradient edge.

For all simulations, a water tank (10 × 10 × 20 cm³) was placed at the exit of each phantom in order to simulate a range detector. The energy deposited in the water tank was measured every 0.4 mm along the beam propagation axis. The generated signal was analyzed and the intensities and ranges of all the BPs were measured for each considered lateral distance. The simulations were performed using the Geant4 Monte Carlo simulation code (v4.9.6.p02) [14]. Standard processes include energy loss and straggling, Multiple Coulomb Scattering using the model from [15] and elastic/inelastic ion interactions from dedicated packages of ion interactions provided in Geant4 [16]. The irradiation beam consisted of 10⁵ carbon ions, with an initial energy normally distribution with a mean of 400 MeV/u and an energy spread of 0.4 MeV/u.

This section presents the results for the first two phantoms, i.e. the cubic and the semi-cylindrical inserts. Figure 4 shows the absolute difference between the Y distance found by the algorithm and the known lateral distance, as a function of the later. Figure 4(a) shows the results for the cubic insert while figure 4(b) for the semi-cylindrical insert.

Zoom In Zoom Out Reset image size

The BPs intensities and consequently the number of carbons ions that led to each BP dropped greatly when the edge was positioned farther from the beam propagation axis. For a distance to the edge superior to the half of the FWHM of the beam ±1 mm, no multiple BPs were detected. Hence, in order to guarantee a reasonable signal/noise ratio and that multiple BP were detected a Y of FWHM/2 mm was set as the maximal distance from the propagation beam axis that could be detected.

Figure 5(a) shows a representation of two beam directions considered, as explained earlier. Figure 5(b) shows the obtained results for the beam direction Beam₁, while figure 6 shows the results obtained for Beam₂. Figure 6(a) displays the signal detected on the range detector for a single beam direction.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to determine which BPs belonged to the most/least dense medium, the change in the carbon fractions, W_i, was calculated (section 2.2). The results showed a maximal error of 1 mm. A total dose of 23.63 μGy was delivered to retrieve the edge position below a millimiter precision, no matter the direction.

This study showed a theoretical approach to confirm and detect the position of high-gradient edges. The detection of two BPs in the measured signal is limited by three factors, the detector quenching and its depth spatial resolution, and the minimum difference between two BP. To minimize the effect of the later a threshold of 3 mm was chosen as a conservative minimum difference between two BPs to be separable. Although BP identification is affected by the detector range resolution, Rinaldi et al [18] have shown a potential resolution of 0.8mm, which is below the imposed limit (3 mm). Nonetheless, these effects might affect the decomposition method and the detector response should be considered in future studies. In addition, respiratory motion effects were also not considered. In a real clinical scenario this effect should be tackled by combining the proposed method with gating strategies. Furthermore, the results obtained in this study may be limited by the two following underlying assumption: (1) the tail of secondary particles in the deposited energy distribution can be neglected and (2) carbon ions can be assumed to travel linearly through the patient. The second assumption has been investigated in Collins-Fekete et al [17]. These limitations and their impacts will be the subject of a subsequent experimental work. Future work could also consider the applicability of this method in imaging reconstruction as means to increase spatial resolution.