Bioheat Transfer with Thermal Memory and Moving Thermal Shocks

: This article investigates the effects of thermal memory and the moving line thermal shock on heat transfer in biological tissues by employing a generalized form of the Pennes equation. The mathematical model is built upon a novel time-fractional generalized Fourier’s law, wherein the thermal ﬂux is inﬂuenced not only by the temperature gradient but also by its historical behavior. Fractionalization of the heat ﬂow via a fractional integral operator leads to modeling of the ﬁnite speed of the heat wave. Moreover, the thermal source generates a linear thermal shock at every instant in a speciﬁed position of the tissue. The analytical solution in the Laplace domain for the temperature of the generalized model, respectively the analytical solution in the real domain for the ordinary model, are determined using the Laplace transform. The inﬂuence of the thermal memory parameter on the heat transfer is analyzed through numerical simulations and graphic representations.


Introduction
In the field of biological sciences, bioheat transfer is the intriguing and interdisciplinary study of heat transport inside living systems. This intriguing area of research focuses on a variety of processes, including heat transfer inside cells and tissues and the thermoregulatory mechanisms that enable animals to control their body temperatures. Bioheat transfer examines the principles of thermodynamics, fluid mechanics, and physiological processes in order to comprehend the intricate relationship between heat and life.
The way that living tissues work is fundamentally impacted by heat as an energy type. Numerous physiological functions are impacted, including the maintenance of enzymatic reactions and the control of metabolic activity. One must have a complete understanding of the generation, transport, and dissipation of heat within biological systems in order to understand the behavior and operation of such systems. Researchers in this rapidly evolving field employ theoretical models, experimental techniques, and computer simulations to better understand the intricate dynamics of heat transport in biological systems. By combining knowledge from several domains, these interdisciplinary techniques increase the bonds between biologists, physicists, engineers, mathematicians, and medics [1][2][3][4].
Thermal shock can occur when a biological system's temperature changes significantly and fast. In the context of bioheat transport, the term "thermal shock" refers to the detrimental effects that sudden temperature fluctuations have on biological tissues and organisms. Understanding the physics and impacts of thermal shock is crucial in many fields, including biomedical engineering, cryobiology, thermal medicine, and thermal therapy [5,6].
The fractional integral is a mathematical concept that broadens the definition of the term "integral" to include orders other than integers. It offers a potent tool for exploring and modifying functions across a range of mathematical and scientific fields. Due to its versatility and broad variety of applications, the study of fractional integrals has attracted a lot of attention in recent years [7,8].
The area of bioheat transfer investigates how heat moves and is distributed within biological systems. The Pennes bioheat equation and other conventional partial differential equations have been used to model heat transfer in biological tissues. Instantaneous heat conduction and local thermal equilibrium are presupposed in these models [9][10][11][12][13].
The concept of differentiation and integration is broadened to encompass non-integer orders through the use of fractional calculus. Because fractional calculus takes into account memory-like effects, non-locality, and long-range interactions, it may be used to represent complex physical processes like bioheat transmission. The fractional integral operator, a crucial tool in fractional calculus, enables the description of heat transfer processes with memory and long-term repercussions. The typical non-local heat conduction behavior observed in biological tissues results from a combination of factors including blood perfusion, capillary networks, and metabolic heat sources. By including non-local heat conduction components into bioheat transfer equations using fractional integral models, more precise estimates of the temperature distribution inside tissues may be achieved. Fractional integral operators can be used to explain events in biological tissues that mimic memory. Understanding these effects, which are the result of prior temperature exposure, is necessary to comprehending thermal responses and adaptive processes. Models can better accurately depict the long-term effects of previous heat exposure on tissue temperature distribution by employing fractional integral operators. Thermal treatments and diagnostic procedures may benefit from the use of fractional integral operators in bioheat transfer models. Treatment planning and optimization may be strengthened by taking non-local and memory-like effects into account. This will result in thermal therapies such as cancer hyperthermia treatments being more successful. Models that take into consideration memory effects can also help in the interpretation of temperature readings and the creation of diagnostic tools [14][15][16][17][18].
The mathematical model of heat transfer proposed by Pennes [19] has been widely studied in different circumstances. Also, different fractional forms of the heat transfer equation in biological tissues have been investigated. We must mention that most of the time, the studied fractional equations were obtained from the classical form in which, the temperature derivative in relation to time was artificially replaced with a fractional derivative. Even if from a mathematical point of view things are correct, there is no physical argument in favor of this construction.
In the present work, the fractional model is built starting from the constitutive equation of the thermal flow. The dimensionless form of Fourier's law is generalized by proposing a constitutive relationship in which the thermal flow is determined not only by the temperature gradient but also by its history. The choice of the relaxation (damping) function of the heat flux is very important in the formulation of the bio-heat equation with memory.
The solution corresponding to a problem with a fractional operator contains as a parameter the fractional order of the operator. The presence of this parameter in the solution could lead to the optimal choice of the parameter value in order to achieve a solution with a certain property.
On the other hand, fractional models could be useful tools in choosing the mathematical model that leads to values close to the values obtained through experiments. In this way, an analytical model could be associated with a set of experimental values.
In this paper we develop a mathematical model by fractionalization of the temperature gradient by the fractional integral Riemann-Liouville operator. The proposed mathematical model is new in the specialized literature. Another interesting aspect is given by the property of the heat source we considered. The proposed form allows the application, at Using the Laplace transform (LT), the analytical solution in the Laplace domain for the temperature distribution in the biological tissue is obtained. Numerical values, in the real domain for temperature are determined using the Stehfest numerical algorithm for the inversion of the LT. The analytical solution for the ordinary case (α = 0) is also determined. The influence of the memory parameter on the temperature distribution is investigated.

Statement of the Problem
Heat transfer in a one-dimensional biological tissue with length L, subjected to the moving thermal shocks is considered (Figure 1).
In this paper we develop a mathematical model by fractionalization of the temperature gradient by the fractional integral Riemann-Liouville operator. The proposed mathematical model is new in the specialized literature. Another interesting aspect is given by the property of the heat source we considered. The proposed form allows the application, at certain moments of time, of a thermal shock in a specified area of the tissue. This property could have applications in various medical treatments.
Using the Laplace transform (LT), the analytical solution in the Laplace domain for the temperature distribution in the biological tissue is obtained. Numerical values, in the real domain for temperature are determined using the Stehfest numerical algorithm for the inversion of the LT. The analytical solution for the ordinary case ( 0   ) is also determined. The influence of the memory parameter on the temperature distribution is investigated.

Statement of the Problem
Heat transfer in a one-dimensional biological tissue with length L , subjected to the moving thermal shocks is considered (Figure 1). The governing equations of heat transfer are [18,20]: Fourier's Law where ( )   is the Dirac function. The special form of the heat source 0 ( ) shows that at each time instant t , a constant heat source 0 Q is applied in the position x vt    . In Equation (1) the parameter 0 v  denotes the velocity of the moving heat source.
Along with Equations (1) and (2) the following conditions are considered: The governing equations of heat transfer are [18,20]: Energy balance equation.
Fourier's Law where δ(·) is the Dirac function. The special form of the heat source Q 0 δ(v t − x) shows that at each time instant t, a constant heat source Q 0 is applied in the position x = v t. In Equation (1) the parameter v > 0 denotes the velocity of the moving heat source. Along with Equations (1) and (2) the following conditions are considered: Introducing the dimensionless variables, functions and parameters, Equations (1) and (2) become The non-dimensionless forms of the conditions are ∂θ(x, t) ∂t

Fractional Mathematical Model
We considered the fractional constitutive equation of the thermal flux where is the fractional integral Riemann-Liouville operator [21][22][23]. In Equation (12), It is used as a normalization coefficient. The Laplace transform of the fractional integral operator (11) is given by whereθ(x, t) is the LT of the function θ(x, t) [24]. For α = 0, the Equation (10) reduces to Equation (7). Equation (10) can be written as highlights that time-variation of the thermal flux is significantly influenced by the history of the temperature gradient whose values are scaled by the kernel h α (t) of the power law.

Semi-Analytical Solution for the Generalized Bio-Heat Equation
To determine the solution of Equations (6) and (10), the Laplace transform method is used. Applying the LT to Equations (6) and (10) and using the initial value (8) we obtain Using (16) into Equation (15), we obtain the ordinary differential equation for the transformed temperatureθ(x, s): Satisfy the conditions A particular solution of the non-homogenous Equation (17) iŝ where The general solution of the homogeneous equation associated to Equation (17) iŝ and the general solution of differential Equation (17) iŝ Using the boundary condition (18), we get Finally, we obtain The analytically expressionθ(x, s) of the LT is complicated. Therefore, determining the inverse of this transform using complex analysis is difficult. Even if this inverse transform can be determined, its expression is complicated enough to be useful in numerical calculations. For these reasons, we will determine the values of the inverse Laplace transform using Sthefest's numerical algorithm [25].
According with Sthefest's algorithm, the numerical values of the function θ(x, t) are given by where γ ∈ [q] and [q] denotes the integer value function. Formula (25) that generates the numerical values of θ(x, t) temperature will be used in Section 3 of this work. In this particular case, Functions (20)-(24) have simplified forms. We will analyze these functions in turn.
The inverse LT of Equation (26) is [26] a(t) = − 2Q 0 On the other hand, In the considered particular case, functionθ(x, t) given by Equation (24) has the following expression: We consider: we get the inverse Lapace transform of F m (s) given by Now, using (35) and the formula [26] The last function that must be inverted is The inverse LT of H(s) is Finally, the temperature θ(x, t) is given by In this section, the analytical solution was determined in the Laplace domain for the temperature of the biological tissue, expression given by Equation (24). In particular, in the ordinary case corresponding to α = 0, the expression of the temperature in the real domain given by Equation (40) was determined. The values in the real domain of the temperature (24) will be obtained using Equation (25). The discussion of the numerical results for certain values of the problem parameters in presented in Section 3 of the paper.
Regarding the particular case α = 0, we make a few remarks, namely: There are many published articles in which the Pennes equation has been investigated considering different forms of heat sources or boundary conditions. We only mention two articles [27,28] in which the authors studied the problem of microwave heating of skin tissue together with the effects of surface cooling and blood flow, respectively, three heat transfer models for radiofrequency ablation. The form of the thermal source considered by Wang et al. [27], allowed the authors to use an elegant method for the integration of the studied partial differential equation. They transformed the initial boundary value problem into a standard heat equation problem. The different form of the thermal source considered by us with respect to Wang et al. [27] does not allow the comparison of the results of this article with those obtained in Wang et al [27].
The expression of the heat source with the Dirac distribution in the problem studied by us generates difficulties in the use of classical methods from the theory of partial differential equations. This is one reason why we used the LT to determine the solution. On the other hand, the presence of the fractional integration operator in the equations of the generalized model requires the use of the LT method to determine the analytical solution of the considered problem.
To verify the correctness of the analytical solution (40) we used the form of the general solution (24) and determined its numerical values for alpha tending to zero. The results are presented in Figure 2. It is observed that when α → 0 the graphs of the generalized solution tend to overlap the graph of the analytical solution given by Equation (40). This fact shows the correctness of the solutions determined by us. the thermal source considered by us with respect to Wang et al. [27] does not allow the comparison of the results of this article with those obtained in Wang et al [27]. The expression of the heat source with the Dirac distribution in the problem studied by us generates difficulties in the use of classical methods from the theory of partial differential equations. This is one reason why we used the LT to determine the solution. On the other hand, the presence of the fractional integration operator in the equations of the generalized model requires the use of the LT method to determine the analytical solution of the considered problem.
To verify the correctness of the analytical solution (40) we used the form of the general solution (24) and determined its numerical values for alpha tending to zero. The results are presented in Figure 2. It is observed that when 0   the graphs of the generalized solution tend to overlap the graph of the analytical solution given by Equation (40). This fact shows the correctness of the solutions determined by us.

Numerical Results and Discussions
In this paper, a generalized form of the Pennes equation that models heat transfer in biological tissues is investigated. To analyze the complex thermal diffusion that occurs in such tissues, we have considered the thermal flux defined by a fractional differential equation with the Riemann-Liouville integral operator. In such thermal processes, the thermal flux is influenced not only by the temperature gradient but also by its history. In addition, the heat transfer is studied under the hypothesis of the existence of a mobile heat source that, at every moment, generates a thermal shock in a well-established position in the field of biological tissue. The proposed equation is brought to the non-dimensional form by conveniently choosing dimensionless variables. The analytical solution in the Laplace domain for the dimensional temperature was determined using the LT. Because the expression of the LT of the non-dimensional temperature is complicated, the determination of the inverse LT is difficult to find with complex analysis methods. Therefore, the determination of the temperature values in the real domain was achieved with the help of the Stehfest numerical algorithm.
To analyze the influence of the thermal memory parameter  on the heat transfer process, numerical simulations have been carried out. The obtained results are presented in graphical illustrations. In all the cases analyzed, the particular case of Fourier's law corresponding to the value 0   of the fractional parameter is also presented. In this way, the difference between the classic heat transfer and the one based on the fractional

Numerical Results and Discussions
In this paper, a generalized form of the Pennes equation that models heat transfer in biological tissues is investigated. To analyze the complex thermal diffusion that occurs in such tissues, we have considered the thermal flux defined by a fractional differential equation with the Riemann-Liouville integral operator. In such thermal processes, the thermal flux is influenced not only by the temperature gradient but also by its history. In addition, the heat transfer is studied under the hypothesis of the existence of a mobile heat source that, at every moment, generates a thermal shock in a well-established position in the field of biological tissue. The proposed equation is brought to the non-dimensional form by conveniently choosing dimensionless variables. The analytical solution in the Laplace domain for the dimensional temperature was determined using the LT. Because the expression of the LT of the non-dimensional temperature is complicated, the determination of the inverse LT is difficult to find with complex analysis methods. Therefore, the determination of the temperature values in the real domain was achieved with the help of the Stehfest numerical algorithm.
To analyze the influence of the thermal memory parameter α on the heat transfer process, numerical simulations have been carried out. The obtained results are presented in graphical illustrations. In all the cases analyzed, the particular case of Fourier's law corresponding to the value α = 0 of the fractional parameter is also presented. In this way, the difference between the classic heat transfer and the one based on the fractional model is highlighted. Figure 3 has been sketched to highlight the tissue temperature variation for three values of t and two values of the parameter b 0 . Considering the special shape of the heat source, the dimensionless parameter b 0 has a significant influence on the temperature distribution. Let us note that, at a moment of time t, the thermal shock is applied in the position x = t/b 0 in the biological tissue. This is clearly highlighted in Figure 3 where the maximum temperature is achieved in such positions. It should also be noted that until the thermal shock is applied, the temperature corresponding to the fractional model is higher than that corresponding to the classic model. There is a position in the tissue located after the one in which the temperature has the maximum value, a position after which the temperature corresponding to the classical model becomes higher than that of the generalized model. The convenient choice of the fractional parameter α and the b 0 parameter could lead to optimal heating/cooling of certain tissue areas. b has a significant influence on the temperature distribution. Let us note that, at a moment of time t, the thermal shock is applied in the position in the biological tissue. This is clearly highlighted in Figure 3 where the maximum temperature is achieved in such positions. It should also be noted that until the thermal shock is applied, the temperature corresponding to the fractional model is higher than that corresponding to the classic model. There is a position in the tissue located after the one in which the temperature has the maximum value, a position after which the temperature corresponding to the classical model becomes higher than that of the generalized model. The convenient choice of the fractional parameter  and the 0 b parameter could lead to optimal heating/cooling of certain tissue areas.     The non-dimensional blood perfusion parameter a 0 gives the behavior of the convective term (−a 0 θ(x, t)) in Equation (6), therefore this parameter significantly influences the temperature distribution. Figure 5 shows profiles of the tissue temperature when the values of the a 0 parameter increase. The temperature values are decreasing with the parameter a 0 because the convective term plays the role of a cooling condition.
The influence of heat source parameter b 0 is shown in Figure 6. The temperature profiles are drawn in the middle of the biological tissue. Obviously, for different values of the source parameter b 0 , the position of the thermal shock is changed; therefore, the temperature values are oscillating.  In the final part of this section, we investigate the variation of the thermal flux of the proposed model. It is known that in thermal processes the rate of heat transfer is characterized by thermal flux. Using Equations (16) and (24) Values of the thermal flux in the real domain are determined using Stehfest's algorithm. In the final part of this section, we investigate the variation of the thermal flux of the proposed model. It is known that in thermal processes the rate of heat transfer is characterized by thermal flux. Using Equations (16) and (24) we obtain the LT of the thermal flux given bŷ Values of the thermal flux in the real domain are determined using Stehfest's algorithm. Figure 7 shows the heat flux profiles at three moments of time and two values of the heat source parameter b 0 . For each time value, the curves corresponding to the classic model (alpha = 0) and the generalized model (alpha = 0.1) were drawn. The numerical data chose were the same as in Figure 3 to see the correlation between the evolution in space of the thermal flux and of the temperature. Obviously, there is a very good correlation between the graphs in Figures 3 and 7. In the intervals in which the thermal flux has negative values, the temperature increases and in the intervals in which the thermal flux has positive values, the temperature decreases. Figure 7 shows the heat flux profiles at three moments of time and two values of the heat source parameter b0. For each time value, the curves corresponding to the classic model (alpha = 0) and the generalized model (alpha = 0.1) were drawn. The numerical data chose were the same as in Figure 3 to see the correlation between the evolution in space of the thermal flux and of the temperature. Obviously, there is a very good correlation between the graphs in Figures 3 and 7. In the intervals in which the thermal flux has negative values, the temperature increases and in the intervals in which the thermal flux has positive values, the temperature decreases.

Conclusions
A semi-analytical solution for the time-fractional generalized Pennes bio-heat equation in the presence of a localized moving heat source and a constant metabolic heat generation source has been obtained using the LT and a numerical algorithm for LT inversion.
The mathematical model considered is based on the time-fractional generalized Fourier's law, which highlights certain aspects of anomalous thermal diffusion. The constitutive equation of the thermal flux is described by the time-fractional Riemann-Liouville integral operator. As a result, the thermal flux is influenced not only by the temperature gradient but also by its historical values.
The thermal source generates a linear thermal shock at every instant at a specified position within the tissue.

Conclusions
A semi-analytical solution for the time-fractional generalized Pennes bio-heat equation in the presence of a localized moving heat source and a constant metabolic heat generation source has been obtained using the LT and a numerical algorithm for LT inversion.
The mathematical model considered is based on the time-fractional generalized Fourier's law, which highlights certain aspects of anomalous thermal diffusion. The constitutive equation of the thermal flux is described by the time-fractional Riemann-Liouville integral operator. As a result, the thermal flux is influenced not only by the temperature gradient but also by its historical values.
The thermal source generates a linear thermal shock at every instant at a specified position within the tissue.
Numerical simulations demonstrate that the fractional model is capable of generating higher tissue temperatures than the classical model based on Fourier's law. By choosing the appropriate parameters for the mobile source that generates the thermal shock, it becomes possible to induce heating or cooling in specified areas of the biological tissue.