Simulation study of melanoma detection in human skin tissues by laser-generated surface acoustic waves

Abstract. Air pollution has been correlated to an increasing number of cases of human skin diseases in recent years. However, the investigation of human skin tissues has received only limited attention, to the point that there are not yet satisfactory modern detection technologies to accurately, noninvasively, and rapidly diagnose human skin at epidermis and dermis levels. In order to detect and analyze severe skin diseases such as melanoma, a finite element method (FEM) simulation study of the application of the laser-generated surface acoustic wave (LSAW) technique is developed. A three-layer human skin model is built, where LSAW’s are generated and propagated, and their effects in the skin medium with melanoma are analyzed. Frequency domain analysis is used as a main tool to investigate such issues as minimum detectable size of melanoma, filtering spectra from noise and from computational irregularities, as well as on how the FEM model meshing size and computational capabilities influence the accuracy of the results. Based on the aforementioned aspects, the analysis of the signals under the scrutiny of the phase velocity dispersion curve is verified to be a reliable, a sensitive, and a promising approach for detecting and characterizing melanoma in human skin.

Simulation study of melanoma detection in human skin tissues by laser-generated surface acoustic waves 1 Introduction According to figures of the World Health Organization corresponding to global mortality statistics in 2012, the increasing level of environmental pollution has become a major environmental health risk in big cities worldwide, causing about 7 million deaths annually, nearly one in eight deaths, and twice of those in 2008. Furthermore, according to the statistics of the American Cancer Society, skin cancer accounts for nearly half of all cancers in the United States. Melanoma, the most serious type of skin disease, accounts for more than 9000 of the 12,000-plus skin cancer deaths each year. 1 Based on the current situation, early diagnosis of human skin diseases in the skin surface, as a main medical mean, becomes more and more significant. However, the investigation of human skin tissues has received only limited attention, to the point that there are not yet satisfactory modern detection technologies to accurately, noninvasively, and rapidly characterize and diagnose human skin. For example, the most common methods of skin disease detection so far are still palpation and clinical section analysis. 2,3 The former is limited by subjective experience, and the latter is a destructive testing not proper for preventive inspection.
The application of ultrasonic and photoacoustic detection for medical testing has received considerable attention in recent years due to its advantage of noninvasive application. 4,5 Ultrasonic and photoacoustic detection are usually used for identification, characterization, or monitoring of tissues in vivo, while the investigation of surface tissues like skin has received only limited attention with still no satisfying outcomes so far. In order to meet the increasing demand for noninvasive testing of human skin tissues and for early diagnosis of skin diseases, a finite element method (FEM) simulation of the application of the laser-generated surface acoustic wave (LSAW) technique on a basic study of melanoma detection in human skin is developed in this paper.
The LSAW technique uses pulses of laser beams to generate ultrasound waves by the action of a thermo-elastic effect in the surface of solid media, and uses contact or noncontact detection approaches to acquire acoustic waves propagated in the surface vicinity of the source, with the potential to characterize some mechanical and physical properties of surface tissues like skin. [6][7][8] When a pulse laser beam impinges on the skin, a part of the energy is scattered, but a significant part is converted into heat. With this conversion process, elastic waves are mainly caused by local thermal expansion during temperature rise and spreading in the skin surface. Under this excitation mode, called thermo-elastic mechanism, the energy of the generated surface acoustic waves (SAW) is concentrated at the depth of approximately one wavelength of the elastic wave below the surface. 8,9 Since the characteristics of SAW propagation are directly related to the properties of the surface material, the changes that happened during the propagation process can be used as a detection index, reflecting the local differences of tissue's mechanical behavior, which may be a useful indicative of the presence of skin diseases, as Fig. 1 shows. Indeed, this principle is similar to that of manual palpation.
With more and more attention to skin diseases in the medical field, some research teams have tried to explore the interactions of pulse LSAW on skin by FEM. 6,10,11 A proper numerical model built by FEM may provide helpful guidance in stages previous to experiments, so it is usually a first step in achieving the goal of nondestructive and noninvasive experimental detection of skin diseases. As far as we know, this paper presents the first FEM simulation investigation of the application of LSAW spectroscopy to find and analyze skin cancer such as melanoma. Specifically, FEM is used here not only to build models for SAW generation, propagation, and analysis of the changes in human skin with melanoma, but also for the acquisition of SAW signals containing the information of mechanical properties of tissues. Time domain analysis and an especial emphasis on the utility of frequency domain analysis of the waves are used to analyze SAW dispersion and identify melanoma. Based on the analysis and discussion of the simulation results, the dispersion curve of LSAW is proposed to be a feasible tool to detect melanoma in skin.

Principles of FEM Simulation of LSAW Propagation
The basic procedure for the simulation of the laser-generated thermo-elastic mechanism is shown in Fig. 2. First, the parameters of the laser such as intensity, spatial and time distribution should be defined by using an appropriate input function. Through the absorption and scattering of laser by the skin, a local heat source is generated. Based on the heat conduction equation and the thermal properties of the material, the temperature rise and the thermal expansion cause thermal strain; and finally, the bulk and surface waves are excited by the thermal deformation of the material.

Modeling of Laser Irradiation
Considering that the incident laser beam has a Gaussian spatial distribution nature, and that there is an attenuation of the laser light caused by the optical scattering and the absorption properties of the skin and the energy conversion, the laser fluence ϕðr; zÞ can be described using Beer-Lambert's law: 12,13 ϕðr; zÞ ¼ E 0 fðrÞ exp½−ðμ a þ μ s Þz; (1) where E 0 is the irradiance on the skin surface, fðrÞ is the Gaussian spatial distribution function, μ a and μ s are the absorption and the scattering coefficients of the skin, respectively, z is the co-ordinate that describes the depth below the surface, r is the radial coordinate, and r 0 is the beam radius. If we use the heat source concept to describe the rate of heat absorption of the laser irradiation in skin, then the rate of heat generation per unit area Sðr; zÞ is given as the product of the absorption coefficient of the skin μ a and the laser fluence 13 Sðr; zÞ ¼ μ a ϕðr; zÞ: (3) Taking into account the temporal distribution gðtÞ of the laser irradiance, the heat source described as a volumetric heat generation to be used later in FEM simulation is then written as Sðr; z; tÞ ¼ μ a ϕðr; zÞgðtÞ. (4) Equation (4) represents the behavior of a distributed input, which is feasible to be used as a load in the FEM model. The temporal distribution gðtÞ is written as

Finite Element Analysis
Finite element models provide an approximate solution to the heat conduction equation, which is solved for a series of individual nodes of a mesh using shape functions. For the laser irradiation heat transfer, the heat loss by convection and radiation can be neglected. 13 Then, the classical heat conduction equation can be described as follows ½KfTg þ ½CfṪg ¼ fSg; where ½K, the conductivity matrix, and ½C, the heat capacity matrix, are determined by the material properties. fTg is the temperature vector, fṪg is the temperature rise rate vector, and fSg is the heat source vectorfound in Eq. (4). On the other hand, the finite element equation of wave propagation is described as ½MfÜg þ ½KfUg ¼ fFg; (7)  where ½M, the mass matrix, and ½K, the stiffness matrix, are the generalization of mass and stiffness parameters in generalized coordinates. fUg is the displacement vector, fÜg is the acceleration vector, and fFg is the external force vector.
In the theory of thermo-elasticity, the external force vector is the same as the thermal stress, which can be expressed as well as where V ε is the volume of thermal strain, fεg is the thermal strain vector, ½B T is the transpose of the derivative of shape functions, and ½D is the matrix of the material properties. It should be noted that Eqs. (6) and (7) are both established at the same particular time t. Then, by integrating over the whole simulation time, the entire time history of the generalized displacement U is obtained and can be acquired by picking points in the model to further apply the LSAW method.

Model of Human Skin
Human skin has a multilayer structure; usually it is divided into three lamellar structures in the medical field, the epidermis (0.07 to 0.12 mm), the dermis (0.3 to 4 mm), and the subcutaneous fat (usually more than 4 mm), respectively. 14 Therefore, in this research, the skin model is assumed to consist of three elastic and isotropic homogeneous layers in parallel, each layer corresponding to the above-described structure. To simplify the theory study and the calculating method, the small pores, glands, and capillaries are ignored; only part of their effect is reflected in the skin parameters like density and specific heat, etc. In addition, considering that SAW propagation can be easily analyzed by decomposing its movement only in a plane without transversal component, only a two-dimensional model is built. The model dimensions are 8 mm in length, with an epidermis depth of 0.08 mm, dermis depth of 1 mm, and subcutaneous fat depth of 4 mm, as shown in Fig. 3. In each of the three layers, the optical, the physical, and the thermal properties are different to make it more realistic. The thermal and the mechanical properties of each layer used in the model 14 are written in Table 1.

Model of Skin with Melanoma
The diagnosis of melanoma in the surface or subsurface of the skin is commonly carried out by detecting the changes in the local mechanical behavior of the tumor compared with its surrounding tissue, such as an increase in the Young modulus of squamous cell carcinomas and malignant melanocytes. 15,16 The closer to the core of the melanoma, the more serious and harder the localized sclerosis of the skin is. Based on this condition, the model of an early-stage melanoma is proposed as a concentric five-layer structure in the subsurface of the skin, whose radius is 1 mm and located at 2 mm from the left margin of the model, and with a Young's modulus that increases gradually toward its core. As Fig. 4 shows, the Young's modulus of the five layers of the melanoma from the outermost one to the core are assumed to be 1.5 × 10 5 , 1.8 × 10 5 , 2.1 × 10 5 , 2.4 × 10 5 , and 3 × 10 5 Pa, respectively. Other parameters of the melanoma model, which usually do not change so dramatically as the Young's modulus, are assumed to be the same as the parameters in the epidermis.

Model of the Load
In order to simplify the model and reduce the simulation time, a symmetric model of the load signal is assumed. Only half of a spatial distribution of laser beam load is set on the left boundary of the model, disregarding the other half load and its corresponding part of the skin model. The load signal parameters for the model were taken from a real pulse laser usually used for SAW experiments with the following parameters: Nd:YAG pulse laser, 532-nm wavelength, 50-mJ pulse power, 10-ns rise time, and 0.5-mm beam radius. The absorption and scattering coefficients were 0.05 and 2.5 mm −1 , respectively. Figure 5 shows the spatial distribution and the temporal distribution of laser beam load after the laser beam impinges the surface. Fig. 3 The three-layered skin model.

Meshing and Time Step Set
Temporal and spatial resolution of the finite element model is critical for the convergence and computational efficiency of the numerical results. In particular, meshing is an important step in FEM, since a proper meshing can help achieve more realistic results without making the model heavier in unnecessary calculations. For the skin model, on one hand, taking into account that the minimum thickness is 0.08 mm and the property parameters of the three layers will not change along the sample; and on the other hand, as SAWs propagate mainly at the surface of the highly absorbing skin, carrying significant energy in comparison with other bulk waves, only tiny displacements can pass through the deepest layer; so meshing by a small node meshing size at the epidermis and the dermis layers, and a larger size at the subcutaneous fat layer, would be sufficient to get reliable outcomes and simplify calculations at the same time. For the present model, the mesh size of the uppermost two layers is set as 0.03 mm, and for the deepest layer as 0.5 mm. The meshing of the melanoma is set to have the same meshing size as the uppermost two layers. Figure 6 depicts the mesh result of skin with melanoma.
In order to ensure that most of the frequency components of the generated waves can be detected, the integration time step should be carefully chosen. The needed time step is related to the time that the fastest waves need to propagate between successive nodes in the mesh. A smaller time step will acquire more accurate information; however, it will also require more processing time. To begin with, the time step should obey the Nyquist-Shannon sampling theorem, but considering the accuracy of signal processing, here the sampling time is chosen to contain at least 20 sampling points in the minimal wave package cycle time. The relationship is expressed in Eq. (9) where f max stands for the highest frequency of the surface waves to be detected. The time step for analysis is set as 2 × 10 −9 s during the input pulse duration, then increased to 2 × 10 −8 s during the propagation, so the value of the highest frequency expected to be detected is 2.5 MHz. The total time is 1 × 10 −4 s. An HP ProLiant DL580 G7 computer with 8G memory and two Xeon E7 processors was used for the simulation.

Simulation Results of the Skin Temperature Distribution
Assuming an initial temperature in the skin of 300 K, Fig. 7 shows that the heat gradient is propagating through the model, and the heat affected zone during the laser irradiation is localized. At the present laser simulation, the pattern of the heat penetration corresponds to the case when absorption predominates, and the laser energy is mainly absorbed at the surface of the skin, rather than being deeply penetrated into Fig. 4 The model of skin with melanoma. Fig. 5 The spatial distribution (a) and the temporal distribution of laser beam effect (b). the material, so the excitation efficiency of surface waves is strong enough to generate out-of-plane amplitudes of several subnanometers detectable by the simulation model or by a sensitive interferometer system in an experimental setting. 17

Simulation Results and Frequency Domain Analysis
Characterizing highly changing SAW signals only using time domain is not an easy task; therefore, a more proper analysis method is adopted. This approach uses the dispersion curve analysis in the frequency domain. 18,19 In order to calculate the phase velocity and plot the dispersion curve, two signals containing the information of the vertical displacements of the simulated SAWs in two different locations are needed. The two receiving points for detecting these signals are located at 3 mm (arrow A in Fig. 6) and 3.5 mm (arrow B in Fig. 6) from the left boundary of the model where the laser beam input signal is applied. Figure 8 shows the acquired signals y 1 and y 2 at the receiving points 3 and 3.5 mm, respectively. Frequency domain analysis is then used to obtain basic information of melanoma features in skin by analyzing the change of the dispersion curve. 10,20 The acquired signals y 1 and y 2 in time domain can be transformed into frequency domain signals Y 1 ðfÞ and Y 2 ðfÞ by applying the fast Fourier transform (Fig. 9). Then, these signals are further processed by calculating their corresponding angular phase ϕðfÞ [Eq. (10)], where 2nπ is used to ensure the continuous phase, so the phase velocity vðfÞ can be obtained by using Eq. (11) where Re UðfÞ and Im UðfÞ are the real and imaginary parts of the signal, x 1 and x 2 are the locations of the receiving points, and ϕ 1 ðfÞ and ϕ 2 ðfÞ are the phase of signals Y 1 ðfÞ and Y 2 ðfÞ. It should be noted that every time when the phase difference approaches or passes zero, which means the value ϕ 2 ðfÞ− ϕ 1 ðfÞ may be very close to zero, then based on Eq. (11), the phase velocity will present a big peak, which may influence the analysis of the results. Usually, this happens at the beginning of the measurements at low frequencies (usually below 200 kHz) where the value is very small and white noise exists, as well as at a higher frequency zone (about 2 MHz) where the maximum time step has been surpassed and the calculations are not reliable. The values of the phase velocity in this zone can take any random form, so it is reasonable to filter those zones from the plot and concentrate on the analysis in the useful part of the working bandwidth. After filtering the unreliable frequency zones, the dispersion curve, showing the relation of the phase velocity vðfÞ and the frequency f, is displayed in Fig. 10. From Fig. 10, it can be seen that the dispersion curve has an overall upward trend; a green dotted line is plotted to emphasize it. This is because when the SAWs propagate in the layered structure, if the Young's modulus of the upper layer is bigger than that of the lower structure, the phase velocity will gradually increase. This upward trend is in agreement with the results of previous simulation and experimental works of film materials. 10,20,21 The arched part of the dispersion curve occurred due to the existence of the melanoma, where the increasing Young's modulus pattern toward the melanoma core causes the surface waves to propagate from soft to hard regions, which will cause a significant increase in the dispersion curve. After the surface waves have passed the melanoma core, the Fig. 7 The contour plot of the skin temperature distribution. Fig. 8 The signals acquired at each of the receiving points at 3 and 3.5 mm. decreasing Young's modulus causes a decline in the dispersion curve. The arched pattern may be explained by the trade-off between the gradual rise and decline of the curve, which provides an insight into how the change of the mechanical properties caused by the melanoma will be directly shown in the dispersion curve. After the rise and the fall of the dispersion curve, the curve retakes the smooth rise curve, which reflects the detection of the mechanical properties of the uppermost normal skin layer. The significant changes in the dispersion curve of the skin with melanoma located in the range of 9 × 10 5 to 1.7 × 10 6 Hz (a red dotted line is plotted to make the pattern clear) are the result of the combination of the different melanoma properties. The lower-frequency surface waves are influenced more by the layered skin since low-frequency surface waves have a greater penetration depth; and higher-frequency waves are influenced more by the Young's modulus changes of top layer since their energy is more concentrated in the surface.

Influence of the Meshing Size
In order to get a higher sensitivity of the melanoma detection, a set of dispersion curves with different meshing size has been drawn to select a proper meshing size. While keeping parameters such as the model size, the total time and the time step unchanged, signals are acquired at the same picking positions depicted in Fig. 6, only the size of the meshing is changed, and then the corresponding dispersion curves are obtained, as shown in Fig. 11.
Three representative sizes of the meshing and their results are shown in Fig. 11. When the size of the meshing was 0.02 mm (green line), due to time-consuming calculations and lack of memory space, the program ran extremely slowly and was prone to present stops during the execution, so the size was chosen to be larger than 0.02 mm. Conversely, with a meshing size of 0.05 mm (red line), the result is not sensitive to the mesh size, the changes of the curve are equally obvious as when the size is 0.03 mm, while the curve became less refined, more averaged and rough; some valuable information may be lost from the tiny surface wave displacements and, consequently, may influence the results. Therefore, an appropriate mesh size for the present simulation conditions was set at 0.03 mm (blue line). If the computer allows more powerful processing and memory capabilities, a smaller mesh size will lead to more accurate results.

Minimum Size of Melanoma Detection
In order to evaluate the minimum size that can be detected by the simulation model, the model size, meshing size, and time step were fixed; only the radius of the melanoma was changed in the range of 0 to 0.5 mm. After testing some values of radii, a set of representative dispersion curves resulted as shown in Fig. 12.
The blue line shows the dispersion curve of the model based on a skin sample including a melanoma of 0.5-mm radius. There is a significant increase in the dispersion curve, so a melanoma with such radius can be perfectly detected in the simulation. Decreasing the radius of the melanoma to 0.1 mm (green line), it can be seen that the increase in the dispersion feature becomes roughly half of the former one, but it is still easily identifiable. If we keep decreasing the size of the melanoma radius step by step until the point where the dispersion curve seems to be only a small bump, finally we realize that this happens with a radius of 0.03 mm (red line). It is worthwhile to note that this value is the same as the size of meshing. After this threshold Fig. 10 The arched dispersion curve in the skin with melanoma. Fig. 11 The changes of dispersion curve with different meshing size values. Fig. 12 The changes of dispersion curve with different melanoma feature radii. radius value, there was already no difference with that of a skin model without melanoma (yellow line). Therefore, 0.03 mm appeared to be the minimum size of melanoma that the simulation model is able to detect. If we compare the resulting dispersion curves with the curve for the case without melanoma (yellow line), all the dispersion curves detecting melanoma have the same tendency to detect an arched pattern, with slightly small bandwidth differences for each case. Furthermore, because of damping loss, energy conversion, and the influence of reflections from bulk waves, small vibrations are shown during the whole process, with more dramatic changes at the beginning of the process or when the melanoma size is not large enough, which means the detection method is in its limit of sensitivity. Another reason for the high sensitivity to noise in the low-frequency region is due to the calculation method, which utilizes the difference of small values of the phase signals in the denominator of the phase velocity formula, values that can change dramatically due to the changes in the detected signals.

Discussion and Conclusion
In this paper, the research team developed a finite element model to simulate pulse laser beam irradiation on human skin, SAWs generation and propagation, as well as the computational procedure to process the changes of the dispersion curve of the surface wave phase velocity influenced by its pass through a melanoma model. We not only simulated the propagation of LSAWs in ideal pure skin, but also tried to go further and find an effective way of detection of the melanoma features in skin. The simulation results point out a promising potential use of the surface wave dispersion curve for early diagnosis of the existence of melanoma in skin. In a broader sense, the dispersion curve also characterizes some basic melanoma features since it is directly related to the mechanical and the physical properties of the tissue. Hence, the main conclusions are the following: 1. A three-layer model of human skin was established. The simulation model includes the effects of the excitation and propagation of SAWs, finding that the dispersion curve has a clear upward pattern, presenting a dramatic increase and decrease of the curve when the melanoma exists. This proves the simulation model was functional, and that the method of SAW spectroscopy applied to obtain the dispersion curve of SAW phase velocity as a tool for detecting and characterizing melanoma in skin is reliable.
2. It was found that the size of the FEM model meshing influences the simulation results, and depending on the conditions of the computing power, the meshing resolution can be improved, leading to more accurate simulation results. In the current research, due to computing restrictions, the recommended minimal mesh size was set as 0.03 mm.
3. By comparing the simulation results of the different dispersion curves, the current simulation model and the detection method could detect a minimum size of melanoma of 0.03-mm radius in skin model, which corresponded to the value of the size of the meshing, validating this way the reliability of the method. Therefore, it is reasonable to believe that the minimum size of melanoma that can be detected is restricted by the size of the FEM model meshing, meaning that the more powerful capabilities of a computer running the simulation model, the better resolution of the detected features in the human skin. The recommendation when having a break-even compromise between step time and meshing size is to set the meshing as the smallest real potential size the melanoma may be, and then, set the step time larger to compensate for the reduction of the meshing, but small enough to accomplish Nyquist-Shannon sampling theorem.
As further research, a broadband, sensitive, and accurate LSAW measurement system will be designed and built to verify present results on real animal skin such as pork skin, as a previous investigation stage before dealing with human skin. In addition, a parallel investigation will be conducted on pulse laser parameters to optimally generate SAWs on skin.