## 2 Materials and Methods
### The Spectral Analysis (SA) Approach  
The spectral analysis approach employs an input function and exponential response function to identify kinetic components of a tracer [47,48]. The SA approach was initially developed for PET kinetics analysis, where it indicates a time-invariant single input/single output model used for the data quantification. During proton therapy, different radioisotopes of positron emitters such as \({}^{11}\)C, \({}^{15}\)O and \({}^{13}\)N are produced by inelastic nuclear interactions of protons with the target elements. Considering that the radioactivity of each radionuclide is one component at time \(t\) and the time-activity \(A(t)\) of each voxel of the PEMGRAPH image, the concentration \(C_{\text{\tiny p}}(t)\), can be represented as  
\[C_{\text{\tiny p}}(t)=\sum_{j=1}^{M}A(t)\otimes\alpha_{j}e^{-\beta_{j}t} \tag{1}\]  
where \(\otimes\) is the convolution operation and \(A(t)\) is the time-activity input function of the system, presumed to be the impulse function in this study. The limit \(M\) represents the maximum number of terms to be included in the model and this is, in general, set to a large number, usually 1000. The values of \(\beta_{j}\) are predetermined and fixed in order to cover an appropriate spectral range. For the studies involving short-lived positron emitting isotopes, this range needs to extend to the slowest possible event up to a value appropriate to transient phenomena. Each \(A(t)\otimes exp\) (\(-\beta_{j}\)\(t\)) (impulse response function (IRF)) is pre-calculated. Subsequently, sets of \(\alpha_{j}\) values were solved using a non-negative least squares (NNLS) estimator. Once the IRFs are pre-calculated, the estimation sets of \(\alpha_{j}\) are entirely linear, and SA can instantly determine groups of \(\alpha_{j}\) without iterated computation.  
Therefore, SA can be used for large datasets. To extract the short half-life (large \(\beta=0.693/\text{T}_{1/2}\)) \({}^{15}\)O radioisotope, a range of \(\beta\) was used from 0.005 (athresh in spectral_analysis.py script) to 0.006 (betacut in spectral_analysis.py script) s\({}^{-1}\), and for the \({}^{13}\)N radioisotope, a range from 0.001 to 0.002 s\({}^{-1}\) was used, respectively. On the other hand, for the longer half-life (small \(\beta=0.693/\text{T}_{1/2}\)) radioisotopes such as \({}^{11}\)C, the range of \(\beta\) from 0.0005 to 0.0006 s\({}^{-1}\) was used. We therefore computed the numerical \(S_{v}\) values for each voxel \(v\) with the threshold of \(\alpha\beta>1.5\) to remove the background region:  
\[S_{v}=\sum_{j=1}\alpha_{j}\beta_{j} \tag{2}\]  
As a result, the SA approach can quickly calculate sets of \(\alpha\) and \(\beta\) in each voxel. Few positive sets of \(\beta\) would be obtained that match the decay constants of the produced radioisotopes.  
The SA technique was used to extract \({}^{11}\)C, \({}^{15}\)O and \({}^{13}\)N positron emitting radioisotopes These were extracted separately from the merged images of 77 frames in different regions (see Figure 4). Firstly, the SA technique was applied to 3D merged images that contained 77 frames of 40 min images (from 15 to 55 min). Secondly, the time activity was calculated from the 3D data from 40 frames with 1 min intervals over 40 min using four separate regions of interest (ROIs): (1) whole, (2) edge, (3) plateau and (4) Braggpeak region, as shown in Figure 4. Then, the time activity data from 15 to 55 min (i.e., total 40 frames) were made into a single voxel image for each ROI. Figure 4 shows these four regions, together with their corresponding heights, widths and depths, respectively. The presence of proton-induced radioisotopes for each voxel was then confirmed using the obtained results.