<figure><img src="image_1.png" /></figure>

HAL
open science

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

Antoine Capart, Khaled Metwally, Chiara Bastiancich, Anabela da Silva

► To cite this version:

Antoine Capart, Khaled Metwally, Chiara Bastiancich, Anabela da Silva. Multiphysical numerical study of photothermal therapy of glioblastoma with photoacoustic temperature monitoring in a mouse head. Biomedical optics express, 2022, 13 (3), pp.1202. 10.1364/BOE.444193. hal-03588832

HAL Id: hal-03588832

https://hal.science/hal-03588832

Submitted on 25 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
To be published in *Biomedical Optics Express*:

**Title:** Multiphysical numerical study of photothermal therapy of the glioblastoma with photoacoustic temperature monitoring in the mouse head

**Authors:** Antoine Capart, Anabela Da Silva, Chiara Bastiancich, khaled metwally

**Accepted:** 03 January 22

**Posted:** 05 January 22

**DOI:** https://doi.org/10.1364/BOE.444193

© 2022 Optica under the terms of the [Optica Open Access Publishing Agreement](#)

OPTICA
PUBLISHING GROUP
Formerly OSA
# Multiphysical numerical study of photothermal therapy of the glioblastoma with photoacoustic temperature monitoring in the mouse head

ANTOINE CAPART<fnref n="1" />,*, KHALED METWALLY<fnref n="1" />,<fnref n="2" />, CHIARA BASTIANCICH<fnref n="3" /> AND ANABELA DA SILVA<fnref n="1" />

**Abstract:** 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.

© 2022 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. 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–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, 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 treatment (including surgery, radiotherapy and chemotherapy) the median overall

<footnote n="1">¹Aix Marseille Univ, CNRS, Centrale Marseille, Institut Fresnel, Marseille, France</footnote>

<footnote n="2">²Aix Marseille Univ, CNRS, Centrale Marseille, LMA, Marseille, France</footnote>

<footnote n="3">³Aix Marseille Univ, CNRS, UMR 7051 Institute Neurophysiopathol, INP, CNRS, Aix-Marseille University, 13005 Marseille, France</footnote>

<footnote>*antoine.capart@fresnel.fr</footnote>
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–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. These atlases are built with different imaging methods and are available online [26–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.

## 2. 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.

### 2.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 $\mathbf{r}=(x,y,z)$ be the position, $\mathbf{s}$ a direction, $t$ the time and $\Omega$ a solid angle, this equation is expressed as:

$$ (\mathbf{s} \cdot \nabla + \mu_a(\mathbf{r}) + \mu_s(\mathbf{r})) L(\mathbf{r}, \mathbf{s}, t) = \mu_s(\mathbf{r}) \int_{4\pi} p(\mathbf{s}, \mathbf{s}') L(\mathbf{r}, \mathbf{s}', t) d\Omega' + q(\mathbf{r}, \mathbf{s}, t). \quad (1) $$

with g, the anisotropy coefficient defined as:

$$ g = \int_{4\pi} p(\mathbf{s}, \mathbf{s}') \mathbf{s} \cdot \mathbf{s}' d\Omega. \quad (2) $$

$p(\mathbf{s}, \mathbf{s}')$ represents the phase function. The Henyey Greenstein phase function was used for $p$.

The fluence is obtained via:

$$ \Phi(\mathbf{r}) = \int_{4\pi} L(\mathbf{r}, \mathbf{s}) d\Omega. \quad (3) $$

### 2.2. 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...).
<table><caption>Table 1. Glossary of notations</caption><thead><tr><th>Symbol</th><th>Description</th><th>Unit</th></tr></thead><tbody><tr><td>μ<sub>a</sub>/μ<sub>s</sub></td><td>Absorption/Scattering coefficient</td><td>cm<sup>-1</sup></td></tr><tr><td>g</td><td>Anisotropy coefficient</td><td>dimensionless</td></tr><tr><td>L</td><td>Specific intensity</td><td>W.cm<sup>-2</sup>.sr<sup>-1</sup></td></tr><tr><td>q</td><td>Source term</td><td>W.cm<sup>-3</sup>.sr<sup>-1</sup></td></tr><tr><td>p</td><td>Phase function</td><td>dimensionless</td></tr><tr><td>Φ</td><td>Fluence</td><td>W.cm<sup>-2</sup></td></tr><tr><td>β</td><td>Volumetric thermal expansion coefficient</td><td>K<sup>-1</sup></td></tr><tr><td>C<sub>p</sub>(Cp<sub>b</sub>)</td><td>Specific heat capacity (of the blood)</td><td>J.kg<sup>-1</sup>.K<sup>-1</sup></td></tr><tr><td>ρ(ρ<sub>b</sub>)</td><td>Density (of the blood)</td><td>kg.m<sup>-3</sup></td></tr><tr><td>T(T<sub>b</sub>)</td><td>Temperature (of the blood)</td><td>°C</td></tr><tr><td>w<sub>b</sub></td><td>Blood perfusion rate</td><td>s<sup>-1</sup></td></tr><tr><td>Q<sub>H</sub></td><td>Volumetric rate of heat deposition</td><td>W.cm<sup>-3</sup></td></tr><tr><td>Q<sub>PA</sub></td><td>Volumetric rate of light's power deposition</td><td>W.cm<sup>-3</sup></td></tr><tr><td>Γ</td><td>Grüneisen parameter</td><td>dimensionless</td></tr><tr><td>k</td><td>Thermal conductivity</td><td>W.m<sup>-1</sup>.K<sup>-1</sup></td></tr><tr><td>α<sub>0</sub></td><td>Attenuation prefactor</td><td>dB.MHz<sup>-1</sup>.cm<sup>-1</sup></td></tr><tr><td>b</td><td>Attenuation power</td><td>dimensionless</td></tr><tr><td>α</td><td>Power law acoustic attenuation</td><td>dB.cm<sup>-1</sup></td></tr><tr><td>ν<sub>s</sub></td><td>Speed of sound in the medium</td><td>m.s<sup>-1</sup></td></tr></tbody></table>

The temperature in a biological system subject to a heat source can be described with the Pennes' Bioheat equation [32]:

$$ \rho C_p \frac{\partial T}{\partial t} = \nabla \cdot (k \nabla T) - w_b \rho C_{p,b} (T - T_a) + Q_H. \quad (4) $$

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 $\Phi$ using:

$$ Q_H = \mu_a \Phi(\mu_a, \mu_s, g). \quad (5) $$

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:

$$ CEM43 = \int_{0}^{t} R^{43-T} dt. \quad (6) $$
$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].

## 2.3. 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:

$$p_0 = \Gamma \mu_a \Phi(\mu_a, \mu_s, g) = \Gamma Q_{PA} \quad \text{with} \quad \Gamma = \frac{\beta v_s^2}{C_p}. \qquad (7)$$

The fluence $\Phi$ is accessed through the resolution of the RTE (Eq.(3)). $Q_{PA}$ represents the light's energy deposition rate. The Grüneisen parameter $\Gamma$ 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 $\beta$ 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]:

$$\frac{\partial u(\mathbf{r}, t)}{\partial t} = -\frac{1}{\rho_0(\mathbf{r})} \nabla p(\mathbf{r}, t); \qquad (8)$$

$$\frac{\partial \rho(\mathbf{r}, t)}{\partial t} = -\rho_0(\mathbf{r}) \nabla \cdot u(\mathbf{r}, t); \qquad (9)$$

$$p(\mathbf{r}, t) = v_s^2 \rho(\mathbf{r}, t) [1 - \tau(\mathbf{r}) \frac{\partial}{\partial t} (-\nabla^2)^{\frac{b}{2}-1} - \eta(\mathbf{r}) (-\nabla^2)^{\frac{b+1}{2}-1}], \qquad (10)$$

with initial conditions:

$$p(\mathbf{r}, t=0) = p_0(\mathbf{r}), \quad u(\mathbf{r}, t=0) = 0. \qquad (11)$$

$p(\mathbf{r}, t)$ and $u(\mathbf{r}, t)$ are respectively the acoustic pressure and the acoustic particle velocity at position $\mathbf{r}$ and time $t$, $p_0$ represents the medium density at ambient temperature. Absorption and dispersion in the medium are integrated in the equation through $\tau$ and $\eta$ coefficients. They are given by:

$$\tau(\mathbf{r}) = -2\alpha_0 v_s (\mathbf{r})^{b-1}, \quad \eta(\mathbf{r}) = 2\alpha_0 v_s (\mathbf{r})^b \tan(\frac{\pi b}{2}). \qquad (12)$$

$\alpha_0$ is obtained, in biological tissues, through the use of a power law:

$$\alpha_0 = \alpha f^b, \qquad (13)$$

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]:

$$\frac{\Delta p}{p_i} = \gamma \Delta T = \frac{\Delta \Gamma}{\Gamma_i}, \qquad (14)$$

where $\Delta p$ and $\Delta \Gamma$ represent the variations of the pressure and the Grüneisen parameter, respectively, when the temperature rises by $\Delta T$ and $\gamma$ is a constant depending on the considered tissue. $p_i$ and $\Gamma_i$ represent respectively the initial pressure and initial Grüneisen parameter of the medium before the PTT begins. The calibration factor $\gamma$ has to be determined experimentally.
### 3. 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.

#### 3.1. 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³ ball. This value was chosen to represent a 10 days GBM and is within the order of values found in the literature [40–42]. The phantom model is surrounded by water, as it is a good acoustic impedance matching medium for biological tissues [19, 22].

<figure><img src="image_3.png" /><figcaption>Fig. 1. Phantom of a mouse head. Top: 2D cortical (left) and sagittal (right) views of the phantom (Tumor's contour is indicated by the dashed purple line). Bottom: source and sensor configurations (left): sensor position is represented in black. Laser's illumination configuration is in red. The bottom right panel shows an image of the curved sensor used (top) as well as a 3D representation (bottom). Reused from [43].</figcaption></figure>

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.
* "$\mu_{aT} \times 2$": the case of a tumorous brain with added NP, so that the absorption coefficient of the tumor is multiplied by 2.

* "$\mu_{aT} \times 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.

## 3.2. 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⁻², 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⁻² [8, 44–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⁻² 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⁻² 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.

## 3.3. 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.

<figure><img src="image_4.png" /><figcaption>Fig. 2. Flow chart of the computed simulations.</figcaption></figure>

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.

<table><caption>Table 2. Thermal thresholds for irreversible damage used in each medium of the phantom [36]. GM: Gray matter, WM: White matter</caption><thead><tr><th></th><th>Water</th><th>Skin</th><th>Skull</th><th>CSF</th><th>Blood</th><th>GM</th><th>WM</th><th>Tumor</th></tr></thead><tbody><tr><th>Thermal threshold (in CEM43 unit)</th><td>/</td><td>40</td><td>16</td><td>/</td><td>/</td><td>10</td><td>10</td><td>10</td></tr><tr><th>Species</th><td>/</td><td>Mouse</td><td>Rabbit</td><td>/</td><td>/</td><td>Dog</td><td>Dog</td><td>Dog</td></tr></tbody></table>

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.
### 3.4. 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.

<table><caption>Table 3. Phantom optical parameters [29, 51–57], at 808 nm. When not found, the anisotropy factor g was supposed equal to 0.9.</caption><thead><tr><th></th><th>μ<sub>a</sub>(cm<sup>-1</sup>)</th><th>μ<sub>s</sub>(cm<sup>-1</sup>)</th><th>g</th><th>Ref.</th><th>Species</th></tr></thead><tbody><tr><td>Water</td><td>0.02</td><td>0</td><td>1</td><td>[51]</td><td>/</td></tr><tr><td>Skin</td><td>1</td><td>60</td><td>0.9</td><td>[52]</td><td>Mouse</td></tr><tr><td>Skull</td><td>0.16</td><td>160</td><td>0.9</td><td>[53, 54]</td><td>Human</td></tr><tr><td>CSF</td><td>0.02</td><td>1.5</td><td>1</td><td>[29]</td><td>assumed same as water</td></tr><tr><td>Gray matter</td><td>0.27</td><td>74</td><td>0.9</td><td>[55]</td><td>Human</td></tr><tr><td>White matter</td><td>0.88</td><td>360</td><td>0.87</td><td>[55]</td><td>Human</td></tr><tr><td>Blood</td><td>3.9</td><td>765</td><td>0.98</td><td>[56]</td><td>Human</td></tr><tr><td>Tumor</td><td>4.6 (No NP)<br>9.2 (μ<sub>aT</sub> × 2)<br>23 (μ<sub>aT</sub> × 5)</td><td>219</td><td>0.81</td><td>[57]</td><td>Rat</td></tr></tbody></table>

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.

<table><thead><tr><th rowspan="2"></th><th rowspan="2">ρ</th><th rowspan="2">k</th><th rowspan="2">C<sub>p</sub></th><th colspan="3">α = α<sub>0</sub>f<sup>b</sup></th></tr><tr><th>α<sub>0</sub></th><th>b</th><th>w<sub>b</sub></th></tr></thead><tbody><tr><td>Water</td><td>1000</td><td>0.6</td><td>4178</td><td>0.00217</td><td>2</td><td>/</td></tr><tr><td>Skin</td><td>1100</td><td>0.35</td><td>3437</td><td>0.5*</td><td>1.83*</td><td>0.0013*</td></tr><tr><td>Skull</td><td>1990</td><td>0.39</td><td>1300</td><td>8.7</td><td>1.7</td><td>/</td></tr><tr><td>CSF</td><td>1007</td><td>0.6</td><td>4191</td><td>0.00217</td><td>2</td><td>/</td></tr><tr><td>Gray matter</td><td>1039</td><td>0.565</td><td>3680</td><td>0.1**</td><td>1.21**</td><td>0.019*</td></tr><tr><td>White matter</td><td>1043</td><td>0.503</td><td>3600</td><td>0.56**</td><td>1.27**</td><td>0.019*</td></tr><tr><td>Blood</td><td>1060</td><td>0.53</td><td>3840</td><td>0.06</td><td>1.56</td><td>0.019*</td></tr><tr><td>Tumor</td><td>1043</td><td>0.5</td><td>3621</td><td>0.1</td><td>1.21</td><td>0.0038*</td></tr></tbody></table>

*: Mouse **: Cow

<caption>Table 4. Thermal and acoustic parameters of the media in the phantom [58–66].</caption>
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 $\beta$, the speed of sound $v_s$ and the specific heat capacity $C_p$ for each tissue of the model. Expressions of $v_s$ and $\beta$ as a function of the temperature for most of the phantom's tissues were collected from the literature [67–72]. The obtained temperature depending curves for $v_s$ and $\beta$ are plotted on Figs. 3 and 4. The constant values for $\beta$ 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 $\beta$ 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–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 $\alpha_0$ and power $b$ were found [58, 60–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 $\alpha_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.

<table><caption>Table 5. Constant thermal expansion coefficient used in the phantom [58, 79].</caption><thead><tr><th></th><th>Skin [79]</th><th>Skull [58]</th><th>Gray matter [58]</th><th>White matter [58]</th><th>Tumor</th></tr></thead><tbody><tr><td>β (10<sup>-4</sup>·K<sup>-1</sup>)</td><td>3</td><td>1.97</td><td>1.6</td><td>1.6</td><td>1.6</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>Fig. 3. Speed of sound in each tissue of the phantom. Dependence with temperature was found for every tissue except for the skin. A constant value of 1540 m.s⁻¹ was used instead, indicated by the magenta line [61, 67, 68, 70–72].</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Thermal expansion coefficients used in the phantom as a function of the temperature [77, 78].</figcaption></figure>

# 4. Simulation's results

## 4.1. 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.

<figure><img src="image_7.png" /><figcaption>Fig. 5. Heat deposition rate $Q_H$ in W.cm⁻³ in the four cases considered at 808 nm.</figcaption></figure>
4.2. *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.

<figure><img src="image_8.png" /><figcaption>Fig. 6. Temperature evolution at the surface of the tumor with and without added NP. The black cross (left) indicates the position where temperature is calculated, with and without tumor.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Damaged tissues in the phantom at different times of the PTT.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 8. Percentage of tumor burnt (left) and percentage of surrounding healthy tissues burnt (right) in the three tumorous cases out of the four considered cases. Percentage of healthy tissue burnt is computed over the total volume of the brain, which is approximately 395 mm³.</figcaption></figure>

### 4.3. 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><img src="image_11.png" /><figcaption>Fig. 9. Initial pressure maps $p_0 = \Gamma Q_{PA}$ at 37°C (in Pa) of the tumor and its surrounding for the four considered cases.</figcaption></figure>

Fig. 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⁻¹). 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 bandwidth ratio 80%. We performed this operation to give an insight on what a realistic signal would look like.

<figure><img src="image_12.png" /><figcaption>Fig. 10. Top left: PA signals calculated in all different cases at 37°C (temperature at the initial time of the therapy). Top right: PA signal in all different cases at 37°C after applying a gaussian filter of central frequency 2.1 MHz and bandwidth ratio of 80%. Bottom left: PA signals evolution as the PTT progresses in the $\mu_{aT} \times 5$ case. The legend presents the duration of the PTT at the time of acquisition of the corresponding signal. Bottom Right: Zoom in the third peaks of the bottom left plot, corresponding to the tumor's response.</figcaption></figure>

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 "$\mu_{aT} \times 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 $\Delta 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.

<figure><img src="image_13.png" /><figcaption>Fig. 11. Differential signal $\Delta 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 "$\mu_{aT} \times 5$" case.</figcaption></figure>
We computed the ratios $\Delta p/p_i$ of the tumor's response for the three tumorous cases during the PTT in Fig. 12.

<figure><img src="image_14.png" /><figcaption>Fig. 12. Evolution of the PA response of the tumor $\frac{\Delta p}{p_i}$ as a function of the heating time (left) and the temperature rise at the tumor surface (right). Dashed lines in the right figure represent a linear fit for each case. Peak to peak values were used.</figcaption></figure>

Linking the signal evolution with the temperature rise can be done with the Eq.(14). The value of the calibration factor $\gamma$ was calculated in each studied case then summarized in Table 6. It represents the slope of the linear fit in Fig. 12 (right).

<table><caption>Table 6. Values of the constant $\gamma$ calculated from the variations of the received PA signals.</caption><thead><tr><th></th><th>No NP</th><th>μ<sub>aT</sub> × 2</th><th>μ<sub>aT</sub> × 5</th></tr></thead><tbody><tr><td>mean γ (10<sup>-4</sup> °C<sup>-1</sup>)</td><td>-5.8</td><td>4.1</td><td>6.6</td></tr></tbody></table>

It is worth noting, that the value of $\gamma$ 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 $\gamma$ values in the table. Results of these 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.

<figure><img src="image_15.png" /><figcaption>Fig. 13. Temperature evolution calculated from PA signals variations using $\gamma$ compared to simulated evolution for the blood (center) and the tumor's surface (right). Positions of the temperature monitoring is indicated on the left figure. Dashed lines correspond to simulated evolution of $\Delta T$ in the tumor, and solid lines correspond to the approximated evolution obtained with $\gamma$.</figcaption></figure>
# 5. 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.

## 5.1. 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 "$\mu_{aT} \times 5$" case.

## 5.2. 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 "$\mu_{aT} \times 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 "$\mu_{aT} \times 5$", the duration of heating is 139 s. The corresponding volumes of burnt healthy tissues in the "No NP" and "$\mu_{aT} \times 5$" cases are respectively 7.27 mm³ and 6.76 mm³. 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³ of burnt healthy tissues for the "$\mu_{aT} \times 5$" case and 2 mm³ 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].

### 5.3. 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 "$\mu_{aT} \times 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 "$\mu_{aT} \times 5$" case in the Fig. 10 but is also clearly present in the "No NP" and "$\mu_{aT} \times 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 "$\mu_{aT} \times 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 "$\mu_{aT} \times 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 $\Delta p/p_i$, presented in Fig. 12. In
the "$\mu_{aT} \times 2$" and "$\mu_{aT} \times 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 $\Delta 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 "$\mu_{aT} \times 2$" and "$\mu_{aT} \times 5$" cases but not in the "No NP" case. This fact impacts directly the efficiency of the temperature monitoring in the tumor: as $\gamma$ 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 "$\mu_{aT} \times 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 "$\mu_{aT} \times 5$". In the "$\mu_{aT} \times 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 $\Delta 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.

## 5.4. 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 \text{ cm}^{-1}$ at 808 nm and a five-fold increase reach $23 \text{ 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.
## 6. 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} \times 2$" and "$μ_{aT} \times 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.

## Acknowledgement

This project is developed in the field of the "Glioblastoma Risk Attrition by Vectorized Ther-motherapy - GRAVITY" project (PC201613) funded by the ITMO "Plan Cancer 2014–2019" INSERM program. Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high performance computing resources. Authors would like to thank Benoit Larrat-CEA Neurospin and Serge Mensah-LMA, for fruitful discussions.

## Disclosure

The authors declare that there are no conflicts of interest related to this article.

## Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

## References

1. W. He, K. Ai, C. Jiang, Y. Li, X. Song, and L. Lu, "Plasmonic titanium nitride nanoparticles for in vivo photoacoustic tomography imaging and photothermal cancer therapy," *Biomaterials* **132**, 37–47 (2017).

2. A. A. Popov, G. Tselikov, N. Dumas, C. Berard, K. Metwally, N. Jones, A. Al-Kattan, B. Larrat, D. Braguer, S. Mensah, A. Da Silva, M.-A. Estève, and A. V. Kabashin, "Laser-synthesized TiN nanoparticles as promising plasmonic alternative for biomedical applications," *Sci. Reports* **9**, 1194 (2019).

3. N. Harris, M. J. Ford, P. Mulvaney, and M. B. Cortie, "Tunable infrared absorption by metal nanoparticles: The case for gold rods and shells," *Gold Bull.* **41**, 5–14 (2008).

4. P. Wang and G. Zhou, "Frost-heaving pressure in geotechnical engineering materials during freezing process," *Int. J. Min. Sci. Technol.* **28**, 287–296 (2018).
5. M. R. L. Zanetti, A. M. Percebom, T. A. Ribeiro, M. L. Dias, A. P. Oliveira, E. R. Júnior, A. L. Rossi, C. A. G. Soares, and P. H. S. Picciani, "Improving in vitro biocompatibility of gold nanorods with thiol-terminated triblock copolymer," Colloid Polym. Sci. **297**, 1477-1487 (2019).

6. J. F. Hainfeld, M. J. O'Connor, P. Lin, L. Qian, D. N. Slatkin, and H. M. Smilowitz, "Infrared-Transparent Gold Nanoparticles Converted by Tumors to Infrared Absorbers Cure Tumors in Mice by Photothermal Therapy," PLoS One **9**, e88414 (2014).

7. H. Bahadar, F. Maqbool, K. Niaz, and M. Abdollahi, "Toxicity of Nanoparticles and an Overview of Current Experimental Models," Iran. Biomed. J. **20**, 1-11 (2016).

8. C. Bastiancich, A. Da Silva, and M.-A. Estève, "Photothermal Therapy for the Treatment of Glioblastoma: Potential and Preclinical Challenges," Front. Oncol. **10**, 610356 (2021).

9. L. Lavaud, S. Pascal, K. Metwally, D. Gasteau, A. Da Silva, Z. Chen, M. Elhabiri, G. Canard, D. Jacquemin, and O. Siri, "Azacalixphyrins as NIR photoacoustic contrast agents," Chem. Commun. **54**, 12365-12368 (2018).

10. H. Sung, J. Ferlay, R. L. Siegel, M. Laversanne, I. Soerjomataram, A. Jemal, and F. Bray, "Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries," CA: Cancer J. Clin. **71**, 209-249 (2021).

11. A. B. Mariotto, K. Robin Yabroff, Y. Shao, E. J. Feuer, and M. L. Brown, "Projections of the Cost of Cancer Care in the United States: 2010-2020," J. Natl. Cancer Inst. **103**, 117-128 (2011).

12. R. Stupp, M. E. Hegi, W. P. Mason, M. J. van den Bent, M. J. Taphoorn, R. C. Janzer, S. K. Ludwin, A. Allgeier, B. Fisher, K. Belanger, P. Hau, A. A. Brandes, J. Gijtenbeek, C. Marosi, C. J. Vecht, K. Mokhtari, P. Wesseling, S. Villa, E. Eisenhauer, T. Gorlia, M. Weller, D. Lacombe, J. G. Cairncross, and R.-O. Mirimanoff, "Effects of radiotherapy with concomitant and adjuvant temozolomide versus radiotherapy alone on survival in glioblastoma in a randomised phase III study: 5-year analysis of the EORTC-NCIC trial," Lancet Oncol. **10**, 459-466 (2009).

13. Y. Chen and L. Liu, "Modern methods for delivery of drugs across the blood-brain barrier," Adv. Drug Deliv. Rev. **64**, 640-665 (2012).

14. A. H. Mesiwala, L. Farrell, H. Wenzel, D. L. Silbergeld, L. A. Crum, H. Winn, and P. D. Mourad, "High-intensity focused ultrasound selectively disrupts the blood-brain barrier in vivo," Ultrasound Med. Biol. **28**, 389-400 (2002).

15. M. A. Lewis, R. M. Staruch, and R. Chopra, "Thermometry and ablation monitoring with ultrasound," Int. J. Hyperth. **31**, 163-181 (2015).

16. N. W. Lutz and M. Bernard, "Contactless Thermometry by MRI and MRS: Advanced Methods for Thermotherapy and Biomaterials," iScience **23**, 101561 (2020).

17. J. Kosir, D. Vella, M. Lukac, and M. Jezersek, "Towards personalized and versatile monitoring of temperature fields within heterogeneous tissues during laser therapies," Biomed. Opt. Express **12**, 4530 (2021).

18. L. Yin, M. S. R. Gudur, Y.-S. Hsiao, R. E. Kumon, C. X. Deng, and H. Jiang, "Tomographic Reconstruction of Tissue Properties and Temperature Increase for High-Intensity Focused Ultrasound Applications," Ultrasound Med. Biol. **39**, 1760-1770 (2013).

19. J. Shah, S. Park, S. Aglyamov, T. Larson, L. Ma, K. Sokolov, K. Johnston, T. Milner, and S. Y. Emelianov, "Photoacoustic imaging and temperature measurement for photothermal cancer therapy," J. Biomed. Opt. **13**, 034024 (2008).

20. I. V. Larina, K. V. Larin, and R. O. Esenaliev, "Real-time optoacoustic monitoring of temperature in tissues," J. Phys. D Appl. Phys. **38**, 2633-2639 (2005).

21. S.-H. Wang, C.-W. Wei, S.-H. Jee, and P.-C. Li, "Photoacoustic temperature measurements for monitoring of thermal therapy," in *Photons Plus Ultrasound: Imaging and Sensing 2009*, vol. 7177 A. A. Oraevsky and L. V. Wang, eds., International Society for Optics and Photonics (SPIE, 2009), pp. 451-461.

22. K. Metwally, C. Bastiancich, F. Correard, A. Novell, S. Fernandez, B. Guillet, B. Larrat, S. Mensah, M.-A. Estève, and A. Da Silva, "Development of a multi-functional preclinical device for the treatment of glioblastoma," Biomed. Opt. Express **12**, 2264 (2021).

23. V. Cuplov, F. Pain, and S. Jan, "Simulation of nanoparticle-mediated near-infrared thermal therapy using GATE," Biomed. Opt. Express **8**, 1665 (2017).

24. J. C. G. Jeynes, F. Wordingham, L. J. Moran, A. Curnow, and T. J. Harries, "Monte Carlo Simulations of Heat Deposition During Photothermal Skin Cancer Therapy Using Nanoparticles," Biomolecules **9**, 343 (2019).

25. N. Manuchehrabadi and L. Zhu, "Development of a computational simulation tool to design a protocol for treating prostate tumours using transurethral laser photothermal therapy," Int. J. Hyperth. **30**, 349-361 (2014).

26. Y. Ma, P. R. Hof, S. C. Grant, S. J. Blackband, R. Bennett, L. Slatest, M. D. McGuigan, and H. Benveniste, "A three-dimensional digital atlas database of the adult C57BL/6J mouse brain by magnetic resonance microscopy," Neuroscience **135**, 1203-1215 (2005).

27. E. Chan, N. Kovacevic, S. K. Y. Ho, R. M. Henkelman, and J. T. Henderson, "Development of a high resolution three-dimensional surgical atlas of the murine head for strains 129S1/SvImJ and C57Bl/6J using magnetic resonance imaging and micro-computed tomography," Neuroscience **144**, 604-615 (2007).

28. A. Dorr, J. Sled, and N. Kabani, "Three-dimensional cerebral vasculature of the CBA mouse brain: A magnetic resonance imaging and micro computed tomography study," NeuroImage **35**, 1409-1423 (2007).

29. Y. Liu, S. L. Jacques, M. Azimipour, J. D. Rogers, R. Pashaie, and K. W. Eliceiri, "OptogenSIM: a 3D Monte Carlo simulation platform for light delivery design in optogenetics," Biomed Opt Express **6**, 4859-4870 (2015).

30. The MathWorks, Inc., "Matlab image processing toolbox," https://www.mathworks.com/products/
image.html.

31. C. Jean-Quartier, F. Jeanquartier, I. Jurisica, and A. Holzinger, "In silico cancer research towards 3R," BMC Cancer **18**, 408 (2018).

32. H. H. Pennes, "Analysis of Tissue and Arterial Blood Temperatures in the Resting Human Forearm," J. Appl. Physiol. **1**, 93-122 (1948).

33. S. Chandrasekhar, *Radiative transfer* (Oxford University, 1950).

34. S. A. Sapareto and W. C. Dewey, "Thermal dose determination in cancer therapy," Int. J. Radiat. Oncol. **10**, 787-800 (1984).

35. M. W. Dewhirst, B. L. Viglianti, M. Lora-Michiels, M. Hanson, and P. J. Hoopes, "Basic principles of thermal dosimetry and thermal thresholds for tissue damage from hyperthermia," Int. J. Hyperth. **19**, 267-294 (2003).

36. G. C. van Rhoon, T. Samaras, P. S. Yarmolenko, M. W. Dewhirst, E. Neufeld, and N. Kuster, "CEM43°C thermal dose thresholds: a potential guide for magnetic resonance radiofrequency exposure levels?" Eur. Radiol. **23**, 2215-2227 (2013).

37. B. E. Treeby and B. T. Cox, "Modeling power law absorption and dispersion for acoustic propagation using the fractional Laplacian," J. Acoust. Soc. Am. **127**, 2741-2748 (2010).

38. P. Subochev, E. Smolina, E. Sergeeva, M. Kirillin, A. Orlova, D. Kurakina, D. Emyanov, and D. Razansky, "Toward whole-brain in vivo optoacoustic angiography of rodents: modeling and experimental observations," Biomed. Opt. Express **11**, 1477 (2020).

39. T. K. Bhattacharyya and J. R. Thomas, "Histomorphologic Changes in Aging Skin: Observations in the CBA Mouse Model," Arch Facial Plast Surg **6**, 21 (2004).

40. L. Lemaire, F. Franconi, B. Siegler, C. Legendre, and E. Garcion, "In vitro expansion of U87-MG human glioblastoma cells under hypoxic conditions affects glucose metabolism and subsequent in vivo growth," Tumor Biol. **36**, 7699-7710 (2015).

41. A. Clavreul, M. Pourbaghi-Masouleh, E. Roger, N. Lautram, C. N. Montero-Menei, and P. Menei, "Human mesenchymal stromal cells as cellular drug-delivery vectors for glioblastoma therapy: a good deal?" J. Exp. Clin. Cancer Res. **36**, 135 (2017).

42. A. Clavreul, E. Roger, M. Pourbaghi-Masouleh, L. Lemaire, C. Tétaud, and P. Menei, "Development and characterization of sorafenib-loaded lipid nanocapsules for the treatment of glioblastoma," Drug Deliv. **25**, 1756-1765 (2018).

43. A. Capart, K. Metwally, C. Bastiancich, and A. D. Silva, "Development of a realistic phantom of a mouse's head: application to photoacoustic temperature monitoring during photothermal therapy," in *European Conferences on Biomedical Optics 2021 (ECBO)*, (Optical Society of America, 2021), p. ETu4D.1.

44. O. Akhavan and E. Ghaderi, "Graphene Nanomesh Promises Extremely Efficient In Vivo Photothermal Therapy," Small **9**, 3593-3601 (2013).

45. W. Lu, M. P. Melancon, C. Xiong, Q. Huang, A. Elliott, S. Song, R. Zhang, L. G. Flores, J. G. Gelovani, L. V. Wang, G. Ku, R. J. Stafford, and C. Li, "Effects of Photoacoustic Imaging and Photothermal Ablation Therapy Mediated by Targeted Hollow Gold Nanospheres in an Orthotopic Mouse Xenograft Model of Glioma," Cancer Res. **71**, 6116-6121 (2011).

46. E. S. Day, P. A. Thompson, L. Zhang, N. A. Lewinski, N. Ahmed, R. A. Drezek, S. M. Blaney, and J. L. West, "Nanoshell-mediated photothermal therapy improves survival in a murine glioma model," J. Neuro-Oncology **104**, 55-63 (2011).

47. C. Li and L. V. Wang, "Photoacoustic tomography and sensing in biomedicine," Phys. Med. Biol. **54**, R59-R97 (2009).

48. D. Marti, R. N. Aasbjerg, P. E. Andersen, and A. K. Hansen, "MCmatlab: an open-source, user-friendly, MATLAB-integrated three-dimensional Monte Carlo light transport solver with heat diffusion and tissue damage," J. Biomed. Opt. **23**, 1 (2018).

49. S. L. Jacques, "Coupling 3D Monte Carlo light transport in optically heterogeneous tissues to photoacoustic signal generation," Photoacoustics **2**, 137-142 (2014).

50. B. E. Treeby and B. T. Cox, "k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields," J. Biomed. Opt. **15**, 021314 (2010).

51. H. Buiteveld, J. Hakvoort, and M. Donze, "Optical properties of pure water," Proc. SPIE 2258, Ocean. Opt. **XII** (1994).

52. Caetano Padial Sabino, A.M. Deana, Tania Mateus Yoshimura, and Martha S Ribeiro, "The optical properties of mouse skin in the visible and near infrared," J. Photochem. Photobiol. B: Biol. **160**, 72-78 (2016).

53. V. Tuchin, *Tissue Optics* (SPIE, 1000 20th Street, Bellingham, WA 98227-0010 USA, 2007).

54. L. Wu, Y. Lin, and T. Li, "Effect of Human Brain Edema on Light Propagation: A Monte Carlo Modeling Based on the Visible Chinese Human Dataset," IEEE Photonics J. **9**, 1-10 (2017).

55. A. N. Yaroslavsky, P. C. Schulze, I. V. Yaroslavsky, R. Schober, F. Ulrich, and H.-J. Schwarzmaier, "Optical properties of selected native and coagulated human brain tissues in vitro in the visible and near infrared spectral range," Phys. Medicine Biol. **47**, 2059-2073 (2002).

56. N. Bosschaart, G. J. Edelman, M. C. G. Aalders, T. G. v. Leeuwen, and D. J. Faber, "A literature review and novel theoretical approach on the optical properties of whole blood," Lasers Med. Sci. **29**, 453 (2014).

57. E. A. Genina, A. N. Bashkatov, D. K. Tuchina, P. A. Dyachenko (Timoshina), N. Navolokin, A. Shirokov,
A. Khorovodov, A. Terskov, M. Klimova, A. Mamedova, I. Blokhina, I. Agranovich, E. Zinchenko, O. V. Semyachkina-Glushkovskaya, and V. V. Tuchin, "Optical properties of brain tissues at the different stages of glioma development in rats: pilot study," Biomed. Opt. Express **10**, 5182 (2019).

58. Francis A. Duck, *Physical Properties of Tissues: A comprehensive reference book* (Elsevier, 1990).

59. S. M. Yacoob and N. S. Hassan, "FDTD analysis of a noninvasive hyperthermia system for brain tumors," Biomed. Eng. Online **11**, 47 (2012).

60. J. M. M. Pinkerton, "The Absorption of Ultrasonic Waves in Liquids and its Relation to Molecular Constitution," Proc. Phys. Soc. B **62**, 129-141 (1949).

61. P. Bhagat, W. Kerrick, and R. Ware, "Ultrasonic characterization of aging in skin tissue," Ultrasound Med. Biol. **6**, 369-375 (1980).

62. V. R. Amin, "Ultrasonic attenuation estimation for tissue characterization," Master of Science, Iowa State University (1989).

63. B. E. Treeby, E. Z. Zhang, A. S. Thomas, and B. T. Cox, "Measurement of the Ultrasound Attenuation and Dispersion in Whole Human Blood and its Components From 0-70 MHz," Ultrasound Med. Biol. **37**, 289-300 (2011).

64. E. R. Muir, Q. Shen, and T. Q. Duong, "Cerebral blood flow MRI in mice using the cardiac-spin-labeling technique," Magn. Reson. Med. **60**, 744-748 (2008).

65. W. T. Stott, M. D. Dryzga, and J. C. Ramsey, "Blood-flow distribution in the mouse," J. Appl. Toxicol. **3**, 310-312 (1983).

66. S.-G. Kim and J. J. H. Ackerman, "Quantitative determination of tumor blood flow and perfusion via deuterium nuclear magnetic resonance spectroscopy in mice," Cancer Res. **48**, 3449-3453 (1988).

67. J. Bamber and C. Hill, "Acoustic properties of normal and cancerous human liver—I. Dependence on pathological condition," Ultrasound Med. Biol. **7**, 121-133 (1981).

68. R. M. Souza, R. P. B. Costa-Felix, and A. V. Alvarenga, "The influence of temperature on the speed of sound of cortical bone phantom: a metrological view," J. Physics: Conf. Ser. **975**, 012021 (2018).

69. G. Clement and K. Hynynen, "Ultrasound phase-contrast transmission imaging of localized thermal variation and the identification of fat/tissue boundaries," Phys. medicine biology **50**, 1585-600 (2005).

70. C. H. Seo, Y. Shi, S.-W. Huang, K. Kim, and M. O'Donnell, "Thermal strain imaging: a review," Interface Focus. **1**, 649-664 (2011).

71. C. R. Hill, J. C. Bamber, and G. R. ter Haar, eds., *Physical Principles of Medical Ultrasonics* (Wiley, 2004), 1st ed.

72. W. D. Wilson, "Speed of Sound in Distilled Water as a Function of Temperature and Pressure," J. Acoust. Soc. Am. **31**, 1067-1072 (1959).

73. S. M. Nikitin, T. D. Khokhlova, and I. M. Pelivanov, "Temperature dependence of the optoacoustic transformation efficiency in ex vivo tissues for application in monitoring thermal therapies," J. Biomed. Opt. **17**, 061214 (2012).

74. C. W. Song, "Effect of local hyperthermia on blood flow and microenvironment: A review," Cancer Res. **44**, 4721s-4730s (1984).

75. H. Masuda, A. Hirata, H. Kawai, K. Wake, S. Watanabe, T. Arima, F. Pouletier de Gannes, I. Lagroye, and B. Veyret, "Local exposure of the rat cortex to radiofrequency electromagnetic fields increases local cerebral blood flow along with temperature," J. Appl. Physiol. **110**, 142-148 (2011).

76. E. Moriyama, "Cerebral Blood Flow Changes During Localized Hyperthermia," Neurol. Med.-Chir. **30**, 923-929 (1990).

77. G. S. Kell, "Precise representation of volume properties of water at one atmosphere," J. Chem. Eng. Data **12**, 66-69 (1967).

78. H. Hinghofer-Szalkay, H. Leopold, T. Kenner, and H. Holzer, "Über den thermischen Ausdehnungskoeffizienten des Blutes und seiner Komponenten," Biomedizinische Tech. **25**, 151-157 (1980).

79. C. H. Li, S. A. Li, Z. H. Huang, and W. B. Xu, "Skin Thermal Effect by FE Simulation and Experiment of Laser Ultrasonics," Appl. Mech. Mater. **24-25**, 281-286 (2010).

80. F. J. van der Meer, D. J. Faber, M. C. G. Aalders, A. A. Poot, I. Vermes, and T. G. van Leeuwen, "Apoptosis- and necrosis-induced changes in light attenuation measured by optical coherence tomography," Lasers Med. Sci. **25**, 259-267 (2010).

81. Q. Chen, Y. Ren, Y. Yin, and H. Qi, "Anisotropic scattering characteristics of nanoparticles in different morphologies: improving the temperature uniformity of tumors during thermal therapy using forward scattering," Biomed. Opt. Express **12**, 893 (2021).