Acoustic Bragg Peak Localization in Proton Therapy Treatment: Simulation Studies

: A full chain simulation of the acoustic hadron therapy monitoring for brain tumors is presented in this work. For the study, a proton beam of 100 MeV was considered. In the first stage, Geant4 was used to simulate the energy deposition and to study the behavior of the Bragg peak. The energy deposition in the medium produced local heating that can be considered instantaneous with respect to the hydrodynamic time scale producing a sound pressure wave. The resulting thermoacoustic signal was subsequently obtained by solving the thermoacoustic equation. The acoustic propagation was simulated by the Finite Element Method (FEM) in the brain and the skull, where a set of piezoelectric sensors were placed. Lastly, the final received signals in the sensors were processed in order to reconstruct the position of the thermal source and, thus, to determine the feasibility and accuracy of acoustic beam monitoring in hadron therapy.


Introduction
The use of heavy particles as protons in the treatment against solid tumors in areas not accessible with standard photon radiation has been increasing in recent decades [1]. This is shown in the number of new accelerators in new centers of proton therapy [2]. The main advantage of protons over photons is the dose distribution, minimizing the damage to tissue out of the region near the Bragg peak region. The location of the Bragg peak profile depends mainly on the initial energy of the beam and its width. This deposition is particularly important in clinical conditions with a total energy deposit, such as melanoma of the eye or tumors in brain tissue, among others, where the possibility of damage to nearby tissue increases due to the size of the tumor tissue. Although monitoring and control systems are carried out before and after each radiation session, different alternatives for the monitoring process have been proposed based on physical processes produced in the irradiated medium such as prompt-gammas, charged particles, βþ-emitter, and PET-imaging [3]. This article proposes an alternative method based on the pressure produced by the proton beam, the amplitude of which was evaluated from a set of piezoelectric sensors. In particular, it is based on simulated studies of the absorbed dose distribution for a set of energies from the Monte Carlo methods. Thus, the thermoacoustic model calculates the pressure received at a point in space that will be one of the input parameters for the Finite Element Method (FEM) that solves the behavior of the piezoelectric effect in a piezoelectric (PZT) material as the set of sensors for brain tumor applications.

Monte Carlo Simulation
In order to calculate the energy deposited in the water due to the Bragg peak, the Monte Carlo methods were applied using the Geant4 hadron therapy simulation package. The model features a 100 MeV monoenergetic proton beam using the Electromagnetic Quark-Gluon String model with Binary Cascade for primary protons (QGSP BIC EMY) as Physics List [4]. The absorbent material was G4Water with a density of 1003 kg m ⁄ and a number of simulations of about 10 histories. The detector volume consisted of a solid G4Box with a voxel size of 1 mm. To evaluate the boundary conditions of bone tissue, the Anthropomorphic model was used within the Geant4 libraries. The models were simulated in the volume of water including the bone structure. Generally, in medical physics applications, the results of Geant4 are evaluated in dose deposit per unit area, however, in this study, the energy of the beam deposited in the tissue was evaluated. Thus, the behavior of the Bragg peak was used in its energy part in the thermoacoustic model described in the following section.

Thermoacoustic Model
According to the Thermoacoustic model [5], the deposition of energy from particles that pass through a liquid produces local heating in the medium that can be considered instantaneous with respect to the hydrodynamic time scale. Due to the temperature change, the medium either expands or contracts according to its coefficient of volumetric expansion . The accelerated movement of the heated medium forms a pressure pulse that propagates through the fluid. The model is described in Equation (1), with the particular solution given in Equation (2).
where ( ⃗, ) is the pressure at a snapshot in time , and spatial position ⃗, represents the speed of sound in the middle; is the specific heat; ϵ( ⃗, ) is the energy density deposited in the medium; and is the coefficient of thermal expansion. When a pulse of protons radiates in a homogeneous volume, it corresponds to a pressure source proportional to the energy deposition ϵ( ⃗, ). Each differential volume element of the irradiated source pressure emits a micro bipolar pressure wave. The pressure measured in the sensor is the sum of the micro-pressure waves emitted by each differential element of the volume. Due to the constant speed of sound throughout the medium, from the perspective of the sensor, the pressure that arrives after a certain time t is related to the pressure waves emitted in a common radius r. As shown in Equation (1), the derivative of these pressure waves translates a bipolar behavior that is eventually the signal received by the sensor.

FEM Simulation
To solve the problem of propagation of the pressure calculated at a certain distance from the Bragg peak, the FEM model was used for acoustic-structure iteration, which combined the Pressure Acoustics, Transient, and Solid Mechanics interfaces to connect the acoustics pressure variations in the fluid domain with structural deformation in the solid domain. In essence, the default sound pressure function models harmonic sound waves in the water domain using Equation (3) for sound pressure.
where is the pressure, and represents the density. is the monopole domain source that corresponds to a mass source on the right-hand side of the continuity equation. The dipole domain source corresponds to a domain force source on the right-hand side of the momentum equation. The combination is called the adiabatic bulk modulus, commonly denoted . The bulk modulus is equal to one over the adiabatic compressibility coefficient = 1⁄ = 1 ⁄ [6]. The pressure generated by the thermoacoustic model is distributed in a sphere within the skull corresponding to the range of the Bragg peak deposition. Figure 1 shows the geometric model implemented in COMSOL Multiphysics and the positions of the sensors. The pressure value propagated to a spherical volume of 2 mm in diameter, which corresponded to the Bragg peak for 100 MeV, as shown in Figure 2b.

Localization Method
For the reconstruction of the position of the Bragg Peak, a system of nonlinear equations of the form ( ) = 0, was solved using Newton's method [7]. Consider a system of equations and unknowns: ( , , , … , ) = 0.
This system can be written in vector form as ( ) = 0, where is a vector of dimensions, and is a vector of dimensions. To solve this system of equations, we have to find a vector such that the function ( ) equals the null vector. If we call to the solution of the system and to an approximation of it, we can develop in Taylor series around this approximation as: Since ( ) = 0, then we get, as an approximation, the following: defining the vector as this approximation, which is closer to the root than . We can continue with the iterative method to obtain approximations closer and closer to the solution. To write the iterative method, the term ∇ ( ) is replaced by the Jacobian of the function , that is: which is a × matrix. Then, it is possible to obtain a new value of by solving the following relationship: Then, iteratively, we can approximate the to more and more until a solution error | − | previously fixed is reached. This method has been used in different experiments of acoustic source localization [7] and compared with algorithms for solving systems of nonlinear equations.

Energy Deposition
As a result of the interaction of protons with matter following the behavior of the Bragg peak, it is possible to determine the energy deposited in the medium for a 100 MeV proton beam. For this, the deposition in the middle was evaluated directly and with a 1 cm bone layer that simulates the skull following the geometric models of the Geant4 libraries. Thus, Figure 2 shows the Bragg profile as a function of the range in water, with the bone layer described and the distribution of this energy in a plane. Figure 2 shows the range in water for 100 MeV with and without the skull material.

Thermoacoustic Signal
Because most of the energy was deposited in a volume that had a depth of approximately 2 mm, it was necessary to propagate the pressure from the simulation of different points in the space on the beam emission axis. For this, the pressure received simulated by a 5 mm to 40 mm sensor in 5 mm steps was simulated. For this, the expression = • was used where is a known reference pressure, the distance between the source and the sensor, and the water absorption coefficient. With these, Figure 3 shows the pressure signal received by a sensor 20 mm from the Bragg peak, the maximum simulated pressure values for different distances and a setting that establishes a pressure of 1.66 Pa at a distance of 2 mm.
Once the pressure behavior was simulated and established over time, the pressure value at a distance of 2 mm was set in the FEM model for wave propagation at a volume that corresponds to the peak maximum of Bragg shown in Figure 2b.

Acoustic Propagation
Once the acoustic pressure was calculated using the thermoacoustic model and obtaining the pressure propagated to a source whose volume coincides with the energy deposition, the propagation in the brain was performed using the FEM model. For this, a uniform density was taken in the region inside the skull with a value of 1003 kg m ⁄ that corresponded to the density of the water since the brain tissue and the cerebrospinal fluid have a similar density [9]. Figure 4 shows an image of the propagation in a horizontal plane in the instant of time = 53 μs. In addition, the pressure signal received on the bone inside the skull and the acceleration on the outside surface of the skull are shown. The received signal is in terms of the acceleration on the surface of the skull. This signal can be converted to pressure to have a voltage reference due to the Receiving Voltage Response (RVR) of the sensor attached to an adaptation layer. This signal is evaluated in time and correlated with the signal emitted by the source to obtain the time differences of arrival among the sensors (TDOA) which are referenced to one of the sensors. The signal received by the sensors corresponds to an acceleration analysis that delivers the FEM method obtained in one of the outputs in the solid mechanic module. Therefore, each of these signals was evaluated to get the arrival times from the source.

Acoustic Localization
In the first approach, the signals arrival times were obtained by cross-correlation between the emitted signal ( Figure 3a) and the set of received signals. Thus, it was possible to determine the TDOA and solve the system of equations for the 6 sensor locations tested. Table 1 shows the positioning of the sensors and the results of the reconstructed position by the location algorithm for energy deposition. It can be seen that there was a deviation between the reconstructed position and the simulated source of about 1.8 mm. However, in a closer approach, this reconstructed position was recalculated by correcting the different sound speed velocity by the skull with respect to the water. The results are also shown in Table 1, noting that the last deviation reduced to 1.0 mm.
The spatial reference for the sensors was given by the software (COMSOL), whereby the position of the source within the skull was known.

Conclusions
Studies were conducted on the deposition of energy in brain tissue through simulations with the Monte Carlo methods for the energy of 100 MeV. It was possible to determine the pressure at a point in space by discretizing the thermoacoustic model and propagating the pressure to a size similar to the Bragg peak. Therefore, through sensors located on the surface of the skull, the position of the source was reconstructed from the moment of detection by the different sensors, while also considering the corrections in time due to the change of impedance among the propagation media. The position of the source was solved with few iterations, using a robust resolution method for nonlinear equations. Therefore, this technique can be considered as a good alternative for monitoring the location of the Bragg peak in the treatment of hadron therapy.