\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
International Conference on Applications of Mathematics to Nonlinear Sciences,\newline
\emph{Electronic Journal of Differential Equations},
Conference 24 (2017), pp. 75--83.\newline
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu}
\thanks{\copyright 2017 Texas State University.}
\vspace{8mm}}

\begin{document} \setcounter{page}{75}
\title[\hfilneg EJDE-2017/conf/24\hfil
Heat regulation in human body with dermal tumor]
{2D model on heat regulation in human body with dermal tumor}

\author[K. Nazir, M. A. Khanday, B. A. Ganai \hfil EJDE-2017/conf/24\hfilneg]
{Khalid Nazir,  Mukhtar A. Khanday,  Bashir A. Ganai }

\address{Khalid Nazir  \newline
Department of Mathematics,
University of Kashmir,
Srinagar, Jammu and Kashmir, 190006, India}
\email{khldnzr99@gmail.com}

\address{Mukhtar A. Khanday (corrresponding author) \newline
Department of Mathematics,
University of Kashmir,
Srinagar, Jammu and Kashmir, 190006, India}
\email{khanday@gmail.com}

\address{Bashir A. Ganai \newline
CORD, University of Kashmir,
Srinagar, Jammu and Kashmir, 190006, India}
\email{bbcganai@gmail.com}


\thanks{Published November 15, 2017.}
\subjclass[2010]{92B05, 92C10, 92C42}
\keywords{Thermoregulation; mathematical model; dermal tumor}

\begin{abstract}
 To predict the thermal fluctuations of a finite biological tissue in presence
 of tumor, an attempt has been made to formulate a 2D mathematical model based
 on the cross sectional temperature distribution in the tissues of the human
 limbs and variational finite element approach has been employed to establish
 the solution of the model. It is assumed that the dermal region of the human
 body is hosting tumor/cancerous cells. Thermal distribution at the tumor
 region with respect to different input parameters has been computed
 using MATLAB software. The physiological and bio-physical parameters
 like metabolic  heat generation, blood mass flow rate and thermal conductivity
 are assumed   to vary in the sub regions independently. The model describes
 the exchange of heat between the internal biological tissues and other
 surrounding media.  Thermal fluctuations at the targeted regions were
 obtained with respect to  various power densities of the heating sources.
 The results obtained may be  helpful for various cases of practical interest
 especially in the treatment of cancerous tumors and in local hyperthermic
 therapies.
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\allowdisplaybreaks

\section{Introduction} \label{s1}

Skin is an important element in the mechanism of temperature regulation,
it acts as a barrier between internal system of the body and the surroundings
by preventing excessive loss of water, organic and inorganic materials
that are necessary for the maintenance of internal homoeostasis and normal
functioning of the cells. In 1962, the structure and the functioning of the
skin was first studied by Montagna \cite{mon}. The physiology of the skin
also plays a vital role for various practical interests. Depending upon the
physiological properties of the skin, epidermis and dermis are its main layers.
However, below dermis subcutaneous tissues also plays an important role in
 thermoregulation.

Skin temperature also plays a vital role to keep the body temperature within
certain stable range, even when the surrounding temperature is unstable.
 The hypothalamus regulates the body temperature by thermoregulatory
mechanisms such as vasomotion, shivering and regulatory sweating.
It receives inputs from the central and peripheral temperature receptors,
situated at the core and the outer shell respectively. The control of body
temperature depends upon the balance between the cooling and heating effects.
The normal temperature of the body is conducted in blood stream to the skin,
 where it is released into the environment by means of conduction, convection
and radiation and evaporation. Conduction takes place within the tissue governed
by the Fourier law. Radiation takes place at the outer most surface governed
by Stefan's law. Convection transfer by the blood circulation through the
vessels is the main mode of transfer of heat from the core to the periphery
 and thermal conductivity solely depends upon rate of blood flow through the
 cardiovascular system. Since densities of blood vessels are very high, it is
hard to compute the detailed temperature in even a small part of the body.
Fortunately, the effects of blood vessels on body tissue temperature has been
 described collectivily by various researchers with some success.

 In 1948, Pennes \cite{pen} devised a bio-heat equation, were he  described the
effect of blood perfusion and metabolic heat generation on heat transfer within
the living tissue. The purpose of Pennes' study \cite{pen} was
``to evaluate the applicability of heat flow theory to the forearm of the
human body in basic terms of local rate of tissue heat production and volume
flow of blood''. In 1958, the paper by Pennes \cite{pen} was revisited
by Wisler \cite{wis} because the experimental data seem to be at a
variance with the theoretical results.  In 1986, Knudsen and Overgaard \cite{over}
identified the new thermal model for human tissue. They have calculated
temperature profiles within the human forearm which were in agreement with
the results obtained by Pennes \cite{pen}. In 1987, Saxena and Bindra \cite{saxbin}
 used quadratic shape functions in variational finite element method to study
heat distribution in cutaneous and subcutaneous tissues (SST).
Recently Khanday et al.\ \cite{khalid} has thoroughly discussed and established
few models on bio-heat equation through eigenvalue expansion approach.
They have discussed the effect of atmospheric temperature on heat regulation
 at the peripheral regions of the human tissue.

Hyperthermia is the procedure of excessive tissue temperature of the body organ
or the whole body. The effectiveness of hyperthermia depends upon the value
of the evaluated temperature and exposure time. In order to determine accurately
the temperature field over the entire region, many numerical and experimental
 methods have been developed in order to solve the bio-heat equation.
Gupta et al.\ \cite{gup} obtained approximate analytic solution of Pennes'
 bio-heat equation in thermal therapy. Ng et al.\ \cite{ng} did the
parametric analysis of thermal profiles within heated human skin using
boundary element method. Kengne et al.\ \cite{ken} gave analytical and numerical
solution of temperature distribution with oscillatory surface and spatial heating.
 In the last few decades attempts were made further by Saxena and Pardasani
\cite{Sap2},\cite{Sap1} to study the effect of dermal tumors on temperature
distribution in skin and subcutaneous tissue with variable blood flow.
They have discussed the radial heat flow in skin and underlying tissue
layers of spherical regions of human or animal body. The model was further
extended and studied by Khanday et al. \cite{kma1}. They established a
mathematical model for the treatment of cancerous tumors based on local
hyperthermia with external heat source. Agrawal et al. \cite{mam3} also
worked on the thermal disturbances in dermal regions of human limbs
 involving metastasis of tumors.

 The proposed work in this paper has many applications in various heat-transfer
related problems. The behavior of tumor tissue temperature with respect to
ambient conditions and in relation with external heat source has been studied.
It is assumed that dermal regions host tumor cells in the human body.
Variational finite element method has been employed for the numerical
solution of Pennes' bio-heat equation with appropriate boundary conditions
to understand the thermal processes in the tumor tissues.

\section{Mathematical formulation}\label{s2}

The bio-heat equation developed by Pennes \cite{pen} is one of the earliest
models for energy transport in biological tissues. If $T(x,y,t)$ denotes
the temperature at time $t$ and at a position $(x,y)$, then the heat
flow in tissues for two dimensional unsteady state is given below:
\begin{equation} \label{e2.1}
\rho c \frac{\partial T}{\partial t}
= k\Big( \frac{\partial ^ 2 T }{\partial x ^ 2}
+\frac{\partial ^ 2 T }{\partial y ^ 2}\Big)
+ \rho _b m _ b c_b (T_a - T ) + S+Q_{\rm ex}\,.
\end{equation}

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.98\textwidth]{fig1} % title1
\label{Tissue}
\end{center}
\caption{Triangular element arrangement for annular cross section of human limb.}
\end{figure}

 The left-hand side represents the storage of heat through the tissue and the first
terms on the R.H.S represents the conduction of heat through the 2D
structure of tissue, the middle term represents the blood perfusion
 and the term $S$ represents the heat generation due to metabolism.
The last term denotes the external heat source  $Q_{\rm ex}=Q(x,y,t)$.
The parameters $\rho$, $c$ and $k$ respectively denote density,
specific heat and thermal conductivity of the tissue. $m_b$, $c_b$ and $T_a$
are the blood mass flow rate, specific heat of the blood and arterial blood
temperature respectively. The outer surface of the tissue (skin) is exposed
to the environment so there is continuous exchange of heat flux between
the two media. Thus the boundary condition at the outer surface of the body
is governed by Newton's cooling law given by
\begin{equation} \label{e2.2}
-k\frac{\partial T}{\partial n} =h (T - T_A)+LE\,,
\end{equation}
where $h$ is the heat transfer coefficient, $T_A$ is atmospheric temperature
and $L, E$ are the latent heat and evaporation rates at the skin surface
respectively.

 Human beings maintain relatively constant temperature at the inner core of
the limb, therefore the inner boundary conditions are given as
\begin{gather}
T_b=37^\circ C,\quad  \text{at } t \geq 0, \label{e2.3} \\
T(x,y)=T_b,\quad  5d \leq x \leq 13d,\quad 0 \leq y \leq 2d, \label{e2.4}\\
\frac{\partial T}{\partial y}=0,\quad  0 \leq x \leq 5d,\quad
 13d \leq x \leq 18d, \label{e2.5}
\end{gather}
where the constant $d$ denotes the length of the domain, and it can assume
any value between $(0\leq d < 1)$ depending on the sample of the limb under study.

We assume that initially the temperature distribution is relatively equal
to the body core temperature, i.e.,
\begin{equation} \label{e2.6}
 T(x,0)=T_b\,.
\end{equation}


\section{Solution of the model}\label{s3}

The domain of the underlying tissue is assumed to be annular in geometry.
Keeping the complex structure of the dermal regions into account, the
discretization of the domain has been approximated by the assembly of
triangular elements of different sizes. Different types of parameters
were considered for subregions such as epidermis, dermis and hypodermis.
Also due to the complex structure of the domain, the exact solution of
the formulated model is impossible to obtain. Therefore, the numerical
solutions of boundary value problems are more effective and reasonable
to get the solution of the formulated model. The variational finite
element (VFEM) is one of the relevant numerical methods for handling
irregular geometrical domains.

On comparing equation \eqref{e2.1} with the Euler-Lagrange differential equation,
the variational integral is given below
\begin{equation} \label{e3.1}
\begin{aligned}
I&=\frac{1}{2}\iint_{{A}}\Big[k\Big\{\Big(\frac{\partial T}{\partial x}
\Big)^2+\Big(\frac{\partial T}{\partial y}\Big)^2\Big\}
+\rho_b c_b \Big(T-T_b\Big)^2 +\rho c\frac{\partial}{\partial t}T^2 \\
&\quad -2SQ_{\rm ex} T\Big]\,dx\,dy
 +\int_{\Omega}\Big[ h\Big(T-T_a\Big)^2+2LE T\Big] d\Omega\,,
\end{aligned}
\end{equation}
where $A$ and $\Omega$ is the cross sectional area and the boundary of 
the cylindrical limb respectively. The region under study has  
been divided into $60$ triangular elements based on the anatomy of 
the dermal region. The element-wise variational integrals are given by
\begin{equation} \label{e3.2}
\begin{aligned}
I_i&=\frac{1}{2}\iint_{{A}}\Big[k_i\Big\{\Big(
\frac{\partial T^i}{\partial x}\Big)^2
+\Big(\frac{\partial T^i}{\partial y}\Big)^2\Big\}+\rho_b c_b
\Big(T^i-T_b\Big)^2 +\rho_i c_i\frac{\partial}{\partial t}
\Big(T^i\Big)^2 \\
&\quad -2S^iQ_{\rm ex}^i T^i\Big]\,dx\,dy
 +\int_{\Omega}\Big[ h\Big(T^i-T_a\Big)^2+2LE T^i\Big] d\Omega\,,
\end{aligned}
\end{equation}
where $i=1,2,3,\dots ,60$.

To solve the variational integral $I_i$, the following shape functions were used
\begin{equation} \label{e3.3}
T^i (x,y)=C_1 ^i +C_2 ^i x +C_3 ^i y
\end{equation}
In the matrix form, these equations can be written as
\begin{equation} \label{e3.4}
T^{(i)}=P^{(u)}C^{(i)},
\end{equation}
where
$P^{(u)}= (1 \;  x \; y)$.

These values of $C_1 ^i,C_2 ^i,C_3 ^i$ can be obtained from the 
system of equations
\[
\begin{pmatrix}
T_i \\
T_j \\
T_k
\end{pmatrix}=
\begin{pmatrix}
1 & x_i & y_i \\
1 & x_j & y_j\\
1 & x_k & y_k
\end{pmatrix}
\begin{pmatrix}
C_1^i \\
C_2^i \\
C_3^i
\end{pmatrix},
\]
where $T_i,T_j,T_k$ are the nodal temperatures of the $i^{th}$ element.

Or,
\begin{equation} \label{e3.5}
U=P^{(i)}C^{(i)},
\end{equation}
where
\[
U=\begin{pmatrix}
T_i \\
T_j \\
T_k
\end{pmatrix}, \quad
P^{(i)}=\begin{pmatrix}
1 & x_i & y_i \\
1 & x_j & y_j\\
1 & x_k & y_k
\end{pmatrix}, \quad
C^{(i)}=\begin{pmatrix}
C_1^i \\
C_2^i \\
C_3^i
\end{pmatrix}.
\]

On solving  equation \eqref{e3.5} for $C^{i}$, we get
\begin{equation} \label{e3.6}
C^{(i)}=(P^{(i)})^{-1} U\,.
\end{equation}
Now from equation \eqref{e3.4} we have
\begin{equation} \label{e3.7}
T^{(i)}=P^{(u)}(P^{(i)})^{-1} U\,.
\end{equation}
Using this equation to evaluate $I_i$ over each element and later assembling
these integrals to obtain $I=\sum_{1}^{60}I_i $. $I$ can be optimized with
respect to the nodal values $ T_i$ to get the system of algebraic equations,
 written in matrix form as
\begin{equation} \label{e3.8}
XT+Y \dot T=Z,
\end{equation}
where $X$ and $Y$ are the matrices of order $(44\times 44)$, $Z$, $T$ and
$\dot T$ are the vectors of order $(44\times 1)$ defined as
\[
X=\begin{pmatrix}
a_{1,1} & a_{1,2} &\dots  & a_{1,44}\\
a_{2,1} & a_{2,2} & \dots  & a_{2,44}\\
. & . & . & .\\
. & . & . & .\\
. & . & . & .\\
a_{44,1} & a_{44,2} & \dots  & a_{44,44}
\end{pmatrix}, \quad
Y=\begin{pmatrix}
b_{1,1} & b_{1,2} &\dots  & b_{1,44}\\
b_{2,1} & b_{2,2} & \dots  & b_{2,44}\\
. & . & . & .\\
. & . & . & .\\
. & . & . & .\\
b_{44,1} & b_{44,2} & \dots  & b_{44,44}
\end{pmatrix},
\]
$T=\big(
T_1 \; T_2 \; \dots \; T_{44}\big)'$,
$\dot{T} =\big(\dot{T_1}  \; \dot{T_2} \; \dots  \; \dot{T_{44}}\big)'$,
$Z=\big(Q_1 \; Q_2 \; \dots \; Q_{44} \big)'$,
where primes denotes the transpose of the vectors and dot denotes the
derivative of temperature with respect to time.

Now to solve the system of ordinary equations, the Crank Nicholson implicit
 scheme has been employed to obtain
\begin{equation} \label{e3.9}
\Big( X+\frac{\varDelta t}{2}Y \Big)T^{i+1}
=\Big( X-\frac{\varDelta t}{2}Y \Big)T^{i}+\varDelta t Z\,.
\end{equation}
On solving this equation, we obtain the values of the nodal temperatures
for the different values of the ambient conditions $t_{\rm ex}=30 s,40 s$
 as given in Table \ref{table3}.

\section{Numerical results}\label{s4}

The numerical and physiological values of the parameters used in the solution
of the model are given in the Tables \ref{table1} and \ref{table2}.

\begin{table}[htb]
\caption{Variants of Heating}
\label{table1}
\begin{center}
\begin{tabular}{|c |c |c|}
\hline
 Variant No. & Power density$[MW/m^2]$ & Heating duration$[s]$ \\
\hline
1 & $7.0$ & $30,40$ \\
\hline
2 & $3.5$ & $30,40$\\
\hline
3 & $1.0$ & $30,40$\\
\hline
\end{tabular}
\end{center}
\end{table}

\begin{table}[htb]
\caption{Numerical and Physiological values of Parameters}
\label{table2}
\begin{center} \renewcommand{\arraystretch}{1.1}
\begin{tabular}{|l| l| l|}
\hline
\qquad Parameters & Value & Units \\
\hline
Thermal conductivity of Epidermis$(k_1)$ & $0.030$ & $cal/cm-min^{0} C$ \\
\hline
Specific heat  of Epidermis$(c_1)$ & $0.83$ & $cal/g$\\
\hline
Density of Epidermis$(\rho_1)$ & $1.05$ & $g/cm^{3}$\\
\hline
Thermal conductivity of Dermis$(k_2)$ & $0.060$ & $cal/cm-min^{0} C$ \\
\hline
Specific heat  of Dermis$(c_2)$ & $0.95$ & $cal/g$\\
\hline
Density of Dermis$(\rho_2)$ & $1.09$ & $g/cm^{3}$\\
\hline
Thermal conductivity of SST$(k_3)$ & $0.090$ & $cal/cm-min^{0} C$ \\
\hline
Specific heat  of SST$(c_3)$ & $0.95$ & $cal/g$\\
\hline
Density of SST$(\rho_3)$ & $1.15$ & $g/cm^{3}$\\
\hline
Heat transfer coefficient$(h)$ & $0.009$ & $cal/cm^{2}-min^{0} C$\\
\hline
Body core temperature$(T_4=T_b)$ & $37$ & $^{0}C$\\
\hline
Artrial temperature$(T_a)$ & $37$ & $^{0}C$\\
\hline
Density of blood$(\rho_b)$ & $1060$ & $kg/m^{3}$\\
\hline
Specific heat of blood$(c_b)$ & $3500$ & $Ws/kg/m^{3}$\\
\hline
Latent heat$(L)$ & $2.42$ & $J/g$\\
\hline
Rate of sweat evaporation$(E)$ & $0.001$ & $g/cm^3/min$\\
\hline
\end{tabular}
\end{center}
\end{table}

\begin{table}[htb]
\caption{Nodal temperatures of the tissue at different atmospheric
temperatures when exposed to external heat for $t_{\rm ex}=30s$.}

\label{table3} 
\begin{center} \renewcommand{\arraystretch}{1.1}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|}
\hline
 & $T_1$ & $T_2$ & $T_3$  & $T_4$ & $T_5$ & $T_6$ & $T_7$ & $T_8$  \\
$T_A = 25^\circ$C & 37.78 & 35.50 & 29.91 & 24.96 & 37.89 & 34.59 & 28.99 & 24.99 \\
$T_A = 30^\circ$C & 37.80 & 35.62 & 29.93 & 29.91 & 37.97 & 34.96& 33.98 & 29.97\\
\hline

&$T_9$  & $T_{10}$ & $T_{11}$ & $T_{12}$ & $T_{13}$ & $T_{14}$ & $T_{15}$ & $T_{16}$ \\
$T_A = 25^\circ$C &   36.99 & 34.41 &28.98 & 25.01& 37.00 & 44.01 & 25.00 & 37.00 \\
$T_A = 30^\circ$C &  37.00 & 35.49 & 34.01 & 29.89 & 37.03 & 44.50 & 29.76 & 37.23 \\
\hline

 & $T_{17}$ & $T_{18}$ & $T_{19}$ & $T_{20}$ & $T_{21}$ & $T_{22}$ & $T_{23}$ & $T_{24}$ \\
$T_A = 25^\circ$C & 43.02 & 41.50 & 24.78 & 36.98 & 36.98 & 35.23 & 29.99 & 24.75 \\
$T_A = 30^\circ$C & 44.01 & 43.50 & 29.99 & 37.01 & 37.01 & 35.60 & 33.98 & 29.98\\
\hline

&$T_{25}$ & $T_{26}$ & $T_{27}$ & $T_{28}$ & $T_{29}$ & $T_{30}$  & $T_{31}$ & $T_{32}$ \\
$T_A = 25^\circ$C &   36.92 & 35.00 & 29.01 & 24.88 & 37.01 & 35.32 & 28.34 & 24.67 \\ 
$T_A = 30^\circ$C & 36.98 & 35.72 & 33.76 & 29.99 & 37.03 & 35.14 & 33.23 & 29.98  \\
\hline

&$T_{33}$ & $T_{34}$ & $T_{35}$ & $T_{36}$ & $T_{37}$  & $T_{38}$ & $T_{39}$ & $T_{40}$ \\
$T_A = 25^\circ$C &  37.03 & 34.45 & 28.54 & 24.93 & 36.68 & 34.76 & 28.93 & 24.95 \\
$T_A = 30^\circ$C & 37.15 & 35.07 & 32.92 & 29.97 & 36.78 & 35.01 & 32.73 & 29.95\\
\hline

&$T_{41}$ & $T_{42}$ & $T_{43}$ & $T_{44}$ &&&&\\
 $T_A = 25^\circ$C &  36.95 & 35.35 & 29.76  & 24.98 &&&&\\
$T_A = 30^\circ$C & 37.00 & 35.01 & 32.27 & 29.98  &&&&\\
\hline
\end{tabular}
\end{center}
\end{table} 


\section{Discussion and conclusion}\label{s5}

\begin{figure}[htb]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig2} % pro311
\end{center}
 \caption{{Thermal distribution of tumor tissue ($21^{th}$ element)
at $t_{\rm ex}=30s$.}}
\label{Tissue2}
\end{figure}

\begin{figure}[htb]
\begin{center}
 \includegraphics[width=0.6\textwidth]{fig3} % pro441
\end{center}
 \caption{Thermal distribution of tumor tissue ($21^{th}$ element) at
$t_{\rm ex}=40s$.}
\label{Tissue3}
\end{figure}

The thermal behaviour of biological tissues hosting tumor cells at the dermal
regions has been analysed using local hyperthermic therapy. Variational
finite element method has been employed to solve the 2D mathematical model
based on Pennes' bio-heat equation. Laser therapy is one of the most commonly
used therapy in which the local part of the underlying tissue is being heated
using some spherical concave lenses and all the energy radiations can be
focused on the target area of the tissue. In this connection the more
effective heating source is the point heating source which is commonly
written in terms of dirac delta function as
\[
Q(x,y,t)=P(t)\delta(x-x_0,y-y_0)
\]
where $P(t)$ is the point heating power for the treatment of local tumors,
 $\delta$ is the Dirac delta function and $(x_0,y_0)$ is the position of
the point heating source. This type of heating gains much more
importance in the hyperthermic treatment. In this type of heating style
only the selected tissue (tumor tissue) is subjected to the heating source,
where as the temperature of the surrounding healthy tissue remains under
the threshold temperature. In this study, we have taken different values of
$P(t)$ (see Table-(1)). Also, it has been assumed that the heating region
is the $21^{th}$ element of our discretized domain. The variation of
temperature at different values of power density and for different exposing
times on the tumour region is shown in Figures $2\,\text{and}\,3$.
It is clear from figures that during heating of the tissue all curves are
almost similar up to the peak value, as soon the time crosses the exposing
time the curves show a little difference and in case of long heating,
the maximum tissue temperature  also differs.

 The main application of the study is that it gives some realistic values
of temperature profiles as compared to the numerical results established
by various researchers including Weinbaum et al.\ \cite{wein1}, Saxena
and Pardasani \cite{Sap1}.
The estimation of the temperature profiles may be helpful for various
clinical purposes especially during the treatment of cancer through
radiotherapy and other local hyperthermic approaches. The proposed
model can contribute to medical application by understanding the
macro-scale tissue response under the external heating. These results
are also useful in the design and characterization of hyperthermic
treatment. Moreover the model may be useful for optimization of a
treatment by maximizing the therapeutic effect and minimizing side effect,
proposing new treatment strategies and evaluation of their outcomes.


\subsection*{Acknowledgements}
The authors are highly thankful to the UGC, New Delhi, for their financial
support under the grant of Major Research Project No. 42-11/2013(SR).


\begin{thebibliography}{99}

\bibitem{mam3} Agrawal, M.; Adlakha, N.; Pardasani, K. R.;
\emph{Thermal disturbances in dermal regions of human limbs involving metastasis
 of tumors}. International Mathematical Forum, Vol. (39), (2010), 1903-1914.

\bibitem{gup} Gupta, P. K.;  Singh, J.;  Rai, K. N.;
\emph{Numerical simulation for heat transfer in tissues during thermal therapy}.
Journal of Thermal Biology, Vol. 35(6), (2010), 295-301.

\bibitem{ken} Kengne, E.; Lakhssassi, A; Vaillancourt, R.;
\emph{Temperature distribution in living biological tissue simultaneously
subjected to oscillatory surface and spatial heating: analytical and numerical
 analysis}. International Mathematical Forum, Vol. 7, (2012), 2373-2392.

\bibitem{kma1} Khanday, M. A.; Mir, Aijaz; Aasma, Rafiq;
\emph{Mathematical analysis on the treatment of cancerous tumors based
on local hyperthermia}. Journal of Energy, Heat and Mass Transfer,
Vol. 35(2), (2013), 0970-9991.

\bibitem{kma2} Khanday, M. A.; Nazir, Khalid;
\emph{Temperature distribution in biological tissue under the influence of
external heat source}, Bull. Cal. Math. Soc., Vol. 107, (2015), 145-152.

\bibitem{khalid}  Khanday, M. A.; Nazir, Khalid;
\emph{Eigenvalue expansion approach to study bio-heat equation}.
J. Multiscale Modelling, Vol. 07, (2016), 1650002-8 pages.

\bibitem{mon} Montanga, W.;
\emph{The structure and function of skin}, 3rd edition, Academic Press,
 Vol. 3,(1962), 448.

\bibitem{ng} Ng, E. Y.; Tan, H. M.;  Ooi, E. H.;
\emph{Prediction and parametric analysis of thermal profiles within heated human
skin using the boundary element method}. Philos. Trans. A. Math phy.
Eng Sci., Vol. 368, (2010), 655-678.

\bibitem{over} Kundsen, M.; Overgaard, J.;
\emph{Identification of thermal model for human tissues},
IEEE. Trans. Biomed. Eng. Vol. 33, (1986), 477-485.

\bibitem{pen} Pennes, H. H.;
\emph{Analysis of tissue and arterial blood temperature in the resting human
 forearm}. J. Appl. Physiol., Vol. 1, (1948), 93-122.

\bibitem{saxbin} Saxena, V. P.; D. Bindra, J.;
\emph{Quadratic shape functions in variational finite element approach
to heat distribution in cutaneous and subcutaneous tissues},
Ind. J. of Pure and Appl. Math., Vol. 18 (1987), 846-855.

\bibitem{Sap2} Saxena, V. P.; Pardasani, K. R.;
\emph{Steady state radial heat flow in skin and underlying tissue layers
of spherical regions of human or animal body}. Ind. J. Tech., Vol. 25, (1978),
501-505.

\bibitem{Sap1} Saxena, V. P.; Pardasani, K. R.;
\emph{Effect of dermal tumor on temperature distribution in skin and variable
blood flow}, Bull., Math  Biol., Vol 53, (1991), 1251-1265.

\bibitem{wein1} Weinbaum, S.;  Jiji, L. M.; Lemons, D. E.;
\emph{Theory and experiment for the effect of vescular micro-structure on
surface tissue heat transfer - Part 1. Anatomical foundation and model
 conceptualization}, J. Biomech. Eng., Vol 106, (1984), 321-330.

\bibitem{wein} Weinbaum, S.; Jiji, L. M.;
\emph{A new simplified bio-heat equation for the effect of blood flow on local
 average tissue temperature}, J. Biomech. Eng., Vol. 107, (1985), 131-139.

\bibitem{wis} Wissler, E. H.;
\emph{Pennes' paper revisited}, Journal of Applied Physiology,
 Vol. 85, (1998), 35-41.

\end{thebibliography}

\end{document}