A Study on Non-Linear DPL Model for Describing Heat Transfer in Skin Tissue during Hyperthermia Treatment

Abstract

The article studies the simulation-based mathematical modeling of bioheat transfer under the Dirichlet boundary condition. We used complex non-linear dual-phase-lag bioheat transfer (DPLBHT) for analyzing the temperature distribution in skin tissues during hyperthermia treatment of infected cells. The perfusion term, metabolic heat source, and external heat source were the three parts of the volumetric heat source that were used in the model. The non-linear DPLBHT model predicted a more accurate temperature within skin tissues. The finite element Runge–Kutta (4,5) (FERK (4,5)) method, which was based on two techniques, finite difference and Runge–Kutta (4,5), was applied for calculating the result in the case of our typical non-linear problem. The paper studies and presents the non-dimensional unit. Thermal damage of normal tissue was observed near zero during hyperthermia treatment. The effects of the non-dimensional time, non-dimensional space coordinate, location parameter, regional parameter, relaxation and thermalization time, metabolic heat source, associated metabolic heat source parameter, perfusion rate, associated perfusion heat source parameter, and external heat source coefficient on the dimensionless temperature profile were studied in detail during the hyperthermia treatment process.

1. Introduction

Tumors or cancerous cells are a classical sign of inflammation and can be benign or malignant (cancerous). In America, nearly 606,880 people were anticipated to die from cancer in 2019, which translated to about 1660 deaths each day. Cancer is the second most common cause of death in the U.S., exceeded only by heart disease [1]. Therefore, the study of tumor treatment is required to save human lives in the world. Several researchers [2,3,4,5] studied the therapeutic treatment of bioheat transfer in skin tissue with the help of mathematical modeling. Mathematical modeling of heat transfer in biological systems has been a broad field of study for various biologists, physicians, mathematicians, and engineers [6]. An efficient clarification of the physiological relation between the vascular system and tissue is necessary in medical science for treating fatal diseases like tumors. Currently, mathematical models are commonly used to describe the process of hyperthermia, cryosurgery, and many other techniques for the treatment of tumors. It is mandatory to know the thermal effect in skin tissue during the hyperthermia treatment process. The size, shape, and location of tumors are important factors for the treatment process [7].

Several bioheat transfer-based models have assumed the physiological properties of human beings to be constant, which are not described accurately for hyperthermia treatment of tumors or cancer. However, because the inner structure of the human body is inhomogeneous, the physiological parameters depend on local tissue temperature. Some researchers [8,9,10] considered the perfusion term to be a function of the temperature in local tissue. Similarly, some authors [11,12] assumed a metabolic heat source in their model, and this was also a function of the temperature in local tissue. The perfusion term and metabolic heat source are both considered a function of temperature in local tissue, i.e., a realistic-type function, and the external heat source is taken as electromagnetic radiation [13]. However, the location and shape parameters are not derived very well.

Modeling of the tumor treatment is done by the study of the heat transfer in the biological system. The treatment of tumors has been broadly studied in pre-clinical models with human clinical trials [14]. The treatment techniques such as hyperthermia, thermal ablation, cryoablation, and cryosurgery are used for selectively destroying the tumor in living skin tissue. Thermal therapy is an ideal modality for the treatment of infected cells using different types of external heat sources like electromagnetic irradiation [6,10,12,13], magnetic nanoparticles (MNPs) [15,16], etc.

Andreozzi et al. (2019) [17] performed a sensitivity analysis of the hyperthermia effects on a typical transient percolation process in a tumor. In this process, the temperature was raised in a tumor region according to different categories of thermal therapy. Bioheat transfer was analyzed by Andreozzi (2019) [18], who took many mechanisms into account, such as thermal conduction in tissues, convection and blood perfusion, metabolic heat generation, vascular structure, and the change of tissue properties depending on the physiological condition. A numerically investigated bioheat transfer model has been used for hyperthermia treatment with the convection term instead of the perfusion term in the energy conservation equation for tissue and blood [19]. Wang et al. (2015) [20] studied the temperature distribution within biological organs for therapeutic aspects related to hyperthermia treatments such as radiofrequency ablation. The accuracy of temperature-based treatment depends on accurate prediction and control of the temperature in skin tissue [12]. A quantitative analysis of the relationship between arterial blood and tissue temperature was done by Pennes (1948) [21]. There are many bioheat transfer models for studying the heat transfer in skin tissue in the existing literature, and it was found that the commonly used bioheat transfer model for analyzing the temperature distribution is the Pennes bioheat transfer model [21], which is based on the classical constitutive relation that was introduced by Fourier, i.e.,

$$q(r,t)=▿T(r,t).$$

(1)

Penne’s bioheat transfer (PBHT) model predicts the temperature with the infinite speed of propagation, which is incompatible in the real domain. To unify this, consecutively, Cattaneo [22] and Vernotte [23] introduced it in 1958, independently; so, the heat flux and temperature gradient are combined with a constitutive relation, which is given as:

$$q(r,t+{\tau }_q)=▿T(r,t).$$

(2)

Equation (2) is known as the single-phase-lag (SPL) constitutive relation. Relaxation time $\left({\tau }_q\right)$ indicates the lag time due to heat flux. In 1995, Tzou [24,25] introduced his concept in the generalization of the SPL model by assuming thermalization time due to the temperature gradient, called the dual-phase-lag (DPL) constitutive relation, i.e.,

$$q(r,t+{\tau }_q)=▿T(r,t+{\tau }_T).$$

(3)

where ${\tau }_T$ is known as the thermalization time, and the combination of the DPL constitutive relation and energy balance equation is known as the DPLBHT model. For the study of micro-scale responses in time and to capture the micro-scale responses in space, the DPL bioheat transfer model has been used. Therefore, in the existing literature, the DPL bioheat transfer model is the most realistic in comparison with others [26].

Thermal-probes and cryoprobes are used for tumor or cancer treatment. In the DPLBHT model, the thermal correlation between the cylindrical cryoprobe and skin tissues was studied by Mochnaki and Machrzak [27]. A relation was developed between the heat transfer in perfused skin tissue with the thermal-probe and a local symmetric component of the vascular system [28]. Many authors [8,29,30,31,32,33] have assumed a non-linear PBHT model with a physiological property such as perfusion rate for finding the temperature in skin tissues. In reality, These types of models do not give realistic data of the temperature in skin tissue because they do not consider the relaxation and thermalization time in the heat flux and temperature gradient, respectively. This is a drawback of this type of nonlinear bioheat transfer model.

The thermal behaviors of a perfused tissue with two co-current and counter-current vascular networks were investigated numerically under an interstitial hyperthermia process using both local thermal equilibrium (LTE) and local thermal nonequilibrium (LTNE) assumptions [34]. Zhang et al. [30] analyzed the PBHT model under the steady-state condition with the perfusion rate in skin tissues varying linearly, quadratically, and exponentially with local skin tissue temperature and solved it using the boundary reciprocity method. The smoothed particle hydrodynamic method was used for the results of the nonlinear PBHT model with space coordinate-dependent thermal conductivity. Several researchers studied the temperature distribution in skin tissues using the DPL bioheat transfer model, undertaking different types of volumetric heat sources. They solved the DPLBHT model using the finite element Legendre wavelet Galerkin method [15] and used the FERK (4,5) method [10,12,13] and finite difference-decomposition method [6]. The development of the reconfigurable distributed multiple-input multiple-output technique in a practical communication environment was proposed by Do and Haas [35].

In this paper, we propose the highly non-linear DPLBHT equation under the constant boundary condition, which consists of temperature-dependent metabolism and blood perfusion heat generation, as well as a Gaussian heat source. This model is very useful for hyperthermia treatment because the physiological properties of biological skin tissue are considered as a realistic function of local tissue temperature. These types of physiological properties with a Gaussian heat source have been used in the DPLBHT model till now. Due to the lower computational complexity and less data storage, combined with the high accuracy, the problem is converted into a system of ordinary differential equations with initial conditions using the finite difference scheme. This system of ordinary differential equations with initial conditions is solved using the RK (4,5) scheme. All parameters such as the location parameter, regional parameter, and relaxation and thermalization time provide a better understanding of the control temperature in the hyperthermia condition. The metabolic heat source, associated metabolic heat source, and external heat source increase as the local skin tissue temperature increases. The perfusion rate, associated with the perfusion heat source parameters, decreases as the local skin tissue temperature increases as well.

This paper is organized into seven sections. In the first section, the introduction of the bioheat transfer models, some methods, and also the nomenclature are given, which support our proposed work. Our mathematical problem is formulated in Section 2. In the third section, our mathematically formulated problem is converted into a dimensionless form. Section 4 describes the solution of the proposed problem using the FERK (4,5) method. In Section 5, we propose the exact solution of our problem for a particular case to verify the FERK (4,5) method. The results and discussion are given in Section 6. The last section consists of the conclusions of the proposed work.

2. Formulation of the Problem

Hyperthermia is a treatment process of tumors and cancer. In this treatment process, the temperature of the tumor region is kept between 41 °C and 46 °C, with approximately a time period of 15 to 60 min [36]. The outer surface of the skin tissue is kept at a fixed temperature $T_0$ = 37 °C, initially, during hyperthermia treatment. The outer surface is heated with a Gaussian heat source externally. The inner surface of the skin, i.e., $r=0$, is insulated, and the temperature of the outer surface is maintained constant with the help of a cooling pad, which is shown in Figure 1. The tumor region is indicated by the schematic geometry of the skin tissue, which depends on the probe region parameter $r_p$. If this parameter increases or decreases, then the location of the tumor in skin tissue changes. This is seen in the external heat source. We used a heated metal disc with temperature control and approximated by the one-dimensional non-linear DPLBHT model in the Cartesian coordinate system under the first kind (constant) boundary condition. This is a combination of the DPL constitutive relation and also the one-dimensional energy balance equation. The energy balance equation in one-dimensional form is written as [21]:

$$\rho c\frac{\partial T(r,t)}{\partial t}=-\nabla q(r,t)+Q_b+Q_m+Q_r,$$

(4)

where the left-hand side denotes the conduction term in the skin tissue and the first term of the right-hand side denotes the convection term in the skin tissue. $Q_b$, $Q_m$, and $Q_r$ are the perfusion heat source, the metabolic heat source, and the externally applied heat source term, which is taken as a Gaussian-type heat source. The metabolic heat source is generated in the body by the intake of food, and the perfusion heat rate is the heat source that is spent in blood circulation. Still, the Gaussian heat source is externally applied on the outer surface.

The blood perfusion $Q_b$ indicates convection in the blood. This term removes heat due to the flow of blood. It was defined by several researchers [10,13,30].

$$Q_b=w_b\left(T\right){\rho }_b{c}_b(T_b-T),$$

(5)

The perfusion coefficient $w_b\left(T\right)$ depends on the tissue temperature due to the anatomical structure of the skin tissue containing blood vessels. $w_b\left(T\right)$ was taken from Zhang et al. [30]:

$$w_b\left(T\right)=w_{b0}\times e^{a\left(\frac{T-T_0}{T_0}\right)},$$

(6)

where $w_{b0}$ and a are the assumed blood perfusion coefficient and associated blood perfusion coefficient initially.

In Equation (4), $Q_m$ indicates the metabolism of the human body. It is the factor that increases the local tissue temperature. The realistic function $Q_m$ of the local tissue temperature was taken from Mitchell et al. [11]:

$$Q_m=Q_{m0}\times 2^{\beta \left(\frac{T-T_0}{10}\right)},$$

(7)

where $Q_{m0}$ and $\beta$ are the initial metabolic heat source coefficient and the associated metabolic constant, respectively. T and $T_0$ are the local skin temperature and the initial temperature of the skin tissue, respectively.

A Gaussian heat source-type expression is curve fitted based on the experimental measurements of the specific absorption rate distribution in the target region [37]. The selection of a Gaussian distribution source termed as a spatial heating source helps to determine the hyperthermia position. Therefore, that Gaussian-type heat source expression is applied on the outer surface, and it was taken from Kumar et al. [13]:

$$Q_r=Q_{r0}{e}^{\eta {(\overline{r}-r_p)}^2},$$

(8)

where $Q_{r0}$ is the reference value of external heat generation in the tissue; $\overline{r}(=L-r)$ is the measured length from the outer surface; L is the length of the tissue temperature. $r_p$ is the length of the probe region.

Now, the energy balanced equation is combined with the approximation of the DPL constitutive relation of the first order in one-dimensional form, then we obtain the non-linear DPLBHT equation [13], i.e.,

$$\begin{array}{c}\hfill \rho c{\tau }_q\frac{{\partial }^{2}T}{\partial t^2}+\rho c\frac{\partial T}{\partial t}=k\frac{{\partial }^{2}T}{\partial r^2}+k{\tau }_T\frac{{\partial }^{3}T}{\partial t\partial r^2}+Q_b\\ \hfill +Q_m+Q_r+{\tau }_q\frac{\partial Q_b}{\partial t}+{\tau }_q\frac{\partial Q_m}{\partial t}+{\tau }_q\frac{\partial Q_r}{\partial t},\end{array}$$

(9)

initially subjected to:

$$T(r,0)=T_0,\frac{\partial T(r,0)}{\partial t}=0.$$

(10)

Physically, the heating/cooling condition at the outer surface $r=L$ was given in [15], i.e.,

$$T(L,t)=T_w,$$

(11)

and the inner surface is adiabatic [10,12], so that it is defined as:

$$\frac{\partial T(0,t)}{\partial r}=0.$$

(12)

In the existing literature, the DPLBHT model with a Gaussian heat source in the presence of a realistic function of the metabolic and perfusion heat generation terms has not been studied till now. We considered the Gaussian heat source because it helps with the control of temperature during hyperthermia treatment.

3. Conversion of the Problem into Dimensionless Form

One way of breaking away from the quantitative features of the corresponding model is to rewrite the equation in terms of non-dimensional quantities. The proposed DPLBHT equations reduces in terms of non-dimensional quantities. These types of studies allow an easy and direct comparison of the proposed mathematical model. Therefore, the non-dimensionless variables with similarity criteria was introduced in [13]:

$$\begin{array}{c}\hfill x=\frac{r}{L},F_o=\frac{kt}{\rho c{L}^2},F_{oq}=\frac{k{\tau }_q}{\rho c{L}^2},F_{oT}=\frac{k{\tau }_T}{\rho c{L}^2},\\ \hfill \theta =\frac{T-T_0}{T_0},P_f=\sqrt{\frac{w_{b0}{c}_b{\rho }_b}{k}}L,P_m=\frac{Q_{m0}{L}^2}{k{T}_0},\\ \hfill P_r=\frac{Q_r}{k{T}_0}{L}^2,\alpha =0.1\times T_0\times \beta ,{\eta }_1=\eta \times L,\\ \hfill \overline{x_p}=1-x_p.\end{array}$$

(13)

Upon using Equation (13), then Equations (9)–(12) can be reduced as follows:

$$\begin{array}{c}\hfill F_{oq}\frac{{\partial }^2\theta }{\partial F_o^2}+\Theta \frac{\partial \theta }{\partial F_o}=\frac{{\partial }^2\theta }{\partial x^2}+F_{oT}\frac{{\partial }^3\theta }{\partial F_o\partial x^2}+P_f^2({\theta }_b-\theta )e^{a\theta }\\ \hfill +P_m\times 2^{\alpha \theta }+P_{r}exp\left({\eta }_1{(\overline{x}-x_p)}^2\right),\end{array}$$

(14)

where $\Theta =\left[1-(a({\theta }_b-\theta )-1)F_{oq}{P}_f^2{e}^{a\theta }-F_{oq}{P}_m\alpha log\left(2\right)·2^{\alpha \theta }\right]$, initially subjected to:

$$\theta (x,0)=0\mathrm{and}\frac{\partial \theta (x,0)}{\partial F_o}=0.$$

(15)

The outer boundary condition:

$$\theta (1,Fo)={\theta }_w,$$

(16)

and the inner boundary condition is taken as insulated:

$$\frac{\partial \theta (0,F_o)}{\partial x}=0.$$

(17)

Our problem transformed into dimensionless form is solved using the finite element Runge–Kutta(4,5) method. The detail of this method is described in the following section.

4. Finite Element Runge–Kutta (4,5) Method

The FERK (4,5) method consists of two techniques: (i) discretization in space coordinates using the finite difference scheme [38]; then our proposed problem is reduced to a system of ordinary differential equations (ODEs); and (ii) another scheme, RK (4,5) [39], is used for the solution of the ODEs. These are described in the following subsection.

4.1. Discretization Technique in Space Coordinates

We discretize the finite range of space coordinates [0,1] into $l+1$ discrete points, i.e., $0=x_0<x_1<x_2<x_3<\cdots <x_i<\cdots <x_l<x_{l+1}=1$, where $x_{i+1}=x_i+h$, and h is the constant length of the subinterval of domain $\left[\mathrm{0,1}\right]$, which is equal between any two nodes. The approximation of the second order derivative is defined as:

$$\frac{{\partial }^2\theta (x,F_o)}{\partial x^2}{|}_{x=x_i}=\frac{{\theta }_{i+1}\left(F_o\right)-2{\theta }_i\left(F_o\right)+{\theta }_{i-1}\left(F_o\right)}{h^2},1\le i\le l.$$

Using the central finite difference formula, Equations (14)–(17) are reduced to the following form:

$$\begin{array}{c}\hfill F{o}_q\frac{d^2{\theta }_1}{dF{o}^2}=-{\Theta }_1\frac{d{\theta }_1}{d{F}_o}+\frac{\left(-29{\theta }_1+38{\theta }_2-9{\theta }_3\right)}{21h^2}\\ \hfill +\frac{F_{oT}}{h^2}\frac{d}{d{F}_o}\left(-29{\theta }_1+38{\theta }_2-9{\theta }_3\right)+P_f^2({\theta }_b-{\theta }_1)\times e^{a{\theta }_1}\\ \hfill +P_m\times 2^{\alpha {\theta }_1}+P_r\times exp\left({\eta }_1({\overline{x}}_1-x_p)\right),\end{array}$$

(18)

$$\begin{array}{c}\hfill F{o}_q\frac{d^2{\theta }_i}{dF{o}^2}=-{\Theta }_i\frac{d{\theta }_i}{d{F}_o}+\frac{\left({\theta }_{i+1}-2{\theta }_i+{\theta }_{i-1}\right)}{h^2}\\ \hfill +\frac{F_{oT}}{21h^2}\frac{d}{d{F}_o}\left({\theta }_{i+1}-2{\theta }_i+{\theta }_{i-1}\right)+P_f^2({\theta }_b-{\theta }_i)\times e^{a{\theta }_i}\\ \hfill +P_m\times 2^{\alpha {\theta }_i}+P_r\times exp\left({\eta }_1({\overline{x}}_i-x_p)\right),1<i<n,\end{array}$$

(19)

$$\begin{array}{c}\hfill F{o}_q\frac{d^2{\theta }_n}{dF{o}^2}=-{\Theta }_n\frac{d{\theta }_n}{d{F}_o}+\frac{\left({\theta }_w-2{\theta }_n+{\theta }_{n-1}\right)}{h^2}\\ \hfill +\frac{F_{oT}}{h^2}\frac{d}{d{F}_o}\left({\theta }_w-2{\theta }_n+{\theta }_{n-1}\right)+P_f^2({\theta }_b-{\theta }_n)\times e^{a{\theta }_n}\\ \hfill +P_m\times 2^{\alpha {\theta }_n}+P_r\times exp\left({\eta }_1({\overline{x}}_n-x_p)\right).\end{array}$$

(20)

with initial conditions:

$${\theta }_i\left(0\right)=0\mathrm{and}\frac{d{\theta }_i\left(0\right)}{d{F}_o}=0,0\le i\le n,$$

(21)

where ${\Theta }_i=\left[1-(a({\theta }_b-{\theta }_i)-1)F_{oq}{P}_f^2{e}^{a{\theta }_i}-\alpha F_{oq}{P}_{m}log\left(2\right)\times 2^{\alpha {\theta }_i}\right],1\le i\le n.$

4.2. Runge–Kutta (4,5) Scheme

Let:

$$\frac{d}{dFo}\left(F{o}_q\frac{d{\theta }_i\left(F_o\right)}{dFo}\right)=\frac{d{\Psi }_i\left(F_o\right)}{dFo}.$$

(22)

Thus, Equations (18)–(21) are reduced to the following form:

$$\frac{d{\theta }_i\left(F_o\right)}{dFo}=\frac{{\Psi }_i\left(F_o\right)}{F{o}_q},1\le i\le n,$$

(23)

$$\begin{array}{c}\hfill \frac{d{\Psi }_1\left(F_o\right)}{d{F}_o}=-{\theta }_1\frac{{\Psi }_1\left(F_o\right)}{F_{oq}}+\frac{\left(-29{\theta }_1+38{\theta }_2-9{\theta }_3\right)}{21h^2}\\ \hfill +\frac{F_{oT}}{21h^2{F}_{oq}}\left(-29{\Psi }_1+38{\Psi }_2-9{\Psi }_3\right)+P_f^2({\theta }_b-{\theta }_1)\times e^{a{\theta }_1}\\ \hfill +P_m\times 2^{\alpha {\theta }_1}+P_r\times exp\left({\eta }_1({\overline{x}}_1-x_p)\right),\end{array}$$

(24)

$$\begin{array}{c}\hfill \frac{d{\Psi }_i\left(F_o\right)}{dFo}=-{\Theta }_i\frac{{\Psi }_i\left(F_o\right)}{F_{oq}}+\frac{\left({\theta }_{i+1}-2{\theta }_i+{\theta }_{i-1}\right)}{h^2}\\ \hfill +\frac{F_{oT}}{h^2{F}_{oq}}\left({\Psi }_{i+1}-2{\Psi }_i+{\Psi }_{i-1}\right)+P_f^2({\theta }_b-{\theta }_i)\times e^{a{\theta }_i}\\ \hfill +P_m\times 2^{\alpha {\theta }_i}+P_r\times exp\left({\eta }_1({\overline{x}}_i-x_p)\right),\\ \hfill 1<i<n,\end{array}$$

(25)

$$\begin{array}{c}\hfill \frac{d{\Psi }_n\left(F_o\right)}{dFo}=-{\theta }_n\frac{{\Psi }_n\left(F_o\right)}{F{o}_q}+\frac{\left({\theta }_w-2{\theta }_n+{\theta }_{n-1}\right)}{h^2}\\ \hfill +\frac{F_{oT}}{h^2{F}_{oq}}\left(-2{\Psi }_n+{\Psi }_{n-1}\right)+P_f^2({\theta }_b-{\theta }_n)\times e^{a{\theta }_n}\\ \hfill +P_m\times 2^{\alpha {\theta }_n}+P_r\times exp\left({\eta }_1({\overline{x}}_n-x_p)\right),\end{array}$$

(26)

with initial conditions:

$$\begin{array}{c}\hfill {\theta }_i\left(0\right)=0,\mathrm{and}\end{array}$$

(27)

$$\begin{array}{c}\hfill {\Psi }_i\left(0\right)=0,\end{array}$$

(28)

where ${\Theta }_i=\left[1-(a({\theta }_b-{\theta }_i)-1)F_{oq}{P}_f^2{e}^{a{\theta }_i}-\alpha F_{oq}{P}_{m}log\left(2\right).2^{\alpha {\theta }_i}\right],$$i=1,2,3,\cdots ,n.$

The typical type of ODEs are computed with the help of the RK (4,5) technique. The proposed technique is very efficient for calculating the non-dimensional temperature of our problem. MATLAB-2018 and DEV-C++ software were used for all computational work.

5. Exact Solution

In order to justify the validity of the results, the exact solution is required. In the present non-linear DPLBHT model in Equation (13), when we assumed the associated blood perfusion constant $a=0$, associated metabolism parameter $\alpha =0$, and external heat source $P_r=0$, then this model was reduced to the following form:

$$\begin{array}{c}\hfill F_{oq}\frac{{\partial }^2\theta }{\partial F_o^2}+\left[1+F_{oq}{P}_f^2\right]\frac{\partial \theta }{\partial F_o}=\frac{{\partial }^2\theta }{\partial x^2}+F_{oT}\frac{{\partial }^3\theta }{\partial F_o\partial x^2}+P_f^2({\theta }_b-\theta )+P_m,\end{array}$$

(29)

subject to the initial conditions:

$$\theta (x,0)=0\mathrm{and}\frac{\partial \theta (x,0)}{\partial F_o}=0.$$

(30)

boundary conditions:

$$\theta (1,Fo)={\theta }_w,$$

(31)

and symmetric conditions:

$$\frac{\partial \theta (0,F_o)}{\partial x}=0.$$

(32)

Taking the Laplace transform and then its inversion technique, the solution of Equation (29) under the initial condition (30) and boundary conditions (31) and (32) turns out to be:

$$\begin{array}{l}\theta \left(x,F_o\right)={\theta }_w\left(\frac{cosh\left(P_{f}x\right)}{cosh\left(P_f\right)}+{\displaystyle \sum _{n=1}^{\infty }}\frac{{\mathrm{e}}^{S_{n1}Fo}cosh\left(M_{Sn1}x\right)}{S_{n1}\left(S_{n1}-S_{n2}\right)}+{\displaystyle \sum _{n=1}^{\infty }}\frac{{\mathrm{e}}^{S_{n2}{F}_o}cosh\left(M_{Sn2}x\right)}{S_{n2}\left(S_{n2}-S_{n1}\right)}\right)\hfill \\ +\left({P_f}^2{\theta }_b+P_m\right)\left({P_f}^{-2}-\frac{cosh\left(P_{f}x\right)}{{P_f}^{2}cosh\left(P_f\right)}-{\displaystyle \sum _{n=1}^{\infty }}\frac{{\mathrm{e}}^{S_{n1}{F}_o}cosh\left(M_{Sn1}x\right)}{S_{n1}\left(F_{oq}{S_{n1}}^2+\left(1+{P_f}^2{F}_{oq}\right)S_{n1}+{P_f}^2\right)\left(S_{n1}-S_{n2}\right)}\right)\hfill \\ -\left({P_f}^2{\theta }_b+P_m\right)\left({\displaystyle \sum _{n=1}^{\infty }}\frac{{\mathrm{e}}^{S_{n2}{F}_o}cosh\left(M_{Sn2}x\right)}{S_{n2}\left(F_{oq}{S_{n2}}^2+\left(1+{P_f}^2{F}_{oq}\right)S_{n2}+{P_f}^2\right)\left(S_{n2}-S_{n1}\right)}\right),n=1,2,3,\cdots ,\hfill \end{array}$$

(33)

where $r_n=1+{P_f}^2{F}_{oq}+{\left(2n-1/2\right)}^2{\pi }^2{F}_{oT},$$S_{n1}=\frac{-r_n+\sqrt{{r_n}^2-4F_{oq}\left({P_f}^2+{\left(2n-1/2\right)}^2{\pi }^2\right)}}{2F_{oq}},$$S_{n2}=\frac{-r_n-\sqrt{{r_n}^2-4F_{oq}\left({P_f}^2+{\left(2n-1/2\right)}^2{\pi }^2\right)}}{2F_{\mathit{oq}}},$$M_{Sn1}=\frac{F_{oq}{S_{n1}}^2+\left(1+{P_f}^2{F}_{oq}\right)S_{n1}+{P_f}^2}{1+F_{oT}{S}_{n1}}$ and $M_{Sn2}=\frac{F_{oq}{S_{n2}}^2+\left(1+{P_f}^2{F}_{oq}\right)S_{n2}+{P_f}^2}{1+F_{oT}{S}_{\mathit{n}\mathit{2}}}.$

6. Results and Discussion

In the proposed computational work, the temperature profile in the skin tissue was computed from the highly non-linear DPLBHT model whenever the outer boundary was kept at a fixed initial temperature. The proposed mathematical model consisted of a temperature-dependent perfusion term and also a temperature-dependent metabolic heat source, the physical function of both terms being experimentally validated. The results are shown graphically in Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11. Those parameters whose values differed from the reference values of the non-dimensional parameters to calculate the non-dimensional temperature profile in the skin tissue with a finite length are given in the following

Table 1. Physiological properties of skin tissues.

Parameters	Units	Numerical Values	References
$\rho$	kg/m${}^3$	1000	[6]
c	J/kg °C	4000	[32]
$c_b$	J/kg ${}^{\circ }$C	3340	[6]
k	W/m °C	0.5	[34]
L	m	0.05	[10]
$T_0$	°C	37	[31]
$T_w$	°C	37	[31]
${\tau }_q$	s	28	[10]
${\tau }_T$	s	20	[10]
$T_b$	°C	37	[32]

,

Table 2. Metabolic properties of skin tissues.

Parameters	Units	Numerical Values	References
$Q_{m0}$	W/m${}^3$	0.17	[11]
$\beta$	°C ${}^{-1}$	2.15	[13]

,

Table 3. Blood perfusion properties.

Parameters	Units	Numerical Values	References
$w_{b0}$	s${}^{-1}$	$3.075\times {10}^{-3}$	[13]
a	–	2.15	[13]

,

Table 4. External heat source properties of the targeted region.

Parameters	Units	Numerical Values	References
$Q_{r0}$	W/m${}^3$	$5.17\times {10}^5$	[13]
$\eta$	m${}^{-1}$	−127	[10]
$r_p$	m	0.005	[10]

and

Table 5. Numerical value of the parameter used in the thermal damage function.

Parameters	Units	Numerical Values	References
A	s${}^{-1}$	$3.11\times {10}^{98}$	[15,40]
$\delta E$	J mol${}^{-1}$	$6.75\times {10}^5$	[15,40]
R	J mol${}^{-1}$ K${}^{-1}$	$8.314$	[15,40]

.

6.1. Comparison of the Exact and Numerical Solutions

The FERK (4,5) method was applied to find the numerical results of the proposed highly non-linear DPLBHT model. The proposed technique unified the essence of the RK (4,5) method with more efficiency and less local error [38,39]. The accurate feasibility of the present numerical scheme was shown by comparing it with the exact solution under a particular case of the non-linear DPL bioheat transfer equation in skin tissues. The comparison of the FERK (4,5) method with the exact solution is shown in Figure 2, and we observed high accuracy with less computational complexity, being merits of the proposed method.

6.2. Effect of $F_o$ and x

Figure 3 shows the temperature profile with respect to the dimensionless space coordinates and also with dimensionless time. Figure 3a is considered a targeted point at $x=0.5$. The temperature profile was approximately 0.24 ≈ 45.88 ${}^{\circ }$C, where the hyperthermia temperature lied between 41 °C and 46 °C for 15 to 60 min. Therefore, the infected or tumorous cells died for different values of $F_o$. The graphs are drawn for $x=0.9$ when dimensionless time $F_o=0.01$. Then, the temperature profile was less and $F_o=0.015$, and then, the temperature profile reached the hyperthermia temperature. Figure 3b demonstrates the temperature profile along $F_o$ for different values of the non-dimensional distance coordinates. The region $x=0.8\mathrm{and}0.9$ reached the hyperthermia region. In this case, $x=0.9$ was the targeted region, so that the temperature at that position was $0.24\approx 45.88$. Therefore, the temperature in the skin tissue achieved the hyperthermia temperature (41 ${}^{\circ }$C to 46 ${}^{\circ }$C) in the tumor region.

6.3. Effect of $x_p$

Figure 4 demonstrates three-dimensional graphs such as $a,b,c,d,\mathrm{and}e$. These graphs were drawn for different values of location parameter $x_p$. This shows that if $x_p$ were changed, then the location of the tumor also changed in the skin tissue. Figure 4a–e are drawn for different values of $x_p=0.1,0.3,0.5,0.7,\mathrm{and}0.9.$ From this parameter, we heated the appropriate targeted point. In a similar way, we also clearly see in Figure 5 that the position of parameter $x_p=0.1,0.3,0.5,0.7,0.9$ was heated with the hyperthermia temperature approximately in the nbd.of 0.25 ≈ 46 ${}^{\circ }$C. From Figure 4 and Figure 5, we see that the position of hyperthermia depended on the location parameter $x_p$.

6.4. Effect of $\eta$ in the Contour Plot

Figure 6 shows the effect of the regional parameter $\eta$ when the values of $\eta$ decreased, then the width of infected or cancerous cells decreased. From regional parameter, we confirmed accordingly the width of tumorous cells (see Figure 6). We chose the value of $\eta$ according to the width of the tumor area. The regional parameter helped in the accurate heating of the cancerous cells, and the region was controlled by it. The regional parameter was very important for the hypothermia treatment of the tumor.

6.5. Effect of $F{o}_q$ and $F{o}_T$

In Figure 7, the effect of the relaxation time in heat flux is shown. The value of the temperature in the local tissue increased as dimensionless time increased from zero to 0.03 and, after that, constant at the targeted region at 0.9. Therefore, we can say that a significant effect of $F{o}_q$ was shown in the targeted region. In Figure 8, we show the effect of lag time $F{o}_T$ due to the temperature gradient. The value of the temperature profile decreased with the lag time $F{o}_T$ in between dimensionless time zero to 0.3 at $x=0.9$. However, in comparison to $F{o}_T$, $F{o}_q$ was more effective in the targeted region. Therefore, $F{o}_T$ and $F{o}_q$ in the hyperthermia position played a vital role because the internal cells were very sensitive, and it was very effective with a small amount of temperature increase.

6.6. Effect of $Pm$ and $\alpha$

The effects of the metabolic heat source and its associated parameter are shown in Figure 9. In Figure 9a, we find that the value of the temperature profile increased with the metabolic heat source, and also, we observe that $P_m$ increased, then the temperature profile was constant. Similarly, the value of the associated metabolic coefficient increased as the temperature profile remained constant. Therefore, both $Pm$ and $\alpha$ affected the hyperthermia temperature.

6.7. Effect of $P_f$

Figure 10 shows the effect of the blood perfusion coefficient and the associated blood perfusion parameter. The value of $P_f$ decreased as the temperature profile increased (seen in Figure 8a), and the value of $\alpha$ increased as the temperature profile decreased. The effects of $P_f$ and $\beta$ were both meaningful to maintain the hyperthermia temperature at the hyperthermia position.

6.8. Effect of $P_r$ and $\eta$

In Figure 11, we see the effect of the external heat source, as well as the cancerous cells’ region parameter. We observed the value of $P_r$ and $\eta$ to increase as the temperature profile increased (seen in Figure 9a,b). The external heat source in the cancerous cell region, the location parameter, was present in the external heat source, which was a type of Gaussian heat source. Therefore, all these parameters had a more significant effect at the position of hyperthermia.

6.9. Effect of Damage Integral Function $(\Omega )$

We can see that the temperature of the skin tissue increased rapidly along with the increase in the value of reference heat source $Q_{r0}$. Therefore, the thermal damage in normal tissue may have occurred by the heating with the external heat source. The thermal damage function $\Omega$ was established by Henriques and Moritz [41], which is given as follows:

$$\Omega \left(t\right)={\int }_0^{t}A\exp \left(\frac{-\delta E}{R(T+273)}\right)ds,$$

(34)

where $\delta E$ is the energy of the initiation of irreversible ablation; A is the frequency constant; R is the universal gas constant; and $T(=T_0+T_0\times \theta )$ is the local tissue temperature, which is taken from Equation (13).

Jiang et al. [40] solved Equation (34) using the difference method, and it is given as follows:

$$\Omega \left(i\right)=\sum _{i=0}^{n}A\exp \left(\frac{-\delta E}{R(T_0(1+\theta \left(i\right))+273)}\right)t,$$

(35)

where i indicates the time step when the node temperature was above the threshold; and n indicates the final time step.

Xu et al. [42] studied a formula of the probability of the thermal damage of normal tissue surrounding the tumor region, which is given as:

$$P\left(i\right)=1-e^{-\Omega \left(i\right)},$$

(36)

where P is the probability of the thermal damage of normal tissue surrounding the tumor region.

The result was calculated when the value of the time step was 0.001 s [40], which is presented in

Table 6. Thermal damage of the normal tissue during the hyperthermia treatment.

Value of ${\mathit{Q}}_{\mathit{r}0}$	Thermal Damage of Normal Tissue	Approximate Value of P
$5.17\times {10}^5$	$2.4985\times {10}^{-10}$	0
$5.17\times {10}^5$	$2.8665\times {10}^{-11}$	0
$5.17\times {10}^5$	$3.2429\times {10}^{-12}$	0

. The numerical result showed that the thermal damage of normal tissue in the target region increased as the value of the external heat source increased. We also observe from Equation (35) that when the value of thermal damage was in the neighborhood of zero, then the probability of the normal tissue damage was near zero, i.e., no thermal damage of the normal tissue during the hyperthermia treatment.

Due to the above description, we can say that the method is very beneficial for medical science and will benefit the user.

7. Conclusions

A highly non-linear DPLBHT model was analyzed under the constant boundary condition in the presence of perfusion and metabolic heat sources, which were a realistic function of temperature. The external heat source was considered as a Gaussian heat source-type function, which consisted of a region of cancerous cells as a parameter, the location parameter. The results of this study are given as follows:

• The finite element RK(4,5) method was applied for the solution of the highly non-linear DPLBHT model and showed high accuracy with less computational complexity.
• The effect of the temperature profile on the different points was analyzed for a different times for the hyperthermia time.
• The effect of regional parameter $\eta$ was that when the values of $\eta$ decreased, then the width of the infected or cancerous cells decreased. The value of $P_r$ and $\eta$ increased as the temperature profile increased.
• We observed that the relaxation time $F{o}_q$ was more effective in the targeted region in comparison with the thermalization time $F{o}_T$.
• As the value of the temperature profile increased, the metabolic heat source increased, but the perfusion rate decreased.
• Graphs were drawn for different values of location parameter $r_p$, then the situation of the targeted region was identified with the help of this parameter.
• We also calculated the probability of the thermal damage of the normal tissue surrounding the tumor and found no thermal damage of the normal tissue during hyperthermia treatment.

From the above observations, we concluded that the presented highly non-linear DPLBHT model played a vital role in the hyperthermia treatment of infected cells.