Open Access
Add to BibSonomy
Abstract
Diffusion approximation (DA) of the radiative transport equation allows derivation of enclosed solutions for diffuse reflectance from multi-layer scattering structures, such as human skin. Although the DA is known to be inadequate near tissue boundaries and light sources, analytical tractability makes such solutions very attractive for use in noninvasive characterization of biological organs based on measured diffuse reflectance spectra (DRS). For the presented three-layer model of human skin, which enables a good match with DRS in visible spectral range measured with an integrating sphere, the DA solutions systematically overshoot numerically simulated DRS (using Monte Carlo approach) by 1–2 percentage points. However, using the former in inverse analysis of the latter can result in much larger artifacts, most notably overestimations of the melanin and blood contents by up to 15%, which must be considered when analyzing experimental DRS. Despite such systematic errors, the described approach allows simple and robust monitoring of physiological changes in human skin, as demonstrated in tests involving temporary obstruction of blood circulation and seasonal variations due to extensive sun exposure.
© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Diffuse reflectance spectroscopy is a practical and affordable technique for noninvasive characterization of biological organs in vivo. Diffuse reflectance spectra (DRS) depend on absorption and scattering properties of the involved tissues, and thus contain information about their chemical composition and structure. Their analysis was thus applied to diagnosis of various pathologies in bladder [1], colon [2], brain [3], breast [4,5], esophagus [6], and skin [7,8]. Specifically, DRS from human skin reflect the contents of epidermal melanin, oxy- and deoxy-hemoglobin, carotenoids, lipids, etc. [9–12], with potential applications ranging from monitoring of skin inflammation [10], changes induced by exposure to ultraviolet light [11,13], aging of traumatic bruises [14], and evaluation of scars [15].
In order to correlate DRS with optical and structural properties of the target organ, a mathematical model of light propagation in heterogeneous optically scattering structures is required. For most practical purposes, this is described with sufficient accuracy by the radiative transport equation. Due to its mathematical complexity, however, enclosed analytical solutions for practically relevant geometries can often be derived only within the so-called diffusion approximation (DA). Such DA solutions are well known to be inaccurate in the vicinity of tissue boundaries and light sources [16,17]. Nevertheless, their mathematical tractability makes them very attractive, especially when solving the inverse problem – assessment of selected properties of a target organ from experimental DRS.
In contrast with the above, numerical modeling of light transport using the Monte Carlo (MC) method [18,19] is very versatile and provides much more accurate predictions of DRS. Because MC results are intrinsically affected by stochastic noise, a large number of “photon” histories must be computed and averaged out in order to obtain reliable data. Moreover, independent MC runs must be performed for each wavelength of interest, so building up a broadband DRS is computationally very intensive. Luckily, this task is amenable to massive parallelization and can thus be performed in a reasonable time by employing a powerful graphics card (see Sect. 2.3). However, applying this approach to extraction of tissue properties from experimental DRS by means of iterative optimization (a.k.a. inverse MC) remains challenging and time consuming, especially when involving simultaneous assessment of multiple parameters.
To overcome this limitation, some researchers employed so-called look-up tables (LUT), which include pre-calculated reflectance values for a large number of tissue property combinations [20,21]. While such assessment of optical properties using the LUT is very fast and was proven to be useful under certain circumstances, it is practically limited to extraction of only a few free parameters (typically 2–3). Populating a LUT with significantly higher dimensionality, as would be required for characterization of multilayer structures, would namely require a prohibitive amount of time – or else the accuracy and reliability of the inverse solution would be compromised.
Most recently, several groups investigated also extraction of optical properties using artificial neural networks (ANNs), which also avoids the computationally intensive multidimensional optimization [22,23]. The main advantage of this approach is that it should be possible to train ANNs with a significantly smaller number of reflectance values as compared to filling the LUT. However, analysis of multi-layered media using ANNs have so far resulted in large prediction errors. This can be attributed to the intrinsic non-uniqueness in mapping of reflectance values to optical properties and inability to capitalize on prior knowledge such as spectral properties of the involved absorbers.
In this study we analyze first the accuracy of the DA solutions for diffuse reflectance in a three-layer model of human skin [24] which enables a good match with DRS measured using an integrating sphere, while keeping the number of variables sufficiently small. We find that the DA solutions overshoot numerically simulated DRS (using the MC approach) in visible spectral range by 1–2 percentage points, in good agreement with similar previous reports. However, as we demonstrate further on (to the best of our knowledge for the first time), much larger systematic errors can occur in the inverse analysis, i.e., assessment of skin properties by fitting the DA solutions to simulated DRS. Such artifacts most likely occur also in analyses of experimental DRS using the same approach, which should be considered in interpretation of the assessed values. Nevertheless, the same approach allows simple and robust monitoring of physiological changes in human skin in vivo when assessment of accurate absolute values is unnecessary, as demonstrated in our tests involving temporary obstruction of blood circulation using a blood-pressure cuff and seasonal variations due to extended sun exposure.
2. Materials and methods
2.1 Skin models
We model human skin with up to three optically homogeneous layers (representing the epidermis, dermis, and subcutis, respectively) separated by flat boundaries. Each layer is characterized by spectrally dependent absorption and scattering coefficients, and anisotropy factor.
The absorption coefficients of the epidermis and dermis, which vary with the fractional contents of melanin and blood, respectively, are computed as reported earlier [25]. In a nut shell, we apply the baseline absorption coefficient for both layers and absorption spectrum of the epidermal melanin from Jacques [26]:
The reduced scattering coefficient of the epidermis and dermis is assumed to be the same and is described as a combination of Rayleigh and Mie contributions [21,26,29,30]. Scattering properties of the subcutis can vary considerably between different subjects and anatomical locations due to differences in the adipocyte size and the amount of connective tissue septa [7,31]. In agreement with this, our preliminary analyses of DRS measured in vivo also indicated that we should allow for some variability with regard to the scattering properties of subcutis. To that end, we introduce a model function with varying amplitude A [24,32]:
Fig. 1. Reduced scattering coefficient spectra of subcutaneous fatty tissue from different authors: (a) Salomatina [7], (b) Bashkatov [33], (c) Bashkatov [34], and (d) Simpson [28]. Dashed lines represent our model function (Eq. (3)) for different values of amplitude A (see the labels).
Refractive indices of the characteristic skin layers vary between different sources [31]. In order to simplify the analytical (DA) expressions, we set its value to n = 1.4 for all skin layers. Our dedicated numerical simulations (MC) namely showed that differentiating these values between the three model layers has a negligible effect on the predicted DRS.
The parameters describing the structure and chromophore contents in our three-layer model of human skin are thus the epidermal and dermal thickness (depi, dder), melanin volume fraction (m), epidermal and dermal blood contents (bepi, bder), oxygen saturation of blood (S), and subcutis scattering amplitude (A).
2.2 Diffusion approximation (DA)
Treating light propagation in strongly scattering media within the DA leads to familiar diffusion equation for fluence rate distribution inside the tissue. However, a key condition for validity of the DA is nearly isotropic fluence rate throughout the treated volume. While this requirement is usually satisfied within the dermis, this is not the case near the skin surface, especially when irradiated with a collimated and/or spatially constrained light beam. Therefore, a special form of boundary condition at the air-tissue interface and dedicated “source functions” were derived to enable treatment of such situations [35–37]. The corresponding analytical solutions for fluence rate distribution and diffuse reflectance values were derived earlier for one-, two-, and three-layered tissue structures [37–40].
It is worth noting that these analytical solutions express diffuse reflectance relative to the effective light sources, which reside inside the tissue. In all experimental situations, however, diffuse reflectance is invariably considered as the ratio between the light flux exiting the tissue surface and the power incident on the same from an external light source. In order to remove this inconsistency, we thus account for partial reflection of the incident light at the air-tissue interface by multiplying the diffuse reflectance expressions in the reported DA solutions with (1 – RF), according to Fresnel’s formula for perpendicular incidence:
2.3 Numerical simulations
Numerical simulations using the MC technique are generally accepted as a valid model of light transport in strongly scattering media. We have applied the customary MCML software [18], where all energy packets (a.k.a. “photons”) are launched from the same point on the skin surface in the direction perpendicular to it. For realistic representation of our DRS measurements with an integrating sphere (see Sect. 2.4), only photons exiting the skin within a certain distance from their origin are accounted for in the computation of diffuse reflectance.
All presented examples involve launching of 106 photons at every considered wavelength, which provides sufficiently low numerical noise while not presenting an excessive computational load. (In our preliminary tests with 2×104 photons per run [24,32], the inverse analyses may have been affected by the noise clearly visible in simulated DRS). Upon massive parallelization of the calculations utilizing a graphics processing unit (NVIDIA GTX 770) and CUDA technology [19], computation of a typical DRS with an increment of 1 nm took ∼160 s on a personal computer with Intel Core i7, 3.4 GHz and 16 GB of DDR3 RAM.
2.4 Measurements of DRS from human skin in vivo
All DRS were measured using an integrating sphere (IS) with an integrated white-light source and a sample opening diameter of 10.3 mm. The reflected light was collected in a nearly perpendicular direction (8° from the normal) by an optical fiber connected to a compact diffraction spectrometer (ISP-REF and USB4000, respectively, by Ocean Optics, Dunedin, FL). Spectral response of the entire setup was calibrated using a white standard (Spectralon© by Labsphere, North Sutton, NH) with Lambertian reflectance exceeding 0.99 within the visible spectral range. The measurement artifact which arises from inevitable change of the light field inside the IS upon the substitution of the white standard with the sample (a.k.a. single-beam substitution error) was removed by performing additional measurements at the reference port of the IS [41].
DRS data were acquired from healthy volunteers with skin of fair complexion (Fitzpatrick types I–II; age 20–37), involving primarily the authors and a couple of close colleagues who have provided an informed consent. Tests involving temporary obstruction of cutaneous blood upon application of a blood-pressure cuff (inflated to 200 mm Hg for 5 minutes) were performed in three volunteers. In addition, monitoring of seasonal changes in human skin ran on a weekly basis from late June to early December in two volunteers (dorsal and volar side of the forearm, forehead). Dedicated transparent masks, marked with specific visual features of each test area enabled us to perform several measurements at the same test spot.
A study protocol allowing for such measurements (and some additional ones) in 30 subjects, was approved by the Medical Ethics Committee of the Republic of Slovenia.
2.5 Extraction of skin characteristics
The selected skin characteristics were assessed from either numerically simulated (MC) or measured DRS by seeking the best match with predictions obtained using the analytical (DA) solutions. This was performed by means of multidimensional minimization of the residual norm using the nonlinear least-squares algorithm, implemented as function lsqnonlin in the Matlab Optimization Toolbox (Mathworks, USA). The fitting increment was always 1 nm.
3. Diffusion approximation solutions vs. numerical simulations of DRS
3.1 Single- and two-layer skin models
Figure 2(a) shows the DRS computed using the DA solutions for a semi-infinite dermis with fractional blood content of bder = 1.0% and oxygen saturation level S = 75%. The results obtained by applying the isotropic or δ-Eddington source functions are presented by dashed and dotted line, respectively. Similarly, Fig. 2(b) presents the DRS computed for a two-layer skin model, where a 100 μm thick epidermis with fractional melanin content of m = 1.5% and a small amount of blood (bepi = 0.2%) is supported by a semi-infinite dermis with the same properties as above.
Fig. 2. DRS as predicted by DA solutions with two different source functions (navy lines; see the legend and text for details) and MC simulations (red solid lines) for: (a) semi-infinite dermal layer (with bder = 1.0%, S = 75%); and (b) two-layer skin model imitating the epidermis (depi = 100 μm, m = 1.5%, bepi = 0.2%) and dermis with the same properties as in (a).
By considering the corresponding MC results as a reference (MC; red solid lines), we can see that in both examples, the DA solutions based on the δ-Eddington source function (dashed lines) provide significantly more accurate results than those based on isotropic one (dotted). Even the former solutions, however, overshoot the numerically computed diffuse reflectance values by 1–2 percentage points across the entire spectral band considered in our analysis.
Similar comparisons were presented earlier by Spott et al. [38] and Randeberg et al. [39] for single- and two-layer skin models, respectively. However the discrepancies between the DA predictions and MC results in our data are somewhat smaller in comparison, due to the appropriate attenuation of the source functions in accordance with Eq. (4).
3.2 Three-layer skin model
In our three-layer skin model, which accounts for the finite dermal thickness and adipose subcutis under it, the discrepancies between the DRS predicted based on the DA solutions and MC simulations are very similar as in the previous examples. The DA solutions slightly, yet consistently, overestimate diffuse reflectance across the entire spectral range, both in the case of lower and higher melanin contents (Figs. 3(a) and (b), respectively). Between the two DA approaches, δ-Eddington source function again provides more realistic DRS values than the isotropic source function. The latter is therefore excluded from subsequent analyses.
Fig. 3. DRS as predicted by three-layer DA solutions with two different source functions (navy lines; see the legend) and MC simulations (red solid line). The epidermal melanin content is m = 0.5% (a) and 1.5% (b). The dermal thickness is dder = 1.0 mm and subcutis scattering amplitude A = 1.5; the remaining parameter values are the same as in Fig. 2(b).
Since all presented comparisons between the DA and MC results considered the same skin properties for both computational approaches, the observed discrepancies result exclusively from the intrinsic inaccuracy of the DA approach near the air-tissue boundary; specifically, for the case of collimated tissue irradiation treated using the approach derived by Haskell et al. [36].
4. Fitting diffusion approximation solutions to numerically simulated DRS
Because of the close similarity between the DRS predicted using the DA solutions (with the δ-Eddington source function) and numerical results for the same skin structure and composition (see Figs. 2 and 3), it is tempting to disregard the small systematic difference and assess skin properties by fitting experimental DRS with analytical DA solutions.
In this section, we thereby analyze the artifacts which may arise from such an approach. We perform this by extracting skin properties from numerically simulated DRS for various skin models with known structure and composition.
For realistic simulation of experimental DRS, we take into account also the finite sample opening in the IS. The applied DA solutions namely assume one-dimensional light transport, which corresponds to an infinitely large sample opening.
As evidenced by a dedicated MC simulation (Fig. 4), DRS collected using our IS with a 10.3 mm wide sample opening will indeed deviate from the result representing a (hypothetical) infinitely large opening. Because the discrepancy increases with the wavelength, some researchers refer to this effect as the “red loss” artifact. In order to prevent its adverse influence on analysis of experimental DRS (Sect. 5) we limit the fitting with DA solutions to λ ≤ 600 nm. The same protocol is therefore applied also in the following analyses of simulated DRS.
Fig. 4. Numerically simulated DRS, taking into account the finite diameter of the IS sample opening (orange line) vs. the hypothetical case of infinitely large opening (red).
4.1 Two-layer skin model
Figure 5 presents the results of fitting the two-layer DA solutions (with δ-Eddington source function) to numerically simulated DRS (MC) for the same two-layer skin model with either low or high melanin content (Figs. 5(a) and (b), respectively). Five skin model parameters are optimized in this example: depi, m, bepi, bder and S; whereby the latter applies to both epidermal and dermal blood.
Fig. 5. Numerically computed DRS for a finite IS aperture (orange solid lines) and best fitting analytical predictions (DA; dashed) for a two-layer skin model with (a) m = 0.5%, (b) m = 1.5%. Fitting was performed in the interval λ = 400−600 nm.
Within the fitting interval of 400–600 nm, the match between the MC spectra and best fitting DA prediction is evidently excellent. At longer wavelengths, the DA solutions increasingly overshoot the MC values, just as anticipated based on their one-dimensional nature (see Fig. 4).
The input parameter values used in the MC simulation are presented in Table 1 (columns MC). The assessed, best fitting values (fit) are complemented with the algorithm-provided estimates of the respective confidence intervals, and r marks the ratio of the latter and the former.
Table 1. Analysis of the skin parameter values as assessed by fitting numerically simulated DRS for a two-layer skin model with the corresponding DA solutions (see Fig. 5).
The last column displays relative deviations of the actual input values from the corresponding fitting results:
This metric will become relevant in interpretation of the results obtained when fitting the DA solutions to experimental DRS. Namely, inasmuch as our MC simulations produce valid representations of experimental DRS, δMC can then be used as a predictor for deviation of the true parameter value from the assessed one.The results in Table 1 indicate the largest estimates of the relative error (r) for epidermal parameters depi, m, and bepi. One could attribute this to the fact that the epidermis is rather thin and thus has a limited effect on DRS in the considered spectral range, which prevents reliable assessment of these parameters. In addition, it is plausible that the chromophore contents m and bepi are (inversely) correlated with thickness depi, since the DRS will depend primarily on the corresponding products (i.e., optical densities).
The results in Table 2 illustrate the effect of fixing depi to a set value; in this case the correct one. The estimated relative errors for the extracted epidermal parameters m and bepi are now significantly reduced, indicating enhanced stability of our inverse analysis. Meanwhile, the consequent reduction in the number of free parameters (from 5 to 4) has no adverse effect on the quality of the fits.
Table 2. Skin parameter values assessed by fitting numerically simulated DRS for a two-layer skin model with the corresponding DA solutions with four free parameters.
Nevertheless, the assessed parameter values still deviate from the respective input values by considerable margins. After having analyzed many similar examples we can safely state that the chromophore contents (m, bepi, and bder) systematically overshoot the input values, while the oxygen saturation (S) is usually somewhat underestimated.
4.2 Three-layer skin model
Adding the third model layer, representing the subcutis, introduces two additional parameters, dermal thickness (dder) and scattering amplitude (A, Eq. (3)). Figure 6 presents two examples of fitting MC-generated DRS for such skin model (with dder = 1.0 mm and A = 1.5) using the three-layer DA solutions. The remaining input values are the same as in Fig. 5.
Fig. 6. Numerically predicted DRS for a three-layer skin model (orange solid lines) and best DA fits (dashed). Fitting was performed in the interval λ = 400−600 nm. The epidermal melanin content is (a) 0.5%, (b) 1.5%.
Just as in the previous example, an excellent match between the simulated DRS and best fitting DA solutions is obtained within the fitting range (λ = 400−600 nm), and the latter overshoot the MC predictions at λ > 600 nm.
Despite the good match, however, the assessed values of parameters m, bepi, and bder (see Table 3) are again significantly overestimated, as evidenced by the large and negative values of the corresponding metric δMC. This is a direct consequence of the intrinsic overshoot of diffuse reflectance values as predicted by the DA solutions in comparison with the MC results (see Fig. 3). Meanwhile, the oxygen saturation values (S) are assessed rather accurately.
Significantly larger relative errors (r) were assessed for dder and A, as opposed to other model parameters. This indicates that the DRS obtained using an IS have lower sensitivity to these two skin properties. This is plausible since the boundary between dermis and adipose lies relatively deep inside the skin, and the absorption spectrum of the latter is rather featureless in the considered spectral range.
In addition, the larger values of r could also result from correlation between the parameters dder and A, similar to that between depi and m discussed in Sect. 4.1. Specifically, a thinner dermis will allow more light to reach the subcutis, where absorption in the involved spectral range is relatively weak and light is strongly backscattered. The resulting effect on DRS is thus similar to that of increased scattering amplitude A.
Table 3. Skin parameter values as assessed by fitting numerically simulated DRS for a three-layer skin model with the corresponding DA solutions.
5. Fitting diffusion approximation solutions to experimental DRS
5.1 Two-layer DA solutions
Figure 7 presents DRS as measured at a dorsal side of the forearm in a healthy subject before and after sun tanning (orange solid lines). In this subsection, we attempt extraction of key skin properties from these DRS using the DA solutions corresponding to our two-layer skin model. As discussed earlier (Sect. 4.1), the fitting process is limited to the spectral range of λ = 400–600 nm, and the DRS below 450 nm, where the SNR is rather low, is smoothed by sliding averaging (window width: 9 nm). Based on our experience from simulation tests (Sect. 4.1), we fix the epidermal thickness to depi = 101 μm, which is the average value for this anatomical location reported in literature [42].
Fig. 7. DRS as measured in a healthy subject with fair skin (orange solid lines) and best fitting two-layer DA solutions (dashed): (a) before, and (b) after sun tanning. Fitting was performed in the interval λ = 400−600 nm.
All assessed skin parameters fall within physiologically plausible ranges (Table 4). The values for S are in agreement with those reported by other groups using similar methodologies [15,30,43]. Perhaps most importantly, the results indicate a considerable increase of the melanin content upon extensive sun exposure, from m = 1.5% to 3.3%. The last line (ɛ) presents the quadratic norm of the residuum (i.e., difference between the experimental and best fitting model DRS vectors), an objective indicator of the quality of the match.
Table 4. Skin parameter values as assessed by fitting two-layer DA solutions to DRS acquired from human skin before (left column) and after sun tanning (right). ɛ marks the quadratic norm of the residuum.
However, Fig. 7 reveals a concerning mismatch between the measured and fitted DRS above 600 nm. Namely, the latter values are consistently lower than the former, which conflicts with our numerical simulation of the “red loss” artifact (Fig. 4) and experience with similar analysis of numerically simulated DRS (Fig. 5). This discrepancy warrants a closer inspection.
5.1.1 Influence of the epidermal thickness and presence of pheomelanin
For human skin on dorsal side of the forearm, analyzed in the discussed example, literature data indicate a considerable inter-personal variation of the epidermal thickness, 65−137 μm [42]. We have thereby tested first whether the qualitative discrepancy mentioned just above could arise from fixing the epidermal thickness in our inverse analysis to an incorrect value.
To that end, we have repeated the analysis of the experimental DRS presented in Fig. 7(a) with depi set to 65 μm. As can be seen in Fig. 8(a), however, the discrepancy above 600 nm was practically unaffected. We attribute this effect to the strong correlation between depi and m (see Sect. 4.1). Consequently, setting the depi to a lower value results primarily in a corresponding increase of m, with only marginal influence on the remaining variables. This was confirmed also by using depi = 137 μm, which resulted in very similar best-fitting DRS (graph not presented here) [32]. Based on such evidence, we can exclude the inaccurate epidermal thickness value as the primary origin of the investigated discrepancy.
Fig. 8. DRS as measured in a healthy subject before sun tanning (orange solid lines) and best fitting two-layer DA solutions (dashed navy lines): (a) with epidermal thickness set to depi = 65 μm; and (b) with the pheomelanin-to-eumelanin ratio of 0.50.
Our skin models use the customary melanin absorption spectrum [26], which corresponds to the prevalent form of melanin found in human skin, eumelanin. Some researchers, meanwhile, considered a mix of eumelanin and up to 30% of pheomelanin [44], which is prevalent in red hair, lips, etc. However, we have found that even by replacing 50% of eumelanin with pheomelanin, the peculiar underestimation of diffuse reflectance above 600 nm in best-fitting DA solution is only partially removed (Fig. 8(b)). Because such a high pheomelanin content is implausible for our test subject with brown hair (Fitzpatrick skin type II), potential presence of some pheomelanin in his skin can also be excluded as the source of the discussed discrepancy. In the interest of simplicity (and in line with most of related literature) such a possibility is not considered in the following analysis.
5.2 Three-layer DA solutions
As demonstrated in Fig. 9(a), applying our three-layer model of skin eliminates the implausible behavior of the best-fitting DA solution discussed just above. Instead, the latter features increasingly overestimated reflectance values at λ > 600 nm, exactly as seen in fitting of numerically simulated DRS (see Fig. 6). In addition, the match between the measured and fitted DRS is considerably improved as compared to the two-layer approach, especially between 540 and 580 nm (Fig. 7(a)). This is reflected also in the considerable drop of the residual norm, from ɛ = 30.8 to 11.2 (see left columns in Tables 4 and 5, respectively).
Fig. 9. DRS as measured in a healthy subject with fair skin (orange solid lines) and best fitting three-layer DA solutions (dashed): (a) before, and (b) after sun tanning. Fitting was performed in the interval λ = 400−600 nm.
Table 5. Skin parameter values as assessed by fitting three-layer DA solutions to DRS acquired from human skin before (left column) and after sun tanning (right). ɛ marks the quadratic norm of the residuum.
A marked improvement with respect to the two-layer approach can be seen also in analysis of the same test site after sun tanning (Fig. 9(b); compare with Fig. 7(b)), which also leads to a significant reduction of the residual norm (Tbl. 5, right column).
The assessed values for both examples fall within physiologically and anatomically plausible intervals. In perfect analogy with our analysis of simulated DRS (Sect. 4), the confidence intervals in Table 5 are generally somewhat larger as compared to Table 4 (due to the increased number of free parameters), and the relative error estimates (r) for dder and A are significantly larger than for the remaining parameters.
5.2.1 Modification of the fitting interval and linking of selected model parameters
Despite the evident improvement gained by application of our three-layer skin model, the best fitting DA solution in Fig. 9(b) still features the implausible undershoot of reflectance values at λ > 600 nm with respect to experimental DRS, albeit considerably reduced as compared to that seen in Fig. 7(b). In the following, we thus present three augmentations of the analysis approach, which help remove this deviation and at the same time lead to improved quality of the fit.
First, we introduce limited sensitivity of the fitting procedure to wavelengths beyond 600 nm. However, because the fitted DA solutions should be allowed to overshoot the experimental DRS in this spectral range, any difference between the two is accounted for in the merit function only when the model prediction is lower than the measured value. In this way, only the interim iterative solutions which undershoot the reflectance above 600 nm are penalized, while an overshoot doesn't influence the optimization process. To be conservative with this augmentation, only five points (from 610 to 650 nm) are treated in this manner.
Secondly, the fitting increment in the range of 400–450 nm is also increased from 1 nm to 10 nm in order to reduce the impact of this spectral range. The main rationale behind this step were low SNR in this part of measured DRS and the questionable validity of DA in skin at such short wavelengths, due primarily to their low penetration depths.
Finally, we take advantage of the fact that certain skin properties won't vary significantly over time. We can therefore analyze in parallel two or more DRS acquired from the same test site, while requiring that selected model parameters assume the same values across all respective solutions.
Figure 10 shows the same experimental DRS as in Fig. 9, fitted with three-layer DA solutions according to the augmented approach described just above. Specifically, the dermal thickness dder and subcutis scattering amplitude A were forced to obtain the same values for both examples. The relations between the experimental DRS and best fitting DA predictions above 600 nm now resemble more closely those predicted by our MC simulations (Fig. 6) as opposed to Fig. 9.
Fig. 10. The same experimental DRS as in Fig. 9 (orange solid lines) fitted with three-layer DA solutions (dashed navy lines). Fitting was performed with the modified spectral sampling (see the text) and the parameters dder and A linked to the same values for both cases.
As is evident from Table 6, the described augmentation of the fitting approach has dramatically reduced the estimates of relative error (r) for the linked parameters dder and A, and to a smaller extent also for most of the remaining ones. This demonstrates how simultaneous analysis of two related experimental DRS can alleviate the problem of non-uniqueness in our inverse analysis, arising from intrinsic correlations between certain model parameters.
Table 6. Skin parameter values as assessed from the same experimental DRS as in Table 5 but using the augmented fitting approach (see text for details).
This effect is illustrated in Fig. 11, which presents the dependence of residual norm (ɛ) on the values of model parameters dder and A, with the remaining parameters fixed at their optimal values. When analyzing two experimental DRS independently, the presented “landscapes” feature extended banana-shaped minima, characteristic for strongly correlated parameters (Figs. 11(a) and (b)). Upon linking A and dder (and simultaneously modifying the spectral sampling), however, the minimum of the common residual norm is much better localized and the “valley walls” considerably steeper (Fig. 11(c)), which results in more robust assessment of these parameters.
Fig. 11. Residual norm (ɛ) as a function of the dermal thickness and subcutis scattering amplitude (A): (a) and (b) - when fitting the DRS in Figs. 9(a) and (b) independently; (c) when fitting both datasets simultaneously with linked values for dder and A (see text for details).
6. Monitoring of changes in human skin in vivo
The results in Table 6 clearly reveal the prominent increase in epidermal melanin content upon extensive sun exposure (from m = 1.6% to 3.3%). A smaller increase of the dermal blood content (bder) and oxygen saturation (S) are also indicated, indicative of the related erythema (skin redness). This matches earlier reports, which demonstrated elevated cutaneous blood contents for 3 weeks after UV irradiation [13].
In order to test more thoroughly the reliability of presented methodology and its suitability for monitoring of changes in human skin in vivo, we submit it in the following to two tests, involving temporary obstruction of cutaneous blood circulation by application of a pressurized arm cuff and a series of measurements performed over the course of several months.
6.1 Obstruction of cutaneous blood circulation
Figure 12 presents a comparison of the most relevant skin parameter values assessed from DRS measurements on the dorsal side of the forearm in three subjects (one female and two males) before and after application of a blood-pressure cuff for 5 minutes. A substantial increase of the dermal blood content (on average by a factor of 2) and a dramatic drop of oxygen saturation level after the application of the pressure cuff are evident in all three subjects (Figs. 12(a) and (b), respectively). This matches the expected effect of blood accumulation in the skin, as pressurization of the cuff initially obstructs blood flow through the veins, and only later also through the arteries [45]. The epidermal blood content consistently remains much smaller in comparison (bepi < 0.1%, not plotted).
Fig. 12. (a) Dermal blood contents, (b) oxygen saturation levels, and (c) melanin contents as obtained by fitting three-layer DA solutions to DRS acquired before and after the application of the blood–pressure cuff in three healthy subjects (denoted A, B and C).
It is also important to notice that only minimal changes of the epidermal melanin content are indicated in all three test subjects (Fig. 12(c)), despite the fact that the DRS acquired before and after blood obstruction were analyzed independently. Moreover, the subject with considerably lower m (marked with C) is indeed a male with very light skin (Fitzpatrick type I) and red hair, while the volunteers A and B have skin type II and brown hair. This test thus demonstrates the sensitivity of the discussed approach for detection of both inter-personal and temporal variations of skin composition, as well as robustness in terms of minimal correlations (cross-talk) between the assessed parameter values.
6.2 Seasonal variations of skin composition
In Fig. 13 we present our analysis of DRS acquired from a dorsal side of the forearm in a healthy male subject over a time span of several months, from late June through early December. All DRS from the same test site were analyzed together, as a single optimization problem with linked values for dder and A (see Sect. 5.2.2). In the example presented below, involving 16 measurement sessions, this leads to simultaneous optimization of 64 independent parameters.
Fig. 13. Seasonal changes of (a) melanin content, and (b) dermal blood content (red solid dots) and oxygen saturation (blue open dots) as assessed by our analysis of DRS acquired from dorsal side of the forearm in a healthy volunteer (LV).
In the graph showing the obtained seasonal variations of the epidermal melanin (Fig. 13a), dashed vertical lines indicate three time points of particular interest, namely the measurements performed after extensive sun exposure (during a seaside vacation, hiking trip, etc.). A pronounced jump of the melanin content is evident at all such time points, followed by a gradual decline. The decline following the time point #3 (Aug. 31) is monotonic and smooth, which indicates that randomness in these results is minimal.
The dermal blood content (Fig. 13(b), red solid dots) features less prominent, yet still observable peaks at the same time points. This is consistent with skin erythema during the recovery from sun damage [13], seen also in Table 6. Note that all DRS data were acquired on Fridays, while sun exposure occurred primarily over (or ended during) the previous weekends. This was a deliberate decision, aiming at limiting the influence of acute erythema on our analysis.
The oxygen saturation values (Fig. 13(b), blue open dots) exhibit considerable variations, but no evident correlation with the behavior of m or bder. This is not uncommon in analyses of skin using DRS [14,15] and could be related, e.g., to uncontrolled physiological processes or variations in pressure applied to the skin during the measurements.
The epidermal blood content values, while fitted independently at each time point, were negligible throughout the season (bepi ≤ 0.1%).
The assessed dermal thickness and subcutis scattering amplitude A (representing the average values for the entire time series) are dder = 0.72 ± 0.02 mm and A = 2.00 ± 0.07, both within the respective ranges of plausible values.
The analyses for the volar side of the forearm and the second volunteer led to the same observations and correlations as discussed above, with inevitable differences in the specific timelines and parameter values [32]. In contrast, the measurements on the forehead could not be fitted well with the present skin model [32]. We tentatively attribute this effect to a too thin subcutis layer, thus requiring explicit accounting for optical properties of the underlying cartilage and/or bone.
7. Conclusion
Analytical solutions derived within the DA enable simple analysis of DRS acquired from human skin in vivo using an IS. A very good match between the measured and best fitting model DRS is obtained by using the presented three-layer skin model (accounting for the epidermis, dermis and subcutis), unlike the two-layer model which assumes a semi-infinite dermis. However, the spectral range included in the analysis must be carefully selected, e.g., by considering the diameter of the IS sample opening.
Mathematical simplicity of the analytical DA solutions allows analysis of multiple DRS acquired form the same test site in parallel, with selected parameters linked to the same value for the entire data set. This reduces the ambiguity in the assessed values due to intrinsic correlations between the model parameters, thus enabling robust monitoring of physiological changes in skin upon various stimuli (e.g., obstruction of cutaneous blood circulation and sun exposure). Applying the same approach in direct inverse MC analysis is also possible, but leads to a considerably higher computational load due to increased dimensionality of the resulting optimization problem confounded with the stochastic noise inherent in this numerical technique [46].
After accounting for the Fresnel reflection at the skin surface, the diffuse reflectance values predicted using the DA solutions (with the δ-Edington source function) systematically overshoot the numerical results (Monte Carlo), typically by 1–2 percentage points across the visible spectral range. Nevertheless, our simulation study shows that applying the former to analyze the latter leads to surprisingly large systematic errors, most notably an overestimation of the melanin and dermal blood contents by 5–15% and 10–15%, respectively, depending on the remaining skin properties. Such artifacts most likely occur also in analysis of experimental DRS using the same approach, and should be duly considered in interpretation of the assessed values.
Funding
Javna Agencija za Raziskovalno Dejavnost RS (ARRS) (J2-7211, P1-0192, P2-0232, PR-04360).
Acknowledgment
The methodology and some provisional results were presented at Biophotonics Riga 2013 [24].
References
1. J. R. Mourant, I. J. Bigio, J. Boyer, R. L. Conn, T. Johnson, and T. Shimada, “Spectroscopic diagnosis of bladder cancer with elastic light scattering,” Lasers Surg. Med. 17(4), 350–357 (1995). [CrossRef]
2. G. Zonios, L. T. Perelman, V. Backman, R. Manoharan, M. Fitzmaurice, J. Van Dam, and M. S. Feld, “Diffuse reflectance spectroscopy of human adenomatous colon polyps in vivo,” Appl. Opt. 38(31), 6628–6637 (1999). [CrossRef]
3. F. Bevilacqua, D. Piguet, P. Marquet, J. D. Gross, B. J. Tromberg, and C. Depeursinge, “In vivo local determination of tissue optical properties: applications to human brain,” Appl. Opt. 38(22), 4939–4950 (1999). [CrossRef]
4. I. J. Bigio, S. G. Bown, G. Briggs, C. Kelley, S. Lakhani, D. Pickard, P. M. Ripley, I. G. Rose, and C. Saunders, “Diagnosis of breast cancer using elastic-scattering spectroscopy: preliminary clinical results,” J. Biomed. Opt. 5(2), 221–228 (2000). [CrossRef]
5. G. M. Palmer, C. Zhu, T. M. Breslin, F. Xu, K. W. Gilchrist, and N. Ramanujam, “Monte Carlo-based inverse model for calculating tissue optical properties. Part II: Application to breast cancer diagnosis,” Appl. Opt. 45(5), 1072–1078 (2006). [CrossRef]
6. M. B. Wallace, L. T. Perelman, V. Backman, J. M. Crawford, M. Fitzmaurice, M. Seiler, K. Badizadegan, S. J. Shields, I. Itzkan, R. R. Dasari, et al., “Endoscopic detection of dysplasia in patients with Barrett’s esophagus using light-scattering spectroscopy,” Gastroenterology 119(3), 677–682 (2000). [CrossRef]
7. E. Salomatina, B. Jiang, J. Novak, and A. N. Yaroslavsky, “Optical properties of normal and cancerous human skin in the visible and near-infrared spectral range,” J. Biomed. Opt. 11(6), 064026 (2006). [CrossRef]
8. G. Zonios, A. Dimou, I. Bassukas, D. Galaris, A. Tsolakidis, and E. Kaxiras, “Melanin absorption spectroscopy: new method for noninvasive skin investigation and melanoma detection,” J. Biomed. Opt. 13(1), 014017 (2008). [CrossRef]
9. T. M. Bydlon, R. Nachabé, N. Ramanujam, H. J. C. M. Sterenborg, and B. H. W. Hendriks, “Chromophore based analyses of steady-state diffuse reflectance spectroscopy: current status and perspectives for clinical adoption,” J. Biophoton. 8(1–2), 9–24 (2015). [CrossRef]
10. N. Kollias, I. Seo, and P. R. Bargo, “Interpreting diffuse reflectance for in vivo skin reactions in terms of chromophores,” J. Biophoton. 3(1-2), 15–24 (2009). [CrossRef]
11. D. Yudovsky and L. Pilon, “Retrieving skin properties from in vivo spectral reflectance measurements,” J. Biophoton. 4(5), 305–314 (2011). [CrossRef]
12. I. V. Ermakov and W. Gellermann, “Dermal carotenoid measurements via pressure mediated reflection spectroscopy,” J. Biophoton. 5(7), 559–570 (2012). [CrossRef]
13. G. N. Stamatas and N. Kollias, “Blood stasis contributions to the perception of skin pigmentation,” J. Biomed. Opt. 9(2), 315–322 (2004). [CrossRef]
14. L. L. Randeberg, O. A. Haugen, R. Haaverstad, and L. O. Svaasand, “A novel approach to age determination of traumatic injuries by reflectance spectroscopy,” Lasers Surg. Med. 38(4), 277–289 (2006). [CrossRef]
15. S.-H. Tseng, C.-K. Hsu, J. Y.-Y. Lee, S.-Y. Tzeng, W.-R. Chen, and Y.-K. Liaw, “Noninvasive evaluation of collagen and hemoglobin contents and scattering property of in vivo keloid scars and normal skin using diffuse reflectance spectroscopy: pilot study,” J. Biomed. Opt. 17(7), 0770051 (2012). [CrossRef]
16. A. Ishimaru, Wave Propagation and Scattering in Random Media (Academic, 1978), Vol. 2.
17. K. Furutsu, “Diffusion equation derived from space-time transport equation,” J. Opt. Soc. Am. 70(4), 360–366 (1980). [CrossRef]
18. L. Wang, S. L. Jacques, and L. Zheng, “MCML—Monte Carlo modeling of light transport in multi-layered tissues,” Comput. Meth. Prog. Bio. 47(2), 131–146 (1995). [CrossRef]
19. E. Alerstam, T. Svensson, and S. Andersson-Engels, “Parallel computing with graphics processing units for high-speed Monte Carlo simulation of photon migration,” J. Biomed. Opt. 13(6), 060504 (2008). [CrossRef]
20. R. Hennessy, S. L. Lim, M. K. Markey, and J. W. Tunnell, “Monte Carlo lookup table-based inverse model for extracting optical properties from tissue-simulating phantoms using diffuse reflectance spectroscopy,” J. Biomed. Opt. 18(3), 037003 (2013). [CrossRef]
21. M. Sharma, R. Hennessy, M. K. Markey, and J. W. Tunnell, “Verification of a two-layer inverse Monte Carlo absorption model using multiple source-detector separation diffuse reflectance spectroscopy,” Biomed. Opt. Express 5(1), 40 (2014). [CrossRef]
22. Q. Wang, D. Le, J. Ramella-Roman, and J. Pfefer, “Broadband ultraviolet-visible optical property measurement in layered turbid media,” Biomed. Opt. Express 3(6), 1226 (2012). [CrossRef]
23. S.-Y. Tsui, C.-Y. Wang, T.-H. Huang, and K.-B. Sung, “Modelling spatially-resolved diffuse reflectance spectra of a multi-layered skin model by artificial neural networks trained with Monte Carlo simulations,” Biomed. Opt. Express 9(4), 1531 (2018). [CrossRef]
24. P. Naglič, L. Vidovič, M. Milanič, L. L. Randeberg, and B. Majaron, “Applicability of diffusion approximation in analysis of diffuse reflectance spectra from healthy human skin,” Proc. SPIE 9032, 90320N (2013). [CrossRef]
25. M. Milanic and B. Majaron, “Three-dimensional Monte Carlo model of pulsed-laser treatment of cutaneous vascular lesions,” J. Biomed. Opt. 16(12), 128002 (2011). [CrossRef]
26. S. L. Jacques, “Skin Optics,” http://omlc.org/news/jan98/skinoptics.html, January 1998.
27. W. G. Zijlstra, A. Buursma, and O. W. van Assendelft, Visible and Near Infrared Absorption Spectra of Human and Animal Haemoglobin: Determination and Application (VSP, 2000).
28. C. R. Simpson, M. Kohl, M. Essenpreis, and M. Cope, “Near-infrared optical properties of ex vivo human skin and subcutaneous tissues measured using the Monte Carlo inversion technique,” Phys. Med. Biol. 43(9), 2465–2478 (1998). [CrossRef]
29. S. L. Jacques, “Optical assessment of cutaneous blood volume depends on the vessel size distribution: a computer simulation study,” J. Biophoton. 3(1-2), 75–81 (2009). [CrossRef]
30. D. Yudovsky and L. Pilon, “Rapid and accurate estimation of blood saturation, melanin content, and epidermis thickness from spectral diffuse reflectance,” Appl. Opt. 49(10), 1707–1719 (2010). [CrossRef]
31. A. N. Bashkatov, E. A. Genina, and V. V. Tuchin, “Optical properties of skin, subcutaneous, and muscle tissues: A review,” J. Innov. Opt. Heal. Sci. 04(01), 9–38 (2011). [CrossRef]
32. P. Naglič, “Determination of structure of human skin using diffuse reflectance spectroscopy,” Master Thesis, University of Ljubljana (2013).
33. A. N. Bashkatov, E. A. Genina, V. I. Kochubey, and V. V. Tuchin, “Optical properties of the subcutaneous adipose tissue in the spectral range 400–2500 nm,” Opt. Spectrosc. 99(5), 836–842 (2005). [CrossRef]
34. A. N. Bashkatov, E. A. Genina, V. I. Kochubey, and V. V. Tuchin, “Optical properties of human skin, subcutaneous and mucous tissues in the wavelength range from 400 to 2000 nm,” J. Phys. D: Appl. Phys. 38(15), 2543–2555 (2005). [CrossRef]
35. R. A. J. Groenhuis, H. A. Ferwerda, and J. J. Ten Bosch, “Scattering and absorption of turbid materials determined from reflection measurements. 1: Theory,” Appl. Opt. 22(16), 2456–2462 (1983). [CrossRef]
36. R. C. Haskell, L. O. Svaasand, T.-T. Tsay, T.-C. Feng, M. S. McAdams, and B. J. Tromberg, “Boundary conditions for the diffusion equation in radiative transfer,” J. Opt. Soc. Am. B 11(10), 2727–2741 (1994). [CrossRef]
37. L. O. Svaasand, L. T. Norvang, E. J. Fiskerstrand, E. K. S. Stopps, M. W. Berns, and J. S. Nelson, “Tissue parameters determining the visual appearance of normal skin and port-wine stains,” Laser. Med. Sci. 10(1), 55–65 (1995). [CrossRef]
38. T. Spott and L. O. Svaasand, “Collimated light sources in the diffusion approximation,” Appl. Opt. 39(34), 6453–6465 (2000). [CrossRef]
39. L. L. Randeberg, A. Winnem, R. Haaverstad, and L. O. Svaasand, “Performance of diffusion theory vs Monte Carlo methods,” Proc. SPIE 5862, 58620O (2005). [CrossRef]
40. T. Spott, “Characterization of layered tissue structures with diffusely propagating photon-density waves,” Doctoral Thesis, Norwegian University of Science and Technology (1999).
41. L. Vidovič and B. Majaron, “Elimination of single-beam substitution error in diffuse reflectance measurements using an integrating sphere,” J. Biomed. Opt 19(2), 027006 (2014). [CrossRef]
42. Y. Lee and K. Hwang, “Skin thickness of Korean adults,” Surg. Radiol. Anat. 24(3-4), 183–189 (2002). [CrossRef]
43. G. Zonios and A. Dimou, “Modeling diffuse reflectance from semi-infinite turbid media: application to the study of skin optical properties,” Opt. Express 14(19), 8661–8674 (2006). [CrossRef]
44. C. Donner and H. W. Jensen, “A spectral BSSRDF for shading human skin,” in Proceedings of the 17th Eurographics Conference on Rendering Techniques (Eurographics Association, 2006), pp. 409–417.
45. T. Strömberg, F. Sjöberg, and S. Bergstrand, “Temporal and spatiotemporal variability in comprehensive forearm skin microcirculation assessment during occlusion protocols,” Microvasc. Res. 113, 50–55 (2017). [CrossRef]
46. N. Verdel, M. Milanič, and B. Majaron, “Physiological and structural characterization of human skin in vivo using combined photothermal radiometry and diffuse reflectance spectroscopy,” Biomed. Opt. Express 10(2), 944–960 (2019). [CrossRef]
