Multiphysical numerical study of photothermal therapy of glioblastoma with photoacoustic temperature monitoring in a mouse head

: This paper presents a multiphysical numerical study of a photothermal therapy performed on a numerical phantom of a mouse head containing a glioblastoma. The study has been designed to be as realistic as possible. Heat diffusion simulations were performed on the phantom to understand the temperature evolution in the mouse head and therefore in the glioblastoma. The thermal dose has been calculated and lesions caused by heat are shown. The thermal damage on the tumor has also been quantified. To improve the effectiveness of the therapy, the photoabsorber’s concentration was increased locally, at the tumor site, to mimic the effect of using absorbing contrast agents such as nanoparticles. Photoacoustic simulations were performed in order to monitor temperature in the phantom: as the Grüneisen parameter changes with the temperature, the photoacoustic signal undergoes changes that can be linked to temperature evolution. These photoacoustic simulations were performed at different instants during the therapy and the evolution of the photoacoustic signal as a function of the spatio-temporal distribution of the temperature in the phantom was observed and quantified. We have developed in this paper a numerical tool that can be used to help defining key parameters of a photothermal therapy.


Introduction
PhotoThermal Therapy (PTT) is a promising therapeutic treatment for cancer relying on optical radiation with Near InfraRed (NIR) light. This technique relies on the conversion of light absorbed in the irradiated tissue into heat leading to increased temperatures, which in turn induces cell death and lead to thermal damage. Through the use of NIR optical absorbers, the efficacy of the treatment can be improved due to the increase in energy conversion. Therefore, regions of interest such as tumors can be targeted and locally treated. An important number of studies tackles the question of which absorber to use, proposing different materials but also different forms and sizes. Metallic nanoparticles (NP) (e.g. Silver, Gadolinium) were the first kinds of NP used, showing strong tunable absorption properties but also important cytotoxicity [1][2][3][4][5][6][7][8].
Recently, the attention of the research community has turned to the development of organic chromophores, presenting some advantages over their metallic counterparts: a lower toxicity and a simpler vectorization into biological systems [8,9].
In 2020, 308 102 patients were diagnosed with central nervous system (CNS) tumors, which include brain tumors. CNS tumors are one of the deadliest kind of cancers with a mortality rate of approximately 81.6% [10]. In the United States in 2010, the treatment of such brain cancer had a national cost of more than 4 billions US dollars [11]. Regarding brain tumors, glioblastoma multiforme (GBM) is the most common and aggressive primary brain tumor in adults. Despite standard treatments (including surgery, radiotherapy and chemotherapy) the median overall survival of GBM patients is 14.6 months with a five-years survival rate of 9.8% [12]. Therefore developing innovative therapeutic strategies for the effective treatment of such tumors is essential in cancer research. The delivery of drugs into the brain is limited by the blood-brain barrier (BBB). Several strategies have been described to cross the BBB and obtain therapeutic concentrations of NP and chemotherapeutics within the tumor. Some methods exploit biological processes allowing for cell-mediated transport of the drugs to the tumor site [13]. Other methods induce transient permeabilization of the BBB via the use of microbubbles and high intensity focused ultrasound (HIFU): microbubbles' vibrations created by the use of HIFU induce a transient physical separation between the cells of the BBB, allowing drugs or NP to cross the barrier [14]. Another strategy to increase the treatment efficacy is to deliver drug or drug-loaded NP directly by intratumoral injection in the brain. Local, intracerebral drug delivery allows to bypass the BBB, obtaining a higher drug concentration at the tumor site and reduced systemic side effects. In this paper, we did not consider the NP delivery process and supposed they were already successfully delivered to the tumor site.
During the PTT, a safe real-time measurement of the temperature is needed to ensure the efficacy of the treatment of the tumor tissue while controlling damages to the surrounding areas. Many methods exist and have been used in previous studies. Ultrasound imaging showed good penetration depth and spatial resolution, but also presents low sensitivity to temperature variations for temperature in the ablation range [15]. Infrared thermography allows a better accuracy but is limited to measurements of the superficial layers [16]. However, methods exist to improve the penetration depth [17,18]. Magnetic resonance imaging has a good sensitivity and penetration depth. But it has also low temporal resolution when high spatial resolution is needed [16]. The PhotoAcoustic (PA) technique has been shown to be a good solution for in depth, non invasive real time technique to monitor temperature [19][20][21][22]. A relation between temperature and the amplitude of the PA signal exists and allows this technique to be used to control the temperature evolution in the tumor and its surrounding.
With the constant improvement of PTT, new numerical tools appeared: Cuplov et al. [23] developed an extension of the open source software GATE (Geant4 Application for Emission Tomography) to perform NP-mediated PTT; Jeynes et al. [24] presented a numerical study of PTT on a two layers skin cancer model showing results that can be faced with experimental results. Manuchehrabadi and Zhu [25] presented a numerical PTT study for the tumor in the prostate. The prostate was simulated as an ellipsoid containing a tumor ball. The numerical study allowed to identify the treatment duration and irradiance power needed to remove the entirety of the tumor with limited damage to healthy surrounding tissues. Such in silico approaches give access to additional information inaccessible in experimental setups which could be used to improve the therapeutic protocol.
Following these approaches and with the aim of enlightening experimental results in a more complex media, a numerical study simulating the effect of PTT on a GBM murine phantom is presented in this paper. A synthetic phantom of a mouse head was implemented. The morphology (skull, skin, blood vessels network, gray and white matter) of the head was obtained from different atlases. Theses atlases are built with different imaging methods and are available online [26][27][28][29]. They were assembled with a control point registration method [30]. Additionally, to ease the interpretation of future experimental results, a thorough literature review on the physical parameters was performed to design this realistic phantom.
In this paper, we present a numerical phantom that can be used to simulate preclinical studies in the field of photothermal therapy for GBM. Orthotopic mice models are generally used to evaluate the efficacy and safety of innovative therapeutic strategies in the early development phases. The development of in silico models able to mimick the tumor microenvironment are highly recommended in order to screen therapeutics and limit the number of animal studies to be performed in the respect of the replacement, reduction and refinement (3Rs) rule [31]. Indeed, this numerical phantom could be used as a model to select the NP type and dose (thus reducing the number of in vivo studies to be performed) or which configuration of laser illumination or sensor is best suited (thus reducing the number of animals in the in vivo studies).
The purpose of the present study is i) to model the temperature evolution in the phantom and more precisely within the tumor during PTT, ii) to quantify the thermal damages caused to the tumor and the surrounding tissues and iii) to perform PA simulations at different stages of the PTT process and link the different signals obtained with the temperature evolution in the phantom in order to evaluate the possibility of monitoring the temperature through PA signals evolution.
To this end, Monte Carlo simulations were performed to model light propagation in the complex numerical phantom. The temperature evolution due to a continuous wave (CW) illumination was modeled with heat diffusion simulations, by solving the Pennes' bioheat equation [32], accurate to describe heat transfer in biological tissues. In the GBM phantom, NP were considered to locally increase the temperature at the tumor site. Following the thermal simulations, PA simulations were performed. PA signals depend on the Grüneisen parameter Γ. This parameter gathers the mechanical properties of the probed medium and undergoes changes when the temperature evolves. This dependency was implemented in the phantom, and its influence on the PA signals was studied. The results support the fact that PA signals can be used to monitor temperature in real time during PTT, provided a proper calibration is performed.
The governing optical and thermal equations are presented in Section 2.. The phantom geometry, its physical properties and the simulation methods are discussed in the Section 3.. In Section 4., the results of the different simulations performed are presented and discussed. Finally, thermal damages in the brain are quantified and PA signals measured during the PTT process are discussed.

Theoretical background
The important equations describing the multiphysical phenomena happening during a PTT and during the PA signal generation are presented in this section. All the principal parameters and the associated notations used in this paper are reported in Table 1.

Light propagation in tissues
Light propagation in a tissue illuminated by a CW source (e.g. laser) is described by the radiative transfer equation (RTE) [33]. Let r=(x, y, z) be the position, s a direction, t the time and Ω a solid angle, this equation is expressed as: (s · ∇ + µ a (r) + µ s (r))L(r, s, t) = µ s (r) ∫ 4π p(s, s')L(r, s', t)dΩ ′ + q(r, s, t).
(1) with g, the anisotropy coefficient defined as: p(s, s') represents the phase function. The Henyey Greenstein phase function was used for p. The fluence is obtained via:

Heat diffusion in biological tissues
Heat diffusion in biological tissues is a complex theory as it must take into account biological processes that regulate temperature inside tissues (convection, conduction, blood perfusion. . .). The temperature in a biological system subject to a heat source can be described with the Pennes' Bioheat equation [32]: Cooling of a biological system by the blood circulation is accounted for with the value of the perfusion rate w b of the blood in the tissue. The deposited heat rate Q H is calculated from the fluence Φ using: Temperature rises in tissues lead to thermal damages once the received thermal dose reaches a threshold. Sapareto and Dewey [34] proposed the cumulative equivalent minute at 43 • C (CEM43) to calculate the thermal dose. Following their work, a heating process realized for a certain duration t produces a thermal dose equal to: R is an experimentally determined value and depends on the tissue and species. This parameter R is related to the rate of cell survival in function of the temperature. When temperature rises, this survival rate decreases linearly but with two different slopes depending if temperature is below or above 43 • C. R is directly related to these slopes [34]. In the case of the mouse, R is usually set at 0.5 for T > 43 • C and at 0.25 for T < 43 • C [35,36] .

Photoacoustic equations and temperature monitoring
Using a pulsed laser satisfying the assumptions of thermal and stress confinements (few nanoseconds), the pressure rise p 0 generated at initial time in an irradiated tissue is equal to: The fluence Φ is accessed through the resolution of the RTE (Eq. (3)). Q PA represents the light's energy deposition rate. The Grüneisen parameter Γ defines the proportion of optical energy which will be transferred to the acoustic pressure and is temperature dependent, as the speed of sound v s , the thermal expansion coefficient β and the specific heat capacity C p depend on the temperature. The pressure created at initial time p 0 in a heterogeneous and acoustically absorbing medium propagates according to the following coupled equations [37]: with initial conditions: p(r, t) and u(r, t) are respectively the acoustic pressure and the acoustic particle velocity at position r and time t, ρ 0 represents the medium density at ambient temperature. Absorption and dispersion in the medium are integrated in the equation through τ and η coefficients. They are given by: α 0 is obtained, in biological tissues, through the use of a power law: where f is the frequency of the acoustic pressure waves created by the pulsed laser. Temperature variations in a biological tissue lead to pressure variations as well. In the case of a linearly varying Grüneisen parameter, a linear relation between pressure and temperature can be expressed [19]: where ∆p and ∆Γ represent the variations of the pressure and the Grüneisen parameter, respectively, when the temperature rises by ∆T and γ is a constant depending on the considered tissue. p i and Γ i represent respectively the initial pressure and initial Grüneisen parameter of the medium before the PTT begins. The calibration factor γ has to be determined experimentally.

Building a realistic numerical phantom of a mouse head
The phantom in which the PTT was simulated is presented in this section: here, the geometry and parameters used, as well as the simulation methods and configurations are presented.

Phantom geometry
The mouse head model used in this paper was gathered from different sources and assembled using a coregistration method. The brain atlas of an adult C57BL/6J mouse [26,29], the skull atlas of an adult 129S1/SvIMJ mouse [27] and the cerebral vasculature atlas of a CBA mouse [28] were co-registered to obtain the 3D model, as shown in Fig. 1. In total, the head's dimension was approximately 15 mm × 10 mm × 10 mm digitalized into 512 × 512 × 512 voxels, with a discretization size of 60 µm in every dimension. This discretization corresponds to the skull atlas resolution we used to create the phantom [27]. The used brain atlas was segmented into gray and white matter, allowing more accurate representation of the morphology. Skin was added manually in the form of an additional layer on top of the skull whose thickness was set to 180 µm, which is the closest to values from the literature we can reach with our grid discretization method [38,39]. Moreover, the cerebrospinal fluid (CSF) was also assigned manually between the skull and the brain. The tumor position and its contour was indicated by the dashed purple line in Fig. 1. It was placed on the right hemisphere of the mouse brain and modeled by a simple 3 mm 3 ball. This value was chosen to represent a 10 days GBM and is within the order of values found in the literature [40][41][42]. The phantom model is surrounded by water, as it is a good acoustic impedance matching medium for biological tissues [19,22]. Effect of using additional thermal agents as NP was numerically implemented by an increase of the absorption coefficient of the tumor at the source's wavelength. Four different cases were considered: • "No tumor": the case of a healthy brain, the phantom contains no tumor.
• "No NP": the case of a tumorous brain without added NP, so that the absorption coefficient is due to the local increase of vascularity within the tumor.
• "µ aT × 2": the case of a tumorous brain with added NP, so that the absorption coefficient of the tumor is multiplied by 2.
• "µ aT × 5": the case of a tumorous brain with added NP, so that the absorption coefficient of the tumor is multiplied by 5.
The first case considered is the case of a healthy brain which serves as a reference for comparison with the three other cases which present a tumorous brain. Three types of coupled simulations were performed: light propagation simulations were performed first to compute the source terms Q H and Q PA needed in the heat diffusion and PA simulations.

Source and sensor configuration
A 320 mW top-hat CW-laser was used for the PTT heating simulations. The diameter of the laser beam is 3 mm and its wavelength is 808 nm. This wavelength was chosen in the therapeutic window, in which the penetration depth is the highest. Power per area of the PTT laser we used is approximately 4.5 W.cm −2 , which is a typical value for laser power per area in PTT: in the literature, values for such laser range from 0.1 to 16 W.cm −2 [8,[44][45][46]. It should be pointed out that the power per area used far exceeds the maximum permitted exposure value for skin illumination, which is approximately 329 mW.cm −2 at the wavelength used [47]. Most PTT studies go beyond this safety limit to achieve an efficient thermal treatment with limited and reversible damages, as it can be seen in the reference [8], and no consensus actually exists on what is the optimal laser power acceptable for an efficient and relatively safe PTT. The purpose of this numerical tool is precisely to help testing various parameters to define the most appropriate protocol.
A 50 µJ pulsed top-hat laser with a pulse duration of 5 ns is used for PA simulations, in order to be consistent with the thermal and stress confinement conditions. The diameter of the laser beam and its wavelength is the same as the PTT laser: 3 mm and 808 nm respectively. These chosen values for the pulsed laser define a power per area of approximately 0.7 mJ.cm −2 and therefore respect maximum permitted exposure value for skin illumination [47].
In both heat and optical simulations, the laser is illuminating the skin surface above the tumor position. As the source configuration and wavelength are the same for both the PTT and the PA signal generation, we calculated a normalized light deposition, which correspond to light deposited by a laser of power 1 W. We then had to multiply this normalized light deposition by the power of the PA laser and the CW-laser to obtain respectively Q PA and Q H .
The sensor for the PA simulations was modeled to mimic a curved sensor presented in Fig. 1. This sensor has a centrally pierced hole of 0.8 cm in diameter and is focused at 2 cm below its surface. In the configuration chosen, the center of the tumor is defined as the focal point. These source and sensor configurations are taken from reference [22]. The simulated sensor is supposed perfect.

Simulation methods
The simulations' flow chart is presented in Fig. 2. Monte Carlo simulations were first performed to access information on the fluence rate and absorbed power using the Monte Carlo program MCmatlab [48]. This program allows simulations in complex media in a voxel-designed simulation grid. Therefore, Monte Carlo simulations allow us to calculate the source terms for both thermal and photoacoustic simulations [49]. These simulations were performed at a wavelength of 808 nm. 10 8 photons were launched for each Monte Carlo simulations.
The versatile k-wave toolbox [50] was used for heat diffusion and PA simulations. For the thermal simulations, the toolbox allows also to solve the Penne's bioheat equation in the 3D phantom. For PA simulations, it solves the first-order coupled acoustic propagation equations for a highly heterogeneous medium with acoustic attenuation. Both computations use a k-space pseudospectral method: spatial gradients are calculated using the Fourier collocation spectral method and temporal gradients are calculated using a k-space corrected finite difference scheme. The numerical phantom is heated continuously for 400 s with a CW-laser. The initial temperature of the phantom is set to 37 • C in every medium. Blood perfusion is taken into account. Blood ambient temperature in the brain is set to 37 • C. Considering the evolution of the temperature in heat diffusion experiments, a small time step of 0.1 s was chosen for the first 100 s of heating, as the temperature slope varies the most in this interval. This time step is then set at a much wider value of 10 s for the remaining heating duration, when the slope is weaker. Heat diffusion simulations were performed in the four scenarios described in Sec. 3.1 to obtain the temperature evolution during the 400 s PTT. It must be noted that the simulation grid used for thermal simulation was extended: the k-space methods impose a periodic boundary condition for the simulation grid. Therefore, the initial 512×512×512 voxels grid was reduced to a 256×256×256 voxels grid with interpolation methods and placed in a 640×640×640 voxels grid. Boundary conditions of continuity between each medium in the heat transfer simulations were imposed.
Thermal damages were then examined for all cases. They were obtained by using the thermal threshold values found in the literature [36] and summarized in Table 2. This threshold presents the limit amount of thermal dose a tissue can withstand before sustaining damage and is a value which depends on the tissue considered. More specifically, the table presents values for which irreversible damages were inflicted to the corresponding tissues [36]. The value for the thermal threshold of the tumor was assumed here to be the same as the brain as a suitable value was not found in the literature. Therefore, in the continuation of the paper, a tissue will be considered burnt, once the thermal dose it received reaches the corresponding thermal threshold.
As the initial pressure depends on temperature and fluence, the results of light propagation and thermal diffusion simulations were used to compute the PA signals measured by the acoustic sensor. PA simulations were performed at different stages of the heating process in order to observe the effect of temperature increase on the amplitude of the PA signals during the PTT.

Multiphysical parameters of the phantom
In order to simulate the PTT process, multi-physical parameters of the phantom need to be defined. On this end, optical, thermal and acoustic parameters have been gathered from multiple sources in the literature. The optical properties of the different tissues in the phantom are presented in Table 3. The optical properties of the tumor are from reference [57]. Optical parameters of a 10 days GBM were chosen, corresponding to a classical tumor growth prior imaging. When specific parameters for the mouse are not found, values from human tissues were considered.
Thermal and acoustic parameters needed for the simulations are presented in Table 4.  Values presented are mostly from human origins. When not from human, origin is specified. CSF's acoustic attenuation parameters and thermal conductivity were assumed the same as water. Tumor's acoustic attenuation parameters are supposed the same as gray matter.
To model the dependence on temperature of the Grüneisen parameter, an important effort has been made to find temperature dependent values of the thermal expansion β, the speed of sound v s and the specific heat capacity C p for each tissue of the model. Expressions of v s and β as a function of the temperature for most of the phantom's tissues were collected from the literature [67][68][69][70][71][72]. The obtained temperature depending curves for v s and β are plotted on Figs. 3 and 4. The constant values for β are presented in Table 5. C p was supposed non dependant on the temperature for the simulations. No temperature dependant parameter was found for the skin, its Grüneisen parameter was therefore supposed constant in the simulations. However, as the skin thickness in the phantom is relatively small, the effect on the acoustic signal is expected to be small. For most of the tissues considered in this study, the expansion coefficient β was supposed constant. However, this parameter is the main source of temperature sensitivity in PA [19]. Nevertheless, the assumption of a constant value was made in the simulations and it is noteworthy to mention that this is a strong approximation in this model and, therefore, the temperature dependency of the Grüneisen parameter is misestimated in this study, but could be more precisely taken into account with new measured values.   The speed of sound evolution of the tumor is calculated using reference [67]. This reference paper presents a comparison of the speed of sound between a healthy and a tumorous liver, showing a decrease of approximately 1.4% of the ultrasounds speed in the tumorous liver. We assumed a similar behavior in the brain tumor: the speed of sound in the tumor was supposed equal to the one in the brain decreased by 1.4%.
It is well known that non linear effect of temperature on the Grüneisen parameter, but also the optical parameters, can occur. These effects appear because of several different biological processes such as coagulation or dehydration [73]. We limited our consideration of the dependencies with temperature in the case of the Grüneisen parameter. However, it is important to note that we did not take into account the non-linear effect on the optical properties in the present study. These effects may be implemented in future improvements of the phantom.
References for the perfusion rate value for the mouse were taken from [64,65]. Blood perfusion of the brain, the skin and the tumor have been shown to undergo significant change with the temperature change in several studies [74][75][76]. However we supposed them constant in the framework of this study. White matter, gray matter and blood share the same blood perfusion value as a mean value for the entire brain has been set [64]. The tumor's perfusion rate was determined using reference [66]. The paper presents blood flow in brain tumor in mice as a function of their sizes. Considering the tumor's volume, a value of 0.0038 s −1 has been deduced.
Acoustic attenuation in each tissue is expressed in the form of a power law. Values for the acoustic attenuation prefactor α 0 and power b were found [58,[60][61][62][63]. The k-wave toolbox allows for a unique acoustic attenuation power b in the computational grid. Therefore, b was set at 1.21, which is the value for gray matter, the most present medium in the phantom, excluding water. The acoustic attenuation prefactors α 0 of the other media were then changed according to the new value of b and chosen by minimizing the L1 norm between the new power law, and the original one.

Monte Carlo simulations
Heat deposition rates Q H obtained from Monte Carlo simulations are presented in Fig. 5. These results serve as heating source in the heat diffusion simulations presented hereafter.

Heat diffusion simulations
Temperature evolution in the phantom can be seen in Fig. 6 in the four considered cases. Temperature curves are extracted at a position corresponding to the surface of the tumor (position indicated in this same figure). Visualization of the damaged area is presented in Fig. 7. The figure presents damaged volumes in the tumor and its surrounding at different times of the heating process. Damages in tumorous and healthy tissues have been evaluated: Fig. 8 shows the percentage of the tumor burnt and the volume of surrounding healthy tissues burnt during the process.

Photoacoustic simulations
PA simulations were performed every 20 s from the beginning of the heating process up until 100 s of heating. A final simulation was done after 400 s of heating. In all considered cases, the whole tumor is burnt at this time. Following the heat diffusion and light propagation simulations, the initial pressure p 0 generated by the laser can be calculated with Eq. (7). For the four cases considered, the initial pressure maps when the phantom was at 37 • C (before the beginning of the PTT) are presented in Fig. 9. Figure 10 shows the obtained PA signals associated with these initial pressure maps. The impact of the tumor on the acoustic signals can be clearly seen and is pointed out on the figure. We can verify that these peaks are the tumor's response by multiplying their time position (approximately 12.2 µs) by the mean speed of sound in the head (approximately 1529 m.s −1 ). The result yields a value of 1.87 cm which corresponds to the tumor's surface (see Fig. 9). In the figure, the top right figure represents the signal after filtering with a gaussian filter of central frequency 2.1 MHz and bandwith ratio 80%. We performed this operation to give an insight on what a realistic signal would look like. Once the PTT begins, temperature rises in the phantom. As temperature rises, the amplitudes of these signals vary due to the evolution of the Grüneisen parameter. PA signals as temperature increases for the case "µ aT × 5" are presented in Fig. 10. As the speed of sound also varies with temperature during the PTT, a shift in time of the signals can be observed. Differential signals ∆p are presented for all three tumorous cases in Fig. 11, showing a clearer representation of the influence of the PTT duration on pressure variations. A time-shift compensation introduced by the temperature change in the phantom was taken into account for the post-processing of the signals.
We computed the ratios ∆p/p i of the tumor's response for the three tumorous cases during the PTT in Fig. 12.
Linking the signal evolution with the temperature rise can be done with the Eq. (14). The value of the calibration factor γ was calculated in each studied case then summarized in Table 6. It represents the slope of the linear fit in Fig. 12 (right).
It is worth noting, that the value of γ is totally dependent on the properties of the medium and configuration of the simulations. By considering a linear behavior of the tumor's PA signals as the temperature rise, temperature can be recovered from γ values in the table. Results of these   11. Differential signal ∆p(t) as PTT process progress for the three tumorous cases. Legend show how much time of PTT has passed in order to obtain the associated differential. The tumor's influence on the signal is indicated on the "µ aT × 5" case.  recovered temperatures are presented in Fig. 13 for the 3 tumorous cases and compared to the actual evolution of the temperature. Additionally, monitoring of the temperature was performed in the blood vessel above the tumor for comparison. Both positions are specified in the figure.

Discussion
The numerical study developed in this paper had for objectives 1) to model the temperature evolution in the phantom, and more precisely in the tumor, 2) to quantify the thermal damages caused by the PTT and 3) to link the PA signals with the temperature evolution in the tumor in order to evaluate the possibility of using the PA technique as a temperature monitoring method.
In keeping with these objectives, simulation results are discussed in this section.

Temperature evolution in the phantom
The effect of the presence of the tumor can be clearly seen in Fig. 6: the temperature increases with a steepest slope at early times allowing to reach a higher temperature plateau. Indeed, the temperature in the "No tumor" case reaches 47 • C vs 51 • C for the tumorous case without NP ("No NP"). This is naturally due to the fact that, considering the values found in the literature, the tumor is more absorbing than the brain tissues (white and gray matter) at 808 nm. This behavior is enhanced by the use of NP in the tumorous brain that increase the absorption properties of the tumor: after 400 s of PTT, the total increase in the "No NP" case is approximately of 14 • C vs 19 • C for the "µ aT × 5" case.

Evaluation of the thermal damages
As PTT progresses, the thermal dose received by certain parts of the phantom reaches their thermal threshold, producing irreversible damages. From Fig. 7, the progression of the thermal damages can be clearly observed: the tumor's surface is burnt first, followed by the close surrounding area, this area widening with time. Indeed, the temperature at the surface of the tumor increases the fastest (see Fig. 6). Because of its strong absorption properties, the majority of the light coming from the laser is absorbed by the surface of the tumor, as shown in Fig. 5.
Consequently, at this position, the thermal threshold is reached early, and as heat diffuses, the temperature in the surroundings of the surface of the tumor rises as well, and reaches the thermal threshold of the peripheral tissues, which leads to the destruction of the lower part of the tumor but also to the damage of the healthy brain tissues above the tumor (see Fig. 7). In order to deliver a sufficient thermal dose to destroy the whole tumor, prolonged exposure is needed, which inevitably causes irreversible damage to healthy tissues. The quantification of the damages done to the tissues has been performed and enables to understand the effect of the NP and their possible advantages in the PTT. The NP can be considered advantageous if, compared to the case without NP, for the same amount of damage on the tumor, the volume of healthy tissues damaged is reduced. For example, we can compare the two cases "No NP" and "µ aT × 5", with an objective to burn 90% of the tumor (see Fig. 8). In the "No NP" case, the phantom needs to be heated for approximately 228 s to reach this objective, whilst for the case "µ aT × 5", the duration of heating is 139 s. The corresponding volumes of burnt healthy tissues in the "No NP" and "µ aT × 5" cases are respectively 7.27 mm 3 and 6.76 mm 3 . Therefore, the added NP show an advantage as they slightly reduce the total volume of unwanted damages. This conclusion remains the same for a lower portion of tumor burnt: If the same reasoning is performed for an objective of 50% of tumor burnt, the collateral damages to healthy tissues are reduced when NP are added, with approximately 1.52 mm 3 of burnt healthy tissues for the "µ aT × 5" case and 2 mm 3 for the "No NP" case. A partial destruction of the tumor can be wanted in a PTT with repeated therapy sessions. In this case, the tumor is burnt gradually thanks to multiple heating phases interspersed with resting phases, in order to avoid any damages to surrounding tissues. In conclusion, the use of NP allows for a reduction of the amount of unwanted damage to healthy tissues surrounding the tumor.
It is important to note that some biological mechanisms were not considered to simplify the simulations. Firstly, we assumed the optical properties of the tumor did not change as the therapy advances but the burnt tissues naturally undergo changes in their optical properties. Reference [80] measures variations in the absorption coefficient in necrotic cells, meaning the dynamic of the temperature evolution could change as PTT is performed. Secondly, the values used for defining the thermal thresholds in the phantom might not be accurate to represent the actual threshold needed for a complete destruction of tumorous cells. The values chosen and presented in Table 2 were values for irreversible damages to the tumor. As irreversible damages does not necessarily lead to a complete necrosis of the cells, the values chosen may be underestimated. However, as we didn't find appropriate and complete reference in literature, these values were used. Thirdly, the regrowth properties of the tumor were not considered. Indeed, if the tumor is not burnt completely during the PTT, a quick regrowth occurs and the tumor recovers from the damages done. It is important to note that, in the case of a repeated PTT process with multiple heating and resting phases, this biological reaction of the tumor can compromise the therapy. Finally, only the effect of the NP on absorption was considered: scattering properties changes induced by the NP delivery were not considered and might change the obtained results [81].

PA measurements as a temperature monitoring method
In order to link temperature and the PA signals received during the PTT, the tumor's effect on the PA signals need to be identified (see Fig. 9). In the PA signals received in the tumorous cases and shown in Fig. 10, three peaks are present. The third one can be clearly identified as the tumor's response. As expected, this peak increases with the amount of added NP. The two other peaks can be identified: the first peak corresponds to the skull and skin's conjoint response and the second peak corresponds to the blood vessels just above the position of the tumor, which gives a strong response because of its strong absorption coefficient at 808 nm. Note that in the non-tumorous brain case, a peak is also seen near the tumor response's location in the tumorous cases and is due to a blood vessel present at the position the tumor is placed (see Fig. 10).
As PTT is performed, the temperature changes and the consequences in the PA signals is clearly shown (see Fig. 10 for the "µ aT × 5" case). With the increase of temperature, the amplitude of the second peak in the signals rises as PTT advances, which is directly linked to the increase of the Grüneisen parameter of the blood with temperature. However, despite the fact that the Grüneisen parameter of the tumor increases with temperature, the maximum amplitude of its response does not show a constant growth. This phenomenon is presented in the "µ aT × 5" case in the Fig. 10 but is also clearly present in the "No NP" and "µ aT × 2" cases. To understand this phenomenon, it is important to note that the third peak of the PA signal is not solely the tumor's response: the contributions from its close environment alters this response. Therefore the PA signals are the results of different media responding together, interfering and inter-meddling with each other. In the considered geometry, the blood vessels positioned above the tumor is the principal cause of interference of the tumor's acoustic response. This blood medium casts a shadow on the lower areas and add therefore a negative contribution on the tumor's acoustic response. This contribution also varies with temperature, leading to the complex variations seen in the figure. Differential signals at different times of the PTT presented in Fig. 11 allow us to see the variations of the tumor's response in the three considered cases. These variations are not clearly visible in the "No NP" and the "µ aT × 2" cases, as they are small compared to the variations of the blood's response, which represent the highest peak of these differential signals. However, in the "µ aT × 5" case, the tumor's response variations can be identified and is presented in the figure. In this case, the peak representing the tumor's response variations is shown to take negative values, which is the consequence of the negative contribution of the blood vessels above the tumor in the considered geometry.
In order to study the feasibility of the PA temperature monitoring, variations of the tumor's response with temperature has been quantified using the ratio ∆p/p i , presented in Fig. 12. In the "µ aT × 2" and "µ aT × 5" cases, the curves show a similar behavior as temperature rises in the tumor. In the "No NP" case, the curve reaches negative values. This is also due to the blood vessel casting a shadow on the tumor, adding therefore a negative contribution on the tumor's response and leading to negative values of ∆p/p i for the "No NP" case. However, the variations of this curve match the variations of the other two cases: it increases with the heating duration which is a fact that agrees with the increase of the Grüneisen parameter of the tumor with the temperature. In the right graph of Fig. 12, the linear fits match with the actual behaviors of the curves in the "µ aT × 2" and "µ aT × 5" cases but not in the "No NP" case. This fact impacts directly the efficiency of the temperature monitoring in the tumor: as γ reflects the slope of the linear fit, the better the linear fit match the actual behavior of the curves, the better the accuracy of the temperature monitoring will be. Therefore, tumor's monitoring obtained for the "µ aT × 5" case show accurate results, with approximation close to the original evolution. A small drift of the approximation can be seen for long heating duration, showing a drift of 2 • C "µ aT × 5". In the "µ aT × 2" case, the result still agrees well with the actual variation despite a larger drift of 4 • C after 400 s. However, for the "No NP case", the approximated ∆T obtained shows significant drifts from the actual evolution, showing that without NP, in the configuration we placed the tumor, monitoring of the temperature is harder as the pressure evolution due to the temperature rise is not linear. Non linear perturbations occur in the signals due to environment influence (e.g. blood's shadow casting). They hinder the assumed linear hypothesis on the relationship between pressure rise and temperature rise. These non linear perturbations are strong in the "No NP" case, however, with added NP, the perturbations are undermined by the amplifying effect of the NP: the linear relationship assumption holds in these cases. In order to test the possibility to monitor the surrounding of the tumor at least, the monitoring of the temperature of the blood vessel above the tumor was performed. Results obtained for the blood show an accurate monitoring of the temperature in all cases (see Fig. 13). Therefore, monitoring temperature at the tumor position without NP is not always possible, however the temperature monitoring of its close environment can be used as a landmark for understanding the temperature evolution in the tumor itself, as the temperature monitoring of the blood vessel above the tumor shows accurate results in all considered cases.
In conclusion, NP allows for a better monitoring of the temperature in the tumor. Without them, the tumor response can be hidden behind strong absorbers such as the blood. In the present case, the blood alters the response of the tumor in the PA signals, and consequently, prevents it from being used for accurate temperature monitoring. Indeed, in the "No NP" case, the temperature monitoring shows important drift compared to the actual temperature evolution. Adding NP allows for a stronger response of the tumor. The stronger the contrast between the tumor and its surroundings, the more accurate the monitoring of the temperature on the tumor will be. However, in comparison, temperature monitoring of the blood vessel above the tumor was possible for all studied cases, showing that monitoring for upper media such as the blood can be possible and used as another way to gather information on temperature in the tumor.

Nanoparticles
NP effect in this numerical phantom has been simulated by an increase of the absorption coefficient of the tumor. Increased values used in this paper are exaggerated as a doubled absorption of the initial tumor means an absorption coefficient of 9.2 cm −1 at 808 nm and a five-fold increase reach 23 cm −1 . These absorption levels are not reachable with actual NP without reaching a cytotoxic level. We aimed with this study to have a better understanding of the needed absorption of tumor to define optimal PTT and temperature monitoring, thus the unrealistic values used.

Conclusion
This paper presents a numerical study performed on a realistic phantom of a mouse head showing the complexity of realizing PTT with and without NP. Temperature in the phantom, and more precisely, in the tumor during the PTT was obtained thanks to heat diffusion simulations. Thermal damages caused were evaluated and NP effect has been studied in the numerical phantom. Added NP in the tumor leads to higher temperature reached in the tumor. Thermal damages are caused faster with NP, but surrounding healthy tissues burn also faster. The study shows that adding NP in the tumor is slightly advantageous as damages to the surrounding healthy tissues are reduced. Temperature rises and PA signal were linked, showing good monitoring capability for the "µ aT ×2" and "µ aT × 5" cases. However, in the "No NP" case, monitoring accurately temperature was not possible with the value of γ calculated. Therefore, in this numerical study, NP are needed for temperature monitoring, as without NP, the complex geometry and strong absorbers surrounding the tumor prevent accurate monitoring. Mainly, we presented here a precious numerical platform that has been used to perform a study on PTT with PA technique as a real time temperature monitoring method in mice. Specific treatment configurations could be simulated accurately and the different key parameters of the photothermal therapy could be optimized to performed an effective treatment with controlled thermal damages. This numerical phantom could also be used as a tool to test NP and their efficiency on the therapy. However, an important point would be to improve the knowledge of some tissues properties and biological mechanisms. Constant values were used in this paper for the dependency with temperature of biological tissues for some media because of the lack of references. Information on accurate thermal threshold for necrosis were also lacking and some natural mechanisms were not implemented such as the modification of the optical properties and blood perfusion with the progression of the PTT. Future researches may improve the model.