A Fibonacci Wavelet Method for Solving Dual-Phase-Lag Heat Transfer Model in Multi-Layer Skin Tissue during Hyperthermia Treatment

: In this article, a novel wavelet collocation method based on Fibonacci wavelets is proposed to solve the dual-phase-lag (DPL) bioheat transfer model in multilayer skin tissues during hyperthermia treatment. Firstly, the Fibonacci polynomials and the corresponding wavelets along with their fundamental properties are brieﬂy studied. Secondly, the operational matrices of integration for the Fibonacci wavelets are built by following the celebrated approach of Chen and Haiso. Thirdly, the proposed method is utilized to reduce the underlying DPL model into a system of algebraic equations, which has been solved using the Newton iteration method. Towards the culmination, the effect of different parameters including the tissue-wall temperature, time-lag due to heat ﬂux, time-lag due to temperature gradient, blood perfusion, metabolic heat generation, heat loss due to diffusion of water, and boundary conditions of various kinds on multilayer skin tissues during hyperthermia treatment are brieﬂy presented and all the outcomes are portrayed graphically.


Introduction
Cancer is the leading cause of death in the human population with an annual death rate of 10 million approximately. Consequently, ample research work has been devoted to study the cause and treatment of cancer in recent years [1]. In particular, hyperthermia treatment has proved to be very important and successful for cancer treatment [2]. During hyperthermia treatment, the temperature at the tumor site is increased by certain external heat sources including ultrasound waves, hot-water tubing, infrared radiation, radio-frequency, and several others [3,4]. As is well-known, the affected region (tumor) is always surrounded by healthy tissue, therefore, accurate information of the temperature distribution in the entire affected region plays a vital role during the treatment process. Recently, Romano et al. [5] investigated the effects of temperature and its potential role in pars plana vitrectomy and found that the variations in temperature during vitreoretinal surgery are clinically significant, as the rheology of tamponades can be better manipulated by modulating intraocular pressure and temperature. Moreover, they showed that rapid circulation of fluid in the vitreous cavity reduces the heat produced by the retinal and choroidal surface, bringing the temperature toward room temperature. However, the situation in multilayer skin tissue is even more complicated as the process involves different mechanisms such as heat conduction, water diffusion, vaporization, convection, and blood perfusion. Therefore, the proper modeling and analysis of treatment play an indispensable role in optimizing the temperature distribution in the affected region. This lead to the birth of several heat transfer models including the Wulff model [6], the Klinger model [7], the Cao-Yang-Srivastava model [8], the Chen-Holmes model [9], the Weinbaum-Jiji model [10], the Nakayama-Kuwahara model [11], and the Pennes bioheat transfer model [12].
Among all of the heat transfer models, the Pennes bioheat transfer model has received considerable attention mainly due to its simple and lucid nature, and is governed by Fourier's law of heat conduction q(r, t) = −k∇T(r, t), where ∇T(r, t) and q denote the temperature gradient and heat flux at position r, respectively. The model (1) is somewhat unrealistic due to the assumption that the thermal disturbance propagates at an infinite speed. To overcome this infeasibility, Cattaneo and Vernotte [13,14] proposed a single-phase-lag (SPL) model by adding a phase-lag term to the time-variable t as q(r, t + τ q ) = −k∇T(r, t), where τ q is the thermal relaxation time between heat flux and temperature gradient, which captures the microscale responses in time due to heat flux. One of the major drawbacks of the SPL model (2) is that it captures only the heat flux microstructural response but fails to capture the temperature gradient microstructural response. To overcome this issue, Tzou [15] proposed a dual-phase-lag model by introducing another phase-lag term τ t to the space variable of the SPL model (2) to get a substantial bioheat transfer model of the following form: q(r, t + τ q ) = −k∇T(r, t + τ t ), where the parameters k, τ q , τ t , r, t, and T represent the thermal conductivity of the tissue, phase-lag due to heat flux, phase-lag due to temperature gradient, space coordinate, time, and tissue temperature, respectively. The significance of this model is that the heat propagates with a finite speed with two different relaxation times: one due to heat flux and the other one due to temperature gradient. Since its origin, this model has been proved to be of substantial importance for acquiring an insight into the mechanism underlying the process of heat transfer in living tissues. Finding the analytical and numerical solutions of the DPL model (3) is one of the prime research areas in mathematical, biological, and thermal sciences. Some of the commonly employed numerical techniques for obtaining the solution of model (3) include the finite volume method [16], the Laplace transform method [17], the boundary element method [18], the finite difference-decomposition method [19], the Adomian decomposition method [20], and the B-spline method [21]. Some of the compelling alternatives of the aforementioned numerical techniques include the wavelet-based numerical methods. Due to their simple procedure, easy computation, and rapid convergence, these methods have achieved a significant place in numerical analysis and approximation theory [22,23]. Some of the commonly used wavelet families for solving different physical, engineering, and biological problems include Haar wavelets, Legendre wavelets, Laguerre wavelets, Chebyshev wavelets, harmonic wavelets, Euler wavelets, Bernoulli wavelets, and ultraspherical wavelets. A recent addition to the class of wavelet families is the Fibonacci wavelets, which have attained considerable attention from researchers working across various disciplines in physical and engineering sciences, mainly due to their several distinguishing features [24,25]. For instance, these wavelets are not based on orthogonal polynomials, however, one can differentiate and integrate these wavelets to obtain the corresponding operational matrices. In addition to this, the Fibonacci wavelets are generated from the Fibonacci polynomials, which contain fewer terms in com-parison to the shifted Legendre polynomials, which in turn speeds up the computational process and minimizes the computational errors [26].
With the pleasant characteristics of the Fibonacci wavelets in hindsight, we are profoundly motivated to introduce a new wavelet method based on nonorthogonal Fibonacci polynomials to solve the dual-phase-lag model of the bioheat transfer DPL model in multilayer skin tissue stratified as epidermis, dermis, and subcutaneous under the general boundary conditions. Moreover, we studied the effect of various parameters such as blood perfusion, tissue wall temperature, metabolic heat generation, and heat loss due to diffusion of water in multilayer skin tissue at the location of tumor. The obtained results shows that the proposed method is an appropriate tool for the approximation of smooth and piecewise smooth functions. Following the strategy of Chen and Haiso [27], we construct operational matrices of integration associated with the Fibonacci wavelets, which have less errors in comparison to the operational matrices based on Legendre wavelets [24]. The operational matrices of integration are then employed for converting the models at hand into a system of algebraic equations, which are subsequently solved by the usual Newton method. We studied the DPL model under different boundary conditions using the classical Haar wavelets in [28], however, the spirit of this study is completely different from the aforementioned study in the sense that (i) our method completely relies on the nonorthogonal Fibonacci wavelets generated via Fibonacci polynomials of degree m, (ii) multilayer skin tissues have been considered, (iii) investigation of the effect of volumetric heat loss due to water vaporization has also been studied.
The remainder of the article is structured as follows. In Section 2, we present the mathematical formulation of the model under consideration. Section 3 is completely devoted to the construction of Fibonacci wavelets and their operational matrices of integration. The solution of the DPL model is presented in Section 4. In Section 5, we discuss the obtained numerical results. Finally, a conclusion is drawn in Section 6.

Mathematical Formulation of the Model
As is known, the human skin tissue constitutes three different layers: epidermis, dermis, and subcutaneous layer (see Figure 1). Throughout this study, we shall consider the skin tissue of length L = 0.012. The initial temperature inside the living biological skin tissue is taken as a constant temperature, i.e., T 0 = 37 • C. The Fourier's model of heat transfer for different layers of skin tissue is given by [29,30].
Dermis Layer: Subcutaneous Layer: After applying the phase-lag in both the time and space of the Fourier's model of multilayer tissue, we shall obtain the following DPL models of bioheat transfer for multilayer skin tissue as follows: Epidermis Layer: Dermis Layer: Subcutaneous Layer: The initial condition is and the generalized boundary condition is where in the BC of the first kind- BC of the third kind-A 0 = −k, B 0 = h, f (t) = hT s ; and the corresponding symmetric condition reads The external heat source term, which is considered as a Gaussian distribution heat source term is given by [31,32] where Q ro = ρSP, P, S, r * , and a o represent the transmitted power, antenna constant, location of tumor, and scattering variable, respectively. Moreover, the metabolic heat generation in skin tissue is defined by [32] where Q mo , Q d , Q v , ∆m, and ∆H represent the constant metabolic heat, volumetric heat loss due to water diffusion, volumetric heat loss due to water vaporization, water vaporization rate of skin surface, and enthalpy of water vaporization, respectively [30]. In addition to this, some thermophysical parameters of the multilayer skin tissue are listed in Table 1 of Section 5 with their given values [30]. Table 1. Some parameters of multilayer skin tissue used in this paper [30].

Parameters Epidermis Dermis Subcutaneous
Consider the dimensionless variables given by Subsequently, all of the models (7) to (12) are converted into the following dimensionless form: Epidermis Layer: Dermis Layer: Subcutaneous Layer: The initial condition is given by the generalized boundary condition is given by where we have the following conditions: In the BC of the first kind: In the BC of the second kind: In the BC of the third kind: together the following corresponding symmetric condition: It is worth noticing that the DPL models of bioheat transfer for multilayer skin tissue (17)- (19) reduce to SPL models for multilayer skin tissue by setting τ t = 0 and to classical Pennes bioheat models for multilayer skin tissue by setting τ q = τ t = 0.

Fibonacci Wavelets and Operational Matrices of Integration
The aim of this section is two-fold: First, to provide an exposition to the construction of Fibonacci wavelets via Fibonacci polynomials; second, to derive the corresponding operational matrices of integration by adopting the strategy of Chen and Haiso [27].

Fibonacci Wavelets and Function Approximation
For any η ∈ R + , the Fibonacci polynomials are defined in terms of the recurrence relation asP with initial conditionsP 0 (η) = 0,P 1 (η) = 1 [24]. Equivalently, they can be defined by using the following general formula: and the closed-form formulaP where α and β are the roots of the companion polynomial x 2 − ηx − 1 of the recursion. Furthermore, the power-form representation of the Fibonacci polynomials follows as in [26]: where · denotes the well-known floor function. The Fibonacci polynomialsP m (η) can also be expressed in the matrix form as follows: It is worth noticing that the Fibonacci polynomials satisfies the following properties: The Fibonacci wavelets belong to a special class of the compactly-supported wavelets generated by the Fibonacci polynomials on [0, 1]. Mathematically, the Fibonacci wavelets are defined on the interval [0, 1] as follows (see [24]): whereP m (η) is the Fibonacci polynomial of degree m given by (28), and k and n denote the level of resolution k = 1, 2, · · · and translation parameter n = 1, 2, · · · , 2 k−1 , respectively. The coefficient 1/ √ w m appearing in (33) is a normalization factor and can be computed via (32) as follows: The Fibonacci wavelets (33) can also be represented as given below: where χ I n,k (η) is the characteristic function defined on n−1 2 k−1 , n 2 k−1 . For the choice k = 2 and M = 3, we have the following family of Fibonacci wavelets: Therefore, any square integrable function f ∈ L 2 [0, 1] can be approximated via the Fibonacci wavelets as follows: where h n,m denotes the Fibonacci wavelet coefficients. The matrix analogue of (37) reads as given below: where F is the discrete version of the continuous function f (η) and H is the row vector of the form The matrix Φ(η) appearing in (38) is the Fibonacci wavelet matrix of order 1 × 2 k−1 M and is given by The coefficient vector H can be obtained as follows: where For obtaining the Fibonacci wavelet approximations, we consider the following collocation points:

Operational Matrices of Fibonacci Wavelets
In this subsection, our endeavour is to construct operational matrices of integration for the Fibonacci wavelets (33) by following the strategy of Chen and Haiso [27]: where S is referred to as the Fibonacci wavelet operational matrix of order 2 k−1 M × 2 k−1 M.
Using the power representation (28) and (35) of Fibonacci polynomials and wavelets along with the binomial expansion for 2 k−1 η − n + 1 m−2i , we have Upon integration in (45), we obtain where The function R j (η) can be expressed in terms of Fibonacci wavelets as follows: Therefore, Equation (46) becomes where Γ n,m r,s = Therefore, the operational matrix of integration S for Fibonacci wavelets is formulated as given below: For instance, if we consider k = 2, M = 3, we can integrate (36) at the collocation points given by (43) so that Therefore, (44) takes the following form: where

Solution of the Dual-Phase Model
The present section is devoted to the solving of the dual-phase-lag model of heat transfer in multilayer skin tissue during hyperthermia treatment by employing the Fibonacci wavelet operational matrices of integration constructed in Section 3. However, we shall first employ the finite difference scheme for discretizing the space variable ζ so that the given problem is converted into a system of ordinary differential equations.

Discretizing the Space Variable ζ
In order to curtail the complexity in the DPL models of multilayer skin tissue (17)- (19), we discretize the space coordinate ζ with the aid of the finite difference scheme by taking grid size ∆ζ = h = 1/k, so that the points ζ = h , = 1, 2, · · · , k, represent the grid points of the space interval [0, 1]. Subsequently, the discrete values of the function θ(ζ) at the grid points ζ becomes θ (ζ) = θ(ζ , η), for η 0. Further, the application of the four-point formula at the boundary points yields (see [33]) Similarly, the implementation of the central difference scheme for second-order derivative reduces the DPL models (17) to (19) to the following matrix form: Epidermis Layer: Dermis Layer: Subcutaneous Layer: with initial condition θ(0) = 0, 0, · · · , 0 T , dθ(0) dη = 0, 0, · · · , 0 T , and (59) for BC of the second kind i = 3, for BC of the third kind, where M i and N i , i = 1, 2, 3 are the matrices of order k × k and k × 1, given by and where for BC of the third kind. (67)

Implementation of the Fibonacci Wavelets
The aim of this subsection is to solve the dual-phase-lag models of multilayer tissue by implementing the Fibonacci wavelet operational matrix method described in the previous section.
The highest derivative d 2 θ(η) dη 2 appearing in (56) to (58) can be expanded in terms of the Fibonacci wavelet as where y = 1, 2, and 3 are indicated for multilayer tissues such as the epidermis, dermis, and subcutaneous layers, respectively. Integration of Equation (68) between the limits 0 to η and using the initial condition (59) yields Once again, integrating (69) and using the same condition, we get Substituting the values of d 2 θ(η)/dη 2 , dθ(η)/dη, and θ(η) in Equations (56)-(58), we obtain the following systems of algebraic equations for epidermis, dermis, and subcutaneous layers, respectively: Epidermis Layer: Dermis Layer: Subcutaneous Layer: Solving the above systems of algebraic Equations (71)-(73), we can obtain the unknown Fibonacci wavelet coefficient vector H T y . Finally, substituting the values of H T y in (70), we obtain the approximate solution of the dual-phase-lag model of multilayer tissue.

Numerical Results and Discussion
During hyperthermia treatment, the temperature is raised between 41-46 • C for 15-30 min at the location of the tumor so that the cancer cells become dead without causing any damage to the healthy tissue. To check that the temperature distribution does not surpass 46 • C in the multilayer skin tissue, some thermophysical properties of skin tissue are required for treatment of the tumor. Some of the selected properties of skin tissue are T 0 = T b = T w = 37 • C, c w = 4200 J kg −1 • C, R a = 8.314 J mol −1 , c b = 3344 J kg −1 • C −1 , ∆H = 2405 J kg −1 , M w =18 g mol −1 , Q mo = 1091 W m −3 , L = 0.012 m, ζ * = 0.0082 m, a 0 = −127 m −1 , S = 12.5 kg −1 , and P = 30 W.
The temperature distributions for boundary conditions of the first, second, and third kinds in multilayer skin tissue are shown in Figure 2. It is quite evident from Figure 2 that the tissue temperature at the tumor site for the first kind of boundary condition is 45 • C, whereas for the second and third kind of boundary conditions, it decreases to 42.9 • C and 42 • C, respectively. Thus, we can say that the temperature distribution at the tumor location during hyperthermia treatment is also affected by the boundary conditions. A comparison for different bioheat models such as Penne's, SPL, and DPL models under boundary conditions of the first kind is depicted in Figure 3. The difference in temperature distribution among all the three models is near about 0.5 • C, which is of great importance in the hyperthermia treatment (see Figure 3). Moreover, the temperature distributions for different values of tissue wall temperature θ w in skin tissue are plotted in Figure 4. We can infer from Figure 4 that the temperature of the skin tissue increases with the increase in θ w at the site of a tumor, while in the epidermis layer, the temperature does not show any change with respect to the values of θ w , which indicates that the surface temperature plays an important role in hyperthermia treatment.     Keeping in view that the values of relaxation time τ q (phase-lag due to heat flux) and τ t (phase-lag due to temperature gradient) are both significant in the DPL model during the treatment process, we have plotted the temperature profiles for various values of τ q and τ t in skin tissue in Figures 5 and 6. From Figure 5, a decrease in temperature distribution is observed at the tumor site for both the dermis and subcutaneous layers with an increase in the relaxation time τ t and keeping τ q fixed. Hence, we conclude that the temperature predicted by the SPL bioheat transfer model will always be greater than that of the DPL bioheat transfer model. Likewise, from Figure 6, we demonstrate that the temperature profiles at the area of the tumor increase as the values of relaxation time τ t decrease by taking the fixed value of τ q . Subsequently, from these figures, we infer that the role of thermal relaxation time due to heat flux and temperature gradient is of great importance for an effective treatment.  The effect of blood perfusion on multilayer skin tissue under the boundary conditions of the first kind is depicted in Figure 7. It is evident from Figure 7 that as the values of w b increase, the temperature distribution at the site of the tumor decreases, mainly due to the heat loss caused by blood flow. However, by taking the value of blood perfusion equal to zero, the temperature at the site of the tumor increases by 0.65 • C, which is also effective for precise prediction of the temperature profile during hyperthermia treatment. In Figure 8, the effect of metabolic heat generation P mo on multilayer skin tissue is shown. From this figure, we infer that the temperature profile at the site of the tumor increases with the increase in the metabolic heat source in the tissue. Moreover, the influence of the metabolic heat term on temperature profiles is quite visible in the dermis and subcutaneous layers in comparison to the epidermis layer.   The temperature profiles in skin tissue for the effect of heat loss due to water diffusion Q do are described in Figure 9. It is worth noticing that if we kept Q do = 0 in the dermis layer, the temperature distribution at the site of the tumor remains small, which is extraordinary for the subcutaneous layer at Q do = 0 due to the fact that the thickness of the subcutaneous layer is more prominent than that of the dermis layer.

Conclusions
In this study, the DPL models of bioheat transfer in multilayer skin tissues during hyperthermia treatment under generalized boundary conditions are studied by introducing an efficient wavelet method based on nonorthogonal Fibonacci wavelets. The facilitation has been carried out by constructing the Fibonacci wavelets and the corresponding operational matrices of integration along with a finite element scheme. The effect of different parameters such as tissue wall temperature, time-lag due to heat flux and temperature gradient, blood perfusion, metabolic heat generation, heat loss due to diffusion of water, and boundary conditions of different types on multilayer skin tissues during hyperthermia treatment has been studied and discussed in detail. Nevertheless, it has been observed that the DPL models of bioheat transfer show a lower temperature profile at the site of a tumor in comparison to the SPL and Fourier models. Thus, we conclude that these models give us a more precise prediction of temperature distribution at the tumor site. Finally, it is hoped that the proposed wavelet method can be employed for solving a wide range of physical and biological problems.

Conflicts of Interest:
The authors declare that they have no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: Arterial blood temperature ( • C)