Influence of natural convection on gold nanorods-assisted photothermal treatment of bladder cancer in mice

Abstract

Background

The thermally-induced urine flow can generate cooling that may alter the treatment outcome during hyperthermic treatments of bladder cancer. This paper investigates the effects of natural convection inside the bladder and at skin surface during gold nanorods (GNR) - assisted photothermal therapy (PTT) of bladder cancer in mice.

Methods

3D models of mouse bladder at orientations corresponding to the mouse positioned on its back, its side and its abdomen were examined. Numerical simulations were carried out for GNR volume fractions of 0.001, 0.005 and 0.01% and laser power of 0.2 and 0.3 W.

Results

The obtained results showed that cooling due to natural convection inside the bladder and above the skin depends on the mouse orientation. For a mouse positioned on its back, on its side or on its abdomen, the maximum temperature achieved inside the tumour at 0.001% GNR volume fraction and 0.2 W laser power was 55.2°C, 50.0°C and 52.2°C, respectively compared to 56.8°C when natural convection was not considered. The average thermal gradients when natural convection was considered were also lower, suggesting a more homogenous temperature distribution.

Conclusions

Natural convection inside the bladder can be beneficial but also detrimental to GNR-assisted PTT depending on the level of heating. At low levels of heating due to low GNR volume fraction and/or laser power, flow inside the bladder may dissipate heat from the targeted tissue; making the treatment ineffective. At high levels of heating due to high GNR volume fraction and/or laser power, cooling may prevent excessive thermal damage to surrounding tissues.

Keywords:

• Gold nanorods
• plasmonic photothermal ablation
• near-infrared laser
• convection heat transfer
• Monte Carlo method

1. Introduction

Nanoparticle-assisted photothermal therapy (PTT) is a thermal ablation technique that uses nanoparticles as photoabsorbers to generate heat within the tumor tissue upon laser irradiation [1]. Gold nanorods (GNR) represent a suitable choice of photoabsorbers due to their tunable optical properties. By altering the aspect ratio, GNR can be tuned to exhibit peak absorbance over a broad spectral range [2]. This is particularly useful in biological applications as GNR can be tuned to allow peak absorption in the near-infrared (NIR) region (600–1300 nm), also known as the biological window, where light penetration through tissue is at its maximum [3].

Of interest in this study is the treatment of bladder cancer using GNR-assisted PTT. Bladder cancer is the most common form of cancer of the urinary tract and is the ninth most common cancer worldwide [4,5]. Treating bladder cancer is difficult mainly due to an incomplete understanding of the disease biology and the limited availability among current therapeutic systems that can effectively eradicate the disease [6]. Recent studies have confirmed GNR-assisted PTT as an effective method for treating bladder cancer [7–13]. Nevertheless, these studies were carried out in vitro and using animal models as a substitute for the human bladder, with research leading to potential clinical trials still ongoing. Under these circumstances, numerical simulations can play an important role in supplementing vital information toward the advancement of this technique, such as the optical and thermal responses of tissues during treatment.

Models of nanoparticle-assisted PTT have been developed in the past to elucidate the role of various parameters influencing the treatment. For instance, Soni et al. [14,15] developed computational models to investigate the influence of optical coefficients and inhomogeneous GNR distribution inside the tumor on the tissue heating characteristics. Singh et al. [16] investigated the effects of different laser parameters during the treatment of vascularized tissues using gold nanoshells-mediated PTT. Jeynes et al. [17] examined the effectiveness of GNR-assisted PTT for eradicating skin cancer, while Manuchehrabadi and Zhu [18] focused on the development of computational models that facilitated the protocol design for treating prostate cancer.

To ensure that the results obtained from the numerical simulations are reliable, the physics describing the different processes that are involved during the treatment must be given careful considerations. In the case of GNR-assisted PTT of bladder cancer, natural convection inside the bladder due to the thermally-driven flow of urine can potentially affect the accuracy of the model predictions if they are not properly accounted for during the simulations. Studies carried out on regional hyperthermic treatment of the bladder have demonstrated how natural convection arising from the thermally driven flow of urine helps the distribution of heat across the entire bladder [19–21]. However, unlike regional hyperthermia where the whole organ undergoes temperature elevation, heating during GNR-assisted PTT is localized to the region of tumor targeted by the laser beam. As such, the role of natural convection during GNR-assisted PTT may be different from that of regional hyperthermia.

Motivated by this, the present study seeks to develop numerical models of mouse bladder to investigate the role of natural convection during bladder cancer treatment with GNR-assisted PTT. The focus on mouse bladder is deliberate since majority of the experimental studies on GNR-assisted PTT at the tissue level were carried out in mice [10,12,13]. As such, the results obtained from the present study may contribute directly to the ongoing development of the treatment. Different bladder orientations due to the different positions of the mouse during treatment, such as laying on its abdomen, on its back, and on its side, are investigated. Since the different positions of the mouse during treatment also affect the flow of air around the skin, the present study takes into account the thermally-induced air flow and its impact on bladder cancer treatment using GNR-assisted PTT. For each orientation, the effects of different GNR volume fraction and laser power on the treatment outcome are examined. The role of natural convection is determined by monitoring the temperature distribution and the formation of thermal damage during the treatment.

2. Theory and methods

2.1. Model geometry

The model of mouse bladder was constructed based on dimensions obtained through visual examination of the photographic images reported by Reis et al. [22]. The bladder was modeled as an ovoid with dimensions shown in Figure 1 with a uniform thickness of 0.35 mm [23]. It was assumed to be filled with urine at the time of treatment. The bladder model has a volume of 0.1 ml, which is similar to the actual bladder capacity in mice [22]. The tissues surrounding the bladder were modeled as an anatomically homogeneous cuboid of 16 mm width (x-direction), 5.7 height (y-direction), and 10 mm depth (z-direction). The distance from the top surface of the bladder to the skin is approximately 1.2 mm, which is within the range of mouse skin thickness [24]. The natural convection above the skin is modeled using a cuboid with the following dimensions: 16 × 4.05 × 10 mm³.

Figure 1. Geometry of the mouse bladder and the three orientations investigated. (a) The 3D computational domain developed in the present study. The bladder is shown in cyan, while the tumor is shown in red. (b) View of the computational model across the z = 0 plane. The cyan domain represents the urine. (c) The three different bladder orientations considered in the present study with respect to direction of gravity and laser irradiation. All dimensions are in mm.

The 3 D model used is shown in Figure 1(a). A tumor with dimensions of 1.96 and 2 mm in the x- and z-directions, respectively, and with a thickness equivalent to that of the bladder was assumed to grow at the top side of the bladder. This is shown in Figure 1(b). The resulting volume of the tumor domain was 0.991 mm³. A continuous wave laser with a flat top beam profile of radius 3 mm irradiates the skin surface along the centerline of the bladder.

The effects of natural convection during GNR-assisted PTT were investigated for three cases that differ by the bladder orientation. The first case (Case I) has the mouse laying on its back. The second case (Case II) has the mouse laying on its abdomen. The third case (Case III) has the mouse laying on its side. These different bladder orientations are illustrated in Figure 1(c). A fourth case (Case IV), which represented the model without natural convection (orientation-independent), was also considered.

2.2. Monte Carlo simulations

In this study, the Monte Carlo method, considered to be the benchmark algorithm for simulating light propagation in biological tissues [25], was used to simulate the tissue optical response during laser irradiation. Details on the implementation of the Monte Carlo algorithm are presented briefly as the implementation has been reported elsewhere [26]. In a typical Monte Carlo algorithm, packets of photons representing the laser beam propagate into the computational domain. Physical processes such as absorption, scattering, transmission, and reflection of the photons at the boundaries and interfaces are described as probability distributions. Typically, millions of photons are used in the algorithm in order to produce an accurate representation of light propagation inside the tissue.

In this study, the Monte Carlo simulations were carried out using ValoMC, an open source mesh-based Monte Carlo algorithm (https://inverselight.github.io/ValoMC/) [27]. With ValoMC, the simulation geometry can be defined using unstructured mesh (triangles in 2D and tetrahedrons in 3D). This allows for a straightforward transfer of information obtained from the Monte Carlo simulations to the finite element heat transfer simulations (see Section Heat transfer model). In order to carry out the Monte Carlo simulations using ValoMC, the algorithm requires four inputs, namely the absorption coefficient (µ_a), the scattering coefficient (µ_s), the scattering anisotropic factor (g), and the refractive index (n), for each domain of the computational model. The algorithm provides the spatial distribution of light fluence, which is used to calculate the amount of heat generated inside the tissue due to laser energy absorption.

2.3. Heat transfer model

Heat transfer inside the tumor, the bladder and the surrounding tissue was described using the Pennes bioheat equation [28,29]: (1) $${\rho }_{\mathrm{t}}{c}_{\mathrm{t}}\frac{\partial T_{\mathrm{t}}}{\partial t}=\nabla ·\left(k_{\mathrm{t}}\nabla T_{\mathrm{t}}\right)+{\rho }_{\mathrm{b}}{c}_{\mathrm{b}}{\omega }_{\mathrm{b}}\left(T_{\mathrm{b}}-T\right)+Q_{\mathrm{m}}+{\mu }_{a,t}{\mathrm{\Phi }}_t{P}_{\mathrm{laser}}$$(1) where $T$ is temperature, $t$ is time, $\rho ,$ $c,$ and $k$ are density, specific heat and thermal conductivity, respectively,${\omega }_b$ is blood perfusion rate, $T_{\mathrm{b}}$ is temperature of the arterial blood that is assumed to be at body temperature of 37 °C, $Q_{\mathrm{m}}$ is tissue metabolic heat generation, ${\mu }_{\mathrm{a},\mathrm{t}}$ is optical absorption coefficient, $\mathrm{\Phi }\mathrm{ }$ is the light fluence distribution obtained from the Monte Carlo simulations and $P_{\mathrm{laser}}$ is laser power. The subscripts ‘$\mathrm{t}$’ and ‘$\mathrm{b}$’ represent tissue and blood, respectively. In Equation (1)(1) $${\rho }_{\mathrm{t}}{c}_{\mathrm{t}}\frac{\partial T_{\mathrm{t}}}{\partial t}=\nabla ·\left(k_{\mathrm{t}}\nabla T_{\mathrm{t}}\right)+{\rho }_{\mathrm{b}}{c}_{\mathrm{b}}{\omega }_{\mathrm{b}}\left(T_{\mathrm{b}}-T\right)+Q_{\mathrm{m}}+{\mu }_{a,t}{\mathrm{\Phi }}_t{P}_{\mathrm{laser}}$$(1) , the second and third terms on the right-hand side represent the contributions due to blood flow and tissue metabolism, respectively.

Heat transfer inside the urine and bladder was described using: (2) $${\rho }_{\mathrm{f}}{c}_{\mathrm{f}}\frac{\partial T_{\mathrm{f}}}{\partial t}+\mathbf{u}·\nabla T_{\mathrm{f}}=\nabla ·\left(k_{\mathrm{f}}\nabla T_{\mathrm{f}}\right)+{\mu }_{\mathrm{a},\mathrm{f}}{\mathrm{\Phi }}_{\mathrm{f}}{P}_{\mathrm{laser}},$$(2) where u = (u, v, w) is the vector describing the velocity of fluid flow in bladder and air in the x-, y-, and z-directions, and the subscript ‘$\mathrm{f}$’ represents either the urine or the air. In the air domain, the last term on the right-hand-side of Equation (2)(2) $${\rho }_{\mathrm{f}}{c}_{\mathrm{f}}\frac{\partial T_{\mathrm{f}}}{\partial t}+\mathbf{u}·\nabla T_{\mathrm{f}}=\nabla ·\left(k_{\mathrm{f}}\nabla T_{\mathrm{f}}\right)+{\mu }_{\mathrm{a},\mathrm{f}}{\mathrm{\Phi }}_{\mathrm{f}}{P}_{\mathrm{laser}},$$(2) was set to zero.

By assuming the flow to be laminar, urine and air to be Newtonian fluids, and ignoring viscous dissipation effects, the velocity vector u in Equation (2)(2) $${\rho }_{\mathrm{f}}{c}_{\mathrm{f}}\frac{\partial T_{\mathrm{f}}}{\partial t}+\mathbf{u}·\nabla T_{\mathrm{f}}=\nabla ·\left(k_{\mathrm{f}}\nabla T_{\mathrm{f}}\right)+{\mu }_{\mathrm{a},\mathrm{f}}{\mathrm{\Phi }}_{\mathrm{f}}{P}_{\mathrm{laser}},$$(2) can be obtained from solving the Navier-Stokes equations: (3) $${\rho }_{\mathrm{f}}\frac{\partial \mathit{u}}{\partial t}+{\rho }_{\mathrm{f}}\left(\mathit{u}·\nabla \mathit{u}\right)=-\nabla p+{\mu }_{\mathrm{f}}{\nabla }^2\mathbf{u}+{\rho }_{\mathrm{f}}\left(T_{\mathrm{f}}\right)\mathbf{g},$$(3) (4) $${\rho }_{\mathrm{f}}\nabla ·\mathbf{u}=0,$$(4) where ρ_f and ${\mu }_{\mathrm{f}}$ are the density and dynamic viscosity, respectively, of urine and air, $p$ is pressure and g is a vector representing gravitational acceleration. The last term on the right-hand side Equation (3)(3) $${\rho }_{\mathrm{f}}\frac{\partial \mathit{u}}{\partial t}+{\rho }_{\mathrm{f}}\left(\mathit{u}·\nabla \mathit{u}\right)=-\nabla p+{\mu }_{\mathrm{f}}{\nabla }^2\mathbf{u}+{\rho }_{\mathrm{f}}\left(T_{\mathrm{f}}\right)\mathbf{g},$$(3) represents the thermally induced buoyant forces causing convective flow, which can be described by the Boussinesq approximation: (5) $${\rho }_{\mathrm{f}}\left(T_{\mathrm{f}}\right)={\rho }_{\mathrm{ref}}\left[1-{\beta }_{\mathrm{f}}(T_{\mathrm{f}}-T_{\mathrm{ref}})\right],$$(5) where ${\beta }_{\mathrm{f}}$ is the fluid thermal expansion coefficient and ${\rho }_{\mathrm{ref}}$ is the reference density evaluated at reference temperature $T_{\mathrm{ref}},$ which in the present study, is equal to body temperature in the urine domain and equal to ambient temperature in the air. Equation (5)(5) $${\rho }_{\mathrm{f}}\left(T_{\mathrm{f}}\right)={\rho }_{\mathrm{ref}}\left[1-{\beta }_{\mathrm{f}}(T_{\mathrm{f}}-T_{\mathrm{ref}})\right],$$(5) shows a decrease in urine density of ∼11 kgm⁻³ when temperature increases from 37 °C to 70 °C (see Section Results).

2.4. Thermal damage model

The formation of thermal coagulation inside both healthy and tumor tissues were described using the thermal damage model of Henriques and Moritz [30]. The model assumes damage to biological tissues due to heat is similar to thermal denaturation, which follows the first-order approximation of the Arrhenius equation. The rate of the thermal denaturation process is given by: (6) $$\mathrm{\Omega }(x,y,z,t)=\underset{0}{\overset{t_{\mathrm{f}}}{\int }}\mathit{Aexp}\left(-\frac{\Delta E}{R_{c}T(x,y,z,t)}\right)dt,$$(6) where $\Omega$ is a dimensionless parameter that describes thermal damage of healthy and tumor tissues at a given point in space ($x,y,z$) and time, $A$ is frequency factor, $\Delta E$ is the activation energy for irreversible thermal denaturation, and $R_{\mathrm{c}}$ is the universal gas constant.

Once Equation (6)(6) $$\mathrm{\Omega }(x,y,z,t)=\underset{0}{\overset{t_{\mathrm{f}}}{\int }}\mathit{Aexp}\left(-\frac{\Delta E}{R_{c}T(x,y,z,t)}\right)dt,$$(6) is solved at all the given points in space, the probability of observing thermal damage, $\mathrm{PD}$ can be estimated using: (7) $$\mathrm{PD}=100\%\times \left[1-\exp (-\Omega )\right].$$(7)

According to Equation (7)(7) $$\mathrm{PD}=100\%\times \left[1-\exp (-\Omega )\right].$$(7) , values of $\Omega$ = 0.7, 1.1, and 4.6 would yield a probability of observing thermal damage at a given point inside the tissue of 50, 67, and 99%, respectively. In the present study, it was assumed that all the points within the tumor must have values of $PD$ above 99%, which is equivalent to $\Omega$ > 4.6, in order to eradicate the tumor.

2.5. Initial-boundary conditions

Since the timescale for light propagation was assumed to be much smaller than the timescale for heat transfer, steady state Monte Carlo simulations were employed. The initial temperature and flow profile inside the bladder were obtained by solving Equations (1)–(5)(1) $${\rho }_{\mathrm{t}}{c}_{\mathrm{t}}\frac{\partial T_{\mathrm{t}}}{\partial t}=\nabla ·\left(k_{\mathrm{t}}\nabla T_{\mathrm{t}}\right)+{\rho }_{\mathrm{b}}{c}_{\mathrm{b}}{\omega }_{\mathrm{b}}\left(T_{\mathrm{b}}-T\right)+Q_{\mathrm{m}}+{\mu }_{a,t}{\mathrm{\Phi }}_t{P}_{\mathrm{laser}}$$(1) at steady-state, i.e., by setting the time derivative term in the left-hand side of Equations (1)–(3)(1) $${\rho }_{\mathrm{t}}{c}_{\mathrm{t}}\frac{\partial T_{\mathrm{t}}}{\partial t}=\nabla ·\left(k_{\mathrm{t}}\nabla T_{\mathrm{t}}\right)+{\rho }_{\mathrm{b}}{c}_{\mathrm{b}}{\omega }_{\mathrm{b}}\left(T_{\mathrm{b}}-T\right)+Q_{\mathrm{m}}+{\mu }_{a,t}{\mathrm{\Phi }}_t{P}_{\mathrm{laser}}$$(1) to be zero and by letting ${\Phi }_{\mathrm{t}}$ = ${\mathrm{\Phi }}_{\mathrm{f}}$= 0, where the latter specifies the condition prior to laser irradiation.

An open thermal boundary condition was applied to the surfaces bounding the air domain (except for the solid surface in Case III, see Figure 1(c)): (8) $$\{\begin{array}{c}T=T_{\mathrm{amb}},\mathit{\text{ if }}\mathit{u}·\mathit{n}<0\\ -k_{\mathrm{f}}\frac{\partial T_{\mathrm{f}}}{\partial n}=0,\mathit{\text{ if }}\mathit{u}·\mathit{n}\ge 0\end{array},$$(8) where $T_{\mathrm{amb}}$ is ambient temperature set in this study to 20 °C. Equation (8)(8) $$\{\begin{array}{c}T=T_{\mathrm{amb}},\mathit{\text{ if }}\mathit{u}·\mathit{n}<0\\ -k_{\mathrm{f}}\frac{\partial T_{\mathrm{f}}}{\partial n}=0,\mathit{\text{ if }}\mathit{u}·\mathit{n}\ge 0\end{array},$$(8) assumes the surfaces of the bounding air domain to be at ambient temperature if air flows into the domain and to have a zero heat flux if air flows out of the domain. In Case III, zero heat flux condition was prescribed at the surface on which the mouse was placed (see Figure 1(c)). The surface at the opposite side of the skin surface was assumed to be sufficiently far from the laser focal zone such that normo-thermoregulation can maintain the temperature at basal level. Therefore, the boundary condition here is given by: (9) $$T_{\mathrm{t}}=T_{\mathrm{f}}=37°\mathrm{C},$$(9)

For the hydrodynamic model, open boundary condition that allows air to flow in and out of the domain was prescribed across the surfaces that bound the air domain except for the solid surface representing the solid side on which the mouse lay in Case III (see above). Across the solid surface (in Case III) and the inner surfaces of the bladder, the no slip boundary condition was prescribed, such that: (10) $$\mathbf{u}=0.$$(10)

2.6. Material properties

2.6.1. Thermal properties

The thermal properties employed in this study are listed in Table 1. Except for the tumor domain, the values presented in Table 1 were obtained from IT’IS database [31]. As there was lack of information on the thermal properties of bladder cancer, thermal properties of liver cancer, which are available in the literature, were used.

Table 1. Thermal properties employed in this study.

Parameter	Value	Source
Bladder
Thermal conductivity, $k$ (W/m^2K)	0.52	[31]
Density, $\rho$ (kg/m^3)	1086	[31]
Specific heat, $c$ (J/kg K)	3581	[32]
Blood perfusion rate, ${\omega }_{\mathrm{b}}$ (1/s)	0.0014	[31]
Frequency factor, $A$ (1/s)	3.3 × 1038	[33]
Activation energy, $\Delta E$ (J/mol)	2.57 × 105	[33]
Tumor
Thermal conductivity, $k$ (W/m^2K)	0.49	[32]
Density, $\rho$ (kg/m^3)	1090	[32]
Specific heat, $c$ (J/kgK)	3421	[32]
Blood perfusion rate, ${\omega }_{\mathrm{b}}$ (1/s)	0.00165	[32]
Frequency factor, $A$ (1/s)	1.46 × 1052	[34]
Activation energy, $\Delta E$ (J/mol)	3.43 × 105	[34]
Surrounding tissue
Thermal conductivity, $k$ (W/m^2K)	0.21	[31]
Density, $\rho$ (kg/m^3)	911	[31]
Specific heat, $c$ (J/kg K)	2348	[31]
Blood perfusion rate, ${\omega }_{\mathrm{b}}$ (1/s)	0.0005	[31]
Frequency factor, $A$ (1/s)	1.61 × 1045	[35]
Activation energy, $\Delta E$ (J/mol)	3.06 × 105	[36]
Urine
Thermal conductivity, $k$ (W/m^2K)	0.59	[31]
Density, $\rho$ (kg/m^3)	1024	[31]
Specific heat, $c$ (J/kg K)	4178	[31]
Thermal expansion coefficient, ${\beta }_{\mathrm{f}}$ (1/K)	3.374 × 10^-4	[37]
Dynamic viscosity, ${\mu }_{\mathrm{f}}$ (Pa s)	0.0007	[37]
Constants
Universal gas constant, $R_{\mathrm{c}}$ (J/mol K)	8.314
Ambient temperature, $T_{\mathrm{amb}}$(^oC)	20
Ambient convection coefficient, $h_{\mathrm{amb}}$ (W/m^2K)	3

Except for thermal conductivity and blood perfusion rate, all other thermal properties employed in the model were assumed to be constant and homogeneous. The thermal conductivity was assumed to increase linearly with temperature at a rate of 1.5% per °C [38]. This assumption was applied to all the domains including the urine.

Studies carried out on porcine kidneys have demonstrated a nonlinear variation in blood perfusion with temperature [33]. At temperatures below 45 °C, tissue responds to temperature rise by increasing blood flow through vasodilation [39]. As temperature continues to increase, blood vessels begin to lose their caliber and vascular stasis starts to occur, thus leading to a decrease in blood perfusion. Cessation of blood flow occurs when complete thermal damage is sustained by the tissue. This temperature-dependent behavior of blood perfusion can be expressed as a piecewise homogeneous function of $\mathrm{PD}$ given by [33,40]: $${\omega }_b\left(PD\right)={\omega }_{b,\mathit{ref}}\times \{\begin{array}{c}0.3PD+1,\\ \begin{array}{c}-0.13PD+1.86,\\ \begin{array}{c}-0.0079PD+0.884,\\ -0.0387PD+3.87,\end{array}\end{array}\end{array}\begin{array}{c}\begin{array}{c}PD\le 2\mathrm{\%}\\ 2\mathrm{\%}<PD\le 8\mathrm{\%}\end{array}\\ \begin{array}{c}8\mathrm{\%}<PD\le 97\mathrm{\%}\\ 97\mathrm{\%}<PD\le 100\mathrm{\%}\end{array}\end{array},$$ where ${\omega }_{\mathrm{b},\mathrm{ref}}$ represents the blood perfusion rate at body temperature.

In this study, Equation (11)(1) $${\rho }_{\mathrm{t}}{c}_{\mathrm{t}}\frac{\partial T_{\mathrm{t}}}{\partial t}=\nabla ·\left(k_{\mathrm{t}}\nabla T_{\mathrm{t}}\right)+{\rho }_{\mathrm{b}}{c}_{\mathrm{b}}{\omega }_{\mathrm{b}}\left(T_{\mathrm{b}}-T\right)+Q_{\mathrm{m}}+{\mu }_{a,t}{\mathrm{\Phi }}_t{P}_{\mathrm{laser}}$$(1) was employed to describe the blood perfusion inside the bladder and the surrounding tissue, but not inside the tumor. This is due to the limited ability of tumor tissues to increase blood flow in response to hyperthermia [39]. As an alternative, a step function given by: (12) $${\omega }_{\mathrm{b}}\left(\mathrm{PD}\right)={\omega }_{\mathrm{b},\mathrm{ref}}\times \{\begin{array}{cc}1,& \mathrm{P}\mathrm{D}<100\%\\ 0,& \mathrm{P}\mathrm{D}=100\%\end{array},$$(12) was used to express the variation in tumor blood perfusion with thermal damage.

2.6.2. Optical properties

The absorption and scattering coefficients of the bladder and the surrounding tissue were obtained from the literature and they are listed in Table 2. The scattering anisotropic factors and the refractive indices of all tissues were chosen to be 0.9 and 1.4, respectively [44]. For the urine, the absorption coefficient, scattering coefficient, and refractive index were assumed to be the same as those of water. The scattering anisotropic factor for the urine was chosen to be 1, implying that photons propagate inside the fluid domain but do not undergo scattering [43].

Table 2. Optical parameters and the peak absorption and scattering coefficients for $\varphi$ = 0.001, 0.005 and 0.01%.

Parameter	Value	Source
Free electron density, $n_o$ (1/m^3)	5.9 × 1028	^a
Permittivity of vacuum, ${\epsilon }_o$ (F/m)	8.85 × 10^-12	^a
Electron charge, $e$ (C)	1.6 × 10^-19	^a
Mass of electron, $m_{\mathrm{e}}$ (kg)	9.1 × 10^-31	^a
Fermi velocity, $v_{\mathrm{F}}$ (m/s)	1.39 × 106	^a
Absorption coefficient, ${\mu }_{\mathrm{a}}$(1/mm)
Bladder	0.02	[41]
Surrounding tissue	0.105	[42]
Urine	0.002	[43]
Tumor embedded with GNR
$\varphi$ = 0.001%	4.29	^b
$\varphi$ = 0.005%	21.47	^b
$\varphi$= 0.01%	42.93	^b
Scattering coefficient, ${\mu }_{\mathrm{s}}$ (1/mm)
Bladder	4.41	[41]
Surrounding tissue	11.16	[42]
Urine	1	[43]
Tumor embedded with GNR
$\varphi$ = 0.001%	0.052	^b
$\varphi$ = 0.005%	0.261	^b
$\varphi$ = 0.01%	0.521	^b

^aConstants.

^bEstimated using Mie–Gans theory.

The absorption and scattering coefficients of tumor with GNR were calculated using the Mie–Gans theory [45], which solves the Maxwell equations based on an electrostatic assumption. This assumption is valid for $D$ ≤ $\lambda$/10, where $D$ is the diameter of the GNR and $\lambda$ is the wavelength of the laser. Accordingly, the absorption and scattering coefficients of a medium embedded with GNR can be calculated using [46,47]: (13) $${\mu }_a=\frac{2\pi \varphi }{\lambda V_{np}}\mathrm{imag}\left(\frac{{\alpha }_1}{3}+\frac{{\alpha }_2}{3}+\frac{{\alpha }_3}{3}\right),$$(13) and (14) $${\mu }_s=\frac{16{\pi }^3\varphi }{18{\lambda }^4{V}_{np}}\left({\left|{\alpha }_1\right|}^2+{\left|{\alpha }_2\right|}^2+{\left|{\alpha }_3\right|}^2\right),$$(14) where $\mathrm{imag}()$ represents the imaginary component, $\varphi$ is the volume fraction of GNR, $V_{\mathrm{np}}$ is the volume of a single nanoparticle calculated by assuming the GNR to be a cylinder capped by hemispheres at both ends and ${\alpha }_i$ (for $i$ =1, 2, 3) is the polarization given by [46,47]: (15) $${\alpha }_i=V_{\mathrm{np}}\left[\frac{\epsilon \left(\omega ,D\right)-{\epsilon }_{\mathrm{m}}}{P_i\left(\epsilon \left(\omega ,D\right)-{\epsilon }_{\mathrm{m}}\right)+{\epsilon }_{\mathrm{m}}}\right],$$(15) where $\epsilon$ is the size- and frequency-dependent dielectric function of the GNR, ${\epsilon }_{\mathrm{m}}$ is the dielectric constant of the surrounding medium and $P_i$ is the geometrical factor. In the present study, the Kreibig–Vollmer model for calculating the size- and frequency-dependent dielectric function of GNR was adopted, which is given by [47]: (16) $$\epsilon \left(\omega ,D\right)={\epsilon }_{\mathrm{bulk}}+{\omega }_p^2\left(\frac{1}{{\omega }^2+{\Gamma }_o^2}-\frac{1}{{\omega }^2+{\Gamma }_d^2}\right)+i\frac{{\omega }_p^2}{\omega }\left(\frac{{\Gamma }_d}{{\omega }^2+{\Gamma }_o^2}-\frac{{\Gamma }_o}{{\omega }^2+{\Gamma }_d^2}\right),$$(16) where ${\epsilon }_{\mathrm{bulk}}$ is the frequency-dependent dielectric constant of bulk gold, and ${\omega }_p$ and ${\Gamma }_d$ are given by: (17) $${\omega }_p={\left(\frac{n_o{e}^2}{{\epsilon }_o{m}_e}\right)}^{\frac{1}{2}},$$(17) (18) $${\Gamma }_d={\mathrm{\Gamma }}_o+\frac{A_s{v}_F}{r_{\mathrm{eff}}},$$(18) (19) $$r_{\mathit{eff}}=\frac{4V_{np}}{S_{np}}$$(19) where $n_o$ is the electron density number, ${\epsilon }_o$ is permittivity of vacuum, $e$ and $m_{\mathrm{e}}$ are the charge and mass of the electrons, respectively, $v_{\mathrm{F}}$ is Fermi velocity, $A_{\mathrm{s}}$ is the surface scattering parameter with an empirical value of 0.3, $r_{\mathrm{eff}}$ is effective radius, $S_{\mathrm{np}}$ is the surface area of a single nanoparticle, and ${\Gamma }_o$ is the free damping coefficient defined as the ratio of Fermi velocity and the mean free path of electrons, which for gold is 42 nm [48]. The values of these constants are listed in Table 2. It is important to point out here that the Equation (16)(16) $$\epsilon \left(\omega ,D\right)={\epsilon }_{\mathrm{bulk}}+{\omega }_p^2\left(\frac{1}{{\omega }^2+{\Gamma }_o^2}-\frac{1}{{\omega }^2+{\Gamma }_d^2}\right)+i\frac{{\omega }_p^2}{\omega }\left(\frac{{\Gamma }_d}{{\omega }^2+{\Gamma }_o^2}-\frac{{\Gamma }_o}{{\omega }^2+{\Gamma }_d^2}\right),$$(16) is expressed in terms of electron volts (eV). The dielectric constant of bulk gold was obtained from data published by Johnson and Christy [49].

The Mie-Gans theory for calculating ${\mu }_a$ and ${\mu }_{\mathrm{s}}$ is valid only for spheroidal nanoparticles. In the present study, in order to evaluate $P_i,$ the coefficients for a quadratic fit were employed, which were derived by Prescott and Mulvaney [50,51]. The quadratic fit relates the geometrical factors and the diameter of the nanorods for spherically capped cylindrical nanorods with different radii. The detailed steps are presented in Appendix A.

2.7. Numerical implementation

A mesh convergence study was carried out to determine the optimum number of elements that yield numerical solutions that are independent of the mesh size. The steps carried out in the mesh convergence study are presented in Appendix B. The setting that resulted in 113,718 tetrahedral elements was found to produce mesh-independent results. This mesh was imported into ValoMC. Monte Carlo simulations were carried out by launching 10 million photons into the computational domain to obtain the fluence distribution across the entire model. Once the fluence distribution has been computed, it was imported into the commercial finite element software COMSOL Multiphysics, where heat transfer simulations were carried out. Boundary layer elements were created across the inner walls of the bladder to facilitate convergence of the fluid flow simulations. First-order approximations were used for both heat transfer and fluid flow simulations, while second-order approximations were employed when solving the thermal damage model.

3. Results

A model verification study was carried out to determine the accuracy of the developed numerical model. Details and the results from the model verification study are presented in Appendix C.

3.1. Optical coefficients

Figure 2 plots the absorption and scattering coefficients against the irradiation wavelength for a monodisperse GNR distribution with diameter and aspect ratio of 10 nm and 3.8, respectively, and volume fractions of 0.001, 0.005, and 0.01%. Peak absorbance occurred at 778 nm wavelength, which matched well the Discrete Dipole Approximation simulations and the experimental results recently reported [52]. Increasing the volume fraction of the GNR led to larger absorption and scattering coefficients; however, the wavelength at peak absorbance was not affected. The absorption and scattering coefficients corresponding to peak absorbance for the different GNR volume fractions investigated are presented in Table 2. It is worth noting that the absorption coefficient of the GNR is a few orders of magnitude higher than the absorption coefficients of the bladder and surrounding tissue.

Figure 2. Plots of the (a) absorption and (b) scattering coefficients at different wavelengths calculated for ϕ = 0.001, 0.005 and 0.01%.

3.2. Fluence distribution

Figure 3 shows the fluence distribution across the z = 0 plane normalized against the laser power for GNR volume fractions of 0.001, 0.005, and 0.01%. The contours are presented in the logarithmic scale for better visualization. The incident irradiation is indicated by the red arrows. The strongly absorbing tumor domain due to the presence of GNR shows very low fluence distribution (see black arrows), which represents one of the distinct features of the contours shown in Figure 3. As GNR volume fraction increases, the absorption coefficient also increases (see Figure 2), which results in greater photon absorption and an even lower fluence distribution.

Figure 3. Fluence distribution across the z = 0 plane. Simulation was performed for three GNR volume fractions (a) $\varphi$ = 0.001%, (b) $\varphi$ = 0.005% and (c) $\varphi$ = 0.01%. The black arrows, indicating low fluence, point to the position of the tumor. An increase in volume fraction leads to lower fluence magnitude due to higher photon absorption.

3.3. Effects of natural convection

Figure 4 illustrates the contours of the temperature distribution across the z = 0 plane and the flow pattern inside the air domain obtained 10 min after the start laser irradiation. The temperature distribution inside the bladder is visualized using isothermal surfaces in order to show the flow-induced variation. The red arrows indicate the air flow in response to heating of the skin surface and the different orientations of the bladder. Regions of high temperatures were observed around the tumor and in the surrounding tissue immediately above the tumor. The maximum and minimum temperatures inside the tumor and the maximum temperature on the skin surface are presented in Table 3. Peak temperature was found at the surface of the skin, which may be explained by the shortest distance to the radiation source; hence the least attenuation.

Figure 4. Contours of the temperature distribution at x = 0 and z = 0 planes 10 min after laser exposure. (a) Case I, (b) Case II, (c) Case III, and (d) Case IV. The temperature distribution inside the bladder is visualized using isothermal surfaces in order to show the flow-induced variation. Results were obtained for $\varphi$ = 0.001% and $P_{\mathrm{laser}}$ = 0.2 W. The yellow arrows show the direction of gravity and the red arrows show the flow field. Figures indicate significant changes in the isotherms inside the bladder due to natural convection under different bladder orientation.

Table 3. Results for the different configurations for $\varphi$ = 0.001% and P_laser = 0.2 W .

Parameter	Maximum tumor temperature (°C)	Minimum tumor temperature (°C)	Maximum skin temperature (°C)	Average convective heat transfer (W/m^2)	Heterogeneity coefficient, $HC$
Case I	55.20	52.27	56.83	3840	0.527
Case II	49.95	46.18	51.14	34,882	0.339
Case III	52.02	48.04	54.25	29,490	0.465
Case IV	55.70	52.76	57.26	—	0.535

The different air flow pattern due to the different bladder orientation was found to affect the degree of convective heat transfer from skin surface to the ambient. From Table 3, in can be seen that Case II produced the smallest maximum tumor and skin surface temperature. As shown in Figure 4(b), gravity in Case II acts in the +y direction. As the skin temperature rises due to laser irradiation, the temperature of air next to the skin surface also increases. Boussinesq approximation states that warm fluid rises due to its lower density, while cold and denser fluid descends. Hence, warm air flows in opposite direction to gravity (−y direction) and impinges on the surface of the skin. The air stream gets diverted across the xz plane, while air at ambient temperature is drawn to the skin surface. This was not observed in Case I because air rises away from the skin surface (see Figure 4(a)), while in Case II, the solid surface prevented cool air to be drawn across the skin surface (see Figure 4(c)).

The effects of natural convection of the urine are better understood if one compares the average temperature gradients inside the bladder. A lower temperature gradient implies a more homogeneous temperature distribution inside the bladder, which from the point of view of GNR-assisted PTT, may be disadvantageous. A homogeneous temperature distribution indicates that heat, which is supposed to be concentrated around the tumor region, is dissipated into other regions of the bladder. The average temperature gradient magnitude inside the bladder in Case I was 18.07 °C/cm, which was not significantly different from Case IV (18.55 °C/cm). Case II had the lowest temperature gradient magnitude of 10.59 °C/cm, while Case III, it was 13.06 °C/cm. These results suggest that the ability for natural convection to produce a more homogeneous temperature distribution inside the bladder during GNR-assisted PTT is dependent on the orientation of the bladder. The heating heterogeneity is further evaluated by quantifying the heterogeneity coefficient $\mathrm{HC},$ which is defined as [53]: (20) $$HC=\frac{T_{10}- T_{90}}{T_{90}- T_{\mathrm{core}}}$$(20) where 10 or 90% of all calculated temperatures at the mesh points are greater or equal to $T_{10}$ and $T_{90},$ respectively, and $T_{\mathrm{core}}$ is the initial temperature of 37 °C. A larger value of $\mathrm{HC}$ implies a more heterogenous system. Values of $\mathrm{HC}$ obtained from Equation (20)(20) $$HC=\frac{T_{10}- T_{90}}{T_{90}- T_{\mathrm{core}}}$$(20) are summarized in Table 3. The $\mathrm{HC}$ value of Case I is approximately 1.52 and 1.1 times that of Cases II and III, respectively, which supports the analysis above.

Figure 5 illustrates the iso-surfaces of the velocity magnitude (for z < 0) and the velocity vector of the thermally driven urine flow inside the bladder at t = 10, 20, 30, 60, and 600 s. During treatment, the urine that is next to the tumor is warmer than in the other regions due to absorption of heat inside the tumor. In Case I, the direction of gravity in -y meant that the warmer fluid that is already at the top of the bladder can only move sideways and downwards, which created two vortices that circulate at either side of the bladder. In Case II, gravity acts in the +y-direction. During the initial stages of heating (see t = 10 s in Figure 5(b)), the warmer fluid next to the tumor moves in the direction opposite to gravity. As the fluid reaches the opposite wall of the bladder, the flow is diverted sideways and upwards to create two vortices (see t = 20 and 30 s). The asymmetrical shape of the bladder between the x > 0 and x < 0 region led to vortex in the x > 0 region to be larger than that in the x < 0 region. As heating continues, the weaker vortex is integrated into the stronger vortex that ultimately resulted in one large vortex in the counter-clockwise direction around the z-axis. In Case III, the direction of gravity along the z-axis meant that the warmer fluid flows in the +z-direction to generate a clockwise vortex around the x-axis.

Figure 5. Contours of the velocity magnitude in mm/s and the velocity vector inside the bladder plotted for z < 0 at t = 10, 20, 30, 60, and 600 s. (a) Case I, (b) Case II, and (c) Case III.

The velocity magnitude in Case I was the lowest, which explains the high urine domain temperature (see Figure 4(a)) and the large values of $\mathrm{HC}$ (see Table 3). In Cases II and III, the maximum velocity magnitude was approximately 5 times larger than in Case I, which explains the lower and near homogeneous temperature inside the bladder for these two cases. The increase in the velocity magnitude found in Cases II and III reflect an increase in convective heat transfer compared to Case I. This can be seen in Table 3 when one compares the average convective heat flux inside the bladder between Cases I, II, and III. As expected, the average convective heat flux in Case II was the highest, followed by Cases III and I.

Contours depicting the probability of thermal damage, PD, across the z = 0 and x = 0 planes, and across the inner surface of the bladder 10 min after laser irradiation are shown in Figure 6. Cooling due to natural convection inside the bladder led to smaller values of PD inside the tumor. Case I, which has the weakest cooling effect, produced the highest probability of observing thermal damage, with the highest value of 77.8%. In Cases II and III, the highest values of PD were 25 and 36.7%, respectively. Additionally, there was a shift in the distribution of PD in these two cases in the direction of the vortex flow, as indicated by the black arrows in Figure 6. In all the cases considered, the targeted PD > 99% thermal damage across the entire tumor domain was not attained.

Figure 6. Contours of PD across the tumor, the bladder and the surrounding tissue at the z = 0 plane, the skin surface (for z < 0) and the inner wall of the bladder, 10 min after the start of laser irradiation. (a) Case I, (b) Case II, (c) Case III, and (d) Case IV, for ϕ = 0.001% and laser power of 0.2 W. Inset: enlarged view of the tumor and the tissue above it.

Figure 7(a) plots the transient variation of the maximum PD values inside the tumor and the bladder. The higher PD in Case I due to the weaker convective heat transfer is clearly demonstrated. Similarly, the highest average convective heat flux (see Table 3) produced the lowest PD values. This suggests different orientations of the bladder can produce different magnitudes of convective heat transfer; some that are sufficiently strong to limit the formation of thermal damage during GNR-assisted PTT of bladder cancer in mice.

Figure 7. Plots of maximum PD values against time inside the tumor and the bladder for (a) $\varphi$ = 0.001%, (b) $\varphi$ = 0.005% and (c) $\varphi$ = 0.01% for $P_{\mathrm{laser}}$ = 0.2 W.

3.4. Effects of GNR volume fraction

The failure to obtain PD > 99% across the entire tumor may be overcome by increasing the amount of heat generated inside the tumor. One way this may be achieved is by increasing the volume fraction of GNR embedded inside the tumor. Simulations were repeated for $\varphi$ = 0.005 and 0.01%, with the values of all other parameters maintained to be the same as those used in Section 3.3. As shown in Figure 2, increasing the value of $\varphi$ resulted in significant elevations in the absorption and scattering coefficients, with the former contributing to larger heat generation during laser irradiation. The temperature distribution and the flow profile inside the bladder and in the air domain (results not presented here) were found to be very similar to those presented in Figures 4 and 5, with the exception of higher temperature and larger velocity magnitude as a result of the increase in the heat generated inside the tissue.

Contours of PD obtained for $\varphi$ = 0.005% are presented in Figure 8(a). As expected, increasing the value of $\varphi$ led to greater thermal damage inside the tumor. At $\varphi$ = 0.005%, the highest PD approached 99.9% in Case I, while the highest PD for Cases II and III were 59.0 and 80.5%, respectively, which were more than doubled compared to the case with $\varphi$ = 0.001%. Although there was a clear elevation in the values of PD, the targeted PD > 99% across the entire tumor domain was not achieved. In Case I, the threshold of PD > 99% was observed only at the top half of the tumor, while the bottom half remained below 80%.

Figure 8. Contours of PD across the tumor, the bladder and the surrounding tissue across the z = 0 plane, the skin surface (for z < 0) and the inner wall of the bladder after 10 min of 0.2 W laser irradiation. (a) ϕ = 0.005% and (b) ϕ = 0.01%. Inset: enlarged view of the tumor and the tissue above it.

The distribution of PD obtained for $\varphi$ = 0.01% are presented in Figure 8(b), where complete tumor damage was found for Case I, as indicated by PD > 99% across the entire tumor domain. Although Case III showed values of PD that approached 99%, these were restricted to only the top half of the tumor. The maximum PD value in Case II was the lowest at 97.5% and was concentrated only at the top surface of the tumor. It is noteworthy that the increase in tumor thermal damage due to the increase in $\varphi$ was also accompanied by increases in the damage to the bladder and the surrounding tissue, especially for Case I. This is indicated by the white arrows in Figure 8. The role of natural convection inside the bladder in the dispersion of heat from the ablated region is further demonstrated by the low PD values in the surrounding tissue for Cases II and III.

Figures 7(b,c) plot the transient variation in maximum PD for the tumor and bladder for Cases I–III obtained for $\varphi$ = 0.005% and 0.01%, respectively. As heating increases due to larger $\varphi$ values, PD also increases in both the tumor and the bladder. Unlike the case when $\varphi$ = 0.001% (see Figure 7(a)), the values of maximum PD in the bladder were not as high as in the tumor. In this case, cooling due to convective heat transfer may have been insufficient to transfer the intense heat, due to the higher GNR concentrations, from the targeted region to the adjacent bladder tissue.

3.5. Effects of laser power

In addition to GNR volume fraction, the heat generated inside the tissue may be elevated by increasing the laser power. To investigate the effects of laser power on thermal damage sustained by tumor, $P_{\mathrm{laser}}$ was increased to 0.3 W and simulations were repeated for $\varphi$ = 0.001, 0.005, and 0.01%. The temperature distribution and velocity profiles obtained were similar to those presented in Figures 5 and 6, apart for their higher magnitudes. Hence, these results are not presented here.

Contours of PD are shown in Figure 9 for $P_{\mathrm{laser}}$ = 0.3 W. The results show significantly larger thermal damage inside the tumor in all cases. The threshold of PD > 99% across the entire tumor was achieved for Case I even at $\varphi$ = 0.001%, while Cases II and III failed to achieve the target. The target was attained in Cases I and III at $\varphi$ = 0.005, while Case II was not far off. At $\varphi$ = 0.01%, it would appear that all cases achieved the target of PD > 99% across the entire tumor domain. The combination of high laser power and GNR volume fraction produced very intense heating, such that significant portion of the bladder in Case I experienced thermal damage (PD > 99%) due to the weak convection inside the bladder. Similar effects were observed in Case III with $\varphi$ = 0.005 and 0.01%. On the other hand, such damage to the bladder was not found in Case II likely due to the strong convective heat transfer inside the bladder. It is important to note that the attainment of PD > 99% across the entire tumor may not be clinically acceptable if the nearby bladder tissue is also compromised.

Figure 9. Contours of PD across the tumor, the bladder and the surrounding tissue at the z = 0 plane, the skin surface (for z < 0) and the inner wall of the bladder after 10 min of 0.3 W laser irradiation. (a) ϕ = 0.001%, (b) ϕ = 0.005% and (c) ϕ = 0.01. Inset: enlarged view of the tumor and the tissue above it.

Significant amount of thermal damage was also found in the surrounding tissue above the tumor and on the skin surface in Case I for all the values of $\varphi$ and in Case III for $\varphi$ = 0.01%. However, such damage to the surrounding tissue and the skin surface was not observed in Case II regardless of the GNR volume fraction. This may be explained not only by the strong convective heat transfer inside the bladder, but also by the larger cooling induced by the air flow on the skin side, as explained in Section 3.3.

3.6. Sustained thermal damage inside the tumor, bladder, and surrounding tissues

The contours of PD presented in the preceding sections describe the probability of observing thermal damage and their visual examinations provide a qualitative indication on the distribution of thermal insult across the tissue. In order to quantify the amount of sustained thermal damage in each tissue region, the following parameters were introduced: $R_{\mathrm{T}},$ $R_{\mathrm{B}},$ and $R_{\mathrm{ST}},$ which represent the ratio of volume of thermal damage sustained by the tumor, bladder, and surrounding tissues relative to their total volume. These parameters are defined as: (21) $$R_{\mathrm{T}}=\frac{V_{\mathrm{TD}}}{V_{\mathrm{T}}},$$(21) (22) $$R_{\mathrm{B}}=\frac{V_{\mathrm{BD}}}{V_{\mathrm{B}}},$$(22) (23) $$R_{\mathrm{ST}}= \frac{V_{\mathrm{STD}}}{V_{\mathrm{ST}}},$$(23) where $V_{\mathrm{TD}},$ $V_{\mathrm{BD},}$ and $V_{\mathrm{STD}}$ are, respectively, the volumes of tumor, the bladder and the surrounding healthy tissues that have been thermally destroyed, i.e., PD > 99% (see Section 2.4), and $V_{\mathrm{T}},$ $V_{\mathrm{B}},$ and $V_{\mathrm{ST}}$ are the total volume of the tumor, bladder and surrounding healthy tissues, respectively. A value of R = 1 implies that the whole tissue is thermally destroyed, while a value of R = 0 implies that there is no sustained thermal damage inside the tissue. An optimum treatment should ideally maximize $R_{\mathrm{T}}$ while minimizing $R_{\mathrm{B}}$ and $R_{\mathrm{ST}}.$

Values of $R_{\mathrm{T}},$ $R_{\mathrm{B}}$ and $R_{\mathrm{ST}}$ calculated using Equations (21)–(23)(1) $${\rho }_{\mathrm{t}}{c}_{\mathrm{t}}\frac{\partial T_{\mathrm{t}}}{\partial t}=\nabla ·\left(k_{\mathrm{t}}\nabla T_{\mathrm{t}}\right)+{\rho }_{\mathrm{b}}{c}_{\mathrm{b}}{\omega }_{\mathrm{b}}\left(T_{\mathrm{b}}-T\right)+Q_{\mathrm{m}}+{\mu }_{a,t}{\mathrm{\Phi }}_t{P}_{\mathrm{laser}}$$(1) are presented in Table 4. When using a laser power of 0.2 W, complete ablation of the tumor (R_T = 1) can be attained only for Case I for $\varphi$ = 0.01%, with a 0.1% damage to the bladder (R_B = 0.001). When using a laser power of 0.3 W, complete tumor ablation was attained at all GNR volume fractions for Case I and at $\varphi$ = 0.01% for Case III. However, these outcomes may not be acceptable due to the increase in the amount of damage sustained by the bladder and the surrounding tissue.

Table 4. Values of $R_{\mathrm{T}},$ $R_{\mathrm{B}},$ and $R_{\mathrm{ST}}$ obtained for Cases I, II, and III for different laser power and GNR volume fraction.

Parameters	Case I			Case II			Case III
	$R_{\mathrm{T}}$	$R_{\mathrm{B}}$	$R_{\mathrm{ST}}$	$R_{\mathrm{T}}$	$R_{\mathrm{B}}$	$R_{\mathrm{ST}}$	$R_{\mathrm{T}}$	$R_{\mathrm{B}}$	$R_{\mathrm{ST}}$
Plaser = 0.2 W
$\varphi =0.001\%$	0	0	0	0	0	0	0	0	0
$\varphi =0.005\%$	0.125	0	0	0	0	0	0	0	0
$\varphi =0.01\%$	1^a	0.001	0	0	0	0	0.066	0	0
Plaser = 0.3 W
$\varphi =0.001\%$	1	0.005	0	0	0	0	0.002	0	0
$\varphi =0.005\%$	1	0.049	0.004	0.389	0	0	0.783	0.001	0
$\varphi =0.01\%$	1	0.112	0.012	0.97	0.005	0	1^a	0.017	0.001

^aValues greater than 0.998 and thus rounded to 1.

For Case III, only the combination of $P_{\mathrm{laser}}$ = 0.3 W and $\varphi$ = 0.01% produced complete tumor ablation None of the combination of laser power and GNR volume fraction led to complete tumor ablation for Cases II; clearly indicating the significance of orientation on the cooling induced by natural convection inside the bladder.

4. Discussions

Natural convection inside the bladder during GNR-assisted PTT of bladder cancer in mice was found to produce cooling effect that can affect the outcome of the treatment. The degree of cooling was found to depend on the orientation of the bladder with respect to the laser beam and gravitational field. Placing the mouse on its back, as shown for Case I in Figure 1(c), minimizes the cooling effect, while having the mouse on its abdomen, as shown for Case II in Figure 1(c), maximizes it. Initially, there were concerns on whether the localized heating of GNR-assisted PTT is sufficiently powerful to produce any meaningful urine circulation inside the bladder. The results obtained indicated that for the case of mice, the effects of natural convection cannot be ignored. This suggests that future experiment in mice of GNR-assisted PTT for bladder cancer treatment of mice must take into consideration the cooling effects due to natural convection inside the bladder.

Depending on the orientation, cooling induced by natural convection inside the bladder can be both detrimental and beneficial to the outcome of bladder cancer treatment. Since heating during GNR-assisted PTT is localized to the targeted tissue, natural convection can be detrimental by distributing heat to other parts of the bladder. This not only causes insufficient heating to the targeted tumor tissue, as can be seen in Figure 4 for Cases II and III, but also raises the temperature of the bladder and the surrounding tissue to levels that may be sufficiently high to induce some degree of thermal damage (see Figure 6(b)). On the other hand, in cases where heating is intensive, either due to the selection of large GNR volume fraction or the use of high-powered laser, natural convection can help to reduce the temperature rise in bladder and surrounding tissues in order to minimize the thermal damage to healthy tissue. Achieving a balance between the detrimental and beneficial effects of natural convection during GNR-assisted PTT of bladder cancer is a matter that is worthy of future investigations.

Different orientation of the bladder also affects the degree of convection heat transfer between the air and the skin surface. It was found that the flow pattern of air generated when the mouse is oriented as in Case II, i.e., when laying on its abdomen, resembles that of an impinging jet, which induces stronger cooling effect than those generated by the air flow in Cases I and III. The ability to cool the skin surface is seen as beneficial as it helps to prevent the risk of over-heating that may arise from the use of larger GNR volume fraction or higher laser power. Unfortunately, this orientation also led to strong convective currents inside the bladder that render the treatment ineffective due to the failure in attaining complete tumor ablation.

The findings obtained from the present study may be used to identify an optimal combination of parameters and treatment orientation for the GNR-assisted PTT of bladder cancer in mice. Based on the data in Table 4, the combination of $\varphi$ = 0.01% and $P_{\mathrm{laser}}$ = 0.2 W with the mouse laying on its back (Case I) is deemed to be optimum as it ensures complete tumor ablation, while minimizing damage to the non-cancerous tissues. In a clinical setting, the orientation pertaining to Case I is akin to the patient in the more natural state of laying down. It is noteworthy that this study does not consider the feasibility of GNR volume fraction (cost, safety etc) in real applications. Should the situation prevent the use of 0.01% GNR volume fraction, then the combination of $\varphi$ = 0.001% and $P_{\mathrm{laser}}$ = 0.3 W a with Case I orientation can be used, though with a slightly higher to the bladder.

It is important to note that the conclusions drawn in the present study were based on investigations in mice. As such, caution must be exercised when interpreting and translating the results from the present study to human patients. The differences in both the size and the shape of human and mouse bladders mean that the effect of natural convection in GNR-assisted PTT of bladder cancer may not scale linearly from mice to humans. If the laser beam diameter remains the same, cooling due to natural convection inside larger bladders, i.e., human bladders, may be more localized and not as significant as in mice bladders. Nevertheless, computational studies remain an integral part of the development of GNR-assisted PTT, particularly at the present stage, where the majority of the experimental studies were carried out in mice [10,12,13].

The two methods considered in this study for overcoming the incomplete tumor destruction due to cooling from natural convection were by increasing the volume fraction of GNR inside the tumor and by increasing the laser power. The former increases the amount of heat generated inside the tumor by increasing the absorption coefficient, while the latter employs a higher intensity laser beam into the tissue. The thermal responses observed were very different for the two methods. Increasing the GNR volume fraction increases the absorption coefficient inside the tumor. Therefore, the tumor experiences higher heat generation and temperature elevation during laser irradiation. The increase in temperature and thermal damage to the bladder and the surrounding tissue were due to conduction of heat from the tumor. Increasing the laser power results in higher energy absorption and heat generation inside the tissue. Since the laser beam is introduced as an external source, the photons will have to travel through the surrounding tissue before arriving at the tumor. The absorption of the photons by the surrounding tissue results in significant elevation in temperature that could lead to severe thermal damage where it is not intended. These results suggest that it may be more acceptable to increase the GNR volume fraction than to increase the laser power when attempting to increase the amount of thermal damage inside the tumor.

The present study identifies some limitations that need to be addressed in future investigations. Firstly, the tumor was assumed to have the same thickness as the bladder, which represents bladder cancer that is at an advanced stage. These cases may be better treated through transurethral resection of bladder tumor [54] or cystectomy [55], although both options are invasive. Secondly, the present study has assumed a spatially homogeneous distribution of GNR. The distribution of GNR is likely to be heterogeneous and the actual distribution will depend on the method employed for introducing the GNR into the tumor. Thirdly, the dimensions of the GNR for a given population are likely to vary over a range instead of being uniform in size. In such cases, one may have to assume a polydisperse distribution of GNR when evaluating the absorption and scattering coefficients using the Mie-Gans theory. Nevertheless, the effects of a polydisperse GNR distribution is expected to affect only the absorption and scattering coefficients, whose effects were investigated in the current study by varying the GNR volume fraction. As such, the use of polydisperse GNR distribution may lead to outcomes that are similar to those found in this study. Finally, conclusions in the present study are based on investigations in mice. While the results obtained are informative, they may not translate well to humans due to differences in anatomy.

5. Conclusions

The present study has demonstrated the effects of natural convection inside the bladder during GNR-assisted PTT of bladder cancer in mice. Although heating is localized, the thermal gradient produced is sufficiently high to induce convective flow inside the bladder causing additional cooling during treatment. The cooling effect can be both favorable and unfavorable, depending on the required level of heating inside the tissue during treatment. The degree of cooling was also found to depend on the orientation of the bladder. Different treatment orientations were found to affect the heat transfer between the skin surface and the surrounding air. The orientation corresponding to the mouse laying on its abdomen induced the strongest cooling between the skin surface and the air, which prevented the formation of thermal damage on the skin, especially when higher laser power was employed. Investigations were also carried out to determine the role of GNR volume fraction and laser power on the level of thermal damage inside the tumor. Increasing the laser power was found to be less favorable due to the extensive thermal damage to the tissue through which the laser beam travels before reaching the tumor.

Disclosure statement

No potential conflict of interest was reported by the author(s).

.

Additional information

Funding

This study has received funding from the European Union’s Horizon 2020 research and innovation program under grant agreement No 801126 (EDIT).

Appendix A

Geometrical factors for spherically capped nanorods

The Mie–Gans theory for calculating the absorption and scattering coefficients assumes the GNR to be spheroids. As such, the geometrical factor $P_i$ (for $i$ = 1, 2, 3) in Equation (15) can be calculated using: (A1) $$P_1=\frac{1-e^2}{e^2}\left[\frac{1}{2e}ln\left(\frac{1+e}{1-e}\right)-1\right],$$(A1) (A2) $$P_2=P_3=\frac{1}{2}(1-P_1),$$(A2) where $e={\left[1-{\left(D/L\right)}^2\right]}^{0.5},$ $P_1$ is oriented in the longitudinal axis of the spheroid, while $P_2$ and $P_3$ are in the transverse directions, with $P_2 = P_3$ due to symmetry. Experimental studies characterizing the morphology of GNR have shown them to resemble cylinders that are capped at the ends by hemispheres, which is far from the spheroidal shape assumed by the Mie-Gans theory. According to Prescott and Mulvaney [50,51], the difference between the actual shape of the GNR and the assumed shape in the Mie–Gans theory can result in the failure to obtain results that match the data from experimental studies.

To overcome this, the same authors derived alternative values of P_i for non-spheroidal nanoparticles that were expressed as quadratic functions of the nanoparticle aspect ratio, given by: (A3) $$\frac{1-P_1}{P_1}=a{\left(\mathrm{AR}\right)}^2+b\left(\mathrm{AR}\right)+c,$$(A3) where $\mathrm{AR}$ is the aspect ratio of the nanoparticle and $a,$ $b,$ and $c$ are the coefficients of the quadratic function. For a nanorod with morphology that is described by a cylinder capped at the ends by hemispheres, the coefficients of the quadratic function for a nanorod of diameter 10 nm, which is the size considered in this study, are $a$ = 0.3504, $b$ = 2.0809, and $c$ = −0.3035 [50,51]. It is important to point out that Equation (A3)(3) $${\rho }_{\mathrm{f}}\frac{\partial \mathit{u}}{\partial t}+{\rho }_{\mathrm{f}}\left(\mathit{u}·\nabla \mathit{u}\right)=-\nabla p+{\mu }_{\mathrm{f}}{\nabla }^2\mathbf{u}+{\rho }_{\mathrm{f}}\left(T_{\mathrm{f}}\right)\mathbf{g},$$(3) is purely empirical and no physical meaning should be inferred from it.

Appendix B

Mesh convergence study

Mesh convergence was examined by systematically decreasing the maximum element size of each domain, starting from the tumor domain, followed by the bladder, the urine, and the surrounding tissue. Convergence was examined for the velocity magnitude and temperature at several points sampled across the solution domain. The solutions were deemed to be mesh-independent when differences in the solutions across the sampled points between two successive mesh settings were less than 5% for velocity magnitude and 1% for temperature. Although the criterion for velocity was less strict, the choice of 5% was deemed to be acceptable due to the nonlinear nature of the Navier–Stokes equations.

Results from the mesh convergence study are shown in Figure B1, which plots the percentage difference of the velocity magnitude and temperature sampled at various points between two successive mesh settings. The mesh settings employed led to total number of elements ranging from 14,857 to 113,718 elements. The results in Figure B1 showed that the model with 85,501 elements was sufficiently fine to achieve mesh convergence. Given that the CPU times required for solving the models with 85,501 elements or higher were not significantly different from one another, the mesh employing 113,718 elements was selected for the current study.

Figure B1. Results from mesh convergence study, where the different bars represent the values obtained at different sampled points.

Appendix C

Model verification study

The in vitro experimental study of Terentyuk et al. [56] was adopted for the purpose of the model verification study. In their study, 1.5 ml solutions containing GNR at concentrations of 8 and 100 µg/ml were placed inside Eppendorf tubes and were irradiated with a laser beam of 0.36 W (laser radiation of 1.2 W/cm² and laser spot area of 0.3 cm²). The GNR has mean diameter and length of 10.2 and 41 nm, respectively.

A 3D model of a standard Eppendorf tube containing 1.5 ml of GNR solution was developed. The governing equations were modified in order to be in agreement with the in vitro conditions of the experiments; in this case, the absence of metabolic heat generation and blood perfusion. A Robin boundary condition similar to Equation (8)(8) $$\{\begin{array}{c}T=T_{\mathrm{amb}},\mathit{\text{ if }}\mathit{u}·\mathit{n}<0\\ -k_{\mathrm{f}}\frac{\partial T_{\mathrm{f}}}{\partial n}=0,\mathit{\text{ if }}\mathit{u}·\mathit{n}\ge 0\end{array},$$(8) was prescribed along the outer surface of the model. The ambient temperature and the initial temperature of the GNR solution inside the Eppendorf tube were assumed to be 27.5 °C. Since no information was available on the ambient convection coefficient h_amb for implementation of the Robin condition, the verification exercise was carried out by first varying h_amb to match the numerical results with the experimental results for the solution with 8 µg/ml GNR. The value of h_amb that led to good match (less than 5% difference) was then used to verify the model with solutions of 100 µg/ml GNR.

compares the maximum temperature obtained with the computational model and the experimental results from the first to the fifth minute, for the case of 8 µg/ml concentration of GNR. One may observe that a value of h_amb of 5 W/m²K yielded percentage difference of less than 10% between the numerical model and the experimental results. When h_amb of 5 W/m²K was employed, the percentage difference ranged from 2.1 to 15% for the case of GNR concentration of 100 µg/ml, see Figure C1. This range is considered to be acceptable for the model verification effort in this study given the uncertainty involved in some of the parameters.

Figure C1. Results from the model verification study for GNR concentration of 100 µg/ml.

Table C1. Results from the model verification study against the work of Terentyuk et al. [45].

Time (min)	Maximum temperature (^oC)
	$h_{\mathrm{amb}}$ (W/m^2K)			Terentyuk et al. [45]
	15	10	5
1	35.94 (1.2)	36.21 (0.4)	36.51 (0.4)	36.36
2	40.83 (6.4)	41.56 (4.8)	42.34 (3.0)	43.64
3	44.28 (9.8)	45.63 (7.0)	47.02 (4.2)	49.09
4	45.15 (16.2)	48.90 (9.3)	50.98 (5.4)	53.86
5	45.90 (20.9)	51.59 (11.1)	54.42 (6.2)	58.00

Numbers in brackets indicate the percentage difference between the numerical and experimental results.