Multidiameter single-fiber reflectance spectroscopy of heavily pigmented skin: modeling the inhomogeneous distribution of melanin

Abstract. When analyzing multidiameter single-fiber reflectance (MDSFR) spectra, the inhomogeneous distribution of melanin pigments in skin tissue is usually not accounted for. Especially in heavily pigmented skins, this can result in bad fits and biased estimation of tissue optical properties. A model is introduced to account for the inhomogeneous distribution of melanin pigments in skin tissue. In vivo visible MDSFR measurements were performed on heavily pigmented skin of type IV to VI. Skin tissue optical properties and related physiological properties were extracted from the measured spectra using the introduced model. The absorption of melanin pigments described by the introduced model demonstrates a good correlation with the co-localized measurement of the well-known melanin index.


Introduction
The knowledge of the optical properties of biological tissues is invaluable for medical applications such as optical diagnosis and therapeutic laser procedures. [1][2][3] The optical properties are also the essential inputs when using a light transport model or Monte Carlo simulation to study the propagation of light within a tissue. 4,5 For many years, research has focused on using spectroscopic techniques in combination with dedicated mathematic models to extract optical properties of biological tissues. The model of diffuse reflectance spectroscopy (DRS) is generally derived from diffusion theory, in which the sampling depth is relatively deep (millimeter to centimeter) and path lengths of the detected light are relatively long. 6,7 To detect local, smallscale changes (<mm) in tissue, subdiffusive techniques, such as single-fiber spectroscopy (SFR), 8,9 are needed. In contrast to DRS, SFR employs a semiempirical model derived from Monte Carlo simulations. 10 SFR enables quantitative determination of the absorption coefficient without prior knowledge of the scattering properties. 11 In order to extract the scattering properties of tissue, two or more fiber diameters are needed. 3,12 Multidiameter single-fiber reflectance (MDSFR) spectroscopy measurements consist of SFR measurements using different fiber diameters d fib on the same tissue location. The measured reflectance spectrum of each fiber diameter R SF;d is the product of a reflectance spectrum without absorption R 0 SF;d and a Lambert-Beer term to account for the effect of absorption: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 6 3 ; 1 2 7 R SF;d ¼ R 0 SF;d ⋅ e −μ a hL d i ; where μ a stands for the absorption coefficient and hL d i stands for the effective pathlength of the detected photons using a fiber of diameter d fib . 10 To utilize the model mentioned above, it is assumed that the absorbers are distributed homogeneously within the sampling volume. With this assumption, the absorption coefficient μ a is independent of fiber diameters. Biological tissues, however, are far from homogeneous ( Fig. 1): blood is mainly confined in blood vessels; 13 melanin pigments are only present in the epidermal layer in the skin. 14 When the absorbers are distributed inhomogeneously, the absorption coefficient μ a is dependent on the effective pathlength and thus on the fiber diameter. Furthermore, the discrepancy between the homogeneity assumption and exact tissue architecture might result in misinterpretation of tissue optical properties and related biological parameters. Earlier work addressing a similar issue focused on the effects of the inhomogeneous distribution of blood. Correction factors have been derived based on an effective blood vessel diameter that allows accurate evaluation of tissue optical properties. 13,15,16 Using MDSFR to extract scattering and vascular properties of tissue has been intensively validated and utilized on less pigmented skin tissue, 3,9,15,17 however, applying MDSFR on heavily pigmented skin tissue might encounter the effect of the inhomogeneous distribution of melanin pigments, which is unaccounted for. Especially in darker skin types, melanin is the dominant absorber in the skin. 18 An additional complication with melanin in analyzing MDSFR spectra using the model described [Eq. (1)] is the fact that the absorption spectra of melanin pigments have the same shape as the commonly used model of the reduced scattering coefficient. 19 Both the reduced scattering coefficient and the absorption coefficient of melanin pigments decrease smoothly with increasing wavelengths in the visible wavelength range. This similarity leads to a competition between the contribution of the scattering of light and the absorption of melanin pigments when fitting the model to measured data. The competition may lead to unstable fit procedures with highly variable or unrealistic outcomes.
The accurate determination of the contribution of the absorption of melanin pigments in the skin is of great value for skin disease detection and skin color modeling. 20,21 To correctly account for the contribution of the inhomogeneously distributed melanin pigments on the measured skin MDSFR spectra and work toward the quantitative analysis of heavily pigmented skin, we introduce a layer model to describe the contribution of the absorption of melanin pigments [Eq. (2)]. The model assumes that all detected photons are independent of their eventual pathlengths and pass through a superficially positioned melanin layer twice: a first time when entering the tissue and a second time just before leaving the tissue. Thus the effect of the absorption of melanin pigments is the same for all detected photons by all fibers used in the MDSFR measurements and is thus independent of the pathlength of all detected photons: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 6 3 ; 2 3 where T layer is the transmission of light through the melanin layer. In this study, the introduced model is applied to the reflectance spectra obtained from in vivo MDSFR measurements on volunteers with heavily pigmented skin (skin type IV to VI). We demonstrate that this model allows the extraction of optical properties of the skin. A strong correlation is observed between the absorption of melanin pigments described by the introduced model and an independent measurement of the epidermal melanin concentration (EMC) on the same location, using mexameter, a commonly used EMC meter in the field. 20 2 Methodology

MDSFR Reflectance Model
The MDSFR reflectance model used as the basis of this study was derived previously by Kanick et al. 22 from Monte Carlo simulations. In their paper, the effect of the absorption on the measured reflectance R SF;d was expressed by a form of the Beer-Lambert law as in Eq. (1) but without the T 2 layer term included in Eq. (2). In Eq. (1), R 0 SF;d stands for the reflectance without absorption measured by the fiber of diameter d fib and can be expressed as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 3 2 6 ; 6 3 6 The lower limit of the single-fiber collection efficiency, η limit , is described in Eq. (4), in which NA is the numerical aperture of the measurement fiber and n medium is the refractive index of the medium in contact with the surface of the measurement fiber: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 3 2 6 ; 5 2 7 η limit ¼ NA n medium 2 : (4) η limit is ∼2.7% for a fiber of NA ¼ 0.22 immersed in water (n medium ≈ 1.33).
The parameter γ is related to the phase function that describes the angular dependence of scattered light. It is defined as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 3 2 6 ; 4 3 8 The parameters g 1 and g 2 are the first and second Legendre moments of the scattering phase function. g 1 is usually referred to as the scattering anisotropy factor. The reduced scattering coefficient μ 0 s is predefined using Eq. (6), in which μ 0 s ðλ 0 Þ is the reduced scattering coefficient at 700 nm and b is the scattering slope (b < 0): E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 3 2 6 ; 3 2 9 The effective path length hL d i of the detected photons by the fiber of diameter d fib is described by Eq. (7). The absorption coefficient μ a;d is dependent on the fiber diameter when chromophores are distributed inhomogeneously: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 3 2 6 ; 2 4 0 This MDSFR model was derived based on Monte Carlo simulations, in which a modified Henyey-Greenstein phase function (MHG PF) was assumed [Eq. (8)], where θ is the scattering angle. The first and second Legendre moments of an MHG can be expressed using Eqs. (9) and (10). 23 Combining Eqs. (9), (10), and (5), the parameter γ for the MHG PF is expressed by Eq. (11), where α ∈ ½0;1 was a factor that weights the contribution of the Henyey-Greenstein phase function (HG PF) P HG and an added Rayleigh-like scattering part ð1 − αÞ 3 4π cos 2 θ. In this study, α was fitted as a wavelength independent parameter. g HG is the first Legendre moment of the HG PF: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 6 3 ; 7 0 9 g 1 ¼ αg HG ; E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 6 3 ; 6 8 8 g 2 ¼ αg HG 2 þ 2 5 ð1 − αÞ; E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 6 3 ; 6 5 7 γ ¼ In the previous MDSFR studies, γ was fitted as a free parameter at each wavelength without a predefined model. 3,24 This approach significantly increased the degrees of freedom during the spectral fit, which might lead to overfitting. A predefined model of γ can avoid the possible overfit. In our study, γ was predefined with the model shown in Eq. (11) combined a predefined model of g HG . The predefined model of g HG was derived from a mathematic fit on the results of Mie calculation of a phantom introduced by Gamm et al. 9 to mimic a MHG PF, where λ is the wavelength (nm) and c is the fit parameter determining the exact shape of g HG (-): (12)

Melanin Layer Model
Based on the morphology of the skin, we have incorporated the following features in the representation of transmission of light through the melanin layer: (1) all melanin pigments are restricted to a superficially positioned layer of a thickness d mel ; (2) all detected photons propagate twice through the melanin layer, a first time when entering the tissue, and a second time just before leaving the tissue, furthermore, the transmission of light through the melanin layer was the same for all detected photons by all fiber diameters; and (3) only a surface fraction f of the layer contains melanin. Consequently, the incident photons can pass the other 1 − f fraction of the layer without encountering any melanin. It is assumed that the surface fraction sensed by fibers of different diameters is identical. When the photons propagate through the fraction of the layer containing melanin pigments, the transmission can be intuitively calculated using Beer-Lambert law as shown in E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 6 3 ; 2 6 6 T melanin ¼ e −μ a_melanin ⋅d mel ; where μ a_melanin stands for the absorption coefficient of the melanin containing fraction of the layer and d mel stands for the thickness of the melanin layer. When light propagates through the fraction of the layer not covered by melanin pigments, the attenuation of light due to melanin is zero, thus T melanin_free ¼ 1.
The effective transmission of the melanin layer can be calculated using the following equation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 6 3 ; 1 5 μ a;melanin is taken as the sum of the contribution of the absorption from two main melanin pigments: eumelanin (black and brown) and pheomelanin (red): 19 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 3 2 6 ; 7 5 2 μ a;melanin ¼ ðε eumelanin C eumelanin þ ε pheomelanin C pheomelanin Þ lnð10Þ; (15) in which C eumelanin and C pheomelanin are the concentration of eumelanin and pheomelanin (mg mm −3 ); ε eumelanin and ε pheomelanin are the corresponding extinction coefficients (mm −1 mg −1 mm 3 ). 25 The transmission of the melanin layer is then rewritten in Eq. (16), where the surface density of eumelanin and pheomelanin (mg mm −2 ) κ eumelanin and κ pheomelanin are the product of the melanin layer thickness d mel (mm) and the concentration C eumelanin and E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 3 2 6 ; 6 2 3 The surface density of melanin pigments is used as the fit parameter rather than thickness and concentration because fitting both the layer thickness and the concentrations independently would lead to a dependency between these parameters.

Combination of the MDSFR Model and the Melanin Layer Model
To correct for the inhomogeneous distribution of melanin pigments, the melanin layer model and the MDSFR model were combined in our study. The absorption of the other components in the tissue was still described using Beer-Lambert law with an effective path length described by Eq. (7). The reflectance R SF;d obtained from a fiber of diameter d fib was then described in E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 3 2 6 ; 3 9 2 T layer was the effective transmission of the melanin layer described in Eq. (16). R 0 SF;d and hL d i were different for different fiber diameters while T layer was independent of the fiber diameter. Note that μ 0 a;d the absorption coefficient of the sampling volume by the fiber of diameter d fib was fiber diameter-dependent because μ 0 a;d accounted for the absorption of blood, which was distributed inhomogeneously and concentrated in discrete cylindrical vessels. In this study, μ 0 a;d was defined as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 8 ; 3 2 6 ; 2 7 0 where μ a;d was the absorption coefficient of whole blood. StO 2;d was the blood oxygen saturation. μ HbO 2 a and μ Hb a were the micromolar absorption coefficients of oxygenated and deoxygenated hemoglobin (mm −1 μM −1 ). C hemoglobin;d was the total concentration (μM) of hemoglobin and F cor;d was a correction factor that accounts for the influence of the inhomogeneous distribution of the blood, which was determined by μ a;d and effective blood vessel diameter D v;d (μm) using the following equation: 11 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 3 2 6 ; 1 1 2

Experimental Setup
An MDSFR spectroscopy system, depicted in Fig. 2, utilized two measurement fibers of 400 and 1000 μm. White light was emitted by a halogen light source (Ocean Optics, HL-2000). Light from the source was guided to two digitally controlled optical shutters (Ocean Optics Inline Shutter) and connected, together with two fibers from the two-channel spectrograph (Ocean Optics SD200), with the two measurement fibers using SMA-SMA connectors. The 400-and 1000-μm fibers were glued next to each other in the channel of a cylindrical aluminum tube of 5-mm diameter. The fiber and tube surfaces were aligned and polished at 15-deg angle to eliminate the internal reflection induced by the refractive index mismatch between the fiber and sample in contact. The whole system was automated electronically and controlled with a LabVIEW program.

Data Acquisition and Analysis
The probe was positioned on the skin of each volunteer with gentle contact to avoid optical property changes due to the pressure applied. Water was applied between the probe and the skin to ensure a good optical contact. The absolute reflectance of the skin was calculated using E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 6 3 ; The absolute reflectance of the undiluted Intralipid 20% sample R IL (the reference sample) was determined using the Fresnel reflection method 26 with flat polished 400-and 1000-μm fibers from the same batch of the fiber used in the probe. The absolute reflectance spectra of undiluted Intralipid 20% of the 400-and 1000-μm fiber were saved in the LabVIEW program prior to the in vivo skin measurements. The reflected intensity of undiluted Intralipid 20% I IL and background intensity I b were measured prior to the measurements on each volunteer to account for the fluctuation of the output of the light source and the ambient light. I b was measured by immersing the probe in water in a black container.
This study was approved by the University Medical Center Groningen (METc M16.199148, UMCG, Groningen, The Netherlands). In vivo MDSFR measurements were performed on the inner forearm skin of 12 volunteers (November 9 to 11, 2016). Three spectra were measured from three close locations on the skin of each volunteer. Each spectrum was the average of a series of ten measured spectra of the same location. The standard deviation was calculated based on the same 10 measured spectra and used as an indication of the measurement error of the averaged spectra. Co-localized measurements of the reflected intensity of the skin I skin were performed on the same skin area of all volunteers where MDSFR measurements were performed using a commonly used EMC meter, Mexameter ® (Courage + Khazaka electronic GmbH, Cologne, Germany). The Mexameter ® emits three specific wavelengths (568, 660, and 870 nm). The melanin index (MI) was presumably determined using a proprietary algorithm. 20 We analyzed the MDSFR spectra from 400 to 900 nm in step of 1 nm. The nonlinear least-squares fit based on MDSFR model [Eqs. (3)- (19)] was performed simultaneously on the spectra obtained from the 400-and 1000-μm fibers. The following parameters were fitted: • the reduced scattering coefficient at 700 nm: μ 0 s ðλ 0 Þ (mm −1 ), • weighing factor of the contribution of Henyey Greenstein phase function (HG PF) and added Rayleigh like scattering part: α (-), • a parameter determining the exact shape of g HG : c (-), • the surface density of eumelanin and pheomelanin: κ eumelanin and κ pheomelanin (mg mm −2 ), • surface fraction of the melanin layer: f (-).
As mentioned previously, the vascular parameters were fiber diameter-dependent and fitted separately from the spectra of the 400-and 1000-μm fibers: • the concentration of hemoglobin: C hemoglobin;d (μM), • the oxygen saturation: StO 2;d (-), • the vessel diameter: D v;d (μm).
The measurement error of the spectra at each wavelength obtained from each fiber was used as a weighing factor in the fit so that spectral regions with large standard deviations have less influence on the fit. 17 The fit algorithm minimizes the reduced chi-square χ 2 red the difference between measurement and model, weighed by the measurement error: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 1 ; 6 3 ; 5 0 1 χ 2 red ¼ and R i;j , σ i;j , and M i;j stand for the reflection measurement, the measurement error and the model value at wavelength i of the spectrograph with fiber number j, υ for the number of degrees of freedom calculated from υ ¼ n − m − 1, where n is the total number of measurements (i.e., wavelengths in all spectra of all fibers) and m for the total number of model parameters in the fit. The model was used to fit the measurement spectra by varying the model parameters until the lowest possible value for χ 2 red was obtained. The confidence interval of each fitted parameters was estimated using the approach described previously by Amelink et al. 27 The 95% confidence intervals represented the statistical error as the square root of the diagonal elements of the covariance matrix multiplied by 1.96, obtained by multiplying the inverse of the second derivative matrix of χ 2 with respect to its free parameters by χ 2 ∕v.

Results
A typically measured spectrum of volunteer 8 (skin type VI) is plotted in Fig. 3 together with the respective best fit with χ 2 red ¼ 1.009. The measured reflectance spectra of all volunteers fitted well, resulting in a fit residue of less than two times the standard deviation of 10 measurement spectra at each wavelength (skin type IV to VI).
The reduced scattering coefficient at 700 nm μ 0 s ðλ 0 Þ and the scattering slope b were fitted from each measurement. The average of μ 0 s ðλ 0 Þ and b of each volunteer was calculated. The average μ 0 s ðλ 0 Þ ranges from 1.09 to 2.45 mm −1 , whereas b ranges from −3.32 to −1.53. The co-localized MI and the respective skin type are also shown in Table 1.
The reduced scattering coefficient was predefined as Eq. (6). The average reduced scattering coefficient spectra of each volunteer, calculated based on the average μ 0 s ðλ 0 Þ and b shown in Table 1, are depicted in Fig. 4 The reduced scattering coefficient ranges from 3.36 to 15.70 mm −1 at 400 nm and it ranges 0.43 to 0.95 mm −1 at 900 nm.
The averaged g HG spectrum of each volunteer at each wavelength was calculated and plotted [see Fig. 5(a)] based on the average of the fitted c and Eq. (12). g HG ranges from 0.92 to 0.99 at 400 nm and 0.60 to 0.95 at 900 nm. The fitted parameter α ranges from 0.87 to 1.00. Combining the anisotropy factor g, α, and Eq. (11), the spectra of γ were calculated and plotted [see Fig. 5(b)]. Both the amplitudes and the slopes of γ spectra vary among volunteers. The value of γ increases with the increase of the wavelength except in volunteer 6. Gamma values range from 0.79 to 1.93 at 400 nm and from 1.14 to 1.80 at 900 nm.
The parameters related to melanin pigments, the surface density of eumelanin κ eumelanin and pheomelanin κ pheomelanin (mg∕mm 2 ) and surface fraction f, were generated from the fit of the measured spectra. The average was calculated based on the fit results of the three measurements on each volunteer (see Table 2). Surface fraction f of the melanin layer varies from 0.70 to 1.00. Furthermore, the transmission spectra of the melanin layer T layer were calculated using Eq. (16) and plotted (see Fig. 6). The one-way transmission of the melanin layer at 700 nm was demonstrated together with the MI values measured by the mexameter in Table 2. The 95% confidence intervals (CI) of the fitted surface density of eumelanin and pheomelanin and surface fraction are also shown. The CI values of the fitted eumelanin concentration are very low (0.00 represents a value much smaller than 0.01), whereas the CI values for pheomelanin concentration and surface fraction are relatively high for some volunteers. The average fitted vascular parameters of each volunteer of both fibers: oxygen saturation (-), hemoglobin concentration (μM), and vessel diameter (μm) are shown in Table 3. For the 400-μm fiber, the fitted oxygen saturation of all volunteers is mostly below 0.47 except volunteer 11 (0.95); the fitted hemoglobin concentration ranges from 6 to 222 μM; the fitted vessel diameter ranges from 4 to 200 μm. For the 1000-μm fiber, the fitted oxygen saturation of all volunteers is mostly below 0.39 except volunteer 2 (0.72) and 5 (0.74); the fitted hemoglobin concentration ranges from 4 to 21 μM; the fitted vessel diameter ranges from 8 to 132 μm.  19 The model generates the single-fiber reflectance without absorption R 0 and the single-fiber reflectance R abs , which considers the absorption within the top layer. The transmission of the layer for different fiber diameters and scattering properties are plotted in Fig. 7. It is observed that the layer transmission as a function of the fiber diameter is saturating at higher Fig. 4 The average reduced scattering coefficient of all volunteers.  (12) fiber diameters. The difference between the layer transmission measured by the 400-and 1000-μm fiber is not significant within the range of simulated optical properties; thus the transmission of the melanin layer can be approximated as fiberdiameter-independent in our study. We do notice that the fiber-diameter-dependence of the transmission gets stronger for smaller fiber diameters, which indicates that the performance of the introduced layer model improves when using relatively large fiber diameters.

Strong Correlation between MI and the Transmission of Melanin Layer
Melanin is the dominant factor that determines skin color. The skin type is often categorized according the MI measured by the mexameter. In this study, the skin types of the volunteers range from Fitzpatrick skin type IV to VI. In this study, a strong correlation was found between the MI values and the total melanin surface density (the sum of the surface density of eumelanin and pheomelanin) ðR ¼ 0.90Þ. Furthermore, a strong correlation was found between the MI and the transmission of the melanin layer T layer ðR ¼ 0.97Þ as shown in Fig. 8. This correlation indicates that the layer model describes the absorption of light caused by melanin pigments in skin tissue as well as the mexameter. The fitted surface density of eumelanin of most volunteers was orders of magnitude higher compared to pheomelanin, which was expected since the skin colors of all volunteers were presented as black/brown rather than red. The correlation between the surface fraction f and the MI was low (R ¼ 0.47).  Journal of Biomedical Optics 127001-7 December 2019 • Vol. 24 (12) Quantifying melanin content using MDSFR with the introduced model does have advantages over DRS since the diffuse approximation works only when the scattering of light is dominant. For heavily pigmented skin tissues, the absorption is very likely to be dominant due to the presence of relatively high concentration of melanin pigments in a strongly inhomogeneous geometry. When the diffusion approximation does not work and/or absorption is more dominant, alternative methods are required to process the data, such as look-up tables. 29,30

Melanin Pigments Contribute to Scattering
The introduced model gives the transmission of the melanin layer (Fig. 6), a spectral shape that is different from the reduced scattering coefficient model (Fig. 4), which reduces the competition between the transmission spectra of the melanin layer and the reduced scattering coefficient during the fitting procedure. As a result, a significant decrease of the reduced chi-square is observed when analyzing the spectra using the introduced model (2.86 in average) compared with the fit when assuming a homogeneous distribution of melanin pigments (39.5 in average), which indicates a significant improvement of the fit quality. Jacques 19 reviewed previously conducted skin measurements: using the same model in Eq. (6), the fitted reduced scattering coefficients at 500 nm range from 2.97 to 6.87 mm −1 and the fitted scattering slope b ranges from −0.705 to −2.453. In this study, the average of the reduced scattering coefficients at 500 nm, which were calculated to be from  2.14 to 7.49 mm −1 , matches well the literature data, while the average of fitted scattering slope of volunteer 7 (−3.14), volunteer 8 (−3.32), and volunteer 11 (−3.32) were lower compared to the fitted scattering slope range reviewed by Jacques. 19 The low scattering slope b fitted in this study compared to the literature data might be caused by the scattering from the high concentration of submicron tissue components, 31,32 which, in our case, is the abundant melanin pigments within the optical sampling volume. Riesz 33 measured the scattering coefficient of a eumelanin solution. The scattering coefficient of the eumelanin solution can be well fitted with a predicted Rayleigh scattering coefficient obtained with a particle radius of 38 nm, 33 which is much smaller compared to the wavelengths of the halogen light source. 34 As a result, Rayleigh scattering (b ¼ −4), or small particle scattering, is more pronounced and presented as a low-scattering slope b. This hypothesis is supported by the strong correlation observed between the fitted scattering slope b and the sum of the surface density of melanin pigments (R ¼ 0.86) as shown in Fig. 9(b). In addition, we also found a correlation between the sum of the surface density of melanin pigments and the fitted reduced scattering coefficient at 700 nm, μ 0 showing the reduced scattering coefficient increases as the increase of the melanin surface density as shown in Fig. 9(a). A similar correlation was also found by Bashkatov et al. 35 who measured the reduced scattering coefficient of phantoms of different melanin concentrations. Their study showed the reduced scattering coefficient increases with the increase of the melanin concentration. The parameter α weights the contribution of Henyey-Greenstein and Rayleigh-like scattering. The smaller α is, the bigger the contribution of Rayleigh-like scattering is, or the smaller the scattering particle size is. However, we did not find a strong correlation between the fitted α and the sum of the surface density of melanin pigments ðR ¼ 0.36Þ in our study.

Fit γ with a Predefined Model
In previous MDSFR studies, the parameter γ at each wavelength was fitted individually without a predefined model. 3,24 This approach grants significant degrees of freedom during the fitting procedure and might lead to overfitting (i.e., the fitting of system noise). To avoid possible overfits, in this study, γ was predefined with Eq. (11) by assuming an MHG PF of the measured skin tissue. Earlier, MHG PF was also assumed when deriving the MDSFR model, 23 because it was believed to describe the back-scattering more accurately compared to HG PF. 10 To limit the degrees of freedom during the fit, g HG was also predefined with Eq. (12), which was derived from Mie calculations on a phantom introduced by Gamm et al. 9 to mimic an MHG PF. Gamm et al. found that a polystyrene beads suspension containing 10 different sphere sizes and a fractal dimension of 4.1 yielded the best match between the true phase function and best fit MHG PF. Using Mie calculations, the wavelength dependence of g HG was derived and expressed in Eq. (11). We would like to stress that the use of the model of g HG assuming an MHG PF is motivated by the reduction of the number of degrees of freedom of the fit. The exact phase function of the measured skin tissue, however, is not known, and we do not claim that we are able to measure the anisotropy factor of the skin tissue. Van Leeuwen-Van Zaane et al. 3 performed MDSFR spectroscopy measurements on murine skin. In their study, γ values were fitted to be roughly 1.2 at 400 nm and 1.4 at 850 nm and increased with increasing wavelength. Brooks et al. 24 utilized the same technique and performed in vivo measurements on Fig. 7 The dependence of the transmission of the melanin layer on the fiber diameter. Fig. 8 The correlation between MI and the average transmission T layer of the melanin layer at 700 nm.
Journal of Biomedical Optics 127001-9 December 2019 • Vol. 24 (12) human skin (skin type II to III) and γ was fitted and averaged. The averaged γ ranges roughly from 1.0 to 1.5. The fitted values of γ in our study range from 0.76 to 1.93. The fitted γ values mostly increase with increasing wavelength for all volunteers except volunteer 6. The descending trend of the fitted γ values for volunteer 6 was accompanied by an extremely high α value fitted (α ¼ 1). In this study, α was fitted as a wavelength-independent parameter. The fitted values range from 0.87 to 1.00, which indicates strong Mie scattering. However, Bashkatov et al. 36 measured a more dominant Rayleigh scattering in the shorter wavelength regime and the fraction of the Rayleigh scattering decreased with the increase of wavelength, which suggested that the ratio of the contribution of HG PF and Rayleigh-like scattering is very likely to change with the wavelength. 36 4.5 Higher Hemoglobin Concentration Fitted from the Reflectance Spectra of the 400-μm Fiber Vascular parameters: hemoglobin concentration, oxygen saturation, and vessel diameter were fitted individually for the 400and 1000-μm fiber to account for the inhomogeneous distribution of blood. Lower hemoglobin concentration and vessel diameter were fitted on the 1000-μm fiber spectra compared to the fitted values from the 400-μm fiber spectra, which is surprising since we expected the photons detected by the 1000-μm fiber to interact with a larger amount of blood compared to the 400-μm fiber. The CI values of the vascular fit parameters were relatively large compared to other fit parameters. The large CI values were due to the superficial sampling depth and heavy absorption of the melanin pigments, the vascular structures that interact with the detected photons were limited. Compared to SFR measurements in Caucasian skin, the absorption dips in the region 400 to 600 nm observed in this study are not well visible due to the dominant melanin absorption in this wavelength region, which may have interfered with the ability of the fit to accurately recognize the hemoglobin features in the spectrum.

Conclusion
A model is introduced to account for the inhomogeneous distribution of melanin pigments during MDSFR measurements. The model improves the quality of the fit and describes the MDSFR spectra of heavily pigmented skin well. We recommend utilizing the model with MDSFR measurements aiming at extracting tissue optical properties in human skin tissue abundant in melanin pigments, preferably with relatively big fiber diameters to minimize the fiber-diameter-dependence of the transmission of the melanin layer. We observe a correlation between the fitted scattering parameters and the surface density of melanin pigments, indicating that melanin contributes to the overall scattering properties of skin tissue. This work is the first step toward quantitative analysis of melanin pigment-related information and accurate extraction of other tissue optical properties in heavily melanin pigmented skin tissue. Future work will focus on full characterization of optical properties and physiological parameters related melanin pigments in skin types IV to VI.

Disclosures
No conflicts of interest, financial or otherwise, are declared by the authors.