A Combination of Laplace Transform and Meshless Method for Analysing Thermal Behaviour of Skin Tissues

A novel numerical model is developed for evaluating the temperature distribution on a twodimensional skin tissue. The method is based on a combination of Laplace transform and meshless method as well as on the Pennes bio-heat equation. It is the first time to combine Laplace transform and RBF-MFS to study the performance of heat transfer in skin tissue governed by Pennes equation. First, Laplace transform of Pennes bioheat equation is applied to handle the derivative of time variable. Then, the particular solution is approximated by a linear combination of radial basis functions, and the homogeneous solution is approximately determined by the method of fundamental solution. The multi-subdomain RBF–MFS technique is implemented for analyzing problems containing different materials and/or multiconnected regions. The last step of the algorithm is the inverse Laplace transforms from Laplace domain to time domain. The efficiency and availability of the proposed method is assessed by several examples including healthy tissue, tissue with tumour and burned tissue.


Introduction
Biology and medicine have come into a new numerical region since human beings explored continuously and the mathematical and chemical systems are consummately set up in 21st century. The computers come up to accelerate the pace of exploration about the numerical medicine and reduce tremendously onerous calculation by hand as well. At present, the bioengineers try to find the most appropriate way to simulate the functions, feelings or consciousness of a human body. The mathematical models are built up to describe the process of functions of cells or tissues. Ideally, these mathematical models will be solved to obtain exact solutions. However, in practice, it is usually not possible to get exact solutions due to physiological complexity and variety of human's bodies. In some cases, only an approximate solution is required. The approximate solution can be evaluated instead of an exact solution if the error between them is acceptable or under the estimation for bioengineering application. The advantage for approximate solutions is that it can save a large number of time, expenditure and labour. To achieve that, meshless method is developed and utilized for scientists and bioengineers for analyzing heat transfer in skin tissues. In this report, meshless method is the main part of the algorithm to solve the bio-heat transfer problems.
The bio-heat transfer is the heat exchange between the blood vessels and the surrounding tissues. The blood perfusion throughout the vascular network system has a considerable impact on the local temperature distribution and the process of bio-heat transfer in living tissues. When the blood flows through the tissue, if the temperature between the blood and the tissue is somewhat different, it will induce convective heat exchange [1]. The effect of the bio-heat transfer is to moderate and balance the temperatures in vivo and in vitro to adapt human beings to the diverse environments.
Kreithand Goswami [1]indicated that there are several factors which will dramatically affect thermal interaction between blood and tissue, including the perfusion rate and the vascular anatomy. These factors may be considerably different among tissues, organs of the body and pathology. In addition, perfusion-based heat transfer interaction is very important to a large amount of physiological processes such as thermoregulation and inflammation.
Therefore, it is important to understand the phenomenon of heat transfer due to blood perfusion. It is also critical for bioengineers to know the process and mechanisms of heat transfer in the human skin. If the mathematical model for such bio-heat process can be built up, it will be extremely helpful and useful to predict the results of the surgical treatment; then, to assist doctors or nurses in practical clinical diagnosis by providing those results.
However, the process of heat transfer in living tissues or organs is very complex. Firstly, the perfusion rate of blood within different vessels, tissues and organs varies from time to time due to different activities, depending on factors such as physical activity, physiological incentive and environmental conditions [1]. Secondly, the vasculature is quite complicated because of the architectural and dimensional variety, which will cause the difficulty of modelling. Thirdly, the heat transfer will be intervened by the diseases such as tumour or therapeutics.
Skin covers almost the whole body of human beings. It protects bodies from the cold or hot environment and against the viruses or bacteria. It can be regarded as the first frontier which can be susceptible to damage. For this account, it is important to understand the physiological and physical properties of skin.

Objectives
This paper aims to develop a new meshless method for predicting the temperature distribution in the skin tissues by solving Penne equation due to the importance of bio-heat transfer in tissues. The new algorithmis derived from the actually mathematical model. The main techniques employed in this article include the Laplace transform and numerical inversion, Radial basis function (RBF) [2] and method of fundamental solutions (MFS) [3][4][5]. Firstly, the time dependence of the Pennes equation is removed by Laplace transform and then the system is replaced by a set of inhomogeneous modified Helmholtz equations. Then, RBF approximation and MFS are employed to construct the particular and the homogeneous solution of the modified Helmholtz equation, respectively. The muti-subdomain method is employed to extend this model to problems with two inhomogeneous domains, such as skin with tumour, which can induce different frequencies in the modified Helmholtz equation system. The algorithmis tested and proved by comparison with several examples which have already been validated before. By comparison with others, the advantages or disadvantages of the new algorithm will be evaluated for effectiveness and efficiency.

Assumption
The skin of human being is divided into three parts: epidermis, dermis and subcutaneous fat. They have different physical and chemical properties so that they can be seen as composite materials which are not homogenous or isotropic. However, if only one single region such as epidermis or dermis is chosen for numerical computation, the skin can be considered as a homogenous and isotropic material. For simplicity, we assume that the skin is only composed of one single region which has the same physical and chemical properties. Therefore, the model can be mathematically set up based on the isotropic region of the skin. The analysis assumes that no interfacial resistance exists between the heating source and the skin surface.

Fundamentals of Transient Heat Conduction
The thermal behavior of biological tissue is approximately governed by Pennes' bio-heat equation, including the effects of blood perfusion and metabolic heat generation [6], as follows: From the Pennes equation, it can be seen that on the right side the first expression describes the conduction of heat in the tissue which is induced by the temperature gradient. The second term stands for the heat flow transmission between the tissue and microcirculatory blood perfusion. The third term on the right hand side represents the internal heat generation due to metabolism and the last one is the spatial heating which is also called special heating. The internal tissues of the human beings have the ability to remove heat by both passive conduction and blood perfusion, which the Pennes equation depicts. The perfusion is defined as the non-vectorial volumetric blood flow per tissue volume in a region that contains sufficient capillaries that an average flow description is considered reasonable. Most tissues, including much of the skin and brain, are highly perfused, with a perfusion coefficient denoted byω b [7] Due to the single isotropic region of the skin tissue being analyzed, the thermal conductivity k is constant during the whole process. Moreover, the model is regarded as twodimensional for the following numerical computation. The Pennes' bio-heat equation can be then rewritten as: If the transient heat conduction problem occupies an arbitrarily shaped region Ω bounded by its boundary Γ shown in Fig.1 [13], = ( 1 , 2 ) ∈ Ω ⊂ 2 . Henceforth, the bold letter x indicates the two-dimensional Cartesian coordinate system. In general, Q m and Q r (x,t) are assumed to be the function of x 1 ,x 2 ,t. Besides the above bio-heat governing equation, the corresponding boundary conditions and initial condition should be provided to make the system solvable [8][9] 1) Dirichlet boundary condition related to unknown temperature field is 2) Neumann boundary condition for the boundary heatflux is 3) Convection or Robin boundary condition is or in a general form Where q represents the boundary heat flux defined as And n is the unit outward normal to the boundary Г of the interesting domain Ω. u Γ , q Γ , c Γ are the boundaries on which the temperature, the heat flux and the convection are known respectively. ( , ) u t x and ( , ) q t x present specified temperature and heat fluxon the related boundaries. The parameter h e is the convection coefficient and u e represents the environmental temperature. Moreover, in general form, β 1 , β 2 , and β 3 are known coefficients which could be expressed as below. .

Laplace Transform
Since the Penne's equation contains the derivative of the skin temperature with respect to time, we should apply some techniques here remove the derivative. Therefore the Laplace transform is utilized to change the problem from time domain to Laplace domain.
Assuming Q r (x,t) is independent of t in this paper to simplify the arithmetic, Laplace transform of Pennes equation is carried out as below Based on the differential theorem of Laplace transform:L(F'(t)) =sf(s)−F(0), Eq. (10) could be written as Where s is the Laplace parameter and 0 ( , 0) u u = x , which presents initial temperature distribution.
For simplicity, Eq. (11) can be written as Where λ and f(x) are known as following and (14) Note that the parameter s in (12) could be set at different values to see the difference of the temperature distributions.

RBF-MFS Scheme for Modified Helmholtz Equations
The solution of (12) can be written as a sum of the particular solution U p and the homogeneous solution U h , (15) where U p is the nonhomogeneous solution (16) and U h homogeneous solution  ( , ). Hence the variable s is not included in the expressions. The other functions and variables derived from U p (x) and U h (x) follow the same guideline in the following parts of the report. The particular solution U p (x) can be achieved by using RBF approximation [4] when it is approximated as (19) Where N I denotes the number of interpolation points in the domain of interest andφ i (x)=φ(r)=φ(│x-x i │)represents radial basis functions with reference point x i and unknown coefficients α i . In a similar way, the particular solution U p (x)could be obtained Where ψ i (x) are corresponding approximated particular solutions which could be acquired by solving the differential equations: The unknown coefficients α i can be uniquely determined by substituting (20) into (16) to obtain a set of equations.ψ is difficult to be obtained for Helmholtz-type operator though achievable by repeated integration for the Laplace operator [10]. Chen and Rashed [11]have given analytic formulae for ψ when φ was a TPS [11] 2 ln r r ϕ = There is another reverse way in obtaining approximate particular solutions [12] whereψ is first chosen directly andφ is derived from (21). e.g. the particular solution ψ is directly chosen as follows [12]: Though it seems to work well so far for many problems, it is difficult to prove mathematically how this approach is reliably correct [13].
Next, the homogeneous solution U h (x) is approximated in a standard collocation fashion by the using of MFS: Where β i are the coefficients to be calculated.U * j (x)=U * (│xx * j │) are the fundamental solutions of the modified Helmholtz operator (∇ 2 -λ 2 ).N S is the number of the source points on the virtual boundary where the source points{x * j } are located outside the solution domain.
A typical fundamental solution for a two-dimensional problem is whereK 0 is a modified second kind Bessel function with order zero.
A system of linear algebraic equations in matrix form can be obtained by collocation method: Once {β} is obtained, U h can be determined by using (26).
Note that the boundary conditions should be transformed from time domain into Laplace domain when solving (28).
In this work, the same strategy is employed to generate source points outside the domain as in [14,15] ( ) Where x j are boundary nodes,x c is the geometric center of the domain and δ is a dimensionless parameter. Then the complete solution U(x) for the modified Helmholtz equation can be written as Then U(x) could be converted to u(x) by inverse Laplace transform.
Noting the normal heat flux can be achieved after U(x) being calculated as ( ) To better illustrate the RBF-MFS algorithm, Fig.2 shows a simple example of the location and distribution of different boundaries and nodes.

Multi-subdomain Method for Multi-material Problem
When it is a multi-material problem, the problem becomes more complicated because λ in (12) is different for different material. So it is necessary to divide the solution domain into several subdomains with different material. For better illustration and self-contained purpose, a figure is taken from [13]. If there is a tumour situated underneath the skin (see Fig.3), the material properties for the tumour and for the skin tissue are different. So we have to divide the whole domain into two subdomains, which aretumour domain (simply connected domain in Fig.3) and tissue domain (multi-connected domain in Fig.3) and solve two different modified Helmholtz equations (12) in each domain respectively. The difficulties in solving this problem are that we only know boundary conditions on the whole domainand don't have any information about the inner boundary of the skin tissue. Moreover, the interface between tissue and tumour belongs to two subdomains at the same time which are solved independently, but continuity of temperature and heat flux on the interface still has to be satisfied. To conquer these difficulties, the multisubdomain method for multi-material problem in Ref. [13] is applied here. Following the method presented in [13], MFS formulation for the tissue domain Ω 1 , we have For the tumour domainΩ 2 , we also have where the superscripts 1 and 2 stand for variables associated with the subdomains Ω 1 and Ω 2 , respectively. On the interface boundary Г I , the following conditions must be satisfied: where subscript 'I' represents the interface boundary between the tissue domain and the tumour domain.
Rearranging the elements of the matrices and vectors in (32) and (33) and making use of (34) and (35) yields { } where the subscripts U and Q represent the boundary conditions respectively where temperature U and heat flux Q are known, B N is the number of collocation points on the boundary.
Collecting the terms with the same variable in (36) and (37) above and rearranging them, we have

Inverse Transform
In modelling many physical processes described by differential or partial differential equations, solutions can be obtained analytically using Laplace transforms. Many applied Maths/Physics subjects involve partial differential equations that may be solved by using numerical Laplace transforms. However the inversion of the expressions in the transform domain is sometimes a formidable task or even impossible and hence a numerical inversion technique may be employed to obtain an approximate solution. Fortunately, there are some efficient numerical methods available for computing the inverse transform. At present, there are several numerical inversion methods available, each with pros and consregarding accuracy and efficiency. An excellent survey and comparison of numerical inversion technique has been done by Davies and Martin [16]. Of these methods, one of the best was presented by HaraldStehfest [17] and has known as the Stehfest Algorithm. The algorithm itself is documented in the reference, but contains a typographical error that was later corrected. Stehfest points out that the algorithm requires an even number of coefficients, and is limited to functions that do not oscillate rapidly and that have no discontinuities. Fortunately, solutions to many practical equations rarely suffer from those problems [18].
All the results obtained in the computations above are in the Laplace domain, which should be transformed into the normal time domain by the inverse Laplace transform. As the results are calculated by RBF-MSF method, it is too complicated to perform the inverse Laplace transform analytically. Hence the numerical inversion is employed, named Stehfest Algorithm, as follows [17]: If f(s) is the Laplace Transform of F(t), an approximate value F a of the inverse F(t) for a specific time t=T is given by Note here f(s)=f(i × ln2/T). For a given time T and iteration number N, f(s) could be inversely transformed approximately to F a in the time domain. The algorithm will be validated in the next calculations.
After the numerical inversion, the temperature distribution could be solved in the skin modelled in 2D assumption. Then the distribution could be displayed in graphs to better illustrate the validity of this algorithm.

Numerical Examples
Four examples are included in this work which cover normal tissue, tissue with tumour and burned tissue. It is worth noting that these examples are only to testify the performance and applicability of the proposed numerical method.

Rectangular Domain of Skin Subjected to Temperature Impact
A benchmark example of is considered for validating the program. Referring to Fig.1, a rectangle of0.03m(x 1 )by 0.08 m(x 2 ) is chosen while the boundary conditions are: The initial condition is ( ) Note that all the boundary conditions above could not be applied into the computation directly. The Laplace transform should be performed to invert the boundary conditions from time domain to Laplace domain. As all the boundary conditions here are constant, it is easy to make Laplace transform following: Where a is constant. The steady-state analytical solution without metabolic heat ofthe rectangular domain is given by [19]: Here ω b =0.0005m 3 /s/m 3 and Q m = 420W/m 3 for healthy tissue. As stated above, the purpose of this example is to assess the efficiency of the proposed scheme. Comparisons are made between different values of N (used in the numerical Laplace inversion) and interpolation node (IN) for thin plate spline (TPS).The results are showed in Table 1.  Table 1, it is obvious that the number of interpolation nodes has almost no effect on the accuracy of this algorithm when the number of loop time N is 2, 4, 6, 8 or 10 in the numerical Laplace inversion. Although when IN = 120 and N=2, the absolute maximum error is the minimum, the outcome is more stable when IN = 221 and the differences between best results of the three IN are neglectable. Furthermore, multi-subdomain problems would be considered in the next part of this paper, then the interpolation nodes would be divided into different domains. The number of interpolation nodes in each subdomain should be around 100 at least. As shown in The analytical solution is a description of steady state temperature of the given area, so firstly, the time for the skin to enter in the steady state should be calculated. For a given number of interpolation nodes (IN=120) and loop time (N=2), the comparison is also made at different time (T). The results are shown in Table 2. From Table 2 we know that although the AME reaches its minimum at 10 4 seconds, the results become stable between 10 5 s and 10 6 s. Hence in the further examples, the time is set at 10 6 s to ensure the stability of the results.
The inverse method to calculate the particular solutions (ψ i and φ i ) is also tested and compared with standard method in Table 3. Here the number of interpolation nodes is 221 for both methods. Even if the inverse method works well in this problem, it might not suit other cases. Thus, for every different instance, the inverse method should be verified by the normal RBF method or analytical or theoretical outcomes. In order to evaluate the efficiency of the proposed method, the calculation is also carried out by using timestepping method in [13] on the same IntelCeleron 32 PC with 2.0GB memory. The CPU time for the proposed method is less than 2 minutes but around 6 minutes for the time-stepping method. We can see that the proposed method could save considerable time because there are only two loop times in the proposed algorithm. But the timestepping method would involve more iteration to achieve the final results. It should be mentioned that the more complicated the problem is, the more time can be saved.

Tumour Hyperthermia
It is know that the presence of a tumour would affect the distribution of both blood perfusion and metabolic heat generation [9].There is a difference of surface temperature distribution between the normal and tumoured skin. This difference canbe employed for non-invasive diagnostics of the physiological stage ofthe biological body. The calculation modelis taken from Liu and Xu [9].The tumour has a 20 by 20 mm size located 5 mm below the skin surface and equal distance(30 mm) to the insulating edges as in Fig.1. Refer to [9] and [13] for detailed illustration.
The boundary conditions are: .015m] is prearranged as the tumour domain for tissue with a tumour.
The temperature distribution in both healthy tissue and the tumour domain can be calculated from (38). The algorithm is carried out in Matlab codes. In the calculation, 56 field nodeson the real boundary and 221 uniform distributed interpolation nodes within the solution domain are used. Further, 56 field nodes are placed on the interface, 112 source nodes over the healthy tissue domain, and 56 source nodes over the tumour domain.See [13] for more details. Fig.4 shows the spatial temperature profiles for both healthy tissue and tumour tissue. It can be seen from Fig.4 that the temperature of tumour is higher than that of healthy tissue, this can be explained by the increased blood perfusion and metabolic heat generation caused by the tumour. Fig. 5 shows the difference of the epidermis temperature of two types of tissue. The trend matches well with the results presented by Liu and Xu [9] and Cao [13]. Finally, the inverse method is tested and compared to standard approach in both healthy and tumour tissues.The same temperature profiles as Fig.4 can be obtained.

Detection of Skin Burn Injury
The blood flow in skin tissue as well as heat conduction will alter to some extent due to skin burns. The degree of burns is responsible for the range of the change for the blood flow and heat conduction. In general, there are three degrees of burns related to the injury depth.
1) First -degree burn, also known as superficial burn, involves only the epidermis. Basal membrane (epidermisdermis interface) is preserved. An example of a first-degree burn is sunburn.
2) Second -degree burn or called partial thickness burn involves both the epidermis and the dermis. Second-degree burns are divided into superficial second-degree burn (superficial partial thickness burn) and deep second-degree burn (deep partial thickness burn), depending on the depth of damaged tissue. The first case is related to the epidermis and the basal membrane destruction. The second case concerns damages observed on superficial dermis (in addition to the basal membrane and the epidermis destruction).
3) Third-degree burn of full thickness burn concerns the destruction of epidermis layer and the entire (or nearly all) dermis layer." [20].
In order to evaluate the performance of the proposed approach on evaluating the harmful skin burn injury, in this example, we consider the skin tissue imposed under a low power of laser radiation. The function for the laser radiation The spatial temperature profiles for normal tissue and burned tissue are shown in Fig.7. The loss of water in the tissue and destruction of the vascular bed in burned tissue make larger temperature amplitude in it than in normal tissue.
Finally, the inverse method is tested and compared to standard approach in both normal and burnt skins. The same temperature profiles as Fig. 7 can be obtained.

Discussions
The algorithm in this paper is verified by some numerical examples and compared with the previous results obtained by Cao et al. [13] to confirm the efficiency and applicability. After comparison, the results match well.
In Section 4.1, the standard approach of TPS is tested by different parameters, N from 2 to 16, IN from 77 to 475. According to the results, the number of loop times (N) should be selected at 2 and the number of interpolation nodes (IN) should be set at no less than 100. The computational time is also dependent on N. It could save much time spent on the program running. The inverse method is also tested and compared to the standard approach. The accuracies of both algorithms are almost the same, only tiny differences, but the inverse method could save considerable time as the calculation of Bessel function is very time-consuming in the standard approach.
In Section 4.2, a special multi-domain MFS formulation was presented to handle multi-material and multi-connected domain problems using RBF-MFS formulation. It should be mentioned that, unlike the traditional domain decomposition method in which the boundary conditions are known on the whole multi-connected boundary, the boundary conditions are here known on the skin surface only, and are unknown on the inner boundary of the skin tissue. In hyperthermia treatment for tumours, two different cases are compared. One case is for healthy tissue and another case is for tissues with tumours. The results seem obviously distinguishable for healthy tissues and tumours in Fig.4. The difference in tumour area is due to a very high blood perfusion and metabolic heat generation. For healthy tissue, the general temperature will range from 37 ℃ to 37.15 ℃ while in tumour region, the temperature increases to about 37.35 ℃. When tackling with the tumour problem, cautions mush always be taken in every step of the calculations, especially the interface between tumour and healthy tissue. Different parts of the tissues have varied properties. Both healthy and tumour tissues have interpolation nodesat the interface. Although the nodes have the same coordinates, they have different parameters because they belong to different parts.
In Section. 4.3, the injury of skin burn is simulated to analyze the different temperature change on the surface of the burned skin. Unlike the hyperthermia instance, the boundary condition on the surface of the skin is not set as constant and spatial heating is also imposed on the whole domain. The spatial heating in this case utilizes the low power laser radiation. Actually, the specific absorption rate (SAR) of heat in the tissue is determined very difficult due to high nonlinearity. The SAR is a function of both location and time. For simplicity, only location dependent is considered in this case. Two occasions for normal and burned tissue are regarded as well. Compared with the normal tissue, the temperature in burned tissue is 1 ℃ higher due to a very low heat conduction and blood perfusion shown in Fig.7. The loss of water in the tissue and destruction to the vascular bed when burning account for the higher values of heat conduction and blood perfusion used for normal tissues.
In the validation examples, some simplifications are used. For all the cases, the tissues are simplified from a realworld 3D problem into a 2D problem which is not too complicated to establish the mathematical models and to apply the algorithm. For the case of tissue with a tumour, three boundaries out of four are assumed to be insulated with heat flux equal to zero. More accurate models could be created in the future study to simulate the problem, taking more practical details into account.
After all, this thesis proposed an efficient mesh-free method, which is easy to understand, has a simple solution procedure, high accuracy, and stability, and can analyze thermal behavior of inhomogeneous materials. It could also save considerable time when compared with the timestepping method.

Conclusions
From the examples discussed above, it can be revealed that the geometry and the properties of biological bodies induce the biological heat transfer significantly complex. Analysis on the heat transfer of biological bodies consists of determination of system parameters, such as thermal conductivity, blood perfusion, specific heat or metabolic generation. This can be done by using invasive or noninvasive experiments. The algorithm proposed in the thesis is capable of solving these problems mentioned above. In addition, it can deal with various problems occurring in the skin tissues and provide a generalized numerical approach to identify the temperature distribution and the heat flux as well. Laplace transform and inversion used in this paper is time independent and proved accurate. Moreover, speed of this algorithm could be explored in further study. Then this method could be applied to produce corresponding equipment for medical examinations.
This paper provides theoretical basis for noninvasive detection of skin abnormity. However, the difference is no so obvious so that professional equipment is required in practice to find the mere 0.1 ℃ difference when scanning the whole skin around the body. But technology would develop constantly and should not be the bottleneck. What's more, the examples included in this report only predict the temperature difference when the tumour is 5mm deep and only one case for skin burn. More study could be done to calculate the maximum depth recognizable with modern technology. More cases with different degrees of skin burns could be included in further research to compare the results. Ultimately, all the results obtained in the simulations should be verified by practical examples to see the efficiency of the method and the difference between the real examples and the models. After that, it could be started to introduce this method to medical practitioners and facility manufacturers.