Dose reconstruction with Compton camera during proton therapy via subset-driven origin ensemble and double evolutionary algorithm

Zhi-Yang Yao ·Yong-Shun Xiao ·Ji-Zhong Zhao

Abstract Compton camera-based prompt gamma (PG) imaging has been proposed for range verification during proton therapy.However, a deviation between the PG and dose distributions, as well as the difference between the reconstructed PG and exact values, limit the effectiveness of the approach in accurate range monitoring during clinical applications.The aim of the study was to realize a PG-based dose reconstruction with a Compton camera, thereby further improving the prediction accuracy of in vivo range verification and providing a novel method for beam monitoring during proton therapy.In this paper, we present an approach based on a subset-driven origin ensemble with resolution recovery and a double evolutionary algorithm to reconstruct the dose depth profile (DDP) from the gamma events obtained by a cadmium-zinc-telluride Compton camera with limited position and energy resolution.Simulations of proton pencil beams with clinical particle rate irradiating phantoms made of different materials and the CT-based thoracic phantom were used to evaluate the feasibility of the proposed method.The results show that for the monoenergetic proton pencil beam irradiating homogeneous-material box phantom,the accuracy of the reconstructed DDP was within 0.3 mm for range prediction and within 5.2% for dose prediction.In particular, for 1.6-Gy irradiation in the therapy simulation of thoracic tumors, the range deviation of the reconstructed spreadout Bragg peak was within 0.8 mm, and the relative dose deviation in the peak area was less than 7% compared to the exact values.The results demonstrate the potential and feasibility of the proposed method in future Compton-based accurate dose reconstruction and range verification during proton therapy.

Keywords Prompt gamma imaging ·Dose reconstruction ·Range verification ·Origin ensemble ·Compton camera ·Evolutionary algorithm

1 Introduction

Proton therapy has been rapidly developed and widely used in clinical cancer treatment over the past decades [1-3].The Bragg peak of the proton beam makes it possible to accurately deliver enough dose to the target tumors and reduce the radiation damage to healthy tissues.However, the uncertainties of the in vivo range and dose, which are caused by the treatment plan, patient positioning, tumor movement,etc., restrict the full clinical potential of proton therapy [4,5].Range verification is key for further improving the clinical effectiveness of proton therapy.The gamma rays derived from proton-induced excited nuclei (e.g.,12C*and16O*) are almost prompt emission (less than 10-11s); they are called prompt gamma (PG), whose distribution is highly correlated with in vivo dose distribution [6].Using PG for range verification is a feasible approach that has been proven in clinical applications [6].

A Compton camera (CC) exploits the electronic collimation principle to realize imaging of radioactive sources.It has a higher detection efficiency and abilities of multienergy reconstruction and three-dimensional imaging compared to gamma detectors with a mechanical collimator such as single-photon emission computed tomography(SPECT).Because of its unique advantages, CCs have been proposed for PG imaging during proton therapy.Several studies on CCs for PG imaging have been conducted,including on Compton system design [7-10] and optimization [11-13], as well as reconstruction algorithm optimization [14-18] and novel algorithms [19-22].Besides, range verification based on CC imaging has been experimentally verified [23, 24].However, the current range verification based on CCs can only predict the peak value and distal falling-off of the PG distribution; the prediction accuracy of the actual dose coverage area and its relative value is low.If dose reconstruction can be realized directly, it will boost the effect of this technology in practical clinical applications.However, there are few studies on Comptonbased dose prediction for proton therapy.The difficulty of this problem mainly involves two aspects: (1) the accuracy of PG reconstruction is limited by the resolution metrics of the detectors and reconstruction algorithm; (2) calculating the dose distribution from the PG distribution is a challenging problem.

Parodi and Bortfeld proposed a filtering approach to obtain the positron emission computed tomography (PET)images from dose distribution, assuming that the distribution of the positron emitter is the convolution of a pre-defined filter function and the dose distribution [25].The convolution formalism based on functions was introduced to deduce the dose distribution from the secondary radiation distribution.Besides, several modified deconvolution methods have been investigated and evaluated [26-28].Schumann et al.proposed a novel deconvolution method based on the evolutionary algorithm (EA) to deduce the dose depth distribution from the PG depth profile [29].Both deconvolution approaches require a pre-defined filter function, which depends on the projectile and target.Another approach is deep learning [30-32].Given that deep learning requires a large number of training samples, the results could be unpredictable for unknown patients, bodies and different organizations.By contrast, the deconvolution approaches could be more portable.

The goal of the study was to realize a PG-based dose reconstruction with a non-ideal CC, thereby furtherimproving the prediction accuracy of in vivo range verification and providing a novel method for beam monitoring during proton therapy.A modified subset-driven origin ensemble with resolution recovery (SD-OE-RR) is proposed to realize the PG reconstruction.We also propose a double evolutionary algorithm (DEA) to reconstruct the dose depth profile (DDP) from the reconstructed PG.We simulated proton pencil beams irradiating into the phantoms made of different materials and the CT-based thoracic phantom, respectively.A non-ideal two-layer cadmium-zinc-telluride (CZT)CC was used to detect the proton-induced PG.Finally, the proposed approach was evaluated by comparing the reconstructed dose distribution with the exact values obtained by the treatment planning system (TPS).

Table 1 Relationship between the particle rate of protons and true coincidence yield of the simulated two-layer CZT CC

Table 2 Relationship between the total number of delivered protons,dose, and irradiation time with the beam time structure used in the simulations

2 Methods

2.1 Monte Carlo simulation

Geant4 version 10.03.p01 and GATE version 9.0 with the QGSP_BERT_HP_EMY physics list were used for protoninduced 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×1010s-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 transportedNpand the delivered dose is given by Eq.(1) [23].

Fig.1 (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

whereDis the delivered dose expressed in Gy(J·kg-1) andAris the beam area expressed in cm2.For a 145-MeV proton pencil beam with a 2D Gaussian broadening ofσ=5 mm,Arwas approximately 1.08 cm2.Note also thatS∕ρis the mass proton stopping power expressed in MeV·cm2·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×108.

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 timeseries 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-3PG 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×108s-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 SDOE-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×107protons 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 107to 109, 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].

2.2 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-OERR proposed in a previous study [38], the calculations of the resolution correction factor and initial guess of the source distributionf0given by Eqs.(2) and (3) were modified for CZT CC-based PG reconstruction.

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

where cosθidenotes the cosine of the Compton scattering angle calculated by the deposited energies for eventi; cosθtiis the theoretical cosine of scattering angle determined by the scattering position, absorbed position, and voxelvjin the field of view (FOV);Ei2is the deposited energy in the absorber;ηPiis the deviation caused by the spatial resolution of the detectors; andηEiandηDiare 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 eventiwere stored as a one-dimensional vectorOi=Σjoij, where the subset of the origin ensemble Σjoijis 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 whenx≠0 and equals 1 whenx=0.Moreover, k+1 and k represent iteration times,β1,k+1represents the randomly selected subset of the origin ensemble in the (k+1)th iteration, andβ2,k+1represents 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 off(β1,k+1) andfof theβ1,k+1-centered four adjacent pixels.Finally, sgn(x) denotes the sign function.

SD-OE-RR was programmed with CUDA C++ for parallel acceleration.The FOV size was set as 200 mm, and 107iterations were implemented for each reconstruction.

2.3 Double evolutionary algorithm

The DEA iteration process used in this study is described next:

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 probabilityp1=0.3 and multiplied with a uniformly distributed random factor between [1-r,1+r] (r= 0.2) with probabilityp2=0.3.It could be also varied around a randomly chosen point with probabilityp3=0.4 and a Gaussian curve of random height (between ±10% ) and width of 10 mm (2σ).

Table 3 Simulated protoninduced PG reconstruction results of the three OE-RR algorithms for 107 incident protons

Fig.3 (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

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.,Niter.

2.4 Dose depth profile reconstruction framework

3 Results

3.1 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., from12C,14N,15O, and16O 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-OERR 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 s for the ordered OE-RR, ROERR, 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 SDOE-RR algorithm to obtain the estimatedwith good agreement in the positions of peak and distal fall-off.

3.2 Reconstruction of dose depth profiles

As shown in Fig.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 Fig.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 Fig.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%.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 Fig.4(d) and (f)that the reconstructed DDPs for the in vivo proton beamhad a deviation of less than 0.6 mm for the distal falloff 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 Fig.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) shows 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.4 (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.( 107 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

Fig.5 (Color online) Simulation results of proton therapy for thoracic mediastinal tumors with 131.2～137 MeV proton pencil beams ( 107 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

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.6 (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

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 107to 109, 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 107incident protons.Moreover, the accuracy of the dose prediction will be better if the number of protons increases.As shown in Fig.6(b) and (c), for a total 1.6-Gy dose, which corresponds to 3×108incident 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 109(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%.

4 Discussion

The difficulty of dose reconstruction based on CC comes from two aspects: the rapid and accurate PG reconstruction with non-ideal CC measured data and the calculation of dose from the reconstructed PG distribution.In this paper,we first propose a modified SD-OE-RR algorithm to realize a faster and more accurate PG reconstruction.As shown in Table 3, the SD-OE-RR algorithms provided peak estimations with an accuracy of approximately 1.0 mm for all the PGs.The accuracy of the distal fall-off positions was within 2 mm except for the 50% fall-off position of14N.Compared with the previous OE-RR algorithm, the proposed SD-OE-RR algorithm additionally considered the correction of the Doppler broadening effect calculated by the CZT extranuclear electron momentum, which further improves the speed of the convergence and accuracy of the relative intensity peak of the reconstructed PG [39].Besides,the parallel SD-OE-RR algorithm exploits the advantages of GPU multi-thread simultaneous computing and greatly reduces the reconstruction time, meeting the requirements of rapid reconstruction with a large number of events.For the proposed SD-OE-RR algorithm, in order to reduce the image blur caused by the large variability between the average states for computing the reconstructed images during OE iterations [17], we selected the imaging field of view with a number of voxels of 256 at most in each dimension while ensuring that it covered at least 20 cm in the depth direction to achieve beam monitoring.Besides, we only used the total energy window to select the effective events of PGs.However, we identified sources of background noise, such as secondary protons or neutrons generated during the irradiation, that could lead to an accidental coincidence event of wrong origin in the actual irradiation; this could hardly be culled by the energy windows [40].Better event selection methods such as those based on neural networks [41] or the distance-of-closest approach [42] can be used for PG event selection before reconstruction by the proposed method in future clinical applications.After PG reconstruction, the corresponding profile (i.e., GDP) with SD-OE-RR for16O was used to obtain theby curve local fitting with～Q.Selecting the fitting interval above 80% of the peak value was due to the good agreement between the reconstructed and exact values, as shown in a previous study for OE-RRbased PG reconstruction.

Based on the original evolutionary algorithm and～Qfunction fitting proposed in Parodi and Bortfeld [25], we propose a DEA algorithm to realize the dose calculation from the reconstructed PG distribution.The main modifications were the calculation approach of the convolution kernel and a novel convolution kernel calculation approach for complex phantoms.Instead of using the analytical method, we used a similar EA to obtain the vectors～kwith the same length of elements to replace the continuous kernel function.We implemented the protons with an energy range from 131 to 159 MeV while irradiating different phantoms.In this case, the main peak positions of these vectors remained almost unchanged when the energy of the proton beams or the density of targets did not vary significantly.Given that the SOBP generally uses an energy range within 30 MeV, the filter kernel of the corresponding proton energy in the distal falling edge or median depth can be used for dose prediction, as verified by Schumann et al.[29].The density of most human tissues is between 1.0 and 1.1 g/cm3, except for the cortical bone, which is approximately 1.82 g/cm3.Thus, the average convolution kernel could be calculated from the convolution vectors obtained by the monoenergetic proton beams irradiating the box phantoms composed of a single material.Moreover, the corresponding estimated convolution vector～kwere used as the pre-defined filter kernel in multi-layer materials and polyenergetic protons, as well as the complex thoracic phantom.In particular, because the specific material proportion of the thoracic phantom and their convolution kernel vectors were often not known in advance, its average convolution vector～kwas calculated using the convolution kernel vectors of materials with similar densities (i.e., water and cortical bone), to simulate the practical application of the proposed method.

The results shown in Figs.3, 4, 5, and 6 verified the feasibility of the proposed method in the range and dose prediction during proton therapy.The DDPs reconstructed by PGs with a non-ideal CC were in good agreement with the exact values of dose distribution calculated by TPS.In the simulation of the box phantom, the reconstructed DDP had an accuracy of less than 0.3 mm for range prediction and within 5.2% for dose prediction of monoenergetic protons irradiating.The prediction deviation of the range was less than 2 mm for a polyenergetic proton beam irradiating the multi-layer phantom (even for a large spread-out Bragg peak of approximately 2 cm).Besides, for the proton pencil beam spot scanning simulations with SOBP of approximately 3 mm with a total dose (approximately 2 Gy ～ 5 Gy) commonly used in clinical applications, the reconstructed DDPs for the in vivo proton beam had less than 0.8 mm deviation for the distal fall-off position and were in good agreement with the exact values in the region above 80% maximum intensity.The accuracy of the reconstruction depends on the number of incident protons, which determined the statistical properties of the measured projection data.The deviations of the range prediction were less than 1 mm for a number of incident protons of at least 107, which corresponds to a delivered dose of 5.3 cGy.Moreover, the accuracy of the dose prediction is better as the number of incident protons increases.

Compared with the previous range prediction methods based on CC-based PG imaging, the proposed method achieved higher accuracy in terms of range verification by directly reconstructing the dose depth distribution.This is because the distal falling edge of the PGs and that of the dose distribution usually had a deviation of approximately 2 mm or more [43].Therefore, the indirect range prediction based on the PGs had systematic errors.The proposed method simultaneously realized the dose reconstruction of the Bragg peak region of the proton beam in different materials, thereby providing a means for online monitoring and control of the beam hot spot.Compared with other offline dose verification methods, the proposed dose prediction method based on CC is promising in terms of optimization of the parameters of the proton beam and therapy plan according to the patient, s condition after a short-time feedback within several seconds, and to further improve the curative effectiveness of proton therapy.In future work, after completing the construction of the CC prototype, we will conduct range and dose prediction experiments based on CC applied on the proton beam under clinical conditions to further optimize the proposed method and apply it in future rapid beam monitoring.

5 Conclusion

In this paper, we propose an approach to realize the dose depth profile reconstruction with a CC for proton therapy.A modified SD-OE-RR algorithm and double evolutionary algorithm are presented to address the two main challenges in CC-based dose reconstruction.The simulation results of the box phantom irradiation and thoracic tumor treatment demonstrate the feasibility of the proposed method in accurate dose depth distribution reconstruction.Moreover, the simulation results also demonstrated that the algorithms proposed in this paper are feasible for more accurate prediction of beam range based on CC during proton therapy.The proposed method may be used in future rapid dose monitoring and more accurate range verification during proton therapy.

Author ContributionsAll authors contributed to the study conception and design.Material preparation, data collection and analysis were performed by Zhi-Yang Yao, Yong-Shun Xiao and Ji-Zhong Zhao.The first draft of the manuscript was written by Zhi-Yang Yao, and all authors commented on previous versions of the manuscript.All authors read and approved the final manuscript.

Data availabilityThe data that support the findings of this study are openly available in Science Data Bank at https:// www.doi.org/ 10.57760/ scien cedb.07866 and http:// resol ve.pid21.cn/ 31253.11.scien cedb.07866.