Application Of Fractional Calculus In Modeling And Solving .

Application of fractional calculus in modelingand solving the bioheat equationR. Magin1, Y. Sagher2 & S. Boregowda1,31Department of Bioengineering, University of Illinois at Chicago, USADepartment of Mathematical Sciences, Florida Atlantic University, USA3University of Chicago Graduate School of Business, USA2AbstractFractional calculus provides novel mathematical tools for modeling physical andbiological processes. The bioheat equation is often used as a first order model ofheat transfer in biological systems. In this paper we describe formulation of bioheat transfer in one dimension in terms of fractional order differentiation withrespect to time. The solution to the resulting fractional order partial differentialequation reflects the interaction of the system with the dynamics of its responseto surface or volume heating. An example taken from a study involving pulsating(on-off) cooling of a peripheral tissue region during laser surgery is used toillustrate the utility of the method. In the future we hope to interpret the physicalbasis of fractional derivatives using Constructal Theory, according to which, thegeometry biological structures evolve as a result of the optimization process.Keywords: bioheat transfer, diffusion, fractional calculus, modelling, lasersurgery, fractals, temperature.1IntroductionThe present paper considers the application of fractional calculus to the analysisof problems in bioheat transfer. The methods of fractional calculus, reviewedrecently by Magin [1], are developed as the basis for formulation and solution ofthe bioheat transfer problem in peripheral tissue regions. Investigators havestudied bioheat transfer using mathematical models for more than 50 years [5-7].In these models tissue cooling (or warming) is approximated by coupling tissueperfusion to the bulk tissue temperature through Newton’s law of cooling (orheating). In addition to full body models, there are numerous models in literatureDesign and Nature II, M. W. Collins & C. A. Brebbia (Editors) 2004 WIT Press, www.witpress.com, ISBN 1-85312-721-3

208 Design and Nature IIthat describe heat transfer mechanisms in a single organ or a portion of the body.In this regard, an analytical model developed by Keller and Seiler examines bioheat transport phenomena with heat generation (metabolism) occurring in theperipheral tissue regions. The Keller and Seiler [8] model was solvednumerically using parallel computers to simulate all possible modes of bioheattransfer by Boregowda et al. [9].Recently a number of investigators [10-12] have applied the bioheat transfermodel to periodic diffusion problems in localized tissue regions such as thatwhich occurs in the skin when laser heating and/or cryogen cooling is applied.Fractional calculus is ideally suited to address this kind of periodic heating orcooling, but to our knowledge has not been used in modeling bioheat transfereither at the tissue, organ or whole body level. The present study demonstratesthat fractional calculus can provide a unified approach to examine periodic heattransfer in peripheral tissue regions. For example, in an experimental studyconducted by Pikkula et al. [13], cryogen spray cooling is utilized to cool theskin surface during the laser skin surgery. A generalized fractional calculusapproach developed by Kulish and Lage [14 -16] is adopted to model thelocalized periodic bioheat transfer problems similar to the one posed by Pikkulaet al. [13].The one-dimensional heat flow problem can be completely solved for welldefined surface temperature or thermal flux boundary conditions by applyingLaplace transforms [17,18]. The solution can also be expressed as a fractionaldifferential equation for the semi-infinite peripheral tissue region [14,15].Further, the fractional differential equation can be solved to compute the heatflux at the boundary for different periodic or on-off boundary conditions thatclosely represent the heating and cooling of skin surface during laser surgery.The approach offered by fractional calculus models a large class ofbiomedical problems that involve localized pulse heating and/or cooling. Oneadvantage of this approach is that there is no need to solve first for thetemperature in the entire domain.2General formulationThe approach used in this study is an approximation to the physical modeldeveloped in the study by Deng and Liu [10]. The region of interest is theboundary and its vicinity, and the total thickness is assumed to be large, so thatrectangular coordinates in one dimension can be used for the analysis. Note thatthe outermost portion, the skin, is considered to be thin so that its thickness is notexplicitly incorporated into the model. The localized tissue region that isrepresented by this approximate physical model is shown in fig. 1.The generalized one-dimensional bioheat transfer equation for thetemperature T ( x, t ) in the tissue developed by Pennes [2] can be written as: T ( x, t ) 2T ( x, t )ρc K ωbρb cb (Ta T ( x, t )) Qm Qr ( x, t ) t x 2Design and Nature II, M. W. Collins & C. A. Brebbia (Editors) 2004 WIT Press, www.witpress.com, ISBN 1-85312-721-3(1)

Design and Nature II209where ρ, c, and K are the density, specific heat and thermal conductivity of thetissue and ρb , cb the density and specific heat of the blood, ωb is the bloodperfusion, Ta is the arterial blood temperature (assumed to be constant), Qm isthe metabolic heat generation and Qr ( x, t ) the heat generation due to spatialheating in the medium.Figure 1:Assumed physical model of the localized tissue region. T ( x, t )represents the temperature in the tissue, while Φ (t ) describes thesurface thermal flux at x 0 .We assume that the problem has the following boundary conditions:T ( x,0 ) Ti ( x,0) T (0, t )Φ (t ) K xlim T ( x, t ) TCx initial temperature distributionsurface fluxconstant core temperatureIf we initially assume Qr to be zero we can solve this problem following Liu et al. [12] in terms of T ( x, t ) T ( x, t ) Ti ( x,0) where we have subtracted theDesign and Nature II, M. W. Collins & C. A. Brebbia (Editors) 2004 WIT Press, www.witpress.com, ISBN 1-85312-721-3

210 Design and Nature IIinitial temperature distribution Ti (x,0) which is just the solution of steady stateproblem. Applying the Laplace transformation to eqn. (1) for the given boundaryconditions we obtain for h ωb ρb cb ρc and k K ρc . t ( x, s )k ( s h) t ( x, s ) 02 x t ( x, s ) φ( s ) K x x 0 T ( x,0 ) 0 ,lim t ( x, s ) 0 ,x This second order ordinary differential equation has the following solution forthe specified boundary conditionsk φ( s )e x ( s h ) k t ( x, s ) .K s h(2)If we consider only the relationship between the flux and the temperature atthe x 0 boundary, then the result can be written in terms of a Laplaceconvolution integral ask T (0, t ) Kt 0e hτke htΦ (t τ)dτ Φ (t ) ,Kπτπt(3) 1 e ht.where we have used the Laplace transform pair L πt s h 1Thus, if the surface flux is modelled by Φ (t ) Φ 0u (t ) , where u (t ) is the unitstep function then the surface temperature will increase ask T (0, t ) Kt 0Φ 0 e hτk Φ0 2dτ K h ππτht 2e u du 0k Φ0erfK hwhere u hτ and the error function is defined by erf ( x ) 2( ht )2 x u 2 e du.π 0However, the convolution integral, eqn. (3) can also be written ask T (0, t ) Kt 0k e htΦ ( τ) e h ( t τ )dτ K ππ(t τ)Design and Nature II, M. W. Collins & C. A. Brebbia (Editors) 2004 WIT Press, www.witpress.com, ISBN 1-85312-721-3t 0Φ (τ)e hτdτ.t τ(4)

Design and Nature II211which can be written in terms of the Riemann-Liouville fractional integral[19-22] defined by 1 t F (τ)dτ , where Γ(α) u α 1e u du .0 D F (t ) 1 αΓ(α) 0 (t τ)0 αtThus, eqn. (4), can be simply expressed in terms of fractional integration by[]k ht 1 2 hτ T (0, t ) e 0 Dt e Φ (τ) .KIf we assume a step input in flux at x 0, Φ (t ) Φ 0 u (t ) , we can write[]k ht 1 2 T (0, t ) e 0 Dt Φ 0 e ht ,Kwhich since fractional integral is a linear operator and the fractional derivative[ ] 1 2 hte 0 Dte hterfh( ht )[21], gives the same result for the surfacetemperature as that obtained above by inversion of the Laplace transform.In the case of a specified surface temperature at the surface x 0 , a parallelanalysis gives the surface flux in terms of the fractional semiderivative of thesurface temperature, which can be writtenΦ (t ) []K ht 1 2 hte 0 Dt e T (0, t ) ,k(5)where the fractional derivative of order 1 2 is defined [21] as120 Dt F (t ) 1 d t F ( τ)dτ .Γ(1 2) dt 0 (t τ)1 2This result can also be obtained using Babenko’s method [22,23].Thus, for the case where T (0, t ) T0u (t ), a step in surface temperatureof T0 at x 0 , and using the semiderivative of eDesign and Nature II, M. W. Collins & C. A. Brebbia (Editors) 2004 WIT Press, www.witpress.com, ISBN 1-85312-721-3ht[21] we obtain

212 Design and Nature IIΦ (t ) KT0 e ht h erf k πt( ht ) .(6) A graph of this result is shown in fig. (2).Figure 2:Graph of the flux Φ (0, t ) necessary to establish a step input intemperature T0u (t ), assuming A K k 1. Two cases areplotted: one for the bioheat equation with h 1, and a second fornormal diffusion without blood flow cooling, e.g., h 0 .Since the relationship between flux and temperature is assumed to followfrom the Fourier law for heat flux it is valid at any point in the domain, not onlyat the x 0 surface. Therefore, for the one-dimensional problem of heatingwith linear surface cooling this allows us to write our fractional integral andderivative results asΦ ( x, t ) and T ( x, t ) []K ht 1 2 ht e 0 Dt e T ( x, t ) ,k[]k ht 1 2 hte 0 Dt e Φ ( x, t ) .KDesign and Nature II, M. W. Collins & C. A. Brebbia (Editors) 2004 WIT Press, www.witpress.com, ISBN 1-85312-721-3(7)(8)

Design and Nature II213Thus, given the flux or temperature profiles at a specific location we can usethis information to determine the corresponding temperature or flux. Thisapproach could be useful in experimental situations where the half-orderfractional integrals or derivatives of known functions could be used to determinethe required input conditions needed for desired temperature or flux outputs [2426]. A few examples are listed in table 1, which is adapted from Oldham andSpanier [21].Note that, if the initial temperature distribution Ti (x,0) is assumed to be uniform and constant, e.g., Ti ( x,0) T0 , then T ( x, t ) T ( x, t ) T0 and theflux expression, eqn. (7), becomes[]K ht 1 2 htKTe 0 Dt e T ( x, t ) 0Φ ( x, t ) kk e ht h erf πt( ht ) . This equation simplifies for h 0 toΦ ( x, t ) K 12T0 0 Dt T ( x, t ) k πt which was previously derived by Kulish and Lage [14]. Kulish and Lage haverecently applied fractional-diffusion theory to thermoreflectance measurementsof the thermal properties of thin films under pulsed laser heating. The currentbioheat model under conditions of volumetric as well as surface heating extendsKulish’s results, eqn. (10) in [16], to yieldΦ ( x, t ) [[K ht 1 2 ht e 0 Dt e T ( x, t )k] ] e ht 0[] P Dt1 2 e ht P ( x, t ) K xwhere P ( x, t ) P ( x, t ) P ( x,0) and P ( x, t ) represents the particularsolution to the Laplace domain inhomogeneous ordinary differential equation.In this short paper we have described a fractional calculus approach to theformulation of the bioheat equation. This method provides a simple expressionfor either the temperature or flux under experimental conditions often specifiedby laser heating and cryogen-cooling procedures. Additional studies are neededto develop a connection between the fractional order of the operators and thematerial structure and properties of the tissue or substate under study. Recentwork by West et al. [27] and others [28-30] is directed toward establishing astronger role for fractional calculus in describing dynamic phenomena incomplex materials.Design and Nature II, M. W. Collins & C. A. Brebbia (Editors) 2004 WIT Press, www.witpress.com, ISBN 1-85312-721-3