Proton range monitoring using 13N peak for proton therapy applications.

The Monte Carlo method is employed in this study to simulate the proton irradiation of a water-gel phantom. Positron-emitting radionuclides such as 11C, 15O, and 13N are scored using the Particle and Heavy Ion Transport Code System Monte Carlo code package. Previously, it was reported that as a result of 16O(p,2p2n)13N nuclear reaction, whose threshold energy is relatively low (5.660 MeV), a 13N peak is formed near the actual Bragg peak. Considering the generated 13N peak, we obtain offset distance values between the 13N peak and the actual Bragg peak for various incident proton energies ranging from 45 to 250 MeV, with an energy interval of 5 MeV. The offset distances fluctuate between 1.0 and 2.0 mm. For example, the offset distances between the 13N peak and the Bragg peak are 2.0, 2.0, and 1.0 mm for incident proton energies of 80, 160, and 240 MeV, respectively. These slight fluctuations for different incident proton energies are due to the relatively stable energy-dependent cross-section data for the 16O(p,2p2n)13N nuclear reaction. Hence, we develop an open-source computer program that performs linear and non-linear interpolations of offset distance data against the incident proton energy, which further reduces the energy interval from 5 to 0.1 MeV. In addition, we perform spectral analysis to reconstruct the 13N Bragg peak, and the results are consistent with those predicted from Monte Carlo computations. Hence, the results are used to generate three-dimensional scatter plots of the 13N radionuclide distribution in the modeled phantom. The obtained results and the developed methodologies will facilitate future investigations into proton range monitoring for therapeutic applications.

Introduction

Considering the currently employed radiation therapy techniques, two types of radiotherapy modes based on photons and protons are used extensively. The ultimate goal of radiation therapy is to deliver a certain amount of radiation dose to the targeted organs while not affecting healthy organs and cells. In this regard, the use of high-energy proton beams has garnered significant attention worldwide [1] owing to their low lateral scattering, no exit dose, and high dose deposition in the Bragg peak region. Previously, we performed extensive comparisons between photon and proton radiation therapies for pediatric applications [2] and discovered that proton beams can significantly reduce the off-target dose to healthy organs and cells.

However, because of the significant gradient of the dose fall-off primarily after the Bragg peak, the proton therapy technique is sensitive to spatial uncertainties. In other words, the uncertainties in the estimated position of the tumor region can result in excessive dose deposition in non-targeted organs and reduce dose deposition in targeted organs [3,4]. These uncertainties primarily originate from approximations associated with dose calculations, unanticipated anatomical changes, and mispositioning errors during accelerator setup for irradiation. In clinical trials, a setup margin is generally allocated to the target volume to circumvent the effects of these uncertainties.

Several techniques for proton range monitoring have been introduced and discussed. Most of these methods are based on the byproduct of proton beam irradiation on patients. Other typically employed techniques include proton radiography and tomography [5,6], which primarily deliver protons of sufficient energy to the patient to reconstruct planar (two-dimensional, 2D) or tomographic (three-dimensional, 3D) images. In this transmission imaging technique, radiography images are created through the proton’s entrance and exit coordinate information provided by a sensitive detector. The primary disadvantage of proton radiography and tomography is the scattering effect, which reduces the resolution of the obtained images [7]. Another direct and cost-effective proton range monitoring technique is the ionoacoustics technique [8,9], which measures acoustic pressure waves for proton range verification. In this technique, the irradiated volume is heated as a result of the deposited radiation dose, and pressure waves are emitted consequently. The acoustic pressure waves are characterized by their amplitude, frequency, and shape, which are governed by the absorbed dose and target material. The ionoacoustics technique offers a direct approach for proton-range verification. However, for the relatively small amplitudes of acoustic signals, this task becomes more challenging. In addition, complexities associated with the coupling between acoustic sensors and human skin exist, rendering this technique laborious. In addition to the abovementioned methods, secondary electron bremsstrahlung can be used for proton range verification [10,11], which uses bremsstrahlung photons generated via charge particle deceleration in matter. Because these photons are of low energy, the method is only applicable to the irradiation monitoring of superficial tumors (i.e., shallow depths). Furthermore, the continuous energy spectrum of bremsstrahlung photons renders it difficult to detect and separate from the background radiation, unlike positron annihilation photons, which have discrete energies. Prompt gamma imaging is another widely employed method for verifying the proton range [12,13]; it uses prompt gamma rays emitted from excited nuclei during the inelastic interactions of incident protons with the target. One significant disadvantage of this approach is its low detector efficiency. Auto-activation positron emission tomography (PET) is another interesting and non-invasive technique that can be used for the range verification of protons; it focuses on measuring photons annihilated from generated positron emitters such as 11C, 15O, and 13N as a result of nuclear interaction between protons and tissues in the body of the patient. The applicability of PET imaging in proton therapy monitoring was previously investigated by several groups [14-20]. In addition, the generated positron emitting radionuclides and their production channel was listed in previous studies; this was reported for proton interaction with human tissues [21].

It is well-known that the 16O(p,2p2n)13N reaction has a relatively low threshold energy (5.660 MeV) [22]. Therefore, by computing the gradient between early and late PET scans, one can extract the 13N creation, which is discovered to be associated closely with the Bragg peak. Considering this property and the high sensitivity and spatial resolution of some previously developed PET systems, it would be useful to extensively investigate the underlying mechanism and feasibility of the 16O(p,2p2n)13N nuclear reaction and the generated 13N peak. The spatial locations of positron emitting radionuclides can be precisely measured using PET or PEM (positron emission mammography) systems around the patient after or during proton irradiation. Furthermore, the correlation between the 13N peak and the actual Bragg peak should be discussed more comprehensively for different incident proton beam energies. These studies would be useful for proton range monitoring, particularly for therapeutic applications. In the present study, we employed the Monte Carlo (MC) method to simulate the proton irradiation of a homogeneous water-gel phantom. The correlation between the 13N peak and the actual Bragg peak is discussed in terms of the offset distance for various incident proton energies. The spectral analysis (SA) approach was used to reconstruct the 13N peak, and the 3D distribution of 13N radionuclide was obtained. In addition, a standalone open-source computer program, i.e., PeakCalib, was developed to precisely calibrate the 13N peak with the actual Bragg peak. The obtained results, introduced methodology, and developed computer program would facilitate future developments in the field of proton therapy based on using the 13N peak for proton range verification. The overall objective of the present work is to investigate the effect of incident proton energy on the production of 13N positron emitting radionuclides, which in turn can be used to estimate the location of the Bragg-peak. The 13N and the Bragg-peak found to have an offset distance, and this was computed for wide range of incident proton energy. The current findings, developed tools and introduced approach would lay the pavement for future investigations and advancement in the field of proton range verification.

Materials and methods

MC method

In the present study, we used the Particle and Heavy Ion Transport Code System (PHITS) code version 3.25 [23]. It is a general-purpose MC simulation code that uses the Jet AA microscopic transport model (JAM) [24] and JAERI quantum molecular dynamics (JQMD) [25] to describe intermediate and high-energy nuclear reactions. Both the JQMD and JAM physical models can be used to describe the dynamic stages of the reactions. In the present study, a water-gel phantom measuring 10 cm × 10 cm × 40 cm was modeled. A schematic illustration of the modeled geometry is shown in Fig 1. The material composition of the modeled water-gel is presented in Table 1. The composition of water-gel phantom that was reported in previous investigations have rather very low nitrogen content [26,27]. Generally, some water-gel phantoms were produced experimentally by mixing agar powder (C14H24O9) with water (H2O), with the ratio of 1/100 (i.e., agar powder/water). Therefore, we have not considered nitrogen (14N) in the modelled water-gel. The modeled incident proton beams considered were monoenergetic pencil-like beams of energies 80, 160, and 240 MeV with protons measuring 1 cm in diameter emitted along the positive z-axis. The location where the incident proton beam is irradiated is additionally shown in Fig 1. In the modelled irradiation setup, 25 cm of air gap was considered between the proton beam and the phantom. To reduce statistical uncertainties associated with the MC method, we launched 109 protons from the modeled beam. The Monte Carlo method is well-established in simulation of radiation transport. The stochasticity of interaction of radiation with matter can be conveniently considered using the Monte Carlo method; this is mainly accomplished by using pseudo random numbers at which determines the interaction with different nuclei and sampling the angular and energy distribution. Considering such stochasticity, the statistical analysis of the results would be important, in a way that low relative error in the estimated results would be desired. More details regarding the MC simulation and modeling are available in our previous publications and the references therein [28-31].

Fig 1

Schematic diagram of water-gel phantom in three dimensions.

Table 1

Density and material composition (fraction by weight for each nuclei) used in present model.

Material	Density (g/cm3)	1H (%)	12C (%)	16O (%)
Water-gel	1.010	11.00	4.650	84.35

Upon the interaction of the protons with the target elements, different positron emitters were produced, which was primarily due to inelastic nuclear interactions. The modeled water-gel was primarily composed of oxygen (see Table 1). In addition, it is noteworthy that hydrogen does not produce stable positron emitters; therefore, it was not considered in the present discussion. In present work, we have employed a homogeneous water-gel phantom mainly to eliminate any complex geometrical effect that might arise from the heterogeneities, such as those can be found in human tissue; this is to investigate the correlation between the production of 13N radionuclides and different incident proton energies. A list of primary nuclear reactions and positron emitters generated as a result of this particular reaction is summarized in Table 2. The created isotope, half-life, reaction channel, and threshold reaction energy are listed in Table 2, these data were taken from Ref. [21].

Table 2

Created isotopes, half-life, reaction, and threshold reaction energy.

Isotope	Half-life (min)	Reaction	Threshold (MeV)
11C	20.39	12C(p,pn)11C	20.61
15O	2.037	16O(p,pn)15O	16.79
13N	9.965	16O(p,2p2n)13N[(a)]	5.660

(a): (p,2p2n) is inclusive of (p,α).

The absorbed dose vs. depth and the spatial distribution of the positron-emitting radionuclides (i.e., 11C, 15O, and 13N) were compared along the z-axis of the incident proton beam. The obtained results were normalized to the primary incident proton (see Ref. [2]). The tally results obtained from Monte Carlo simulations are mostly normalized per primary source particle by the Monte Carlo simulation package, as the absolute values would have no physical meaning. Therefore, the obtained results were normalized to the primary incident proton. Subsequently, the obtained data were converted into activity by considering the physical half-life of each positron-emitting radionuclide. The constructed dynamic frames were generated at 1 min intervals until 75 min. Our previously developed PyBLD software [32] (link: http://www.rim.cyric.tohoku.ac.jp/software/pybld/pybld.html) was used to analyze the output data from the PHITS. The obtained images were analyzed using the software, A Medical Image Data Examiner (version 1.0.4) (link: http://amide.sourceforge.net) [33]. The one-dimensional (1D) depth dose and activity profiles were obtained, which provided information regarding the location of the activity and its distribution. The 3D data were analyzed by considering four different regions of interest (ROIs): (1) whole, (2) edge, (3) plateau, and (4) Bragg-peak region. These four regions, with their respective heights, widths, and depths, are shown in Fig 2.

Fig 2

Side view of four regions of interests (ROIs) with their respective height (h), width (w), and depth (d) values.

Subsequently, the results from these four regions were used in the spectral analysis calculations for the data obtained 60 min after irradiation. In addition, 1D time profiles from the whole region data for 11C, 15O, and 13N radionuclides were calculated for 15, 20, 30, 60, and 75 min after irradiation. The presence of proton-induced radionuclides was confirmed from the obtained results. Finally, 3D scatter plots were generated by performing spectral analysis on each voxel. Additionally, the 3D distribution of 13N was obtained and visualized.

SA approach

SA is widely performed to identify kinetic components (i.e., tracers) in each voxel of a PET image in the field of nuclear medicine [34]. SA does not require non-linear optimization for compartmental modeling. More details regarding compartmental modeling are available in our previous publications and our recently developed compartmental software [35-37]. Compartmental modelling refers to the of modelling substance transport in a system consisting of multiple compartments (i.e., distinct regions/voxels), which is characterized by the transfer rates among the relevant compartments. The variations of certain substances or, more generally, the radionuclides in different compartments could be explained using sets of differential equations. SA requires relatively low computational resources; nonetheless, it can yield voxel-by-voxel functional images. Each voxel of the PET image contains several positron-emitting radionuclides because of the interaction between the incident proton beam and the target elements. In this study, we performed SA to distinguish the 13N component from other positron-emitting radionuclides in each voxel. The counts as a function of time (t) in the voxel (v) of the PET image, denoted as C(t), as a function of the incident proton beam profile of A(t), can be expressed as where ⊗ is the convolution operation; M represents different types of radionuclides produced, numbered from j = 1 to j = M; α and β are the initial radioactivity and decay constant of radionuclide j, respectively. In this study, we assumed an impulse function for A(t). Datasets for β were first prepared (by default, the range of β is from 10−4 to 0.1 s-1 and is logarithmically divided, with M = 1000), and each A(t) ⊗ exp (-βt) (impulse response function) was pre-calculated. Subsequently, sets of α were solved using a non-negative least squares estimator, which was used to solve Eq (2).

The estimated sets of α were linear; hence, SA can determine groups of α without requiring any iterations and therefore promptly calculate the sets of α and β in each voxel. Ideally, we wish to obtain a few positive sets of β that correspond to the decay constants of the produced radionuclides. However, practically, several peaks of β will appear owing to numerical errors arising from the discreteness of β. Therefore, we computed the numerical value S = ∑α for each voxel (v). S enhances the production of short half-life radionuclides (e.g., 13N with a large β) and suppresses that of longer half-life radionuclides (e.g., 11C with a small β). The threshold value of αβ was set to > 1.5, which removes the background region.

In realistic measurements using PET system, the measured signal could be weak and therefore it generates noisy images. There are various ways to circumvent the issue with weak signal and in turn denoise PET images. Recently, Guo et al. [38] introduced a novel kernel graph filtering method that could significantly tackle the issue with noisy PET images as a result of weak signal. The study performed by Guo et al. [38] was tested extensively using simulated and real life in-vivo dynamic PET datasets. The authors showed that the proposed method significantly outperforms the existing methods in sinogram denoising and image enhancement of dynamic PET at all count levels, and especially at low counts which measured signal from isotopes are weak. Therefore, the issue with weak signals that may create difficulties in realistic measurements could be solved rather effectively using denoising methods. In addition, the total body PET scanner is another system that can be used to solve the issue with weak signals; this scanner has 200 cm axial field of view (FOV) and 40 times higher sensitivity than conventional PET systems [39].

Results and discussion

MC computations

The 1D profile of dose vs. depth was obtained for energies of 80, 160, and 240 MeV. Similarly, the relative distributions of positron-emitting radionuclides (i.e., 11C, 15O, and 13N) were calculated along the z-axis (i.e., along the incident proton beam). The obtained results are shown in Fig 3, where the sum of the activities of all three radionuclides are shown in the same plot. The average relative errors were 0.169, 0.072, and 0.147 for energies 80, 160, and 240 MeV, respectively. The sum represents the combination of radionuclides possessing radio activities of 11C (T1/2 ≈ 20 min), 15O (T1/2 ≈ 2 min), and 13N (T1/2 ≈ 10 min), immediately (at time t = 0) after proton irradiation. The dose shown in Fig 3 was scaled such that it can be plotted in the same graph as the radionuclide activities. As shown in Fig 3, the production of 11C and 15O decreased before the Bragg peak and in the fallout region. However, the production of 13N generated a peak near the Bragg peak region. The primary reason causing the earlier decline of 11C and 15O as compared with 13N was that the threshold energy for the production of 11C (20.61 MeV) and 15O (16.79 MeV) was higher than that for 13N (5.660 MeV). It is noteworthy that upon the interaction of protons with matter, the proton will lose energy; therefore, lower-energy protons are to be expected in the deeper region of the water-gel phantom. At superficial depths, more high-energy protons will be present; therefore, the required threshold energy for 11C and 15O production will be satisfied. However, as the depth increases and the proton energy decreases, the dominant production reaction will be 16O(p,2p2n)13N, which has a relatively lower threshold energy. Hence, the byproduct of this reaction (i.e., 13N) is expected to be closer to the Bragg peak region. Considering the Bragg peak and the peak at which the 13N radionuclide was created, the distance offset was discovered to be 2.0, 1.9, and 2.0 mm for 80, 160, and 240 MeV, respectively. In other words, the depths at which the Bragg-peak and the peak from 13N were observed were 49.8 and 47.8 for 80 MeV, 171.8 and 169.9 mm for 160 MeV, and 345.9 and 343.9 mm for 240 MeV. The deviation or the distance offset between the Bragg peak and the 13N peak was likely due to the threshold energy for the 16O(p,2p2n)13N reaction. Based on the definition of the Bragg peak, it is clear that the dose reaches its maximum value at a depth near the end of the particle range, which implies that the incident particle energy will reach its minimum and be lower than 5.660 MeV (i.e., 13N produces the reaction threshold energy). Therefore, the peak from 13N and the actual Bragg peak would be located at different depth positions in the water-gel phantom. However, our calculations show that the offset distance was insignificant. This is similarly indicated in Fig 3(a)-3(c) for incident proton energies of 80, 160, and 240 MeV, respectively. For reference, the range of protons for three different incident energies in the water-gel based on the PHITS and SRIM is shown in Table 3 [40] (link: http://www.srim.org/). The PHITS Monte Carlo package computes the average stopping power for the charged particles and nuclei either using the ATIMA package [41].

Fig 3

1D dose and relative distribution of positron-emitting radionuclides obtained along incident beam direction, immediately after proton irradiation (i.e., t = 0) for (a) 80, (b) 160, and (c) 240 MeV incident proton energies. Statistical uncertainties associated with Monte Carlo computation is shown for sum curve.

Table 3

Proton range comparison for three different incident energies in water-gel from PHITS and SRIM [40].

Energy (MeV)	Range from PHITS	Range from SIRM	Absolute deviation
80	59 mm	58 mm	1 mm
160	199 mm	197 mm	2 mm
240	389 mm	395 mm	6 mm

The computations were performed using three different energies of 80, 160, and 240 MeV emitting along the positive z-axis from a circular source with a diameter of 1 cm. The source was used to irradiate the water-gel phantom, and the results were obtained using the PHITS MC package. Table 3 shows a comparison of the estimated proton range in the water-gel phantom based on our computations using the PHITS and the standard and widely used SRIM. The deviation between the estimated ranges is shown in Table 3. The proton ranges in the water-gel phantom estimated from the PHITS and SRIM showed good agreement. The estimated proton range between the PHITS and SRIM differed slightly owing to the different models and tabulated data used to explain proton straggling and interaction with matter. However, the deviation was relatively small compared with the overall average range for each incident proton energy. For example, considering the 240 MeV incident beam energy, the overall average range based on the PHITS and SRIM was 392 mm, and the deviation was only 1.53% of the overall average range—this can be considered negligible. In addition, the comparison between the estimated range values serves as a good benchmark for our developed MC model.

Because an offset was present between the generated 13N peak and the actual Bragg peak, a wider incident proton energy range should be considered to precisely verify the distance offset. Therefore, we used our developed model to investigate the distance offset for incident proton energies ranging from 45 to 250 MeV with an interval of 5 MeV. This energy range encompassed the most widely used proton energies used in therapeutic applications. For simplicity, they were obtained immediately after proton irradiation (t = 0). The obtained numerical results for the actual Bragg peak, 13N peak location, and distance offset (Bragg-peak location– 13N peak location) are shown in Table 4.

Table 4

Comparison between actual Bragg peak and 13N peak location in water-gel phantom with offset distance (Bragg-peak location– 13N peak location).

Energy (MeV)	Bragg peak (mm)	13N peak (mm)	Offset (mm)	Energy (MeV)	Bragg peak (mm)	13N peak (mm)	Offset (mm)
45	17.0	15.0	2.0	150	153.0	152.0	1.0
50	21.0	19.0	2.0	155	163.0	161.0	2.0
55	25.0	23.0	2.0	160	172.0	170.0	2.0
60	29.5	28.0	1.5	165	181.0	180.0	1.0
65	34.0	32.0	2.0	170	191.0	189.0	2.0
70	39.0	37.0	2.0	175	201.0	199.0	2.0
75	44.0	42.0	2.0	180	211.0	209.0	2.0
80	50.0	48.0	2.0	185	221.0	219.0	2.0
85	56.0	54.0	2.0	190	232.0	230.0	2.0
90	62.0	60.0	2.0	195	242.0	241.0	1.0
95	68.0	66.0	2.0	200	253.0	251.0	2.0
100	75.0	73.0	2.0	205	264.0	262.0	2.0
105	81.5	80.0	1.5	210	275.0	274.0	1.0
110	88.5	87.0	1.5	215	287.0	285.0	2.0
115	96.0	94.0	2.0	220	298.0	296.0	2.0
120	103.5	102.0	1.5	225	310.0	308.0	2.0
125	110.5	109.0	1.5	230	322.0	320.0	2.0
130	119.0	117.0	2.0	235	334.0	332.0	2.0
135	127.0	126.0	1.0	240	345.0	344.0	1.0
140	136.0	134.0	2.0	245	358.0	356.0	2.0
145	145.0	143.0	2.0	250	371.0	369.0	2.0

Based on the obtained results shown in Table 4 for various incident proton energies, it is clear that the offset distance ranged from 1.0 to 2.0 mm, which is within the acceptable range. The obtained data can be used for the future calibration of the measured 13N peak to the actual Bragg-peak location. It was observed that the distance offset values did not fluctuate significantly for different incident proton energies. This is primarily due to the approximately flat energy-dependent cross-section data for the 16O(p,2p2n)13N reaction in the incident proton energy range of 37.5–250 MeV. The energy-dependent cross-section data for the 16O(p,2p2n)13N nuclear reaction reported by (1) ICRU 63 [42], (2) Del Guerra [16], and (3) Litzenberg [20], these were taken from Ref. [1] and are shown in Fig 4.

Fig 4

Energy-dependent cross-section data for 16O(p,2p2n)13N nuclear reaction reported by ICRU 63 [42], Del Guerra [16], and Litzenberg [20], these were taken from Ref. [21].

The current computations were performed at intervals of 5 MeV. Considering the intermediate energies that might be used in certain irradiation facilities, we developed a standalone open-source peak calibration computer program, i.e., PeakCalib, which reports offset distance values with an energy interval of 0.1 MeV using linear and non-linear spline interpolation techniques. Details regarding the PeakCalib program and the obtained results are provided in Appendix A. The PeakCalib program is distributed and the program can be freely downloaded (from a free public repository), recompiled, and redistributed without any restrictions.

The 2D time-dependent images and their respective intensity profiles are shown in Figs 5–7 for 80, 160, and 240 MeV, respectively. They were predicted from the radioactive decay curve using the PHITS MC computer program. The obtained results were for time ranges from 15 to 75 min, which translates to a 60 min dataset. A scaling factor similar to that used to obtain the results shown in Fig 3 was used in this case, and the y-axis of the 1D profiles was labeled as the relative intensity. The relatively short-lived 15O (T1/2 = ~2.0 min) nuclei spectrum vanished almost completely after 30 min in the time-dependent profile data. The primary observable peak that was similar to the Bragg peak originated from the 13N nuclei. This trend was observed for all three simulated incident proton beam energies (i.e., 80, 160, and 240 MeV). In fact, the peak from the 13N nuclei was present for most of the simulated time values. However, owing to its relatively short half-life, the 13N peak disappeared for longer time values, such as 75 min. It is arguable that such long time durations (e.g., 75 min) will not benefit therapeutic applications; however, it is interesting to observe the presence of a 13N peak up to 60 min intervals and the dominance by the long-lived 11C radionuclides for long time durations. The creation of a 13N peak or any other positron-emitting radionuclides are affected primarily by two factors: (1) their production rate and (2) their decay rate. The two main controlling parameters are the incident proton energy and the half-life for the production and decay of these radionuclides. For example, 13N has a shorter half-life than 11C; however, it has a lower threshold energy for its creation compared with 11C. It is important to account for these two factors simultaneously when analyzing the results. Considering these two factors, for relatively long durations, only the 11C spectrum can be observed owing to its relatively long half-life. In fact, the 11C spectrum dominates regions in the phantom with proton energies equal to or higher than its nuclear reaction threshold energy. Comparing the results of different incident proton beam energies, we observed a distinct 13N peak in the time range of interest in medical imaging. Regarding the clinical significance of our study, it needs to be noted that intensity of the positron emitting radionuclides can be measured as long as they are being produced (i.e., when annihilation photon is emitted from the patient’s body). Furthermore, measuring annihilation photon provides mobility of the patient. The patient does not have to stay on the treatment couch for the measurements. The measurements can be performed in another place. In addition, by measuring annihilation photons using PET system, relative time trend would be measured rather than absolute photon counts. Therefore, it would not be necessary to start the measurements immediately after proton irradiation. Having such flexibility would in fact be beneficial in realistic clinical treatment conditions.

Fig 5

2D images and time-dependent activity of three positron-emitting nuclei for 80 MeV incident proton energy.

Fig 7

2D images and time-dependent activity of three positron-emitting nuclei for 240 MeV incident proton energy.

SA

SA was performed on the dynamic time-dependent activity results. The results shown in Fig 8(a) are those of 2D images with their respective 1D profiles obtained via SA. SA was performed for three different incident proton energies of 80, 160, and 240 MeV. The peak positions of SA-extracted 13N were consistent with the simulated 13N peak positions for the three energies presented in Figs 9–11. The production of the 13N peak via SA indicates a promising application of SA for separating the 13N peak from other positron-emitting radionuclides; this approach would be useful for analyzing the experimentally obtained data.

Fig 8

(a) SA images for 80, 160, and 240 MeV, respectively; (b) spectral analysis for different ROIs having 80 MeV incident proton beam energy.

Fig 9

3D scatter visualization and 1D profiles of SA image; 13N production and dose for 80 MeV incident proton energy.

Fig 11

3D scatter visualization and 1D profiles of SA image; 13N production and dose for 240 MeV incident proton energy.

The selected ROIs at which SA was performed are shown in Fig 2. The ROIs contained whole, edge, plateau, and peak regions (see Fig 2), and the results obtained via SA for the 80 MeV incident proton beam energy are shown in Fig 8(b). For all ROIs, the x-axis represents the half-life (i.e., log(2)/β) for the extracted radionuclides, and the y-axis represents the concentration of radionuclides, labeled as α. Based on the SA results for the whole region shown in Fig 8, it was clear that the contributions from 11C and 13N were similar. However, the contribution from the 15O radionuclides was approximately one-half those of 11C and 13N. The SA results of the edge ROI primarily comprised those of relatively long-lived radionuclides; in other words, the contribution from the long-lived radionuclides (e.g., 11C) was greater than those of 15O and 13N. Considering the plateau ROI, it was discovered that 11C offered the greatest contribution, whereas 15O and 13N indicated similar levels of contribution. Finally, the Bragg-peak ROI indicated the greatest contribution from the 13N radionuclides, whereas the contributions from 11C and 15O were negligible.

3D visualizations

For a better visualization of the Bragg peak and the peak from 13N, a 3D SA image was generated. The 13N production, dose, generated 3D SA image, and 1D profiles are shown in Figs 9–11 for incident proton beam energies of 80, 160, and 240 MeV, respectively. The 3D plots from the SA for incident proton energies of 80, 160, and 240 MeV showed a distinct creation of the Bragg peak; this is another promising approach for verifying the results, particularly those obtained experimentally. Similarly, the 3D plots for the 13N yield indicated the creation of a peak near the end of the range of the primary particles. Comparing the abovementioned two plots based on a 1D profile, it was observed that the 13N peak was created near the Bragg peak for all three different incident proton energies. In addition, the 3D dose distributions for the three different incident beam energies indicate that the dose value increased with depth in the water-gel phantom. 3D visualizations would benefit the investigation of inhomogeneous organs (those with irregular geometries) such as the lungs, head, and neck. In fact, the calculation of dose distribution for treatment planning and proton beam positioning are more complex for inhomogeneous organs. Based on the 1D profiles shown in Figs 9–11 for incident proton beams of energies 80, 160, and 240 MeV, respectively, it was observed that the prediction of the 13N peak based on SA was consistent with the computed 13N peak from MC computations (see Table 4 for the numerical values). Therefore, SA is effective for analyzing experimental data obtained from PET systems.

Conclusions

Herein, the concept of proton range monitoring using the 13N peak was discussed for various incident proton energies. The MC method using the PHITS package was used to obtain the production of positron-emitting radionuclides, namely 11C, 15O, and 13N, in the simulated water-gel phantom. Subsequently, the generated 13N peak was compared with the actual Bragg peak for various incident proton energies, i.e., those from 45–250 MeV, which is within the range of interest for therapeutic applications. The offset distance between the 13N peak and the actual Bragg peak was primarily due to the threshold energy of the 16O(p,2p2n)13N nuclear reaction. The fluctuations in the offset distance, which were relatively mild for the energy range investigated, were correlated with the energy-dependent cross-section data for the 16O(p,2p2n)13N nuclear reaction. In addition, we developed an open-source computer program to perform linear and non-linear cubic spline interpolation; the program can obtain the offset distance with an energy interval of 0.1 MeV. In addition, SA was performed to analyze the results, which indicated significant 13N production when compared with other radionuclides (11C and 15O) in the Bragg ROI. SA will benefit future experimental studies as it can separate the 13N peak from other positron-emitting radionuclides for proton range monitoring. Additionally, the obtained results and the tools developed in the present study will benefit future investigations. In future works, we aim to investigate the production of 13N and other positron emitting radionuclide by irradiating heterogenous phantoms with monoenergetic and spread-out Bragg-peak (SOBP) proton beams.

(DOCX)

Click here for additional data file.

22 Dec 2021

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: https://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols.

We look forward to receiving your revised manuscript.

Kind regards,

Mohammadreza Hadizadeh

Academic Editor

PLOS ONE

Journal Requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and

https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. Thank you for stating the following financial disclosure:

This study was supported by Agency for Medical Research and Development (AMED) under Grant Number 20he2202004h0402.

At this time, please address the following queries:

a) Please clarify the sources of funding (financial or material support) for your study. List the grants or organizations that supported your study, including funding received from your institution.

b) State what role the funders took in the study. If the funders had no role in your study, please state: “The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.”

c) If any authors received a salary from any of your funders, please state which authors and which funders.

d) If you did not receive any funding for this study, please state: “The authors received no specific funding for this work.”

Please include your amended statements within your cover letter; we will change the online submission form on your behalf.

3. Thank you for stating the following in the Acknowledgments Section of your manuscript:

[copy in statement]

We note that you have provided additional information within the Acknowledgements Section that is not currently declared in your Funding Statement. Please note that funding information should not appear in the Acknowledgments section or other areas of your manuscript. We will only publish funding information present in the Funding Statement section of the online submission form.

Please remove any funding-related text from the manuscript and let us know how you would like to update your Funding Statement. Currently, your Funding Statement reads as follows:

The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Please include your amended statements within your cover letter; we will change the online submission form on your behalf.

4. We note that you have stated that you will provide repository information for your data at acceptance. Should your manuscript be accepted for publication, we will hold it until you provide the relevant accession numbers or DOIs necessary to access your data. If you wish to make changes to your Data Availability statement, please describe these changes in your cover letter and we will update your Data Availability statement to reflect the information you provide.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: No

Reviewer #2: Partly

Reviewer #3: Yes

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: N/A

Reviewer #2: Yes

Reviewer #3: Yes

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: No

Reviewer #2: Yes

Reviewer #3: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: The authors describe a method of utilizing the auto-activation of nitrogen-13 positron emitters to monitor the location of Bragg peaks in proton therapy. The work builds upon a paper by Cho et al. published in 2017. The novel contributions of the work are to 1) develop a N-13 peak - Bragg peak offset vs proton energy function and 2) develop a method of spectral analysis to determine the N-13 component of a PET voxel signal.

With regard to contribution 1), I don’t believe a fitting function represents a significant contribution to the scientific community. Further, the precision of the data presented (whole integers) is on the same scale as the authors are attempting to fit to, which renders the fitting function obsolete.

With regard to contribution 2), the spectral analysis results appear promising, but no superior to the results presented by Cho et al. using an alternative method. Furthermore, as Cho et al. point out, N-13 peaks are only clearly present in monoenergetic Bragg peaks, which is of little practical use in proton therapy. Without performing the spectral analysis with spread-out Bragg peaks, the value of the work is questionable.

Considering the points above, I cannot recommend the work for publication. If the authors are able to demonstrate the spectral analysis method is superior to that of Cho et al. when applied to spread out Bragg peaks, the work should be resubmitted for consideration.

Reviewer #2: This is a very novel study by considering the generation of the positron emitters, especially the N13 peak proximal to the Bragg peak of proton therapy. The PHITS Monte Carlo simulation results show that the N13 peak is about 2 mm before the Bragg peak; therefore, the authors concluded that the N13 peak can be used to determine the Bragg peak location. However, I think it is still very challenging to apply the method in this study to the actual clinical treatment cases to monitor the proton ranges. My concerns are listed below:

(1) Can you elaborate the clinical significance of your study? For example, if it is a real patient being irradiated, is it practical to let the patient stay in the treatment couch waiting for 15 minutes to monitor the intensity of the positron emitters?

(2) The irradiated phantom is water-gel. I would say the reaction within the human tissues can be different from the water-gel. Have you considered the impact of the heterogeneities in human tissues on the N13 peak location?

(3) Only single energy beams have been considered in this study. Can you evaluate if this method is also applicable to using the SOBP, or even more complicated multi-field Intensity-modulated proton therapy cases?

(4) How will you accurately measure the spatial locations of these positron emitters? Using PET imagers around the patient during the proton therapy?

Reviewer #3: The paper contains interesting results in the reduction of the range uncertainties in proton therapy, and include useful tools developed that could help future research in this field, however, from my view, some parts of the work should be carefully reviewed before the publication.

General comments

In the Introduction to the work, various techniques used for the verification of the range in proton therapy collected in other works are summarized, however, references to previous papers from other authors based on the production of beta (-) emitters are not mentioned. From my perspective, these references have been used, but have not been properly referenced. In addition, it would be important to establish in a clear, explicit and summarized way, at the end of the Introduction, the main objective of the work and the methods followed.

As an example of other works with beta (-) emitters used but no mentioned in this paper would be: Simulation of Proton Therapy Treatment Verification via PET imaging of Induced Positron-Emitters, (2003), J.J. Beebe-Wang et al. The Figure 4 in your paper is exactly the same as Figure 3 of this paper.

In Materials and Methods part, however, the work does not have a continuity, and references are constantly made to other previous works by the same authors. In this way, to understand and follow the paper, it would be necessary to consult up to seven additional papers, as collected by the authors. It is understandable to refer to works and methodologies developed in previous works, but the paper should contain the information and data necessary to be fully understood without having to consult so many additional references. The paper could be organized as a compendium of some of your previous work. However, this is not the main goal stated by you.

As an example of this lack of continuity and permanent references to another works from the same authors would be:

Line 92 - More details regarding the MC simulation and modeling are available in our previous publications and the references therein [18-20].

Line 110 - the obtained results were normalized to the primary incident proton (see Ref [2]).

Line 135 - More details regarding compartmental modeling are available in our previous publications and our recently developed compartmental software [24-26].

Finally, in Conclusions, considering that in some of previous works, the use of 13N to reduce the uncertainty of the range in proton therapy have been already discussed, and it was concluded that although the signal from that isotope can be characterized by Monte Carlo techniques, however, there are difficulties in real measurements due to its low intensity, it would be important that the authors justify this fact in this work, and the real usefulness of this technique could be established in a more obvious way.

Some detailed comments

Line 3 – It is discovered. Sentence unclear. Did you discover that nuclear reaction? Please rewrite

Line 26 – Low exit dose. Sentence incorrect. Protons have not exit dose. Please rewrite

Line 63 – Recently, it was reported the reaction..Sentence incorrect. The reaction and the threshold energy is well known and stated long time ago. Please rewrite.

Line 98 - Table 1. Please justify why you choose this composition of the water-gel and why there is not Nitrogen in the phantom. In the same way, the percentage of Carbone is very high.

Line 107 - Table 2. Please include the correct reference where you have picked these data.

Line 144 - M is the number of radionuclides produced. Sentence unclear. From my view, M represents the different types of radionuclides produced, numbered from j=1 to j=M. Please rewrite.

Line 193 - From SRIM and Continuos Slow Down Approximation for protons it would be possible to obtain the results of Table 3. Please justify more widely how these data was obtained with PHITS.

Line 212 - My recommendation would be to include the average relative errors in the main paper, better than the caption of the Figure.

Line 237 - As already mentioned above, please include the correct reference.

Figure 1 - Please clarify in the paper that there is a gap of 25 cm of air from the proton beam to the phantom

In conclusion, considering that the results could be interesting in the field of range monitoring for proton treatment applications, my advice to the authors would be they should take the time to justify better and more extensively these results, and that the writing of the paper must contain references in a proper and kind way.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

13 Jan 2022

We appreciate the editor and all reviewers's comments, and we revised our manuscript accordingly. We believe the quality of the manuscript has been significantly by this updates. See the response letter for our responses.

Submitted filename: replies_20211228F.docx

Click here for additional data file.

21 Jan 2022

Proton range monitoring using 13N peak for proton therapy applications

PONE-D-21-31024R1

Dear Dr. Watabe,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Mohammadreza Hadizadeh

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #1: All comments have been addressed

Reviewer #2: All comments have been addressed

Reviewer #3: All comments have been addressed

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: N/A

Reviewer #2: Yes

Reviewer #3: (No Response)

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: (No Response)

Reviewer #2: The authors have addressed all of my concerns. Thank you!

I would recommend the acceptance of its publication.

Reviewer #3: The points made in the initial review have been considered and corrected by the authors, so the paper is recommended for publication.

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: No

7 Feb 2022

PONE-D-21-31024R1

Proton range monitoring using 13N peak for proton therapy applications

Dear Dr. Watabe:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Mohammadreza Hadizadeh

Academic Editor

PLOS ONE