## Highlights

- •Protons undergo multiple Coulomb scattering, so require Bayesian path estimation.
- •Most likely path currently calculated assuming a water scattering environment.
- •We account for inhomogeneity of tissues using scattering power ratios.
- •Water assumption is valid for majority of medical scanning scenarios.

## Abstract

### Purpose

### Methods

### Results

### Conclusions

## Keywords

## 1. Introduction

*Bragg peak*in energy deposition, a phenomenon whereby most of a heavy charged particle’s energy per unit length is deposited toward the end of its trajectory in a sharp peak. This quality is advantageous in radiotherapy treatment using such particles as beam energy and intensity can be manipulated to deposit a highly conformable dose to the tumor volume, with a low dose on entry and no exit dose.

*relative stopping power*(RStP). This is simply the stopping power of a material relative to that of water at a given energy. For human tissues, the RStP is approximately independent of proton energy for therapeutic energies.

*most likely path*(MLP) formalism. In this paper we use the matrix-based MLP formula proposed by Schulte et al. [

*relative scattering power*(RScP) hereon. Using the scattering power formulae of Gottschalk [

*a priori*. Furthermore, we are able to investigate a relationship between RScP and RStP in different materials. In a 2017 study, Collins-Fekete et al. [

*TOPAS*tool [

*GEANT4*simulation toolkit [

## 2. Methods

where ${R}_{0}$ and ${R}_{1}$ are change-of-basis matrices,

and ${\mathrm{\Sigma}}_{1}$ and ${\mathrm{\Sigma}}_{2}$ are the covariance matrices,

The elements of the covariance matrices in (3), known as the

*scattering moments*, are given by (for $i=1,2$) [

where $T\left(u\right)$ is the

*scattering power*, which is defined as the rate of increase, with depth

*u*, of the mean square of the projected scattering angle $\theta $. That is,

The variance in lateral displacement is given by the (1,1) matrix element of (8);

This variance reflects the square of the standard error in the MLP and may be used to define a probability envelope surrounding the most likely path. This probability envelope is defined such that it contains at least a specified proportion of all possible MCS influenced paths. For example, given a pair of entry and exit information, only 0.3% of all possible particle trajectories should exit a $3\sigma $ envelope surrounding the MLP. Note that (8), and therefore the size of the envelope, is dependent on the material composition of the medium traversed by the particle.

### 2.1 Scattering power calculations for inhomogeneous materials

*scattering length*, $1/{X}_{s}$. For a scattering medium with mass density $\rho $, atomic number

*Z*, and mass number

*A*, the scattering length is defined as

Here, $\alpha $ is the fine structure constant, ${N}_{A}$ is Avogadro’s number,

*e*is the charge of an electron, and ${m}_{e}{c}^{2}$ is the rest mass of an electron. Using this scattering length, it can be shown that the scattering power obeys

where $E\left(u\right)$ is the depth-dependent kinetic energy of the proton and its

*reduced kinetic energy*, $\tau $, is defined as

where ${m}_{p}{c}^{2}$ is the rest mass of a proton. In addition to Gottschalk’s work [

*n*elements, each with fractional weight per volume $0<{w}_{k}\u2a7d1,k=1,\dots ,n$ then

where $\rho $ is the mass density of the composite material.

*relative scattering power*(RScP), $\widehat{T}$, as the ratio of the scattering power for the material to the scattering power for water at the same energy. Using Gottschalk’s formula (10), the energy dependent terms cancel and the result is simply the ratio of the scattering lengths. It can be shown that

with weights ${w}_{k}$ and mass density $\rho $ (in units of g/cm

^{3}) as defined in (12). The factor of 0.01636 is the reciprocal of the summation term applied to liquid water. In practice, we can now calculate the scattering power at any depth, and hence the covariance matrices in (3) using

where the subscript ‘w’ refers to the value for water.

### 2.2 Stopping power calculations for inhomogeneous materials

*u*. This can be calculated through an appropriately weighted combination of the forward and backward Euler methods, using stopping powers.

to find the energy lost per unit length if the stopping power is known. In discretizing the depth we have ${u}_{j}$ for $j=0,\dots ,N$ for some positive integer

*N*. The spacing between successive depths is constant and equal to $\delta u={u}_{j}-{u}_{j-1}$ for $j=1,\dots ,N$. The forward Euler method gives

with a boundary condition of ${E}_{0}^{F}={E}_{\mathrm{in}}$, the incoming energy of the proton as it enters the medium, assumed to be equal to the exiting energy from the accelerator. The backward Euler method gives

with a boundary condition of ${E}_{N}^{F}={E}_{\mathrm{out}}$, the outgoing energy of the proton as it exits the medium, as measured by a detector.

with boundary conditions ${E}_{0}={E}_{\mathrm{in}}$ and ${E}_{N}={E}_{\mathrm{out}}$.

where ${\rho}_{e}$ is the electron number density of the material,

*I*is the

*mean excitation potential*, ${\u220a}_{0}$ is the vacuum permittivity, and $\beta \equiv v/c$ where

*v*is the speed of the proton and

*c*is the speed of light in a vacuum. Experimental values of the mean excitation potential for a single-element substance (e.g. O

_{2}gas) can be obtained from a database, but for more complicated composite materials consisting of

*n*elements, one may use the Bragg additive rule,

where

${w}_{k}$ is the fractional weight per volume of the

*k*-th element, as in (12).

*relative stopping power*(RStP), $\widehat{S}$, as the stopping power in the material,

*S*, divided by the stopping power in water, ${S}_{\mathrm{w}}$, at the same energy. That is,

Note that the RStP is weakly energy dependent, however the effect for most materials and energies encountered clinically can be assumed negligible. This is explored in more detail in Section 2.4. Now, if the material through which the proton is passing is known, then one may calculate the stopping power by simply multiplying the stopping power in water by the RStP for that material,

In the proposed use of this method, the current estimate of RStP from the reconstructed pCT image will be used to approximate RStP in each voxel. The forward and backward Euler methods can now be alternatively written, respectively, as

where ${\widehat{S}}_{j}$ is the relative stopping power at depth ${u}_{j}$.

### 2.3 MLP-Spline-Hybrid method

### 2.4 Catalogue of anatomical materials

^{3}has been used as the electron number density.) Mean excitation values used in the calculation of stopping power were obtained, for the same set of materials, from ICRU Report 49 [

*GEANT4*RStP values at various proton energies to that at 300 MeV [

*Gammex 467*tissue characterisation phantom, for which the material contributing to the greatest RStP variation (approximately 1%) was the cortical bone equivalent,

*Gammex SB3*. Our calculations show that our nominated RStP value does not differ from the true value by more than 1% for 30 to 300 MeV protons, with the largest error manifested in cortical bone. Further, the error is less than 2% over the entire range of 10 to 300 MeV, and the residual energy of a proton exiting a patient is generally larger than both of these lower limits.

### 2.5 Correlation of relative scattering and stopping powers in human tissues

*calibration curve*that determines RScP from RStP could provide an improvement to the convergence of the algorithm and the final accuracy of the image.

where ${p}_{1},{p}_{2}\in \mathbb{R}$ are constant parameters, given in Table 1, and ${p}_{2}$ is constrained to zero for low density materials.

${\widehat{p}}_{1}$ | ${\widehat{p}}_{2}$ | ${r}^{2}$ | ${N}_{\mathrm{s}}$ | |
---|---|---|---|---|

Solid tissues | 0.3905 | 0.6448 | 0.9945 | 45 |

Low density materials | 1.2127 | 0 | 0.9958 | 4 |

### 2.6 Monte Carlo simulations

*TOPAS*, version 3.1.2 using the standard environment and physics lists [

*t*-

*u*and

*u*-

*v*planes (refer to Fig. 2) differed from their mean values by more than a chosen number of standard deviations. Once the data cuts were performed, the MLP was only calculated for the projection of the path onto the

*t*-

*u*plane, as the scattering in orthogonal directions form two independent probabilistic processes [

*u*, the material composition of a realistic phantom is voxelised and must be treated as such in the implementation of ${\mathrm{MLP}}_{\mathrm{x}}$ and ${\mathrm{MLP}}_{\mathrm{x}}$SH. For each proton registered at the prior and posterior detectors, the entry and exit positions and directions were used to fit a three-dimensional cubic spline. This approximate path was subsequently used to interpolate the RStP in the traversed voxel at depth

*u*. The RStP was pre-calculated for each CT coordinate using the HU-RStP calibration of Schneider et al. [

## 3. Results

### 3.1 Water phantom

tracks outside 3$\sigma $ envelope (%) | ||||
---|---|---|---|---|

${\mathit{MLP}}_{\mathrm{H}2\mathrm{O}}$ | ${\mathit{MLP}}_{\mathrm{x}}$ | ${\mathit{MLP}}_{\mathrm{x}}$SH | ||

Water phantom (200 MeV) | No data cuts | 2.47 | 2.44 | – |

3$\sigma $ cuts | 0.444 | 0.419 | – | |

2$\sigma $ cuts | 0.251 | 0.214 | – | |

Slab phantom | 2$\sigma $ cuts (210 MeV) | 9.63 | 0.223 | 0.208 |

2$\sigma $ cuts (215 MeV) | 3.91 | 0.264 | 0.289 | |

2$\sigma $ cuts (220 MeV) | 2.34 | 0.183 | 0.204 | |

2$\sigma $ cuts (225 MeV) | 1.75 | 0.250 | 0.254 | |

2$\sigma $ cuts (230 MeV) | 1.42 | 0.215 | 0.215 |

### 3.2 Slab phantom

*MATLAB*(MATLAB R2017b, The MathWorks, Inc., Natick, Massachusetts, United States) on a standard

*Intel Core i7*processor. The spline-hybrid method, ${\mathrm{MLP}}_{\mathrm{x}}$SH, regularly achieved very similar results up to 3 times faster on average. It can be seen in Fig. 5 that at lower energies, at which scattering is more pronounced, ${\mathrm{MLP}}_{\mathrm{x}}$SH retains greater accuracy than ${\mathrm{MLP}}_{\mathrm{H}2\mathrm{O}}$. After applying 2$\sigma $ statistical cuts, the percentage of tracks falling outside the 3$\sigma $ envelope, listed in Table 2, was consitently less than 0.3% using the inhomogeneous formalism, indicating that the probabilistic path estimation is suitably accurate. On the contrary, significantly more tracks escaped the ${\mathrm{MLP}}_{\mathrm{H}2\mathrm{O}}$ envelope, with only 90.4% being contained for the 210 MeV beam. An example proton track is shown in Fig. 6(c) along with all three MLP estimates. It is clear that ${\mathrm{MLP}}_{\mathrm{x}}$ and ${\mathrm{MLP}}_{\mathrm{x}}$SH give a larger and more right-skewed probability envelope surrounding the proton path when compared to ${\mathrm{MLP}}_{\mathrm{H}2\mathrm{O}}$ for this particular track.

### 3.3 Human head phantom

## 4. Discussion

*TOPAS*[

*u*. Namely,

We have adopted the subscript “dM” (for “differential Molière”) from Gottschalk’s work. Gottschalk showed that this nonlocal formula outperformed the differential Highland formula presented by Kanematsu [

- Khellaf F.
- Krah N.
- Rinaldi I.
- Létang J.M.
- Rit S.

*Phys Med Biol.*2019; 64065003https://doi.org/10.1088/1361-6560/ab02a8

## 5. Conclusions

## Declaration of Competing Interest

## Acknowledgments

## Appendix A. Implementation of inhomogeneous most-likely-path algorithm

## References

- The calibration of CT Hounsfield units for radiotherapy treatment planning.
*Phys Med Biol.*1996; 41: 111https://doi.org/10.1088/0031-9155/41/1/009 - Vision 20/20: Proton therapy.
*Med Phys.*2009; 36: 556-568https://doi.org/10.1118/1.3058485 - Effects of Hounsfield number conversion on CT based proton Monte Carlo dose calculations.
*Med Phys.*2007; 34: 1439-1449https://doi.org/10.1118/1.2715481 - The application of protons to computed tomography.
*IEEE Trans Nucl Sci.*1978; 25: 657-660https://doi.org/10.1109/TNS.1978.4329389 - Proton Computed Tomography.
*IEEE Trans Nucl Sci.*1979; 26: 1635-1640https://doi.org/10.1109/TNS.1979.4330455 - Computed tomography using proton energy loss.
*Phys Med Biol.*1981; 26: 965-983https://doi.org/10.1088/0031-9155/26/6/001 - Conceptual design of a proton computed tomography system for applications in proton radiation therapy.
*IEEE Trans. Nucl. Sci.*2004; 51: 866-872https://doi.org/10.1109/TNS.2004.829392 - Multiple Coulomb scattering and spatial resolution in proton radiography.
*Med Phys.*1994; 21: 1657-1663https://doi.org/10.1118/1.597212 - A maximum likelihood proton path formalism for application in proton computed tomography.
*Med Phys.*2008; 35: 4849-4856https://doi.org/10.1118/1.2986139 - The most likely path of an energetic charged particle through a uniform medium.
*Phys Med Biol.*2004; 49: 2899https://doi.org/10.1088/0031-9155/49/13/010 - On the use of a proton path probability map for proton computed tomography reconstruction.
*Med Phys.*2010; 37: 4138-4145https://doi.org/10.1118/1.3453767 - The Effect of Tissue Inhomogeneities on the Accuracy of Proton Path Reconstruction for Proton Computed Tomography.in: AIP Conference Proceedings. vol. 1099. AIP, 2009: 476-480https://doi.org/10.1063/1.3120078
- On the scattering power of radiotherapy protons.
*Med Phys.*2010; 37: 352-367https://doi.org/10.1118/1.3264177 - Multiple scattering with energy loss.
*Phys Rev.*1948; 74: 1534-1535https://doi.org/10.1103/PhysRev. 74.1534 - Extension of the Fermi-Eyges most-likely path in heterogeneous medium with prior knowledge information.
*Phys Med Biol.*2017; 62: 9207-9219https://doi.org/10.1088/1361-6560/aa955d - Alternative scattering power for Gaussian beam model of heavy charged particles.
*Nucl Instrum Methods Phys Res B.*2008; 266: 5056-5062https://doi.org/10.1016/j.nimb.2008.09.004 - Reconstruction for proton computed tomography by tracing proton trajectories: a Monte Carlo study.
*Med Phys.*2006; 33: 699-706https://doi.org/10.1118/1.2171507 - Bragg peak prediction from quantitative proton computed tomography using different path estimates.
*Phys Med Biol.*2011; 56: 587https://doi.org/10.1088/0031-9155/56/3/005 - Developing a phenomenological model of the proton trajectory within a heterogeneous medium required for proton imaging.
*Phys Med Biol.*2015; 60: 5071https://doi.org/10.1088/0031-9155/60/13/5071 - TOPAS: An innovative proton Monte Carlo platform for research and clinical applications.
*Med Phys.*2012; 39: 6818https://doi.org/10.1118/1.4758060 - Geant4 - a simulation toolkit.
*Nucl Instrum Methods Phys Res A.*2003; 506: 250-303https://doi.org/10.1016/S0168-9002(03)01368-8 - Approximations to multiple Coulomb scattering.
*Nucl Instrum Methods Phys Res B.*1991; 58: 6-10https://doi.org/10.1016/0168-583X(91)95671-Y - Multiple Scattering of Electrons.
*Phys Rev.*1940; 57: 24-29https://doi.org/10.1103/PhysRev. 57.24 - Multiple scattering of electrons. II.
*Phys Rev.*1940; 58: 36-42https://doi.org/10.1103/PhysRev. 58.36 - The Scattering of Fast Electrons by Atomic Nuclei.
*Proc Royal Soc A.*1929; 124: 425-442 - LXXIX. The scattering of
*α*and*η*particles by matter and the structure of the atom.*Philos. Mag.*1911; 21: 669-688https://doi.org/10.1080/14786440508637080 - High-energy particles. Prentice-Hall physics.Prentice-Hall, New York, NY1952: 63-69
Bethe H, Ashkin J. Passage of Radiations through Matter. In: Experimental Nuclear Physics. 1953;ed: E. Segre:253.

White DR, Griffith RV, Wilson IJ. Report 46. J. ICRU. 1992;os24(1):NP. doi: 10.1093/jicru/os24.1.Report46.

Berger MJ, Inokuti M, Andersen HH, Bichsel H, Powers D, Seltzer SM, et al. Report 49. J. ICRU. 1993;os25(2):NP. doi: 10.1093/jicru/os25.2.Report49.

- Monte Carlo comparison of x-ray and proton CT for range calculations of proton therapy beams.
*Phys Med Biol.*2015; 60: 7585-7599https://doi.org/10.1088/0031-9155/60/19/7585 - Physics settings for using the Geant4 toolkit in proton therapy.
*IEEE Trans Nucl Sci.*2008; 55: 1018-1025https://doi.org/10.1109/TNS.2008.922816 - The UF family of reference hybrid phantoms for computational radiation dosimetry.
*Phys Med Biol.*2010; 55: 339-363https://doi.org/10.1088/0031-9155/55/2/002 - Reconstruction of organ dose for external radiotherapy patients in retrospective epidemiologic studies.
*Phys Med Biol.*2015; 60: 2309-2324https://doi.org/10.1088/0031-9155/60/6/2309 Schultze B, Witt M, Censor Y, Schulte RW, Schubert KE. Performance of Hull-Detection Algorithms For Proton Computed Tomography Reconstruction. 2014. arXiv:1402.1720.

- Density resolution of proton computed tomography.
*Med Phys.*2005; 32: 1035-1046https://doi.org/10.1118/1.1884906 - Development of proton computed tomography detectors for applications in hadron therapy.
*Nucl Instrum Methods Phys Res A.*2016; 809: 120-129https://doi.org/10.1016/j.nima.2015.07.066 - PRaVDA: The first solid-state system for proton computed tomography.
*Phys Med.*2018; 55: 149-154https://doi.org/10.1016/j.ejmp.2018.10.020 - Expected proton signal sizes in the PRaVDA Range telescope for proton computed tomography.
*J Instrum.*2015; 10: P05013https://doi.org/10.1088/1748-0221/10/05/P05013 - Effects of transverse heterogeneities on the most likely path of protons.
*Phys Med Biol.*2019; 64065003https://doi.org/10.1088/1361-6560/ab02a8

## Article info

### Publication history

### Identification

### Copyright

### User license

Creative Commons Attribution (CC BY 4.0) |## Permitted

- Read, print & download
- Redistribute or republish the final article
- Text & data mine
- Translate the article
- Reuse portions or extracts from the article in other works
- Sell or re-use for commercial purposes

Elsevier's open access license policy