^{1}Department of Applied Mathematics, University of Calcutta, , Kolkata-700 009, India

^{1}Corresponding author

Mathematical Models in Engineering, Vol. 2, Issue 2, 2016, p. 151-157.
https://doi.org/10.21595/mme.2016.18024

Received 15 November 2016; received in revised form 19 December 2016; accepted 20 December 2016; published 31 December 2016

Copyright © 2016 JVE International Ltd.

In this paper, a new mathematical model for Pennes’ bioheat equation using the new memory-dependent derivative is established. The one-dimensional thermal behavior in living tissue subject to instantaneous surface heating is investigated. Numerical calculations are performed to study the temperature transients in the skin exposed to instantaneous surface heating. Numerical results are plotted in the form of two-dimensional graphs and discussed. In this novel model, the time delay parameter is a new indicator of bio-heat efficiency in living tissues.

**Keywords: ** Pennes’ bioheat transfer equation, Fourier’s law of heat conduction, memory-dependent derivative, time-delay parameter, Kernel function.

The memory-dependent derivative is defined in an integral form of a common derivative with a Kernel function on a slipping interval [1]. Thus, this kind of definition is better than the fractional one for reflecting the memory effect (instantaneous change rate depends on the past state). Its definition is more intuitionistic for understanding the physical meaning and the corresponding memory dependent differential equation has more expressive force.

The fractional order calculus is more important in many applications in different areas of engineering and modern physics. We can use the fractional calculus to modify the classical Fourier law of heat conduction. The commonly used definition of fractional derivative is the Caputo type [2, 3]:

(1)

${D}_{b}^{\alpha}f\left(t\right)=\underset{b}{\overset{t}{\int}}{K}_{\alpha}(t-\xi ){f}^{\left(m\right)}\left(\xi \right)d\xi ,$
with:

(2)

${K}_{\alpha}\left(t-\xi \right)=\frac{(t-\xi {)}^{m-\alpha -1}}{\mathrm{\Gamma}\left(m-\alpha \right)},$
where $b$ is a fixed real number, $m$ is an integer which satisfies $m-1<\alpha \le m$, $\mathrm{\Gamma}$ is the Euler’s Gamma function and ${f}^{\left(m\right)}$ is the common $m$-order derivative which has specific physical meaning and ${K}_{\alpha}(t-\xi )$ is the Kernel of function.

From Eq. (1), we see the $\alpha $-th order fractional derivative at time $t$ is not defined locally, it relies on the total effects of the commonly used $m$-order integer derivative on the interval [$a$, $t$]. So, it can be used to describe the variation of a system in which the instantaneous change rate depends on the past state, which is called “memory effect” in a visualized manner [4]. This is why the fractional theory is widely used.

The first order memory-dependent derivative of a function $f\left(\xi \right)$ was introduced by Wang and Li [1]. They defined it in an integral form of a common derivative with a Kernel of function on slipping interval in the form:

(3)

${D}_{\omega}f\left(t\right)=\frac{1}{\omega}\underset{t-\omega}{\overset{t}{\int}}K(t-\xi )f\text{'}\left(\xi \right)d\xi ,$
where $\omega $ is the time-delay parameter and $K(t-\xi )$ is the Kernel of function which can be chosen freely [5]. The Taylor theorem in terms of memory-dependent derivatives can found in [5] as:

(4)

$f\left(x,t+\omega \right)=f\left(x,t\right)+\sum _{n=1}^{\infty}\frac{{\omega}^{n}}{n!}{D}_{\omega}^{n}f\left(x,t\right),$
where:

(5)

${D}_{\omega}^{n}f\left(x,t\right)=\frac{1}{\omega}\underset{t-\omega}{\overset{t}{\int}}K\left(t-\xi \right){f}_{\xi}^{\left(n\right)}\left(x,\xi \right)d\xi .$
As in [5], we can expand the heat flux vector $\overrightarrow{q}(x,t+\omega )$ by Taylor’s series of memory-dependent derivative and neglected the terms up to the first order in time-delay $\omega $ to get:

(6)

$\overrightarrow{q}\left(\overrightarrow{x},t+\omega \right)=\left(1+\omega {D}_{\omega}\right)\overrightarrow{q}\left(\overrightarrow{x},t\right).$
Recently, an interesting application of memory-dependent derivative is given by Yu et al. [6]. They introduced the memory-dependent derivative (MDD) instead of fractional calculus into the rate of heat flux in Lord-Shulman theory [7] of generalized thermoelasticity to denote memory-dependence as:

(7)

${q}_{i}+\omega {\dot{q}}_{i}=-k{\theta}_{,i}.$
Eq. (7) has more clear physical meaning.

In this paper, a new mathematical model for Pennes’ bioheat equation using the new memory-dependent derivative is established. The one-dimensional thermal behavior in living tissue subject to instantaneous surface heating is investigated. Numerical calculations are performed to study the temperature transients in the skin exposed to instantaneous surface heating. Numerical results are plotted in the form of two-dimensional graphs and discussed. In this novel model, the time delay parameter is a new indicator of bio-heat efficiency in living tissues.

Heat transfer in biological systems is usually modelled by the Pennes’ bioheat equation [8] based on the classical Fourier’s law (which specifies a linear relationship between heat flux and temperature gradient):

(8)

$\overrightarrow{q}\left(\overrightarrow{x},t\right)=-k\overrightarrow{\nabla}T\left(\overrightarrow{x},t\right),$
as:

(9)

$\rho C\dot{T}\left(\overrightarrow{x},t\right)=-\overrightarrow{\nabla}\cdot \overrightarrow{q}\left(\overrightarrow{x},t\right)+Q\left(\overrightarrow{x},t\right),$
where $T(\overrightarrow{x},t)$ is the temperature of living tissue and $Q(\overrightarrow{x},t)$ is the volumetric heat generated by metabolism and blood perfusion, given by:

(10)

$Q\left(\overrightarrow{x},t\right)={G}_{B}{C}_{B}\left[{T}_{B}-T\left(\overrightarrow{x},t\right)\right]+{Q}_{m},$
where ${G}_{B}$ is the blood perfusion, ${C}_{B}$ is the volumetric specific heat of blood, ${T}_{B}$ is the artery temperature and ${Q}_{m}$ is the metabolic heat source. A well-known problem with Fourier’s law is that it yields infinitely fast propagation of thermal signal, incompatible with physical reality and physiological considerations in a transient process. This equation implies an instantaneous thermal energy deposition in medium, i.e. any local temperature disturbance causes an instantaneous perturbation in temperature at each point in medium.

As explained above, heat pulses obtained by the classical bio-heat conduction equation propagate at infinite speed. Much attention has been devoted to modifying the classical heat conduction equation to ensure finite speed pulse propagation. In mathematical terms, the governing partial differential equation is transformed from parabolic to hyperbolic type [9]. A general form of bioheat transfer model in living tissues based on the generalized Fourier’s law [7]:

(11)

$\overrightarrow{q}\left(\overrightarrow{x},t\right)+\tau \dot{\overrightarrow{q}}\left(\overrightarrow{x},t\right)=-k\overrightarrow{\nabla}T\left(\overrightarrow{x},t\right),$
can be put in the form:

(12)

$\rho C\left[\dot{T}(x,t)+{\tau}_{0}\ddot{T}(x,t)\right]-\left[Q\left(x,t\right)+{\tau}_{0}\dot{Q}\left(x,t\right)\right]=k{\nabla}^{2}T\left(x,t\right),$
where $\tau $ is a relaxation time parameter, $k$ is the tissue thermal conductivity and $C$ is the specific heat.

From a mathematical viewpoint, the Fourier’s law Eq. (8) for the considered new memory-dependent derivative with time-delay parameter $\omega $ is given by:

(13)

$\left(1+\omega {D}_{\omega}\right)\overrightarrow{q}\left(\overrightarrow{x},t\right)=-k\overrightarrow{\nabla}T\left(\overrightarrow{x},t\right).$
Taking the memory-time derivative of Eq. (9) (suppressing $\overrightarrow{x}$ for convenience), we get:

(14)

$\rho C{D}_{\omega}\dot{T}=-\overrightarrow{\nabla}\cdot {D}_{\omega}q+{D}_{\omega}Q.$
Multiplying Eq. (14) by $\omega $ and adding to Eq. (9), we obtain:

(15)

$\left(1+\omega {D}_{\omega}\right)\rho C\dot{T}=-\overrightarrow{\nabla}\cdot \left(1+\omega {D}_{\omega}\right)q+\left(1+\omega {D}_{\omega}\right)Q.$
Substituting from Eq. (13) into the Eq. (15), we finally get:

(16)

$k{\nabla}^{2}T=\left(1+\omega {D}_{\omega}\right)\left(\rho C\dot{T}-Q\right).$
Eq. (16) is the new bioheat transfer equation with memory- dependent derivative, taking into account the time-delay parameter$\omega $. The Pennes’ bioheat equation is follows as a limit case when $\omega \to 0$. This model is more intuitionistic for understanding the physical meaning and the corresponding memory dependent differential equation is more expressive.

The temperature distribution in human skin with instantaneous surface heating for more understanding to the thermal wave propagation behaviour of fractional bioheat transfer equation is considered. One dimensional model of the FPBE in a finite medium was established. The predictions of the FPBE equation are then compared those of other representations of the Pennes’ equation, to reveal the distinct differences in heat transfer modelling processes.

The one-dimensional form of Eq. (16), given ${Q}_{m}$= constant, is written as:

(17)

$k\frac{{\partial}^{2}T}{\partial {x}^{2}}=\left(1+\omega {D}_{\omega}\right)\left[\rho C\dot{T}-{G}_{B}{C}_{B}({T}_{B}-T)-{Q}_{m}\right].$
Then the initial steady state temperature ${T}_{0}(x,0)$ in skin is derived from Eq. (17) as:

(18)

$k\frac{{\partial}^{2}{T}_{0}}{\partial {x}^{2}}=-\left[{G}_{B}{C}_{B}\left({T}_{B}-{T}_{0}\right)+{Q}_{m}\right].$
Subtracting Eq. (18) from Eq. (17), yields:

(19)

$k\frac{{\partial}^{2}\theta}{\partial {x}^{2}}=\left(1+\omega {D}_{\omega}\right)\left({G}_{B}{C}_{B}+\rho C\frac{\partial}{\partial t}\right)\theta ,$
where the evaluation temperature $\theta (x,t)$ is defined as $\theta (x,t)=T(x,t)-{T}_{0}$.

In numerical solutions, we assume that heat flux approaches zero deep within the tissue $(x-d)$ as expected in biological body.

The homogeneous initial conditions are:

(20)

$\theta \left(x,0\right)=\dot{\theta}\left(x,0\right)=0,0\le x\le d.$
The appropriate boundary conditions are taken as:

(21)

$\theta (0,t)={\theta}_{0},{\left.\frac{\partial \theta}{\partial x}\right|}_{x=d}=0,0\le x\le d,t0.$
From now on, the kernel function form $K(t-\xi )$ can be chosen freely as:

(22)

$K(t-\xi )=1-\frac{2b(t-\xi )}{\omega}+\frac{{a}^{2}(t-\xi {)}^{2}}{{\omega}^{2}}=\left\{\begin{array}{ll}1,& a=b=0,\\ 1-\frac{(t-\xi )}{\omega},& a=0,b=\frac{1}{2},\\ 1-(t-\xi ),& a=0,b=\frac{\omega}{2},\\ {\left[1-\frac{(t-\xi )}{\omega}\right]}^{2},& a=b=1,\end{array}\right.$
where $a$ and $b$ are constants.

Appling the Laplace transform with parameter $p$ defined by the formula:

(23)

$L\left\{h(x,t)\right\}=\stackrel{-}{h}(x,p)=\underset{0}{\overset{\infty}{\int}}h\left(x,t\right){e}^{-pt}dt,Re\left(p\right)0,$
to both sides of Eq. (19) taking into account the initial conditions Eq. (20) and to the boundary conditions Eqs. (21), the equations become:

(24)

$\left({D}^{2}-{\eta}^{2}\right)\stackrel{-}{\theta}(x,p)=0,$
with the boundary conditions:

(25)

$\stackrel{-}{\theta}(0,p)=\frac{{\theta}_{0}}{p},{\left.\frac{d\stackrel{-}{\theta}}{dx}\right|}_{x=d}=0,$
where:

$D=\frac{d}{dx},L\left\{\omega {D}_{\omega}h\right(x,t\left)\right\}=H\left(p\right)\beta \left(p\right),\beta \left(p\right)=(1-{e}^{-p\omega})\left(1-\frac{2b}{\omega p}+\frac{2{a}^{2}}{{\omega}^{2}{p}^{2}}\right)$

$-{e}^{-p\omega}\left({a}^{2}-2b+\frac{2{a}^{2}}{\omega p}\right)=\left\{\begin{array}{ll}\left(1-{e}^{-p\omega}\right),& a=b=0,\\ \left(1-\frac{1-{e}^{-p\omega}}{p\omega}\right),& a=0,b=\frac{1}{2},\\ \left(1-{e}^{-p\omega}\right)\left(1-\frac{1}{p}\right)+c{e}^{-p\omega},& a=0,b=\frac{\omega}{2},\\ \left(1-\frac{2}{p\omega}\right)+\frac{2\left(1-{e}^{-p\omega}\right)}{{p}^{2}{\omega}^{2}},& a=b=1,\end{array}\right.$

$H\left(p\right)=L\left\{\left({G}_{B}{C}_{B}+\rho C\frac{\partial}{\partial t}\right)\theta \right\}=\left({G}_{B}{C}_{B}+\rho Cp\right)\stackrel{-}{\theta}\left(x,p\right),{\eta}^{2}=\left(\frac{{G}_{B}{C}_{B}+\rho Cs}{k}\right)\left[1+\beta \left(p\right)\right].$

$-{e}^{-p\omega}\left({a}^{2}-2b+\frac{2{a}^{2}}{\omega p}\right)=\left\{\begin{array}{ll}\left(1-{e}^{-p\omega}\right),& a=b=0,\\ \left(1-\frac{1-{e}^{-p\omega}}{p\omega}\right),& a=0,b=\frac{1}{2},\\ \left(1-{e}^{-p\omega}\right)\left(1-\frac{1}{p}\right)+c{e}^{-p\omega},& a=0,b=\frac{\omega}{2},\\ \left(1-\frac{2}{p\omega}\right)+\frac{2\left(1-{e}^{-p\omega}\right)}{{p}^{2}{\omega}^{2}},& a=b=1,\end{array}\right.$

$H\left(p\right)=L\left\{\left({G}_{B}{C}_{B}+\rho C\frac{\partial}{\partial t}\right)\theta \right\}=\left({G}_{B}{C}_{B}+\rho Cp\right)\stackrel{-}{\theta}\left(x,p\right),{\eta}^{2}=\left(\frac{{G}_{B}{C}_{B}+\rho Cs}{k}\right)\left[1+\beta \left(p\right)\right].$

Imposing the boundary conditions Eq. (21) upon the heat flow Eq. (24), the exact solution for the temperature distribution $\stackrel{-}{\theta}(x,p)$ in the Laplace transform domain is obtained as:

(26)

$\stackrel{-}{\theta}\left(x,p\right)=\frac{{\theta}_{0}}{p}\frac{\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}\eta \left(d-x\right)}{\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}\eta d}.$
The Laplace transform in the above equation is inverted numerically by Zakian algorithm [10, 11]. Biological skin is a very complex heterogeneous medium, containing many different layers, each having their own unique thermal and blood perfusion properties.

Typical values of thermal properties for skin tissue and other parameters have been chosen as $\rho =$ 1000 kg/m^{3}, $C={C}_{B}=$ 4200 J/kg °C [12]. The distance between the skin and body core is $d=$ 0.01208 m and the surface temperature is taken as ${\theta}_{0}=$ 12 °C.

Fig. 1 shows that the effect of the time-delay parameter $\omega $ on the temperature distribution in living tissue. The time-delay parameter plays an important role in bioheat transfer. The heat energy taken away by the blood is directly proportional to the perfusion rate. Also, the skin temperature is maximum at $x=$ 0 for $\omega =$ 0.09. This research also provides some new points for studying transient heat transfer processes in biological systems.

**Fig. 1.**
Temperature distributions versus skin depth temperature for different values of time-delay $\omega $ and $K(t-\xi )=1-(t-\xi )$ with $k=$ 0.2 W/m °C, ${G}_{B}=$ 8.0 kg/m^{3}s and at $t=$ 250 s

Fig. 2 depicts the effect of the perfusion rate of blood on the temperature distribution in living tissue. The perfusion rate of blood shows an important role in bioheat transfer. Since the skin temperature exceeds the arterial temperature, blood perfusion performs a cooling function. The heat energy extracted by the blood is proportional to the perfusion rate. It is clear from the Fig. 2 that the skin temperature for ${G}_{B}=$ 8.0 kg/m^{3}s is lower than that for ${G}_{B}=$ 4.0 kg/m^{3}s and ${G}_{B}=$0.0 kg/m^{3}s. In this case the skin temperature also attains its maximum value at $x=\text{0}$.

**Fig. 2.**
Temperature distributions versus skin depth temperature for time-delay $\omega =$ 0.0009 and $K(t-\xi )=1-(t-\xi )$ with $k=$ 0.2 W/m °C and at $t=$ 250 s for different blood perfusion rate ${G}_{B}$

The main goal of this work is to introduce a new mathematical model for the Pennes’ bioheat transfer equation memory-dependent derivative involving time-delay parameter $\omega $. In this paper, the memory-dependent derivative is successfully incorporated into a bioheat transfer model. In an attempt to reconcile the novel and classical approaches to bioheat transfer, results of our model were compared with those of the classical and hyperbolic bioheat transfer equations

The Pennes’ bioheat transfer model with memory dependent derivative is established and applied to investigate the thermal response in living tissue with instantaneous surface heating. It is observed that thermal wave propagation may provide realistic prediction of temperature distribution in living tissue. The blood perfusion develops the cooling function to prevent the tissue rising, but does not affect the thermal propagation velocity. In this novel theory, we propose that the time-delay parameter$\omega $ becomes a new indicator of bioheat transfer efficiency in living tissues.

**Wang J. L., Li H. F.**Surpassing the fractional derivative: Concept of the memory-dependent derivative. Computers and Mathematics with Applications, Vol. 62, 2011, p. 1562-1567.**Diethelm K.**Analysis of Fractional Differential Equation: An Application-Oriented Exposition Using Differential Operators of Caputo Type. Springer-Verlag, Berlin, Heidelberg, 2010.**Caputo M.**Linear models of dissipation whose Q is almost frequency independent II. Geophysical Journal of the Royal Astronomical Society, Vol. 13, 1967, p. 529-539.**Mishura Y. S.**Stochastic Calculus for Fractional Brownian Motion and Related Processes. Springer-Verlag, Berlin, Heidelberg, 2008.**Ezzat M. A., El-Karamany A. S., El-Bary A. A.**Generalized thermoelasticity with memory-dependent derivatives involving two temperatures. Mechanics of Advanced Materials and Structures, Vol. 23, 2016, p. 545-553.**Yu Y.-J., Hu W., Tian X.-G.**A novel generalized thermoelasticity model based on memory-dependent derivative. International Journal of Engineering Science, Vol. 811, 2014, p. 123-134.**Lord H. W., Shulman Y. A.**Generalized dynamical theory of thermoelasticity. Journal of the Mechanics and Physics of Solids, Vol. 15, 1967, p. 299-309.**Pennes H. H.**Analysis of tissue and arterial blood temperatures in the resting human forearm. Journal of Applied Physiology, Vol. 1, 93, p. 122-1948.**Ozisik M. N., Tzou D. Y.**On the wave theory in heat conduction. Journal of Heat Transfer, Vol. 116, 1994, p. 526-535.**Zakian V.**Numerical inversions of Laplace transforms. Electronics Letters, Vol. 5, 1969, p. 120-121.**Bachher M., Sarkar N., Lahiri A.**Fractional order thermoelastic interactions in an infinite voids material due to distributed time-dependent heat sources. Meccanica, Vol. 50, 2015, p. 2167-2178.**Liu J., Zhang X., Wang C.**Generalized time delay bioheat equation and preliminary analysis on its wave nature. Chinese Science Bulletin, Vol. 42, 1997, p. 289-292.