On One Mathematical Model of Cooling Living Biological Tissue

When cooling living biological tissue (active, non-inert medium), cryomedicine uses cryo-instruments with various forms of cooling surface. Cryoinstruments are located on the surface of biological tissue or completely penetrate into it. With a decrease in the temperature of the cooling surface, an unsteady temperature field appears in the tissue, which in the general case depends on three spatial coordinates and time. To date, there are a large number of scientific publications that consider mathematical models of cryodestruction of biological tissue. However, in the overwhelming majority of them, the Pennes equation (or some of its modifications) is taken as the basis of the mathematical model, from which the linear nature of the dependence of heat sources of biological tissue on the desired temperature field is visible. This character of the dependence does not allow one to describe the actually observed spatial localization of heat. In addition, Pennes' model does not take into account the fact that the freezing of the intercellular fluid occurs much earlier than the freezing of the intracellular fluid and the heat corresponding to these two processes is released at different times. In the proposed work, a new mathematical model of cooling and freezing of living biological tissue are built with a flat rectangular applicator located on its surface. The model takes into account the above features and is a three-dimensional boundary-value problem of the Stefan type with nonlinear heat sources of a special type and has applications in cryosurgery. A method is proposed for the numerical study of the problem posed, based on the use of locally one-dimensional difference schemes without explicitly separating the boundary of the influence of cold and the boundaries of the phase transition. The method was previously successfully tested by the author in solving other two-dimensional problems arising in cryomedicine.

Abstract When cooling living biological tissue (active, non-inert medium), cryomedicine uses cryo-instruments with various forms of cooling surface. Cryoinstruments are located on the surface of biological tissue or completely penetrate into it. With a decrease in the temperature of the cooling surface, an unsteady temperature field appears in the tissue, which in the general case depends on three spatial coordinates and time. To date, there are a large number of scientific publications that consider mathematical models of cryodestruction of biological tissue. However, in the overwhelming majority of them, the Pennes equation (or some of its modifications) is taken as the basis of the mathematical model, from which the linear nature of the dependence of heat sources of biological tissue on the desired temperature field is visible. This character of the dependence does not allow one to describe the actually observed spatial localization of heat. In addition, Pennes' model does not take into account the fact that the freezing of the intercellular fluid occurs much earlier than the freezing of the intracellular fluid and the heat corresponding to these two processes is released at different times. In the proposed work, a new mathematical model of cooling and freezing of living biological tissue are built with a flat rectangular applicator located on its surface. The model takes into account the above features and is a three-dimensional boundary-value problem of the Stefan type with nonlinear heat sources of a special type and has applications in cryosurgery. A method is proposed for the numerical study of the problem posed, based on the use of locally one-dimensional difference schemes without explicitly separating the boundary of the influence of cold and the boundaries of the phase transition. The method was

Introduction
Prediction of the effects of cold on biological tissues leads to the study of Stefan's problems with sources that significantly depend on the desired temperature field. A distinctive feature of such problems is the existence of stationary solutions.
Physically, the time evolution of the temperature field u = u (P, t) is stabilized to the solution to the corresponding stationary problem u = u (P). This is achieved due to the fact that during cryogenic cooling of biological tissue, the removed heat flux is compensated by the sources of heat, blood and lymph flow, metabolism, and oxidative chemical reactions that arise in it.
It should be noted that the stationary limiting state does not always exist in problems with phase transitions. An example is Stefan's classical problem of freezing or thawing of a semi-bounded medium. From the exact solutions of such problems, it follows that the plane of phase separation propagates in the medium according to the law of the square root of t and, therefore, a stationary state is impossible.
As is known, in an inert medium with constant thermophysical characteristics, the rate of propagation of a thermal disturbance is instantaneous, which is mainly due to the linearity of the heat transfer equation. The question is different for media with temperature-dependent thermophysical characteristics, when the heat transfer equation becomes quasilinear. For example, in the case of a power-law dependence of the thermal conductivity coefficient on temperature, self-similar solutions of one-dimensional problems are obtained that describe heat waves propagating in an unperturbed medium with a finite velocity. However, here one can speak only about the spatial localization of the thermal disturbance only at a given moment of time, since at t → ∞ the disturbance region becomes unlimited. Neither linear nor quasilinear approaches to modeling allow describing the evolution of actually observed heat conduction processes stabilizing to limiting stationary states with a clear spatial localization of the temperature field. The latter include the processes of cryogenic cooling and freezing of biological tissues. A characteristic feature of such processes is the special dependence of heat sources on temperature. The models considered below take into account the effect of spatial localization of the thermal field actually observed in practice.
In the majority of literary sources (see, for example, the literature given in the review article [5]), the mathematical model is based on the Pennes equation [9], which, due to the linear dependence of heat sources on the desired temperature field, fundamentally does not allow describing this effect spatial localization of heat.
In this paper, we present a formulation and a method for the numerical study of a new three-dimensional three-phase boundary value problem of the Stefan type, corresponding to heat removal from the surface of biological tissue by a flat rectangular applicator. The desired temperature field, the boundaries of the influence of cold and phase separation are determined using the "through counting" technique [3]. The model takes into account the effect of spatial heat localization actually observed in practice and has applications in cryomedicine.
For the numerical solution of the problem, a locally one-dimensional difference scheme is constructed [2].

Formulation of the Problem
In cryomedicine, cooling and freezing of biological tissue is performed with instruments with various forms of cooling surface, which can either penetrate into the tissue (probes) or be located on its surface (applicators).
The mathematical model considered below corresponds to the process of cooling living biological tissue with a flat rectangular applicator.
In For numerical calculations, the approximate values of these constants can be found using the exact solution of a one-dimensional stationary problem [1].
When the temperature of the cooling surface of the cryoinstrument decreases, to about -3℃, a zone of frozen biological tissue appears (intercellular fluid freezes), and with a further decrease in temperature, to about -30℃, a zone of cryopreserved biological tissue appears (intracellular fluid freezes). To determine the boundaries of these zones, additional conditions are used on the corresponding isothermal surfaces: where 1 u is the biological tissue freezing temperature, 2 u is the cryogenic injury temperature. Surfaces 1 2 , z z -to be determined.
Equation (1) in this case has the form: In the case of equation (3), the coefficients Further, adding the initial and boundary conditions, the problem can be solved numerically using the locally one-dimensional method [2].
In the case when the freezing process is described by equation (2), the function of "heat content" will be written as follows: where η(x) is the Heaviside function; 1 2 , P P -known constants. In this case, equation (4) remains unchanged, and conditions (2) have been already taken into account by the equation. In this case, after "smoothing" the discontinuous functions , , , c H ρ λ [3], the problem can also be solved using locally one-dimensional difference schemes.

Difference Scheme
In the area under consideration, we introduce a grid (for simplicity, we consider a grid that is uniform in all directions): In this case, we will assume that the segments To find the value of the grid function ν at the =̄.
Depending on the mesh node in the boundary condition,  u ( w satisfies this condition in its properties, which cannot be said about it in the negative temperature interval. However, practice shows [7] that sweeping can also be used for 1 u u < .
As noted above, a characteristic feature of the heat conduction process during cooling of living biological tissue is stabilization to the limiting stationary state with a clear spatial localization of the temperature field. This feature is modeled by a specially selected nonlinear functional dependence of heat sources on temperature [1,6,7].

Sources of Heat
As in ordinary Stefan's tasks, the freezing of biological tissue is accompanied by the release of heat during crystallization, first of the extracellular and then of the intracellular water. This is modeled by the discontinuity of the specific internal thermal energy at points 1 2 , u u u u = = (see equation (3) Further, due to the fact that the biological tissue is permeated with a branched network of capillaries, supplying blood to the cooled non-frozen and frozen, but not cryopreserved areas of the biological tissue, internal sources of heat arise in it. According to the thermophysical meaning, the functional dependence of these heat sources on temperature should be limited, continuous and monotonically decreasing in the range of positive temperatures, and in the range of negative temperatures -limited and monotonically increasing. A condition is necessary for the existence of spatial localization ( ) w u ′ = −∞ . It has the following physical meaning: with an arbitrarily small perturbation of the initial temperature of biological tissue, arbitrarily small sources appear in it, the rate of growth of which is unlimited [4]. In addition, a necessary condition is the convergence of the integral: Here is the minimum temperature of biological tissue. When carrying out numerical calculations on a computer [6,7], the following functions of heat sources have proven themselves well: The parameters contained in the indicated dependences 0 and 0 < β < 1 should be determined experimentally.
The simulations were performed for the following values of thermophysical properties of biological tissue [6,8]: Comparison of the results shows a similar character of dynamics for freezing and cryodestruction isotherms for all three types of heat sources. Note that the second functional dependency contains only one experimentally determined parameter.

Conclusions
(1) The effect of spatial localization of heat, associated with the presence of nonlinear heat sources of a special type in biological tissue, makes it possible to formulate the problem in a limited area, which is important for numerical calculations. (2) The considered models take into account the fact that the freezing of the intercellular and intracellular fluid occurs at substantially different temperatures and at different times. (3) Nonlinear heat sources depend on numerical parameters, which must be determined experimentally, which in itself is an independent rather complex problem. (4) The used method of end-to-end counting with preliminary smoothing of discontinuous functions is suitable for the numerical study of models of the considered type of any dimension. (5) The connection between the steps in time and space, found when solving two-dimensional problems [6,7,10], indicates that the constructed difference schemes are conditionally stable. In addition, numerical experiments have shown the possibility of choosing a "large" time step (about two orders of magnitude larger than the space step), which is especially important when solving multidimensional problems. (6) The iterative process of Newton's method converges quickly enough: 5-7 iterations for the choice of