text:
# Laser Induced Heat Diffusion Limited Tissue Coagulation: Problem and General Properties

## 1 What the heat diffusion limited tissue coagulation is: a phenomenon and a theoretical problem

Thermal coagulation of living tissue caused by local heating due to laser light absorption is one of the novel thermotherapy techniques of tumor treatment which is currently undergoing clinical trials (see, e.g., ). Thermal coagulation is used to form a necrosis domain of desired form for the removal of the malignant tissue. So the mathematical modeling of the necrosis growth is required, first, to find out the physical limitations and the basic characteristics of the treatment and, second, to optimize the therapy course. However, living tissue is extremely complex in structure, thereby, for the adequate theoretical model to be developed and for the mathematical modeling of the given process to be implemented reliably typical limit cases should be singled out and studied individually. Such a case is the subject of the present paper.

When the tissue region affected directly by laser light is sufficiently small is size, the necrosis formation is mainly governed by heat diffusion. Namely, in this case we actually deal with the following physical process (Fig. 1). Absorption of laser light delivered into a small internal region of living tissue causes the temperature to attain such high values (about or above $`70^{}`$C) that lead practically to immediate thermal coagulation in this region. Heat diffusion into the surrounding live tissue causes its further coagulation, giving rise to the growth of the necrosis domain. In this case heat diffusion plays a significant role in the necrosis growth because the necrosis size $`\mathrm{}`$ exceeds the depth of laser light penetration into the tissue. Therefore, the temperature distribution inevitably has to be substantially nonuniform and for the tissue to coagulate at peripheral points heat diffusion should cause the temperature to grow at these points. The latter property singles out the specific mode of thermal coagulation under discussion from other possible types of thermotherapy treatment and that is why we refer to the necrosis growth under the given conditions as to thermal tissue coagulation limited by heat diffusion. In particular, as will be seen below, the optimal implementation of this thermotherapy mode is characterized by the formation of necrosis domains of size $`\mathrm{}5`$$`10`$ mm and by the treatment duration of $`t_{\text{course}}\mathrm{}^2/D2`$$`8`$ min (where $`D210^3`$ cm2/sec is the tissue temperature diffusivity). Laser lighted is considered to be absorbed within a layer of thickness less than $`\mathrm{}`$.

Heat diffusion in the live tissue is affected substantially by blood perfusion causing the heat sink . Thus, the characteristics of the spatial distribution and the dynamics of the blood perfusion rate should have a substantial effect on the necrosis growth limited by heat diffusion. Therefore, in modeling mode one has to take into account the tissue response to the temperature growth which can locally give rise to a tenfold increase in the blood perfusion rate . The latter effect, in particular, is responsible for a substantially nonuniform distribution of the blood perfusion rate and visually manifests itself in a red ring (“hyperemic ring”) appearing around the necrosis region during the thermotherapy treatment. Besides, when heated the living tissue will inevitably exhibit spatial nonuniformities in the temperature due to the vessel discreteness . Because of the extremely strong temperature dependence of the thermal coagulation rate such temperature nonuniformities can substantially perturb the necrosis form and, so, this effect should be also taken into account.

In the previous papers  basing on the free boundary description we have developed a model for the heat diffusion limited thermal tissue coagulation and studied the properties of the corresponding necrosis growth. This model has been developed in order to take into account the effects mentioned above. The aim of the present paper is to survey the obtained results, to outline the general properties of the necrosis growth under the aforementioned conditions, and to justify the point of view on the laser induced heat diffusion limited tissue coagulation as a individual mode of laser therapy.

However, before stating the free boundary model let us, first, recall the main features of living tissue affecting heat transfer, discuss the approach previously used by other authors and analyze the corresponding problems in the theoretical description of the local thermal coagulation. In this way we will make the key points of the theory we have developed more clear.

### 1.1 Background: the main features of bioheat transfer

#### 1.1.1 Living tissue as a heterogeneous medium

Blood flowing through vessels in living tissue forms paths of fast heat transport and under typical conditions it is blood flow that governs heat propagation on scales exceeding several millimeters (for an introduction to this problem see, e.g., ). Blood vessels make up a complex network being practically a fractal. The larger is a vessel, the faster is the blood motion in it and, so, the stronger is the effect of blood flux in the given vessel on heat transfer. The blood flux in the smallest vessels, capillaries, practically does not affect heat propagation . Thus, there should be blood vessels of a certain length $`\mathrm{}_v`$ that are the smallest among the vessels wherein the blood flux affects heat transfer remarkably. The value of $`\mathrm{}_v`$ may be different at various points (under nonuniform blood perfusion) and can be found from the expression  (see also ):

$$\mathrm{}_v(𝐫)\sqrt{\frac{\kappa }{c\rho fj_v(𝐫)L_n}}.$$
(1)
Here $`\kappa `$ is the thermal conductivity of the cellular tissue, $`c`$, $`\rho `$ are the density and heat capacity of the tissue, $`j`$ is the blood perfusion rate (the volume of blood going through tissue region of unit volume per unit time) and $`j_v`$ is its value averaged over scales of order of $`\mathrm{}_v(𝐫)`$. The factor $`f<1`$ accounts for the counter-current effect. Initially it was phenomenologically introduced in the bioheat equation as a certain renormalization of the blood perfusion rate  and its theoretical estimate will be discussed below in the present section. The factor $`L_n`$ playing a key role in the theory of bioheat transfer is given by the formula:

$$L_n=\mathrm{ln}\left(\frac{l}{a}\right),$$
(2)
where $`l/a`$ is the mean ratio of the individual length to radius of blood vessels forming peripheral systems of blood circulation. Since expression (1) estimating the length $`\mathrm{}_v`$ contains the blood perfusion rate $`j_v`$ itself averaged over the scale $`\mathrm{}_v`$, it is implicit and has to be completed with an expression relating the value $`\mathrm{}_v(𝐫)`$, the averaged perfusion rate $`j_v(𝐫)`$ and the true one $`j(𝐫)`$. At first, we write it in a symbolic form

$$\widehat{}\{\mathrm{}_v,j_v,j\}=0.$$
(3)
Nevertheless, to estimate the length $`\mathrm{}_v`$ we may ignore the difference between the averaged and true blood perfusion rates in expression (1). Then for the typical values of the ratio $`l/a20`$–40 , the thermal conductivity $`\kappa 710^3`$ W/cm$``$K, the heat capacity $`c3.5`$ J/g$``$ K, and the density $`\rho 1`$ g/cm3 of the tissue, as well as setting the blood perfusion rate $`j_vj0.3`$ min-1 and the factor $`f0.5`$ we get from (1) and (2):

$$\mathrm{}_v4\mathrm{mm}\mathrm{and}L_n34.$$
(4)
The theory of heat transfer in living tissue (bioheat transfer theory) starts from the system of microscopic equations governing heat propagation in the cellular tissue and with blood inside the vessels individually. However, on one hand, the complex structure of the vascular network and, on the other hand, the lack of the detailed information about the arrangement of individual vessels necessitate the development of the macroscopic description of heat transfer in living tissue regarded as an effective continuum with certain, may be, anomalous properties. Moreover, the vascular network can vary in particular details from tissue to tissue or even from patient to patient for one tissue. In this case the development of an adequate macroscopic theory is the only way to model a thermotherapy treatment. The macroscopic governing equation for the tissue temperature can be obtained by averaging the microscopic equations over certain scales that, first, are small enough so the tissue temperature does not exhibit great variation on these scales. Second, these scales should be large enough for the heat exchange between the cellular tissue and blood could be treated within a continuous approximation. This heat exchange is directly controlled by the vessels of length $`\mathrm{}_v`$. So the length $`\mathrm{}_v`$ is a natural scale of the averaging in the theory of bioheat transfer. Besides, it turns out that the characteristic scale of the temperature variations is

$$\mathrm{}_T\sqrt{L_n}\mathrm{}_v\sqrt{\frac{\kappa }{c\rho fj}}.$$
(5)
The latter estimate demonstrates us that the value $`1/L_n`$ plays the role of small parameter in the bioheat transfer theory, exactly which has made it feasible to construct a mean field approximation regarding living tissue as a continuum .

In particular, averaging the microscopic equations in this way, we have got the following macroscopic bioheat equation relating the tissue temperature $`T_v`$ averaged over scales of order of $`\mathrm{}_v`$ and the averaged perfusion rate $`j_v`$:

$$c\rho \frac{T_v}{t}=\left(\kappa _{\text{eff}}T_v\right)fc\rho j_v(T_vT_a)+q_v,$$
(6)
where we have ignored the difference in the density and heat capacity of the tissue and blood, respectively, $`cc_b`$, $`\rho \rho _b`$, the blood temperature $`T_a`$ in large arteries of systemic circulation is considered to be fixed, and $`q_v`$ is the heat generation rate averaged over the same scales. The factor $`f`$, accounting for the counter-current effect, and the ratio $`F=\kappa _{\text{eff}}/\kappa `$ are estimated by the expressions:

$$f\frac{1}{\sqrt{L_n}}0.5,F=F(L_n)2.$$
(7)
It should be pointed out that the factors $`f`$ and $`F`$ are determined solely by the geometry of the vascular network, which directly results from the vascular network being fractal in structure 111The theoretical estimate of $`f0.5`$ has been obtained for the first time in book , announced in our papers  and has been later presented independently in paper . The effect of blood perfusion on the heat propagation rate in terms of the heat conductivity renormalization was quantitatively studied actually for the first time in paper . Moreover, due to the logarithmically weak dependence (2) these factors can be approximately treated as constants determined beforehand.

It should be noted that blood flow through the large arteries can give rise to an anomalously fast heat transport over scales much greater than $`\mathrm{}_T`$. Under certain conditions it may appear that the effective thermal conductivity $`\kappa _{\mathrm{eff}}`$ exceeds the thermal conductivity of the cellular tissue by tenfold due to this effect. However, such a fast heat transport cannot be described in terms of the mean field theory and deserves an individual consideration. Besides, on the average, its role is not too essential because of the sufficiently small relative volume of the large arteries. Therefore dealing with local thermal coagulation we can ignore this effect at the first approximation. We also point out that equation (6) is similar in its form to the phenomenological generalized bioheat equation  except for the latter contains the true perfusion rate $`j`$ rather than the averaged one $`j_v`$. This fact can be essential for substantially nonuniform blood perfusion  as it is the case in local thermal coagulation .

The mean field approximation leading to equation (6) does not allow for the temperature nonuniformities caused by the vessel discreteness  because it regards living tissue as a certain continuum. So, to take into account these nonuniformities we have to go beyond the scope of the mean field theory. Since the particular details of the vessel arrangement on scales about several millimeters are practically unknown and, moreover, alter in various tissues and, may be, at different points of one tissue it is reasonable to regard the vessel arrangement and the corresponding temperature nonuniformities as random . In particular, under uniform heating, i.e. when $`T_v`$ is constant, these nonuniformities will be characterized by the mean amplitude $`\sigma `$ and the correlation length $`\lambda `$:

$`\delta T(𝐫,t)\delta T(𝐫,t)`$ $`=`$ $`\sigma ^2\mathrm{and}`$ (8)
$`\delta T(𝐫,t)\delta T(𝐫^{},t)`$ $``$ $`\sigma ^2\mathrm{for}\left|𝐫𝐫^{}\right|\lambda ,`$
where the symbol $`\mathrm{}`$ stands for averaging. As might be expected, the correlation length $`\lambda `$ is about $`\mathrm{}_v`$ ($`\lambda \mathrm{}_v`$) and depends on the local value of the perfusion rate. By contrast, due to the fractal structure of the vascular network the ratio of the amplitude $`\sigma `$ to the overheating $`T_vT_a`$ is estimated as :

$$\frac{\sigma }{T_vT_a}\frac{1}{L_n}$$
(9)
and can be regarded as a predetermined constant. Taking into account additional numerical factors  we get $`\sigma (0.1`$$`0.2)(TT_a)`$ (see also ). Since the vessel discreteness effect is mainly caused by the vessels of lengths about $`\mathrm{}_v`$ it can be described in terms of a bioheat equation similar to (6) which, however, contains the true tissue temperature $`T`$ and where the averaged perfusion rate $`j_v`$ is replaced by term:

$$j_vj_v+\delta j(𝐫),$$
(10)
involving the random component $`\delta j(𝐫)`$ meeting the conditions:

$`\delta j(𝐫)`$ $`=`$ $`0,`$ (11)
$`\delta j(𝐫)\delta j(𝐫^{})`$ $`=`$ $`j_v^2g^0\left({\displaystyle \frac{\left|𝐫𝐫^{}\right|}{\mathrm{}_v}}\right).`$ (12)
Hear the function $`g^0(x)`$ is such that $`\left|g^0(x)\right|1`$ for $`x1`$ and $`g^0(x)1`$ for $`x1`$.

In this subsection we have discussed the characteristic features of living tissue as a heterogeneous medium. However it is also an active medium because it responds to temperature variations by increasing the blood vessel radius, which gives rise to an increase in the blood perfusion rate $`j`$. The blood perfusion rate can locally grow by tenfold , so, this effect is significant and is the subject of the next subsection.

#### 1.1.2 Living tissue as an active medium

Living tissue tries to keep its temperature from exceeding a certain vital boundary $`T_+44`$$`46`$ C ($`T<T_+`$) to prevent thermal damage (see, e.g., ). So, the main increase in the blood perfusion rate $`j`$ should fall on the temperature variations from $`T_a`$ to a certain value $`T_{\text{vr}}T_+`$ and after the temperature exceeds the value $`T_{\text{vr}}`$ the blood perfusion rate $`j`$ is likely to depend weakly on temperature because the blood vessels exhaust their ability to expand.

In order to analyze the temperature evolution in living tissue the dependence $`j(T)`$ obtained experimentally under practically uniform heating of large tissue regions is typically used. However, whether such a dependence holds when the temperature distribution becomes substantially nonuniform is a question. Indeed, in this case due to the vessel extent the dependence $`j(T)`$ can alter not only its particular form but also become functional. In other words, when the tissue is heated substantially nonuniform the blood perfusion rate $`j(𝐫)`$ at a given point $`𝐫`$ is determined, in general, by the whole temperature distribution $`\{T(𝐫^{})\}`$ over a certain neighborhood of this point (concerning a similar behavior of the tissue response to variations in CO2 concentration see, e.g., ). Specific details of the mechanism by which living tissue responds to local temperature variations on scales about 1 cm is also a question . However, under strong heating localized on such scales thermal selfregulation can be effectively implemented through the response of microcirculatory bed to reduced O2 partial pressure or increased concentration of metabolism products (e.g. CO2, H+, ADP, etc.) because higher temperatures result in a metabolism intensification with a higher O2 consumption. In addition, this explains a possible delay of the tissue response to temperature variations and enables the corresponding delay time to be estimated at several minutes (see, e.g., ).

In order to describe the local thermoregulation we need to specify how each vessel responds to the corresponding piece of information characterizing the state of the tissue, in particular, its temperature. So, in deriving the equation relating the blood perfusion rate $`j(𝐫,t)`$ to the tissue temperature $`T(𝐫,t)`$ we make use of the following physiological data . The local selfregulation of blood perfusion in living tissue on scales about one centimeter is due to expansion or contraction of blood vessels making up a single microcirculatory bed and blood redistribution over this vascular network is mainly controlled by a large group of arteries differing in length significantly. The reaction of the microcirculatory bed is governed by receptors responding to variations in CO2 partial pressure, H+ concentration, or other metabolism products. However, as mentioned above, under strong heating such receptors can play the role of effective thermoreceptors supplying the microcirculatory bed with the needed information. Obviously, none of the vessels receives the whole information on the tissue state, so there must be a certain cooperative mechanism of information self-processing by which the behavior of different vessels is consistent with each other in the manner enabling the tissue to respond properly.

We have shown  that such a cooperative mechanism of self-regulation can be implemented through the vessel response to the blood temperature in the corresponding veins. The receptors mentioned above are located directly in the cellular tissue as well as embedded into the walls of vessels, including veins . Those embedded into the vessel walls are able to measure the concentration of the metabolism products directly in blood and, thus, to measure effectively its temperature. For small vessels (arteries and veins) the position of the receptors governing their behavior is not a factor. This allows us to make use of the proposed cooperative mechanism of self-regulation in the description of tissue response to local and strong heating.

Under the adopted assumptions it turns out that for the normal tissue the dependence of blood perfusion rate $`j`$ on the tissue temperature $`T`$ can be approximately described by a local equation relating the values of $`j(𝐫)`$ and $`T(𝐫)`$ taken at the same point $`𝐫`$ until the temperature comes close to the vital boundary $`T_+`$. This is due to the cooperative mechanism of selfregulation which involves the individual response of each vessel to the corresponding piece of information and the coordination of their behavior by the selfprocessing of information. The adequate selfprocessing of information is implemented through heat conservation during blood motion inside relatively large veins of the microcirculatory bed. The vascular network whose behavior is governed by this mechanism can supply each region of the cellular tissue with blood at such a rate that is required of its individual demand and different regions of the cellular tissue do not interfere with one other. The obtained equation for the tissue response is of the form :

$$\tau _v\frac{j}{t}+j\frac{T_+T}{T_+T_a}=j_0$$
(13)
where $`\tau _v`$ is the delay time of tissue response and $`j_0`$ is the blood perfusion rate provided the tissue is not affected. This result practically holds for living tissue containing also a small necrosis domain  because the temperature of blood in veins whose length exceeds its size by several times is not sensitive to the presence of necrosis domain.

In the next subsection we analyze the models of local thermal coagulation used previously by other authors and discuss the basic theoretical problems met in such approaches.

### 1.2 Distributed description of thermal coagulation

The effect of heat diffusion on local thermal coagulation has been considered and numerically simulated by a number of authors (see, e.g.  and references therein). Leaving aside the description of laser light propagation in the tissue we can generalize their models for the necrosis formation to the following coupled equations for the temperature field $`T(𝐫,t)`$ and the field $`\zeta (𝐫,t)`$ determining the fraction of undamaged tissue at a given point $`𝐫`$ and time $`t`$:

$`c\rho {\displaystyle \frac{T}{t}}`$ $`=`$ $`\left(\kappa _{\text{eff}}T\right)fc_b\rho _bj(TT_a)+q,`$ (14)
$`{\displaystyle \frac{\zeta }{t}}`$ $`=`$ $`\zeta \omega (T).`$ (15)
The thermal coagulation rate $`\omega (T)`$ depends strongly on the tissue temperature $`T`$ and for typical values of temperature attained during the treatment course can be approximated by the expression

$$\omega (T)=\omega _0\mathrm{exp}\left[\frac{TT_0}{\mathrm{\Delta }}\right].$$
(16)
Here $`\omega _0=\omega (T_0)`$, where $`T_0`$ is a certain fixed temperature and $`\mathrm{\Delta }`$ is a constant. Expression (16) is actually a convenient approximation of the Arrhenius dependence $`\omega (T)\mathrm{exp}\{\frac{E}{T}\}`$ inside a not wide temperature interval governing the tissue coagulation. The available experimental data  for the temperature dependence of the exposure time enable us to estimate the value of $`\mathrm{\Delta }`$ as $`\mathrm{\Delta }35^{}`$C ($`\mathrm{\Delta }3.6^{}`$C for pig liver at $`T_0=65^{}`$C). Below in numerical calculations the dependence (16) will be taken in the form

$$\omega (T)=0.2\mathrm{exp}\left[\frac{T60}{3.6}\right]\text{min}^1,$$
(17)
where the temperature $`T`$ is in degrees Celsius.

In what follows the system of equations (14), (15), and expression (16) will be referred to as the “distributed model” for local thermal coagulation. This model in fact is rather a phenomenological one and should be modified in order to allow for sophisticated effects revealing themselves in local coagulation.

Let us consider the properties of local thermal coagulation that, in fact, make the mathematical description of this process stubborn and give rise to the inconsistency of the distributed model. During a typical course of thermotherapy treatment high temperatures of order of $`T_{\text{max }}100^{}`$C are attained at the necrosis center, whereas at distant points $`T=T_a37^{}`$C. Under such conditions the damage rate $`\omega (T)`$ or its reciprocal value $`t_{\text{cg}}(T)=1/\omega (T)`$ called the threshold exposure time of thermal coagulation varies extremely sharp in space on scales of the necrosis size $`\mathrm{}`$. Indeed, at points where the temperature attains, for example, the values of $`T=60^{}`$C, $`65^{}`$C, $`70^{}`$C, and $`75^{}`$C the values of $`t_{\text{cg}}`$ are $`8`$ min, $`100`$, $`20`$, and $`5`$ sec, respectively . The typical duration $`t_{\text{course}}`$ of the thermotherapy course is about several minutes and estimated by the inequality, $`t_{\text{course}}\tau `$ , where

$$\tau \stackrel{\text{def}}{=}\frac{1}{fj}36\text{min},$$
(18)
is the characteristic time scale of the necrosis growth and $`j`$ is the mean value of the blood perfusion rate attained near the necrosis boundary. In estimate (18) we have set $`f=0.5`$ and used the typical values of the blood perfusion rate $`j0.3`$, $`0.5`$, and $`0.7`$ min-1 for stomach, intestine, and spleen, respectively . So the layer wherein the thermal coagulation is under way at a given moment of time is characterized by a narrow temperature interval of thickness $`\mathrm{\Delta }`$ located in the vicinity of the temperature of $`T_{\text{cg}}65^{}`$C. This layer $`𝕃_{\text{pd}}`$ of partially damaged tissue separates the necrosis region ($`\zeta 1`$) and practically live tissue ($`\zeta 1`$ ). Thus, in order to describe the necrosis growth due to local thermal coagulation in the framework of the distributed model one has to consider in detail the temperature field inside the layer $`𝕃_{\text{pd}}`$, because exactly it governs the necrosis growth directly.

Now we estimate the thickness $`\delta _{\text{pd}}`$ of this layer. As follows from the numerical analysis  the temperature distribution is characterized by a single spatial scale, so, in the necrosis domain the mean temperature gradient is about $`T(T_{\text{max}}T_a)/\mathrm{}`$. In addition, the necrosis size attaining typically the values of $`1020`$ mm meets the inequality :

$$\mathrm{}\mathrm{}_T10\text{mm}$$
(19)
by virtue of (4) and (5). Then setting $`\delta _{\text{pd}}T\mathrm{\Delta }`$ we get the desired estimate:

$$\delta _{\text{pd}}\frac{\mathrm{\Delta }\sqrt{L_n}}{(T_{\text{max}}T_a)}\mathrm{}_v0.2\mathrm{}_v1\text{mm }.$$
(20)
So the layer $`𝕃_{\text{pd}}`$ of partially damaged tissue is sufficiently thin, its thickness is less not only than the mean necrosis size $`\mathrm{}`$ but also, what is essential for the theory development, than the characteristic length $`\mathrm{}_v`$ of the vessels directly governing the heat exchange between the cellular tissue and blood. Therefore, in order to describe rigorously the temperature distribution on scales about $`\delta _{\text{pd}}`$ we have to go beyond the framework of the mean field approximation. However, the distributed model is, in fact, a mean field one, because it regards living tissue in terms of an effective continuum. Thereby the distributed model in its own right cannot ensure the reliability for the mathematical modeling of local tissue coagulation.

In particular, as an aspect of this problem we note that the mean field approximation ignores the difference between the true temperature distribution in the tissue and the temperature field $`T(𝐫,t)`$, i.e. does not take into account the random component of temperature distribution caused by the vessel discreteness. This is justified if the mean amplitude $`\sigma `$ of such temperature nonuniformities is sufficiently small. By virtue of (4) and setting the temperature inside the layer $`𝕃_{\text{pd}}`$ of partially damaged tissue equal to $`T_{\text{cg}}65^{}`$C we get the estimate

$$\sigma |_{𝕃_{\text{pd}}}\mathrm{\Delta }.$$
(21)
In other words, inside the layer $`𝕃_{\text{pd}}`$ the amplitude of the temperature nonuniformities due to the vessel discreteness turns out to be about the characteristic temperature drop across the whole layer $`𝕃_{\text{pd}}`$. So such temperature nonuniformities are bound to affect the necrosis formation.

Another problem inherent in the distributed model is the specification of the main kinetic coefficients characterizing the partially damaged tissue as a continuous medium that effectively approximates the real cellular tissue with the embedded vascular network, with both of them having undergone partial damage. Namely, we should specify the effective thermal conductivity $`\kappa _{\text{eff}}`$ and the factor $`f`$ as functions of the blood perfusion rate $`j`$ and the fraction $`\zeta `$ of undamaged tissue:

$$\kappa _{\text{eff}}=\kappa _{\text{eff}}(\zeta ,j),f=f(\zeta ,j).$$
(22)
However, these coefficients being the characteristics of the mean field theory can adequately describe heat transfer in living tissue only on scales about or greater than $`\mathrm{}_v`$. So we may write the corresponding relationships only for the live tissue where $`\zeta 1`$. (Clearly, inside the necrosis domain $`\kappa _{\text{eff}}=\kappa `$ and we may set $`f`$ equal to any value of order unity because of $`j=0`$ in this region.) Inside the layer of partially damaged tissue any dependence of the type $`\kappa _{\text{eff}}=\kappa _{\text{eff}}(\zeta )`$ or $`f=f(\zeta )`$ is nothing more than a formal phenomenological approximation of a real complicated phenomenon. In particular, as follows from the results to be obtained the particular details of the dependence $`f=f(\zeta )`$ are of no concern and below we will treat the value $`f`$ as a constant.

Besides, in order to complete the distributed model one should describe how the blood perfusion rate depends on the undamaged tissue fraction $`\zeta `$ and the temperature. It has been proposed  to set

$$j(\zeta ,T)=\zeta j_t(T),$$
(23)
where $`j_t(T)`$ is the perfusion rate that would occur in tissue without damage. Clearly, it is also a pure phenomenological approximation.

The only fact, that is able to justify the feasibility of applying the distributed model to the analysis of local thermal coagulation can be the independence of the necrosis growth from the particular details of heat transfer in the layer of partially damaged tissue. In this case, however, it would be more consistent to use a free boundary model which ignores the thickness of this layer, i.e. treats it as the boundary of the necrosis domain. The motion of such an interface must be governed by the boundary values of the temperature and its gradient. Particular details of heat transfer in the given layer may be taken into consideration by a certain collection of parameters.

In the next section we will show that this independence is the case and derive the corresponding free boundary model. In particular, formally assuming equation (14) to hold in the whole region under consideration we will reduce the system of equations (14), (15) to a certain equivalent free boundary model dealing with integrated characteristics of the tissue damage rate inside the layer $`𝕃_{\text{pd}}`$. In other words, the model to be obtained aggregates the particular details of the functions $`\kappa _{\text{eff}}(\zeta )`$ and $`j(\zeta ,T)`$ to certain contains of order unity, which demonstrates the desired independence.

## 2 Interface dynamics

In order to obtain the free boundary description of local thermal coagulation it is suffice to consider a certain neighborhood $`Q`$ of the partially damaged tissue layer $`𝕃_{\text{pd}}`$ (Fig. 2). The size of the region $`Q`$ is assumed to be, on one hand, much smaller than the characteristic size $`\mathrm{}`$ of the necrosis domain and, on the other hand, much larger than the thickness $`\delta _{\text{pd}}`$ of the layer $`𝕃_{\text{pd}}`$. In other words, we analyze the case when the tissue in the region affected directly by laser light has already coagulated and the necrosis growth is caused by the heat diffusion into the surrounding tissue. To describe the dynamic of local thermal coagulation we introduce the interface $`\mathrm{\Gamma }`$ specified by the condition

$$\zeta (𝐫,t)|_{𝐫\mathrm{\Gamma }}=\zeta _0,$$
(24)
where $`\zeta _00.5`$ is a certain fixed value, and keep track of how it moves. Then the necrosis growth is represented by the motion of the interface $`\mathrm{\Gamma }`$ at the normal velocity $`\vartheta _n`$.

There are two ways of deriving the desired expression for the velocity $`\vartheta _n`$. One of them is as follows. Let us, first, assume equations (14), (15) to hold at all the points of the tissue. Then the derivation is to reduce the distributed model to certain equations dealing with the layer $`𝕃_{\text{pd}}`$ of partially damaged tissue in terms of an interface $`\mathrm{\Gamma }`$ whose velocity $`\vartheta _n`$ is determined by the local characteristics of the temperature field that practically do not vary over this layer and so remains unchanged at the nearest points of the live tissue and the necrosis region. This way has been implemented in our papers  using the singular perturbation technique in the small parameter $`\mathrm{\Delta }/(T_{\text{max}}T_a)`$. In particular, we have obtained the following expression relating the velocity $`\vartheta _n`$ of the interface $`\mathrm{\Gamma }`$, the value $`T_{\text{cg}}`$ of the temperature at this interface, and the boundary value of the temperature gradient, for example, on the necrosis side $`_nT|_{\mathrm{\Gamma }0}`$:

$$\vartheta _n=\mathrm{}_0\frac{\mathrm{\Delta }}{\left|_nT\right|_{\mathrm{\Gamma }0}}\omega (T_{\text{cg}}),$$
(25)
where the function $`\omega (T)`$ is the thermal coagulation rate given by formula (16) and the constant $`\mathrm{}_0`$ takes into account the formal dependence of the effective thermal conductivity $`\kappa _{\text{eff}}(\zeta )`$ on the undamaged tissue fraction $`\zeta `$

$$\frac{1}{\mathrm{}_0}=\underset{\zeta _0}{\overset{1}{}}𝑑\zeta \frac{\kappa }{\kappa _{\text{eff}}(\zeta )\zeta }1.$$
(26)
The other way of obtaining expression (25) is to state it basing on the general properties of thermal coagulation and without using the distributed model at all. Doing so we are not be able to get an expression for the constant $`\mathrm{}_0`$ and, thereby, should regard it as a phenomenological one. However, first, expression (26) is as rigorous as the former approach using the dependence $`\kappa _{\text{eff}}(\zeta )`$. Second, following the latter way we show that the free boundary model based on expression (25) is more general than the distributed model, i.e. it is not restricted to the validity area of the distributed model.

Let thermal coagulation be under way in a layer $`𝕃_{\text{pd}}(t)`$ at a given moment of time $`t`$ and the mean temperature and the mean temperature gradient inside this layer be $`T_{\text{cg}}`$ and $`_nT_{\text{pd}}`$, respectively. Then the characteristic thickness $`\delta _{\text{pd}}`$ of the layer $`𝕃_{\text{pd}}`$ can be estimated by expression

$$\delta _{\text{pd}}\frac{\mathrm{\Delta }}{_nT_{\text{pd}}}.$$
(27)
In fact, when the necrosis growth due to thermal coagulation is limited by the heat diffusion the temperature distribution is substantially nonuniform in space. Therefore, on one hand, in the necrosis region, where the temperature exceeds $`T_{\text{cg}}`$ by values greater than $`\mathrm{\Delta }`$: ($`TT_{\text{cg}}\mathrm{\Delta }`$) and, thus, $`\omega (T)\omega (T_{\text{cg}})`$, the tissue coagulation has to be completed. On the other hand, in the region of undamaged tissue the temperature is sufficiently low, $`T_{\text{cg}}T\mathrm{\Delta }`$, and, so, $`\omega (T)\omega (T_{\text{cg}})`$, thereby the tissue has not enough time to coagulate at such temperatures. Therefore tissue coagulation can be under way only in the region where $`\left|TT_{\text{cg}}\right|\mathrm{\Delta }`$ whence we immediately get estimate (27). After a lapse of the time interval $`t_{\text{ thr}}1/\omega (T_{\text{cg}})`$ the tissue in the layer $`𝕃_{\text{pd}}`$ has to coagulate practically completely. This is equivalent to the displacement of the layer $`𝕃_{\text{pd}}`$ over the distance $`\delta _{\text{pd}}`$. Thereby, if we keep watch on the points at which $`\zeta 0.5`$ then we will see that these points move at the velocity

$$\vartheta \frac{\delta _{\text{pd}}}{t_{\text{ thr}}}\frac{\mathrm{\Delta }}{_{\text{pd}}T}\omega (T_{\text{cg}}),$$

which exactly coincides with expression (25) within a factor of order unity.

The other relation completing the free boundary approximation of the necrosis growth reflects the heat flux conservation through the layer $`𝕃_{\text{pd}}`$ of partially damaged tissue. Obviously, it can be represented in the form

$$\kappa _nT|_{\mathrm{\Gamma }0}=\kappa _{\text{eff}}_nT|_{\mathrm{\Gamma }+0}$$
(28)
relating the temperature gradients $`_nT|_{\mathrm{\Gamma }0}`$ and $`_nT|_{\mathrm{\Gamma }+0}`$ at the interface $`\mathrm{\Gamma }`$ on the necrosis side and the side of live tissue, respectively.

Expressions (25), (28) are actually the essence of the free boundary model for the dynamics of local thermal coagulation. It should be pointed out that the stated free boundary approximation in its turn can be regarded as the initial point of modeling the necrosis growth due to local thermal coagulation. Indeed, all the information required of specifying the dependence $`\vartheta _n(T)`$ can be obtained based on the experimental data for the temperature dependence of the threshold exposure time at fixed temperature. The only parameter of this model that contains the particular information about the properties of heat transfer in the real layer of partially damaged tissue is the numeric factor $`\mathrm{}_0`$ of order unity, $`\mathrm{}_01`$. Keeping the latter in mind we may regard the analysis presented in the given section as the substantiation of the fact that the necrosis growth due to local thermal coagulation limited by heat diffusion is insensitive to the particular details of heat transfer inside partially damaged tissue. Thus, it actually justifies the distributed model rather than the free boundary model.

Now the obtained results (25), (28) enable us to state the desired free boundary model for heat diffusion limited thermal coagulation.

## 3 Free boundary model

Following the previous section the tissue is considered as comprising two regions (Fig. 3), the necrosis domain $`Q_\text{n}`$ where the blood perfusion rate is equal to zero,

$$j(𝐫,t)=0\text{for}𝐫Q_\text{n},$$
(29)
and the undamaged tissue $`Q_\text{t}`$ responding to temperature variations by increasing locally the blood perfusion rate $`j(𝐫,t)`$. Besides, we allow for that the tissue response can be delayed. Blood vessels can expand only to a certain extent as the temperature grows. So when it becomes high enough, $`T>T_{\text{vr}}`$, the blood perfusion rate $`j(𝐫,t)`$ attains a large but finite value $`j_{\text{max}}`$ and remains approximately constant. Taking into account expression (13) we describe this behavior of the undamaged tissue by the equation:

$$\tau _v\frac{j}{t}+j\mathrm{\Phi }(T)=j_0\text{for}𝐫Q_\text{t}.$$
(30)
Here $`\tau _v`$ is the delay time of the tissue response and the function $`\mathrm{\Phi }(T)`$ is of the form

$$\mathrm{\Phi }(T)=\{\begin{array}{cc}\alpha +(1\alpha )\frac{T_{\text{vr}}T}{T_{\text{vr}}T_a}\hfill & \text{for}T<T_{\text{vr}}\hfill \\ \alpha \hfill & \text{for}T>T_{\text{vr}}\hfill \end{array}$$
(31)
where $`\alpha =j_0/j_{\text{max}}`$ is the ratio of $`j_0`$ and the maximum $`j_{\text{max}}`$ of the blood perfusion rate that can be attained in living tissue due to the vessel expansion caused by the tissue response to temperature increase, and $`T_{\text{vr}}45`$$`46^{}`$C is the temperature at which the blood vessels exhaust their ability to expand. Inside the necrosis domain $`Q_\text{n}`$ the tissue temperature obeys the heat diffusion equation for solids:

$$c\rho \frac{T}{t}=\kappa ^2T+q$$
(32)
where $`\kappa `$ is the cellular tissue conductivity and $`q`$ is the rate of heat generation due to the laser light absorption. Inside the undamaged tissue the temperature is governed by the equation taking into account also the effect of the vessel discreteness (replacement (10)):

$$c\rho \frac{T}{t}=F\kappa ^2Tfc_b\rho _b(j_v+\delta j)(TT_a)+q.$$
(33)
Here, as before, $`j_v`$ is the regular component of the blood perfusion rate averaged over spatial scales of order $`\mathrm{}_v`$, the value $`\delta j`$ is its component due to the vessel discreteness treated as random, and the constants $`F=\kappa _{\text{eff}}/\kappa >1`$ and $`f<1`$ are of order unity. At the interface $`\mathrm{\Gamma }`$ between the necrosis domain and the undamaged tissue the temperature and the heat flux are assumed to have no jumps, i.e. the temperature distribution meets the following boundary conditions

$$T|_{\mathrm{\Gamma }+0}=T|_{\mathrm{\Gamma }0}\stackrel{\text{def}}{=}T_{\text{cg}},F_nT|_{\mathrm{\Gamma }+0}=_nT|_{\mathrm{\Gamma }0}.$$
(34)
By virtue of (16) and (25) the normal velocity of the interface $`\mathrm{\Gamma }`$ is given by the expression

$$\vartheta _n=\frac{\mathrm{}_0\omega _0\mathrm{\Delta }}{\left|_nT\right|_{\mathrm{\Gamma }0}}\mathrm{exp}\left[\frac{T_{\text{cg}}T_0}{\mathrm{\Delta }}\right]$$
(35)
For points distant from the necrosis domain

$$T|_{\mathrm{}}=T_a.$$
(36)
Finding the relationship between the averaged and true blood perfusion rates, $`j_v`$ and $`j`$, we have to take into account that the scale $`\mathrm{}_v`$ of averaging in its turn depends on the local value of $`j_v`$ as demonstrated by expression (1). This dependence enables us to specify functional (3) in the form :

$$j_v\frac{\lambda _v\kappa }{c\rho }^2\mathrm{ln}j_v=j,$$
(37)
where $`\lambda _v1/\sqrt{L_n}`$ is also a constant of order unity. Equation (37) should be subjected to a certain boundary condition at the interface $`\mathrm{\Gamma }`$ because it makes no sense to average the blood perfusion rate over the necrosis domain impermeable to blood. The physical layer separating the necrosis domain and the undamaged tissue where the local vascular network is not damaged is complex in structure and contains a spatial increase of the blood perfusion rate from zero to the value in the undamaged tissue. In order to avoid the problem of analyzing the blood perfusion rate in this layer we take into account the following simplifying circumstance. On one hand, the typical size of the necrosis domain formed during a thermal therapy course and the characteristic length of temperature variations are of the same order about 1 cm. So, particular details of spatial variations in the blood perfusion rate on scale much less than 1 cm are not the factor. On the other hand, the damaged part of the vascular network located inside the necrosis domain is most probable to be made up of an artery and vein having supplied previously this region with blood as a whole and of shorter vessels formed by their branching. Indeed, the mean volume of living tissue falling on a single small artery (or vein) of a fixed length $`l`$ is about $`l^3`$ . Thus, the regions where the total blood perfusion is directly controlled by different vessel of fixed length do not cross each other substantially and the architectonics of microcirculatory bed can be approximately represented as the binary tree embedded uniformly into the cellular tissue . Therefore, the region containing the vascular network part in which blood flow is remarkably disturbed because of the necrosis formation does not exceed substantially the necrosis domain. The latter enables us not to make difference between the given layer and the interface $`\mathrm{\Gamma }`$ and to deal with a sharp jump of the blood perfusion rate at the necrosis interface. The desired boundary condition imposed on the averaged blood perfusion rate $`j_v`$ must obey the law of blood conservation, which in this case can be written as

$$\underset{Q_\text{t}}{}𝑑𝐫j_v(𝐫,t)=\underset{Q_\text{t}}{}𝑑𝐫j(𝐫,t).$$
(38)
So in order to fulfill identity (37) we have to set the normal gradient of the averaged blood perfusion rate equal to zero at the interface $`\mathrm{\Gamma }`$:

$$_nj_v|_{\mathrm{\Gamma }^+}=0.$$
(39)
We note that the adopted assumption will not hold if a large vessel goes through the necrosis domain. However, the probability of this event is small enough and this case should be analyzed individually.

The random component $`\delta j`$ of the perfusion rate allowing for the vessel discreteness on scales $`\mathrm{}_v`$ is assumed to meet the conditions:

$`\delta j(𝐫,t)`$ $`=`$ $`0,`$ (40)
$`\delta j(𝐫,t)\delta j(𝐫^{},t)`$ $``$ $`j_v(𝐫,t)j_v(𝐫^{},t)g(𝐫,𝐫^{}).`$ (41)
Here the correlation function $`g(𝐫,𝐫^{})`$ is symmetric with respect to the exchange of the variables $`r`$, $`r^{}`$, i.e. $`g(𝐫,𝐫^{})=g(𝐫^{},𝐫)`$ and has to obey the general law of blood conservation:

$$\begin{array}{c}\delta j(𝐫,t)\underset{Q_\text{t}}{}𝑑𝐫^{}\delta j(𝐫^{},t)=\hfill \\ \hfill j_v(𝐫,t)\underset{Q_\text{t}}{}𝑑𝐫^{}j_v(𝐫^{},t)g(𝐫,𝐫^{})=0.\end{array}$$
(42)
For the living tissue without necrosis it is natural to assume:

$$g(𝐫,𝐫^{})=g^0\left(\frac{\left|𝐫𝐫^{}\right|}{\mathrm{}_v}\right),$$
(43)
where the function $`g^0(x)`$ is such that $`g(x)1`$ for $`x1`$, $`\left|g(x)\right|1`$ for $`x1`$, and

$$\underset{0}{\overset{\mathrm{}}{}}𝑑xx^2g^0(x)=0.$$
(44)
The latter reflects the blood conservation in three-dimensional space. A typical form of the function $`g^0(x)`$ is shown in Fig. 4

The system of equations (29), (30), (32), (33), (37) subject to the boundary conditions (34)–(36), (39) with expressions (40), (41) makes up the desired free boundary model.

Keeping in mind that the distributed model was widely used let us, now, compare the dynamics of the necrosis growth predicted with this model and with the stated free boundary model under different physical conditions. In this way, first, we will verify the basic question as to whether the value $`\mathrm{\Delta }`$ of the order of $`35^{}`$C is small enough so the ratio $`\mathrm{\Delta }/(T_{\text{max}}T_a)`$ may be regarded as a sufficiently small parameter for the free boundary model to hold. Second, we will actually justify the distributed model because, as will be seen, it does be independent of the particular details of heat transfer inside the layer of partially damaged tissue.

### 3.1 Model comparison

To be defined we simulate the necrosis growth in the tissue phantom shown in Fig. 5. The applicator presented by the hatched circle locally heats the tissue up to high temperatures of the order $`T_b100`$ C, which causes immediate tissue coagulation in the nearest neighborhood of the applicator. This can be the case, for example, because of the direct heat exchange between the tissue and the applicator heated to such temperatures or due to irradiation of laser light and its absorption in an adjacent thin layer (dotted region in Fig. 5). In the latter case it is typically used additional internal cooling of the applicator boundary, so the maximum $`T_b`$ of the temperature attained just near the applicator does not practically depend on heat diffusion into the surrounding tissue. So keeping in mind the possible vaporization control over the temperature maximum we will treat the value $`T_b`$ as a boundary temperature fixed at the interface of a certain radius $`r_0`$. Generalizing both these situations let us confine ourselves to the analysis of local coagulation assuming that:

* at the initial time $`t=0`$ the necrosis region under consideration is layer of zero thickness whose finite radius is $`r_0`$: $`\zeta (r)=0`$ and $`T(r)=T_a37`$ $`{}_{}{}^{}C`$ for $`r>r_0`$;
* the following necrosis growth is governed solely by heat diffusion, i.e. we set $`q(𝐫,t)=0`$ for $`r>r_0`$;
* at the boundary $`r=r_0`$ the temperature is a fixed value: $`T(r_0)=T_b100`$ $`{}_{}{}^{}C`$;
* at distant points the tissue temperature is equal to $`T_a`$: $`TT_a`$ as $`r\mathrm{}`$.

In other words, we confine ourselves to the region $`r>r_0`$ where the layer $`r_0<r<r_0+\mathrm{}(t)`$ of thickness $`\mathrm{}(t)`$ represents the necrosis domain whose growth is directly governed by the heat diffusion only. The processes in the region $`r<r_0`$ remain beyond consideration, i.e. we ignore the real dynamics of initial coagulation in the immediate vicinity of the applicator boundary. The typical duration of the latter process can be estimated as $`1/\omega (T_b)1`$ sec for the applicator directly heating the surrounding tissue and as $`\mathrm{max}[1/\omega (T_b),q(r_0)/(c\rho (T_bT_a))]`$ for the laser applicator. In the present analysis this duration is regarded as a small parameter.

Keeping in mind applicators of various forms we study the necrosis growth in the one-, two-, and three-dimensional tissue phantoms. The tissue response to temperature variations is described by equation (30) and expression (31) with $`\alpha =0.1`$ because the tissue response is sufficiently strong and the blood perfusion rate can locally increase by tenfold . For the main tissue parameters we use the values pointed out in Sec. 1.1.1 which give us the following estimates

$$\tau _0=\frac{1}{fj_0}6\text{min}\text{and}\mathrm{}_0=\sqrt{\frac{\kappa }{c\rho fj_0}}10\text{mm}$$

of the characteristic temporal and spatial scales, $`\tau _0`$ and $`\mathrm{}_0`$. We compare the dynamics of the necrosis growth, namely, the time dependence of the thickness $`\mathrm{}(t)`$ of the necrosis layer and the coagulation temperature $`T_{\text{cg}}(t)`$ predicted by the developed free boundary model with that given by the distributed model. In the latter case we use expression (23) and specified the required dependence of the effective thermal conductivity $`\kappa _{\text{eff}}`$ on the fraction $`\zeta `$ of undamaged tissue by the function

$$\kappa _{\text{eff}}(\zeta )/\kappa =(F1)\zeta +1.$$

Besides, for the distributed mode the necrosis interface $`\mathrm{\Gamma }`$ is specified by the condition $`\zeta |_{r\mathrm{\Gamma }}=0.5`$. For the sake of simplicity in the model comparison we also ignore the effect of vessel discreteness and identify the true and averaged blood perfusion rates.

First of all, Fig. 6 demonstrates the typical form of the temperature distribution $`T(r)`$, the distribution of the undamaged tissue fraction $`\zeta (r)`$, and the blood perfusion rate $`j(r)`$ obtained for the one-dimensional tissue phantom in the frames of the distributed model. As seen in Fig. 6 the region wherein the undamaged tissue fraction $`\zeta `$ varies substantially in space is sufficiently thin. So on spatial scales characterizing the temperature decrease such an increase of the value $`\zeta (r)`$ may be treated as a sharp jump. The latter actually justifies using the free boundary model for the given values of the tissue parameters, namely, the capability for assigning to an effective necrosis boundary a certain coagulation temperature $`T_{\text{cg}}`$ and certain values (on both the sides) of the temperature gradient.

Fig. 7 demonstrates the fact that the distributed model leads practically to the same dynamics of the necrosis growth as that predicted by the free boundary model for one-, two- , and three-dimensional tissue phantom (figures ($`a`$, $`b`$), ($`c`$, $`d`$), and ($`e`$, $`f`$), respectively). These results have been obtained for the tissue phantom with the strong ($`j_{\text{max}}=10j_0`$) and delayed ($`\tau _v=2`$ min) response to temperature variations and so actually comprises the characteristic features of the necrosis growth in the tissue without response as well as with a strong immediate response.

The present figures not only clearly justify that the particular details of thermal coagulation in a real partially damaged layer are not the factor but also show that the ratio $`\mathrm{\Delta }/(T_{\text{max}}T_a)`$ for $`\mathrm{\Delta }35^{}`$C can be treated in fact as a small parameter. So the two models may be regarded formally as equivalent, but the free boundary one does not contain self-inconsistent elements. In addition the free boundary model can be applied to constructing a faster numerical algorithm of simulating the necrosis growth because in this model we need not consider the thin layer of partially damaged tissue. In the free boundary model we deal with the temperature field only in the necrosis region and the region of undamaged tissue where it is sufficiently smooth.

In this analysis we have gotten another characteristic feature of local thermal coagulation illustrated in Fig. 8. The free boundary approximation is rigorously justified provided the temperature distribution inside the layer of partially damaged tissue can be regarded as quasistationary. The latter is the case when, in particular, the time variations $`\delta _tT_{\text{cg}}`$ of the coagulation temperature $`T_{\text{cg}}`$ are small enough during the necrosis growth, which in mathematical terms may be stated as the condition $`\delta _tT_{\text{cg}}0`$ as $`\mathrm{\Delta }0`$. Fig. 8 demonstrates that this condition may be fulfilled. Indeed, the smaller is the parameter $`\mathrm{\Delta }`$ of the tissue damaged rate $`\omega (T,\mathrm{\Delta })`$, to a greater extent are smothered the time variations of the coagulation temperature $`T_{\text{cg}}`$ except for a short initial period of the necrosis growth. Besides, the given feature of the local thermal coagulation justifies, at least at the qualitative level, the a certain simplification of the free boundary model  that considers the coagulation temperature $`T_{\text{cg}}`$ fixed during the necrosis growth.

### 3.2 Fixed coagulation temperature approximation

In this approximation expression (35) is replaced by the conditions on the temperature at the necrosis boundary $`\mathrm{\Gamma }`$. Namely, it is assumed that inside the undamaged tissue the temperature cannot exceed the coagulation temperature $`T_{\text{cg}}`$, i.e.

$$T(𝐫,t)<T_{\text{cg}}\text{for}𝐫Q_\text{t}$$
(45)
and at the interface $`\mathrm{\Gamma }`$ the boundary value $`T_i`$ is either equal to the coagulation temperature, $`T_i=T_{\text{cg}}`$, or rigorously less: $`T_i<T_{\text{cg}}`$. The former case takes place when growing near the interface $`\mathrm{\Gamma }`$ the temperature comes closely to the value $`T_{\text{cg}}`$ and the interface has to move in order to keep up the boundary temperature $`T_i`$ inside the interval $`T_iT_{\text{cg}}`$. In the second case the interface is fixed. Both these conditions can be formally described by the expression:

$$(T_iT_{\text{cg}})(\frac{T}{t}|_\mathrm{\Gamma }\frac{T_i}{t})=0,$$
(46)
where the boundary temperature $`T_i(𝐬,t)`$ is treated as a function of the interface coordinates $`𝐬`$ and the time $`t`$.

This approximation, except for simplifying the mathematical analysis of local coagulation, demonstrates the fact that the necrosis growth limited by heat diffusion weakly depends on the particular details of thermal coagulation. All the characteristics of thermal coagulation are aggregated into the coagulation temperature $`T_{\text{cg}}`$ and the coagulation rate $`\omega (T)`$ is only required to depend on temperature sufficiently strong.

This fact is illustrated in the following way. We study the necrosis growth in the one-dimensional tissue phantom within the framework of the distributed model modified so to take into account the effect of the blood perfusion nonuniformities near the necrosis domain. In other words, in equation (14) we replace the true perfusion rate by the averaged one, $`jj_v`$, related by equation (37) subject to the boundary condition (39). As before, the effect of vessel discreteness is ignored. We keep track of the points $`x_{0.2}`$, $`x_{0.5}`$, $`x_{0.8}`$ specified by the equalities $`\zeta =0.2`$,$`0.5`$,$`0.8`$, respectively. In this way the dynamics of the necrosis growth is characterized by the time dependence of the coordinates $`x_{0.2}`$, $`x_{0.5}`$ , $`x_{0.8}`$ (in mm) and the corresponding temperatures $`T_{0.2}`$, $`T_{0.5}`$, $`T_{0.8}`$ (in degrees Celsius). Different conditions representing various possible limiting cases are considered. Namely, the tissue phantom without response to temperature variations, with immediate and delayed response is analyzed. In the first case the blood perfusion rate remains unchanged, $`j=j_v=j_0.`$ For tissue with immediate response ($`\tau _v=0`$) the perfusion rate can attain large values directly at the beginning of the necrosis growth, whereas for tissue with delayed response ($`\tau _v>0`$) this increase will occur only after a lapse of time of order $`\tau _v`$.

Fig. 9 shows the time dependence of the quantities $`x_{0.2}`$, $`x_{0.5}`$, $`x_{0.8}`$, and $`T_{0.2}`$, $`T_{0.5}`$, $`T_{0.8}`$ for the tissue phantom without thermoregulation ($`j_{\mathrm{max}}/j_0=1`$), with immediate strong response to temperature variations ($`j_{\mathrm{max}}/j_0=5`$, $`\tau _v=0`$) and for the tissue where the delay in the temperature response has a pronounced effect on thermal coagulation ($`\tau _v=2`$ min). As seen in Fig. 9($`b1`$$`b3`$) the decrease of the temperature $`T_{0.5}`$ becomes slow sufficiently quickly (within approximately $`12`$ min) and the temperature $`T_{0.5}`$ remains inside the interval $`6070^{}`$C during all the time of order of $`510`$ min which is a typical duration $`t_{\text{tr}}`$ of thermal treatment based on thermal coagulation. The width of this interval is sixth as less as the typical value of the tissue overheating, $`T_bT_a60^{}`$C. So at first approximation, the temperature $`T_{0.5}`$ may be treated as a fixed constant $`T_{\text{cg}}`$. Besides, whence it follows that the coagulation temperature $`T_{\text{cg}}`$ being a phenomenological parameter of this approximation can be estimated from the expression

$$\frac{1}{\omega (T_{\text{cg}})}t_{\text{tr}}.$$
(47)
Indeed, under such conditions thermal coagulation proceeds practically at the temperature $`T_{\text{cg}}`$ and the value of $`1/\omega (T_{\text{cg}})`$ is approximately the time it takes for the live tissue located before in the necrosis region to be damaged. Because of a strong temperature dependence of $`\omega (T)`$ this estimate gives us the value of $`T_{\text{cg}}`$ to sufficient accuracy. In addition, to make the comparison of the two model more clear we have used in the simulation the value $`T_{\text{cg}}=60`$ C found from expression (47) (for $`t_{\text{tr}}5`$ min) rather than the value of $`T_{\text{cg}}`$ approximating the dependence $`T_{0.5}(t)`$ to the best degree.

Another characteristics of thermal coagulation is illustrated in Fig. 9($`a1`$$`a3`$). Comparing the time dependence $`x_{0.5}(t)`$ with the curves “FBM” (describing the motion $`\mathrm{}(t)`$ of the necrosis interface in the free boundary model with the fixed coagulation temperature) we see that there are two stages of the necrosis growth. The former corresponds to the time interval $`(0,t_{\text{cg}}),`$ where $`t_{\text{cg}}1/j_{\text{max}}4`$ min. At this stage the necrosis domain grows fast enough and in the free boundary model the interface $`\mathrm{\Gamma }`$ reaches its limit position $`\mathrm{}_{\text{lim}}`$. This saturation of the interface displacement is due to the temperature distribution becoming stationary. At the latter stage (from $`t_{\text{cg}}`$ to $`t_{\text{tr}}`$) the necrosis domain grows slowly and in the free boundary model it is fixed. In other words, the proposed model makes the difference of these stages more pronounced. If the treatment is continued, the real necrosis domain will grow further and after a lapse of $`2030`$ minutes the necrosis domain will deviate significantly in form from that predicted by the given model. However such a prolonged treatment is typically used to produce a hyperthermia effect (without visible injury) rather than to cause thermal coagulation directly .

It should be noted that the existence of the two stages does not obviously result from the dependence $`T_{0.5}(t)`$ because in fact the necrosis continues to grow slow at the second stage too. However, as follows from the present analysis, these stages differ from each other not only in the necrosis growth rate but also according to the behavior of the temperature distribution. The slow stage is characterized by the quasistationary temperature distribution. In other words, at this stage time scales on which the size $`\mathrm{}`$ of the necrosis domain increases substantially are much larger than time scales on which temperature distribution becomes steady-state provided the necrosis boundary is fixed. The given property is caused by the exponential dependence of the damaged tissue rate $`\omega (T)`$ on the temperature. At the fast stage these scales are of the same order.

Concluding this subsection we note that the fixed boundary temperature approximation predicts the dynamics of thermal coagulation to a sufficiently good accuracy. Indeed, the motion of the partially damaged tissue layer is practically represented in Fig. 9($`a1`$$`a3`$) by the region bounded by the curves $`x_{0.2}(t)`$ and $`x_{0.8}(t)`$. So the free boundary model can be regarded as adequate, if the difference $`\left|x_{0.5}(t)\mathrm{}(t)\right|`$ does not exceed the value $`\left|x_{0.8}(t)x_{0.2}(t)\right|`$ remarkably. As seen in Fig. 9($`a1`$$`a3`$) this is the case except for the beginning of the growth and times larger than 10 min, where we have gone beyond the fixed coagulation temperature approximation.

In the next section basing on the proposed model we will study the specific properties of the necrosis growth limited by heat diffusion.

## 4 Characteristics of the necrosis growth

The stated above model describes several different features of living tissue affecting simultaneously the necrosis growth limited by heat diffusion. So to penetrate deeper in its properties it is reasonable to analyze them individually, which is the subject of the present section.

### 4.1 Effect of the vessel discreteness

To clarify the essence of this effect let us confine our consideration to its qualitative analysis referring interested readers to specific papers  for a detailed analysis.

The vessel discreteness manifesting itself in the random component $`\delta j(𝐫,t)`$ of the blood perfusion rate causes random perturbation $`\delta T(𝐫,t)`$ in the temperature field. Actually $`\delta j`$ is the difference between the true perfusion rate and one averaged over scales about $`\mathrm{}_v`$. So the latter is also the correlation length of field $`\delta j(𝐫,t)`$ (Sec. 1.1.1). In this subsection we mainly consider the formation of the necrosis domain of size $`\mathrm{}\mathrm{}_T\mathrm{}_v\sqrt{L_n}`$, which enables us, at least formally, assume the inequality $`\mathrm{}\mathrm{}_v`$ to hold. The exception is the formation of a three-dimensional necrosis domain whose growth is due to a small applicator or the laser light absorption inside a small region  and it will be discussed below. Therefore in the case under consideration the mean amplitude $`\sigma `$ and the correlation length $`\lambda `$ of the temperature random nonuniformities $`\delta T(𝐫,t)`$ should be approximately the same as for such nonuniformities under uniform heating and be given by formula (9) and the estimate $`\lambda \mathrm{}_v`$. The field $`\delta T(𝐫,t)`$ in its turn affects the form of the necrosis interface $`\mathrm{\Gamma }`$ making it randomly perturbed as shown in Fig. 10. It is obviously, that the correlation length $`\mathrm{}_\mathrm{\Gamma }`$ of these interface perturbations is also about $`\mathrm{}_v`$, i.e.

$$\mathrm{}_\mathrm{\Gamma }\mathrm{}_v\frac{1}{\sqrt{L_n}}\mathrm{}_T\frac{1}{\sqrt{L_n}}\mathrm{}.$$
(48)
The remaining characteristic of the interface perturbations that should be found is their amplitude or, what is the same, the thickness $`\delta _\mathrm{\Gamma }`$ of the layer bounding these perturbations (Fig. 10).

In the fixed coagulation temperature approximation (Sec. 3.2) the form of the necrosis domain perturbed by the random temperature nonuniformities $`\delta T(𝐫,t)`$ due to the vessel discreteness is specified by the condition:

$$\left[T_v(𝐫,t)+\delta T(𝐫,t)\right]|_{𝐫\mathrm{\Gamma }}=T_{\text{cg}},$$
(49)
where $`T_v(𝐫,t)`$ is the regular component of the temperature field (i.e. the solution of equations (32), (33) for $`\delta j=0`$). In a neighborhood of the necrosis boundary $`\mathrm{\Gamma }`$ whose thickness is much less than the mean necrosis radius $`\mathrm{}`$ the spatial variations of the averaged temperature $`T_v(𝐫,t)`$ can be approximated by a linear dependence on the spatial coordinates $`𝐫`$:

$$T_v(𝐫,t)T_{\text{cg}}G\left(r\mathrm{}\right).$$
(50)
Here $`G`$ $`(T_{\text{max}}T_a)/\mathrm{}(T_{\text{cg}}T_a)/\mathrm{}_T`$ is the temperature gradient near the necrosis boundary $`\mathrm{\Gamma }`$ and $`r=\left|𝐫\right|`$ provided the origin, $`𝐫=0`$, is placed at the necrosis center. Then from (49) and (50) we obtain $`G\delta _\mathrm{\Gamma }\sigma `$ whence taking also into account expressions (5), (9) we get the following estimate chain:

$$\delta _\mathrm{\Gamma }\frac{1}{\sqrt{L_n}}\mathrm{}_v\frac{1}{L_n}\mathrm{}_T\frac{1}{L_n}\mathrm{}.$$
(51)
Expressions (48) and (51) substantiate, at least formally, the validity of the following scale hierarchy: $`\delta _\mathrm{\Gamma }<\mathrm{}_\mathrm{\Gamma }<\mathrm{}`$. So, when describing the growth of a necrosis domain as a whole one can ignore the effect of the vessel discreteness and identify the averaged tissue temperature $`T_v(𝐫,t)`$ and the true one $`T(𝐫,t)`$ . However when the attention is focused on the particular form of a necrosis domain, i.e. spatial scales smaller than the mean necrosis size are under consideration such necrosis perturbations are substantial.

To justify the latter statement and to measure the effect of these perturbations on the necrosis form we compare the value $`\delta _\mathrm{\Gamma }`$ with the thickness $`\delta _\zeta `$ of the layer of partially damaged tissue where thermal coagulation is under way. As results from (20) and (51) for typical values of the tissue parameters we get

$$\frac{\delta _\mathrm{\Gamma }}{\delta _\zeta }\frac{(T_{\text{cg}}T_a)}{L_n\mathrm{\Delta }}2.$$
(52)
Therefore, these random perturbations of the necrosis interface exceed sufficiently its physical thickness, i.e. the thickness of the layer of partially damaged tissue. In other words, because of the temperature nonuniformities due to the vessel discreteness the necrosis form should be disturbed substantially, deviating remarkably from that predicted by the distributed model dealing solely with the averaged tissue temperature $`T_v(𝐫,t)`$.

However when, for example, the characteristic size of the applicator is sufficiently small (about a few millimeters) the temperature field becomes substantially nonuniform not only due to heat dissipation caused by blood perfusion but also because of the geometric factor caused by the heat propagation from the applicator into the surrounding tissue. In three-dimensional space the temperature field $`T(𝐫)`$ induced by a small localized source even inside a tissue without perfusion will decrease as $`1/r`$. Under such conditions it is this spatial decrease in temperature that mainly controls the necrosis growth rather than heat diffusion into the surrounding live tissue. So the thickness $`\mathrm{}`$ of the necrosis domain is less then $`\mathrm{}_T`$, namely, can be about $`\mathrm{}_v`$ and the perturbations of its boundary due to the vessel discreteness are depressed considerably .

Expressions (48) and (51) are actually the main results of the present subsection. In the following subsections we will focus out attention on the growth of a necrosis domain as a whole and, thus, will ignore the vessel discreteness effect.

### 4.2 Effects of the tissue response to heating

In this section we analyze how the tissue response to local strong heating affects the necrosis growth. To single out these effects in their own right and to simplify the analysis we confine ourselves to studying the necrosis growth in the one-dimensional tissue phantom described in Sec. 3.1 (Fig. 5) and adopt the fixed coagulation temperature approximation.

First, we consider the tissue with immediate response, corresponding to $`\tau _v=0`$ in equation (30). In this case the thickness $`\mathrm{}_{\text{lim}}`$ of the necrosis domain attained during the fast stage of the coagulation is determined by the stationary temperature distribution. In particular, for the tissue phantom without thermal regulation ($`\alpha =1`$) the blood perfusion rate is constant, $`j=j_0`$, and as follows from equations (32) and (33):

$$\mathrm{}_{\text{lim}}^{\alpha =1}=\sqrt{\frac{F\kappa }{c\rho fj_0}}\frac{T_bT_{\text{cg}}}{F(T_{\text{cg}}T_a)}1\text{cm}.$$
(53)
For the tissue with thermal regulation expression (53), after the replacement $`j_0j_{\text{max}}`$, may be also used to estimate the value $`\mathrm{}_{\text{lim}}`$, thus

$$\mathrm{}_{\text{lim}}\sqrt{\frac{F\kappa }{c_t\rho _tfj_{\text{max}}}}.$$
(54)
The dynamics of coagulation under such conditions is represented in Fig. 11 for different values of the parameter $`\alpha =j_0/j_{\text{max}}`$.

Fig. 11$`a`$ shows the size $`\mathrm{}`$ of the necrosis domain as function of time. The higher is the tissue response, the smaller is the necrosis domain and the shorter is the fast stage of the necrosis growth. The duration of the fast stage is actually the characteristic time $`t_{\text{cg}}`$ during which the size of the necrosis domain attains values about $`\mathrm{}_{\text{lim}}`$. Comparing the numerical values of the corresponding quantities we find that the duration of the fast stage can be estimated by the expression

$$t_{\text{cg}}\frac{1}{fj_{\text{max}}},$$
(55)
which conforms to the general properties of heat transfer in living tissue. Indeed, in the given model, as follows from equations (32) and (33), the time it takes for the temperature distribution to become steady-state is about $`fj_{\text{max}}`$ and the establishment of this steady-state (in reality, quasistationary) temperature distribution is actually the essence of the fast stage.

Fig. 11$`b`$ shows the temperature distribution for the tissue without thermoregulation (curve 1) and for the tissue with strong immediate response (curve 2, $`j_{\mathrm{max}}/j_0=10`$). In order to compare them with each other, lengths are measured in units of the corresponding necrosis size. Whence it follows that, in contract to the time dependence $`\mathrm{}(t)`$, the tissue response to heating practically does not affect the form of temperature distribution.

When the tissue response to temperature variations is sufficiently strong ($`\alpha 1`$) the blood perfusion rate $`j`$ becomes substantially nonuniform. In this case the averaged perfusion rate $`j_v`$ differs significantly from the true one $`j`$, which has a certain effect on the growth of small necroses. The latter feature is illustrated in Fig. 12

Fig. 12$`a`$ compares the dynamics of the necrosis growth described by the given model and by the same model where, however, equation (37) is omitted and the replacement $`j_vj`$ is made. Fig. 12$`b`$ shows the distribution of the tissue temperature $`T`$, the true blood perfusion rate $`j`$, and the averaged one $`j_v`$ that occur when the tissue responds in such an intensive way. In this case, as seen in Fig. 12$`b`$, the averaged blood perfusion rate can be twice as less as the true one. The latter has a certain effect on the necrosis growth because ignoring the difference between $`j_v`$ and $`j`$ gives a more lower estimate of the necrosis size (Fig. 12$`a`$).

Now we consider how the delay in the tissue response can affect thermal coagulation. This effect is remarkable when the delay time $`\tau _v`$ is comparable with the duration of the fast stage $`t_{\text{cg}}`$. So we may confine ourselves to values $`\tau _v`$ $`t_{\text{cg}}`$. The difference in dynamics of the necrosis domain growth for the tissues responding immediately ($`\tau _v=0`$) and with a certain delay ($`\tau _v2`$ min) is illustrated in Fig. 13 and Fig. 14.

As seen in Fig. 13$`a`$ when the tissue response is delayed essentially the necrosis domain, firstly, grows fast enough keeping ahead of one growing in the tissue with the immediate response. Then the growth of the necrosis domain is suppressed and further its form remains unchanged. It should be noted that in this case the saturation of the necrosis domain growth is not due to the temperature distribution becoming stationary. Its reason is the increase of the blood perfusion rate after a lapse of a time about $`\tau `$. Under such conditions there is enough time for the size $`\mathrm{}`$ of the necrosis domain to attain values of order $`\mathrm{}_{\text{lim}}^{\alpha =1}\sqrt{1/j_0}`$ until the blood perfusion rate increases substantially. These values could not be attained if the tissue response were not delayed. So after the blood perfusion rate increases, the following growth of the necrosis domain becomes impossible and the temperature $`T_i`$ at the interface $`\mathrm{\Gamma }`$ must go below $`T_{\text{cg}}`$ and the necrosis domain has to cease to grow. This behavior of the interface temperature is illustrated in Fig. 14. It should be noted that such a saturation of the necrosis growth for real tissues seems to be more pronounced because a real necrosis continues to grow at the slow stage until the blood perfusion rate becomes high enough.

Fig. 13$`b`$ shows the distribution of the temperature and blood perfusion rate at different time moments for the tissue with delayed response. As before, the length is measured in the corresponding values of the necrosis size in order to compare these distributions. As time elapses the perfusion rate increases due to the tissue response. In contrast to this behavior of the perfusion rate, the form of the temperature distribution practically remains unchanged. The comparison of the given result and one obtained for the tissue responding to temperature variations immediately (Fig. 11$`b`$) leads us to the conclusion that the form of the temperature distribution occurring in tissue during the necrosis growth depends weakly on the specific values of the tissue parameters. This conclusion, in particular, forms the basis for applying variational techniques to analysis of the thermal coagulation due to laser-tissue interaction.

It should be noted that the conclusion concerning the universal form of the temperature distribution has been made analyzing the necrosis growth in the tissue phantoms that differ only in the properties of thermoregulation. The other tissue parameters (for example, the tissue thermal conductivity $`\kappa `$ and the initial value $`j_0`$ of the blood perfusion rate) take on particular values fixed in the present analysis. This raises the question of whether the stated conclusion will hold if we change these parameters too. However, choosing the appropriate units of time and length aggregating such parameters we can rewrite the governing equations in the dimensionless form. Thus their particular values are not the factor.

### 4.3 Effect of the applicator geometry

In the previous section we have actually studied the characteristics of the necrosis growth when the applicator is sufficiently large and the necrosis domain can be locally treated as a plane layer. If the applicator is small enough so its size is about or even less than the thickness of the necrosis layer the heat diffusion into the surrounding tissue gives rise to the temperature distribution depending on the applicator form. The latter in turn affects the necrosis formation. This effect is the subject of the present section reviewing the results obtained in .

As in the previous section we study the necrosis growth in the tissue phantom described in Sec. 3.1 (Fig. 5) but applying the free boundary model in its general form (Sec. 3), i.e. assuming the necrosis boundary motion to be governed by expression (35) because under such conditions time variations in the coagulation temperature become remarkable. However, keeping in mind that the latter fact is the main reason of the effects under consideration we will make no difference between the true and averaged perfusion rates in order to simplify the analysis.

Characteristic features of the necrosis growth depending actually on the applicator form are demonstrated in Fig. 15–Fig. 17. Namely, Fig. 15$`a`$ shows how the time dependence of the necrosis layer thickness $`\mathrm{}(t)`$ changes for the applicator of plane ($`d=1`$), cylindrical ($`d=2`$), and spherical ($`d=3`$) form for the tissue phantom with the same properties. Here is illustrated the dynamics of the necrosis growth in the tissue phantom without response to temperature variations for $`r_0=5`$ mm. Fig. 15$`a`$ shows that for one-dimensional tissue phantom the difference between the fast and the slow stages is not well distinctive and the necrosis growth is under way practically throughout the whole course of a typical thermotherapy procedure. For the three-dimensional tissue phantom we meet the opposite situation. In this case the necrosis growth exhibits the clearly evident tendency to saturation after a lapse of several minutes and after 12 minutes the thickness $`\mathrm{}|_{d=3}`$ of the necrosis layer can attain the value only somewhat less than the initial radius $`r_0`$ of the necrosis domain. At the same time in the one-dimensional tissue phantom the necrosis radius exceeds this value by several times. In the two-dimensional tissue phantom the necrosis growth is characterized by the intermediate behavior. It should be noted that in contrast to this the time dependence of the interface temperature $`T_i(t)`$ (Fig. 15$`b`$) is practically of the same form for the three cases.

These results prompt us that the rate of the necrosis growth in the three-dimensional tissue phantom should depend substantially on the initial necrosis radius $`r_0`$ at least for $`r_0\mathrm{}|_{d=1}`$, where $`\mathrm{}|_{d=1}`$ is the thickness that the necrosis layer would attain for an applicator of plane geometry. This fact is directly demonstrated in Fig. 16 showing the dynamics of the necrosis growth in the tissue phantom without response to temperature variations (figure $`a`$) as well as with an immediate strong response (figure $`b`$). At a fixed moment of time the necrosis layers can differ in thickness by a factor of three to four as the initial radius $`r_0`$ changes from 2 mm to 18 mm. As should be expected for the two-dimensional tissue phantom this dependence is smoothed (Fig. 16) and under the same conditions the necrosis thickness can increase not more than two-fold.

Fig. 17 illustrates the dependence of the necrosis growth on the initial necrosis radius in the most clear form. These curves have been obtained by fixing the time $`t=8`$ min and treating the point collection $`\{\mathrm{}|_{t=8\text{min}},r_0\}`$ for different values of $`r_0`$ as the partition points of certain continuous curves. In other words, this figure shows the thickness $`\mathrm{}(t,r_0)|_{t=8\text{min}}`$ of the necrosis layer as a function of its initial radius $`r_0`$ for the fixed moment of time under various physical conditions. We see in Fig. 17 that for the three-dimensional tissue phantom the curve $`\mathrm{}|_{d=3}(r_0)`$ practically goes into the origin of the plane $`\{\mathrm{},r_0\}`$ as formally $`r_00`$. When the initial radius $`r_0`$ becomes large enough ($`r_0>\mathrm{}|_{d=1}`$) the curve $`\mathrm{}|_{d=3}(r_0)`$ (as it must) tends to the value $`\mathrm{}|_{d=1}`$ (figure $`b`$). In other words, until $`r_0<\mathrm{}|_{d=1}`$ it is the initial radius of the necrosis domain that directly controls the size of the necrosis layer which can be attained during a typical course of thermotherapy. This property of the necrosis growth is due to the fact that in the three-dimensional space the temperature field governed by heat diffusion from a local source would remain substantially nonuniform ($`T(r)\frac{1}{r}`$ as $`r\mathrm{}`$) even though we ignored the heat sink effect caused by blood perfusion. So in the three-dimensional tissue phantom the temperature at the necrosis interface $`T_{\text{cg}}`$ inevitably has to decrease substantially as the necrosis layer grows. So due to the strong temperature dependence of the coagulation rate $`\omega (T_{\text{cg}})`$ the necrosis growth will be practically suppressed after a lapse of time it takes for the $`\mathrm{}|_{d=1}(t)`$ to exceed the initial necrosis radius $`r_0`$. Whence it follows that for real applicators of the spherical form whose size does not exceed several millimeters the thickness of the necrosis layer will be directly controlled by its size and, may be, by the depth of laser light penetration into the tissue because in this case heat diffusion standing along leads to the saturation in the necrosis growth. In the two-dimensional space (and more so in the one-dimensional space) heat diffusion tends to make the temperature field uniform. Thus, for the two-dimensional tissue phantom the dependence of the necrosis growth on the initial radius should be smoothed and the thickness of the necrosis layer will be finite even though we formally set $`r_0=0`$ (Fig. 17). In particular, for the case shown in Fig. 17 the change of the initial radius from zero to infinity causes the necrosis layer in the two-dimensional tissue phantom to increase in thickness by 2–3-fold.

The dynamics of the necrosis growth depends substantially on the physical conditions. However, it turns out that the time dependence of the coagulation temperature $`T_{\text{cg}}(t)`$ (i.e. the temperature at the necrosis interface) is practically insensitive to particular details of the growth conditions. This fact is demonstrated in Fig. 18. We have discovered that all the curves $`T_{\text{cg}}(t)`$ obtained numerically for various growth conditions are located inside the dotted region which is thin enough. This property of the necrosis growth is relative to the universality of the one-dimensional temperature distribution demonstrated in Sec. 4.2 and justify the conclusion that the general form of the temperature field can depend slowly on particular values of the main tissue parameters. Such a fact shows the feasibility of applying variational methods to analysis of the necrosis growth in more complicated cases and is undoubtedly worth special consideration.

### 4.4 Self-localization of the necrosis growth in tissue with a tumor due to thermoregulation

In section 4.2 we have demonstrated that the tissue response to heating affects substantially the necrosis growth limited by heat diffusion. It is due to the caused increase in the blood perfusion rate strengthening the heat dissipation in the surrounding undamaged tissue, which suppresses the further heat diffusion into it and, as the result, the necrosis growth. According to the experimental data  only the normal tissue can exhibit a significant increase in the blood perfusion rate in the response to the temperature growth. whereas inside a tumor an increase in the blood perfusion rate is depressed. This leads to some self-localization of the temperature inside tumor and self-localization of necrosis growth due to thermoregulation. The latter feature is illustrated in Fig. 19, 20 where results of mathematical modeling of the two and three dimensional cases are presented.

## 5 Laser induced heat diffusion limited tissue coagulation as a specific therapy mode

One of the main parameters controlling a thermotherapy treatment is the power of irradiation delivered into the tissue and the treatment duration. These parameters as well as a specific therapy mode should be chosen in such a way that gives rise to the necrosis of desired form. Different modes, in principle, can give rise to a necrosis region of the same form and, in this case, a specific mode may be chosen keeping in mind, for example, its optimality and stability. Various thermotherapy modes are singled out by those physical mechanisms that play the leading role and endow them with individual properties . The presented results enable us to regard the laser induced heat diffusion limited tissue coagulation as an individual therapy mode where heat diffusion into the live tissue, i.e. the tissue with active blood perfusion, plays the governing role.

To justify the latter statement let us briefly summarize the aforesaid. Namely, we have considered thermal tissue coagulation induced by local heating due to laser light absorption and limited by heat diffusion into the surrounding live tissue. The irradiation power is assumed to be enough high for the temperature inside the region affected directly by laser light to grow up to values about 100$`^{}`$C. Then this process is characterized by the following.

(I) In the considered case the necrosis growth is governed by the tissue thermal coagulation inside a layer of thickness about 1 mm that separates the necrosis domain and the live tissue and involves the partially damaged tissue. This layer is thin enough so the standard bioheat equation cannot be applied to describing temperature distribution in it because on such scales this equation does not hold. Therefore, the approach used previously by other authors (see, e.g.,  and references therein) which is based on the continuous description of heat transfer and thermal coagulation inside the layer of partially damaged tissue (the distributed model) is inconsistent.

(II) This inconsistency is overcome in the free boundary model of local thermal coagulation which regards the layer of partially damaged tissue as an infinitely thin interface of the necrosis domain whose motion is governed by the boundaries values of tissue temperature and its gradient. In addition, the free boundary model takes into account the perturbations of the tissue temperature caused by the vessel discreteness.

(III) The necrosis growth caused by heat diffusion limited thermal coagulation depends weakly on the particular details of heat transfer inside the layer of partially damaged tissue. This property justifies the distributed model. Nevertheless, we think that using the free boundary model is favored because in addition to being consistent it singles out the characteristic features governing the necrosis growth.

In particular, basing on the free boundary model we have found that:

(IV) Heat diffusion limited thermal coagulation involves two stages, fast and slow. It is the former stage during which the necrosis domain is mainly formed and its size $`\mathrm{}`$ can be estimated by expression (54). Expression (55) estimates the duration $`t_{\text{cg}}`$ of this stage. The latter stage is characterized by a substantially slower growth of the necrosis domain, with this slowdown becoming more pronounced as the applicator dimensions decrease. For the applicators of the spherical form whose radius does not exceed several millimeters the characteristic size of a necrosis domain formed during a typical course of thermotherapy is directly determined by the applicator radius and, may be, by the depth of laser light penetration into the tissue. Under such conditions even heat diffusion itself gives rise to the practical saturation in the necrosis growth and the fast and slow stages are most distinctive.

(V) The blood perfusion rate can become extremely nonuniform in the vicinity of the necrosis domain because of the strong tissue response to temperature variations. In this case in modelling thermal coagulation one should take into account that the temperature distribution is governed by the averaged blood perfusion rate rather than the true one. The delay in the tissue response on time scales of order $`t_{\text{cg}}`$ can affect thermal coagulation essentially. In this case, in particular, the duration of the fast stage is controlled by the delay time $`\tau _v`$. The size $`\mathrm{}`$ attained by the necrosis domain under such conditions is determined by the initial value of the blood perfusion rate $`j_0`$ rather than those ($`j_{\text{max}}j_0`$) occurring near the necrosis domain at the second stage due to the temperature tissue response.

(VI) The discreteness of the vessel arrangement causes perturbations in the tissue temperature, which is described by introducing additional nonuniformities of the blood perfusion rate regarded as random. These temperature nonuniformities it turn give rise to perturbations in the necrosis boundary whose amplitude $`\delta _\mathrm{\Gamma }`$ and the correlation length $`\mathrm{}_\mathrm{\Gamma }`$ are estimated as:

$$\delta _\mathrm{\Gamma }\frac{1}{L_n}\mathrm{},\mathrm{}_\mathrm{\Gamma }\frac{1}{\sqrt{L_n}}\mathrm{}.$$

Here the factor $`L_n\mathrm{ln}(l/a)4`$ ($`l/a40`$ is the characteristic ratio of the individual length to radius of blood vessels forming peripheral circulation systems). It should be noted that, first, these boundary perturbations are remarkable because they exceed in amplitude the thickness $`\delta _\zeta `$ of the layer of partially damaged tissue where thermal coagulation is under way. Second, the universal form of these relations is due to the vascular network being fractal in structure. However, for spherical applicators of small size (about several millimeters) the effect of the vessel discreteness is depressed because in this case blood perfusion does not affect the necrosis growth substantially.

Concerning certain additional aspects in the theoretical description of heat diffusion limited thermal coagulation we can state the following:

(VII) The free boundary model can be the basis of a faster numerical algorithm of simulating the necrosis growth because it deals with solely the necrosis region and the undamaged tissue where the temperature distribution $`T(𝐫,t)`$ is a smooth field.

(VIII) The time dependence of the temperature in the layer where thermal coagulation is under way (the coagulation temperature) is practically insensitive to the particular details of the growth conditions and to the particular values of the tissue parameters. For applicators of large dimensions (exceeding 1 cm) the form of the temperature distribution occurring in tissue during the necrosis growth also depends weakly on the particular values of the tissue parameters. The given properties form the basis for applying different variational techniques to analyzing heat diffusion limited thermal coagulation.

## 6 Acknowledgments

This work was supported by STCU grant #1675.