# NUMERICAL HYPERTHERMIA SIMULATION FOR A 3-D TRIPLE-LAYERED SKIN STRUCTURE WITH EMBEDDED VASCULAR COUNTERCURRENT NETWORK AND NANOPARTICLES

NUMERICAL HYPERTHERMIA SIMULATION FOR A 3-D TRIPLE-LAYERED SKIN STRUCTURE WITH EMBEDDED VASCULAR COUNTERCURRENT NETWORK AND NANOPARTICLES

Casey Orndorff Weizhong Dai*

Mathematics and Statistics, College of Engineering and Science, Louisiana Tech University, Ruston, Louisiana, 71272, USA

Corresponding Author Email:
dai@coes.latech.edu
Page:
S179-S184
|
DOI:
https://dx.doi.org/10.18280/ijht.34S124
|
Accepted:
|
Published:
31 January 2016
| Citation

OPEN ACCESS

Abstract:

In recent years, there has been interest in research relating to hyperthermia skin cancer treatment.  However, obtaining accurate temperature distributions in living tissue without using an intruding sensor is a challenge.  The purpose of this study is to develop a mathematical model that can accurately determine the temperature distribution in the tumor region and surrounding normal tissue.  Our model is based on a modified Pennes’ equation for the bioheat transfer in a 3-D triple-layered skin structure embedded with a vascular countercurrent network and a tumor appearing in the subcutaneous region.  The vascular network is designed based on the constructal theory of multi-scale tree-shaped heat exchangers.  The tumor is injected with golden nanoparticles in order to heat up the tumor quickly.  The proposed model will be solved numerically using an unconditionally stable finite-difference scheme.  We expect our method can be a useful tool to provide the optimal laser irradiation in order to eradicate the tumor through hyperthermia treatment while keeping the damage to the surrounding healthy tissue to a minimum.

Keywords:

Constructal Law, Skin Living Tissue, Finite-Difference Method, Hyperthermia, Bioheat Transfer.

1. Introduction

In the recent years, there has been interest in research related to the hyperthermia combined with other treatments, such as radiation and cytotoxic drugs, to enhance the killing of tumors. Conventional hyperthermia (target temperatures of 42-46°C) in conjunction with radiation has demonstrated increased effectiveness in the treatment of certain types of cancer, such as those of liver metastases.  The objective is the control laser heating of the tumor so that the temperatures of the normal tissue surrounding the tumor remains low enough so as not to cause damage to the tissue.  Thus, the entire temperature field in the treatment region, clinical personnel can potentially control the heating source to deliver energy to the treatment target volume to raise its minimum temperature above 42°C, while limiting the temperatures in the normal tissue to prevent damage.  However, it is not easy to obtain an accurate determination of the temperature field over the entire treatment region during clinical hyperthermia treatments, because the number of invasive temperature probes that can be used is limited due to the pain tolerances of the patients.  Hence, it is important to determine the laser intensity and pattern of laser exposure in order to optimize the temperature distribution in the treated region.  The determinants of temperature distributions during thermal therapy are: the power deposition pattern of the heating source, heat removal by conduction, and heat removal by blood flow forced convection.  In order to determine the  temperature  distribution  in  the  treated  region, numerical methods must be developed to solve the bioheat transfer equation for the human body .

Dai et al.  have developed a numerical method for optimizing laser power in the irradiation of a 3-D triple-layered cylindrical skin structure. The method was obtained by solving numerically the 3-D bioheat equation where the surface of the skin is irradiated by the laser according to a designed pattern. Tang and Dai et al.  developed a model for an optimal temperature distribution in a 3-D triple-layered skin structure embedded with artery and vein vasculature that employs the 3-D Pennes' equation to determine heat distribution and relies on the Crank-Nicolson scheme to solve numerically. Dai et al.  also developed a numerical algorithm that takes into account both the laser intensity and laser exposure pattern to obtain an optimal temperature distribution in a 3-D triple layered cylindrical skin structure based on the Pennes' equation; later they approached the same problem with a blood vessel embedded in the model to further the progress of an optimal temperature distribution in a real system . Dai and Bejan et al.  developed a model for a 3-D triple-layered skin structure with a vascular countercurrent network that employed a modified Pennes' equation that accounted for thermal lag of the tissue. Vera and Bayazitoglu  developed a model that used finite difference time domain (FDTD), an explicit method, to calculate the heat distribution of nanoparticles on a single layer of human biological tissue at varying particle distribution densities in different host mediums.

The purpose of this study is to develop a mathematical model that can accurately determine the temperature in the tumor region and surrounding normal tissue.  Our model is based on a modified Pennes’ equation for the bioheat transfer in a 3-D triple-layered skin structure embedded with a vascular countercurrent network and a tumor appearing in the subcutaneous region.  The vascular network is designed based on the constructal theory of multi-scale tree-shaped heat exchangers.  The tumor is injected with golden nanoparticles in order to heat up the tumor quickly.  The proposed model will be solved numerically using an unconditionally stable finite-difference scheme.  We expect our method can be useful tool to provide the optimal laser irradiation in order to eradicate the tumor through hyperthermia treatment while keeping the damage to the surrounding healthy tissue to a minimum.

2. Mathematical Model

The model under consideration is comprised of a modified Pennes’ equation, a vascular network based on the constructal theory of multi-scale tree-shaped heat exchangers, and a heat source that is comprised of a laser irradiating the area of consideration along with nanoparticle heating in the tumor region as shown in Figure 1.

## 1.png Figure 1. A three-dimensional triple-layered skin structure with embedded countercurrent vasculature network with a tumor in the subcutaneous region.

2.1. Modified Pennes’ equation

The modified Pennes’ equation that governs bioheat transport in a 3-D triple-layered skin tissue and tumor region can be written as follows , , [7-9].

\rho_{l} C_{l} \frac{\partial T_{l}}{\partial t}+W_{b}^{l} C_{b}^{l}\left(T_{l}-T_{o u t}\right)=k_{l}\left(\frac{\partial^{2} T_{l}}{\partial x^{2}}+\frac{\partial^{2} T_{l}}{\partial y^{2}}+\frac{\partial^{2} T_{l}}{\partial z^{2}}\right)+Q_{l}, \quad l=1,2,3, t(1) Here, T_{l} is the temperature of the lth tissue layer or tumor region; T_{o u t} is the blood temperature at the exit or entrance of the 7th level vessel for the artery or vein; \rho_{l}, C_{l} and k_{l} denote the density, specific heat, and thermal conductivity of the lth skin tissue layer, respectively; C_{b}^{1} is the specific heat of blood; W_{b}^{l} is the blood perfusion rate; and Q_{l} is the volumetric heating due to spatial heating , . 2.2. Vascular network To simplify computations, we consider the skin tissue to be a rectangular structure embedded with a seven-level countercurrent vascular network, which is a highly branching and hierarchical network in the subcutaneous layer as described in , as shown in Figure 1. It should be pointed out that only large blood vessels can be seen in the subcutaneous tissue because the dermis is very sparingly supplied with capillaries and the capillary beds of skin lying immediately below the epidermis, and thus, the contribution of these small vessels to the heat transfer can be ignored . In Figure 1, the red color dendritic network represents arteries while the blue color dendritic network represents veins, where all are considered as slim cuboids for simplicity. Levels of arteries are designed such that the first-level artery comes from right to left running lengthwise in the y-coordinate; the second-level artery branches from the left end of the first-level artery running lengthwise in the x-coordinate; the third-level artery has two vessels branching from the two ends of the second-level artery running lengthwise in the z-coordinate; there are four fourth-level arteries branching from the four ends of the third-level arteries running lengthwise in the y-coordinate; the fifth-level artery has eight arteries branching from the eight ends of the fourth-level arteries running lengthwise in the x-direction; there are 16 vessels in the sixth-level artery running length wise in the z-direction; and finally, there are 32 vessels in the seventh-level artery running length wise in the y-direction. The vein network has the same number of blood vessels as its counterpart artery in corresponding levels. In sum, there are 128 blood vessels in total in the considered skin structure . To determine the diameters of the many blood vessels, we follow the constructal theory of multi-scale tree-shaped heat exchangers , [12-14] and assume that the diameters of arteries are decreasing by a constant ratio \gamma between successive levels of branched vessels, given by ,  \gamma=\frac{N L_{b}^{m+1}}{N L_{b}^{m}}=\frac{N W_{b}^{m+1}}{N W_{b}^{m}}=2^{-\frac{1}{3}}, \quad m=1, \cdots, 6(2) where N L_{b}^{m} and N W_{b}^{m} are the length and width of the cross section of a blood vessel in level m, respectively. The length of a blood vessel is assumed to double after two consecutive construction steps, which can be expressed in the length-doubling rule  as follows, L_{b}^{m}=2^{-\frac{1}{2}} L_{b}^{m+1}, \quad m=1, \cdots, 6(3) where L_{b}^{m} is the length of the blood vessel in level m. The mass flow of blood in the mth level vessel, M_{m}=v_{m} F_{m} satisfies  M_{m}=2 M_{m+1}, \quad m=1, \cdots, 6(4) where \mathcal{V}_{m} is the blood flow velocity and F_{m}\left(=N L_{b}^{m} \times N W_{b}^{m}\right) is the area of the cross-section in the mth level vessel . Furthermore, the blood temperature in the cross-section of a vessel is assumed to be uniform. We further assume that a steady-state energy balance in the blood vessel can be reached because the length of the considered blood vessel is relatively short and the blood velocity is relatively high. However, one may use the transient heat transfer equation for a more accurate solution. Hence, the convective energy balance equations, which are used to calculate the artery (levels 1-6) blood temperatures, can be expressed as [2-5], ,  C_{B} M_{1} \frac{\partial T_{b}^{1}}{\partial y}-\alpha P_{1}\left(T_{w}^{1}-T_{b}^{1}\right)=0(5a) C_{B} M_{2} \frac{\partial T_{b}^{2}}{\partial x}-\alpha P_{2}\left(T_{w}^{2}-T_{b}^{2}\right)=0(5b) C_{B} M_{3} \frac{\partial T_{b}^{3}}{\partial z}-\alpha P_{3}\left(T_{w}^{3}-T_{b}^{3}\right)=0(5c) C_{B} M_{4} \frac{\partial T_{b}^{4}}{\partial y}-\alpha P_{4}\left(T_{w}^{4}-T_{b}^{4}\right)=0(5d) C_{B} M_{5} \frac{\partial T_{b}^{5}}{\partial x}-\alpha P_{5}\left(T_{w}^{5}-T_{b}^{5}\right)=0(5e) C_{B} M_{6} \frac{\partial T_{b}^{6}}{\partial x}-\alpha P_{6}\left(T_{w}^{6}-T_{b}^{6}\right)=0(5f) where C_{B} is the heat capacity of blood, a is the heat transfer coefficient between blood and tissue, and P_{m} is the vessel perimeter. In addition, T_{w}^{m} and T_{b}^{m} are the wall temperature and the blood temperature in the mth level vessel. For the smallest, terminal arterial vessels (level 7), a decreased blood flow rate ( \dot{P} ) is included in the energy balance equation [2-5],,. C_{B} M_{7} \frac{\partial T_{b}^{7}}{\partial y}-\alpha P_{7}\left(T_{w}^{7}-T_{b}^{7}\right)-\dot{P} C_{B} F_{7} T_{b}^{7}=0(6) The venous network is assumed to be similar to the arterial network, except that the blood flow direction in each vein is opposite of that in the artery; i.e., countercurrent flow occurs in these two kinds of vessels (see Figure 1). Also, the diameter ratio, length ratio, and mass flow ratio of the blood between the successive levels of the branched veins take the same form, as described in (2)-(4) for the arteries. Moreover, the convective energy balance equations (5a)-(5f) and (6) used to calculate the blood temperature in the artery domain are applied to the vein domain at the corresponding levels . 2.3. Heat source We assume that the laser power, Q_{l} , at layer l is continuous and spatial with a normal distribution , , Q_{1}=\alpha_{1} e^{-\alpha_{1} z} \frac{1}{\sqrt{2 \pi \sigma^{2}}} e^{\frac{\left[x-x_{0}(t)\right]^{2}+\left[y-y_{0}(t)\right]^{2}}{2 \sigma^{2}}} * P_{0}\left(1-\operatorname{Reff}_{1}\right)(7a) Q_{2}=\alpha_{2} e^{-\alpha_{1} L_{1}-\alpha_{2} z} \frac{1}{\sqrt{2 \pi \sigma^{2}}} e^{\frac{\left[x-x_{0}(t)\right]^{2}+\left[y-y_{0}(t)\right]^{2}}{2 \sigma^{2}}} * P_{0}\left(1-\operatorname{Reff}_{2}\right)(7b) Q_{3}=\alpha_{3} e^{-\alpha_{1} L_{1}-\alpha_{2} L_{2}-\alpha_{3} z} \frac{1}{\sqrt{2 \pi \sigma^{2}}} e^{\frac{\left[x-x_{0}(t)\right]^{2}+\left[y-y_{0}(t)\right]^{2}}{2 \sigma^{2}}} * P_{0}\left(1-\operatorname{Reff}_{3}\right)(7c) On the other hand, the heat source for inside the tumor region is defined as \begin{align} & {{Q}_{t}}={{\alpha }_{3}}{{e}^{-{{\alpha }_{1}}{{L}_{1}}-{{\alpha }_{2}}{{L}_{2}}-{{\alpha }_{t}}z}}\frac{1}{\sqrt{2\pi {{\sigma }^{2}}}}{{e}^{\frac{{{[x-{{x}_{0}}(t)]}^{2}}+{{[y-{{y}_{0}}(t)]}^{2}}}{2{{\sigma }^{2}}}}} \\ & \text{ } \\ & \text{ }*{{P}_{0}}(1-\text{Ref}{{\text{f}}_{t}})+{u}'''\text{ } \\ \end{align} (7d) Here, for Eqs. (7a)-(7d), \alpha_{1}, \alpha_{2}, \alpha_{3}, and \alpha_{t} are laser absorptivities of the three layers and tumor region, respectively; Reff1, Reff2, Reff3 and Refft are laser reflectivities of the three layers of the skin and tumor region, respectively; \sigma is the standard deviation of the width of a normally distributed laser beam; and L_{1}, L_{2}, and L_{3} are the lengths of the three layers, respectively. Here, \left[x_{0}(t), y_{0}(t)\right] is the location where the laser is focused at time t. P_{0} is the laser intensity, which will be determined later, so that an optimal temperature distribution can be obtained. Optimality is achieved by minimizing the sum of square deviations between observed and pre-specified temperature elevations at different locations on the skin surface . The u^{\prime \prime \prime} term is related to the heating of nanoparticles and is given by  u^{\prime \prime \prime}=-\frac{d q_{R \lambda}\left(\tau_{\lambda}\right)}{d z}(8) where spectral radioactive heat flux, q_{R \lambda}\left(\tau_{\lambda}\right) , is the sum of the collimated radiation, q_{c \lambda}\left(\tau_{\lambda}\right) , and the heat generation caused by diffuse radiation, q_{d \lambda}\left(\tau_{\lambda}\right) . It should be noted that \tau_{\lambda} is defined as the optical depth which is a function of z. 2.4. Boundary and initial conditions We assume on the laser-facing boundary and the non-laser-facing boundaries, respectively , k_{1} \frac{\partial T_{1}}{\partial z}=h\left(T_{1}-T_{a}\right), \quad z=0(9) here, h is the convective heat transfer coefficient, T_{a} is the ambient temperature. For simplicity, we further assume that the heat flux approaches zero as the tissue depth increases, which is realistic for a biological body . The other boundary conditions in the tissue are assumed to be ,  \frac{\partial T_{l}}{\partial \mathbf{n}}=0(10) where n is the unit outward normal vector on the boundary. At the entrance to the first level vessel, we assume T_{b}^{1}=T_{i n}(11) where T_{i n} is the blood temperature at the entrance of the artery. At the exit of the artery, the blood temperature is equal to the surrounding tissue temperature. T_{b}^{\eta}=T_{o u t}(12) The continuity of heat transfer between the lateral blood vessel and the tissue requires , ,  \frac{\partial T_{b}^{m}}{\partial \mathbf{n}}=B_{i}\left(T_{w}^{m}-T_{b}^{m}\right)(13) The interfacial conditions between three skin tissue layers and the tumor are assumed to be perfectly thermal contact and are given by , ,  T_{1}=T_{2}, \quad k_{1} \frac{\partial T_{1}}{\partial z}=k_{2} \frac{\partial T_{2}}{\partial z}, \quad z=L_{1}(14a) T_{2}=T_{3}^{\text { out }}, k_{2} \frac{\partial T_{2}}{\partial z}=k_{3} \frac{\partial T_{3}^{\text {out}}}{\partial z}, z=L_{1}+L_{2}(14b) T_{2}=T_{3}^{i n}, k_{2} \frac{\partial T_{2}}{\partial z}=k_{t} \frac{\partial T_{3}^{i n}}{\partial z}, z=L_{1}+L_{2}(14c) T_{3}^{\text { out }}=T_{3}^{\text { in }}, k_{3} \frac{\partial T_{3}^{\text { out }}}{\partial z}=k_{t} \frac{\partial T_{3}^{m}}{\partial z}, z=L_{1}+L_{2}+L_{t}(14d) T_{3}^{\text {out}}=T_{3}^{\text { in }}, k_{3} \frac{\partial T_{3}^{\text { out }}}{\partial x}=k_{t} \frac{\partial T_{3}^{\text { in }}}{\partial x}, x=\frac{N X}{3}, \frac{2 N X}{3}(14e) T_{3}^{o u t}=T_{3}^{in}, k_{3} \frac{\partial T_{3}^{o u t}}{\partial y}=k_{t} \frac{\partial T_{3}^{in}}{\partial y}, y=\frac{N Y}{3}, \frac{2 N Y}{3}(14f) where T_{3}^{\text { out }}, T_{3}^{\text { in }}, N X and N Y represent the temperature outside and inside the tumor region, and the length of the plane in x and y directions, respectively. L_{t} is the length of the tumor in the z-direction. Because the blood flow in the vein is oriented against the arterial flow, the entrance of the blood to the vein is located at the seventh level, and the blood temperature is equal to the surrounding tissue temperature. The initial conditions are T_{l}=T_{0}, \quad t=0, \quad l=1,2,3, t(15) where T_{1}, T_{2}, T_{3} and T_{t} are tissue temperatures in skin layer one, two, three, and the tumor region, respectively, and T_{0} is the initial temperature in the tissue . 3. Numerical Method The above model would have to be solved numerically. In this study, we will use the finite difference method to solve due to the complexity of Eq. (1) as follows: \rho_{l} C_{l} \frac{\left(u_{l}\right)_{i j k}^{n+1}-\left(u_{l}\right)_{i j k}^{n}}{2}+W_{b}^{l} C_{b}^{l}\left[\frac{\left(u_{l}\right)_{i j k}^{n+1}-\left(u_{l}\right)_{i j k}^{n}}{2}-\left(u_{b}\right)_{o u t}\right]$=k_{l}\left(\delta_{x}^{2}+\delta_{y}^{2}+\delta_{z}^{2}\right) \frac{\left(u_{l}\right)_{i j k}^{n+1}+\left(u_{l}\right)_{i j k}^{n}}{2}+\left(Q_{l}\right)_{i j k}^{n+1 / 2}$

$l=1,2,3, t$(16)

where $\left(u_{l}\right)_{i j k}^{n}$ and $u_{b}$  denote the numerical approximations of $T_{l}(i \Delta x, j \Delta y, k \Delta z, n \Delta t)$ and $T_{b}$ , where $\Delta x, \Delta y, \Delta z,$ and $\Delta t$  are the spatial and temporal mesh sizes, and $i, j, k$ are integers, $0 \leq i \leq N_{x}, 0 \leq j \leq N_{y}, 0 \leq k \leq N_{z}^{l}$ ; so that $N_{x} \Delta x=N X, \quad N_{y} \Delta x=N y$ and $N_{z}^{l} \Delta z=N Z_{l}, l=1,2,3$ . In this mesh, we assume that $\left(u_{3}\right)_{i j k}^{n}=\left(u_{b}^{m}\right)_{i j k}$  when the grid point $(i, j, k)$  is in the mth level blood vessel. It should also be noted that $\delta_{x}^{2} u_{i j k}=\frac{u_{i+1 j k}-2 u_{i j k}+u_{i-1 j k}}{\Delta x^{2}}$ and so on, for the y and z-directions . Similarly, the finite-difference scheme can be obtained for solving Eqs. (5)-(6) together with Eqs. (8)-(15).

To determine the laser intensity $P_{0}$  so that an optimal temperature distribution can be obtained, we prespecify the temperature elevations to be obtained at the center and some locations in the perimeter on the skin surface. The reason that these locations are chosen is because the highest temperature is assumed to be around the center of the skin surface, and it is necessary to have the temperature in the perimeter below a certain threshold so as not to cause damage to the normal tissue. In addition, the temperature can be easily measured at these locations. By guessing initial laser intensity $P_{0}$ and pre-specifying a laser exposure pattern, one can solve the above equations to obtain a temperature field in the entire 3D skin structure. Once the temperatures, $\mathcal{U}_{c a l}^{i}$  are calculated at given locations ( $i=0,1, \ldots, M$ ) a least-square approach can be employed to minimize the difference between the pre-specified temperature elevation $\mathbf{\theta}_{\text { pre }}$  and the calculated temperature $\mathbf{U}_{\mathbf{c a l}}$  as follows :

$S\left(P_{0}\right)=\sum_{i=0}^{M}\left(\theta_{p r e}^{i}-u_{p r e}^{i}\right)^{2}, \quad i=0,1, \cdots, M$(17)

Minimizing $S\left(P_{0}\right)$ in Eq. (17), one obtains

$\frac{d}{d P_{0}} S\left(P_{0}\right)=-2 \sum_{i=0}^{M}\left[\frac{d u_{c a l}^{i}}{d P_{0}}\right]\left(\theta_{p r e}^{i}-u_{p r e}^{i}\right)=0$(18)

A new $P_{0}$  can be calculated iteratively as follows 

$P_{0}^{(J+1)}=P_{0}^{(J)}+\left(X^{t} X+\alpha^{*} \mathbf{I}\right)^{-1} X^{t}\left(\boldsymbol{\theta}_{\mathrm{pre}}-\mathbf{u}_{\mathrm{cal}}\right)$(19)

where $\alpha^{*}$  is a relaxation parameter,  $\mathbf{I}$ is an identity matrix, and X is the sensitivity coefficient matrix, which is an $(M+1) \times 1$  vector  of $X=\left[\frac{\partial u_{c a l}^{1}}{\partial P_{0}}, \frac{\partial u_{c a l}^{2}}{\partial P_{0}}, \cdots, \frac{\partial u_{c a l}^{M}}{\partial P_{0}}\right]^{t}$ , and $\boldsymbol{\theta}_{\mathrm{pre}}=\left[\theta_{p r e}^{0}, \theta_{p r e}^{1}, \cdots, \theta_{p r e}^{M}\right]^{t}, u_{\mathrm{pre}}=\left[u_{p r e}^{0}, u_{p r e}^{1}, \cdots, u_{p r e}^{M}\right]^{t}$

Thus, the computational procedure for optimizing the temperature distribution for a 3-D triple-layered skin structure with embedded countercurrent vascular network and nanoparticles can be described as follows:

Step 1: Guess the initial laser $P_{0}$  and obtain a temperature distribution based on the above numerical method.

Step 2: Increase the power with a small increment as $P_{0}+\Delta P_{0}$ .  Obtain a temperature distribution based on the above numerical method.

Step 3: Obtain the new power level by Eq. (19).  Continue iterating until an optimal temperature distribution is obtained.

We will test this method and numerical procedure in a 3-D triple-layered skin structure embedded with a vascular countercurrent network and nanoparticles in tumor region.  Currently, we are working on the testing and hope that we will obtain the result soon.

Acknoledgment

This research is partially supported by Louisiana Board of Regents’ Support Funds and the National Science Foundation (NSF EPS 1003897).

References

1. Zhang, L., Dai, W., and Nassar, R., “A numerical algorithm for obtaining an optimal temperature distribution in a 3d triple layered cylindrical skin structure,” Computer Assisted Mechanics and Engineering Sciences, 14:107-125, 2007.

2. Zhang, L., Dai, W., and Nassar, R., “A numerical method for optimizing laser power in the irradiation of a 3-d triple-layered cylindrical skin structure,” Numerical Heat Transfer, Part A, 48:21-41, 2005. DOI: 10.1080/10407780590929865.

3. Tang, X., Dai, W., Nassar, R., and Bejan, A., “Optimal temperature distribution in a three-dimensional triple-layered skin structure embedded with artery and vein vasculature,” Numerical Heat Transfer Part A, 50:809-834, 2006. DOI: 10.1080/10407780600669175.

4. Zhang, L., Dai, W., and Nassar, R., “A numerical method for obtaining an optimal temperature distribution in a 3-d triple-layered cylindrical skin structure embedded with a blood vessel,”   Numerical Heat Transfer, Part A, 49:765-784, 2006. DOI: 10.1080/10407780500506691.

5. Zeng, X., Dai, W., and Bejan, A., “Vascular countercurrent network for 3-d triple-layered skin structure with radiation heating,” Numerical Heat Transfer, Part A, 57:369-391, 2010. DOI: 10.1080/10407781003659599.

6. Vera, J., Bayazitoglu, Y., “Gold nanoshell density variation with laser power for induced hyperthermia,” International Journal of Heat and Mass Transfer, 52:564-573, 2009. DOI: 10.1016/j.ijheatmasstransfer.2008.06.036.

7. Pennes, H., “Analysis of tissue and arterial blood temperatures in the resting human forearm,” Journal of Applied Physiology, 1:93-122, 1948.

8. Liu, J., Chen, X., and Xu, L., “New thermal wave aspects on burn evaluation of skin subjected to instantaneous heating,” IEEE Transactions on Biomedical Engineering, 46:420-428, 1999.

9. Liu, J., “Preliminary survey on the mechanisms of the wave-like behaviors of heat transfer in living tissues,” Forschung im Ingenieurwesen, 66:1-10, 2000.

10. Lorente, S., and Bejan, A., Design with Constructal Theory. Hoboken:Wiley, 2008.

11. Zhou, J., and Liu, J., “Numerical study on 3-d light and heat transport in biological tissues embedded with larger blood vessels during laser-induced thermotherapy,” Numerical Heat Transfer, Part A, 45:415-449, 2004. DOI: 10.1080/10407780490269030.

12. Bejan, A., “The tree of convective heat streams: its thermal insulation function and the predicted 3/4 –power relation between body heat loss and body size,” International Journal of Heat and Mass Transfer, 44:699-704, 2001. DOI: 10.1016/S0017-9310(00)00138-1.

13. Bejan, A., and Lorente, S., “Constructal theory of generation of configuration in nature and engineering,” Journal of Applied Physics, 100:041301, 2006. DOI: 10.1063/1.2221896.

14. Silva, A., Lorente, S., Bejan, A., “Constructal multi-scale tree-shaped heat exchangers,” Journal of Applied Physics, 96:1709-1718, 2004. DOI: 10.1063/1.1766089.

15. Dai, W., Bejan, A., Tang, X., Zhang, L., and Nassar, R., “Optimal temperature distribution in a three dimensional triple-layered skin structure with embedded vasculature,” Journal of Applied Physics, 99:104702, 2006. DOI: 10.1063/1.2199193.

16. Han, J., and Jensen, K.F., “Combined experimental and modelling studies of laser-assisted chemical vapour deposition of copper from copper (i)-hexafluoracetylacetonate trimethylvinylsilane,” Journal of Applied Physics, 75:2240-2250, 1994.