Recovery of Optical Transport Coefficients Using the Diffusion Approximation in Bilayered Tissues: A Theoretical Analysis

1. Introduction

Time domain diffuse optical spectroscopy (TD-DOS) is a quantitative, non-invasive technique that has been widely explored across several biomedical applications [1,2,3,4,5]. The quantitative nature of this technique relies on recovering the wavelength-dependent transport coefficients –the absorption coefficient ${\mu }_a\left(\lambda \right)$ and the reduced scattering ${\mu }_s^{\prime }\left(\lambda \right)$ to optically characterize the medium from which measurements are obtained. In TD-DOS, the measured temporal attenuation and broadening of an input source (which experimentally is usually a laser having pulse durations of 10-100 ps) is mathematically modeled using theoretical or numerical methods to recover the optical coefficients ${\mu }_a\left(\lambda \right)$ and ${\mu }_s^{\prime }\left(\lambda \right)$[6,7,8,9,10]. In human tissues, accurate recovery of ${\mu }_a\left(\lambda \right)$ at multiple wavelengths allows for reconstructions of component chromophore concentrations such as hemoglobin, water, and/or lipids, while ${\mu }_s^{\prime }\left(\lambda \right)$ is influenced by tissue microstructure [6,9,11,12,13,14]. These chromophore concentrations are used to estimate tissue vascular oxygen saturation, total blood volume, hydration, and/or cellular density in vivo [15,16,17].

Theoretical models are critical for recovering optical transport coefficients from measurements, as they provide the quantitative platforms to model light transport in turbid media. A commonly used theoretical approach to quantify TD-DOS measurements is diffusion theory (DT)[18,19,20,21,22]. DT-based approaches typically consider the tissue medium as optically-homogeneous and semi-infinite in extent across in various applications, including tracking and monitoring cerebral oxygenation and cancer development in clinical settings [3,8,19,23,24]. Though modeling biological tissues as a semi-infinite homogeneous medium allows for quantitative analysis, it is inherently insensitive to any inhomogeneities of the medium, which has been shown to lead to large errors in recovering physiological parameters [6,7,9,25,26,27].

Representing tissue as a bilayered medium would be more realistic than as a semi-infinite medium –e.g., a human head can be better modeled as having a top layer (representing scalp and skull) and a bottom layer (for the cerebral spinal fluid, gray and white matter) and limb tissues can be modeled as a layer of muscle buried under a upper-layer of fat and skin [6,28]. The use of layered models in DT would permit reconstruction of optical coefficients for each layer, independently, from the measured TD reflectance [11,17,29]. Bilayered models can thus naturally detect transport coefficients of deeper layers better, and could also be used to suppress cross-talk from superficial layers [21,30,31]. Such depth-dependent reconstructions would enhance the accuracy of cerebral and muscular hemodynamic assessments [29,30,32,33,34,35].

Analytical solutions of TD reflectance using DT in bilayered media have been reported previously [18,36,37,38,39]. However, they are not commonly used as the implementation of these solutions require computing complex analytical expressions that often include hyperbolic functions, leading to numerical instabilities because of finite machine-precision. Thus, these solvers have high computational complexity and costs [7,35,40,41,42,43]. Recently, we have developed an easy-to-use, open-source library - lightPropagation.jl - that is numerically stable and highly optimized (with a single calculation for a forward solution taking less than $0.5$ ms on standard laptop computers) for computing time domain reflectance in multilayered tissue models via DT [44].

Here, we seek to establish the accuracy and precision of using lightPropagation.jl for inverse reconstruction of optical coefficients using time-domain reflectance acquired from using Monte Carlo (MC) simulated of time-resolved reflectance in bilayered tissue models. We assess the performance of the inverse model in recovering the absorption and reduced scattering coefficients. To ensure appropriate relevance to biomedical applications, we selected optical properties for tissue models simulated to match values reported in literature [28]. We use wavelength-dependent reconstructions to derive functional biological markers of fractional oxygen saturation (${\mathrm{SO}}_2$) and hemoglobin concentrations for each tissue layer and quantify the errors associated with reconstructions to establish limits for accuracy and precision of recovered optical coefficients and functional endpoints.

2. Materials and Methods

2.1. Bi-Layered Tissue Models

Two types of tissue models, representing head and limb muscle tissue were used to generate TD reflectance using MC simulations and are depicted by Figure 1(a). Both tissue models were represented as cylinders with large radius (set to 100 cm), had upper layer thickness of 1 cm, and had semi-infinite bottom layer (thickness of 100 cm) to compute MC and DT solutions. Input values of the absorption coefficient for each layer in each tissue model was determined by solely by two chromophores - oxygenated ([HbO]) and deoxygenated hemoglobin ([dHb]) [45]. Table 1 lists values for inputs used to generate optical coefficients in each layer, for each tissue model.

Eq.(1) was used to determine ${\mu }_a$ for each layer by specifying values for $\left[THb\right]$ and ${\mathrm{SO}}_2$ and using molar extinction coefficients ${ϵ}_{HbO}\left(\lambda \right)$ and ${ϵ}_{dHb}\left(\lambda \right)$ from literature [28,46]. The reduced scattering coefficients (${\mu }_s^{\prime }$) for the top and bottom layer were calculated using Eq.(2) using literature values for A and b, for each layer and tissue type [12,28].

$${\mu }_a\left(\lambda \right)=S{O}_2·\left({ϵ}_{Hbo}\left(\lambda \right)·THb+{ϵ}_{dHb}\left(\lambda \right)·(1-THb)\right)$$

(1)

$${\mu }_s^{\prime }\left(\lambda \right)=A{\left({\displaystyle \frac{\lambda }{500}}\right)}^{-b}$$

(2)

64 different models were generated at three commonly used near infrared wavelengths (690, 760, and 850 nm) to span the isosbestic points of hemoglobin [6,14,20]. The full dataset consisted of $384(64\times 3\times 2)$ bilayered tissue models, on which we conducted the analysis. Ranges for optical coefficients at each of the three wavelengths used, for each layer, are listed in Table A1, Table A2 and Table A3 for 690 nm, 760 nm and 850 nm, respectively.

2.2. Monte Carlo Simulation of Reflectance Signals

A previously used Monte Carlo code for photon transport was used to simulate time-resolved reflectance from each of the 384 tissue models [47,48]. MC simulations for all models were computed with $3\times {10}^8$ photons. For each simulation, the input source was a pencil beam (representing a delta-function) incident at the origin of the tissue model and the TD diffuse reflectance was recorded between 1cm to 3 cm in intervals of $0.25$ cm with timing relative to the input delta-function source. Reflectance data was stored between $0-7$ ns with temporal resolution of $0.01$ ns simulations were run on a high-performance computing cluster (with each node of the cluster hosting Intel Xeon Gold processors). Each MC simulation required about 15 hours to execute and simulations were run in parallel (on separate compute nodes) to increase output efficiency. For analysis, each simulated time-resolved reflectance was smoothed (using a moving-window of span 5 or a temporal window of 0.05 ns), normalized (peak intensity $=1$), truncated (by retaining values greater than $80\%$ of peak intensity of the rising edge to eliminate early arriving photons and lesser than than ${10}^{-3}\%$ of the tail to limit statistical noise), and lastly log-transformed to obtain $R_{MC}(\rho ,t)$. Although we simulated reflectance data between 1-3 cm, we only use reflectance at three specific source detector separations (SDS) of 1.5, 2 & 2.5 cm because the noise in MC simulations at SDS of 3 cm was too high for use here.

2.3. Forward-Modeling of Diffuse Reflectance Signal

The numerical DT solver lightPropagation.jl was used to obtain the TD-reflectance (the temporal point spread function or TPSF), for any given bilayered cylindrical tissue model [27]. The SDS and was set to match those from the MC simulations, the refractive index of each layer was fixed to be $1.4$ which is consistent for biological tissues from previous work [49] and the top layer thickness was set to 1 cm. To compute the reflectance signal, for each tissue model required inputs for the four transport coefficients, ${\mu }_{a_1}$, ${\mu }_{a_2}$, ${\mu }_{s_1}^{\prime }$ and ${\mu }_{s_2}^{\prime }$. The computed signal was was then truncated, normalized and log-transformed following the same protocol used for processing the MC reflectance signal, to obtain $R_{DT}\left(t\right)$.

Figure 1(b) shows representative data for TD-reflectance simulated (in symbols) and predicted from DT (solid lines) at three different wavelengths (colors) for a muscle tissue model. The values of optical absorption and scattering were calculated at each wavelength using Eq.(1) and Eq.(2) where $\left[THb\right]$ was $18\mu$M for the top layer and $27\mu$M and for bottom layer, with both layers having ${\mathrm{SO}}_2$ of 50%.

2.4. Inverse Fitting of Reflectance

Inverse modeling sought to recover the four optical transport coefficients (the absorption and scattering for each layer) by fitting each time-resolved reflectance $R_{MC}\left(t\right)$, at each SDS, for each of the 384 tissue models. Fits were obtained iteratively using a Levenberg-Marquardt non-linear optimizer (LVM) (LsqFit package in Julia version 1.9.1) to compute a set of inputs to compute $R_{DT}\left(t\right)$ that best matched the input $R_{MC}\left(t\right)$. For DT calculations, the top layer thickness was assumed to be known (kept fixed for all analysis here at 1 cm), the refractive indices of both layers were held fixed at $1.4$ and the SDS was set to match $R_{MC}\left(t\right)$. During optimization the optical coefficients were constrained to be bounded between $0.01<{\mu }_a<0.5$ for absorption and between $2<{\mu }_s^{\prime }<30$ for scattering.

A normalized mean absolute error ($\zeta$) value was computed to evaluate goodness of fits between $R_{DT}\left(t\right)$ and $R_{MC}\left(t\right)$ using equation 3

$$\zeta ={\displaystyle \frac{\mathrm{MAE}}{\overline{|R_{MC}\left(t\right)|}}}$$

(3)

Here, $\overline{|R_{MC}\left(t\right)|}={\displaystyle \frac{1}{N}}{\sum }_{i=1}^N\left|R_{MC}\left(t_i\right)\right|$ represents the average of the absolute value of the reflectance signal, and $\mathrm{MAE}={\displaystyle \frac{1}{N}}{\sum }_{i=1}^N|R_{MC}\left(t_i\right)-R_{DT}\left(t_i\right)|$. N represents the total number of time bins in $R_{MC}\left(t_i\right)$ and $N\approx 180\pm 100$ across the entire data-set, due to data truncation of $R_{MC}\left(t\right)$ that was noted above.

For each simulated reflectance $R_{MC}\left(t\right)$, inverse fitting was repeated until $\zeta \le 0.03$. If $\zeta$ was higher than this after 25 attempts the solution with lowest $\zeta$ was used as the fit. In Figure 1(b) inverse fits are shown by dashed lines and had $\zeta =0.017,0.015\&0.018$ for 690 nm (indigo), 760 nm (orange) and 850 nm (red) models, respectively.

Lastly, we also used the semi-infinite DT approximation to fit $R_{MC}\left(t\right)$ using lightPropagation.jl [27]. The forward fitting used in lightPropagation.jl replicates the expression for semi-infinite DA reflectance in Ref.[37]. Inverse fitting for the semi-infinite model mirrored the bilayered approach, and used the optimizer algorithm to extract two coefficients ${\mu }_a$ and ${\mu }_s^{\prime }$ to the same goodness-of-fit threshold as used in the bi-layer reconstructions ($\zeta =0.03$) and used the same search ranges for each optical coefficient ($0.01<{\mu }_a<0.05$; $2<{\mu }_s^{\prime }<30$).

2.5. Retrieval of Functional Endpoints

As most real-world applications of TD-DOS involve extracting tissue physiological parameters from optical coefficients we emulated that process to the reconstructed absorption coefficients at multiple wavelengths to obtain the functional parameters of $\left[dHb\right]$ and $\left[HbO\right]$ using 4

$$\left[\begin{array}{c}{\mu }_a\left({\lambda }_1\right)\\ {\mu }_a\left({\lambda }_2\right)\\ {\mu }_a\left({\lambda }_3\right)\end{array}\right]=\left[\begin{array}{cc}{ϵ}_{HbO}\left({\lambda }_1\right)& {ϵ}_{dHb}\left({\lambda }_1\right)\\ {ϵ}_{HbO}\left({\lambda }_2\right)& {ϵ}_{dHb}\left({\lambda }_2\right)\\ {ϵ}_{HbO}\left({\lambda }_3\right)& {ϵ}_{dHb}\left({\lambda }_3\right)\end{array}\right]\left[\begin{array}{c}\left[HbO\right]\\ \left[dHb\right]\end{array}\right]$$

(4)

Here, values of $\left[dHb\right]$ and $\left[HbO\right]$ were estimated using the same optimizer (as used for TD-reconstructions) to ensure only non-negative values of $\left[dHb\right]$ and $\left[HbO\right]$ were allowed as solutions by the optimizer. This eliminated any solutions with negative values for concentrations and thus remain physiologically meaningful. The chromophore concentrations were used to compute $\left[THb\right]$ (as $\left[dHb\right]$+$\left[HbO\right]$) and fraction ${\mathrm{SO}}_2$ (as $\left[HbO\right]/\left[THb\right]$) which were each then compared to input values used for each model (shown in Table 1). This process was repeated for each SDS.

3. Results

The inverse fitting protocol was employed to recover optical coefficients from each of the 384 MC models at three SDS. Reconstructions for the absorption (${\mu }_{a_1}$, ${\mu }_{a_2}$) and reduced scattering coefficients (${\mu }^{\prime }{s}_1$, ${\mu }^{\prime }{s}_2$) at 760 nm are shown in Figure 2. The mean goodness of fits at $1.5$ cm and 2 cm had $\zeta \le 0.01$ , across all models, with head reconstruction performing marginally better (average $\zeta$(head) was $0.005$ lesser than $\zeta$(muscle)). However, for SDS of $2.5$ cm, $\zeta$ increased to $0.02$ which we hypothesize was due to an increase in statistical noise in the MC simulations at longer SDS. Each iteration for convergence took approximately $15s$ on average and most $R_{MC}\left(t\right)$ fits needed $10-15$ iterations to meet the required convergence threshold. In Figure 2 each marker corresponds to some fixed value for the optical coefficient being plotted, while the other three coefficients were varied. For instance, each marker in Figure 2(a) is from 16 different tissue models that shared the same ${\mu }_{a_1}$ value, while error bars are the standard error across all those 16 reconstructed values of ${\mu }_{a_1}$ as ${\mu }_{a_2}$, ${\mu }_{s_1}^{\prime }$ and ${\mu }_{s_2}^{\prime }$ were permuted.

Figure 2(d) clearly show that the largest reconstruction errors were seen bottom layer scattering ${\mu }_{s_2}^{\prime }$ with reconstruction errors, generally $>5c{m}^{-1}$ for all SDS. The difficulty of layered DT in reconstructing deep-layer scattering estimation observed here has been reported previously [18,50,51]. Recovery of ${\mu }_{a_1}$ was inconsistent and showed significant errors for some of the tissue models tested but the mean reconstruction error was lower than $0.03c{m}^{-1}$ at SDS $=1.5$ cm and 2 cm, with good agreement in many cases, as noted in Figure 2(a) by the overlapping markers. Larger errors in reconstructed values for ${\mu }_{a_1}$ were noted at SDS of 2.5 cm growing to nearly $0.1c{m}^{-1}$. Recovery of ${\mu }_{a_2}$ and ${\mu }_{s_1}^{\prime }$ were the highly accurate with errors $<0.01\mathrm{and}<2c{m}^{-1}$, respectively, for all SDS tested. These trends preserved at the other two wavelengths (Supplementary Figure A1 for 690 nm and Figure A2 for 850 nm)

The impact of reconstruction errors on derived physiological parameters was next investigated by computing the reconstructed concentrations of ($\left[dHb\right]$) and ($\left[HbO\right]$), as illustrated in Figure 3. Reconstructed absorption coefficients at 690 nm, 760 nm, and 850 nm were used to derive the hemoglobin concentrations. Due to the sparse and non-uniform distribution of hemoglobin concentrations used in simulated tissue models, the x-axis are scaled to reflect the actual concentration values and thus the values from bar heights need to interpreted in conjunction with their x-axis locations, to assess performance of reconstructions.

Results in Figure 3(a) and Figure 3(b), demonstrate that top-layer hemoglobin retrieval was less accurate than reconstructions in the bottom layer. Most of the reconstructed values in the upper layer either overestimated (for smaller values of $\left[dHb\right]$) or underestimated (for larger values of $\left[dHb\right]$). For instance, the reconstructed value for $\left[dHb\right]$ was $6\pm 1\mu$M when its true value was $1.8\mu$ M for muscle (the third bar in Figure 3(a)). However, the spread (standard error) remained within $10\mu$M across all models which is sufficient for clinical use [52]. The reconstructions for bottom-layer concentrations were accurate overall, with reconstructed values for both $\left[dHb\right]$ and $\left[HbO\right]$ being well within the computed standard error (Figure 3(c) and Figure 3(d)). This trend was true at all SDS used, but data at larger SDS were associated with a larger spread in reconstructed values.

Functionally, the most common biomarkers sought by real-world applications using TD-DOS include the total hemoglobin concentration $\left[THb\right]$ and fractional oxygen saturation ${\mathrm{SO}}_2$ which are computed using wavelength-dependent optical transport coefficients. Figure 4 shows the true (input) values used for ${\mathrm{SO}}_2$ and total hemoglobin concentration $\left[THb\right]$ values for the top and bottom layers vs. those derived from the inverse DT. As expected and seen in the figures, errors in hemoglobin concentrations of the top layer directly affect the accuracy of the derived physiological markers of ${\mathrm{SO}}_2$ and $\left[THb\right]$, leading to significant deviations from their true values.

These appear as incorrect median reconstructed values as well as higher standard error all SDS. However, both precision and accuracy of reconstructed values for the bottom layer (Figure 4(c) and Figure 4(d)) were excellent with the reconstructed median values of both ${\mathrm{SO}}_2$ and $\left[THb\right]$ matching the expected physiological ranges to better than 3% across all SDS. One exception to this was seen for muscle models where $\left[Thb\right]=50\mu$M and ${\mathrm{SO}}_2=50\%$ (for the bottom layer) that exhibited large standard error in derived values (the median value was in agreement to the true value to within 3% again) . These same tissue models shown in Figure 3 also had large standard errors in reconstructed $\left[dHb\right]$ and [$HbO$] concentrations.

As a final point of analysis, we reconstructed optical coefficients from the bilayered $R_{MC}\left(t\right)$ using the semi-infinite DT expression for TD reflectance. Since the SI model reconstructed optical coefficients for a homogeneous model (hence only recovers two optical coefficients), we compared the recovered absorption coefficient to both ${\mu }_{a_1}$ and ${\mu }_{a_2}$ and repeated that for the reduced scattering as well. Figure 5, shows the median percent error (bars) along with the standard error (error bars) for all MC models (head and tissue) simulated at 15 mm SDS. It is clear that even with the difficulties in reconstruction of four optical coefficients using bi-layer DT models, the retrieved layer specific transport coefficients were always more accurate than those estimated by the semi-infinite model with a notable exception for retrieval of ${\mu }_{s_2}^{\prime }$. This is also true for SDS $=2$ cm but at SDS $=2.5$ cm, the reconstruction of ${\mu }_{a_1}$ was significantly affected (Figure A5)

4. Discussion

We verified the accuracy of a recently reported, open-source, numerical solver of DT in cylindrical coordinates, lightpropagation.jl, to function as an inverse solver in bilayer tissue models using TD reflectance simulated at multiple SDS. 384 bilayered models were simulated using MC for a range of tissue optical properties that have been reported for head and muscle tissues. The goodness of fit was established by computing the Normalized Mean Absolute Error ($\zeta$), which on average was (across all models) lower than 0.01 at SDS $=1.5$ and 2 cm while $\zeta \ge 0.02$ at SDS of $2.5$ cm. A change of 0.01 in $\zeta$ in goodness of fit significantly impacted accuracy of reconstructions of the top layer absorption ${\mu }_{a_1}$ as seen in Figure 2(a) (and Figure A1(a) Figure A2(a)). Thus $\zeta$ could potentially signal a measure of confidence in the retrieved value of optical coefficients for each fitted reflectance and could also serve as a measure of signal quality in experimental use.

A threshold of $\zeta \le 0.02$ in inverse fits of the reflectance was sufficient for reconstruction of top layer scattering (${\mu }_{s_1}^{\prime }$) and the bottom layer absorption (${\mu }_{a_2}$) with mean errors of lesser than 5% (for ${\mu }_{a_2}$) and lower than 3% (for ${\mu }_{s_1}^{\prime }$), across all SDS tested and as shown in Figure 2(b)-(c) Figure A1(b)-(c) and Figure A2(b)-(c). However, the scattering of the bottom layer (${\mu }_{s_2}^{\prime }$) could not be reconstructed accurately and showed average error of more than 60% across all SDS, indicating layered DT is largely unperturbed by changes in deep-layer scattering, as reported before [7,18,39,50,51]. Although recovery of top layer absorption ${\mu }_{a_1}$ was achieved, it was inconsistent, with mean errors $\le 15\%$ for all 384 models at SDS $=1.5$ and 2 cm that increased to more than 50 % for SDS of 2.5 cm, which indicates that the reconstruction upper layer absorption was highly sensitive to SNR.

The input coefficients for each tissue model simulated here were obtained by assuming only two chromophores ($\left[dHb\right]$ and $\left[HbO\right]$) were present in each layer, the values are shown in Table 1. In the inverse sense, these chromophore concentrations had to be obtained from reconstructed absorption coefficients, for each layer, at each wavelength. Thus, as expected, errors in recovery of ${\mu }_{a_1}$ impacted the estimation of oxygenated and deoxygenated hemoglobin concentrations for the upper layer, with mean error being higher than $5\pm 2\mu$M but lower than $10\pm 3\mu$M for SDS $1.5\mathrm{and}2$ cm (Figure 3(a)-(b) and Figure A3(a)-(b)). The upper bound on of top-layer concentration errors grew to more than $15\mu$M for both chromophore concentrations at SDS of $2.5$ cm (Figure A4(a)-(b)). The bottom-layer reconstructions demonstrated good agreement with ground truth values (to better than $3\mu$M at all SDS) while also having low inter-model variance.

Accurate recovery of individual chromophore concentrations is important, as clinically relevant parameters such as oxygen saturation (${\mathrm{SO}}_2$) and total hemoglobin concentration ($\left[THb\right]$) are derived from the combined contributions of $\left[HbO\right]$ and $\left[dHb\right]$ as markers of tissue health and metabolic demand. Errors in reconstruction of the top layer chromophore concentrations thus impacted retrieval of top layer ${\mathrm{SO}}_2$ (error > 5% for all SDS). The top layer $\left[THb\right]$ values were retrieved with mean errors greater than 15% for shorter SDS (1.5 and 2 cm), which increased to more than $70\%$ for larger SDS (Figure 4(a)-(b)). Again, this is consistent with the idea that the diffuse reflectance collected at larger SDS would be most insensitive to upper layer optical coefficients.

These results highlight that bilayered DT has reduced sensitivity to the top layer absorption and reconstruction of ${\mu }_{a_1}$ is sensitive to signal quality. However, recovery of bottom layer endpoints was accurate with ${\mathrm{SO}}_2$ and $\left[ThB\right]$ being estimated to better than 3% across all SDS. This result is particularly encouraging since it signifies that real-world applications could enhance the quantitative sensitivity of functional properties recovered from deeper tissue layers, while being immune to changes in superficial layers, by using improved analytical models for reconstruction of TD reflectance [11,53,54].

A practical consideration that remained significant was the computational cost involved with reconstructions using the bilayer inverse model. The inverse fits converged in under 100 milliseconds using the analytical DT expression for a SI model, while for bilayer reconstructions inverse fits were approximately three orders of magnitude slower (taking about 150 s for each model). All inverse calculations were run on nodes of a high-performance computing cluster (with each node having an Intel Xeon Gold processor). For the 384 models analyzed at 3 SDS, the total computation time was almost 40 CPU-hours. This disparity presents a significant bottleneck for real-time or high-throughput analysis and motivates the development of computationally efficient strategies or machine-learning based models for acceleration.

Ultimately, the computational costs were recovered by the layered model as it consistently outperformed the semi-infinite model in terms of accuracy to recover ${\mu }_{a_2}$ and ${\mu }_{s_1}^{\prime }$ even at large SDS (Figure 5 and Figure A5). On average the bilayered modeling performed about 10 - 15% better than SI model, across optical coefficients (except ${\mu }_{s_2}^{\prime }$) at SDS $=1.5\mathrm{and}2$ cm. It is interesting to note that bilayered DT models could extract ${\mu }_{a_2}$ accurately, even for shorter SDS could prove to be practically useful, as shorter SDS channels usually have better SNR.

5. Conclusions

We showed that using multi-layer DT for analysis of TD reflectance in bilayered media allowed for quantitative reconstruction of optical absorption of both layers with errors well within requisite thresholds for clinical utility. 64 different head and limb models at three commonly used near infrared wavelengths were used to simulate a total of 384 MC reflectance signals used for our analysis. For the top layer, absorption coefficient (${\mu }_{a_1}$) was reconstructed with mean errors $<0.03{\mathrm{cm}}^{-1}$ at SDS $=1.5\mathrm{and}2$ cm (rising to $\gtrsim 0.1{\mathrm{cm}}^{-1}$ at 2.5 cm), while reduced scattering coefficient (${\mu }_{s_1}^{\prime }$) showed high accuracy ($<2{\mathrm{cm}}^{-1}$ across all SDS). Bottom layer absorption, ${\mu }_{a_2}$ on the other hand was retrieved exceptionally well ($<0.01{\mathrm{cm}}^{-1}$ error) while bottom layer scattering ${\mu }_{s_2}^{\prime }$ could not be well estimated. These translated to errors in biomarkers with top-layer $\left[THb\right]$ and ${\mathrm{SO}}_2$ values having larger errors relative to the bottom-layer, highlighting the utility of layered DT to assess depth-resolved tissue physiology. Layered DT always outperformed semi-infinite DT, across all retrieved optical coefficients but had significantly higher computational costs. Future work will investigate approaches to reduce computational time for optimization, for e.g., by using SI DT to bootstrap the optimization, or to exclude ${\mu }_{s_2}^{\prime }$ reconstructions. Approaches that reduce dimensionality of the inverse problem could potentially accelerate reconstructions without sacrificing accuracy.