Anomalous properties of heat diffusion in living tissue caused by branching artery network. Qualitative description I. A. Lubashevsky,1,2 V. V. Gafiychuk,2 and B. Y. Datsko2 1Theory Department, General Physics Institute, Russian Academy of Sciences 2Institute for Applied Problems in Mechanics and Mathematics, National Academy of Sciences of Ukraine (Dated: February 1, 2008) 2 Weanalyzetheeffectofbloodflowthroughlargearteriesofperipheralcirculationonheattransfer 0 inlivingtissue. Bloodflowinsucharteriesgivesrisetofastheatpropagationoverlargescales,which 0 is described in terms of heat superdiffusion. The corresponding bioheat heat equation is derived. 2 In particular, we show that under local strong heating of a small tissue domain the temperature distribution inside the surroundingtissue is affected substantially by heat superdiffusion. n a J I. INTRODUCTION. LIVING TISSUE AS A 5 HETEROGENEOUS MEDIUM ] t f Blood flowing throughvessels forms paths of fast heat o transport in living tissue and under typical conditions s it is blood flow that governs heat propagation on scales . t a aboutorgreaterthanonecentimeter(foranintroduction m tothisproblemsee,e.g.,[1,2]). Bloodvesselsmakeupa complex network being practically a fractal. The larger - d is a vessel, the faster is the blood motion in it and, so, FIG. 1: Schematic illustration of the effect of blood flow n the strongeristhe effectofbloodflow inthe givenvessel through the vein tree on heat diffusion imitated by random o on heat transfer. Blood flow in capillaries practically walks. Thefigureshowstrappingofarandomwalkerbecause c does not affect heat propagation whereas blood inside ofgettingtheinternalpointsof alarge veinafter passingthe [ large vessels moves so fast that its heat interaction with vein node. 1 the surrounding cellular tissue is ignorable [1]. Thus, v there should be vessels of a certain length ℓ that are 7 v flow on heat transfer is reduced to the renormalization the smallest ones among the vessels wherein blood flow 5 of the temperature diffusivity, D Deff, [5] and the affects heat transfer remarkably. The value of ℓ can be → 0 v appearance of the effective heat sink fj [1, 3, 6] in the 1 estimated as [3] (see also [1, 2]): bioheat equation: 0 2 D ∂T = (D T) fj(T T )+q . (2) 0 ℓv , (1) ∂t ∇ eff∇ − − a T / ∼sjfLn at Here T(r,t)is the tissue temperature fieldaveragedover m where D = κ/(cρ) is the temperature diffusivity of the scales about ℓ , the parameter T is the blood temper- v a cellular tissue determined by its thermal conductivity κ, ature insider the systemic circulation arteries, and the - d specific heat c, and density ρ, the value j is the blood summand q (r,t) called below the temperature gener- T n perfusion rate (the volume of blood going through tis- ation rate is specified by the heat generation rate q as o sue region of unit volume per unit time), and the fac- q = q/(cρ). The renormalization of the temperature T :c tor Ln ∼ ln(l/a) is logarithm of the mean ratio of the diffusivity is mainly determined by the blood vessels of v individual length to radius of blood vessels forming pe- lengthsaboutℓ andduetothefractalstructureofvascu- v Xi ripheral circulation. For the vascular networks made up lar networks the renormalization coefficient F = Deff/D of the paired artery and vein trees where all the vessels is practically a constant of unity order, F &1 [3]. r a are grouped into the pairs of the closely-spaced arter- Let us imitate the temperature evolution in terms of ies andveins withopposite bloodcurrentsthe coefficient randomwalkswhoseconcentrationis(T T ). Thenthe a f L−1/2 accounts for the counter-current effect [15]. part of the vein tree made up of vessel−s whose lengths n ∼ Forthevascularnetworkswherethearteryandveintrees exceedorareaboutthescaleℓ formsthesystemoftraps. v are arranged independently of each other the factor f Infact, bloodstreamsgoingthroughthe veintree merge shouldbesetequaltounity,f =1. Inparticular,forthe into greater streams at the nodes (Fig. 1). Therefore typical values of the ratio l/a 20–40 [4], the thermal an effective random walker after reaching the boundary conductivity κ 7 10−3 W/∼cmK, the heat capacity of one of these veins inevitably will be moved by blood c 3.5J/g K,a∼ndth·edensityρ ·1 g/cm3 ofthetissue flow into the internal points of large veins. Then, due to as∼well as s·etting the blood perfu∼sion rate j 0.3min−1 relatively fast blood motion inside these veins it will be ∼ from(1)wegettheestimatesℓ 4mmandL 3 4. carried away from the tissue region under consideration, v n ∼ ≈ − In the mean-field approximation the effect of blood which may be described in terms of the walker trapping 2 or, what is the same, the heat sink [3]. Since the mean distancebetweentheseveinsisdeterminedmainlybythe shortestones,i.e. bytheveinsoflengthℓ themeantime v during which a walker wanders inside the cellular tissue before being trapped is [3] ℓ2 τ vL . (3) n ∼ D In obtaining the given expression we have assumed the vascular networkto be embedded uniformly in the cellu- lar tissue, so the tissue volume ℓ3 falls per one vein (and v artery, respectively) of length ℓ . Whence it follows, in FIG. 2: Schematic illustration of the anomalously fast heat v propagationcausedbybloodflowthroughlargearterytreein particular, that the rate at which the walkers are being termsofrandomwalks. Thefigureshowstheeffectivewalker trapped by these vein, i.e. the rate of their disappear- motion with blood from a large artery to a small one where ance is estimated as 1/τ, leading together with expres- thewalker leaves it wandering in thecellular tissue. sion (1) to the heat sink of intensity fj in the bioheat equation (2). The characteristic spatial scale of walker diffusion in the cellular tissue before being trapped is ℓ √Dτ, i.e. example, inside a cellular tissue neighborhoodof a point T ∼ r can be found during the time τ inside a cellular tissue neighborhoodofapointr′ atadistancemuchlargerthan D ℓ ℓ L 1cm. (4) ℓ ,i.e. r r′ ℓ . Thegiveneffectmayberegardedas T v n v v ∼ ∼sjf ∼ anomal|ou−sly f|a≫st heat diffusion in living tissue, i.e. heat p superdiffusion. The scale ℓ gives us also the mean penetration depth T Dealingwithheattransferinlivingtissuewemaycon- of heat penetration into the cellular tissue from a point fineourconsiderationtotheperipheralvascularnetworks source or,what is the same,the widening of the temper- typically embedded uniformly into the cellular tissue, at aturedistributioncausedbyheatdiffusioninthecellular leastatthe firstapproximation[8]. Thelatterstatement tissue. It is the result obtained within the mean-field means, in particular, the fact that for a fixed peripheral approximation. vascular network the vessel collection comprising all the Beyondthescopeofthemean-fieldtheorywemeetsev- arteries of length l meets the condition of the volume eralphenomena. Oneofthemisthetemperaturenonuni- l3 approximately falling per each one of these arteries. formities caused by the vessel discreteness [7] which can Therefore as is seen in Fig. 2 a walker going into a large be described assuming the heat sink in equation (2) to arteryoflength l atinitialtime during the time τ before contain a random component [3]. Another is fast heat beingtrappedby theveinscanbe foundequiprobablyat transport over scales exceeding substantially the length each point of the given artery neighborhood of size l. In ℓ caused by the effect of blood flow through the artery T other words, this walker makes a large jump of length l tree. The latter phenomena is the main subject of the thatexceedssubstantiallythemean-fielddiffusionlength present paper. ℓ . T Inwhatfollowswewillanalyzethetemperaturedistri- II. FAST HEAT TRANSPORT WITH BLOOD bution averaged over all the possible realizations of the FLOW THROUGH LARGE ARTERY TREE vascular network embedding. This enables us to regard a walker entering a large artery of length l as a random Let at a certain time a random walker wandering in eventwhoseprobabilityis independent ofthe walkerini- the cellular tissue get a boundary of a large artery, i.e. tial position. Thereby, the probability Pl for a walker to an artery of length exceeding ℓ . It should be noted make a large jump over the distance l is also indepen- v that such an event is of low probability and cannot be dentofthespatialcoordinatesr. Itshouldbenotedthat consideredwithinthestandardmean-fieldapproximation forafixedrealizationofthe vascularnetworkembedding because the relative number of large arteries is small. the probability Pl depends essentially on the spatial co- Due to the direction of the blood motion from larger ordinatesrandthe heattransferinliving tissueonlarge arteries to smaller ones as well as the high blood flow scales has to exhibit substantial dependence on the spe- rate in the large arteries the walker will be transported cific position of the heat sources. fasttooneofthearteriesoflengthℓ (Fig.2). Theblood Now we estimate the value of P assuming the heat v l flow rate in smallvesselsis nothigh enoughto affect the sourcestobelocalizedinsideadomainQLofsize . Two L walker motion essentially and it has inevitably to leave differentfactorsdeterminethevalueofP . First,itisthe l this artery and wander in the cellular tissue until being process of walker trapping by an artery of length l > L trapped by the veins of length ℓv. Thereby a certain goingthroughthedomainQL. Ifl< bloodflowinthis L not too large number of random walkers generated, for arteryhas practically no effect onheat diffusion. On the 3 average a random walker during the time τ travels the distance ℓ in the cellular tissue until being trapped by T theveinsoflengthℓ . Soforawalkertoenterthisartery v and, thus, to leave the domain QL with blood flow in it thewalker,ononehand,shouldbelocatedatinitialtime inside a cylindrical neighborhood Q of the given artery l whose radius is about ℓ andthe volume is ℓ2. Onthe T L T other hand, it has to avoid being trapped by the veins of length ℓ . The probability of the latter event is about v (ℓ /ℓ )2. Indeed, aveinoflengthℓ maybe treatedasa v T v trap of cylindrical form. Thereby in qualitative analysis FIG. 3: Models of the peripheral artery network embedding the walker trapping can be described in terms of two- into the cellular tissue, (a) four-fold node model used in the dimensional random walks in the plane perpendicular to present analysis and (b) a more realistic dichotomic artery treeuniformlyembeddedintoacellulartissuedomainM. In the artery under consideration where the trapping veins the qualitative description of heat transfer both the models are represented by small circular regions. Their density lead to thesame result. is about 1/ℓ2 which directly leads to the latter estimate. v Thereforethetotalnumberofwalkersleavingthedomain QL with blood flow through the given artery per unit For the infinitely long cylinder time is P(2e) L 1 ℓ 2 D l ∼ l v ℓ2 (T T )= (T T ). (5) (cid:18) n(cid:19) τ · ℓ ·L T · − a L L − a (cid:18) T(cid:19) n and for the plane layer P(1e) 1. Multiplying g (r) by l ∼ l In obtaining (5) we have taken into account expres- the corresponding values of P(e) we get the result of av- l sion(3). Sincethetrappedwalkersspreaduniformlyover eraging the walker transition rate g (r) over the possible l aregionofsize l the resultingdensity ofthe walkertran- realizationsof the vascularnetwork embedding. The ob- sitionrateto apoint rspacedata distanceaboutl from tained result is written as the domain QL is D d g (r) L (T T ), (7) g (r) D L(T T ). (6) h l i∼ Lnln2+d − a l ∼ Lnl3 − a wherethevaluedactuallyplaystheroleofthedimension of the space inside which the temperature field can be It should be noted that the transition rate gl(r), as it considered. Atthenextstepweshouldsumtheterms(7) must, does not depend onthe localvalue of bloodperfu- over all the levels of the large artery tree. However due sion rate. to the strongincreaseofthe terms (7), g (r) 2n(2+d), l h i∝ At the second step we should average the obtained as the level number n increases the arteries of length transition rate gl(r) over the possible realizations of the l r mainly contribute to the value of gl(r) . So the ∼ h i vascular network embedding. Let us adopt a simplified term describing the fast heat transport with blood flow model for the vascular network shown in Fig. 3a where throughlargearteriesfromthe domainQL (locatednear the vessel lengths ln and ln+1 of the neighboring hierar- the origin of the coordinate system) can be written as chylevelsnandn+1arerelatedasl =2l . Figure3b n n+1 demonstrates a more adequate model for the peripheral g(r) D dr′ T(r′) , (8) aprrteesreynttrqeueawlithaicthiv,ehaonwaelvyesris, wmitahyinbethreedfruacmedewtoorkthoefftohre- ∼ Ln ZM |r−r′|5 mer one by combining three sequent two-fold nodes into where is the region containing the peripheral vascu- M one effective four-fold node at all the levels. In this case lar network as a whole and the integration in the three the cubic domain of volume l3 falls per each artery of dimensional space over the domain QL allows for all its n level n. Let us now consider individually three charac- three considered types. teristic forms of the domain QL, a ball or a cube of size Expression (8) together with the mean-filed bioheat (d=3), aninfinitely longcylinder ofradius (d=2), equation (2) enables us to write the following equation Land a plane layer of thickness (d=1). For tLhe ball or governing the anomalous heat transfer in living tissue L the cube, i.e. a region bounded in three dimensions the ∂T D T(r′) T(r) probability that an artery of level n passes through the =D 2T + dr′ − fj(T T )+q , domain QL is about ∂t ∇ Ln ZM |r−r′|5+ℓ5v − − a T (9) 2 P(3e) L . where we have added directly the value ℓv in order to l ∼(cid:18)ln(cid:19) cut off the spatial scales smaller than the length ℓv and 4 ignored the difference between the effective temperature givingusalsotheminimalsize oftheregionwherein mim L diffusivity and the true one of the cellular tissue. Equa- the tissue temperature increase(T T ) is mainly lo- max a − tion (9) is the desired governing equation of the anoma- calized. In the neighboring tissue the blood perfusion lous fast heat diffusion in living tissue for the averaged ratekeepsupasufficientlylowvaluej ,whichmakesthe 0 temperature field. It should be noted that the second heat propagation with blood flow through large arteries term on the right-hand side of equation (9) depends considerable. Indeed, let us estimate the temperature weakly on the blood perfusion rate. Therefore, for the increase caused by this effect using the obtained equa- nonuniformdistributionofthebloodperfusionratej(r,t) tion (9). The temperature increase T(r) T at a point a − equation (9) holds also. spaced at the distance r > from the region (of size L ) affected directly, i.e. the tail of the temperature dis- L tribution is mainly determined by the anomalous heat Anomalous temperature distribution under local diffusion and, so, is estimated by the expression strong heating Hyperthermia treatment as well as thermotherapy of smalltumors ofsize aboutorless1 cmisrelatedto local strong heating of living tissue up to temperatures about 45 ◦C or higher values. In this case the tissue region T(r) T jmax (ℓ∗T)2L3 (T T ) . (11) heated directly, for example, by laser irradiation is also − a ∼ Lnj0 r5 max− a ofasimilarsize. Duetothetissueresponsetosuchstrong heatingthe bloodperfusionratecangrowtenfoldlocally whereas in the neighboring regions the blood perfusion rate remains practically unchanged [9]. The feasibility of such nonuniform distribution of the blood perfusion Asseenfrom(11)forasufficientlylocalandstrongheat- ratemaybeexplainedapplyingtothecooperativemech- ing of the tissue, i.e. when ℓ∗ and j j the anismofself-regulationinhierarchicallyorganizedactive L ∼ T max ≫ 0 temperature increase at not too distant points such that media [3, 10]. Therefore in the region affected directly r & can be considerable. Otherwise the anomalous the blood perfusion rate jmax can exceed the blood per- heat Ldiffusion is ignorable. fusion rate j in the surrounding tissue substantially. In 0 Acknowledgments this case the characteristic length of heat diffusion into the surrounding tissue is about D ∗ ℓ , (10) This work was supported by STCU grant #1675 and T ∼sfjmax INTAS grant #00-0847. [1] M. M. Chen and K. R. Holmes, “Microvascular contri- [7] J.W.Baish,P.S.Ayyaswamy,andK.R.Foster.“Small- butions in tissue heat transfer”. Ann. N. Y. Acad. Sci., scale temperature fluctuations in perfused tissue during 335, 137–154 (1980). localhyperthermia”,ASMEJ. Biomech. Eng.108,246– [2] Heat Transfer in Medicine and Biology. Analysis and 250 (1986). Applications, A. Shitzer and R. C. Ebergard, Editors [8] B.B.Mandelbrot,TheFractalGeometryofNature(Free- (Plenum,New York,1970). man, New York,1977). [3] V.V. Gafiychuk and I.A. Lubashevsky. Mathemat- [9] C. W. Song, A. Lokshina, I. G. Rhee, M. Patten, and ical Description of Heat Transfer in Living Tis- S. H.Levitt. “Implication of blood flow in hyperthermia sue, (VNTL Publishers, Lviv, 1999); e-print: adap- treatment of tumors”. IEEE Trans. Biom. Eng., BME- org/9911001,9911002. 31 (1), 9–15 (1984). [4] G. I. Mchedlishvili. Microcirculation of Blood. General [10] I. A. Lubashevsky and V. V. Gafiychuk. “Cooperative Principles of Control and Disturbances (Nauka Publish- mechanism of self-regulation in hierarchical living sys- ers, Leningrad, 1989) (in Russian). tems”, SIAM, J. Appl. Math. 60(2), 633-663 (2000). [5] S. Weinbaum and L. M. Jiji. “A new simplified bioheat [11] J.C. Chato. “Heat transfer to blood vessels”, ASME J. equation fortheeffect ofblood flowon local averagetis- Biom. Eng. 102, pp.110–118 (1980). suetemperature”,ASMEJ.Biomech.Eng.107,131–139 [12] E.H. Wissler. “Comments on the new bioheat equation (1985). proposed by Weinbaum and Jiji”, ASME J. Biom. Eng. [6] S.Weinbaum,L.X.Xu,L.Zhu,andA.Ekpene.“Anew 109, pp.226–233 (1987). fundamental bioheat equation for muscle tissue: Part I– [13] E.H. Wissler. “Comments on Weinbaum and Jiji’s dis- Blood perfusion term”, ASME J. Biomech. Eng. 119, cussion of their proposed bioheat equation ”, ASME J. 278–288 (1997). Biom. Eng. 109, pp.355–356 (1988). 5 [14] I. A. Lubashevsky, A. V. Priezzhev, V. V. Gafiychuk, renormalizationofthebloodperfusionratecausedbythe and M. G. Cadjan. “Free-boundary model for local counter-currenteffect[11,12,13].Itstheoreticalestimate thermal coagulation”. In: Laser-Tissue Interaction VII, waslaterobtainedindependentlybyWeinbaumet al.[6] S.L. Jacques, Editor, Proc. SPIE 2681 81–91 (1996). andGafiychuk&Lubashevsky[3](announcedforthefirst [15] Initiallythefactor f wasphenomenologically introduced time in [14]). in the bioheat equation to take into account a certain