Gold Nanoparticles as Photothermal Agent in Cancer Therapy: Theoretical Study of Concentration and Agglomeration Effects on Temperature

: One promising cancer therapy is related to the treatment of diseased cells through thermal ablation by an individual or an agglomeration of nanoparticles acting as photothermal agents. The main principle of such a therapy consists in the photo-energy absorption by the nanoparticles and its conversion into heat in order to kill the biological media/cells in the neighboring regions of such a photothermal agent. Nevertheless, such a therapy must preserve the surrounding healthy cells (or biological media). In case of agglomerates of nanoparticles, the local concentrations of nanoparticles may increase the temperature locally. In this paper, we use the ﬁnite element method to calculate the temperature elevation for agglomerations of nanoparticles in a biological medium/cell. The positions of nanoparticles, forming the agglomerates, are randomly generated. The temperature elevation for such agglomerations of nanoparticles is then analyzed. We show that the control of the concentration of nanoparticles can preserve the efﬁciency of the thermal agent, but with limited risk of damage to the surrounding biological media/cells.


Introduction
In the last decades, photothermic research, developments, and applications have been the object of intense activity [1,2], with particular emphasis on the scope of photothermal therapy (PTT). The local heat source produced by nanoabsorbers in the near-infrared region (NIR) has revealed a new and promising non-invasive therapeutic area. The thermal ablation consists in producing local heat, through the illumination of nanoparticles, in order to destroy the tumoral cell [3][4][5][6][7]. Such nanoparticles (or agglomerations of nanoparticles with complex shapes) form photothermal agents (PA). The main specificity of such therapies in NIR regime are to treat very local tumoral cells and preserve the surrounding tissues [8][9][10]. Indeed, NIR light (650-950 nm) induces minimal absorption by skin and tissue, leading to minimal tissue invasion and deeper tissue penetration (up to 10 cm). A laser illumination is absorbed by the PA, producing a local temperature elevation. Such a temperature elevation depends on the material, the shape, and the size of the PA, but also on the distance between the nanoparticles. These local heat sources are expected to satisfy two conditions: first, they will overpass the thermal ablation threshold within a range of a reasonable maximum temperature limit; second, they will avoid exceeding the thermal ablation temperature outside the treated biological media/cells (or the target cancer cells). Such a control of the temperature is essential.
For a specific cancer cell, the only difference is not directly on the temperature elevation as a function of the agglomeration of nanoparticles, but on the material functionalization (or the functionalization of individual nanoparticles) improving the binding to the cancer cell [11][12][13][14]. In vitro, that consists in a transfection of the nanoparticles to cancer cells before the internalization of the nanoparticles in the cancer cells. In vivo, the nanoparticles are injected intravenously or intratumorally before migrating to the tumor or cancer cells [15].
This theoretical study focus on the process of temperature elevation above the ablation threshold (and under a limit of maximum temperature) as function of the nanoparticles concentration in a generic biological medium/cell. Here, we compute and analyze the spatial evolution of the temperature around agglomerations of nanoparticles in a range of parameters (the shape and size of the agglomerates, and the diameters and concentration of the nanoparticles) [16,17]. Therefore, we obtain a range of concentrations of particles, satisfying both the constraints of efficiency of PA and limiting the risk for the surrounding media/cells. The computation of the temperature elevation induced by a single nanoparticle and electromagnetically coupled agglomerates in a biological media or cell are achieved through the finite element method (FEM). The weak-coupled electromagnetic and heat partial differential equations are solved by using a 3D adaptive remeshing scheme with an error estimator controlling the accuracy of the numerical solutions in the multiscale domain of computation [18][19][20][21]. Statistical analysis of the temperatures produced from random positions of nanoparticles reveals the domain of an adequate concentration of nanoparticles.
The organization of the paper is as follows. In Section 2, we present the main equations, the physics parameters, and the numerical method (details are given in Appendices A.1 and A.2). In Section 3, we show the simulation results and discuss the effects of the mass concentration and of the agglomerate configurations and, therefore, on the temperature elevation before concluding.

Models for Photothermic Problem and Agglomerates
Here, we describe the model used to solve the PTT problem in a volumic domain of computation Ω with a surface boundary Γ. We also give the model used to describe agglomerates.

Photothermic Model
When illuminated by electromagnetic waves, absorbing materials are sources of temperature elevation. The gold nanoparticles absorb the light energy and play the role of a heat source. The heat density of such heat sources depends on the local electromagnetic field E(x) as follows: where ω is the incident angular frequency, 0 is the permittivity of a vacuum, and r is the complex relative permittivity of materials. The electromagnetic field E(x) is the solution of the Helmholz equation with a radiation boundary condition [20,22] (see Appendix A.1 for electromagnetic and thermic equations). With such a heat source, the temperature elevation is governed through the heat equation (Equation (2)): with the time detoned by t, the temperature by T, the volumic mass density by ρ, the specific heat capacity by C p , the thermal conductivity by κ, and (∇·) and ∇ denoting the divergence and gradient operators, respectively. The solution of the full problem requires solving the heat equation (Equation (2)) by including the electromagnetic source Q (Equation (1)). The distribution of the temperature is obtained by solving the weak-coupled photothermic problem with an illumination duration t longer than the thermic characteristic time of materials τ thermic (i.e., τ thermic ≈ 1 ms). The solution of Equation (2) converges to the stationary solution (i.e., t ≥ τ thermic , then ∂T(x, t)/∂t −→ 0, corresponding to a continuous illumination). The characteristic electromagnetic time τ electro (e.g., 1.9 ≤ τ electro ≤ 2.9 fs in the NIR domain of 650-950 nm) is much shorter than the characteristic thermic time τ thermic . Therefore, the electromagnetic and thermic problems induce a weak coupling and could be solved sequentially.
Equations (1) and (2) are solved through a 3D FEM with an adaptive remeshing process, ensuring the convergence of the numerical solution to a stable solution [18][19][20]. The remeshing process (see Appendix A.2) ensures the automatic and systematic control of the accuracy of computation in the multiscale domain (e.g., the typical size of the nanoparticles is small, between 10-100 nanometers, and that of the cell is 1-10 microns). Figure 1a shows a schematic of the referenced problem: agglomerates of individual gold nanoparticles (blue circles) on individual volume V np and permittivity np , embedded in a volume of biological media V bm with relative permittivity bm . The computations are achieved in a domain Ω with boundary Γ. Figure 1b shows an example of the mesh domain associated with such an agglomerate of nanoparticles.  With such a numerical model, the temperature variations that are induced by the electromagnetic and thermic coupling in the agglomerate of nanoparticles in the biological media/cells are studied.

Model for Agglomerates of Individual Nanoparticles
For a given concentration of metallic nanoparticles, the process of forming agglomerates of nanoparticles is complex [23][24][25]. Such metallic nanoparticles can agglomerate by 1-15 nanoparticles/vesicles [26,27]. In this paper, we consider individual nanoparticles (each one, depending on its size, can produce a sub-or upper-ablation threshold temperature elevation). The relative volumic concentration C rv and the relative mass concentration C rm can be defined and related as: where ρ np , ρ bm , V = NV np , and V bm are the mass densities and global volumes of the gold nanoparticles and the biological media (or cell), respectively. Therefore, the number of nanoparticles is N ≈ C rv V bm /V np . For an agglomerate of N identical nanoparticles of diameters D np , the inter-distances between the constitutive nanoparticles are additional parameters. Indeed, the electromagnetic coupling between nanoparticles induces an electromagnetic field enhancement and produces a temperature elevation in the biological media/cells. However, such an electromagnetic coupling occurs at a short distance between individual nanoparticles (d electro ≤ 50-70 nm ≈ D np ). Moreover, a thermic coupling appears if the inter-distance is typically lower than d thermic ≈ 1 µm [21]. Consequently, we can distinguish three spatial distributions of nanoparticles: • a dense agglomerate: characterized by strong coupling for d ≤ d electro due to both electromagnetic and thermic couplings; • a non-dense agglomerate: characterized by soft coupling for d electro ≤ d ≤ d thermic only due to thermic coupling; • a non-agglomerate: characterized by no coupling for d ≥ d thermic ; each nanoparticle is independant (i.e., isolated individual nanoparticle).
With this spatial distribution, the number N of nanoparticles can be written as: with N st i , N so i , and N no i being the numbers of particles participating in the strong, soft and no coupling (associated with the nanoparticles for the dense agglomerate, the non-dense agglomerate and the non-agglomerate) with occurences n st , n so and n no , respectively. These numbers are characteristics for each kind of coupling that induces a specific temperature elevation. With such a model for agglomerates and with the photothermic model, the temperature elevation in a biological medium/cell is calculated through FEM. The maximum temperature elevation is directly related to these kinds of agglomerate and is given by: (5) where the N st max is the higher strong-coupling N st -agglomerate present in the N-agglomerate, and T 1,max is the maximum temperature for an isolated nanoparticle with a diameter of D np . First we will consider single dense agglomerates before analyzing statistics on random mixtures of multiple agglomerates.

Numerical Results and Discussion
The following material parameters are used in the numerical simulations. The volume of the computational domain Ω is 100 µm 3 . The illumination of the volume of the biological media/cells is achieved at λ = 830 nm. Such a wavelength is in NIR, where the laser irradiation can both penetrate soft tissues (1-10 cm deep) with low absorption and overlap the absorption band of gold nanoparticles [5,28]. The illumination beam has a power of P w ≈ 1.0 W and a beam diameter of D beam = 0.15 mm. The body temperature at the boundary of the computational domain Ω is T 0 = 37.0 • C. The electromagnetic and thermal parameters are given in Table 1 [29,30]. All constitutive materials of the nanoparticles are considered isotropic and homogeneous. Table 1. Values of the parameters (the volumic mass density ρ, the specific heat capacity C p , the thermic conductivity κ, and the relative permittivity r ) for the photothermic model and for the media (biological media/cells and gold nanoparticles) with i 2 = −1. The ranges of diameters (D np and D bm ) and the relative concentrations (C rm and C rv ) used in this study for both materials are also given. The agglomerates are constituted by individual gold nanoparticles of varying diameters D np , with relative volumic concentrations C rv (or the relative mass concentration C rm ) in the biological media/cells. The ranges of diameters and concentrations used in this paper are given in Table 1.
In the following, we first analyze the relevant constitutive parameters for the agglomerates of identical nanoparticles before investigating the temperature elevation produced by mixtures of agglomerates.

Single Dense Agglomerate of Identical Nanoparticles
The case of a single dense agglomerate corresponds to N = N st (i.e., n st = 1, n so = 0 and n no = 0 in Equation (4)). From Equation (5), the maximum temperature elevation produced by such an agglomerate is given by [21]: Therefore, the maximum temperature of a N st -agglomerate is directly related to the dense agglomerate constituted by N st nanoparticles (i.e., related to the number of nanoparticles participating only in a strong coupling: electromagnetic and thermic coupling). We achieve FEM simulations of the elevation of temperature produced by N = N st identical nanoparticles for 1 ≤ N st ≤ 10 identical nanoparticles and diameters of particles D np varying in (15; 180) nm by a step of 5 nm. We solve, numerically, the photothermic problem of dense agglomerates and analyze the maximum temperature as function of the number of nanoparticles (see Figure 2). The temperature increases and is larger than N st times the temperature elevation produced by one uncoupled nanoparticle (Equation (6)). We can note that the larger the diameter D np , the smaller the N st -agglomerate satisfying a maximum temperature overpassing the thermal ablation threshold, but the higher the maximum temperature that could overpass a reasonnable maximum temperature limit T lim (e.g., 55-60 • C) [31]. On the other hand, overpassing the thermal ablation threshold with small D np requires a larger N stagglomerate, but reduces the risk of overpassing a reasonable maximum temperature limit. See, for example, the temperature elevations for agglomerates of identical nanoparticles with diameters of D np = 100 nm and D np = 80 nm (the two magenta vertical dashed lines in Figure 2). The N st -agglomerate of D np = 100 nm overpasses the ablation threshold with only N st > 1, but overpasses the temperature limit with only N st > 3. This contrasts with a N st -agglomerate of D np = 80 nm, which is overpassing the ablation threshold with N st > 2 and which is still under the temperature limit with N st < 10. Therefore, the calculation of the temperature as a function of N st is useful.
For the dense agglomerate, the electromagnetic and thermic coupling has a tremendous influence on temperature elevation, on the risk of inactive PA (for small concentrations or small agglomerates/particles), or of PA causing damage to the surrounding biological media/cells (for large concentrations or large agglomerates/particles). This is depends mainly on the incident laser power, the size of the nanoparticles, and the concentration of such nanoparticles. For a fixed incident laser power, a balance between the size of the nanoparticles and the concentration of nanoparticles (i.e., associated with the production of the N st -agglomerate that induces a maximum temperature) should be found. This balance can ensure an active PA (overpassing the thermal ablation threshold), limit the risk of maximum temperature, and protect the neighboring biological media/cells. However, in real life, for a given concentration of nanoparticles, the control of agglomerates is impossible. Therefore, it becomes necessary to consider the mixture of multiple agglomerates of various kinds of coupling.

Influence of the Concentration on the Mixture of Multiple Agglomerates
By increasing the concentration of the nanoparticles, the probability of producing a mixture of various agglomerates increases. Such a mixture of agglomerates will induce a temperature elevation that is expected to be larger than the temperature elevation of a single strong coupling N st -agglomerate. In order to protect the neighboring biological media/cells, the control of the maximum temperature, relative to the mixture of agglomerates, must be achieved. Therefore, numerical simulations for various mass concentrations (or volumic concentrations) of nanoparticles are performed in order to analyze the concentration effects on the maximum temperature elevation and on its spatial extension. The range of the relative mass concentrations are given in Table 1. For fixed relative mass concentrations C rm , multiple configurations of agglomerates can occur, mixing strong, soft, and no thermic coupling (i.e., corresponding to various occurences n st ≥ 1, n so ≥ 0, and n no ≥ 0). Figure 3 shows two examples of agglomerate configurations and temperature maps for two relative mass concentrations, C rm = 11.8% (a), (c), and C rm = 1.5% (b), (d), respectively. These two concentrations are relative to the local volume of the biological media/cells V bm (illustrated in Figure 3a,b). In these examples, the agglomerate is constituted by mixtures of N st varying in (1; 6) for the relative mass concentration C rm = 11.8%, and with N st in (1; 4) for the relative mass concentration C rm = 1.5%, respectively. As expected (Equation (5)), for a mixture of agglomerates, the maximum temperature T N,max is larger than the maximum temperature obtained with a single strong-coupling N st -agglomerate (i.e., T N,max > T N st ,max ). Here, the maximum is T N,max = 53.2 • C ≥ T N st ,max = 49.9 • C (with N st = 6) for C rm = 11.8% and T N,max = 45.8 • C ≥ T N st ,max = 44.2 • C (with N st = 4) for C rm = 1.5%.  Figure 3. Agglomerates of nanoparticles (blue disks) for two relative mass concentrations, (a) C rm = 11.8% (i.e., relative volumic concentration C rv = 0.66%) and (b) C rm = 1.5% (i.e., relative volumic concentration C rv = 0.08%), and their respective temperature maps (c,d).
Figures 4a-c and 5a-c show the temperature variations above the ablation threshold as a function of the extensions along the x-, y-, and z-axes for the two mixtures of agglomerates with C rm = 11.8% and C rm = 1.5%, respectively.    The zones of higher temperatures reveal the strong coupling for dense agglomerates and the zones of soft coupling revealed in the vicinity of multiple N st -agglomerates (shown in the magenta circles). The larger the N st value, the higher the maximum temperature. The effects of soft couplings (thermic coupling) are more easily demonstrated due to lower maximum temperatures (i.e., with small N st values). These lower temperatures are related to the fact that the strong couplings (electro-thermic couplings) are less present for smaller relative concentrations. By controlling the relative mass concentration C rm (or relative volumic concentration C rv ) of the agglomerate, it becomes possible to limit the range of maximum temperature T N st ,max produced by such agglomerates.
For both minimizing such a risk of overpassing the upper temperature limit T lim (inducing damage to the neighboring biological media/cells) and ensuring an efficiency to the PA above the ablation threshold, the relevant criteria are on the minimum and maximum numbers of dense agglomerates (i.e., only considering the strong coupling), relative to the mass concentration C rm . Therefore, N st min and N st max must satisfy T N st min ≥ T abl and T N st max ≤ T lim , and their values can be obtained by: where Int(.) is the integer part. With such N st min and N st max , a range of relative mass concentrations can be managed with a probability of succes. The probabilities of producing agglomerates are deduced from statistical analyses obtained from multiple relative mass concentrations and random generations of the spatial distributions of the agglomerates. We use 100 relative mass concentration values, varying in (0.02; 35.0)%, and for each value, 1000 random spatial distributions are produced for the agglomerates before solving the photo-thermic problem. With these random generations, statistical analysis of the N st -agglomerates permit assigning probabilities of producing N st -agglomerates and their uncertainties. Figure 6a shows the probability of producing an N-agglomerate as a function of the relative mass concentration C rm for various strong couplings. As expected, for a given concentration, the larger the N st , the smaller the probability of producing such an N st -agglomerate. Such a probability increases with the relative mass concentration. Figure 6b gives the ranges of temperatures (minimum and maximum) as a function of the relative mass concentration C rm for mixtures of agglomerates and for three diameters D np of nanoparticles. The larger the diameter D np of nanoparticles, the larger the range of temperatures and the higher the maximum temperature (that could overpass the maximum temperature limit).
The obtained results are consistent with the expected behavior. Indeed, the trends of Figure 6a,b are similar in C rm and to the sizes of the nanoparticles forming the mixture of agglomerates: an increasing concentration induces an increase of the probability of producing at least one agglomerate with N st ≥ N st min . The larger the diameter D np of the nanoparticles forming the mixture of agglomerates, the higher the maximum temperature that rapidly overpasses the temperature limit. The main parameter governing the efficiency of the PA and limiting the risk of temperature limit T lim is directly associated with the production of one N st min -agglomerate for a given concentration. A compromise between the sizes of the nanoparticles and the relative mass concentration must be found. Here, D np = (85.0 ± 5.0) nm and C rm = (2.0 ± 1.0)% ensure an active PA with a realistic concentration by limiting the exceeding of the temperature limit.

Discussion of Material Parameters and Configuration
In cancer therapy applications, the laser wavelengths are in the biological window (800 to 1000 nm). In such a window, the imaginary part of the permittivity of the gold nanoparticles constituting the agglomerates slightly increases with the incident wavelength, necessitating the designing of the particle shapes forming the agglomerates or the decreasing of the laser power. For individual nanoparticles, our results are in agreement with the experimental ones in Refs. [3][4][5][6][7] (the absorption efficiencies, due to the shapes of the nanoparticles, are counterbalanced by smaller laser power densities, producing compatible heat sources and temperature elevations). Moreover, the treatment of larger biological media/cells (typically a cell size is varying from 1 to 100 µm [32]) or a mix of smaller/larger concentrations of nanoparticles, a smaller/greater power opens the way to control both the PA efficiency (i.e., the local temperature elevation above the ablation threshold) and the protection of the neighboring region of the biological media/cells (i.e., limiting the risk of higher temperature elevation). Each nanoparticle acts as a local temperature source, and in a mixture of agglomerates of nanoparticles, multiple local thermic sources are induced simultaneously. Therefore, the global thermic field is increased and the main effect in the temperature elevation is directly related to agglomerates inducing a strong photo-thermic coupling. Therefore, as the maximum temperature is directly related to the characteristics of the agglomerates (i.e., to the size of the nanoparticles and to the number of nanoparticles participating in the agglomerate and inducing strong coupling), two configurations are to be favored: either promote the use of large nanoparticles in very low concentrations, or use small nanoparticles in high concentrations. Each of these two configurations has advantages and disadvantages, which can be minimized (e.g., adjusting the laser power) and which depend on the use conditions in real life. The use of large nanoparticles in low concentrations permits the production of an active PA with a limited risk of agglomerates exceeding the temperature limit. Conversely, the use of small nanoparticles in high concentrations limits the risk of producing an inactive PA and avoids the risk of exceeding the temperature limit.

Conclusions
In this study, we focus on an analysis of the efficiency and the safety of the mixture of agglomerates of nanoparticles as photo-thermic agents embeded in biological media/cells (i.e., their ability to produce a temperature elevation above the ablation threshold and to limit the risk for the neighboring biological media/cells). We show that the temperature elevation strongly depends on both the nanoparticle size and concentration, with the dense agglomerates inducing strong photo-thermic coupling between the nanoparticles. Agglomerates of nanoparticles are much more efficient than individual nanoparticles, but induce a potential risk for the surrounding media/cells. By tuning the mass concentration (or volumic concentration), the maximum temperature can be controlled. Here, agglomerates constituted by only spherical particles are considered. Nevertheless, FEM with an adaptive remeshing process is a general method. Therefore, FEM could also be adapted to treat more complex and realistic systems (materials, sizes, shapes), and for time-dependent analysis for temperature evolution under pulse illumination. Moreover, advanced heuristic laws of agglomerate formation, material functionalization improving the binding to the cancer cell, as well as experimental distributions of particles sizes, could also be introduced. Acknowledgments: This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: with the velocity of light denoted by c and the curl operator by (∇×). To solve Equation (A1), a radiation boundary condition (Equation (A2)) is applied on the external border Γ of the computational domain Ω [20,22]: with the normal vector to the surfacen, i 2 = −1, and the incident electromagnetic field E i . Such a boundary condition is used to describe the free propagation of the total electromagnetic field. The amplitude of the illumination field is given by E 0 = 2cµ 0 P s , with the permeability of a vacuum µ 0 and the surface laser power density P s . Such a surface laser power density P s is a function of the laser power P w and the diameter of the laser beam D beam (i.e., P s = 4P w /(πD 2 beam )). For the thermic problem, the temperature elevation is governed by the heat equation (Equation (2)) with the external boundary condition (Equation (A3)): with the initial temperature T 0 and the heat transfer coefficient α. Such an external boundary is the Fourier boundary and the initial condition T(x, 0) = T 0 [33].

Appendix A.2. FEM Numerical Method
The main principle of FEM is to solve the problem on a discrete domain represented by a mesh [22,34,35]. The unknown physics quantities (here, the electromagnetic field E and the temperature T) are computed at the nodes of the mesh through a variational method. The use of an improved 3D method and an iterative remeshing process [18][19][20], ensures the control of the error on the numerical solution. The computed solution accuracy is directly related to the quality of the mesh [36,37]. Such a quality of the mesh is automatically improved by adapting the size of the mesh elements to the solution of the physics through adaptive loops of the remeshing process [19,38]. Such mesh adaptions are the only requirements to ensure the convergence of the numerical solution to a stable solution [18,19]: the approximations of the electromagnetic field E and the temperature T (Equation (2)) are computed, at each step of the adaption process, on the basis of the interpolation polynomials between nodes. The difference between the computed solution at the nodes and the maximum deviation between the solution associated with the mesh is controlled by the interpolation error. The interpolation error is based on an estimation of the discrete Hessian of the solution) [39,40]. The adaptive remeshing is obtained from the OPTIFORM software (generator of isotropic/anisotropic remeshing) [18,19] and, during the iterative process, the domain is entirely remeshed and a new mesh is produced. In the present case, the volume of the computational domain Ω is 100 µm 3 . The remeshing of the domain is only governed by the minimum and maximum element sizes, set to h min = 0.1 nm, h max = 500 nm, and with a maximum tolerance δ = 0.2% for the relative error of the computed unknowns.