Broadband diffuse optical spectroscopy of two-layered scattering media containing oxyhemoglobin, deoxyhemoglobin, water, and lipids

We investigated the relationship between chromophore concentrations in two-layered scattering media and the apparent chromophore concentrations measured with broadband optical spectroscopy in conjunction with commonly used homogeneous medium inverse models. We used diffusion theory togenerate optical data from a two-layered distribution of relevant tissue absorbers, namely, oxyhemoglobin, deoxyhemoglobin, water, and lipids, with a top-layer thickness in the range 1–15 mm. The generated data consisted of broadband continuous-wave (CW) diffuse reflectance in the wavelength range 650–1024 nm, and frequency-domain (FD) diffuse reflectance at 690 and 830 nm; two source-detector distances of 25 and 35 mm were used to simulate a dual-slope technique. The data were inverted using diffusion theory for a semi-infinite homogeneous medium to generate reduced scattering coefficients at 690 and 830 nm (from FD data) and effective absorption spectra in the range 650–1024 nm (from CW data). The absorption spectra were then converted into effective total concentration and oxygen saturation of hemoglobin, as well as water and lipid concentrations. For absolute values, it was found that the effective hemoglobin parameters are typically representative of the bottom layer, whereas water and lipid represent some average of the respective concentrations in the two layers. For concentration changes, lipid showed a significant cross-talk with other absorber concentrations, thus indicating that lipid dynamics obtained in these conditions may not be reliable. These systematic simulations of broadband spectroscopy of two-layered media provide guidance on how to interpret effective optical properties measured with similar instrumental setups under the assumption of medium homogeneity.


Introduction
Biomedical diffuse optics typically utilizes a simplified model to analyze data collected from a complex medium. That complex medium being biological tissue, which has a heterogeneous distribution of optical properties, and that simplified model often being the semi-infinite homogeneous model. 1 Therefore, recovered optical properties are said to be effective homogeneous properties. This meaning that a homogeneous medium with these properties would replicate the behavior of the data. When the measurement of a single region is sought, such as the brain, 2 muscle, 3 or breast lesions, 4 these heterogeneous effects on homogeneous recovered data are considered partial-volume confounds. 5 Therefore, an understanding of how heterogeneous properties are related to effective homogeneous properties may help with data interpretation.
Simulations of optically diffuse media offer a method for generating forward model data where the actual optical properties within the medium are known and may be varied. This allows one to generate forward data from a complex medium and then invert the data using a homogeneous model as one would with experimental data, to compare the actual simulated medium to the recovered values. Highly complex heterogeneous media may be simulated with Monte Carlo methods. These methods may simulate voxelized or meshed media with computationally efficient parallel computations of photon propagation. 6 Monte Carlo simulations have the advantage of being able to simulate complex media but these methods do not match analytical expressions in terms of speed. Due to how fast analytical expressions may be evaluated, many media of differing properties can be simulated in a short computation time by repeating the analytical evaluation. However, the media simulated by analytical methods must be simpler compared to Monte Carlo. It may be said that the simplest heterogeneous medium is the two-layer medium for which analytical expressions for the diffuse reflectance exist. 7 In this work, analytical expressions from diffusion theory for the diffuse reflectance in a two-layer medium are used to generate forward data and expressions for a semi-infinite homogeneous medium are used to invert data.
The optical parameters describing these media can be split into two categories, those affecting the (μ a ) and those affecting the μ s ′ . 8 Absorption parameters are related to chromophore concentration. Considering biological tissue in the Near-InfraRed (NIR) range, the primary chromophores are Oxy-hemoglobin (O), Deoxy-hemoglobin (D), Water (W), and Lipid (L). In this work this notation of one letter symbols is used to represent concentrations of these chromophores. Scattering on the other hand can be described by a power law containing two parameters, an amplitude a and the scattering power law exponent (b).
NIR Spectroscopy (NIRS), whether Continuous-Wave (CW), Frequency-Domain (FD), or Time-Domain (TD), seeks to measure μ a and μ s ′, or combinations of the two, in the NIR range with the ultimate goal of recovering chromophore concentrations, scattering parameters, or both. These measurements are conducted at either discrete or broadband NIR optical wavelengths (λs), with a general trade-off of temporal information for spectral information. 9 Broadband spectroscopy creates a overdetermined system when recovering chromophores, allowing for better understanding of how the chromophore model fits the data. 11 When considering broadband absorption, the O and D parameters can be reexpressed as Total-hemoglobin (T) and oxygen Saturation (S), where T = O + D and S = O/T, which have more spectral meaning b as amplitude and shape, respectively. This work is an extension of a previous investigation utilizing Dual-Slope (DS) broadband diffuse optics on human muscle tissue. 11 This previous work observed the aforementioned complex interaction between a heterogeneous tissue and the recovered parameters, and attempted to apply a two-layer model to the measured data. As an extension of that work, this work utilizes a systematic set of simulations to obtain a deeper understanding of how properties in a two-layer medium affect the recovered values for the instruments and methods in the previous work. 11 Despite the results herein being specific to the simulated instrument and recovery methods, this methodology of careful examination of a two-layer biological medium may provide broader guidance in the field of biomedical diffuse optics.

Simulated instrumental conditions
A specific measurement system was simulated when generating forward data. Thus, this measurement system or instrument was not used in the real-world for this work. The real-world version of the simulated instrument was implemented 10 and utilized for tissue measurements 11 in the past. Only instrument parameters pertinent to the simulations are described here.
The simulated instrument was comprised of both a FD and CW portion. The FD portion measured complex reflectance phasors R s with a modulation frequency of 140.625 MHz at two discrete λs of 690 and 830 nm. The CW portion measured real reflectances (Rs) at a broadband range of wavelengths (650-1024 nm spaced by 0.5 nm). Both the FD and CW reflectances were simulated at two source-detector distances (ρs) of 25 and 35 mm. Therefore, the simulated instrument outputs: four FD reflectance phasors (R(λ, ρ); two λs and two ρs) and 1498 CW reflectances (R(λ, ρ); 749 λs and two ρs).
a The scattering amplitude is expressed as the value of the reduced scattering coefficient μ s ′ at a single optical wavelength (λ).
The real-life version of the instrument relied on a DS optode arrangement to eliminate the need for calibration 12 and provide preferentially deep sensitivity. 13,14 The calibration consideration is not pertinent for these simulations, however the expected preferentially deep sensitivity of DS does apply to the simulated measurements.

Two-layer forward model
The forward model used expressions derived from diffusion theory to simulated a two-layer medium and retrieve diffuse reflectance both R and R. The input of the model described chromophore concentrations and reduced scattering parameters in the two layers as well as the thickness, these 13 parameters are as follows: top layer Total-hemoglobin (T Top ) The first step of the model was to calculate the μ a spectrum c for each layer from the absorption parameters: T Top , T Bottom , S Top , S Bottom , W Top , W Bottom , L Top , and L Bottom .
These layer-dependent absorption spectra were calculated as a linear combination of known extinction coefficients (ϵs) [15][16][17] weighted by chromophore concentrations (Cs): where the ϵs are contained within matrix of extinction coefficients (E), C is the vector of chromophore concentrations, O = ST, and D = (1 − S)T.
Next, the μ s ′ spectrum was calculated for each layer using a power law: for each layer using the scattering parameters: The analytical two-layer model took the μ a and μ s ′ spectra for each layer and z Top as input with a final unmentioned parameter of index of refraction (n) which was assumed to be fixed at 1.4 for both layers. This analytical two-layer model has been previously derived 7 and utilized. 11 Here, the retrieval of either R or R is simulated given the above input parameters for each ρ and λ. The expression for the two-layer diffuse reflectance can be found in Appendix A.1 (Eq. (A.23)).

Semi-infinite homogeneous inverse model
Recovery of absolute optical properties was achieved using an iterative linear fit of a linearized reflectance versus ρ. Expressions for this linearized reflectance were derived from diffusion theory for the semi-infinite homogeneous medium. This method was previously described 10,11 and is explained in Appendix A.2. Using this method, R recovered complex effective attenuation coefficient μ eff , from which μ a and μ s ′ were obtained: μ a λ FD = ℜ 2 μ eff λ FD − ℑ 2 μ eff λ FD 3μ t ′ λ FD , (5) μ s ′ λ FD = μ t ′ λ FD − μ a λ FD , (6) using the parameters of speed of light in vacuum (c), n, angular modulation frequency (ω), and total reduced attenuation coefficient μ t ′ . The FD subscript signifies that the property was recovered for two λs.
A continuum of extrapolated μ s ′ was then be calculated across all CW λs using this recovered b and μ s ′ at 830 nm recovered by FD by inputting these parameters into Eq. (3).
Next, the simulated CW R was used to recover broadband effective attenuation coefficient (μ eff ) using the aforementioned iterative recovery methods 10,11 (Appendix A.2). This was combined with the extrapolated μ s ′ d to give a broadband measurement of μ a : where the CW subscript signifies that a value is calculated for each of the 749 CW λs.
Finally, the recovered CW μ a is converted to effective homogeneous e recovered T, S, W, and L using the MATLAB (MathWorks, Natick, MA, USA) mldivide function to invert Eq.
(1). Since the E is nonsquare, MATLAB's qr decomposition is used: where Q 1 and Q 2 have orthogonal columns and R 1 is an upper triangular matrix. From this decomposition Eq. (1) can be inverted on the recovered μ a spectrum: giving the effective homogeneous recovered chromophores.
Therefore, the output of this inverse model necessary to describe an equivalent homogeneous medium is two scattering parameters and four absorption parameters. The scattering parameters are the effective homogeneous recovered μ s ′ at one λ f and the effective homogeneous recovered b. The absorption parameters are the effective homogeneous recovered T, S, W, and L.

Two-layer sensitivity
Sensitivities to changes in absorption parameters in each layer of a two-layer medium were also considered. These sensitivities were calculated as numerical derivatives, where an independent variable g was changed by a small amount h and the dependent variable's i change was retrieved. The ratio of the change in the dependent variable divided by the independent variable was considered the sensitivity of the dependent variable to a change in the independent variable. In all cases μ s ′ recovered from FD was fixed at the baseline value when recovering the small change in the dependent variable. Thus only changes in absorption parameters were considered, and importantly cross-talk between absorption parameters and scattering parameters was ignored. In practice this would be equivalent to making only one FD measurement at the beginning of an experiment and then measuring absorption dynamics with only CW, therefore assuming constant scattering during an experiment.
The aforementioned independent absorption variables and inputs to the forward model (Sec. 2.2) were T Top , T Bottom , S Top , S Bottom , W Top , W Bottom , L Top , and L Bottom . The dependent absorption variables were the recovered T, S, W, and L. Sensitivities could take the form of any pair between these independent and dependent variables. These sensitivities were e The effective homogeneous properties recovered are also dependent on the specific instrument simulated, particularly the sourcedetector distance (ρs) and optical wavelengths (λs).
f This reduced scattering coefficient μ s ′ at a single optical wavelength (λ) is the scattering amplitude, and the single λ used is 830 nm. g Inputs to the two-layer forward model are considered independent variables. h The small amount used to calculate the numerical derivative was the baseline value of the parameter multiplied by 10 −6 . i Dependent variables are the outputs of the homogeneous inverse model. split into two categories, first co-sensitivities which are between like parameters j and second cross-sensitivities which are between differing parameters. k Taking the example of the co-sensitivity between T Top and T (δT/δT top ), this should be interpreted as how much a change in the top-layer T will affect the effective homogeneous recovered T. These sensitivity values also are affected by the baseline properties, therefore this dependence was also investigated by recalculating sensitivities for various values of baseline two-layer absorption parameters.

Recovery of absolute values
The results in this section focus on how actual properties in a two-layer medium are related to the absolute recovered effective homogeneous properties. This was done by simulating the R and R for a two layer medium (Sec. 2.2) with the particular parameters of the simulated instrument (Sec. 2.1), then inverting these data using a semi-infinite homogeneous model (Sec. 2.3). Table 1 shows the baseline absorption model and the recovered baseline absorption properties for said model. Table 2 shows scattering values for this same baseline model. Unsurprisingly, in this case the recovered parameters take values between the actual values of the two layers. However, the dominant layer is different for different parameters, for example recovered T is closer to the actual bottom-layer value but recovered W is approximately the average of the top and bottom layers. Finally, Fig. 7

Sensitivity
The results in this section focus on how changes in absorption parameters in a two-layer medium are related to the recovered effective homogeneous absorption changes. o This was done by calculating sensitivities according to Sec. 2.4. Tables 1 and 2 as baseline, co-sensitivities p were calculated (Table 3). From this it can be seen that recovered dynamic changes in T, L, and especially S are more sensitive to bottom-layer dynamics of the like chromophore. Conversely, recovered W dynamics are more sensitive to changes in the top-layer. Here it is emphasized that these values are only strictly actual for the particular baseline model chosen, and this consideration will be addressed in the following section. Table 4 shows the full set of sensitivities displayed as fractions for the baseline model. In this case the co-sensitivities (Table 3) can be found in the diagonal of the upper and lower sub-tables and all off-diagonal elements are cross-sensitivities q which can be thought of as cross-talk between parameters. Examining the cross-sensitivities in Table 4, it appears that there is little cross-talk between parameters, with the only parameter displaying notable cross-talking being the recovered changes in L. In the case of measuring recovered L dynamics, its value may be affected by changes in W Top , S Bottom , and W Bottom . Again, all values are valid only for the particular baseline model. However, since little cross-sensitivity was observed r overall, the following section will neglect cross-sensitivity for brevity. p Sensitivities between like parameters are dubbed co-sensitivities, between recovered T and T Top for example. q Sensitivities between dis-similar parameters are cross-sensitivities, for example between recovered T and W Top . r Little cross-sensitivity except in the case of recovered L changes, which typically is assumed to not have dynamics. were held at the values of the baseline model (Tables 1 and 2), simply stated, only one parameter was varied at a time.

Variation sensitivities as a result of varying model parameters-In
These results are displayed in Fig. 8, where each sub-plot contains dual y-axes and dual x-axes. To which y-axis a curve belongs to is designated by the line type that is specified in square brackets of the axis title. All solid lines correspond to T sensitivities, dashed lines to S sensitivities, dotted lines to W sensitivities, and dash-dotted lines to L sensitivities. Which x-axis a curve belongs to is signified by color, which is the same as the color of the variable in the axis title. All yellow lines show the sensitivity while varying the absolute value of S, blue lines while varying the absolute value of W, green lines while varying the absolute value of L, and purple lines while varying the absolute value of T. Each row of subplots examines differing combinations of which layer the sensitivity is to and which layer the absolute value is being varied.
Despite the plethora of information in Fig. 8 the story in most subplots is quite simple. This is because the absolute value of a given parameter does not affect many of the sensitivities in a given subplot. For example, in Fig. 8(a) the sensitivities of recovered T or S to top-layer changes as a function of actual absolute top-layer properties is plotted. In this case it is immediately obvious that only the purple curves significantly change. This means that the sensitivities of recovered T or S to top-layer changes are not affected by the absolute value of S Top , W Top , or L Top . However, since the purple line does change, the T sensitivity to the top-layer decreases and the S sensitivity to the top-layer increases each with greater T Top .
Examining Fig. 8 as a whole, some general trends are observed. In the left column the curves show that T and S sensitivities are not affected by the actual absolute values of W or L, but are somewhat affected by the actual value of S and are strongly affected by the actual value of T. The right column concerns that sensitivities of W and L which appear to be more heavily affected by the actual absolute value of all parameters but to a lesser extent on S. The actual absolute value of T has the most effect on these sensitivities, closely followed by the actual absolute value of W. Finally, the actual absolute value of L does not affect the W sensitivity, but the opposite is not true.
Again, it is noted that all these results are specific to the model in question. Despite the parameters in this experiment being varied, they were not covaried, scattering or thickness were not varied at all, and only the specific measurement method in Sec. 2.1 was considered.

Discussion
This work sought to investigate how actual properties from a two-layer model are related to the properties retrieved from methods that utilize a DS instrument and semiinfinite homogeneous model at their core. Arguably the two-layer model is the simplest heterogeneous model one can have, and yet its interactions with a homogeneous inversion model are still complex. Examining this simple pair of mismatched forward and inverse models may help provide insight into the meaning of retrieved optical properties from tissues, which are in fact much more complex and heterogeneous than even a two-layer model.
First, the case of recovery of effective absolute homogeneous optical properties was considered (Sec. 3.1). General trends observed in those simulations are:

2.
Recovered S is dominated by S Bottom (Table 1 and Fig. 2).

3.
Recovered W and L are not dominated by either the top-or bottom-layer value alone (Table 1 and Figs. 3 and 4).

4.
Actual scattering properties have little effect on recovered T and S but do effect recovered W and L (Figs. 5(a), 5(c), 6(a) and 6(c)).

5.
Recovered b can take values outside of the range between actual b Top and b Bottom ( Fig. 6(b)).
As a whole these results provide some valuable insights. For example, a goal in biomedical diffuse optics is to measure the chromophore concentrations of deep tissues. Points (1) and (2) suggest that for this model and measurement method the goal is feasible, at least for a similar z Top to the ones studied. Additionally, point (5) may provide some insight into apparent scattering powers that are outside of the range of what is physically possible. s Finally, these results show that each individual recovered property is not representative of the same region as the others. For example, point (6) suggests that recovered μ s ′ is dominated by the top layer but point (1) suggests that the opposite is true for T. Since the recovered μ s ′ resulted from the FD measurement and μ a from CW, this may be evidence that FD exhibits more superficial sensitivity compared to CW. This is consistent with the idea that at increasing modulation frequencies sensitivity to optical property changes becomes more superficial. 13 However, this effect only occurs at frequencies above 500 MHz, therefore it is likely that this is not the effect here. Thus, this observed effect is likely simply due to the μ s ′ measurement being dominated by superficial layers.
The results for sensitivity (Sec. 3.2) tell a similar story to those of absolute properties, but with the additional focus on dynamic changes in absorption while assuming scattering is constant. Again, a list of general observations follows:

7.
Recovered dynamics in absorption parameters are significantly affected by each layer except for S whose dynamics are dominated by the bottom layer (Table 4).

8.
There is little cross-talk in dynamics of different parameters except for L which has significant cross-talk with W Top , S Bottom , and W Bottom (Table 4).

9.
The actual absolute value of T or W affects the sensitivity of almost every parameter (Fig. 8).

10.
The actual absolute value of S has little effect, for the most part, on the sensitivity of almost every parameter except the sensitivity recovered S to S Bottom (Figs. 8(c) and 8(g)).

11.
The actual absolute value of L only has significant effect on L sensitivity (Figs.

8(b), 8(d), 8(f) and 8(h)).
As may be expected in most cases, sensitivities are affected by the absolute value of multiple parameters within a two-layer medium. However, point (8) may be a significant result. This being that L dynamics will have significant cross-talk with W and S dynamics due to the high cross-sensitivity (Table 4). Therefore, caution is recommended when allowing L to dynamically change in analysis, especially for the wavelength range examined here. This confound may not be a significant issue since L is not typically expected to dynamically change, thus it is reasonable and, likely safer due to the cross-talk, to assume no L changes from baseline during one's analysis.
These results should be considered when interpreting certain recovered properties or dynamic changes in a heterogeneous medium. Furthermore, since the particular results are dependent on instrumental setup, it may be helpful for one to repeat such simulations for one's own instrument, thus aiding in the interpretation of data from said instrument. The two-layer medium is much simpler than the complex heterogeneous one within biological tissue, but the model does provide insight into how different instruments and methods behave when measuring heterogeneity.

Conclusion
Within this article various simulations were carried out to evaluate a particular simulated instrument (Sec. 2.1) and optical property inversion model (Sec. 2.3). The simulations focused on using a two-layer forward data (Sec. 2.2) which was considered that simplest type of heterogeneous media. In the case of absolute optical property recovery (Sec. 3.1) the significant result was that bottom-layer hemoglobin recovery is feasible, but scattering measurements are dominated by a top-layer. When examining dynamic changes in terms of sensitivities (Sec. 3.2) the main result was to be weary when measuring lipid dynamics due to potential cross-talk.
Overall, these simulation, analysis, and evaluation methods provide valuable guidance when interpreting data measured from a heterogeneous medium. While these results are specific to the model within this work, these types of simulations may be repeated or applied for other instruments and model media. Thus, this work is not only specific to this model but, also may be used as guidance on how other systems may be modeled to determine how heterogeneous media affect said system's measurements.
scattering coefficient μ s ′ , and the top layer thickness (z Top ), the dependent variables in these expressions are shown as functions of them. Other variables are in fact independent but are not included in this notation since they are not varied within this work. t Furthermore, for brevity, this function notation is only used on the left-hand-side of the following expressions.

A.1. Expression for two-layer reflectance
To express the two-layer (cylinder) 7  The approximate extrapolated boundary height (z b ): where n is the index of refraction of the medium (assumed homogeneous here) and A is the index of refraction mismatch parameter. 18 Next, various coefficients are defined as a function of the kth zero of the 0th-order Bessel function of the first kind. First, the optical attenuation parameter a k : where μ a is the absorption coefficient, j 0,k the kth zero of the 0th-order Bessel function of the first kind, B is the radius of the two-layer cylinder (for which the source is in the center; assumed to be 150 mm for this work), and i is the imaginary unit. Then, the source-detector distance parameter (Q k ): t Examples of un-varied and de-emphasized independent variables are index of refraction (n) and angular modulation frequency (ω). where J ν is the νth-order Bessel function of the first kind and ρ is the source-detector distance.
Now, the reflectance for the case when the isotropic source is in the first layer R 1 can be expressed:  Then combinations of these coefficients are defined.

A.2.2. Inversion of semi-infinite reflectance
Given that the measurement method focuses on measuring complex R or real (R) reflectance versus source-detector distance (ρ), x Eq. (A.24) is rearranged to find a linear relationship between the measured reflectance and the effective attenuation coefficient μ eff across ρ: ln (4πR)/ C 1, m − 1 + C 2, m − 1 e μ eff, m − 1 r 1, m − 1 − r 2, m − 1 = − μ eff, m r 1, m − 1 . (A.32) Here the notation of dependent variables being a function of independent is dropped. With an initial guess for μ eff the slope of the left-hand-side of Eq. (A.32) versus r 1 can be fit y to find an updated estimate of μ eff . This procedure is repeated iteratively, where m is the iteration number, until μ eff converges. The initial guess was chosen using previously described methods which make further assumptions to achieve simpler expressions. z ,19 For CW data the same procedure is used, only with the complex variables replaced with the real ones and the real effective attenuation coefficient (μ eff ) as the output instead of the complex one.          Tables 1 and 2, and top layer thickness was 5 mm. Black lines show the co-sensitivities of the baseline model (Table 3). X-axes: Curves belong to the x-axis of like color. Purple curves correspond to the bottom x-axes and show the variation of the actual absolute value of T. Yellow, blue, and green    Co-sensitivities for a two-layer medium with the baseline optical properties of Tables 1 and 2. δT/δT lay δS/δS lay δW/δW lay δL/δL lay All sensitivities for a two-layer medium with the baseline optical properties of Tables 1 and 2.  Table 3, off-diagonal elements are cross-sensitivities. Note 3: All values are fractions none are represented in percent.