2.1.

Monte-Carlo-Based Multi-Layer Inverse Models

Monte Carlo (MC) simulations were performed using Monte Carlo eXtreme²² in the MATLAB environment (Mathworks, Natick, MA, United States) on a standalone workstation (Ryzen 9 3800x, Corsair 64 GB RAM, and Nvidia RTX 2070). LUTs were generated as inverse models based on MC results. MC-based multi-layer DCS LUTs have been described in prior works^19,23 and multi-layer FD-NIRS LUTs are an extension of single-layer LUTs that our group has previously demonstrated.²⁴ Briefly, a semi-infinite geometry was defined with three layers: the skin, adipose, and muscle. A source and detector were positioned on top of the skin layer with a source–detector separation (SDS) of 25 mm [Fig. 1(b)]. Paired single-layer models were also developed for comparison [Fig. 1(a)]. The upper two layers (skin and adipose) of the multi-layer simulations were assigned static optical properties, and the skin layer was assigned a static thickness. Assumed values were informed by prior literature and are shown in Table 1. The skin optical properties were determined from a prior study that used spatial frequency domain imaging to measure skin optical properties and used the Fitzpatrick skin type scale to classify the subject’s skin tones.²⁵ Subjects with a Fitzpatrick score of less than 16 were categorized as having light skin tone, subjects with a score of greater than 26 had dark skin tone, and subjects that scored between 16 and 26 had medium skin tone.²⁵ In addition, prior literature has shown that upper subcutaneous tissue (skin + adipose) has a ${\mathrm{BF}}_{\mathrm{i}}$ baseline value of $0.6{\mathrm{mm}}^2/\mathrm{s}$. For each wavelength, the bottom layer was first assigned as ${\mu }_a=0{\mathrm{mm}}^{-1}$ and 10 simulations were run for each ${\mu }_s^{\prime }$ value with a defined step size of $0.1{\mathrm{mm}}^{-1}$ over a range of 0.1 to $1{\mathrm{mm}}^{-1}$. Photons were launched from the source position and recorded at the detector location. For each detected photon, the total pathlength and the total dimensionless momentum transfer were recorded.²³ For DCS, both the total pathlength and total momentum transfer information were fed into the electric temporal field autocorrelation function equation²³ alongside ${\mu }_a$ and ${\mathrm{BF}}_{\mathrm{i}}$ to obtain the temporal field autocorrelation function (${\mathrm{G}}_1$). The resulting ${\mathrm{G}}_1$ and an instrument-dependent coherence factor (0.45 for phantom measurements and $\sim 0.35$ for in vivo measurements) were then input into the Siegert equation³³ to obtain the normalized temporal autocorrelation function (${\mathrm{g}}_2$). In the case of FD-NIRS, the pathlength and ${\mu }_a$ were input into Beer’s Law to calculate absorption weights.^24,34 The resulting photon absorption weights were then binned based on photon time of flights to create photon time-of-flight distributions. A Fourier transform was applied to the photon time-of-flight distributions function to generate a complex frequency-domain reflectance measurement at a given modulation frequency. This procedure was done for each MC simulation that had a distinct bottom layer ${\mu }_s^{\prime }$. Linear interpolation was performed on both the complex FD reflectance and $g_2$ from each distinct MC simulation to create high-density LUTs, the ${\mu }_s^{\prime }$ interpolation step size was 0.001 for the complex FD reflectance and 0.05 for the $g_2$. Each individual simulation took $\sim 5\min$ and a single LUT took on average 15 min to create. Fitting 16 min of data collected at a sampling frequency of 0.5 Hz took $\sim 25\mathrm{s}$ for either FD-NIRS or DCS dataset.

Fig. 1

(a) Monte Carlo geometry used in single-layer models. (b) Monte Carlo geometry used in multi-layer models. (c) Stackable phantoms used in phantom measurement. From left to right: light skin, medium skin, dark skin, adipose tissue, and muscle. SDS, source–detector separation; ${\mu }_a$, absorption coefficient; ${\mu }_s^{\prime }$, reduced scattering coefficient; ${\mathrm{BF}}_{\mathrm{i}}$, blood flow index.

Table 1

Properties of the skin, adipose, and muscle layers and the corresponding references used in the MC simulations and g2 calculations.

2.2.

Stackable Multi-Layer Silicone Phantoms

Stackable silicone phantoms were fabricated to validate the multilayer LUTs [Fig. 1(c)]. Thin silicone phantoms were created to mimic skin and adipose layers. These phantoms consisted of clear silicone (Eager Polymers, P4 Silicone), nigrosin (Sigma-Aldrich, N4754), and titanium dioxide (Sigma-Aldrich, 248576) and were cast in custom molds to ensure accurate layer thicknesses (1, 2, and 4 mm). The custom mold consisted of two clear acrylic pieces with one plate having an overflow hole and a 3D-printed perimeter piece with the desired thickness. A technical drawing of the mold is shown in Fig. S1 in the Supplementary Material. Thin phantoms were fabricated by pouring the silicone solution on the bottom plate and perimeter then laying the top plate, excess material flowed out from the overflow hole. The full phantom recipes are shown in Table S1 in the Supplementary Material. A matching larger bulk homogeneous phantom ($70\times 108\times 38\mathrm{mm}$) was made for each thin phantom and used to measure the optical properties using a custom FD-NIRS system. The measured optical properties are shown in Table 2. The thin phantoms were stacked to create a multi-layer phantom with the bottom layer being either a silicone phantom mimicking muscle optical properties for the FD-NIRS validation dataset or a liquid phantom for DCS validation. The liquid phantom consisted of intralipid (Fresenius Kabi, Intralipid 20%) and nigrosin. A custom DCS system was used to measure the ${\mathrm{BF}}_{\mathrm{i}}$ from the Brownian motion of the liquid phantom. The liquid phantom’s recipe is shown in Table S1 in the Supplementary Material, and the optical properties and ${\mathrm{BF}}_{\mathrm{i}}$ are shown in Table 2.

Table 2

The optical properties of silicone and intralipid phantoms. Values for g and n were assumed.

2.3.

Sensitivity Analysis

Several assumptions related to the skin and adipose layers in the multi-layer models were made based on prior literature (Table 2). These include the skin thickness as well as the skin and adipose optical properties. Although assigning fixed values to the model parameters reduced computational complexity, any discrepancies between the assumed properties and measured tissue could induce errors in extracted target tissue optical properties. A sensitivity analysis was conducted to quantify how errors in these assumptions would affect measured target tissue properties. Each assumed model parameter was perturbed one at a time in separate MC simulations to evaluate the effect on extracted target tissue ${\mu }_a$, ${\mu }_s^{\prime }$, and BFi. The magnitude of the tested model perturbations is shown in Table 3 and is based on the reported variance of baseline literature values from healthy subjects. The perturbation percentages were defined as the approximate reported standard deviations divided by the mean of each parameter. The exception was the adipose thickness where the perturbation was set to 10%. Sensitivities were determined for the case of 1-mm adipose thickness and light, medium, and dark skin tones. Sensitivity for ${\mu }_a$, ${\mu }_s^{\prime }$, and ${\mathrm{BF}}_{\mathrm{i}}$ was defined as

Table 3

Perturbation values for the model assumptions and the corresponding reference that justifies the magnitude of the perturbation.

2.4.

Healthy Volunteer Study

All measurements were conducted under an institutionally approved protocol (BU IRB 5618E), and all subjects completed a consent form prior to the experiments. A custom system that consisted of FD-NIRS (730 nm modulated at 139 MHz and 830 nm modulated at 149 MHz) and DCS (850 nm) was used to measure the muscle tissue of the sternocleidomastoid muscle during breathing exercises.³⁷ The SDS was 25 mm. The study protocol is explained in detail in a previous work.³⁷ Briefly, 17 healthy young ($26.1\pm 1.8$ years) volunteers performed a breathing exercise that consisted of a baseline, load, and recovery period while the sternocleidomastoid, a neck muscle, was monitored. The load consisted of subjects breathing through a respiratory trainer and each subject performed both a low load and a high load. A colorimeter (PCE Instruments, PCE-CSM1) was used to determine the subject’s individual typology angle (ITA).^38,39 The determined ITA was then used to classify each subject’s skin tone as either light ($n=7$), medium ($n=6$), or dark ($n=4$) according to the thresholds determined by literature. An $\mathrm{ITA}>41\mathrm{deg}$ was classified as light skin tone, an $\mathrm{ITA}<10\mathrm{deg}$ was classified as dark skin tone, and an ITA between 10 deg and 41 deg was classified as medium skin tone.^38,39 In addition, an ultrasound (GE, Vscan Extend) with a broad-bandwidth linear array (3.4 to 8 MHz) was used to determine the adipose thickness ($2.4\pm 0.7\mathrm{mm}$). The skin was assumed to be 1-mm thick for all subjects.

Recovered target tissue ${\mu }_a$ and ${\mathrm{BF}}_{\mathrm{i}}$ values were processed to obtain chromophore values and the metabolic rate of oxygenation (${\mathrm{MRO}}_2$). Oxygenated hemoglobin plus myoglobin (oxy [Hb + Mb]) and deoxygenated hemoglobin plus myoglobin (deoxy [Hb + Mb]) concentrations were calculated using Beer’s Law from ${\mu }_a$ values at 730 and 830 nm with an assumption of 20% lipid fraction and 62.5% water fraction.^40–43 Total hemoglobin plus myoglobin (total [Hb + Mb]) is the addition of oxy [Hb + Mb] and deoxy [Hb + Mb]. Tissue saturation (${\mathrm{S}}_{\mathrm{t}}{\mathrm{O}}_2$) is defined as

The ${\mathrm{MRO}}_2$ is defined by the following equation:⁹

The Wilcoxon signed-rank test was used to evaluate the statistical significance of values derived from the single-layer versus multi-layer LUTs. The test was performed for all muscle parameters at both baseline and during the peak of the perturbation. Baseline values were calculated by finding the mean value of the first 50 s of the measurement. For the perturbation value, the largest absolute maximum difference from the baseline measured in either the low or the high load was selected as the perturbation value. Additional subanalyses were conducted to disentangle the effects of the skin tone and adipose thickness. The difference between the single layer and multi-layer for both baseline and perturbation were compared with both adipose thickness and skin tone. Robust linear regression analysis was conducted to show the strength of the relationship that each parameter difference had with skin tone and adipose thickness, respectively.

3.1.

Stackable Multi-Layer Silicone Phantoms

The multi-layer LUTs greatly improved the accuracy of target tissue optical properties in stackable silicone phantoms across skin tone and adipose thicknesses, as shown in Fig. 2. Here, the recovered ${\mu }_a$ and ${\mu }_s^{\prime }$ at 830 nm, as well as the ${\mathrm{BF}}_{\mathrm{i}}$ are shown for both the simple single-layer LUTs (single-layer experimental values) and multi-layer LUTs (multi-layer experimental values). In addition, the known target tissue properties (target muscle values) are indicated as well as the expected extractions for the single-layer case from MC simulation (single-layer theoretical values). In general, when the mimicked skin tone and adipose thickness were not taken into account (i.e., when a single-layer LUT was used), there was an underestimation of both ${\mu }_a$ and ${\mathrm{BF}}_i$. This underestimation increased as a function of adipose thickness, with errors reaching $-35\pm 4\%$ and $-32\pm 0\%$ at 6-mm adipose thickness, respectively. In addition, ${\mu }_s^{\prime }$ was overestimated by up to $121\pm 2\%$ at 6-mm adipose thickness. When the multi-layer model was used, the extracted properties closely matched the known target values across all the skin tones and adipose thickness. Similar trends were seen in the recovered ${\mu }_a$ and ${\mu }_s^{\prime }$ at 730 nm, as shown in Fig. S2 in the Supplementary Material.

Fig. 2

Recovered target ${\mu }_a$ and ${\mu }_s^{\prime }$ at 830 nm and ${\mathrm{BF}}_{\mathrm{i}}$ from multi-layer phantom measurements. When measurements were processed with multi-layer LUTs, the recovered values (teal line) closely matched those of the known ground truth of the bottom layer (black line). In addition, the recovered values calculated via a single-layer LUT (red line) matched the expected extractions from simulations (yellow line). An experimental measurement error occurred for ${\mathrm{BF}}_{\mathrm{i}}$ for the dark skin phantom at 5-mm adipose thickness, and those data are not shown. ${\mu }_a$, absorption coefficient; ${\mu }_s^{\prime }$, reduced scattering coefficient; ${\mathrm{BF}}_{\mathrm{i}}$, blood flow index.

3.2.

Sensitivity Analysis

The sensitivity analysis results are shown in Fig. 3 for ${\mu }_a$ and ${\mu }_s^{\prime }$ at 830 nm and ${\mathrm{BF}}_{\mathrm{i}}$. The sensitivities are presented as %/% (i.e., unitless) so that they can be compared across parameters. The model was most sensitive to the assumption of skin thickness, which had $\sim 1\%/\%$ sensitivity for ${\mu }_a$ and $-1\%/\%$ for ${\mathrm{BF}}_{\mathrm{i}}$. This means that a 1% mismatch in the skin thickness leads to a $\pm 1\%$ error in ${\mu }_a$ and ${\mathrm{BF}}_{\mathrm{i}}$. Prior literature suggests that skin thickness can vary by $\sim 10\%$ across different anatomic locations,²⁶ which could lead to similar magnitude errors in ${\mu }_a$ and ${\mathrm{BF}}_{\mathrm{i}}$. Although the sensitivity to skin ${\mu }_a$ was relatively low, approximately $\pm 0.2\%/\%$, this parameter is known to vary by as much as 75% across a measured population,²⁵ which would lead to, for example, a $\pm 15\%$ error in target ${\mu }_a$. In general, ${\mu }_s^{\prime }$ was relatively insensitive to mismatches in model assumptions; the largest sensitivity was only $\sim 0.2\%/\%$ which arose from a perturbation in either skin or adipose ${\mu }_s^{\prime }$. Similar trends were seen for both ${\mu }_a$ and ${\mu }_s^{\prime }$ at 730 nm (Fig. S3 in the Supplementary Material). Although sensitivities had minimal dependence on the skin tone for the 830 nm results, this was not the case for 730 nm, in which the dark skin tone on average was less sensitive to model assumption compared with medium and light skin tones, especially for the ${\mu }_a$ results.

Fig. 3

Sensitivity of the multi-layer model to different model assumptions for extractions of ${\mu }_a$ and ${\mu }_s^{\prime }$ at 830 nm, and ${\mathrm{BF}}_{\mathrm{i}}$ model was most sensitive to skin thickness for extractions of both ${\mu }_a$ and ${\mathrm{BF}}_{\mathrm{i}}$. No substantial differences were observed across skin tones. ${\mu }_a$, absorption coefficient; ${\mu }_s^{\prime }$, reduced scattering coefficient; ${\mathrm{BF}}_{\mathrm{i}}$, blood flow index.

3.3.

Healthy Volunteer Study

Example time traces for both single-layer and multi-layer results from a single subject are shown in Fig. 4, and others can be seen in Figs. S4 and S5 in the Supplementary Material. The time traces shown are those of ${\mathrm{MRO}}_2$, and the multi-layer derived values are more biologically plausible than the single-layer values. The comparison between in vivo measurements using the single-layer and multi-layer LUTs is shown in Fig. 5. There were statistical significances between extractions from the two models for both baseline and perturbation values. Notably, oxy [Hb + Mb], total [Hb + Mb], ${\mathrm{S}}_{\mathrm{t}}{\mathrm{O}}_2$, and ${\mathrm{BF}}_{\mathrm{i}}$ were underestimated with the single-layer model at baseline. Oxy [Hb + Mb] values were underestimated by $78\pm 33\mu \mathrm{mol}$, total [Hb + Mb] by $69\pm 33\mu \mathrm{mol}$, ${\mathrm{S}}_{\mathrm{t}}{\mathrm{O}}_2$ by $18\pm 9$ %pt, and ${\mathrm{BF}}_{\mathrm{i}}$ by $7\pm 8\times {10}^{-6}{\mathrm{mm}}^2/\mathrm{s}$. The perturbation magnitudes of oxy [Hb + Mb], total [Hb + Mb], ${\mathrm{BF}}_{\mathrm{i}}$, and ${\mathrm{MRO}}_2$ were also underestimated with the single-layer model. Across most of the six muscle-based parameters for both the baseline and perturbation, the values from multi-layer LUTs had a wider distribution than the dataset derived from a single-layer model.

Fig. 4

Example mean time traces for ${\mathrm{MRO}}_2$ from one subject with an ITA value of $-17.8\mathrm{deg}$ (dark skin tone) and an adipose thickness of 1.5 mm. The values derived from single-layer LUTs are denoted by the red line, and the multi-layer LUTS values are denoted by the teal line. The vertical dashed lines indicate the start ($t=60\mathrm{s}$) and stop ($t=120\mathrm{s}$) of the loading phase. The ${\mathrm{MRO}}_2$ values were substantially underestimated when the simpler single-layer model was used.

Fig. 5

Comparison between single-layer and multi-layer values for both the baseline and the absolute difference of baseline values and the magnitude of perturbations during a breathing exercise. Baseline values overall were increased when a multi-layer LUT was used to recover the values except for deoxy [Hb + Mb] and ${\mathrm{MRO}}_2$. The perturbations from multi-layer LUT were larger than those from single-layer LUT except for deoxy [Hb + Mb] and ${\mathrm{S}}_{\mathrm{t}}{\mathrm{O}}_2$. Each box chart displays the following information: the median which is denoted by the notch, the lower and upper quartiles which are denoted by the shaded area, any outliers which are denoted by circles, and the minimum and maximum values that are not outliers which are denoted by the whiskers. * $p<0.05$, ** $p<0.01$, *** $p<0.001$, **** $p<0.0001$, and ***** $p<0.00001$.

The relationship between skin tone and adipose thickness and the difference in baseline values between the two LUTs are shown in Fig. 6. Adipose thickness had weak relationships for all six parameters, with all $R^2<0.30$. Skin tone was strongly correlated to both differences in extracted deoxy [Hb + Mb] and ${\mathrm{S}}_{\mathrm{t}}{\mathrm{O}}_2$, both having an $R^2>0.50$. Skin tone and deoxy [Hb + Mb] had the strongest relationship and showed that the single-layer LUT overestimated deoxy [Hb + Mb] for dark skin tone subjects while underestimating the lighter skin tone subjects. ${\mathrm{S}}_{\mathrm{t}}{\mathrm{O}}_2$ was underestimated for all skin tones, with the largest underestimates occurring for the darkest skin tones. There were no $R^2>0.2$ for the differences in perturbation magnitude (Fig. S6 in the Supplementary Material).

Fig. 6

Relationship of skin tone (ITA) and adipose thickness on the differences in baseline optical parameter extractions between the models. $Y$-axis values show extractions from the multi-layer LUT subtracted from the single-layer LUT. These values can be interpreted as the extraction errors when using a single-layer LUT that does not take into account adipose thickness or skin tone. Therefore, a negative value means that the use of a simple single-layer model would cause an underestimate in that parameter compared with a multi-layer model. Robust linear regression analysis was used to determine the best-fit line and the $R^2$ value. The adipose thickness relationship with the six tissue parameters was relatively weak, $R^2<0.30$. Skin tone had stronger relationships, especially with deoxy [Hb + Mb] ($R^2=0.70$) and ${\mathrm{S}}_{\mathrm{t}}{\mathrm{O}}_2$ ($R^2=0.56$).