Proposal for a Skin Layer-Wise Decomposition Model of Spatially-Resolved Diffuse Reflectance Spectra Based on Maximum Depth Photon Distributions: A Numerical Study

Abstract

In the context of cutaneous carcinoma diagnosis based on in vivo optical biopsy, Diffuse Reflectance (DR) spectra, acquired using a Spatially Resolved (SR) sensor configuration, can be analyzed to distinguish healthy from pathological tissues. The present contribution aims at studying the depth distribution of SR-DR-detected photons in skin from the perspective of analyzing how these photons contribute to acquired spectra carrying local physiological and morphological information. Simulations based on modified Cuda Monte Carlo Modeling of Light transport were performed on a five-layer human skin optical model with epidermal thickness, phototype and dermal blood content as variable parameters using (i) wavelength-resolved scattering and absorption properties and (ii) the geometrical configuration of a multi-optical fiber probe implemented on an SR-DR spectroscopic device currently used in clinics. Through histograms of the maximum probed depth and their exploitation, we provide numerical evidence linking the characteristic penetration depth of the detected photons to their wavelengths and four source–sensor distances, which made it possible to propose a decomposition of the DR signals related to skin layer contributions.

1. Introduction

Diffuse Reflectance Spectroscopy (DRS) applied to biological tissues is a widely studied [1,2,3,4], non-destructive optical characterization technique. It consists of measuring back-reflected intensity spectra carrying information from light–tissue interactions, from which the optical and structural properties of the probed medium can then be extracted and analyzed. By illuminating the tissue with a broadband light, photons propagating in the medium undergo wavelength-dependent absorption and scattering events, before some of them may move back towards the surface and potentially be captured by a spectrometric detection system. An analysis of the amplitude and shape features of the latter intensity spectrum thus gives access to information of interest such as absorption bands and elastic scattering properties related to tissue-specific chromophores and structures. Its non-invasive nature, its sensitivity to metabolic and structural changes in tissues and its clinical applicability are very well suited to in vivo biomedical diagnostics.

Spatially Resolved DRS (SR-DRS) refers to fibered systems of DRS featuring several Detection Fibers (DF) located at different distances from the Source Fiber (SF) allowing for depth-resolved acquisitions. Indeed, photons detected at a short/large source–detector (SD) distance will be associated with shallower/deeper trajectories into the biological tissues under evaluation. The literature on this subject points to numerous applications of SR-DRS such as the measurement of blood glucose levels [5,6], hemoglobin oxygen saturation [7] and for pre-cancer and cancer diagnosis: lung [8], cervix [9], colon [10], bladder [11], brain [12], breast [13], oesophagus [14] or pancreas [15]. There are also a large number of applications of SR-DRS for the in vivo characterization of cutaneous pathological states, including melanoma [16,17,18] and carcinoma [19,20], and more generally for the characterization of skin layer optical properties [21,22].

In this context of skin optical biopsy, point spectroscopy devices exploiting light–tissue interactions have been designed and used in several clinical studies [11,23,24]. The estimation of diagnostically relevant optical parameters from in vivo SR-DR spectra, and especially their evolution around the lesion, can be linked to the tissue pathological state. As the metabolic modifications associated with carcinogenesis are local, the quantitative estimation of the penetration depth of the photons detected using point spectroscopy is at stake, as reported by Breslin et al. [25], who studied the sensitivity of DR spectra to elastosis (accumulation of elastotic material in the papillary dermis located underneath epidermis). Cancer growth is also related to the enzymatic degradation of collagen fibers in the dermis leading to collagen fluorescence decrease [26] and to the decrease of Nicotinamide Adenine Dinucleotide Hydrogen (NADH) fluorescence in the epidermis [27]. Finally, morphological changes may be significant for part of skin carcinomas such as ulceration and hyperkeratosis, which are respectively characterized by a thinning or skin thickening (∼×10) of the stratum corneum superficial layer.

The need to explain DR signals acquired on the surface, resulting from light–tissue interactions in all skin layers, justifies the use of numerical simulation as a tool to interpret the acquired signals. The literature reports work by several teams using photon transport simulation in a multi-layer skin model to extract features that describe the propagation of the acquired photons [28,29]. First, in a very visual way, absorption density maps of detected photons in a three-layer skin model [30] provide information on a photon’s trajectories. In a more quantitative manner, teams have presented their results on the penetration depth of photons. In a two-layer model, histograms of probed depths have been created at 450 and 650 nm to extract the average penetration depth [31], while the interrogation depth (mean depth of all detected photon scattering events) [32] or measurement depth (median of the cumulative absorption on fluence absorption maps) [30] can also be extracted from simulation to characterize the photon’s penetration depth. Other quantities of interest are exploited to know more about the subcutaneous behavior of backscattered photons, such as the path length distribution with increasing SD distance [33,34]. This quantity was even declined in the mean partial path length [35,36], which consists of dividing the total path length into layer path contributions.

Although these approaches have made it possible to study the behavior of photons (path, depth) in different multi-layer skin models, none of them, to our knowledge, has led to the analytical formulation of an explicit decomposition of the “global” DR spectra, as a function of the SD distances and of the modeled skin layers. This is an essential issue to (i) directly interpret the measured spectra, especially in the case of clinical application and patients’ inter-individual variability in terms of epidermal thicknesses, skin type (healthy, carcinoma), phototype and blood dermal content and to (ii) indirectly estimate the optical properties of the different layers in an efficient way. In order to address this gap, the present contribution deals with the acquisition of simulated DR spectra for a multi-layer medium considering (i) the absorption and scattering of the five main endogenous chromophores (ii) on the near UV-visible range and (iii) for a geometrical model of the probe adapted to the previous wavelengths, namely SD distances smaller than a millimeter. Histograms of the maximal probed depth for photons detected at several SD distances were created and analyzed to quantify light penetration in the media. The study was led for skin media with varying parameters: epidermal thickness, phototype, dermal blood concentration and health status of the tissue (ulcerated and hyperkeratosed carcinomas). Finally, a DR signal decomposition according to the deepest layer probed is offered.

2. Materials and Methods

Section 2.1 describes the geometrical and spectral configuration of the probe, Section 2.2 provides an explanation of the numerical skin model developed, and Section 2.3 refers to the simulation description.

2.1. SR-DRS Probe—Geometrical and Spectral Features Implemented for Simulations

The geometrical and optical features of source and detection fibers implemented in our SR-DRS simulations were based on the technical characteristics of an existing device developed by our team [37], whose purpose is to assist a surgeon in establishing surgical margins during skin carcinoma resection [38]. Work with a similar goal has already been presented in the literature using DR [39] (elastic scattering) or Raman [40] (inelastic scattering) spectra. The multiple fiber optical probe of the clinical device consists of (i) one central 600 $\mathsf{\mu }$m-core diameter Source Fiber (SF) and numerical aperture $N{A}_{SF}=0.22$ delivering broadband 365–765 nm white light illumination to the tissue in contact and (ii) 24,200 $\mathsf{\mu }$m core diameter detection fibers (DF) of $N{A}_{DF}=0.22$ each, equally distributed on 4 concentric circles corresponding to the 4 SD center-to-center distances $D_n$ ($n\in \{1,...,4\}$) from 400 to 1000 $\mathsf{\mu }$m by 200 $\mathsf{\mu }$m steps. The photons reaching each of the six collection fibers of the $D_n$ distance are merged to constitute the diffuse reflectance signal $D{R}_{D_n}\left(\lambda \right)$. To maximize the signal to noise ratio in the simulation, cylindrical symmetry was exploited using a ring detection surface before multiplying the quantity by a factor that ensures the geometrical correspondence with the clinical device. Complementary simulations were conducted to quantify the error related to this geometric approximation (detailed results not shown). In the “worst” case—i.e., for small SD distances—the relative error was less than 0.5%. A schematic representation of the acquisition geometry is shown in Figure 1, while the numerical values of the different parameters used are summarized in

Table 1. Values of the geometrical and spectral source–detector features for the SR-DRS simulation. $SF$: source fiber; $DF$: detection fiber; $D_n$: SD distance; $r_{SF/DF}$: fiber radius.

$[{\mathit{\lambda }}_{\mathbf{min}},{\mathit{\lambda }}_{\mathbf{max}}]$ (nm)	${\mathit{D}}_{\mathit{n}}$ ($\mathsf{\mu }$m)	${\mathit{NA}}_{\mathit{SF}/\mathit{DF}}$	${\mathit{r}}_{\mathbf{SF}}$ ($\mathsf{\mu }$m)	${\mathit{r}}_{\mathbf{DF}}$ ($\mathsf{\mu }$m)
$[365,765]$	400, 600, 800, 1000	0.22	300	100

. The optical probe geometrical configuration considered in our simulation does not take into account the epoxy resin and the polished stainless steel in which optical fibers (whose core and cladding are made of silica) are embedded in the clinical device. As mentioned by Zelinskyi et al., this may limit the accuracy of the results obtained [41].

2.2. Skin Five-Layered Model

2.2.1. Layer Thicknesses and Skin Types Modeled

The numerical model mimicking the skin was designed following a five planar multi-layered medium, representing (from superficial to deepest layers) the stratum corneum (1-SC), the living epidermis (2-LE), the papillary dermis (3-PD), the reticular dermis (4-RD) and the subcutaneous fat (5-SF). The total thickness was fixed to $z^{tot}=6120\mathsf{\mu }$m, while the total radius of the simulation medium was set to $r^{tot}=3000\mathsf{\mu }$m. The results presented hereafter are valid under the assumption that the optical properties of each of the skin layers are homogeneous in the cylinder thus described (cf. Figure 1). Considering the order of magnitude of skin optical properties and the position of the the farthest collecting fiber ($D_4=1000\mathsf{\mu }$m), edge effects (i.e., photons that would be collected by DF after reflection on the side or bottom border of the medium) were negligible. This was verified numerically with absorption maps of detected photons: none of them reached the bottom (in the worst case—i.e., a wavelength in near IR, great SD distance and fair phototype—only 1% of detected photons were scattered in the SF), while less than 0.1% were collected after having reached the side border.

Consequently, five combinations of SC and LE thicknesses were defined, among which three featured a thin, medium and thick healthy epidermis and two represented some specific values encountered in carcinoma. SC and LE thicknesses may significantly vary according to age, gender, and body location [42,43,44] and also for different skin physiological states. For instance, the stratum corneum can either completely disappear or thicken 10-fold [45] in the case of ulceration or hyperkeratosis, respectively. Indeed, other epidermal layers are also affected by cell proliferation due to skin carcinoma. However, for the sake of comparison with healthy conditions, only SC thickness was varied to model these two sub-types of skin carcinomas; i.e., in the case of ulceration, the absence of SC was modeled through the removal this superficial layer. As shown in Figure 1, the thickness values of middle layers PD and RD were fixed (200 and $1800\mathsf{\mu }$m, respectively), while those of the top layers (SC and LE) were varied to simulate different tissue conditions. The thickness of the underlying SF layer was consequently adjusted to keep the total thickness $z^{tot}$ of our simulation model constant. We considered the range of values of SC and LE thicknesses as well as their approximate ratio (∼1/4) reported for sun-exposed skin sites only (head including face, scalp, ear, cheek and shoulders, back, arms, back of hand and legs), which are 10–30 $\mathsf{\mu }$m and 60–120 $\mathsf{\mu }$m for SC and LE layers, respectively. The numerical values associated with the elaboration of these three healthy and two other pathological epithelium (carcinoma) simulation media configurations are summarized in

Table 2. Values of the SC (stratum corneum) and LE (living epidermis) layer thicknesses implemented in our skin model to simulate five different healthy skin and carcinoma conditions.

Skin Categories		SC Thickness ${\mathit{z}}^1$ ($\mathsf{\mu }$m)	LE Thickness ${\mathit{z}}^2$ ($\mathsf{\mu }$m)
Healthy	Smallest thicknesses	10	60
	Mean thicknesses	20	90
	Highest thicknesses	30	120
Carcinoma	Ulcerated	0	90
	Hyperkeratosed	300	90

.

2.2.2. Skin Optical Properties

The optical properties defined for each skin layer in terms of absorption and scattering are synthesized in the following paragraphs. Complementary details can also be found in [46].

Absorption Properties

Absorption coefficients ${\mu }_a^k\left(\lambda \right)$ were defined for each layer $k\in \{1,...,5\}$ of the numerical skin model taking into account the contribution of the five main endogenous absorbers i of human skin, which are blood (decomposed in oxygenated hemoglobin $Hb{O}_2$, deoxygenated hemoglobin $Hb$, and bilirubin contributions), water, melanin (decomposed in eumelanin and pheomelanin), lipids and $\beta$-caroten, respectively, as well as their volume fraction $f_{v,i}^k$, in the mixture. A residual absorption baseline was also added according to [47] (cf.${\mu }_{a,base}$ in Equation (1)). Thus, the total absorption coefficient of a skin layer was defined as the sum of the latter baseline and the contributions ${\mu }_{a,i}^k$ of each of the five chromophores, weighted by their respective volume fraction, in accordance with [48]:

$$\begin{array}{cc}\hfill {\mu }_a^k\left(\lambda \right)=& {\mu }_{a,base}\left(\lambda \right)\hfill \\ \hfill & +f_{v,blood}^k[\gamma {\mu }_{a,Hb{O}_2}^k\left(\lambda \right)+(1-\gamma ){\mu }_{a,Hb}^k\left(\lambda \right)+{\mu }_{a,bilirubin}^k\left(\lambda \right)]\hfill \\ \hfill & +f_{v,water}^k{\mu }_{a,water}^k\left(\lambda \right)\hfill \\ \hfill & +f_{v,mel}^k[{\mu }_{a,eum}^k\left(\lambda \right)+{\mu }_{a,phe}^k\left(\lambda \right)]\hfill \\ \hfill & +f_{v,lipids}^k{\mu }_{a,lipids}^k\left(\lambda \right)\hfill \\ \hfill & +f_{v,\beta -caroten}^k{\mu }_{a,\beta -caroten}^k\left(\lambda \right)\hfill \end{array}$$

(1)

where ${\mu }_{a,i}^k$ are all expressed in cm${}^{-1}$, ${\mu }_{a,base}\left(\lambda \right)=7.84\times {10}^8\times {\lambda }^{-3225}$ ($\lambda$ in nm), and $\gamma$ is the ratio between $Hb{O}_2$ and $Hb$, set to $\gamma =0.9$ in our simulations. Lipids and water absorption spectra as well as extinction coefficients of the other compounds were directly implemented in our model from [49]. It is noticeable that, for $\beta$-caroten and blood (namely both hemoglobins and bilirubin), molar extinction coefficient ${ϵ}_i^{mol}\left(\lambda \right)$ (in cm${}^{-1}$/mol·L${}^{-1}$) was given, while for melanin, the mass extinction coefficients ${ϵ}_i^{mass}\left(\lambda \right)$ (in cm${}^{-1}$/g·L${}^{-1}$) of eu- and pheo-melanin were provided. These are linked to absorption coefficients according to the following formulas:

$$\begin{array}{cc}\hfill {\mu }_{a,i}^k\left(\lambda \right)=\frac{{ϵ}_i^{mol}\left(\lambda \right)}{M_i}{c}_i^k& {\mu }_{a,i}^k\left(\lambda \right)={ϵ}_i^{mass}\left(\lambda \right)c_i^k\hfill \end{array}$$

(2)

where $c_i^k$ is the mass concentration (in g/L) of chromophore i in the mixing model and $M_i$ the respective molar weight with the following values implemented in our simulations:

$$M_{\beta -Caroten}=587\mathrm{g}·{\mathrm{mol}}^{-1}{M}_{bilirubin}=585\mathrm{g}·{\mathrm{mol}}^{-1}{M}_{Hb{O}_2,Hb}=66,500\mathrm{g}·{\mathrm{mol}}^{-1}$$

In order to take into account skin pigmentation-related absorption variations, four different melanosome volume fraction values in the LE layer were used to simulate very light skin (phototype I) to dark skin (phototype IV) types. Considering our targeted clinical application of interest, it was chosen to firstly consider a melanosome volume fraction $f_{v,mel}^2=0.01$, corresponding to phototype I, which is much more prone to the development of skin carcinomas. In Section 3.2.3, darker skins are considered. The correspondence between phototypes and volume fractions of melanosome, adapted from [50], is summarized in

Table 3. Values of the volume fraction of melanosome in the living epidermis used for simulating skin phototypes I to IV.

Skin Type	Very Fair	Fair	Moderately Fair	Dark Skin
Phototype	I	II	III	IV
Melanosome volume fraction $f_{v,mel}^2$	0.01	0.04	0.08	0.11

.

In our skin model, this parameter only influences the LE absorption coefficient, but its impact is significant. The absorption spectra resulting from the different calculations are finally shown in Figure 2. The top right subplot focuses on the absorption coefficient of LE.

Similarly, the impact of the amount of blood in the dermal layers, modified by the volume fraction of blood $f_{v,blood}^3$ and $f_{v,blood}^4$, was considered as an additional parameter of interest [51]. Simulations were then conducted for three pairs of blood volume fractions to represent three media; respectively, weakly, moderately and strongly vascularized. The numerical values of these configurations are shown in

Table 4. Values of the volume fraction of blood in the papillary dermis (PD) and reticular dermis (RD) used for simulating a dermis with more or less blood content.

Blood Content	Weak	Mean	Strong
PD blood volume fraction $f_{v,blood}^3$	0.0144	0.024	0.0336
RD blood volume fraction $f_{v,blood}^4$	0.012	0.02	0.028

. Here again, the couple of mean values—i.e., $f_{v,blood}^3=0.024$ and $f_{v,blood}^4=0.02$—was used for all simulations, except in Section 3.2.4, which focuses on this blood amount effect.

This variation influences mainly the absorption coefficient of the two dermal layers (cf. Equation (1)) but slightly influences the scattering coefficient (cf. Equation (3)). The optical parameters of these three couples are then isolated in the middle-right (${\mu }_a^3$ and ${\mu }_a^4$) and bottom-right (${\mu }_s^3$ and ${\mu }_s^4$) subplot of Figure 2.

The values of the volume fraction and mass concentration of each chromophore applied in every layer of our numerical model are summarized in

Table 5. Values of the geometrical and optical parameters of our five layer skin model. SC, LE and SF thickness ranges, $f_{v,mel}^k$ melanosome fraction range in LE and $f_{v,blood}^k$ blood fraction in dermis layers refer to Table 2, Table 3 and Table 4, respectively.

Name (Unit)	Symbol	SC	LE	PD	RD	SF
Absorption parameters
Volume fraction (dimensionless)	$f_{v,water}^k$	0.35	0.6	0.5	0.7	0.05
	$f_{v,lipids}^k$	0.2	0.15	0.173	0.173	0.753
	$f_{v,\beta -caroten}^k$	0.45	0.24	0.303	0.125	0.141
	$f_{v,mel}^k$	0	See Table 3	0	0	0
	$f_{v,blood}^k$	0	0	See Table 4	See Table 4	0.005
Concentration (g/L)	$c_{\beta -caroten}^k$	$2.1\times {10}^{-4}$	$2.1\times {10}^{-4}$	$7.0\times {10}^{-5}$	$7.0\times {10}^{-5}$	0
	$c_{eum}^k$	0	80	0	0	0
	$c_{phe}^k$	0	12	0	0	0
	$c_{bilirubin}^k$	0	0	0.05	0.05	0.05
	$c_{Hb{O}_2/Hb}^k$	0	0	150	150	150
Scattering parameters
Volume fraction	$f_{v,blood}^k$	0	0	See Table 4	See Table 4	0.05
Vessel diam. (cm)	$d_{vessels}^k$	/	/	$6.0\times {10}^{-4}$	$1.5\times {10}^{-3}$	$7.5\times {10}^{-3}$
Scat. coeff. (cm${}^{-1}$)	${\mu }_{s,577}^k$	400	300	120	120	130
Other optical properties and medium geometry
Layer thick. ($\mathsf{\mu }$m)	$z^k$	See Table 2	See Table 2	200	1800	See Table 2
Anisotropy factor	$g^k$	0.92	0.75	0.72	0.72	0.9
Refractive index	$n^k$	1.55	1.44	1.39	1.38	1.34

. They were chosen following (i) a literature review [47,49,52,53] and (ii) comparisons between simulated spectra and experimental data acquired in clinics [38] (see Figure 3) using our SR-DRS device [37].

Scattering Properties

The values of skin scattering coefficients ${\mu }_s^k\left(\lambda \right)$, expressed in cm${}^{-1}$ in the various layers of our model, were obtained with Equation (3) given from [53] featuring two contributions: a blood scattering component and a bloodless “base” scattering component, such that

$${\mu }_s^k\left(\lambda \right)=(1-f_{v,blood}^k){\mu }_{s,577}^k\frac{577}{\lambda }+f_{v,blood}^k\frac{{\mu }_{s,blood}\left(\lambda \right)}{(1+a{\left(0.5{\mu }_{s,blood}\left(\lambda \right)d_{vessels}^k\right)}^b}$$

(3)

In this first term, ${\mu }_{s,577}^k$ stands for the scattering coefficient value at 577 nm of the kth layer, and $\lambda$ in the denominator is in nm. For the second term, ${\mu }_{s,blood}\left(\lambda \right)$ is the scattering coefficient of pure blood in cm${}^{-1}$ [49]. Moreover, we set $a=1.007$ and $b=1.228$ [53], while $d_{vessels}^k$ is the mean vessel diameter of the layer, expressed in cm. The numerical values of the parameters implemented in our simulation are summarized in Table 5, and the corresponding scattering spectra ${\mu }_s^k\left(\lambda \right)$ used in our model are presented for each layer in Figure 2.

Other Optical Parameters

Both the anisotropy factor g and refractive index n were considered as wavelength-independent in our spectral range of study. For example, concerning near UV wavelengths, this assumption results in a maximum relative error of respectively 3% and 1.5% for the LE and two dermis layers [54]. The values of the former were set to $g^1=0.92$, $g^2=0.75$, $g^3=0.72$, $g^4=0.72$ and $g^5=0.9$ for each of the five layers [55], while the second values were set to $n^1=1.55$, $n^2=1.44$, $n^3=1.39$, $n^4=1.38$ and $n^5=1.34$ in each of the five layers, respectively [52].

2.3. Simulation Configuration, Method and Tools

2.3.1. Simulation Parameters Description

The simulation algorithm developed was based on the CudaMCML code from Alerstam et al. [56]. It was modified to consider the probe geometrical configuration mentioned in Section 2.1. Multi-thread calculations were performed on a DELL T630 device using the Intel Xeon E5-2683v4 2.1 GHz, AVX CPU and the GEFORCE GTX 1080 Ti GPU card cluster available within the High Performance Computing resources center “Explor” at the University of Lorraine. Technical aspects of the simulation are summarized in

Table 6. Technical characteristics of the GPU-based parallel computing and MC simulation configurations.

Photon Number	Block Number	Threads per Block	$\mathbf{\Delta }{\mathit{T}}_{\mathbf{sim}}$ (s)	${\mathit{w}}_{\mathbf{th}}$	${\mathit{t}}_{\mathbf{th}}$ (s)
$81\times {10}^8$	28	320	[900–1500]	$w_0\times {10}^{-4}$	1

. The computation time $\Delta T_{sim}$ was different for the media simulated within an interval of 900–1500 s. This duration corresponds to the launch of the ${10}^8$ photons repeated for each of the 81 wavelengths used to build a DR spectrum between 365 and 765 nm with a spectral step of $\Delta \lambda =5$ nm. This value was chosen to satisfy the best compromise between calculation time and the spectrum appearance. Ten simulations with ${10}^8$ photons were performed to quantify the trial by trial fluctuation for the the healthy mean epidermal thickness medium of phototype I (see Section 3.1). This led to the conclusion that the relative error in DR spectra compared to the average of the 10 simulations was less than 0.5% in the worst case; i.e., when the least photons are detected, in the hemoglobin absorption peak (415 nm) at the $D_4$ SD distance. The spatial resolutions according to radial distance and depth were respectively set to $dz=5\mathsf{\mu }$m and $dr=10\mathsf{\mu }$m. Finally, concerning the photon packet stopping criteria, the weight threshold was set to $w_{th}=w_0\times {10}^{-4}$, with $w_0$ the initial photon packet weight, while the time of flight threshold was established at $t_{th}=1$ s.

Details of the photon propagation algorithm and detection conditions are presented in our previous work [46]. The changes introduced to the original code [56] allow for the acquisition of simulated spectra $D{R}_{D_n}\left(\lambda \right)$, the unit of which is the photon ratio, comprised of photons that have reached the DF located at distance $D_n$.

2.3.2. Probed Depth Histograms

Probed depth histograms were determined to enable the quantitative analysis of the detected photon penetration depth. The algorithm that allowed us to obtain the latter data set is presented below. It can be applied at a given $D_n$ distance and at a given wavelength $\lambda$.

The normalization line in Algorithm 1 was implemented to eliminate the signal level and thus allow for a proper comparison between histograms of different SD distances. The occurrence of histograms is then relative. A very similar approach was presented by Arifler et al. [31] and applied to resolve spectral information from epithelium and stroma layers in SR-DR spectra. The histogram created here has the particularity of considering the deepest coordinate of the trajectory of the detected photon. This can be associated with the deepest scattering event or to a reflection event between two consecutive layers.

Algorithm 1: Maximum depth histogram of weighted detected photons

1 Here, the maximal depth $z_{max}$ is discretized to allow the design of the histogram. The width of the length interval is the spatial resolution in depth i.e., $dz=5\mathrm{\mu }$m.

3. Results and Discussion

Results of DR spectra and probed maximum depth distribution graphs are detailed in Section 3.1 for the healthy skin model with an intermediate thickness, mean dermal blood content and very fair phototype. Then, Section 3.2 focuses on the effect of the epidermis thickness, pathological state, phototype and dermal blood content on photon propagation. Based on these results and their analysis, a quantitative decomposition of the DR signal is then proposed in Section 3.3.

3.1. Detailed Results for Healthy Skin Model with Intermediate Thicknesses

3.1.1. SR-DR Simulated Spectra

Figure 3 (left) shows the DR spectra simulated for a model configuration corresponding to healthy skin of intermediate thickness (cf. Table 2), intermediate blood content (cf. Table 4) and phototype I (cf. Table 3) at each SD distance. DR intensity tends to increase toward the red spectral range due to the lower opacity of the near-infrared (near IR) wavelengths. Local minimums at 415, 540 and 575 nm correspond to hemoglobin absorption.

Figure 3 (right) shows the superposition of the normalized simulated spectra and the clinical spectra for each SD distance. Each of the subplots is composed of the the average clinical spectrum of three acquisitions ${\overline{DR}}^{clin}$ (dashed line) of this same quantity, to which we added the standard deviation (dash-doted line) and finally the simulated spectrum (continuous line), for which model parameters (epidermal thickness, phototype, dermal blood content) were determined by inverse problem solving to fit clinical data. Those experimental DR spectra was acquired on the anterior skin surface of a patient’s wrist. This skin site was clinically considered as healthy by the plastic surgeon in charge of the patient’s inclusion in the SpectroLive clinical trial (ClinicalTrials.gov Identifier: NCT02956265) approved by the Institutional Board for Ethics and Patients Safety (CPP Est III authorization: 16.06.03) and by the French National Agency for Medicines and Health Products Safety (ANSM authorization: DMDPT-RIAL/MM/2016-A00608-43). Comparing spectra from a MC simulation with clinical spectra is not an easy task [57]; thus, both sets of spectra were normalized (i.e., we set the minimum value to 0 and the maximum value to 1 in the 400–650 nm exploitable spectral band). The fitting between acquired and simulated spectra partly asserts the validity of the numerical skin model and more generally the fidelity of the simulation. A more complex skin model including skin intrinsic fluorophore absorption and emission would be interesting to increase the level of media accuracy, even if, for the considered spectral band, their influences seems to be negligible in comparison with other considered chromophores.

3.1.2. Depth Histogram and Exploitation

Depth Histogram

Quantitative information is presented in Figure 4 to describe and analyze the maximum penetration depth of the detected photons. The normalized y-axis of the histograms is the percentage of occurrence of detected photons for a given SD distance. The presence of abrupt transitions is noticeable in the histograms corresponding to a refractive index mismatch (i.e., reflection) between consecutive layers, represented by a black vertical line. Finally, a box-plot approach was chosen to extract a characteristic probed depth associated with the couple $(\lambda ,D_n)$. Thus, the three vertical colored lines for each of the histograms represent, from left to right, the first quartile (dotted) $z_{q1}$, the median (solid line) $z_{med}$ and the third quartile (dotted) $z_{q3}$. The colored area under the histogram, between the two quartiles, then represents 50% of the total area. In other words, half of the photons contributing to the $D{R}_{D_n}\left(\lambda \right)$ value have reached at most a depth between $z_{q1}$ and $z_{q3}$. For the considered wavelength—i.e., $\lambda =565$ nm—a large part of the detected photons traveled through LE and RD for the shortest and the longest distance, respectively. For the second and third distance, the maximum penetration depths of detected photons are distributed between those layers; i.e., in the PD layer. A shorter distance would be of interest for this wavelength to focus on superficial SC and LE layers.

Extraction of a Typical Probed Depth and Spectral Behavior

The histograms presented in Figure 4 are available at all wavelengths and for each of the four SD distances, from which the median and quartiles of the maximal penetration depth can be calculated. The choice to synthesize our depth histograms using a box-plot approach rather than the mean penetration depth [31] of these histograms is justified by its simplicity of interpretation from the point of view of the resulting DR signal. Indeed, by considering the signal captured for the pair $(\lambda ,D_n)$, we can decompose it into four equal contributions of different depths, for which the quartiles and medians represent the boundaries. Moreover, the interquartile range gives a good idea of the depth spread of the photons around the median, thus illustrating the enlargement of the probed volume in the three-dimensional medium.

Other alternatives were developed to extract a penetration depth for collected photons in SR-DR spectroscopy. In their study [32], Tseng et al. introduce a quantity called “interrogation depth”, defined as the z average of all detected photon scattering events in a two-layer skin model. The photon penetration was not the main focus of their study, but the paper revealed that their SD distances in a range of 1440–2880 $\mathsf{\mu }$m allow researchers to probe between 293–949 $\mathsf{\mu }$m within the 550–860 nm wavelength region. Other teams have adopted an original and experimentally transposable technique to extract a sampling depth $Z_S$ in a mono-layer medium defined by its thickness $Z_0$ and its optical coefficients ${\mu }_a$ and ${\mu }_s^{\prime }={\mu }_s(1-g)$. This approach consists of varying $Z_0$ from 0 to a sufficiently large thickness so that it no longer has any impact on the DR signal. For a very thin thickness, photons escape through the bottom of the layer, which is responsible for a very weak signal. Thickening the medium allows more photons to reach the DF, until it becomes wide enough to converge to a maximum (semi-infinite medium). $Z_S$ is then defined as the thickness at which the DR signal reaches 50% of its final maximum value [58,59]. Unfortunately, this method can not be employed for multi-layered tissue. Finally, the depth measurement can be defined as the median of the cumulative absorption of fluence absorption maps. This was applied in the study by Arimoto et al. [30] for a three-layer skin medium considering the water absorption and scattering in the near IR spectral range.

The results obtained from our box-plot approach from maximum depth histograms are presented in Figure 5, with the same graphic conventions (cf. Figure 4). The quasi-linearity of the probed depth can be seen as a function of the two varying parameters—i.e., the SD distance and the wavelength—except at the first hemoglobin absorption peak (415 nm). More importantly, it can be noticed that by choosing a couple $(\lambda ,D_n)$, it is possible to probe depths in a range of values covering the interval 50–900 $\mathsf{\mu }$m. Moreover, one can read in that work that the photons that contributed to the spectrum $D{R}_{D_4}\left(\lambda \right)$ and verified $\lambda >450$ nm mostly visited the RD layer. On the other hand, it is visible that the photons of $D{R}_{D_1}\left(\lambda \right)$ satisfying $\lambda <450$ nm remained globally confined in SC, LE and PD layers. One can also observe that the inter-quartile range widens when $\lambda$ increases for all distances $D_n$. Indeed, the probed volume is very localized for the near UV wavelengths and spreads when it approaches the near IR.

3.2. Analysis of Spectral and Probed Depth Data Related to Medium Configuration

3.2.1. Epidermis Thickness Modifications in Healthy Skin Model

Some studies have focused on the impact of the thickening of the most superficial skin layer on DR spectra acquisition in a skin multi-layered medium [60]. Among them, six levels of epidermis thickness in the range of 0–500 $\mathsf{\mu }$m have been tested in the three-layer model (epidermis and two dermal layers) [36], leading to the conclusion that its impact on the simulated spectra is very significant. In [61], the authors quantified the accuracy of estimating the optical properties of the upper and lower layers in a bi-layer medium when the thickness of the most superficial layer varies. Thus, they showed that the accuracy of estimation of the upper (respectively lower) layer increases (decreases) when the thickness of the surface layer increases over the range 0–1000 $\mathsf{\mu }$m. However, those studies do not show an evolution of the probed depth with the thickening of the superficial layer(s). The skin model configurations considered in this section are those described in Table 2. For the sake of conciseness and clarity, only the DR spectra and the graphs of probed depths are presented (see Figure 6), where only the median of the maximum depth value is shown for each of the media (quartiles are not shown).

By looking at the DR spectra, it is noteworthy to state that (i) at $D_1$, the amplitude of the DR spectrum for the maximum thickness of the epidermis is (for 365–450 nm and 650–765 nm spectral ranges) greater than that for the minimum thickness of the epidermis, while (ii) along with longer SD distances, the amplitude of the DR spectra for the maximum thickness of the epidermis tends to be lower than the one for the minimum thickness of the epidermis. If we observe the probed depth at $D_1$ for the thickest epidermis model, we notice that most of the photons verify $\lambda <550$ nm, wgucg contributes to $D{R}_{D_1}\left(\lambda \right)$ remaining confined in SC and LE layers. For the thinnest epidermis model, those photons largely visited the PD layer. The $D{R}_{D_1}\left(\lambda \right)$ spectra evolution for the different thicknesses may be related to the high scattering capacity of the SC and LE layers at the origin of a short path length return of many photons to the surface (and thus captured at $D_1$). Indeed, the figure of the probed depths shows the well-known fact [62] that, in general, the increase of the SD distance forces photons to dive deeper in the medium: the superficial back-scattering effect is no longer observed. Finally, we can notice that the wide spectral range of the excitation light and the geometry of the device (SD distances available) allow us to probe different depths without redundancy, and this is true for every epidermal thickness.

As indicated in Section 2.2, both dermal layers (PD and RD) were considered of fixed thickness. This choice results from a preliminary study in which the impact of the dermal thickness was considered secondary in comparison with the other parameters that we vary in this paper. As we have treated the impact of the epidermal thickness, we obtained results for three pairs $(z^3=150;z^4=1600)$, $(z^3=200;z^4=1800)$, $(z^3=250;z^4=2000)\mathsf{\mu }$m of PD and RD thicknesses, respectively, representing minimum, mean and maximum dermal thicknesses. The DR and maximal probed depth appear in Figure 7. It can be noticed that the impact of this parameter on the signal (left) and on the detected photon penetration (right) is very small (variations smaller than 2% in DR for $D_1$ and for $\lambda >600$ nm).

3.2.2. Epidermis Modifications in Healthy and Pathological Skin Models

We present in Figure 8 the DR spectra obtained for skin models simulating two types of cutaneous carcinoma: ulcerated and hyperkeratosed (cf. Table 2). By way of comparison, the healthy medium of average thickness also appears in the various results.

DR spectra for healthy and ulcerated skin models are very similar; i.e., they are identical for $D_2$, $D_3$ and $D_4$. Slight distinct amplitudes are observed at $D_1$ but with a globally identical spectral shape. Indeed, the disappearance of the SC layer, which differs from both, can be seen as an overall decrease in the thickness of the epidermis $z^{Epi}=z^1+z^2$. The back-scattering effect for $D_1$ is thus more important when $z^{Epi}$ is greater (healthy medium) than when it decreases (ulcerated medium). The associated signal increases and then fades when the SD distance increases.

Unlike healthy as well as ulcerated skin models, the results obtained for the hyperkeratosis model show that the majority of photons detected, whatever the SD distance, remained confined in SC (for $D_1$ and $D_2$) and LE (for $D_3$ and $D_4$). This information is confirmed by the shape of the DR spectra, since the dermal hemoglobin signature (pits at 415, 540 and 575 nm) appears only slightly. One can again notice the superficial back-scattering effect responsible for much larger amplitudes for every $D_n$, except for $D_1$ for the near IR wavelengths. Such results emphasize the importance of taking into account the skin layer thickness of each pathological case to properly diagnose it.

3.2.3. Phototype Modifications in Healthy Skin Model

Figure 9 shows the DR spectra and probed depth graphs obtained for the healthy skin medium of intermediate thickness (cf. Table 2), mean dermal blood content (cf. Table 4) and for the four skin phototypes considered (cf. corresponding values of volume fraction of melanosome $f_{v,mel}^2$ in Table 3 and LE absorption coefficient ${\mu }_a^2$ in Figure 2, top right insert).

The decay of the signal in the near UV and the beginning of the visible light spectrum, when the skin darkens, was expected and is the most obvious characteristic. Rather than simply diffusing, the LE becomes more absorbent with increasing pigmentation. The optical barrier of the epidermis (including SC+LE) becomes increasingly difficult to penetrate for photons verifying $\lambda <450$ nm. For $D_1$ (shortest paths), the maximum depth probed by the detected photons decreases when $f_{v,mel}^2$ increases. The information carried by $D{R}_{D_1}(\lambda <450)$ appears thus to be very superficial with $z_{max}$ confined in LE layer. The very weak signature of the hemoglobin peak in the spectra for the dark phototypes confirms the previous assertion, since photons that reach the vascularized layers are rare. However, for this same distance and for the rest of the spectrum ($\lambda >500$ nm), it is interesting to note that the probed depth, in the PD, is virtually unchanged with the phototype. The path (i.e., the trajectory) traveled by these photons seems then to be globally identical, and only the absorption during the forward passage—the way back in the LE—justifies the signal drop in the DR spectrum. This is partly explained by the fact that the change in phototype in our model is simply considered by an increase in the melanosome volume fraction. A more complex model taking into account the volume distribution and shape of melanosomes varying with phototype would probably have more impact on the reduced scattering coefficient and thus on the photon trajectories in the medium. This specific behavior observed at short distances ($D_1$), which is different for short ($\lambda <450$ nm) and long ($\lambda >500$ nm) wavelengths, might be of particular interest, namely to improve the robust and precise estimation of the absorption coefficient ${\mu }_a^2\left(\lambda \right)$ (and thus by extension a value for the volume fraction of melanosome) of the LE layer, specifically when solving the complex inverse problem of multiple optical properties in multilayer skin models.

For the other, longer SD distances, one can notice that all the detected photons reach at least the RD layer for $D_2,D_3,D_4$ and $\lambda >500$ nm, and it is possible to observe the sensitivity of $z_{max}$ to the hemoglobin absorption peak (at 415 nm), especially for $D_2$ and $D_3$. The increasing probed depth with phototype, for all longer SD distances $D_n$ (n > 1), can be explained by considering two complementary behaviors during the photon penetrating trajectory from SC to RD layers. The first is to dive as little as possible into the medium and thus reduce the optical path traveled. The second is to quickly reach the PD in order to escape the LE with its high absorbing capacity. For fair phototypes, where the LE absorption is moderate, the first behavior seems to have more impact; that is, detected photons resulted from shallower paths. On the other hand, for dark phototypes, it is more probable to detect photons that cross the LE layer with the shortest path lengths; i.e., paths that are as vertical or straight as possible. These photons thus need more time to make their way back towards the detection fiber: the probed depth is then more important.

3.2.4. Blood Content in Dermal Layer Modifications in Healthy Skin Model

Figure 10 shows the DR spectra and probed depth graphs obtained for the healthy skin medium of intermediate thickness (cf. Table 2) for the fair phototype I (cf. Table 3) and for the three media representing dermal layers with more or less blood content (cf. Table 4), for which absorption and scattering coefficients of the two concerned layers (PD and RD), respectively, appear on the middle-right and bottom-right subplots of Figure 2. The most remarkable and expected observation is the relative signal decrease in the three hemoglobin absorption peaks (415, 540 and 575 nm). This can be linked to the probed depth plot where, for distances $D_2$ to $D_4$, the spectral region where hemoglobin is absorbed (i.e., 365–630 nm) is affected by a decrease in the detected photon penetration with increasing blood concentration.

Concerning $D_1$, the interpretation is even simpler, because of the quasi-non-influence of the quantity of blood on the path of the detected photons. In $D{R}_{D_1}\left(\lambda \right)$, we notice that for the first part of the spectrum (365–600 nm), the blood concentration is responsible for an attenuation of the signal, while it is associated with a gain of intensity for the rest of the spectrum (600–765 nm). The analysis of the optical properties on these spectral bands leads us to the interpretation that the first phenomenon is caused by the absorption of hemoglobin in the dermal layers. For the second phenomenon, the absorption of blood becomes negligible, and only the scattering coefficients remain significantly different (about 15% between the two extremes maximum and minimum $f_{v,blood}$). This means then that the scattering power caused by the increase of the blood concentration of the dermal layers (particular PD for this shorter distance $D_1$) is at the origin of many short optical paths between SF and DF and therefore responsible for this non-obvious signal increase.

3.3. A Proposal of a Model for DR Signal Decomposition by Skin Layer

In order to more quantitatively determine the content of the global DR signal in terms of layer contributions, we propose to exploit the results of analysis of the aforementioned developed maximum probed depth histograms and spectra to decompose the total signal $D{R}_{D_n}\left(\lambda \right)$, as a sum of the DR contributions related to each of the five skin layers considered, such that for each SD distance $D_n$,

$$\begin{array}{cc}\hfill D{R}_{D_n}\left(\lambda \right)& =D{R}_{D_n}^{SC}\left(\lambda \right)+D{R}_{D_n}^{LE}\left(\lambda \right)+D{R}_{D_n}^{PD}\left(\lambda \right)+D{R}_{D_n}^{RD}\left(\lambda \right)+D{R}_{D_n}^{SF}\left(\lambda \right)\hfill \end{array}$$

(4)

where $D{R}_{D_n}^{layer}\left(\lambda \right)$ is composed of photons that reached at deepest the $layer$ appearing in exponent verifying $d_{min}^{layer}<z_{max}<d_{max}^{layer}$ with $d_{min}^{layer}$ and $d_{max}^{layer}$, respectively, the minimum and maximum depth coordinates of the $layer$. In practice, this decomposition was obtained by integrating the histograms of probed depths (cf. Section 3.1.2 and Figure 4) between the upper and lower border of each skin layer. Figure 11 illustrates this decomposition for the four SD distances and the healthy skin model with intermediate thicknesses, mean dermal blood content for SC and LE layers and phototype I. This “deepest layer probed” separation criterion has already been proposed in the literature. In [36], Fredriksson and his team studied the path length distributions for detected photons in a three-layered skin model, although they did not propose any decomposition approach. However, a different “layer”-based approach of DR signal decomposition has been proposed in a few studies. In [31], the contribution of the detected photon associated with the k-layer was defined as the quotient of the number of scattering events undergone in this layer on the total number of scattering events between SF and DF. In [35], the authors proposed to explain the skin layer sensitivity by the mean partial path length, which consists of assigning to the layer k the relative weight $l_k/l_{tot}$ in the DR signal, with $l_k$ the path length in the layer k and $l_{tot}$ the total path length between SF and DF.

The approach proposed in the present contribution differs from the latter by the following points. First, our proposed decomposition provides a straightforward and simple interpretation regarding the layer reached by the detected photons making the global DR signal, compared to less non-trivial interpretations linking the mean path length or the ratio of scattering events to the layer’s weight in the DR signal. Second, as a consequence, the proposed decomposition should help in designing a method to estimate sequential optical properties by studying the layers from the shallowest to the deepest and by providing new, a priori knowledge on the selection of the sensitive distances and wavelength range. This approach will bring a physical and anatomical sense to the estimation and should be less time consuming due to the successive isolation of the parameters.

The analysis to be drawn from Figure 11 is in accordance with that performed qualitatively in the Section 3.1 but allows us to better explain the DR signal in layer decomposition. To further compare the proportion of a signal from a layer, for a fixed couple $(\lambda ,D_n)$, we defined the weight fraction of the layer contribution in percentage, as follows:

$$w_{D_n}^{layer}\left(\lambda \right)=100\times \frac{D{R}_{D_n}^{layer}\left(\lambda \right)}{D{R}_{D_n}\left(\lambda \right)}(\%)$$

(5)

Finally, Figure 12 shows the associated evolution of $w_{D_n}^{layer}\left(\lambda \right)$ defined in Equation (5) as a function of the wavelength and the SD distance. This decomposition turns out to be very useful to extract the spectral band and SD distance of a DR signal that is susceptible to carry information associated with local depth information. For instance, the spectral “signature” of elastosis (accumulation of disordered elastic tissue in the PD due to carcinogenesis) should be looked for in the UV range (365 < $\lambda$ < 400 nm) at the $D_2$ distance, where more than 40% of the DR signal is constituted of photons that have visited the concerned layer.

It is worth noting that all the geometrical and optical parameters considered in the model described above were chosen within standard ranges of values, accounting for anatomical and physiological variability. The aim of the present work was not to perform a complete variation study taking into account those variability intervals but rather to investigate the variation of the dependence of epidermal thickness (see Section 3.2.1), pathological state (see Section 3.2.2), phototype (see Section 3.2.3) and dermal blood content (see Section 3.2.4) considering a unique set of basic but representative values for all other parameters of the model. A possible extended study could be led with other parameters, such as the individually variable scattering coefficient [63], the blood oxygen saturation [60] or the average blood vessel diameter [64], even if those parameter have a more limited impact on the layers’ optical coefficients and thus on the penetration of photons in depth. It would also be interesting to consider 3D tumor models, since the assumption of semi-infinite layers imposed by CudaMCML is a limitation. A study shows that the validity of this assumption depends on the type of tumor simulated (squamous cell carcinoma in epithelium or a basal cell carcinoma in upper dermis), the angular incidence of the exciting fiber, the height of the tumor and its width relative to the SD distance at which the signal is acquired [65].

All the results established so far allow us to extend our discussion towards the targeted clinical application of SR-DRS devices. Firstly, the in-depth elements characterizing the path of the detected photons within the various skin layers allow us to assert that the geometrical layout of the probe and its spectral configuration provide access to the different layers of the skin with a sensitivity depending on the chosen couple $(\lambda ,D_n)$ and on the physiological and morphological skin configuration. In particular, the results showed that the photons contributing to the DR spectra were likely to carry the discriminating local information that distinguishes healthy and pathological tissues. It should be noted, however, that the results presented here are valid for the case of fiber-based contact-type DRS configurations, and under the assumption that the optical properties of each of the skin layers do not significantly vary laterally within the ∼2 mm diameter circular region covered by the multi-fiber optics probe tip. Moreover, for the estimation of diagnostically relevant optical parameters from in vivo SR-DRS measurements through an inverse problem solving method [66,67], such results should be useful to understand this iterative process better by defining a cost function adapted to the layer structure of the skin.

4. Conclusions and Prospects

In order to learn more about the depth behavior of detected photons and thus allow for a better interpretation of the simulated and measured diffuse reflectance signals, we decided to carry out the current numerical study. A photon propagation algorithm was adapted to match the measurement features of the SpectroLive device, which is a spatially-resolved diffuse reflectance spectroscopy device used in clinical applications. To do so, a numerical five-layered skin model, including realistic absorption and scattering properties from main intrinsic chromophores, was defined in order to reach a certain fidelity between the simulation and the experimental acquisition. The histograms of probed depths allowed us to extract the typical depths at which the photons took the return path to the collection fiber. The results dependent on the skin type placed under examination and the fate of photon behavior in terms of depth penetration for different physiological and pathological states were investigated. Indeed, diffuse reflectance spectra and probed depth graphs were analyzed for different sites on the body (by varying the thicknesses of the stratum corneum and the living epidermis), for different pathological states (ulcerated and hyperkeratosed carcinoma) and finally for different phototypes (by varying the volume fraction of melanosome in living epidermis).

Finally, to further the understanding of the diffuse reflectance signal, a layer decomposition model of the latter was proposed and applied to a healthy medium of intermediate thickness, mean dermal blood content and phototype I (most carcinomas occur in the fair-skinned population). By providing quantitative layer contributions in the resulting signal, these results refine the existing relationships between the wavelength, source–detector distance and probed skin layer. This relationship will be further used to interpret spectra recorded during a clinical trial that aims at evaluating the SR-DRS diagnostic accuracy for the differential diagnosis of skin carcinomas. More specifically, the exposed numerical evidence may facilitate the development of analytical photon propagation models that enable the inverse estimation of diagnostically relevant optical parameters or allow researchers to distinguish a healthy from a pathological tissue using spectral criteria, from in vivo clinical spectra. Such a numerical study may also be used to standardize measured spectra taking into account intra- and inter-individual variability. Currently, ratios of measured spectra are proposed to eliminate this variability, but such calculation methods can sometimes be impaired by artifacts of measurements.