A Photoplethysmography Melanin Evaluation System by Modified Boltzmann Transport Equation (BTE)

With the advance in cosmetic medical technology in recent years, more and more people get cosmetic medical treatments, especially skin whitening treatments. Nevertheless, people usually assess the effect of skin whitening products by vision, which is subjective and will be different from each person. To acquire the value of melanin concentration objectively, people need to go to cosmetic medical clinics. This will cause inconvenience to people. This paper develops a novel evaluation platform based on optical assessment methods, which employ different absorption and scattering properties to different wavelengths of light in human tissue to obtain melanin concentration. Moreover, this paper proposes a new method that compensates the interaction between epidermis and dermis to acquire the melanin concentration more accurately. The novel platform designed in this paper is smaller and consumes lower-power and smaller when comparing to other conventional devices in market.


Introduction
The number of medical devices used in the global market is expected to be much more in the future [1], while noninvasive treatments and skin care products are more accepted to the mass population due to their safety and efficiency and then result in the fast growing of the global cosmetic market [2].Among various techniques, photoplethysmography (PPG) is a noninvasive electrooptic, low-cost, and simple technique that can be analyzed to provide the information of physiological status and skin properties, such as blood oxygenation, melanin, blood volume, heart rate, and other features of a human [3][4][5][6].To evaluate the melanin concentration and other tissue properties, most bibliographies are usually based on Boltzmann transfer equation (BTE) [7][8][9][10][11][12], Beer-Lambert law simplified from the BTE [13][14][15][16][17][18], or empirical equations derived from diffuse reflectance model and Boltzmann transfer equation [19][20][21][22].Nevertheless, the light scattered from the dermis (where hemoglobin mainly exists) will interact with the tissue in the epidermis and the main absorber is melanin.This will cause the nonlinearity and nonindependence of contribution from each concentration.To overcome this problem, this study proposes a new method that compensates the nonlinearity of chromophores to calculate the melanin concentration.
Skin consists of three layers: the epidermis, the dermis, and the subcutis; most of the melanin exists in the epidermis [23].There are five layers in the epidermis, including basal layer, spiny layer, granular layer, clear layer, and horny layer [24][25][26][27].Melanin is generated from tyrosine in melanocytes.It has two forms; one is eumelanin, which contributes a brown-black color and is more general; the other one is pheomelanin, which contributes a yellow-red color [28][29][30].
The optical properties of tissue include multiple scattering, absorption, optical pathways in skin, and the penetration depth of light.The electromagnetic spectrum plays a significant role in interacting with biological systems [31].Initially, cells, fibers, and chromophores are specific relevance to the discussion about the propagation of light in skin, although skin is composed of various kinds of cellular level elements   [32,33].Mie scattering is used to represent the scattering property when the size of the constituent is approaching to the wavelength of the incident such as mitochondria, nuclei, and collagen fibers [31,34].Secondly, with a low absorption coefficient of the element, light can propagate longer distance without large attenuation [35].Absorption in skin over the visible and near-infrared spectrum is mainly caused by many kinds of chromophores: melanin, hemoglobin, water, bilirubin, and beta-carotene [13,21,[36][37][38][39][40].The absorption spectra of hemoglobin have been mentioned in many studies [37,38,41,42].Bilirubin and beta-carotene give the yellowish-olive color of an individual's skin due to their absorption band [40,[43][44][45].Thirdly, for optical pathways in skin, literature assumes that the centers of absorption and scattering are uniformly spread across skin layers, and anisotropy induced by preferential alignments of collagen fibers is insignificant [32,37,46,47].Finally, the penetration depth of light depends on the wavelength of light [32,48,49].
PPG has been studied and widely used in medicine because it is a simple, low-cost, and noninvasive optical technique that can be used to detect microcirculation of tissue [6,50,51].Appropriate filtering and amplification can be used to obtain nonpulsatile part and pulsatile part of PPG signal [50] for the following pulse wave analysis.Finally, according to the principles of near-infrared optical tomography, the design approach of the systems is discussed [13,15,48,[52][53][54][55][56].The main three methods are continuous-wave (CW), time-domain (TD), and frequency-domain (FD).CW uses sources with constant intensity to interrogate the tissue.The configuration of the system is the easiest and the cost is also the lowest.
The rest of this paper is organized as follows.In Section 2, we introduce the algorithm for obtaining tissue properties and the method for compensating the chromophore.In Section 3, a novel evaluation platform, Physiskin, is proposed and the algorithm introduced was realized in this platform.In Section 4, the experiments are illustrated and discussed.In Section 5, we conclude this work and point out the future work.

Algorithm
These steps are illustrated in Figure 1 and will be explained in each subsection.Initially, the Boltzmann transport equation is introduced to model the light propagation in the turbid media from Sections 2.1 to 2.4.Secondly, this study describes the method to calculate the tissue properties in Section 2.5.Thirdly, after obtaining the necessary parameter and concerning the constant contribution of water, the tissue properties of a subject are calculated as well as the concentration of each chromophore from a linear equation in Section 2.6.Finally, the Empirical Mode Decomposition method is employed to compensate the nonlinear relationship between the epidermis and the dermis for the melanin concentration measured in Section 2.7.

Modeling Light Propagation in Scattering Media.
In this study, we use the transport theory, also known as radiative transfer theory, from the previous works to describe the light propagation in turbid media.Boltzmann transport equation is applied to describe the scattering and absorbing events of light propagation in the media [8,9,[57][58][59][60].However, in most practical cases, this equation cannot be solved precisely.It is necessary to consider an approximate approach.Diffusion approximation, which is validated in the strongly scattering media, is one of the simplified approaches.This approach is valid in most biological tissues due to the very large coefficient in scattering than absorbing [59].The diffusion approximation is performed in the semi-infinite geometry to obtain the appropriate solution for the interrogation of biological tissues [59,60].

Boltzmann Transport Equation (BTE).
The Boltzmann transport equation can be used to properly describe the movement of particles in scattering and absorbing media [11,59,60].The propagation of electromagnetic waves in scattering media can be described with the Boltzmann transport equation as follows: where the radiance  is in units and is the vector pointing to the interesting direction.The linear scattering coefficient and the absorption coefficient represent the inverses of the mean free paths for scattering and absorption, respectively.Moreover, (Ω  ⋅ Ω), the normalized differential scattering cross section, satisfies The function (Ω  ⋅ Ω) in ( 1) and ( 2) is the normalized probability function of scattering a photon that propagates from direction Ω  into direction Ω.The source term (r, Ω, ) is the power injected into a unit solid angle centered on Ω in a unit volume at r.The parameter  represents the speed of photons in the medium which is equal to  = (3 × 10 8 m/s)/, where  is the refractive index of the medium, which typically equals 1.40 for biological tissue.As shown in (1), the BTE has both integral and differential terms, and its solution needs initial and boundary conditions to obtain (r, Ω, ).This makes the equation sophisticated and hard to be solved.Therefore, the diffusion approximation is used to simplify the BTE to obtain the solution that can be analyzed easier.

Diffusion Approximation in an Infinite
Medium.The diffusion approximation hypothesizes that the radiance can be presented as an isotropic fluence rate plus a directional flux when scattering is much stronger than absorption in the media [57,59,60].The equation of the radiance in (1) will be expressed in the following: where We assume the form of the source term is like the radiance in (3): With the assumption of diffusion approximation, we can take (3) and ( 5) into (1) to obtain a closed set of two equations ( 6) and (7), where  is the diffusion coefficient and  is defined as the average of the cosine of the scattering angle and is the effective transport scattering coefficient which is   = (1 − ): If we assume that the source term is isotropic meaning that q 1 = 0 and ignore the time derivative of J(r, ), which means that the time between photon collisions with scattering particles is much shorter than the variation of j(r, ) with time in the medium, with this assumption, take ( 7) into ( 6), and we rewrite (6) in the form of the photon-diffusion equation as With utilizing differential equation skills and assuming that intensity of source is fixed, defined as , the photondiffusion equation can be solved.The solution of fluence rate (r, ) in the macroscopically homogeneous and infinite medium is shown in the following: 2.4.Light Propagation in a Semi-Infinite Medium.It is important that these optical techniques for clinical, diagnostic, and therapeutic applications are usually noninvasive, which means the light source and the detector must be mounted on the surface of the tissue.This kind of measurement is called semi-infinite geometry [60].Nevertheless, the infinite medium geometry is that the source and detector are both immersed in the medium, which cause the solutions of the diffusion approximation to be violated near the surface.Hence, the boundary of tissue is unavoidable and the diffusion approximation must be concerned with the boundary condition to avoid enormous errors, such as 50% or more in the measurement of optical properties.
In the semi-infinite geometry, when the boundary condition is utilized in Fresnel reflections at the surface, the violation of diffusion approximation is partially recovered.There are two boundary conditions which can be employed to interpret the Fresnel reflections that bring about the refractive index mismatch at the tissue-air interface.One is the extrapolated boundary condition and the other one is the partialcurrent boundary condition.The extrapolated boundary and the partial-current configurations are depicted.

Advances in Electrical Engineering
Sound physical principles are the basis for both the partial-current and the extrapolated boundary approaches.On the one hand, the partial-current boundary condition is logically corresponding with the limitations of the diffusion approximation.On the other hand, the extrapolated boundary condition is based on the mathematical solution to the transport equation in general situation.Each has its own pros and cons.However, the simulations and experiment measurements in both boundary conditions show that they have similar results to the infinite medium results.Moreover, with some simplifications, the configuration of the partialcurrent boundary condition resembles the configuration of extrapolated boundary condition.Hence, the solution of the partial-current boundary condition can be connected with the solution of extrapolated boundary condition.We obtain the solution of fluence rate, shown in (11), with considering the boundary condition.The geometry of the extrapolatedboundary condition to interrogate the tissue is as follows: The variable is the transport mean free path for the first collision for a single photon.As mentioned before, in the strong scattering, the effective transport scattering coefficient   is much larger than the absorption coefficient .Therefore, we simplify  tr = (  + ) −1 ≈ (  ) −1 .Other variables in the equation of fluence rate are where  is the horizontal distance between the source and the detector and  eff is the effective reflection coefficient which is typically equal to 0.493 in the biological tissue.

Reflectance Calibration Model.
The optical properties vary from each part of human tissue.Only using simulation to confirm the measurement is not appropriate.Hence, this study uses the optical properties from other studies as reference data [7,32,61].This study chooses the cheek, similar to the forehead, as the measuring part due to the signal-tonoise ratio, practical using, and clear wave [62].We use the optical properties from bibliographies, the LED source intensity, and source-detector distance of platform to calculate the theoretical reflection intensity from (11).The theoretical reflection intensity depending on these parameters can be defined as where  is the source intensity of LED,  is the source-detector distance of platform,    is scattering coefficient in tissue,   is absorption coefficient in tissue, and  rel is relative refractive index from bibliographies.
In (13), the spectrometer (Ocean optics USB4000) is used to measure the LED source intensity of each wavelength.Due to the difference between the reflection intensity measured from the system in this study and the theoretical reflection intensity calculated from the reference, we fit the reflection difference by where (, ) is called the calibration coefficient.It is the specific wavelength of light. is the source-detector distance.
To calculate , we use the platform in this study to measure the reflection intensity,  ref (, ), of a yellow race human's cheek as the criterion for the calibration coefficient.With obtaining the theoretical reflectance calculated by (13) and  ref (, ) measured by the system, the calibration coefficient (, ) can be obtained.After obtaining the calibration coefficient (, ), all the data   (, ) acquired from the platform in the future will be calibrated as  ,cal (, ) =  (, ) ×   (, ) .
With (15), the tissue properties can be calculated.

Chromophores Concentration.
We assume that the chromophores giving the absorption coefficient are principally melanin, oxyhemoglobin, deoxyhemoglobin, and water in the spectral band window which this study adopts.The residual three unknown chromophores' concentration can be derived from ( 16).This equation has been shown in many studies for the determination of human tissue chromophores concentration in vivo [7,11,14,61,[63][64][65][66]: where   () is the extinction coefficient of a chromophore at a specific wavelength and   is the concentration of a specific chromophore.
In (16), the chromophores are all independent of each other, meaning that the chromophores will not have interaction with each other.Hence, the absorption coefficient can be the summation from each chromophore with its extinction coefficient.Furthermore, it is obviously observed that if we want to derive three chromophores' concentration, we need three different wavelengths which can produce three different absorption coefficients.Based on the absorption spectra, this study chooses 640, 700, and 940 nm of LED wavelength to minimize the influence of the unwanted chromophores and interrogate the tissue.The 640 nm is mainly for melanin, and 700 nm and 940 nm are for the oxyhemoglobin, which is at the left side of the isobestic point (800 nm), and deoxyhemoglobin, which is at the right side of the isobestic point.Therefore, the concentration of each chromophore can be extracted from the following:  ,640 nm =  HbO 2 ,640 nm  oxy +  Hb,640 nm  deoxy +  mel,640 nm  mel  ,700 nm =  HbO 2 ,700 nm  oxy +  Hb,700 nm  deoxy +  mel,700 nm  mel  ,940 nm =  HbO 2 ,940 nm  oxy +  Hb,940 nm  deoxy +  mel,940 nm  mel .

(17)
The experiment for each subject takes about 25 seconds, and the first three seconds will be deleted due to the transient time caused by the low-pass filter designed in the software.Then, take the physiological signal filtered and calibrated to acquire the tissue properties in time series.With these tissue properties, we can obtain the mean melanin concentration and dynamic hemoglobin concentration.Furthermore, the dynamic oxygen saturation in tissue can be calculated as

Compensation and Calibration of Melanin Concentration.
The concentration of each chromophore calculated from ( 16) is based on the independence of each other.From the configuration of optical path way [32], the light incident on the skin surface is propagated into the epidermis and dermis under scattering and absorption.Some portion of the scattered light coming back from the dermis to the epidermis may interact with the melanin in the epidermis again and then finally be emitted to the skin surface that can be detected by the sensor.Therefore, this propagation course causes the nonlinearity of ( 16).The absorbance effects in the epidermis and dermis are not easily separated [67].However, the concentration of oxyhemoglobin and deoxyhemoglobin should not have large variation at the same part of skin interrogated between individuals [68,69].From our experiment results, we have found out that the hemoglobin concentration we acquired does not have large difference between individuals.Furthermore, we find out a critical feature that can be used to compensate the melanin concentration.
We assume that more heartbeat power obtained from the physiological signal gives more blood information.The reason is that when less light is absorbed by the melanin in the epidermis, more light can propagate into the dermis and then the constituents of the physiological signal received contain more pulsatile information.Consequently, the heartbeat power can be used to compensate the melanin concentration.
Nevertheless, the physiological signal that the system received has too much noise and signal unwanted.Although high frequency noise can be filtered out by using lowpass filter, the artifact or some noise that comes from high frequency band with the intermodulation still cannot be removed.Thus, using Fast-Fourier transform (FFT) directly to acquire the heartbeat power may be wrong due to the complexity mentioned before.To obtain the signal without using the filter, we introduce a method called Empirical Mode Decomposition (EMD) [3] that can obtain the signal without high frequency noise and artifact and can acquire the correct heartbeat power.
EMD has been proposed as a method that can analyze time series data.It has extensive applications to extract noisy nonlinear and nonstationary signals.EMD decomposes a signal into a finite set of AM-FM constituents.These constituents are called Intrinsic Mode Functions (IMFs), which is a counterpart to the common harmonic function.An IMF is not limited to be a narrow band signal but can be modulated in both amplitude and frequency.
To implement the needed decomposition of functions, we first specify all the local maxima and the local minima in the signal ().Then, a cubic line is used to link all the local maxima to produce the upper envelope.The lower envelope is also produced by linking all the local minima.Both envelopes should include all the data between them.Define  1 as the mean of the envelope.The difference between the data and  1 is defined as ℎ 1 , a typical IMF, After extracting the first IMF, we specify  1 = () − ℎ 1 , which is called the residual.We continue the procedures to specify all local maxima and the local minima in the signal  1 and then obtain the second IMF, called ℎ 2 .By following the steps mentioned above, we obtain a set of IMFs and a residual.The criterion specified to stop sifting is that when there are less than two extreme values in defining the envelope.Finally, we can express the data () in a series by In all IMFs, we search the IMF in our data analysis that has the most heartbeat power by using FFT.The mean melanin calculated by ( 17) is then compensated with the power of the found IMF calculated by FFT and the compensation equation can be expressed as The square of IMF used to compensate the melanin concentration in (21) is to increase the effect of heartbeat power.With more heartbeat power measured of the individual, this individual is supposed to have less melanin concentration.The block diagram of Empirical Mode Decomposition (EMD) is illustrated in Figure 2. From the processes mentioned above, the average melanin concentration can be calculated and compensated.After obtaining the modified melanin concentration, this study transforms the melanin concentration into melanin index for easier differentiation with the same form like [70,71].The equation of transformation from melanin concentration to melanin index is shown in (22).Three states of the transformation are employed to fit the melanin index from 0 to 999 because the reference subject is not the lightest or the darkest person.Moreover, the Mexameter MX18 is employed to calibrate the melanin index acquired from the Physiskin and four subjects are investigated as the reference.The equation from MX18 to calibrate the melanin index is in (23) where  mel,sub is the average melanin concentration of the subject,  mel,ref is the average melanin concentration of the reference, and MI is melanin index.With melanin index, people can compare their melanin pigmentation with each other or make sure that the skin whitening products are surely effective to them.All of these algorithms are first realized on Matlab 2009b and implemented as an application program for demonstration in microprocessor.If the data is uploaded on the personal computer with 32-bit operating system and i7 CPU@2.67GHz, it takes about 30 seconds to show all final results that consist of average hemoglobin, heart rate, and melanin index.

Implementation
Due to large attenuation of light propagating in the medium, about 100 dB, and the noise that interferes with the signal, an appropriate design of the system is quite important to have good performance in subsequent analysis [72].The signal first passes a transimpedance amplifier (TIA), bandpass filters (BPF), a multiplexer (MUX), and a constant gain amplifier (AMP).Later, the signal is converted to digital signal by an 8-bit analog-to-digital converter (ADC) inside the microprocessor and then sent to the computer for data analysis.The schematic of circuit blocks is shown in Figure 3.The power supply used is 5 V and the source voltage, which is 3.3 V and is not shown in Figure 3, is produced in the system by a low-dropout (LDO) regulator.To implement a high gain and good SNR system, the gain for each stage is arranged and the band of the filter is selected properly.In this study, the circuit of the system is implemented on the PCB.Figures 4(a

Comparison of Devices Performance.
To show the performance of the device implemented in this study, the device is compared with some commercially conventional devices.Although there are still many devices that can interrogate the skin and obtain the tissue properties, this study only takes devices that have similar system characters like Physiskin, such as being portable and easy to use, into account.Table 1 illustrates the comparison between devices [73][74][75].The platform, Physiskin, designed in this study has ultra-low power as compared with other devices.Melanin indices of the two devices are recorded in Table 2.Moreover, the correlation between the MX18 and Physiskin is about 0.987 and the average error is 4.85, only about 0.5% in the range of melanin index.It shows that the Physiskin can match with the MX18.

Conclusion
Based on the diffusion approximation simplified from BTE, this study designs a system that can evaluate the tissue to acquire melanin concentration of an individual.However, the equations used to obtain the concentration of each chromophore are based on the noninteraction between epidermis and dermis.This assumption is not enough because some light propagating through the dermis will scatter and go back to the epidermis and then will be received by the detector.Therefore, a new compensation method based on EMD is proposed in this paper and shows a good operation on compensating the melanin concentration.In addition to the PCB implementation, the IC solution of the system will be implemented to improve the performance of the system.
of platform R ref (, r) (, r) = R theory /R ref (, r)Compensation for C mel and → Diffusion approximation

Figure 1 :
Figure 1: Algorithm procedure to obtain tissue properties.

Figure 2 :
Figure 2: Block diagram of EMD for the compensation of melanin.

Figure 3 :Figure 4 :
Figure 3: Schematic of each block in the apparatus.
) and 4(b) show the front and the back of the PCB.The LEDs and PDs are all on the front.They were turned on and arranged in the same form as in the bread board.Other components and ICs are arranged on the back of PCB to avoid the uncomfortable condition during the experiment.
Concentration.Twenty subjects, 15 males and 5 females, aged from 20 to 29 years are measured in this experiment.Each subject's left cheek is measured by Mexameter MX18 and the Physiskin on the same part.

Table 1 :
Comparison between conventional devices and device in this study.

Table 2 :
Experiment result of melanin index from MX18 and Physiskin.