{'fraction_non_alphanumeric': 0.03235494589758882, 'fraction_numerical': 0.01723183096191531, 'mean_word_length': 4.631402183039462, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dictionary-based protoacoustic dose map imaging for proton range verification

Proton radiotherapy has the potential to provide state-of-the-art dose conformality in the tumor area, reducing possible adverse effects on surrounding organs at risk. However, uncertainties in the exact location of the proton Bragg peak inside the patient prevent this technique from achieving full clinical potential. In this context, in vivo verification of the range of protons in patients is key to reduce uncertainty margins. Protoacoustic range verification employs acoustic pressure waves generated by protons due to the radio-induced thermoacoustic effect to reconstruct the dose deposited in a patient during proton therapy. In this paper, we propose to use the a priori knowledge of the shape of the proton dose distribution to create a dictionary with the expected ultrasonic signals at predetermined detector locations. Using this dictionary, the reconstruction of deposited dose is performed by matching pre-calculated dictionary acoustic signals with data acquired online during treatment. The dictionary method was evaluated on a single-field proton plan for a prostate cancer patient. Dose calculation was performed with the open-source treatment planning system matRad, while acoustic wave propagation was carried out with k-Wave. We studied the ability of the proposed dictionary method to detect range variations caused by anatomical changes in tissue density, and alterations of lateral and longitudinal beam position. Our results show that the dictionary-based protoacoustic method was able to identify the changes in range originated by all the alterations introduced, with an average accuracy of 1.4 mm. This procedure could be used for in vivo verification, comparing the measured signals with the precalculated dictionary.


Introduction
Proton therapy employs proton beams for cancer treatment. This technique is potentially advantageous compared to conventional photon radiotherapy since protons deposit most of their energy at the end of their path, in the so-called Bragg peak region. Besides, protons have a limited range and therefore no dose is deposited beyond the Bragg peak [1]. The depth reached by proton beams is determined by their energy [2], which varies from 50 to 250 MeV in therapeutic situations [3]. These energies give proton ranges in water between 22 and 380 mm according to the NIST PSTAR database [4]. On the other hand, photons do not have a well defined range, instead photon beams decay in intensity exponentially somewhere after entry in the patient. Therefore, the input photon beam energies that are suitable for cancer treatment are different from the ones used in proton therapy (the appropriate energy mega-voltages of photon beams vary from 2 to 18 MV [5]). Consequently, proton therapy offers the possibility to achieve more conformal dose distributions, thus reducing the radiation received by healthy tissues, especially those distal to the target. The first technique developed for proton beam delivery in radiotherapy is passive scattering. In this approach, the beam is spread using scatter foils, and then is shaped to the target volume with collimators and compensators [6,7]. On the other hand, the most modern mode for the delivery of proton beams is active scanning. In this technique, thousands of narrow and quasi-monoenergetic proton 'pencil' beams are magnetically scanned to cover voxel by voxel each of the energy layers in which the tumor is divided [6,7]. All these characteristics make proton therapy especially suitable for the treatment of pediatric, head and neck, liver, pancreatic and prostate cancers [2].
Due to the sharp fall-off of the dose after the Bragg peak (Fig. 1a), in proton therapy it is crucial to know with great accuracy the position in which the Bragg peak is located inside the patient. However, proton range calculation is affected by several factors, such as patient positioning errors, organ motion due to respiratory and cardiac movement [8], interfractional variations in patient's anatomy [6], or uncertainties in obtaining the proton stopping power maps from photon-based CT scans [9][10][11]. Currently, clinical plans extend the irradiated target volume by a variable safety margin to ensure it is fully covered [8]. For instance, some US proton centers apply a margin of 3.5 % + 3 mm, which leads to an overshoot of 1 cm for a proton range of 20 cm in soft tissue [12], limiting the benefits of proton therapy and increasing the dose received by surrounding healthy tissues.
In vivo range verification techniques have been proposed to reduce the uncertainty in the determination of the proton range. Some of these non-invasive methods are PET and prompt-gamma imaging, which are based on the detection of secondary photons resulting from nuclear interactions [6]. An alternative method for in vivo range verification is protoacoustics [13], also named ionoacoustics. This promising technique, which is still under development, is based on the detection of acoustic pressure waves generated as a result of energy deposition by proton beams in tissue (Fig. 1b). This is known as radio-induced thermoacoustic effect. The main advantage of protoacoustics compared to other methods is that the necessary equipment is more compact and less expensive [14,15].
The first experiment in which the acoustic signal induced by proton beams was demonstrated was performed in 1979 [16], and the acoustic pulses generated in a patient during treatment were detected for the first time in 1995 [17]. Recently, simulation studies have characterized the shape of the acoustic signals [18], and have established the relation between the proton pulse characteristics (pulse duration, beam energy and beam diameter) and the amplitude of the generated acoustic wave [19,20]. Measurements have demonstrated the feasibility of determining the proton range with submillimeter accuracy in water [21][22][23], while simulation and experimental studies in heterogeneous media have been able to reconstruct the Bragg peak position using different approaches, such as 3D filtered back projection [24][25][26], time-of-flight (TOF) calculation [27], time-reversal (TR)-based reconstruction [28], model-based inversion reconstruction [29] and iterative reconstruction algorithms [30].
It should be noted that, even though protoacoustics and photoacoustics are based on the same physical process, there are some differences between them. In photoacoustic imaging, a laser of optical photons induces the acoustic waves, so that the different absorption properties of tissues allow to form images of specific structures such as hemoglobin, lipids or melanin, using different image reconstruction algorithms to improve the quality of the image [31][32][33][34]. To the contrary, in protoacoustics the acoustic waves are generated by the proton energy deposition in tissue, being the Bragg peak area the most predominant source of those waves, represented as spherical waves in Fig. 1b. Therefore, the information that the protoacoustic wave carries depends on the shape of the Bragg peak, given by the beam energy and width, although the wave propagation is also affected by the geometry of the studied object.
One of the main difficulties in the implementation of photoacoustic/ protoacoustic is the need to simulate, as realistically as possible, the protoacoustic signal propagation through the patient to the ultrasound detectors. Sophisticated software packages, such as k-Wave [35], provide high quality image reconstruction from acoustic signals, as can be seen in the case of photoacoustic images [36]. However, the accuracy of these models often goes hand to hand with a large computational cost which makes it difficult to implement this technique on-the-fly during the radiotherapy procedures.
In this work we propose a different method for proton range verification in which we use the a priori computation of the dose distribution map, taken from the treatment plan, to create a signal dictionary. This dictionary contains simulated acoustic signals (propagated to the detection point) of all possible proton pencil beams that integrate the dose plan. These precomputations can be done once the plan is approved, at least one day in advance of the treatment, giving the physics team enough time to precalculate the signal dictionary with a realistic model including the full physics of the propagating waves. In contrast to other methods such as filtered back projection, time reversal and model-based inversion, our approach does not involve explicit reconstruction of dose maps. Instead, we compare the measured signals with the precalculated dictionary to identify, for a given signal, the most probable dose map which originated it. This system is very efficient at detecting treatment delivery errors, when the delivered spot differs from the "most-probable" spot identified by the acoustic transducers. Modern proton therapy treatment plans (in its most common modality of pencil beam scanning) are usually formed by thousands of proton pencil beams, which are delivered with a time spacing of the order of 10 ms [37], 'painting' the dose distribution in the planned target volume [7]. This spacing between consecutive beams is sufficient to distinguish the acoustic pressure wave generated by each one of them, as the acoustic wave takes approximately t = d/c = 0.2 ms to travel from its source to a detector placed at d = 30 cm, considering the water speed of sound value c = 1500 m/s. Since the dictionary is calculated before treatment and no reconstruction is required, this method requires less computational time and it is per se regularized, making possible its application for in vivo range verification. The aim of this simulation is to achieve accurate range verification (localizing the dose maximum within 1 mm) with a reduced number of detectors (under 5) placed on the skin of the patient.
We present a proof-of-concept study for the proposed dictionarybased protoacoustic dose reconstruction method. The study consists of four stages. First, we take a two-field proton treatment plan created for a prostate tumor patient, and study 4 coplanar realistic positions to place simulated detectors. The study is limited to a single axial slice of the patient and from a single field of the plan, comprising a total of 426 pencil beam spots. Second, we use the acoustic dose propagation code k-Wave to simulate the thermoacoustic waves created by each of the pencil beams in the plan, as they reach the simulated detectors. The simulation includes a full description of the temporal dose deposition pattern (as delivered by the accelerator), as well as specifics on the frequency response of the simulated transducers. These sets of simulated acoustic waves, one for each pencil beam, form the plan-specific protoacoustic dictionary. Then, we describe an algorithm to infer the position of a pencil beam from the detected acoustic wave in a number of detector positions, using the a priori information from the precalculated dictionary. Finally, we apply a number of tests to the described algorithm, in order to estimate its usefulness as a range verification tool, by measuring its sensitivity to density, geometry and range variations.

Theoretical background
Proton beams deposit energy in tissue mainly via electromagnetic interaction. This energy deposition results in a small increase in temperature which produces a thermoelastic expansion and a local pressure increase, which in turn induces a measurable ultrasonic pressure wave [13]. This process is known as the radio-induced thermoacoustic effect, and it is also used for photoacoustic and optoacoustic imaging [38][39][40][41][42]. Under thermal and stress confinement conditions [29,43], the initial pressure distribution can be calculated from the dose distribution as [29,43]: where Γ(r) is the Grüneisen coefficient, a material-specific dimensionless parameter that indicates the conversion efficiency between the absorbed heat energy and the induced pressure (Γ = c 2 βC p , where c is the speed of sound, β is the isobaric volume expansion coefficient and C p is the specific heat capacity) [43]; D(r) is the dose distribution and ρ(r) is the medium mass density. The proton pulse temporal profile is defined as a gaussian distribution G(t). So in order to take into account the finite proton pulse duration, the time varying initial pressure distribution is calculated as the convolution between eq. (1) and the gaussian distribution [28]: The initial pressure distribution p 0 (r, t) produces an acoustic wave that propagates across human tissue following wave equation [40]: where p(r, t) is the induced pressure in the position r and time t and H(r, t) is the heat function, which is related to the initial pressure in eq. (1), since H(r) = D(r)ρ(r) if thermal and stress confinement conditions are fulfilled [43].

Construction of the protoacoustic dictionary
The CT image of the patient, together with the three dimensional dose distribution provided by the treatment-planning system, can be used to precalculate, for each pencil beam, the initial pressure distribution p 0 (r, t) generated by its corresponding dose distribution D(r) applying eqs. (1) and (2). Then, numerical software k-Wave (version 1.2) for MATLAB [35] is used to solve the wave equation (eq. 3), as well as to simulate the propagation of the initial pressure signal across the patient CT to a set of predefined detector positions, executing the CUDA code on GPU using the MATLAB parallel computing toolbox.
Each pencil beam included in the treatment is identified by the set of signals that it originates at the detector positions and by the position coordinates on the CT slice (i, j) where its dose maximum is located. This 2D position of the dose maximum is used to uniquely identify a pencil beam and retrieve its information when the dictionary is utilized.
For a standard plan slice with a single field and 426 beams, the total size of the dictionary is in the order of 20 Mb. To speed up the calculation of the acoustic signals, k-Wave simulations were performed in a NVIDIA GeForce RTX 2080 Ti GPU. The total computation time to calculate the acoustic signals of the dictionary was 90 min, which is approximately 25 times faster than the time it would have taken if the simulations were carried out on the CPU. The simulation time on GPU is short enough to precompute the dictionary before treatment.
It should be noted that this relatively high computational time arises due to the size of the computational grid, given by the size of the CT image (specified in section 3.1.) and the size of the source, i.e. the dose distribution, which must match; and from the fact that the simulation is performed in a realistic way over heterogeneous human tissue. However, the resulting data, that is the acoustic dictionary (formed by the acoustic signal induced by each pencil beam and measured by all the detectors set, 3000 × 4 × 426), has itself not a large size, which makes its use easy (it does not require high storage capacity) and fast (negligible time to load the data).

Localization of the position of maximum dose
We chose a parameter to evaluate similarity between a given measured signal and precomputed dictionary signals. To quantitatively assess this difference, we use the averaged, absolute difference between them as a metric [44]: where Nt is the total number of time steps (total duration of the signal piece), Ns is the number of detectors, D s (t; i, j) is the dictionary signal at detector s identified by the coordinates (i, j) on the CT slice where the pencil beam that originated it deposits its maximum dose, and s s (t) is the measured signal at detector s. This was found to be the most robust metric for the seeked purpose after testing other widely-used metrics such as the mean square error and the correlation. Before calculating the metrics in eq. (4), we multiplied the signals D s (t; i, j) and s s (t) by the following exponential filter: where t 0 is the trigger that indicates when the filter starts to be applied, and it is defined as 10 μs before the arrival time of the first peak of each measured signal. As it will be explained in the following section, the acquisition time goes from 0 to 250 μs. So, for arrival times t close to the trigger value t 0 , the difference t − t 0 in the exponential in eq. (5) will be in the order of 10 − 6 − 10 − 5 s, while for the final arrival times it will be in the order of 10 − 4 s. The aim of the exponential filter is to give major relevance to the first peak signal, which is the most reliable part of the signal, and attenuate the end of the signals. Therefore, the attenuating factor k in eq. (5) must be in the order of 10 4 . In this case, k = 4⋅10 4 was found to be the optimal value for the function f(t). The use of this exponential filter is suitable since the last part of the acoustic signals may be affected by unwanted echoes, resonances and reflections of the waves in bone structures of immediately upper and/or lower CT slices.
In this way, we made the calculation of ε(i, j) more robust.
Besides, the attenuated signals had to be normalized so that the Bragg peak identification did not depend on the relative amplitude between measured and dictionary signals. All attenuated signals were normalized by its maximum value. Other methods, such as the normalization of the attenuated signal squared by the sum of its values and the normalization of the attenuated signal by the sum of its absolute values, were tested. These two last methods provided acceptable results, but the best and most robust results were obtained with the normalization by the maximum value, and therefore these are the ones presented in this paper.
For each set of signals {s 1 (t), s 2 (t), ..., s n (t)} measured by the detectors 1…n, coming from a given pencil beam, the error value ε(i, j) is calculated for the complete set of signals saved in the dictionary, where the position (i, j) identifies the pencil beam with dose maximum located on coordinates (i, j) on the CT slice. The graphical representation of ε(i, j) over the CT slice (see Fig. 2) clearly indicates the "most-probable" pencil beam originating the set of signals {s 1 (t), s 2 (t), ..., s n (t)}, corresponding to the beam with maximum on (i, j) which minimizes the error ε(i, j). The error value distribution within the mask is smooth and convex, making the determination of the minimum simple and robust. This in-house developed code in MATLAB takes about 100 ms to determine the minimum error, and hence the most-probable pencil beam giving rise to the measured set of signals.

Benchmark study: methods
To evaluate the proton range verification method that we propose, we performed a simulation study over a prostate case, where we precalculated the dictionary of signals and used it to reconstruct the dose position under different conditions. We simulated six different changes from the original plan representing possible sources of uncertainty in proton range determination: proton range variations caused by anatomical changes, patient shifts in transverse position, and changes in medium density. That is, we focused on the most important source of uncertainties, corresponding to changes in the patient. The proton beam (energy, pulse width…) is considered to be well characterized with daily QA activities, and therefore with negligible uncertainty compared to the other factors. We performed the evaluation of the six cases for three different pencil beams of the plan (Fig. 3), depositing their maxima (that we fixed as 1 cGy) in different regions of the target.

Simulated treatment configuration
A 90-slice prostate CT image was taken from the CORT dataset [45]. A 2-field, intensity-modulated proton therapy plan was created with matRad [46] using default quadratic objectives on tumor volume, normal tissue and select organs at risk (bladder, rectum and femoral heads). Proton beams with energies ranging from 120 to 196 MeV were used to build the plan. In the default proton machine provided by matRad, these beams have a gaussian distribution with FWHM (full width at half maximum) ranging from 8.9 to 10.0 mm. The dose deposited by each pencil beam of the plan is approximately 3 mGy. CT slices dimensions are 219 × 219 (voxel sizes of 2.5 mm in X and Y).
Without loss of generality, we evaluated the dictionary system over the CT axial central slice, simulating only the pencil beams whose dose maxima laid on this slice according to the calculated plan. In total, the number of pencil beams in this slice was 426. To take into account the beam total width in Z direction (vertical axis), which increases along its path in human tissue, we performed the acoustic simulation on a 3D grid of 219 × 219 × 12 voxels, that is, including the five immediately lower and the six immediately upper slices of the CT, with a voxel size of 2.5 × 2.5 × 3 mm 3 . An absorbing boundary layer, called perfectly matched layer (PML) in the k-Wave toolbox, of 7 × 7 × 2 voxels is considered inside the grid in order to prevent waves leaving one side of the computational grid to reappear at the opposite side [27,28,35]. To perform the acoustic simulation in k-Wave, we considered for simplicity a uniform Grüneisen coefficient of Γ = 0.8 [47], as in [29]; and medium density and speed of sound values were calculated from the CT Hounsfield Units (HU) based on the works of Schneider et al. and Mast [48,49].
The proton beam temporal structure is given by a Gaussian profile, as it is indicated in eq. (2), with σ = 7 μs, matching the σ value of the IBA S2C2 synchrocyclotron [50]. This is the machine with the most advanced proton technology developed by IBA (one of the leading companies at proton therapy technology), which is already implemented in more than 20 proton therapy facilities in Europe, America and Asia [51,52]. It should be noticed that the time of dose deposition, given by the σ value, is in the same order that the time a perturbation takes to propagate along a single voxel, which is Δt = dx/c = 1.67 μs, with dx = 2.5 mm the voxel size and c = 1500 m/s the approximate value of the speed of sound in soft tissue. Therefore, this choice of the voxel size ensures that the time of proton beam deposition is well supported within  the k-Wave simulation, and so it also is the maximum frequency of the generated signals. The sampling time step was 83 ns, and the total sampled duration was 250 μs.
To record the simulated acoustic signals, four ideal point detectors were located over the patient's skin of the numerical phantom, in positions in which bone structures interfere as little as possible. This could be a possible configuration for a real-case scenario. We compared results with 2-6 detectors, and we found that the accuracy did not improve beyond 4 detectors, so that is the suitable number to get the best accuracy and the most cost-effective and comfortable setup at the same time. Fig. 3 shows the detectors setup over the CT image. Random, uniform noise between -0.1 and 0.1 mPa was added to the acoustic signals to make the simulations as realistic as possible, and a low-pass filter with a frequency band of 100 kHz was introduced to mitigate the effect of highfrequency noise.
Acoustic attenuation can be considered negligible in the frequency range of the protoacoustic signals, which is approximately between 0 and 100 kHz (see Fig. 6c). The absorption of acoustic waves is defined by the absorption coefficient [44]: where α is expressed in units of Np/m, α 0 is the power law pre-factor in Therefore, the value of the absorption coefficient α is very low for frequencies considerably smaller than 1 MHz, as it is in this situation. The attenuation of the signal amplitude is expressed as [44]: where A 0 is the non-attenuated signal amplitude and A is the reduced amplitude when the acoustic wave has traveled a distance d from the initial location. Thus, for frequencies in the order of tens of kHz and distances around 10 cm, the exponential in eq. (7) has values close to 1, and the signal amplitude barely attenuates. Moreover, the low-pass filter eliminates the high frequency part of the signal, which could be more affected by attenuation. Hence, it is a good approximation to consider the attenuation is negligible in this work.
On the other hand, scattering was considered by adding Gaussian variations of the speed of sound (20 m/s FWHM) and density (5% FWHM) in tissues.

Selection of test cases
Six different cases of variations with respect to the original treatment situation were analyzed. In all these cases, the dose distributions of the Beams 1, 2 and 3 and the propagation of the acoustic waves related to them were calculated based on the corresponding, modified CT slices.
Range variations caused by anatomical changes. Three different situations were studied. First (test case #1), an air region (HU = -1000) with a size of 1 × 1 cm 2 was inserted in patient, near the beam entrance channel, to represent anatomical changes (i.e. air bubbles) in thoracic or pelvic cases. Test cases #2 and #3 represent milder variation of patient anatomy, designed to test the sensibility of the proposed method. For test case #2, a region of soft tissue of 2 × 1 cm 2 , with HU = 40-80, was inserted in the proton entrance channel, simulating patient weight gain. Finally, in test case #3, material of a single voxel was changed (from soft tissue to water) in the patient CT. These three material changes cause moderate to small variations in proton range (from 12.0 mm in test case #1 to 0.1 mm in test case #3). Fig. 4 shows the effect of the addition of an air region on the CT and on the proton beam range.
Shift in patient position. A global shift of 1 mm in transverse direction, representing patient misalignment during treatment setup, was simulated. This was labeled as test case #4.
Medium density variations. Variations in medium density were studied by introducing a decrease of 3.5 % in CT Hounsfield Units in the whole patient tissue, which affect the amount and the location of the  (Fig. 4a) and the dose deposited when there is an air region in the beam way (Fig. 4b).
deposited dose. Two different situations were analyzed: on the one hand, the propagation of the generated acoustic signal was performed over the original CT (test case #5), thus studying errors in the CT Hounsfield Units conversion to relative stopping power. On the other hand, the acoustic propagation was calculated over the modified CT (test case #6), studying in this case errors in CT-HU values.

Evaluation of the method for error detection
To evaluate the accuracy of the proposed method for the detection of errors, we calculate the following parameters: where Z 0 is the original position of the Bragg peak, Z 1 is the real position of the Bragg peak when the variation is introduced, Z 0,ac is the original position of the Bragg peak calculated with the dictionary-based method and Z 1,ac is the Bragg peak position estimated by the dictionary-based method when a variation is introduced. That is to say, ΔZ real represents the real change in the Bragg peak position in the different test cases explained in section 3.2, while ΔZ ac represents the change that the dictionary-based method is able to detect. These four Bragg peak positions (Z 0 , Z 1 , Z 0,ac and Z 1,ac ) are interpolated using splines 2D method, with an interpolation step of 0.1 mm, to achieve better accuracy (see Fig. 5). Then, the level of accuracy is given by the absolute difference between ΔZ real and ΔZ ac , and by the relative error: ε r = Accuracy⋅100/ΔZ real (11) In each test case, the average accuracy is calculated as the mean of the values of accuracy (eq. 10) and relative error (eq. 11), respectively, obtained for the three beams. The total average is calculated as the mean of the values of accuracy and relative error obtained in all test cases.

Simulated acoustic pressure
The resulting set of simulated acoustic signals for test beam 1 (Fig. 6  a) is shown in Fig. 6b. The maximum pressure amplitude received by the detectors ranges between 5.0 and 1.3 mPa, for detectors placed between 9.25 and 22 cm away from the Bragg peak and a pencil beam of 5 × 10 7 protons depositing a dose of 1 cGy in the voxel of maximum dose, which is the maximum dose value we set in the simulations of Beams 1, 2 and 3 for the test cases (as previously said in the beginning of Section 3). That is, the nearest simulated detector measures 0.5 Pa/Gy, while the farthest one measures 0.13 Pa/Gy. The frequency spectrum of the signals ranges from 8 to 100 kHz, with a central frequency of 28 kHz (Fig. 6c). Table 1 summarizes the results obtained for the real change ΔZ real (eq. 8) and the detected change ΔZ ac (eq. 9) when the different anatomical changes are introduced. The presence of an air region (test case #1) near the beam entrance channel produces the highest variation in the Bragg peak position (between 11.3 and 12 mm). The acoustic dictionary-based method indicates changes in the maximum dose position between 9.0 and 12.4 mm. That is to say, the method is able to detect that some error in dose delivery is happening. The corresponding values of accuracy (eq. 10) and relative error ε r (eq. 11) are also shown in Table 1. The level of accuracy is found to be between 0.4 and 2.6 mm, giving an average accuracy to estimate the error for the three beams of 1.8 mm. The relative error ranges between 3 % and 22 %, which results in an average value of 15.2 %. The introduction of a soft tissue region causes changes in the Bragg peak location between 1.4 and 1.7 mm for the three test beams, while the acoustic method identifies that the variations are in the range of 2.0-2.6 mm. So again, the disagreement between the planned Bragg peak position and the real position is detected. Moreover, our method estimates how much the Bragg peak is displaced in relation to the plan

Table 1
Bragg peak original position (Z 0 ), Bragg peak original positions calculated with the dictionary-based method (Z 0,ac ), real Bragg peak position in a test case (Z 1 ), Bragg peak position estimated by the dictionary-based method in a test case (Z 1,ac ), real change in Bragg peak position in the test case (ΔZ real ) and change detected by the dictionary-based method in the test case (ΔZ ac ). Accuracy and relative error are calculated using Eqs. (10) and (11) with an accuracy between 0.5 and 0.9 mm (average value of 0.7 mm). The relative error is in the range of 31-53 %. Lastly, the change in the material of a single voxel from soft tissue to water only causes a variation in the dose location in one of the three test beams (Beam 2), being the Bragg peak shifted just 0.1 mm. For Beam 1 and Beam 2, the detected change ΔZ ac perfectly matches the real change ΔZ real , so the relative error is 0 %. For Beam 3, there is a detected change of 0.3 mm when actually no variation has been created in the Bragg peak position. However, this difference of 0.3 mm between ΔZ real and ΔZ ac can be considered within the uncertainty of the method. In this test case, the average accuracy is 0.1 mm. Table 1 collects as well the real and detected changes in maximum dose position when a 1 mm displacement in patient transverse position is simulated, and the corresponding values of accuracy and relative error. This situation originates shifts in the Bragg peak positions of 0.3− 0.6 mm, while the dictionary-based method detects variations of 0.4-1.5 mm. For Beams 2 and 3 the accuracy is 0.1 mm, and the relative error is, respectively, 20 % and 33 %. Therefore, for these two beams the method is able to detect the submillimetric change with great accuracy. On the other hand, for Beam 1 the difference between real and detected change is 0.9 mm, that is, greater than the variation that is caused, which gives a relative error higher than 100 %. So in this case, the method detects that there has been a change in the Bragg peak location, but it overestimates the amount of the variation.

Detection of medium density variations (test cases #5-#6)
The 3.5 % decrease in the CT Hounsfield Units produces changes in the Bragg peak position between 5.3 and 8.1 mm for the different test beams (see Table 1). However, the detected change ΔZ ac presents different values depending on whether the acoustic propagation was simulated over the original CT (test case #5) or the modified CT (test case #6). In test case #5, the detected change values are between 7.3 and 8.0 mm, which indicates that the dictionary-based method senses these variations with an accuracy between 0.4 and 2.7 mm, being the average value 1.6 mm. Therefore, in test case #5 the acoustic method is able to detect the change originated, which represents an error in CT-HU conversion to relative stopping power. Besides, the relative error is 28 %, 5 % and 51 % for Beam 1, 2 and 3, respectively, so for the first two beams the method estimates with acceptable accuracy the value of the error.
On the other hand, the ΔZ ac values in test case #6 are 0.6, 5.1 and 1 mm for Beams 1, 2 and 3, which give an accuracy in the range 3.0-5.1 mm. This represents an average accuracy of 4.1 mm, that is to say, more than two times smaller than in test case #5, and also worse than the values presented in the other four test cases. The relative error is between 37 % and 89 %, which clearly points out that the acoustic method is not optimal for accurately detecting errors in CT-HU values in its current form.

Discussion
According to the results presented in Table 1, the dictionary-based protoacoustic method is able to detect changes in the dose maximum position due to material anatomical variations, which are the most common source of uncertainty in interfractional treatments. These changes are represented here by the presence of an air region and a soft tissue region in the way of the test beams (test cases #1 and #2). In the case of the region of air, the mean value of the relative error is 15 %, which indicates that the proposed method is also able to accurately quantify the deviation of the Bragg peak position when the deviation is in the order of 1 cm. In the case of the soft tissue region, the higher relative error values obtained are due to the fact that the changes originated in the Bragg peak position are relatively small (around 1.5 mm), so the precision needed to achieve a relative error similar to the one in test case #1 would be around 0.25, which is four times better than the accuracy of 1 mm that we aim to get and that has been achieved in other range verification studies. Still, it is remarkable that the proposed dictionary method detects variations of less than 1 mm. Moreover, test case #3 shows that the method is still robust when a material change in some voxel is introduced, since the average accuracy obtained in this case is 0.1 mm.
Regarding displacements in patient transverse position (test case #4), the proposed method can also detect the variation in the Bragg peak location originated by patient positioning errors, and this change can be quantified with an average accuracy of 0.4 mm. In addition, the method is able to identify errors in the conversion from CT Hounsfield Units to relative stopping power, since it accurately detects the changes caused by a 3.5 % decrease in CT Hounsfield Units when the acoustic propagation is performed over the original CT (test case #5). Nevertheless, the method is not so accurate in detecting errors in CT Hounsfield Units themselves (test case #6). This is due to the fact that an error in CT-HU not only affects the dose deposition, but also the acoustic propagation, on which is based the dictionary method. In fact, in light of the results for test cases #5 and #6, the effect of a global change in CT-HU is more significant in the acoustic calculation than in the dose calculation. This would limit the ability of the proposed method to detect this kind of errors, although it should also be considered that in most cases CT uncertainties are related to Hounsfield Units conversion (test case #5), not to CT-HU values.
It should be taken into account that this dictionary-based method can identify changes in the dose maximum location, but in a real scenario it will not be able to distinguish the specific nature of the disturbance that originated the change. This could be considered as a limitation of the method, but to the best of our knowledge, no other state-of-the-art range verification method is able to identify the source of deviations in the Bragg peak positions.
It should be noticed as well that even though in some cases the estimated Bragg peak position (Z 0,ac ) does not perfectly match with the real Bragg peak location (Z 0 ), as it happens for Beam 1, this does not affect the ability of the method for errors detection, since that difference between Z 0,ac and Z 0 is systematic and stems from the modification of the acoustic pulse shape when noise is added. Therefore, the dictionarybased method for errors detection can identify changes in dose delivery independently of its aptitude for locating the Bragg peak exact position, which makes it a stable and robust method.
The pressure values presented in this work (between 1.3 and 5.0 mPa) are in accordance with the results presented in several, previous experimental and simulation studies. In [54], a pressure amplitude of 5 mPa is measured in water with a detector placed at a distance of 5 cm from the Bragg peak after averaging 2048 pulses, each one containing 1.2 × 10 7 protons and depositing a dose of 1.9 cGy, with a temporal pulse of 18 μs FWHM and beam width of 8 mm at the entrance.
A Brüel and Kjaer 8105 hydrophone was used, which has an omnidirectional, uniform response in the range 10 Hz -100 kHz. In [24], simulations in a water phantom of proton beam dose deposition of 1 cGy at the Bragg peak in one pulse of 1 μs and 8 × 10 6 protons (beam width of 1 cm FWHM) resulted in pressure amplitudes between 8.7 and 25.5 mPa for a transducer array placed along the beam axis (closest transducer at 15 cm, farthest transducer at 25 cm), with a detector bandwidth of 0.5 MHz. Another simulation study in a water phantom [55] reported a pressure peak of 140 mPa induced by a proton pulse of 10 μs with a 64.7 cGy dose (protons per pulse in the order of 5 × 10 6 ), 5 mm spot size, using a realistic PZT ultrasound transducer (5 cm diameter) placed 4 cm distal to the Bragg peak in the simulation. The effect of the bandwidth in the signals was simulated with a Butterworth bandpass filter. The maximum frequency of the signals was 100 kHz, centered around 76 kHz. Finally, a similar simulation study [27] to the one presented in this work indicated pressures of 160− 50 mPa in a prostate case for a proton beam of 1.3 cGy at the Bragg peak and an ideal point detector at 1.66 cm, and a 1 cGy beam and a detector at 7.1 cm, respectively (pulse of 1 × 10 7 protons and Gaussian temporal profile of 14 μs FWHM). All these results are consistent with the ones obtained in our work (see Table 2), considering that the pressure amplitude decreases with higher Bragg peak-to-detector distance and with larger temporal pulses. Also, a greater number of protons (or greater proton beam dose) increases the pressure amplitude of the generated acoustic signal, while simulations in homogeneous media (i.e. water) gives higher pressure peaks in comparison to simulation in heterogeneous tissue.
In this study we used ideal point detectors, as in [27]. The main differences between these ideal detectors and realistic detectors are the finite size, bandwidth, and directionality of the laters. The size of the detectors could be represented in the k-Wave simulation by averaging the signals recorded by several point detectors at single voxels. However, the contribution of each point detector to the final signal would have to be determined taking into account the acoustic center and directionality of the real detector to accurately represent that particular real detector. This will have to be considered when carrying out experimental measurements with specific detectors, but it was considered not worthy for the moment since this work aims to present general results for the novel dictionary-based method. Regarding directionality, there are omnidirectional detectors, like the ones employed in [54], which are well represented by the point detectors here simulated. Constraints in directionality could be introduced in the simulation of finite size, but again it would depend on the specific detector that would like to be represented. Finally, the limited bandwidth of realistic detectors is represented in our simulations by the low-pass frequency filter of 100 kHz. Bandwidth could cause changes in the arrival times of the signal, but these will be negligible if detectors with proper sensitivity are used [27].
Just as it was explained in section 3, the dose simulation in matRad and the acoustic simulation in k-Wave were performed with a voxel size of 2.5 mm both in X (transverse direction) and Y (depth). The dose simulation could be carried out with a smaller voxel size; however, that will not change our ability to obtain an accurate Bragg peak positioning based on the acoustic signals simulated with k-Wave software and stored in the dictionary. Hence, we did not consider it necessary.
Compared to previous simulated protoacoustic studies in heterogeneous media [27][28][29], the method we propose needs a smaller number of detectors to determine the proton dose distribution, as it uses the a priori known shape of the proton beamlets as a strong regularizer of the inverse problem. As a consequence, the setup needed to put into practice this method would be very cost effective and more comfortable for the patient. It is expected that the use of more realistic detectors will not imply an increase in the number of devices required to carry out our method, although it should be proved with experimental measurements. Besides, the dictionary is calculated before the treatment, since the dose plan is usually available the day before. Therefore, our method can provide the precalculated dictionary before the treatment, so error values (eq. 4) can be obtained in real time during treatment, allowing for in vivo range verification. In a way, we trade off precomputation time by the ability to perform the analysis on-the-fly.
It is also important to mention that although in this work we have only considered cases in which the detectors are all located in the same plane as the proton beam, different detector distributions could be easily considered (for instance, placing them on different CT slices). The 3D acoustic simulations performed in this work would not change if other positions of the detector elements are considered. How to define the optimal locations of the detectors is out of the scope of this work, and it is a problem that deserves further studies.
To carry out the proposed method in a real-case scenario, it would be necessary to determine with great accuracy the exact position of the detector sensitive area, so that we can ensure that the positions where the detectors have been placed in the simulation of the dictionary match the positions in the real patient setup. This could be measured beforehand by ultrasound means, as it is proposed in [22].
Although acoustic signals generated by proton beams in patients have been already detected in previous studies [17], one of the main drawbacks of protoacoustics is still the low signal-to-noise ratio (SNR). Our results show that the frequency range of interest of the protoacoustic signals is up to 100 kHz, as it has been reported in other studies as well [22,32,54,55]. Thus, detectors with uniform frequency response and maximum sensitivity in this range will provide the best SNR. Experimental measurements with suitable detectors will determine if some signal processing (for example, averaging) is needed to achieve the required SNR.
The method proposed in this work has the highest potential in combination with active scanning beam delivery systems. The dictionary strategy would be different for passive scattering delivery system, since this technique does not employ pencil beams, and the basis of our method is precisely to correlate each proton pencil beam with its corresponding acoustic signal. However, active scanning is the technology of choice in essentially all modern proton therapy facilities [56], under planning or in construction, because it allows to achieve a better dose conformity [6]. Besides, active scanning provides significantly smaller beam spots and pulse length (see Section 1), which favors a higher acoustic amplitude.

Conclusion
This work demonstrated the feasibility of proton range verification in proton therapy employing a precalculated dictionary, which is formed by the simulated, expected acoustic signals that proton beams generate during patient treatment, and a cost-effective equipment. The use of a precomputed dictionary opens the path towards real time, on-the-fly, verification of range, independently of the computational burden that the calculations based on realistic acoustic simulations can require. The accuracy of this novel approach was quantified by the comparison between the real change in the Bragg peak position, ΔZ real , and the Table 2 Summary of the results obtained in similar simulation and experimental studies of the acoustic pulse induced by proton beams and the results obtained in this work. The main characteristics that affect the acoustic pressure amplitude are gathered: Bragg peak-to-detector distance, protons per pulse, dose deposited by the proton beam, proton pulse duration, proton beam width, detector bandwidth and propagation medium. detected change, ΔZ ac , being on average 1.4 mm for the eighteen presented cases. This method is applicable with pencil beam scanning systems, available in all modern proton facilities. Experiments are planned to validate the results of the simulation study.

Funding
Work supported by several public agencies: Spanish Government (FPA2015-65035-P, RTC-2015-3772-1, RTI2018-095800-A-I00), Comunidad de Madrid (B2017/BMD-3888 PRONTO-CM), European Regional Funds and the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 793576 (CAPPERAM). This is a contribution for the Moncloa Campus of International Excellence, "Grupo de Fisica Nuclear-UCM", Ref. 910059. Part of the calculations of this work were performed in the "Cluster de Calculo para Tecnicas Fisicas", funded in part by UCM and in part by EU Regional Funds.

Declaration of Competing Interest
The authors declare that there are no conflicts of interest.