Simulated secondary particle distributions described in detail by Meric et al. [6] were used as input data for this study. The setup and distribution reconstruction are shortly summarized in the following.

The simulations were performed using the toolkit GATE v. 8.2 [7] and a publically available lung cancer patient CT [8] (see Figure 1). The introduced hybrid neutron/gamma detector was a scintillator allowing for position-resolved time and energy determination and gamma/neutron discrimination by pulse-shape discrimination. The detector is referred to as NOVCoDA (NOVO Compact Detector Array), has a volume of 30 cm³ × 20 cm³ × 20 cm³, and consists of optically segmented organic glass bars with dual-ended silicon photomultiplier light readout. A pencil beam scanning (PBS) treatment plan spot was simulated by a posterior proton beam of 85 MeV with a width of 4.7 mm (FWHM) and 2.5 mrad divergence directed centrally at the tumor.

Secondary particles generated in the patient were tracked and the phase-space information of those interacting with the detector volume were scored. Double elastic scatter events on hydrogen-1 nuclei were evaluated for neutrons and triple scatter events were evaluated for prompt gamma-rays. Background due to wrong event sequences was included, whereas random coincidences were not modelled. The applied energy cuts were 100 keV for incident neutrons and 10 keV for incident prompt gamma-rays. The time resolution, position resolution, energy resolution and segmentation effects of a realistic detector were accounted for by smearing out the phase-space distributions accordingly prior to reconstruction (see [6]; Figure 2). Event cones were then kinematically reconstructed from the time-of-flight between detector elements (only for incident neutrons), the deposited energy and the positions of interaction. Finally, the distribution of the particle origin in the patient was estimated by a simple back projection algorithm and projected onto the beam axis.

In order to evaluate the ability of the system to quantify range deviations, this process was repeated 21 times with range-shifted versions of the patient [6]. For this, range deviations of up to ±5 mm were introduced in steps of 0.5 mm in the CT by adding and removing surface tissue, respectively. For each case, the simulation was performed with 10⁹ initial protons. Realistic spot intensities (10⁶ to 10⁸ protons per spot) were achieved by sub-sampling [9] of the projected distributions (see Figure 3).

To estimate the range deviation ΔR from the projected distributions p(z), previous work [6] used the arithmetic mean of the distribution as a range landmark RL = ∑(p(z) ⋅ z)/∑p(z) and considered the difference between the RL of two distributions to be approximately equal to the introduced range deviation between those two (ΔR ≈ ΔRL). The performance of this RL was compared for each particle species individually and for the summed distribution of both particle species.

In this work, we aimed to replace this simplistic range deviation estimation by a multiple linear regression model in order to improve the estimation accuracy. As apparent from Figures 2, 3, the introduction of range deviations does not only shift the arithmetic mean (RL) of the distribution, but leads to changes in the distribution shape. Additionally, the distributions are subject to statistical noise, especially for low spot intensities (see Figure 3), which may affect the arithmetic mean more strongly than more robust measures such as distribution percentiles. Therefore, a set of 28 distribution parameters was tested in this work to estimate the range deviation more precisely, including higher statistical momenta, robust histogram features, percentiles, interquantile ranges, and the integral number of counts.

The method is described in detail in [10]. First, relevant distribution features were calculated from both the fast neutron and the prompt gamma-ray distributions based on the Image Biomarker Standardisation Initiative parameter set [11]. This parameter set contained measures of the distribution width, such as the standard deviation, interquartile ranges, coefficients of variation and mean absolute deviation, measures of position, such as the arithmetic mean and percentiles, and further parameters such as the area under the curve. The most relevant parameters were then selected by forward feature selection and the Least Absolute Shrinkage and Selection Operator (LASSO) [12]. In short, the forward selection method iteratively adds parameters based on the p-value of linear regression, while the LASSO method reduces the number of parameters by an additional penalty term in the cost function. The feature selection was performed once for each particle species individually and once for a merged parameter set containing the parameters of both particle distributions. Based on this, multivariate linear models were generated using the most predictive features as regressors to predict the applied range deviation from the reconstructed distributions.

For each investigated spot intensity, the models were trained on 40 individually re-sampled projected distributions and validated on further 20 individually re-sampled projected distributions. The root mean squared error (RMSE) between estimated and actual range deviation was used as a measure to evaluate the model performance in comparison to the previously used method based solely on RL.

The two studied feature selection methods (LASSO and forward selection) showed a high agreement in the selected parameters and attainable RMSE. The selected parameters are discussed in the following.

For neutrons, the distribution parameter “area under the curve,” i.e., the integral number of detected events per spot, was identified as the single most predictive feature by both selection algorithms for all spot intensities. Here, the model performance did not further increase when adding more parameters to the model.

For gamma rays, the area under the curve was also most often chosen by the feature selection. Here, a clear improvement of accuracy was observed for multivariate models over univariate ones. The forward selection algorithm took less parameters to converge to the final accuracy than LASSO. It converged at three to four parameters per model, without further relevant improvement when adding more parameters. This result is consistent with previous results for the prompt gamma-ray timing system [10]. Next to the area under the curve, the median and the mean absolute deviation of the distributions were chosen as predictive features.

The hybrid models combining information of the neutron and gamma-ray distributions all relied on the area under the curve of the neutron distribution as most predictive parameter, independent of the spot intensity and feature selection algorithm. The area under the curve of the gamma-ray distribution was identified as the second most important feature, confirming that the combined information of neutrons and gamma rays is advantageous. Further often-selected parameters were the median and the coefficient of variation of the prompt gamma-ray distributions. The forward-selected models converged with slightly less parameters than LASSO, again reaching their optimal reconstruction accuracy at three to four parameters.

The distribution features selected by the presented method clearly improved the range reconstruction accuracy, as compared to the previously used RL (see Figure 4). For medium-intensity spots (10⁷ initial protons), the RMSE decreased from 2.8to 1.3 mm, from 1.9 to 1.6 mm, and from 2.0 to 1.3 mm, for neutrons, gamma rays, and hybrid models, respectively, when taking into account only the feature identified as most predictive. When combining the three most predictive features, the RMSE further decreased to 1.2 mm for gamma rays only and 1.0 mm for hybrid models, and stayed at 1.3 mm for neutrons only. This corresponds to a decrease by 37%–54% relative to the previously used method for all distributions.

For all studied spot intensities (10⁶ to 10⁸ initial protons), the uncertainty of the range reconstruction RMSE was reduced by 30%–50% by the new method (see Figure 5). For the lowest (10⁶) and highest (10⁸) intensity spots, an RMSE of 1.9 and 0.4 mm was reached, respectively.

For low and medium-intensity spots, the multivariate models using the neutron signal alone were the least predictive of the range deviation, which confirms previous work [6]. The multivariate models based on the gamma-ray signal alone performed better, possibly due to a higher number of collected events per spot. The highest accuracy could be reached by hybrid models combining features of the gamma-ray and neutron distributions, which underlines the advantage of a detector system capable to evaluate both particle types.