Methods of Finite Differences in the Domain of Time Applied to the Equation of Heat Biotransfer in Ocular Tissues

This article presents a methodology for two-dimensional modeling of the human eye to analyze the heat transfer in biological tissues, in this case, ocular tissues. The technique used was the finite differences method in the time domain. With a proposal for a numerical simulation of the problem, see that this type of work aims to promote the development of tools that help detect errors in biomedical engineering. Several medical procedures, including surgical ones that use laser, whose main limiting factor is the increase in temperature in living tissue, a fact that can cause irreversible damage, even evaporating. For example, refractive eye surgery using a laser. The simulation is intended to calculate the temperature gradient to guide the physician more precisely before an intervention. An example of this is certain types of inoperable tumors that can be irradiated with laser sources. Their function is to cause a rise in local temperature to destroy cancer cells without thermal damage to the healthy region.


Introduction
Over the years, several numerical and applied applications to engineering to analyze and solve problems in medicine and particularly involving heat transfer, which is the area of interest of this article.These tools are employed through numerical simulation methods of physical and biological phenomena.An example of this is laser surgeries, which have been increasing their applications because of their positive results, using the fundamental principles of heat transfer.Hypothermic treatments are used to destroy cancerous lesions, for example.Consequently, many therapies and medical devices operate by destructively heating or cooling diseased tissue while leaving adjacent healthy tissue undisturbed.With the growing advances in biomedical engineering, many medical devices have been developed that aid treatments, such as the destruction of eye tumors.These devices or tools must be developed to control the distribution of temperatures throughout the heat treatment and the distribution of chemical species in chemotherapy (Guimarães, 2003;Incropera, 1976).With the study of the Bioheating Transfer Equation -BHTE, it is possible to analyze and develop tools that

Description of the Mathematical Model
The heat biotransfer equation is a heat conduction equation, with a specific term for heat generation due to added blood perfusion.The analysis considers the effect of the heat pulse on living tissue -estimates of parameters of the Pennes equation (Guimarães, 2003).In his work, Harry H. Pennes (1948) proposed an equation that describes heat transfer in human tissues and organs, including the effects of blood flow through the heat biotransfer equation.The linear mathematical model obtained in this study represents the energy balance within biological tissues through the interaction of blood perfusion and metabolism.A term in the equation will need to be neglected because, according to (Sturesson & Andersson-Engels, 1995), generally, the generation of metabolic heat   is much smaller than the external heat deposited, in our case, the laser.Also, the heat source due to perfusion is given by   ; that is, it is like removing the heat due to the capillary vascularization present in the tissues.This term can be represented by equation 2 (Diller, 1992): Where:  volumetric blood perfusion rate;   blood specific mass;   specific heat of the blood;   temperature of arterial blood entering tissue;   temperature of the venous blood leaving the tissue.
The blood temperature is not the same at each point of the tissue.According to (Guimarães, 2003) the temperature of the blood that enters the capillary region is equal to the temperature of the arterial blood and the temperature of the blood that leaves is the temperature of the venous blood,   , which can be considered the same as the local temperature () of the tissue.The velocity of blood in the capillaries is minimal, with a Peclet number (which expresses the relationship between heat transfer by convection and heat transfer by conduction) much less than unity.This fact justifies the consideration that the temperature of the venous blood leaving the tissue is equal to the temperature of the latter.
Replacing the   by  at the local temperature.Soon   becomes: This blood perfusion rate cannot be measured directly (Charny, 1992), making it difficult to obtain values closer to real values for different types of tissues.Therefore, to arrive at acceptable values, as well as tissue properties in vivo, it is necessary to analyze before removal.Thus separating the effects of blood flow and metabolic heat generation so that the thermal properties of the medium can be deduced.
Assuming that heat transfer occurs by convection between the corneal surface and the external environment; internal driving; and convective heat transfer from the retina to the choroid due to blood flow in this region.Therefore, considering the specifications of these parameters that we will use from now on, the eye is represented by six specific areas with different thermophysical properties: the cornea, the aqueous, the lens or lens, the vitreous, and the retina, as shown in figure 1. (Silva, 2004).
Figure 1 -Cross-section scheme of the human eye.Source: Google Images Ocular biological tissue has a certain degree of sensitivity, as the most of them contain a large percentage of aqueous liquid.However, if they are exposed to high temperatures, thermal damage can occur.That is, it can cause first, second, and third-degree burns.The process of thermal tissue damage occurs with the denaturation of proteins, thus losing their biological activity.Therefore, irreversible destruction of living tissue can occur (Diller, 1992).For example, these effects may be due to significant heat sources, which in the proposed case is laser irradiation.
Several researchers have made the first quantitative assessments of damage or burning in living tissue.
Most of these experiments were done on animal tissue or human skin to determine the temperature-time exposure relationship to create various classifications of damage levels.With these data, it was possible to describe the thermal damage as a process rate dependent on temperature.It can be observed that the increase in damage occurs as an exponential proportion to the absolute temperature (Diller, 1992).
In this work, we will analyze the temperatures in order not to cause irreversible damage.Because, some analyzes, such as the one given by (Rol, 2000), indicate that it only starts to suffer thermal damage and protein denaturation from 37°C onwards.And Second (Diller, 1992), irreversible thermal tissue damage occurs at temperatures between 65°C and 75°C.The model by Henriques and Moritz (1947) was used because they were one of the first researchers to present a quantitative model to assess thermal damage in biological tissues, even if it is an in vitro experiment, with animal and human skin.They related thermal tissue damage as an integral part of a chemical process and used an Arrhenius-type equation to exemplify the damage mathematically.The finite difference method (FDM) is used to calculate the partial derivatives present in differential equations.The technique replaces each derivative or differential of the differential equations with numerical approximations of the same order of approximation using the Taylor Series.(CHAPRA,2013).
Taylor series expansion: (5.0) As the central derivative, obtained through the arithmetic mean of the other two derivatives, it is given by: The numerical methods replace the differential equation of heat conduction with a set of  algebraic equations for unknown temperatures at  points strategically chosen in the middle.The simultaneous solution of these equations results in the temperature values at these points discreet.
The first derivative of () at the midpoint is equivalent to the slope of the line tangent to the curve at this point, defined as: Thus, the ratio between the ∆ increment of the function and the ∆ increment of the independent variable, when ∆ → 0. If we do not take the indicated limit, we will have the following approximate relationship for the derivative: This expression is the finite difference form of the first derivative.The same equation can also be obtained by writing the Taylor series expansion of the function  over the point , And the second derivative of () at the midpoint is equivalent to the slope of the line tangent to the curve at that point, defined as: 2.1.2Finite Difference Method Two-dimensional problems involving geometries or boundary conditions are mostly complex; we can't develop solutions by analytical methods so quickly.Therefore, in this work, we chose to use the numerical methods, in our case, the finite difference method in the time domain.
The numerical solution allows the determination of temperature only discrete points, unlike the analytical solution that determines the temperature at any point of interest in the medium.The first step in numerical analysis is the selection of these points.For this, we should look at Figure 2, which shows how this can be done by subdividing the medium of interest into a series of small regions and specifying for each a reference point located at its center.The reference point is often called a nodal point (or simply a node), and the aggregate of points is called a nodal network (or mesh).The nodal points are identified by a numbering scheme which, for a twodimensional system, can take the form shown in Figure 2. The  and  positions are identified by the indices  and , respectively.The numerical accuracy of the calculations strongly depends on the number of nodal points used.If this number is large, the solutions will be more accurate.
The numerical determination of the temperature distribution requires that an appropriate conservation equation is written for each of the nodal points of unknown temperature.
The resulting set of equations must then solve the resulting equation set simultaneously for the temperatures at each node.
For any interior node in a two-dimensional system with no generation and uniform thermal conductivity, the exact form of the energy conservation requirement is given by the heat biotransfer equation, which is the heat equation with the energy generation term heat by blood perfusion.However, if the system is characterized in terms of a nodal network, it becomes necessary to work with an approximate, or finite difference, form of this equation.Thus, the finite difference equation suitable for the internal nodal points of a two-dimensional system can be deducted directly from the heat biotransfer equation: To obtain the form of finite differences in the time domain, we need to make some considerations.As we saw earlier, we can use the central difference approximations for the spatial derivatives represented by equation (16.0).
Using the subscript indices  and  can designate the positions of discrete nodal points in relation to the  and  axes.However, in addition to being discretized in space, the problem must also be discretized in time.The integer p is introduced for this purpose, where: And the finite difference approximation for the derivative concerning time in equation 16 is represented by: The superscript index p is used to indicate the temporal dependence of the temperature T, and the derivative concerning time is expressed in terms of the difference between the temperatures associated with the new (p + 1) and previous (p) instants time.
Thus, must perform the calculations at successive time instants separated by time interval ∆.Just as the solution of finite differences in space is limited to determining the temperature at discrete points, it also restricts temperature determination to discrete points in time.
Replacing equation 16 into  into equation 18, in the explicit form of the finite difference equation for the internal node  and  is: Convergence means that when x and t tend to zero, the results of the finite difference technique approach the true solution.
Stability means that errors at any stage of the calculation are not amplified but attenuated as the progress of the analysis.
Therefore, the ∆ must be less than a quarter of the ratio of the square of the mesh pitch and the thermal diffusivity of the fabric.

Implementation and Computational Modeling
With the physical-mathematical model developed and with the described equations that govern the domain of the object of study, including knowledge of the constitutive models and parameters of the parts that form it, the computational model to be used in the simulation must be built.And in turn, with the application of heat biotransfer equation in the simulation of biological systems, which will be subject to the most varied thermal procedures, requires great geometric flexibility in all stages of the adopted computational modeling (GUIMARÃES, 2003).
To start such an implementation, we need to discretize the domain.We will use image processing because a digital image is obtained from the sampling and quantization of a function.A function (, ) can be represented, only in this case,  and  are integers.Once the plane is sampled, we have various values representing the color.Each element of this array is called a Pixel.In these images, the pixel value contains a single number corresponding to an index for the table, including the color value in RGB (Red, Green, Blue), according to the table in Figure 3.The three components are reduced to 1. Table 1 shows the index mapping of the Color The pixels are normally square, generating a regular grid due to evenly spaced sampling, we can compare with the nodal network.
Using the Matrix Laboratory -Matlab software, the six main parts of the eye with different thermophysical properties were discretized: the cornea, the aqueous humor, the crystalline or crystalline, the vitreous humor, the sclera and the retina, as shown in Figure 5.After finishing the discretization of the image with the standard colors, a matrix of 350x350 pixels is generated with colors and their respective numbers shown in Table 2.The modeling procedure used in the study of bioheat transfer comprises the following steps: image discretization; geometric modeling and generation of discrete matrices; data pre-processing; numerical analysis program via Finite Differences Method; post-processing and visualization of results.The analysis used the Finite Differences Method implemented in Matlab language with a microcomputer with Intel i7 3 GHz processor with 4GB of RAM.
Thermal properties of various components in the human eye are listed in Table 3.The analytical solution of the Penne bioheat transfer equation is limited to some simple geometries with a high degree of symmetry.Still, using the Finite Element Method, we adopt the following considerations: I. Heat transfer within the eye occurs by conduction;

II.
Not considered a calorie generation   ≪ ; III.
Blood perfusion was not considered; IV.
The initial temperatures in the different regions of the eye were obtained through the study of (Amara, 1995) Blood temperature is considered to be 37°.The eye will be represented by six regions of different thermophysical properties: the cornea, the aqueous humor, the lens or crystalline, the vitreous humor, the sclera, and the retina.The boundary conditions and initial conditions are: I.
Convective heat transfer between the corneal surface and the external environment; II.
Convective heat transfer between the retina and the choroid.The initial temperatures of certain regions will be adopted second (Amara, 1995) with subdivisions in the vitreous, crystalline, aqueous.Also, we consider the sclera part of the vitreous humor as it has similar characteristics.The following table 4 presents the initial temperatures.The blood temperature is taken equal to 37°C.The ambient control temperature was considered constant and equal to 20°C.Source: Amara, 1995 The dimensions of the human eye were obtained from the work of Amara (1995) and Gokul (2013), as shown in Figure 6.We also consider that the sclera is part of the vitreous because it has similarities in its characteristics.Measurements were made from the central axis to its specific regions, as shown in Figure 5b, resulting in Table 5.In this way, we divided the vitreous into three parts, the crystalline into two and the aqueous into two, as shown in Figure 7.The reference model is a physical system described by equivalence differences consisting of the heat equation that includes the laser heat source and the threshold and initial condition equations.
The analytical system is transformed into an integral formulation where a Galerkin function is applied.Finite element method to drive a discrete algebraic system solved by numerical methods.The method used in this work is based on the discretization of an RGB image later transformed into a matrix of integers, transforming it into an algebraic system, thus inserting thermophysical properties and solving by the numerical method of finite differences.We will consider temperature control from 65° to 75°.
The transfer of heat in a biological tissue occurs more slowly because it has a large amount of water in its composition, causing heat propagation in the tissue to have a very long simulation time since their characteristics vary.After all, they have several different types of fabrics and with different temperatures.
We applied an external source at a temperature of 70° and used 1 second since most heat sources for surgery use the time interval from 30  to 1000 .The numerical approach is based on equation ( 20).Numerical simulation parameters are shown in table 6.
Then, spacings of 0.0029  were adopted in both x and y, which implies a time step of ∆ equal to 14.48 .By the stability, criterion ∆ must be less than a quarter of the ratio between the square of the mesh pitch and the thermal diffusivity of the fabric.Source: Author The heat propagation depends on the thermal diffusivity of the fabric because the higher the heat transfer, we can see from equation ( 14).
As for the ocular tissues, the thermal diffusivity varies between 1,11 × 10 −7 to 1,510 × 10 −7  2  ⁄ making the heat transfer occur more slowly, so we need large time steps for the simulation to be more accurate.
Table 7 shows the differences in thermal diffusivity between materials and fabrics.And Figure 8 presents the position of the heat source and the points for the output graph.The graph shown in Figure 9a shows the point in the corneal region, where the temperature only starts to increase after 1  (milliseconds), applying the external source with a temperature of 70°.The graph shown in Figure 9b shows the thermal damage done to the tissue in approximately 0.53 .The simulation presented values very close to the literature, demonstrating the method's functionality.The temperature is the same in the sclera until it reaches 10  because it is further away from the heat source.And because their thermophysical characteristics are different.Also, it did not suffer thermal damage; we can observe that the temperature does not suffer many variations.Aqueous humor is composed of large amounts of water and protein, and the graph shown in Figure 9b shows that conduction happens quickly but does not cause damage due to its composition.In the lens, because it has a smaller thickness than the other eye regions, the temperature variation is low compared to the other parts.Figure 9b showing the thermal damage graph, shows that there is not enough change to cause a burn in this region.
The vitreous humor consists of a transparent gel with hydration of approximately 98%, and it is two to four times more viscous than water.Also, due to its greater thickness and composition in liquid form, the temperature variation only occurs at about 15 .Not suffering any thermal damage because of the lowtemperature variation in the region.The retina is the outer part of the eye, so the heat transfer took place over a more extended period, even though the thickness was small.Most diseases of the human eye affect this specific part.The region in the simulation did not suffer any thermal damage because it was far from the source with a slight temperature variation.Tissue damage is defined as the denaturation of proteins or the loss of biological functions of molecules found in cells or in extracellular fluids (Diller, 1992).This thermal damage can occur when tissue temperature exceeds the range within which normal life processes exist.High and low temperatures can cause irreversible changes in molecules, resulting in thermal damage.Typical examples are burns and skin frostbite.The most successful mathematical models consider that thermal damage is a chemical reaction that depends on the temperature and the time interval to which the living tissue was subjected.The models describe the effect of heat on the chemical reaction rate and consider the variation in the concentration of undamaged molecules using Arrhenius' law.Although the meaning of the thermal damage function is not very well defined, it was considered that when the damage function assumes the value one (Ω  = 1), which corresponds to approximately 63% of the denatured proteins, complete necrosis, or damage, occurred.Irreversible damage to biological tissue (Henriques and Moritz, 1947).

Conclusion
Considering the observations described in this work, we can conclude that it is possible to determine thermal damage in ocular tissues using the time domain differences method with few thermophysical characteristics of the eye.The results presented were very close to the literature.It evidences the possibility of developing a tool that helps in the two-dimensional modeling of the human eye with bioheat transfer analysis.
Although, several works have pointed out challenges such as the lack of information about the natural values of the physical parameters of eye tissues due to the need for in vivo experiments.The simulation was able to show the behavior of temperature variation in the human eye.Given the results achieved for future research, treat the problem in three dimensions and implement an error estimator to control the stability of the proposed method.

Figure 3
illustrates the image discretization process and shows a possible color coding for each element of the matrixDiscretization of an Image Process.Source: Scuri,2002.

Figure 8 .
Figure 8. Position of heat source and points for output graph.Source: Author a) Graph of temperature of the parts of the eye (b) Graph of Thermal Damage Equation 1 represents the Pennes heat biotransfer equation: Metabolic heat is neglected as it is much smaller than the external heat source   ≫ .Then: And applying finite differences, the second derivative of ():

Table . Table 1 .
Index Mapping in Color Table

Table 3 .
Thermal properties of various parts of human eye

Table 4 .
Temperature without heat source