Numerical Study of Photoacoustic Pressure for Cancer Therapy

Abstract

A commonly used therapy for cancer is based on the necrosis of cells induced by heating through the illumination of nanoparticles embedded in cells. Recently, the photoacoustic pressure shock induced by the illumination pulse was proved and this points to another means of cell destruction. The purpose of this study is to propose a model of the photoacoustic pressure in cells. The numerical resolution of the problem requires the accurate computation of the electromagnetism, the temperature and the pressure around the nanostructures embedded in a cell. Here, the problem of the interaction between an electromagnetic excitation and a gold nanoparticle embedded in a cell domain is solved. The variations of the thermal and photoacoustic pressures are studied in order to analyze the pressure shock wave inducing the collapse of the cell’s membrane in cancer therapy.

1. Introduction

In the last decades, the theory and the applications of photoacoustics have been developed, in particular for the photothermal measurement techniques [1,2]. The applications of photoacoustics include spectroscopy, monitoring deexcitation processes, probing physical properties of materials, tomography and applications in photothermal therapy [3,4,5]. Photothermal therapy and photoacoustic techniques are procedures using laser illumination of photoabsorbers. Such therapies are non-invasive when working in the near-infrared spectrum. Metallic nanoparticles embedded in the tumor site absorb the laser irradiation and transfer the thermal energy to the surrounding medium (i.e., cancer cells) [6,7]. Photoacoustics refers to the generation of acoustic waves by modulating an optical radiation. Photoacoustic imaging is applied to measure temperature by tracking the temperature-induced changes in the photoacoustic signal amplitude [8].

In this context, we propose a numerical model for the accurate computation of photoacoustic pressure around gold nanopatricle embedded in a cell and the analysis of the spatio-temporal evolution of the pressure. The numerical method, including a 3D version of the adaptive remeshing scheme with error estimator [9,10], takes into account the local and global variations of the pressure ensuring the convergence of the solution to the physical solution. The organization of the paper is as follows: Section 2 describes the equations of the problem and the numerical method. In Section 3, the results of the simulations are presented and a discussion on the ability of the pressure to destroy a cell is acheived before concluding the paper.

2. Photoacoustic Model

This section presents the physical model and the numerical method used to solve the photoacoustic problem. The photoacoustic or thermoacoustic process refers to the generation of acoustic waves by absorption of electromagnetic energy. Such an absorption of the radiation by a nanoparticle acts like a heat source that depends on the space $\mathbf{x}$ and time t and is given by:

$$Q(\mathbf{x},t)=\frac{\omega }{2}{ϵ}_{0}Im\left({ϵ}_r\left(\mathbf{x}\right)\right){\left|\mathbf{E}(\mathbf{x},t)\right|}^2,\forall \mathbf{x}\in \mathsf{\Omega },$$

(1)

where ω is the angular frequency of the incoming wave, $\mathbf{E}$ the electric field, ${ϵ}_0$ is the permittivity of vacuum, $Im\left({ϵ}_r\right)$ is the imaginary part of ${ϵ}_r$ the relative complex permittivity of materials that are functions of the 3D spatial coordinates $\mathbf{x}$ in the domaine Ω. The laser illumination pulse satisfies the electromagnetic equation:

$$\nabla \times \left(\frac{1}{{ϵ}_r\left(\mathbf{x}\right)}\nabla \times \mathbf{E}(\mathbf{x},t)\right)-\frac{1}{c^2}\frac{{\partial }^2\mathbf{E}(\mathbf{x},t)}{\partial t^2}=\mathbf{0},\forall \mathbf{x}\in \mathsf{\Omega },$$

(2)

where ($\nabla \times$) is the curl operator and c the velocity of light. Equation (2) is solved with a radiation boundary condition [11,12] at the external border Γ of the computationnal domain Ω corresponding to the free propagation of the scattered field and the incident laser pulse of amplitude:

$$E_0={\left[\frac{2c{\mu }_0{E}_{laser}}{\pi R_{beam}^2{\tau }_{pulse}}\right]}^{1/2},$$

(3)

where ${\mu }_0$ is the permeability of the vacuum, $E_{laser}$ is the laser energy, $R_{beam}$ is the radius of the laser beam and ${\tau }_{pulse}$ the duration of the laser pulse. From such a heat source, the temperature elevation is governed by the equation:

$$\left[\rho \left(\mathbf{x}\right)C_p\left(\mathbf{x}\right)\frac{\partial }{\partial t}-\nabla .\left(\kappa \left(\mathbf{x}\right)\nabla \right)\right]T(\mathbf{x},t)=Q(\mathbf{x},t),\forall \mathbf{x}\in \mathsf{\Omega },$$

(4)

where $(\nabla .)$ and ∇ are the divengence and the gradient operators, ρ is the mass density of the material, $C_p$ is the specific heat capacity, κ the thermal conductivity and T is the temperature. These are functions of time and of the 3D spatial coordinates $\mathbf{x}$. Therefore, the resolution of the thermal problem requires solving the heat equation (Equation (4)) with a source Q produced by the incoming electromagnetic illumination (Equation (1)) with the boundary condition $T(\mathbf{x},t)=T_0$ a the boundary of the computational domain (i.e., $\forall \mathbf{x}\in \mathsf{\Gamma }$). Such a thermal expansion induces acoustic displacement and the equation for photoacoustic generation is derived [2]:

$$\left[\frac{1}{\rho \left(\mathbf{x}\right)v^2\left(\mathbf{x}\right)}\frac{{\partial }^2}{\partial t^2}-\nabla .\left(\frac{1}{\rho \left(\mathbf{x}\right)}\nabla \right)\right]P(\mathbf{x},t)=\beta \left(\mathbf{x}\right)\frac{{\partial }^2}{\partial t^2}T(\mathbf{x},t),\forall \mathbf{x}\in \mathsf{\Omega },$$

(5)

with v the sound velocity in the medium, β is the expansion coefficient and P is the acoustic pressure with initial and boundary conditions:

$$P(\mathbf{x},t=0)=p_0\mathrm{and}n.\left(\frac{1}{\rho \left(\mathbf{x}\right)}\nabla \right)P(\mathbf{x},t)=0,\forall \mathbf{x}\in \mathsf{\Gamma },$$

(6)

with $p_0$ being the initial pressure and n the outgoing normal vector. The resolution of the weak coupled photothermal problem gives the spatio-temporal distribution of the pressure in the computational domain.

Equations (1)–(5) are numerically solved with a 3D Finite Element Method (FEM) with adaptive remeshing scheme using an a priori error estimator based on the Hessian of the solution calculated at the previous step of remeshing. The Finite Element Method (FEM) is now well known and has been applied for many years in mechanics, thermodynamics, electromagnetics and in electrical engineering [13,14] and consists of solving systems of partial differential equation with boundary conditions in open or closed domains [11]. The general problem is solved on a discrete mesh of the domain [11] and the unknown physical quantities ($\mathbf{E}$, T, P) are computed at the nodes of the mesh by using a variational method. The control of the error on the solution is achieved through an extended and improved 3D method, including a process of iterative remeshing [9,10,12]. The accuracy of the computed solutions are closely related to the quality of the mesh [15,16]. The quality of the solution is improved by adapting the size of the mesh elements to the physical solution [10,12,17] through the remeshing process with adaptive loops. If strong variation is induced in the unknown physical quantities of interest, the convergence of the solution to a stable solution requires such mesh adaptions [9,10]. At each step of the adaption process, the approximations of the solutions of the thermic and acoustic wave equations T and P, the temperature and the pressure are calculated in the basis of the interpolation polynomials between nodes. The maximum deviation between the solution associated with the mesh and the exact solution at nodes is limited by using the interpolation error (which is based on an estimation of the discrete Hessian of the solution) [18,19] and is obtained from OPTIFORM software (adaptive remeshing generating isotropic or anisotropic meshes) [9], to govern the remeshing of the domain (with minimum and maximum element sizes $h_{\min }=0.1$ nm, $h_{\max }=100$ nm, and a tolerance on the relative error on the computed unknown physical quantities $\delta =0.001$). Therefore, the domain is entirely remeshed each time and a new mesh is constructed.

3. Numerical Results and Discussion

Here, we consider a 3D spherical gold nanoparticle of radius $R_p=50$ nm embedded in a 3D spherical cell of radius $R_c=1.0\mathsf{\mu }$m illuminated by an electromagnetic wave at $\lambda =532$ nm and immersed in a water domain Ω of radius $R_w=22\mathsf{\mu }$m at initial temperature $T_0$ = 37.5 °C (i.e., 310.5 K). This wavelength is chosen to be near the plasmon resonance for such a gold nanoparticle where light absorption is maximum and was already used as a contrast agent for in vivo tumor imaging [20]. The illumination pulse of duration ${\tau }_{pulse}=5$ ns is achieved with a laser with energy $E_{laser}=10$ mJ and a radius of the spotsize $R_{beam}=15$ mm.

The parameters of Equations (4)–(5) for the thermic and acoustic problems are given in

Table 1. Values of the parameters for the photothermal model for the three media (water, cell, gold) and $j^2=-1$.

Material	ρ (kg·${\mathbf{m}}^{-\mathbf{3}}$)	${\mathit{C}}_{\mathit{p}}$ (${\mathbf{m}}^2$·${\mathbf{s}}^{-1}$·${\mathbf{K}}^{-1}$)	κ (kg·m·${\mathbf{s}}^{-3}$·${\mathbf{K}}^{-1}$)	v (m·${\mathbf{s}}^{-1}$)	β (${\mathbf{K}}^{-1}$)	${\mathit{ϵ}}_{\mathit{r}}$
water	1000	4185	0.60	1493	0.000215	1.768
cell	1090	2185	1.20	1535	0.000515	1.768
Au	19,300	129	250	3240	0.000014	−5.58 + j2.25

.

Figure 1a–d show the geometries of the problem with the meshes’ associated materials for the computational domain and the cell, respectively. The relative permittivities of materials (which describe the electromagnetic problem of Equations (1) and (2)) are given in Table 1. Here, the materials of the problem are considered to be isotropic and homogeneous.

The results of the computation of the pressure maps are illustrated in Figure 2. Figure 2a–f shows the associated variations of the pressure at 0.5, 1.5, 2.5, 4.5, 14.5 and 40 ns after the illumination. The variations and the propagation of the pressure wave is shown in the hole computational domain. The numerical scheme not only takes into account the shape and size of the nanoparticle but also the variations induced by the electromagnetism, the temperature and the pressure. The pressure is initiated at the beginning of the temperature shock that is induced by the modulation of the illumination. Such an acoustic pressure wave propagates in the whole media (cell and water, corresponding to a spatial extension of some micrometers the cell membrane has spatial extension from −1 $\mathsf{\mu }$m to 1 $\mathsf{\mu }$m) before vanishing. That contrasts with the temperature propagation (see Figure 3a,b) that is mainly located in the near vicinity of the gold nanoparticle (i.e., only a few nanometers from the gold nanoparticle, corresponding to the spatial expansion of the nanoparticle along the y-axis: −50 to 50 nm). Therefore, only the nearest zone around the nanoparticle sees the temperature elevation. Figure 4 illustrates that the cell’s membrane is receiving a pressure shock wave (i.e., the variation of pressure is non vanishing at the position of the cell membrane). Figure 5 shows the maximum of the pressure at the position of the cell membrane (i.e., at $1\mathsf{\mu }$m) at the initiation of the temperature and pressure shock and the evolution of the pressure at the limit of the membrane. For an illumination at $\lambda =532$ nm, the two curves are obtained for two nanoparticle shapes. The first one consists of a spherical gold nanoparticle of radius $R=50$ nm, the second one is an ellipsoidal particle of semi-axis a = 100 nm and b = 50 nm. Such a pressure variation can induce a deformation of the cell. The deformation is determined by the shear rate of the streaming flow field (computed from the strain tensor) and the cell elasticity [21,22]. Here, with the used materials, the shear stress components are: ${\sigma }_{zz}=\Delta P\approx 300$ Pa, ${\sigma }_{xx}={\sigma }_{yy}=\Delta PR/\left(2e\right)\approx 0.05$ MPa with R being the radius of the cell, e being the width of the cell’s membrane and the associated shear modulus being $\tau \approx 0.01$ to $0.1$ GPa [23]. Moreover, the maximum rate of strain is about $S_{\max }\approx {10}^4$–${10}^5{\mathrm{s}}^{-1}$ (i.e., $S_{\max }\approx \delta R/R/\Delta t$ [24] where $\Delta t$ is the evolution time). That leads to a maximum tension in the membrane ${\tau }_{\max }\approx \mu S_{\max }R$ varying in the range $[{10}^{-4};{10}^{-1}]{\mathrm{N}·\mathrm{m}}^{-1}$, where $\mu \approx 0.1$ to $0.001{\mathrm{kg}·\mathrm{m}}^{-1}$·${\mathrm{s}}^{-1}$ is the dynamic viscosity of the aqueous solution in the cell (i.e., depending on the cell composition). Such a tension is sufficient to stabilize and neglect the thermal fluctuation of the membrane area, but will only expand the actual membrane area by a fraction $\Delta S/S\approx 2.0\times {10}^{-3}$ (with an expansion modulus of about $0.24{\mathrm{N}·\mathrm{m}}^{-1}$ [24]). To achieve the shear of the membrane (or at least opening of pores), a ratio of area $\Delta S/S\ge 0.03$ is necessary [24]. To extend the application in cancer therapy, the illumination must use wavelengths in the biological window (800 to 1000 nm). Such an illumination of a spherical gold nanoparticle at $\lambda =800$ nm is shown in Figure 5. In such a window, the imaginary part of the permittivity of the gold nanoparticle is smaller and requires us to design the particle shape or to increase the energy laser [25]. Moreover, the use of larger cell (typically a cell size is varying from 1 to 100 $\mathsf{\mu }$m) or a mix of larger nanoparticles, greater pulse energy or a higher repetition rate of ultra-short illumination shocks (which will induce an increase in the local pressure variation at the membrane boundary and therefore the maximum rate of strain) can open the way to a process inducing a local stress (shear stress) at the membrane in order to produce a disruption of the cell’s membrane. Here, we focus on the study of one nanoparticle, nevertheless, metal nanoparticles can also be taken up by cells in 1–15 nanoparticles/vesicles [26,27]. The increase of the number of nanoparticles, if aggregated, has a similar effect of increasing the size and shape of the nanoparticle. Indeed, each nanoparticle plays the role of local temperature source inducing a photoacoustic pressure wave. Therefore, multiple pressure waves can be induced at the same time and the pressure shock at the membrane can also be increased due to the aggregated nanoparticles. Such a disruption could be produced by cytolysis (or osmotic lysis) that occurs when a cell bursts due to an osmotic imbalance that has caused excess water being released into the cell.

4. Conclusions

The paper focuses on the analysis of the direct 3D numerical computation of the time-dependant photoacoustic pressure of a gold nanoparticle embedded in a cell. The adaption of the computation process was achieved to optimize both the temperature amplitude and the propagation of the acoustic wave pressure. The variation of the pressure at the limit of the cell’s membrane relative to the time of excitation is presented. We discuss the ability of the use of pressure to disrupt the cell membrane, which can be an alternative way to induce lysis of the membrane. The use of complex particles or cell shapes can extend the domain of its application, especially to the study and the design of a pressure-based non invasive theragnostic for cancer.