A Two-Dimensional Mathematical Model to Analyze Thermal Variations in Skin and Subcutaneous Tissue Region of Human Limb during Surgical Wound Healing

During wound healing, the metabolic activity associated with each phase must occur in the proper sequence, at a specific time, and continue for a specific duration at an optimal intensity. Any disturbance in appropriate thermal environment may complicate the wound healing process and may give rise to wound infection. In the presented paper a transient state two-dimensional mathematical model has been developed to analyse thermal variations in skin and subcutaneous tissue (SST) region of human limb. Due to circular shape of human limb, model has been developed in polar coordinates. The domain of the study consists of two types of tissues: abnormal tissues and normal tissues. The post surgery peripheral tissue of human limb during healing time is considered as abnormal tissues. The effect of variable density of blood vessels in dermal layer of both tissues on the physical and physiological parameters is incorporated in the model. The effect of healing on physiological parameters of abnormal tissue is incorporated by considering the physiological parameters to be function of time “t”. The effect of different climatic conditions is considered in the model. Taking into account the variable core temperature due to anatomy of arteries and variable physiological parameters in dermal layer of peripheral region, the well known Pennes’ bio heat equation is used to analyse the time-dependent temperature distribution of both normal and abnormal tissues. Comparison between temperature profiles of both normal and abnormal tissue has been done using finite element approach with bilinear shape functions in polar coordinates. A computer program in MATLAB has been developed to simulate the results.


Introduction
Enzymes play a key role in the healthy functioning of organs and glands of our body.They work properly in 20˚C -37˚C body temperature [1] [2].Any disturbance in body temperature disturbs the proper functioning of enzymes (metabolic activity).Thus it becomes important to study thermal variations in order to develop a mechanism for establishing the limits of thermal regulation in human body during the healing process.
Pennes [3] in 1948 first gave the well known bio-heat equation of heat transfer in perfused tissue for the analysis of temperatures in the resting human forearm.Perl [4] in 1962 used it to illustrate heat distribution and tissue blood flow in body tissues.Since then it has been used by many researchers and scientists [1] [5]- [11] for analysis in their respective areas of research.A lot of research on wound healing has also been done by considering different parameters like size and shape of the wound, the effect of growth factors, the effect of angiogenesis, the effect of oxygen concentration at the wound site etc. [12]- [18].A few scientists have studied temperature variations caused due to wound healing process after surgery [1] [19] [20].
In this paper, a two-dimensional finite element model has been developed for transient state temperature distribution in normal and operated peripheral tissues of human limb during healing time by using Pennes' bio heat equation [3].Assuming the cross section of human limb to be circular in shape, model has been developed in polar coordinates (r, θ).The effect of the density of blood vessels in different layers of dermal region on temperature distribution is incorporated by considering the physical and physiological parameters to be function of space coordinate "r" only in case of normal tissues and increasing function of time "t" also in case of abnormal tissues [1] [19] [20].The position of major arteries, lying deep inside the skin, carrying blood from the core to the limbs is not symmetrically distributed along the angular direction [21].Variable core temperature along angular direction "θ" [22] has been considered to incorporate the effect of position of major arteries in limbs.The outer surface is exposed to the environment and an appropriate boundary condition for heat loss at the outer surface has been incorporated.( )

Model Description
Here K, S, m b , c b , ρ, c, T and T b are thermal conductivity, metabolic heat generation rate, rate of blood mass flow per unit volume, specific heat of blood, tissue density, specific heat of tissue, unknown tissue temperature and temperature of blood.
The inner core of the human limb is at a variable temperature and hence it is formulated as: Here T α is greater than T β .The boundary condition at outer surface of skin is: ( ) Here h, T a , L, E and T n ∂ ∂ are heat transfer coefficient, atmospheric temperature, latent heat of vaporization, rate of evaporation and partial derivative of unknown tissue temperature along the normal direction to the skin surface.Heat transfer along angular direction is assumed to be negligible.The surface of the limb is insulated initially [1] [9].So initial condition is given by: ( ) ( )

Use of Finite Element Method-Discretisation of Domain
The circular cross section of annular shaped region of peripheral tissues of human limb is discretised into a total of twenty four elements (Figure 1).Skin tissues are naturally divided into three layers viz.epidermis, dermis and subcutaneous tissues.We have further divided dermis into four layers.Each layer is further divided into four coaxial circular sectors at 0, π 2, π and 3π 2 θ = . Thus there are in all twenty four elements and twenty eight nodes.The bio heat equation can be written for each element in discretised variational form as: where Ω (e) and Γ (e) represents the domain and outer boundary along the angular direction of e th element.In the second integral λ (e) is 1 for surface elements and zero for other elements.

Physiological Parameters for Normal and Abnormal Tissues
In this model, it is assumed that there is negligible variation in thermal conductivity, blood mass flow rate and metabolic heat generation rate along angular direction.

For Normal Tissues
There is no blood flow and metabolic heat generation in the epidermis.Hence it is assumed that M 3 = 0 = S 3 in epidermis.Physical and physiological parameters like thermal conductivity (K 1 ), rate of heat transfer due to blood mass flow (M 1 ) and metabolic heat generation rate (S 1 ) have almost constant value in subcutaneous layer.In dermis K 2 , M 2 and S 2 are function of radial coordinate only and expressed by interpolating linearly between the corresponding values at epidermis and subcutaneous tissues.

For Abnormal Tissues
Since surgery results in loss of blood and tissue, it is assumed that initially at time t = 0, blood mass flow rate and metabolic heat generation rate is almost negligible and gradually increase with time to attain normal values.
For abnormal tissues these are increasing function of time too.

Interpolation Function for Circular Coaxial Sector Element
We have chosen the bilinear shape function in polar coordinates for the variation of temperature within each element, , Here Here i, j, k and l are the local node numbers.The interpolation functions, thus, obtained are- , In this paper, the domain has been divided into circular coaxial sector elements of varying thickness i.e. "a" takes three different values: ( ) for subcutaneous tissues, dermis and epidermis respectively.Using the above values and expressions for physical and physiological parameters and expression for temperature variation (12) in (7), integral ( ) e I is evaluated as: Equation ( 16) has been differentiated with respect to each nodal temperature of the element then expressed in matrix form and minimized.The corresponding global matrix equation has been formed and assembled.Crank Nicolson Technique has been employed to solve the system of equations to obtain the temperature profile at each node for progressive time.A computer program has been developed in MATLAB to simulate the results.

Numerical Results and Discussions
In order to enhance the applicability of the present model in different climatic conditions, the numerical calculations have been done for three cases of atmospheric temperatures.The values of physical and physiological parameters taken [1] [10] [19] [22] [23] to obtain the numerical results are given in Table 1 and Table 2.   12-14) have been plotted for three different atmospheric temperatures for both, normal and abnormal tissues.
For plotting Figures 2-8 four nodes chosen are: one at core of the skin, one at the interface of subcutaneous and dermis, one at the interface of dermis and epidermis and one at the surface of the skin for each value of θ = 0(π/2)3π/2.For plotting Figures 9-11 seven nodes chosen are at r i , i = 0(1)6 for each value of θ = 0(π/2)3π/2.For plotting Figures 12-14 the five nodes chosen are at θ = 0, π/2, π, 3π/2 and 2π for each r i ; i = 0, 1, 5 and 6.These graphs give us an idea about the magnitude of the thermal effect of surgical wound on temperature profiles in the skin layers.
All the graphs (Figures 2-14) show that at the skin core nodal temperatures are given by the function f(θ) as assumed.
From Figures 2-8 it is observed that with the increase in time, there is steep decrease in temperature for first 10 minutes and 20 minutes for normal and abnormal tissues respectively.This is due to the fact that on exposure to atmosphere (at lower temperature), heat is lost from the tissues due to radiation and conduction.Heat is also lost due to evaporation of sweat, which is produced on activation of sweat glands when peripheral tissues are at higher temperature.Heat of peripheral tissues is lost in the form of latent heat of vaporization.This fall in temperature is higher in abnormal tissues than that in the normal ones.It is because of temporary vasoconstriction in incised blood vessels of the tissues [2] [24].Results show that, abnormal tissue temperature takes more time to obtain steady state than that of normal tissue.For normal tissues, the temperature becomes almost steady after  20 minutes.In case of abnormal tissues, after attaining the minimum value the temperature begins to rise for about 180 -250 minutes without and with different rates of evaporation (Figures 2-8).It is further observed that during this increment period, the rate of increase in temperature is fast for first 40 -50 minutes approximately followed by further increment but with slower rate in the remaining period.This may be due to the fact that initially the rate of blood mass flow and heat generated due to metabolic activity of the operated tissues are very small in comparison to normal tissues but as time progresses, the value of these parameters also increase which helps the abnormal tissue to return to normal temperature and restart cell mitotic division to accomplish wound healing [25] [26].Thus abnormal tissues take more time to attain steady state temperature in comparison to normal tissues and it can be validated on the basis of experimental results of Gannon [26].In all the graphs, the fall in tissue temperature is more at the skin surface in comparison to the interior tissues because more heat is lost at the surface due to convection, radiation and evaporation.For the same rates of evaporation, the decline in tissue temperature is more at lower atmospheric temperature than that at higher atmospheric temperature and the steady state temperature of the tissue at lower atmospheric temperature is less than that with the higher atmospheric temperature (Figure 2, Figure 3 or Figure 4, Figure 6 or Figure 5, Figure 7).This is because the temperature gradient at the skin surface is more at lower atmospheric temperature than that at higher atmospheric temperature and amount of heat transfer is directly proportional to temperature gradient.But with the slight increase in evaporation rate at higher atmospheric temperature, the temperature profiles become closer in spite of considerable amount of difference in atmospheric temperature.For example, it is observed that the temperature profiles for T a = 15˚C and E = 0.0 gm/cm 2 (Figure 2) and T a = 33˚C, E = 0.24 × 10 −3 gm/cm 2 (Figure 6) are close in spite of a difference of 18˚C in atmospheric temperature and a small rate of sweat evaporation at higher temperature.This highlights the significant role of evaporation of sweat by utilizing latent heat of vaporization in gaining thermal balance.For the same atmospheric temperature, the fall in tissue temperature increases as the rate of evaporation increases (Figures 3-5 or Figures 6-8).This again shows that rate of evaporation has a significant role in temperature regulation in SST region.
From Figures 9-11 it is observed that the slope of the curve changes at the junction of each layer (r = 2.5(0.1)2.9 cm).This is due to the different physiological properties of each layer.The steepness of the curves increases as we move from the inner core towards outer surface.It is because of the maximum temperature gradient at the surface which results in maximum temperature fall at the surface and this effect reduces as we move towards core of the limb.Also it is observed that at a given atmospheric temperature and evaporation rate, the thermal variation along angular direction decreases as we move radialy outwards from the core.It implies that effect of core temperature on thermal distribution of tissues decreases in angular direction with increase in distance from the core.From Figure 12, Figures 13(a)-(c) and Figures 14(a)-(c) it is noted that the variation in tissues temperature in angular direction, in each layer reflects the effect of assumed boundary condition at the core of limb.
On comparing Figure 9 with

Conclusions
The two-dimensional unsteady state finite element model has been developed to analyse thermal variations during wound healing in the peripheral tissues of human limb with bilinear shape function in polar coordinates.It is concluded that the use of coaxial circular sector elements is appropriate to incorporate the thermal variations along the circumference of the human limb.The variation in tissues temperature in angular direction, in each layer due to asymmetric arrangement of arteries is modelled well by assumed boundary condition at the core of the limb.The fall in tissue temperature is more at the skin surface in comparison to the interior tissues.The fall in temperature profile in abnormal tissues is more in comparison to that in normal tissues.After attaining minimum temperature values, abnormal tissues take more time (180 -250 min) to resume steady state temperature profiles in comparison to normal tissues.The time difference to attain steady state temperature profiles and magnitude of difference in temperature profiles between normal and abnormal tissues varies with ambient temperature and evaporation rate.The effect of core temperature on thermal distribution of tissues decreases in angular direction with increase in distance from the core.The present model can be extended to the deep tissues of human limbs where the properties also vary along angular direction.
The findings of the present work are in coordination with the experimental results of Gannon [26].This information can be used to generate the thermal information which may be useful to biomedical scientists in diagnosis and development of treatment regimen for surgical wounds.

Figure 1 .
Figure 1.Cross section of annular shaped region of peripheral tissues of human limb.In all there are twenty four coaxial circular sector elements and twenty eight nodes.

Figure 9 .
Figure 9. Nodal temperature versus distance from core graph for T a = 15˚C, E = 0.0 gm/cm 2 •min, η = ν = 0.01 for normal (black solid line) and abnormal (blue solid line) tissues at thermal equilibrium.Note: Nodal temperature of normal and abnormal tissues is almost same at equilibrium so overlapping of curves has occurred.

Figure 10 (
a), Figure 10(b) with Figure 11(a), Figure 10(c) with Figure 11(b), it is observed again that (for fixed evaporation rate) the fall in nodal temperature along radial direction increases as atmospheric temperatures decreases but lowering of atmospheric temperature has no significant effect on the fall of nodal temperature along angular direction (compare Figure 12 with Figure 13(a), Figure 13(b) with Figure 14(a) and Figure 13(c) with Figure 14(b)).It is because lowering of atmospheric temperature is symmetric about "θ" and increases temperature gradient significantly in radial direction only.On comparing Figures 10(a)-(c) or Figures 11(a)-(c), it is observed that at given atmospheric temperature, the increase in evaporation rate increases the decline in nodal temperature along radial direction but has no ef-fect on the thermal variations along angular direction (compare Figures 13(a)-(c) or Figures14(a)-(c)).It is because the increment in evaporation rate is also symmetric about "θ" and increases temperature gradient significantly in radial direction only.The steady state nodal temperature (SSNT) values shown in Figures2(a)-(c) can also be observed in Figure9for θ = 0 (or 2π), θ = π/2 (or 3π/2) and θ = π respectively.Further same values can also be observed in Figure12at θ = 0 (or 2π), π/2 (or 3π/2) and π respectively.Similarly, SSNT values shown in Figures3(a)-(c), Figures 4(a)-(c) and Figures5(a)-(c) can also be observed in Figures10(a)-(c) respectively and in Figures13(a)-(c) respectively at θ = 2π.The SSNT values shown in Figures 6(a)-(c), Figures 7(a)-(c) and Figures 8(a)-(c) can also be observed in Figures 14(a)-(c) respectively.

Figure 10 .
Figure 10.Nodal temperature versus distance from core graph for T a = 23˚C, (a) E = 0.0 gm/cm 2 •min, (b) E = 0.00024 gm/cm 2 min, (c) E = 0.00048 gm/cm 2 •min and η = ν = 0.01 for normal (black solid line) and abnormal (blue solid line) tissues at thermal equilibrium.Note: Nodal temperature of normal and abnormal tissues is almost same at equilibrium so overlapping of curves has occurred.

Figure 11 .
Figure 11.Nodal temperature versus distance from core graph for T a = 33˚C, (a) E = 0.00024 gm/cm 2 •min, (b) E = 0.00048 gm/cm 2 •min, (c) E = 0.00072 gm/cm 2 •min and η = ν = 0.01 for normal (black solid line) and abnormal (blue solid line) tissues at thermal equilibrium.Note: Nodal temperature of normal and abnormal tissues is almost same at equilibrium so overlapping of curves has occurred.

Figure 12 .
Figure 12.Graph between T(r i ), i = 0, 1, 5, 6 and θ at T a = 15˚C, E = 0.0 gm/cm 2 •min and η = ν = 0.01 for normal (black solid line) and abnormal (blue solid line) tissues at thermal equilibrium.Note: Nodal temperature of normal and abnormal tissues is almost same at equilibrium so overlapping of curves has occurred.

Table 1 .
Value of physical and physiological parameters.

Table 2 .
Value of M 1 and S 1 at different T a and E.