Monte Carlo Simulation

Geant4 version 10.03.p01 and GATE version 9.0 with the QGSP_BERT_HP_EMY physics list were used for proton-induced PG emission and CC detection simulation. Besides, the G4EMLivermorePhysics list was used to simulate the physics processes in the production and interaction effect of PG, including the Doppler broadening effect [33]. The beam time structure was referred to the IBA cyclotron C230 used in the clinical proton therapy process [34]. The intensity of C230 clinical treatment was approximately 2× 10¹⁰ s^-1, in which the current was approximately 3.2 nA. The beam pulse duration was 3.2 ns with a period of 9.4 ns. In this case, the number of protons contained in a single pulse was 217. The approximate relationship between the number of protons transported N_p and the delivered dose is given by Eq. (1) [23]. $N_{\text{p}}=6.24\times {10}^9\frac{D}{S/\rho }{A}_{\text{r}},$ (1) where D is the delivered dose expressed in Gy(J⋅ kg^-1) and A_r is the beam area expressed in cm². For a 145-MeV proton pencil beam with a 2D Gaussian broadening of σ=5 mm, A_r was approximately 1.08 cm². Note also that S/ρ is the mass proton stopping power expressed in MeV⋅ cm²⋅ g^-1. The average energy of the proton beam at the depth of the Bragg peak was approximately 14 MeV, corresponding to a residual range of approximately 2 mm. According to the Report 90 of the International Commission on Radiation Units and Measurements (ICRU) [35], S/ρ was approximately 35.2 MeV⋅ g^-1⋅ cm^-2. When the delivered dose at the Bragg peak was 2 Gy, the number of protons transported calculated by Eq. (1) was 3.8 × 10⁸.

A non-ideal two-stage CZT CC (i.e., presenting limited position and energy resolution) with a coincidence time window of 1.5 μs was used to detect PGs. Each stage of the CZT detectors comprised 44 × 44 crystals, each with a size of 2 mm × 2 mm × 15 mm. The spatial resolutions of the detectors were both 1 mm in the lateral and depth directions. The energy resolution was set as 1.5% at 662 keV. The time resolution of the CC was set as 25 ns and the dead time of each event was approximately 250 ns. The CC recorded the deposited energies and positions of the pixel where the photons interacted in the two layers and output the time-series projection that contained a timestamp and the corresponding interaction pixel label and deposited energy. The list-mode data were obtained by selecting the time-series events in a coincidence time window in which the following requirements were met: (i) the total deposited energy was the characteristic energy of PGs; (ii) the interactions of the events were once with the scatterer first and once with the absorber. The ideal detection efficiency of the simulated CC was approximately 1 × 10^-3 PG coincidence events for one irradiated proton. However, the true coincidence yield, which would determine the detection efficiency of effective events in practice, was affected by the limited coincidence time window of the CC and dose rate of the proton beam. The dose rate could be calculated by the particle rate of protons using Eq. (1). After considering the incorrect coincidence caused by multiple protons in a bunch and different secondary particles in a coincidence time window, the true coincidence yield of the CC for different particle rates of protons was evaluated by simulation; the results are provided in Table 1.

Given that the proton beam in a bunch was too large in the case of the delivered current in the range of several nA, the probability of incorrect coincidence events when using the simulated CC detection was greater than 99.9%. Therefore, it is almost impossible to obtain correct coincidence events that can be used for reconstruction. One feasible method is to reduce the current when delivering the beam and prolong the irradiation time. When the current intensity was reduced by a factor of 10, the number of protons in a single reactor would be reduced to approximately 22 by using a delivered current of approximately 0.32 nA. For the two-layer CZT CC, the true coincidence yield was approximately 28%. After reducing the current intensity by two orders of magnitude to 0.03 nA, the number of protons in a single bunch was 2, and the true coincidence yield was approximately 90%. The particle rate delivered was 2× 10⁸ s^-1, and the time required to deliver a 2-Gy dose was 1.9 s. With this particle rate, the relationship between the delivered dose, number of protons, and irradiation time used in simulations is listed in Table 2.

Figure 1 shows the diagrams resulting from the Monte Carlo simulations. To evaluate the performance of the SD-OE-RR algorithm for PG distribution reconstructions, a 120-MeV proton beam irradiated the water box phantom three times independently. Then, a proton pencil beam with different energies irradiated the box phantom composed of different body-like materials to evaluate the proposed dose reconstruction approach. The materials (e.g., muscle, cortical bone, soft tissue) used in the simulations were referred to the ICRT-46A. To evaluate the performance of the proposed method in a situation close to the clinical proton therapy, we simulated the proton therapy in the thoracic cavity and used the proposed method to reconstruct the dose depth profiles with the convolution vectors obtained by the above box phantom simulations. The thoracic cavity phantom was referred to the Geant4 extended example DICOM [36]. Three different proton pencil beams of 2D Gaussian shape in the transversal plane irradiated the esophagus or mediastinum in different directions to simulate cancer treatment in diverse situations. For the thoracic tumor proton therapy simulation, two treatment modalities were considered. One was a single spot scan of the pencil beam with approximately 0.8 MeV energy broadening and spatial broadening parameters σ equal to 2 mm; the coverage depth of the Bragg peak maximum irradiation dose was approximately 1 mm. The other modality was the spread-out Bragg peak (SOBP) irradiation with 5.8 MeV energy broadening and spatial broadening parameter σ=5 mm for full coverage of the target area. A total of 1 × 10⁷ protons were used in the simulations to evaluate the proposed method. Finally, to investigate its performance with a different dose, various numbers of protons, ranging from 10⁷ to 10⁹, were implemented for the thoracic cavity proton therapy. The corresponding doses of the different numbers of protons delivered to the phantoms varied from 5.3 cGy to 5.3 Gy, covering the region of generally delivered dose in clinical treatments [37].

Fig. 1Full-screen View

(Color online) Diagrams resulting from the Monte Carlo simulations; (a) CC-based beam range and dose monitoring in the multi-layer box phantom; (b) non-ideal CC with limited position and energy resolution; (c) proton therapy simulation of the thoracic cavity based on the CT slice.

Subset-driven origin ensemble with resolution recovery

The subset-driven origin ensemble with resolution recovery (SD-OE-RR) algorithm was used for PG reconstruction with list-mode projection data. Different from the SD-OE-RR proposed in a previous study [38], the calculations of the resolution correction factor and initial guess of the source distribution f₀ given by Eqs. (2) and (3) were modified for CZT CC-based PG reconstruction. $\text{Δ}\left(\cos {\theta }_{\text{q}i}\right)\approx \frac{m_0{c}^2}{E_{i2}}\cdot \left({\eta }_{\text{E}i}+{\eta }_{\text{D}i}\right)+{\eta }_{\text{P}i}\mathrm{,}$ (2) $f_0={\Sigma }_j{\Sigma }_i\delta \left(v_j\mathrm{,}{e}_i\right)\left\{\left|\cos {\theta }_{\text{t}i}-\cos {\theta }_i\right|⩽\text{Δ}\left(\cos {\theta }_{\text{q}i}\right)\right\}\mathrm{,}$ (3) where $\cos {\theta }_i$ denotes the cosine of the Compton scattering angle calculated by the deposited energies for event i; $\cos {\theta }_{ti}$ is the theoretical cosine of scattering angle determined by the scattering position, absorbed position, and voxel vj in the field of view (FOV); E_i2 is the deposited energy in the absorber; ηPi is the deviation caused by the spatial resolution of the detectors; and ηEi and η_Di are the energy deviation percentages caused by the statistical fluctuation and Doppler broadening effect, respectively. For the CZT CC used in the study, the total cosine Δ (cos θ_qi) was approximately equal to 0.05.

The voxels probably containing the sources of event i were stored as a one-dimensional vector Oi=∑joij, where the subset of the origin ensemble ∑joij is given by Eq. (4). Then, the iteration of the SD-OE-RR was based on the Monte Carlo-Markov chain, updating the probability density function by Eq. (5); this is similar to the original origin ensemble algorithm [19]. The defined δ(x) function equals 0 when x ≠ 0 and equals 1 when x=0. Moreover, k+1 and k represent iteration times, β_{1, k+1} represents the randomly selected subset of the origin ensemble in the (k+1)th iteration, and β_{2, k+1} represents another independent random number used in the (k+1)th iteration to determine whether to move the location of the origin. Besides, F(β_{1, k+1}) is the sum of f(β_{1, k+1}) and f of the β_{1, k+1}-centered four adjacent pixels. Finally, sgn(x) denotes the sign function. $o_{ij}=\text{ln}\left[1+\delta \left(\max \left\{\left|\cos {\theta }_{ti}-\cos {\theta }_i\right|-\text{Δ}\left(\cos {\theta }_{\text{q}i}\right)|\mathrm{,0}\right\}\right)\right]/\text{ln}2$ (4) $\begin{array}{c}f\left({\beta }_{\mathrm{1,}k+1}\right)=\max \{\text{sgn}\left(2\left(F\left({\beta }_{\mathrm{1,}k+1}\right)+1-F\left({\beta }_{\mathrm{1,}k}\right)\right)\right)\mathrm{,}\\ \text{sgn}\left(2\left(F\left({\beta }_{\mathrm{1,}k}\right)-F\left({\beta }_{\mathrm{1,}k+1}\right)\right)\right)\cdot \\ \text{sgn}\left(\frac{F\left({\beta }_{\mathrm{1,}k+1}\right)+1}{F\left({\beta }_{\mathrm{1,}k}\right)}-{\beta }_{\mathrm{2,}k+1}\right)\cdot \left[F\left({\beta }_{\mathrm{1,}k+1}\right)+1\right]\mathrm{,}\\ \text{sgn}\left(2\left(F\left({\beta }_{\mathrm{1,}k}-F\left({\beta }_{\mathrm{1,}k+1}\right)\right)\right)\right)\cdot \text{sgn}({\beta }_{\mathrm{2,}k+1}\\ -\frac{F\left({\beta }_{\mathrm{1,}k+1}\right)+1}{F\left({\beta }_{\mathrm{1,}k}\right)})\cdot \left[F\left({\beta }_{\mathrm{1,}k}\right)+1\right]\}\end{array}$ (5) SD-OE-RR was programmed with CUDA C++ for parallel acceleration. The FOV size was set as 200 mm and 10⁷ iterations were implemented for each reconstruction.

Double evolutionary algorithm

In a previous study, the PG depth profiles (GDPs) and dose depth profiles (DDPs) could be fitted by an analytical approximation of the Bragg curve called $\tilde{Q}$ function. The convolution relation between the corresponding fitting curve $\widetilde{GDP}$ and $\widetilde{DDP}$ was given by Equation (6), which is an analytical approximation of the Bragg curve. The $\tilde{Q}$ function is the convolution of a Gaussian function G(x) and a power-law function Pv (x), which was introduced by Parodi and Bortfeld [20] for PET and adapted to PG by Schumann et al [29]. In this study, the evolutionary algorithm was used to obtain the kernel function $\tilde{k}$ of various known simple phantoms and deduce the unknown $\widetilde{DDP}$ of complex phantoms. This approach was abbreviated as double evolutionary algorithm (DEA). $\widetilde{GDP}=\widetilde{DDP}\cdot \tilde{k}$ (6) The DEA iteration process used in this study is described next:

1. Initialization. A specific array of a $\tilde{Q}$ distribution is used to create a set of N_pop individuals. Each individual represents a $\widetilde{DDP}$ or $\tilde{k}$ array for two different evolutionary algorithm applications.

2. Parent selection. Two parental arrays are obtained with a fitness proportional selection. The arrays with higher fitness give rise to new offspring.

3. Crossover. With random points cutting in two arrays, the parental arrays implement the single-point crossover to form two-child arrays.

4. Mutation. The new offspring is mutated to derive a potential improvement in fitness. Each modified array could be shifted in depth by a random integer value in the interval [-2,2] with probability p₁=0.3 and multiplied with a uniformly distributed random factor between [1-r,1+r] (r=0.2) with probability p₂=0.3. It could be also varied around a randomly chosen point with probability p₃=0.4 and a Gaussian curve of random height (between ± 10%) and width of 10 mm (2σ).

5. Replacement. The next generation is created by choosing 50% of individuals with the best fitness value from the parental and newly generated population.

6. Iteration. Repeat steps 1 to 4 until reaching the pre-determined iterations, i.e., N_iter.

The specific array mentioned above, which replaces their corresponding continuous function, represents a vector with 400 elements in the interval [-200 mm, 200 mm] for $\tilde{k}$ and a vector with 200 elements in the interval [0 mm, 200 mm for $\widetilde{DDP}$, respectively. The DEA iterations were implemented 3000 times for the kernel function $\tilde{k}$ of various known simple phantoms and 1000 times for unknown $\widetilde{DDP}$ of complex phantoms.

Dose depth profile reconstruction framework

As shown in Fig. 2, the CC-based dose depth reconstruction framework mainly comprises two parts: the $\widetilde{GDP}$ obtained by the CC and the convolution kernel $\tilde{k}$ obtained by prior known phantoms and final dose estimation via the DEA. Then, the estimated convolution kernel $\tilde{k}$ of multiple materials or multi-energy was obtained by algebraic averaging (i.e., interpolation method) or weighted average of the known kernels. The previously known kernels with proton beams of different energies irradiating different materials were calculated by the DEA. The original GDP was obtained by SD-OE-RR reconstruction with the projection data of a non-ideal CC. Moreover, the estimated $\widetilde{GDP}$ was given by local peak fitting based on the original GDP. The local fitting interval of curves used in the study was the interval above 80% of the peak value.

Fig. 2Full-screen View

Framework of dose depth profile reconstruction from CC detection-based PG by the proposed method.

PG reconstruction

Table 3 shows the PG reconstruction results for three origin ensemble (OE) algorithms (i.e., ordered OE-RR [15], ROE-RR [18], and SD-OE-RR algorithms). The PGs with four characteristic photons (i.e., from ¹²C, ¹⁴N, ¹⁵O, and ¹⁶O de-excitations) were chosen to evaluate their performance. To alleviate the effect due to the incomplete absorption and background radiation, the effective events for reconstructions were selected by using the total energy windows of coincidence events within ± 0.2 MeV of the four known PG energy spectral peaks (i.e., 4.44 MeV, 2.31 MeV, 5.25 MeV, and 6.13 MeV) [15]. As shown in Table 3, compared with previous OE algorithms, SD-OE-RR presents the same or slightly higher reconstruction accuracy. Besides, in a 64-bit Linux computer with a 2.50 GHz Intel i5-7200U CPU and a GTX 1650 Ti Nvidia GPU, for approximately 200000 events, the reconstruction times were 26, 21, and 3 seconds for the ordered OE-RR, ROE-RR, and SD-OE-RR algorithms, respectively. Therefore, the SD-OE-RR algorithm proposed in this study consumes less time while providing reconstruction whose distal-fall position deviation was less than 2 mm. The results also demonstrate that the proposed algorithm can make use of the PG distribution reconstructed by the SD-OE-RR algorithm to obtain the estimated $\widetilde{GDP}$ with good agreement in the positions of peak and distal fall-off.

Reconstruction of dose depth profiles

As shown in Figs. 3(a) and (b), compared to the exact DDP in the two-layer phantom made of cortical bone and muscle, the reconstructed DDP with CC by the proposed method had less than 0.3 mm deviation at the 50% distal fall-off position. Moreover, the relative deviation was within 2.5% in terms of absolute dose in the region above 80% maximum intensity. Besides, the reconstructed DDP had an average relative deviation of 7.7% in the region of the level area behind the dose mutation interface. As shown in Figs. 3(c) and (d), the reconstructed DDP of the three-layer phantom made of cortical bone, muscle, and soft tissue had a deviation below 0.26 mm at 50% distal fall-off position, and within 5.2% in terms of absolute dose around the Bragg peak. As shown in Figs. 3(e) and (f), for the SOBP induced by a 145~160 MeV proton beam irradiating a phantom made of cortical bone, muscle, and soft tissue, the reconstructed DDP could reproduce the peak region with an accuracy within 2 mm and predict the relative dose in the peak area with a deviation less than 4%. Furthermore, the relative deviation of the global curves (until the distal end fall-off) was approximately 4.5%.

Fig. 3Full-screen View

(Color online) Simulation results of proton beam irradiation for multi-layer box phantoms: (a)(c)(e) show the comparison of the DDPs between the exact values obtained by MC simulations and reconstructed values via the proposed method from CC for three different multi-layer phantoms consisting of cortical bone, muscle, and soft tissue, respectively; (b)(d)(f) show their corresponding relative dose deviation, respectively.

Figure 4 shows that, with respect to the exact 2D dose distributions, the Compton-based reconstructed PG distributions obtained by the SD-OE-RR algorithm were in good agreement with the positions of the distal fall-off, but the peak broadening was skewed, similar to the distributions of the initial PGs. Note also in Figs. 4(d) and (f) that the reconstructed DDPs for the in-vivo proton beam had a deviation of less than 0.6 mm for the distal fall-off position. Moreover, the reconstructed DDPs were in good agreement with the exact values at the Bragg peak in the region above 80% maximum intensity, where the relative deviation was within 5.3%. However, as shown in Figs. 4(d) and (e), the reconstruction in the region before the Bragg peak had a large deviation due to the proton beam passing through a more complex structure such as a lung. In contrast, Fig. 4(i) and (j) show that when the proton beam passes through a more homogeneous structure such as mediastinum, most of the relative deviations were within 10% in the region of the level area behind the dose mutation interface.

Fig. 4Full-screen View

(Color online) Simulation results of proton therapy for thoracic mediastinal tumors with 130.2~131 MeV proton pencil beams ((a)~(e)) and for thoracic esophageal tumors with 144.2~145 MeV proton pencil beams ((f)~(j)), respectively. (10⁷ protons). (a)(f) show the exact distributions of dose at the CT slice obtained from MC simulations; (b)(g) show the distributions of the initial PGs induced by the proton beams; (c)(h) show the reconstructed PG distributions obtained by SD-OE-RR with CC data; (d)(i) show the comparison between the exact DDPs in (a)(f) and reconstructed DDPs from (c)(h); (e)(j) show their corresponding relative deviations.

Compared to the exact initial PG distribution shown in Fig. 5(b), the Compton-based reconstructed PG distributions obtained by the SD-OE-RR algorithm shown in Fig. 5(c) were in good agreement with the distal fall-off position and spatial distribution in the area around the peak. As shown in Fig. 5(d), the reconstructed DDP for the in-vivo proton beam could predict the position of the distal fall-off with an accuracy within 0.8 mm. Moreover, as shown in Fig. 5(e), the reconstructed DDP was in good agreement with the exact value at the Bragg peak in the region above 80% maximum intensity, with a relative deviation of less than 4.8% in terms of absolute dose. However, for the region of the level area before the Bragg peak, the reconstructed DDP had larger values than the exact values with a mean relative deviation of approximately 9.2%.

Fig. 5Full-screen View

(Color online) Simulation results of proton therapy for thoracic mediastinal tumors with 131.2~137 MeV proton pencil beams (10⁷ protons) of 2D Gaussian shape in the transversal plane (σ=5 mm): (a) exact distribution of dose obtained by MC simulation; (b) corresponding initial PG distribution; (c) reconstructed PG distribution obtained by SD-OE-RR with CC data; (d) comparison between the exact DDP in (a) and the reconstructed DDP from (c); (e) corresponding relative deviations of dose.

Figure 6 shows the simulation results of proton therapy for thoracic mediastinal tumors with different numbers of incident protons. Given that the reconstruction accuracy of the PGs with the CC was influenced by the number of detected events, the DDP reconstructed by the proposed method was correlated with the number of incident protons. When the number of incident protons increased, the number of effective events available for the reconstruction of PGs also increased by almost the same proportion, under the premise that the particle rate of incident protons (which was proportional to the dose rate) remained unchanged. When the number of incident protons varied from 10⁷ to 10⁹, corresponding to a delivered dose ranging from 5.3 cGy to 5.3 Gy, the number of PG events by the simulated CZT CC ranged from approximately 2000 to approximately 200000. As shown in Fig. 6(a), the accuracy of the reconstructed GDP deteriorated as the particle number decreased, especially in the area before the Bragg peak. This is consistent with the results of previous studies on OE-RR-based PG reconstruction [15, 18]. However, the reconstructed positions of the distal fall-off were almost the same, which means that the distal fall-off was reproduced accurately. Therefore, the reconstructed DDP could provide an accurate range prediction for 10⁷ incident protons. Moreover, the accuracy of the dose prediction will be better if the number of protons increases. As shown in Figs. 6(b) and (c), for a total 1.6-Gy dose, which corresponds to 3 × 10⁸ incident protons, the reconstructed DDPs were in good agreement with the distal falling edge of the Bragg peaks with an accuracy within 0.6 mm, and the corresponding relative deviations in the region above 80% maximum intensity were less than 5.5%. When the total delivered protons reached 10⁹ (i.e., approximately 5.3 Gy), the relative deviations in the region above 80% maximum intensity were less than 4% and the average value of the global deviation was approximately 5%.

Fig. 6Full-screen View

(Color online) Simulation results of proton therapy for thoracic mediastinal tumors with different numbers of incident protons: (a) reconstructed GDPs obtained by SD-OE-RR for different particle rates; (b) relative deviations between the exact and reconstructed DDPs for different particle rates; (c) corresponding reconstructed DDPs.