2.1.

Two-Layer Head Model

We modeled the head as a cylinder with radius $a$ that consists of a homogeneous extracerebral layer of thickness $\ell$ above an infinitely thick homogeneous cerebral layer [Fig. 1(a)].²³ The tissue absorption coefficient, reduced scattering coefficient, and blood flow index of the extracerebral layer are ${\mu }_{a,1}$, ${\mu }_{s,1}^{\prime }$, and $F_1$, respectively. The corresponding properties of the cerebral layer are ${\mu }_{a,2}$, ${\mu }_{s,2}^{\prime }$, and $F_2$. The index of refraction is assumed to be the same for both tissue layers. Note that we used a cylindrical geometry instead of a slab geometry because the numeric approximation of the cylindrical Green’s function solutions for diffusive light transport is more robust than the laterally infinite 2L solution.¹⁶ In the model, a point source is incident on the middle of the cylinder top, and multiple detectors are positioned on the cylinder top at different distances from the source. The distance between the source and the $i$’th detector is ${\rho }_i$. For FD-DOS, the source is radio-frequency intensity modulated light, which produces a diffuse photon density wave in the tissue oscillating at the same frequency.^1,3 At each detector position, the wave’s amplitude and phase are measured, i.e., ${\mathrm{AC}}_{\mathrm{meas}}({\rho }_i)$ and ${\theta }_{\mathrm{meas}}({\rho }_i)$. For DCS, the source is a continuous-wave (CW), long-coherence-length laser. At each detector position, the detected normalized intensity temporal autocorrelation function, $g_2^{(\mathrm{meas})}({\rho }_i,\tau )={⟨I(t)I(t+\tau )/I(t)⟩}^2$, is computed at multiple delay times, $\tau$; here $I(t)$ is the detected light intensity at time $t$, and the angular brackets, $⟨⟩$, represent time averages.

Fig. 1

Geometry used for the 2L model and the simulations. (a) Our 2L model comprised a homogeneous cylindrical extracerebral layer of thickness $\ell$ and radius $a$ (corresponding to extracerebral scalp and skull tissue) above an infinitely thick homogeneous cylindrical cerebral layer (corresponding to the brain cortex). The tissue absorption coefficient, reduced scattering coefficient, and blood flow index of the extracerebral layer are ${\mu }_{a,1}$, ${\mu }_{s,1}^{\prime }$, and $F_1$, respectively. The corresponding properties of the cerebral layer are ${\mu }_{a,2}$, ${\mu }_{s,2}^{\prime }$, and $F_2$. A point source (S1) is incident on the middle of the cylinder top, and multiple point detectors $(D_1,D_2,\dots ,D_n)$ are positioned at different distances from the source. (b) Using the NIRFASTer package, we generated synthetic data for a $10\times 10\times 10{\mathrm{cm}}^3$ 2L cube, with an extracerebral layer thickness of $\ell =1.2\mathrm{cm}$. (c) We additionally used NIRFASTer to generate synthetic data for a realistic adult head geometry, wherein the scalp and skull were combined to form a homogeneous extracerebral layer, and the CSF, white matter, and gray matter were combined to form a homogeneous cerebral layer. The source and detectors were positioned on the right side of the head, and we used the average skin-to-brain distance under the middle portion of the optical probe [i.e., the average thickness of the 2-cm long gray line in (c)], as $\ell$.

To evaluate the 2L fitting algorithm, we used simulated data for 8 FD-DOS source–detector separations (${\rho }_i=0.8$ to 4.0 cm) and 2 DCS source–detector separations (${\rho }_i=0.8$ and 2.5 cm) across a wide range of blood flows and optical properties.

2.2.

Forward Model Simulations

Our first synthetic dataset was generated using the analytical Green’s function solutions of the frequency-domain photon diffusion equation and the correlation diffusion equation for the 2L cylinder geometry (subject to the extrapolated-zero boundary condition). The frequency-domain Green’s function (i.e., ${\mathrm{\Phi }}_{2\mathrm{L}}({\rho }_i)={\mathrm{AC}}_{\mathrm{theo},2L}({\rho }_i)\exp (-i{\theta }_{\mathrm{theo},2L}({\rho }_i))$) in this geometry was derived elsewhere^16,23 and is presented in Appendix A. The corresponding correlation diffusion Green’s function (i.e., $G_1^{(\mathrm{theo},2L)}({\rho }_i,\tau )$) is also presented in Appendix A. It has the same form as the Green’s function solution of the continuous wave photon diffusion equation, with the new feature of the decay constant in the solution depending on $\tau$.¹

The frequency-domain solution depends on the tissue optical properties (i.e., ${\mu }_{a,1}$, ${\mu }_{s,1}^{\prime }$, ${\mu }_{a,2}$, and ${\mu }_{s,2}^{\prime }$), the source–detector separation (${\rho }_i$), the cylinder radius ($a$), the extracerebral layer thickness ($\ell$), the source intensity modulation frequency ($f$), and the tissue index of refraction ($n$). All FD-DOS data were generated using $f=110\mathrm{MHz}$, $a=30\mathrm{cm}$, and $n=1.4$. Note that 110 MHz is a commonly used modulation frequency in FD-DOS instrumentation (e.g., Imagent, ISS), and $a$ was sufficiently large such that the Green’s function solution at the detector positions is not affected by the cylindrical border.

We evaluated the frequency-domain Green’s function at eight different ${\rho }_i$ (see Sec. 2.1) across a wide range of optical properties for four evenly spaced $\ell$ between 1.0 and 1.6 cm (this range approximates the range of thicknesses for adult humans ^24–26). Specifically, at each $\ell$ value, the Green’s function was evaluated for 2030 different combinations of ${\mu }_{a,1}$, ${\mu }_{s,1}^{\prime }$, ${\mu }_{a,2}$, and ${\mu }_{s,2}^{\prime }$. To mimic the range of properties observed in adult humans for the 700 to 900 nm spectral range,^27–29 ${\mu }_{a,1}$ and ${\mu }_{a,2}$ were randomly selected between 0.08 and $0.18{\mathrm{cm}}^{-1}$, and ${\mu }_{s,1}^{\prime }$ and ${\mu }_{s,2}^{\prime }$ were randomly selected between 6 and $15{\mathrm{cm}}^{-1}$ subject to the constraint that the fractional difference between ${\mu }_{s,1}^{\prime }$ and ${\mu }_{s,2}^{\prime }$ was $<20\%$. Of note, this latter constraint is justified by a recent study that observed considerable variations in overall scattering across the near-infrared spectral range but small scattering differences between skin, skull, and brain tissue.²⁸

Random amplitude and phase noise derived from a Gaussian noise model with zero mean was then independently generated for each source–detector separation and each combination of optical properties. For the amplitude, we generated data with a signal-to-noise ratio (SNR) defined as $\mathrm{SNR}\equiv \mu /\sigma =100$. For the phase, we added noise with a standard deviation equal to 0.1 deg. These amplitude and phase noise levels were chosen based on previously published in vivo data in adults.³⁰ We assumed that noise was independent of wavelength and source–detector separation. This roughly resembles the case in practice wherein the detected intensities at short source–detector separations are attenuated to approximately the same scale as the intensities at longer source–detector separations (i.e., to reduce the dynamic range of detection across separations). However, it will be interesting to consider more complex noise models in future work.

The correlation diffusion Green’s function solution depends on the same parameters as the frequency-domain solution and additionally depends on the blood flow indices ($F_1$ and $F_2$), light wavelength ($\lambda$), and delay-time ($\tau$). For $\lambda =785\mathrm{nm}$, we evaluated the solution at two ${\rho }_i$ (0.8 and 2.5 cm) and 100 different $\tau$ (spanning $0.6\mu \mathrm{s}$ to 3.7 ms in a multitau scheme³¹). Specifically, for each FD-DOS optical properties combination, the correlation diffusion solution was evaluated for a randomly selected $F_1$ and $F_2$ combination. To mimic adult humans, $F_1$ was selected between ${10}^{-9}$ and $2\times {10}^{-8}{\mathrm{cm}}^2/\mathrm{s}$, and $F_2$ was selected between ${10}^{-9}$ and ${10}^{-7}{\mathrm{cm}}^2/\mathrm{s}$. The normalized intensity autocorrelation function was then obtained via the Siegert relation,^1,32 i.e., $g_2^{(\mathrm{theo},2L)}({\rho }_i,\tau )=1+\beta {|G_1^{(\mathrm{theo},2L)}({\rho }_i,\tau )/G_1^{(\mathrm{theo},2L)}({\rho }_i,0)|}^2$, where $\beta =0.5$ was assumed.

Intensity autocorrelation noise, derived for three different integration times (i.e., $T=0.1$, 1, and 10 s), was independently added to each $g_2^{(\mathrm{theo},2L)}({\rho }_i,\tau )$ evaluation to obtain synthetic DCS data as a function of $T$, i.e., $g_2^{(\mathrm{meas},T)}({\rho }_i,\tau )$. The autocorrelation noise was derived using a correlation noise model³³ evaluated with DCS photon count rates of 200 and 40 kHz for the short and long source–detector separations, respectively. Note that 40 kHz is on the high end for the 2.5-cm source–detector separation, but it is still within the range observed in previously published in vivo measurements on adults.³⁴ We added random Gaussian noise (with zero mean and a standard deviation based on the correlation noise model described above) independently for each delay-time and source–detector separation and independently for each combination of optical properties and flow indices.

2.3.

NIRFASTer Simulations

We used the open-source finite-element software package NIRFASTer^35,36 to generate additional synthetic datasets for the same set of eight FD-DOS (${\rho }_i=0.8$, 1.2, 1.6, 2.0, 2.8, 3.2, 3.6, and 4.0 cm) and two DCS source–detector separations (${\rho }_i=0.8$ and 2.5 cm) placed in two different geometries. The first geometry was a $10\times 10\times 10{\mathrm{cm}}^3$ 2L cube [Fig. 1(b)] with a node size of 0.07 cm, which provided a final mesh containing 482,460 nodes. The extracerebral layer thickness and absorption coefficient were set to $\ell =1.2\mathrm{cm}$ and ${\mu }_{a,1}=0.1{\mathrm{cm}}^{-1}$, respectively. The reduced scattering coefficients of both layers in the cube were set to the same value, i.e., ${\mu }_s^{\prime }={\mu }_{s,1}^{\prime }={\mu }_{s,2}^{\prime }=10{\mathrm{cm}}^{-1}$, and held constant. NIRFASTer was then used to simulate $\mathrm{AC}({\rho }_i)$ and $\theta ({\rho }_i)$ via a finite-element method for 11 evenly spaced cerebral layer absorption coefficients (${\mu }_{a,2}$) between 0.08 and $0.18{\mathrm{cm}}^{-1}$. Similar to the forward-model simulations, Gaussian noise was added to $\mathrm{AC}({\rho }_i)$ and $\theta ({\rho }_i)$ to obtain 20 different pairs of ${\mathrm{AC}}_{\mathrm{meas}}({\rho }_i)$ and ${\theta }_{\mathrm{meas}}({\rho }_i)$ synthetic data for each value of ${\mu }_{a,2}$ (amplitude $\mathrm{SNR}=100$; phase $\sigma =0.1\mathrm{deg}$).

For each combination of optical properties, NIRFASTer was also used to generate $G_1^{(\mathrm{theor},2\mathrm{L})}({\rho }_i,\tau )$ via a finite-element method for 16 different CBF indices between ${10}^{-9}$ and ${10}^{-7}{\mathrm{cm}}^2/\mathrm{s}$ (the extracerebral flow index was held constant at $F_1={10}^{-8}{\mathrm{cm}}^2/\mathrm{s}$). Then as described in Sec. 2.2, correlation noise was added to $G_1^{(\mathrm{theor},2\mathrm{L})}({\rho }_i,\tau )$ for a given DCS integration time $T$ to obtain a synthetic DCS measurement (i.e., $g_2^{(\mathrm{meas},T)}({\rho }_i,\tau )$) independently for each source–detector separation. For each integration time ($T=0.1$, 1, and 10 s) and each combination of optical properties, flow indices, and noise additions from FD-DOS, 15 synthetic DCS measurements were generated (in total, we generated 300 autocorrelation curves for each source–detector separation at each combination of optical property, flow, and integration time).

The second geometry was a realistic adult head mesh created using an open-source library (brain2mesh, with a Delaunay sphere radius of 0.11 cm, radius-to-edge ratio of 1.24, and maximum element volume of $4{\mathrm{mm}}^3$).³⁷ The head was segmented into the scalp, skull, cerebral spinal fluid (CSF), white matter, and gray matter, containing $\sim 1.4$ million nodes. We removed the nodes further than 10 cm from the simulated source, reducing the final mesh to 663,470 nodes. For our simulations, the scalp and skull were merged to form one homogeneous tissue type (i.e., the extracerebral layer), whereas the CSF, gray matter, and white matter were merged to form a second homogeneous tissue type (i.e., the cerebral layer). The synthetic data for this geometry were generated with NIRFASTer in the same manner as the cube simulations (including the same combinations of extracerebral and cerebral tissue properties and noise additions). Note that, although the scalp and skull blood flows are quite different under normal conditions, the concatenation of the scalp and skull into one layer is closer to reality for applications wherein a transient high probe pressure can be applied against the scalp to reduce the scalp flow closer to levels in the skull.¹² Our results are thus most relevant for these conditions.

We did, however, conduct a pilot test of the algorithm under conditions of the scalp flow being higher than the skull flow. This test used simulated data for the same realistic head geometry, except that the scalp and skull tissues were assigned distinct optical properties and flow indices (i.e., a three-layer realistic head geometry). Specifically, we simulated data in which the absorption coefficients for the scalp and skull were ${\mu }_{\text{scalp}}=0.1{\mathrm{cm}}^{-1}$ and ${\mu }_{\text{skull}}=0.15{\mathrm{cm}}^{-1}$,²⁸ respectively, and the blood flow indices were equal to $F_{\text{scalp}}={10}^{-8}{\mathrm{cm}}^2/\mathrm{s}$ and $F_{\mathrm{skull}}={10}^{-9}{\mathrm{cm}}^2/\mathrm{s}$, respectively. Here we varied the true cerebral absorption coefficient (${\mu }_{a,2,\mathrm{act}}$) between 0.08 and $0.16{\mathrm{cm}}^{-1}$ (in steps of $0.02{\mathrm{cm}}^{-1}$). For each change in cerebral absorption, we also varied the cerebral flow ($F_{2,\mathrm{act}}$) between 4 and $10\times {10}^{-8}{\mathrm{cm}}^2/\mathrm{s}$ (in steps of ${10}^{-8}{\mathrm{cm}}^2/\mathrm{s}$). As with the other simulations, homogeneous scattering was assumed (i.e., ${\mu }_{s,\text{scalp}}^{\prime }={\mu }_{s,\text{skull}}^{\prime }={\mu }_{s,2}^{\prime }=10{\mathrm{cm}}^{-1}$), and we added random Gaussian noise to generate multiple datasets from each simulation.

2.4.

Liquid Phantom Experiments

Finally, we performed in vitro testing of the algorithm in a 2L liquid phantom. To this end, we developed a custom black acrylic aquarium to mimic a 2L cube [Fig. 2(a)]. We used a removable thin plastic film attached to a black frame to separate the two layers. To attach the optical fibers, we drilled holes in one side of the phantom, such that fibers were in direct contact with the liquid. We glued neutral-density filters in some fiber positions to avoid saturation of the shorter source–detector separations. Finally, to allow flow changes in the second layer, we also attached a peristaltic pump to the aquarium’s lateral sides [Fig. 2(a)].

Fig. 2

2L phantom setup. (a) We used a hybrid diffuse optical system composed of an FD-DOS and a DCS module. For the phantom experiment, we developed a black acrylic aquarium capable of mimicking a 2L geometry. To separate each layer, we inserted a black frame with a thin plastic film. The second layer of the phantom was connected to a peristaltic pump. (b) The experiment was performed with an optical sensor that allowed for measurements with eight different SDS ($\rho =0.7$, 1.2, 1.5, 1.7, 2.0, 2.2, 2.5, and 3.0 cm) for FD-DOS and four SDS ($\rho =0.8$, 1.5, 2.0, and 2.5 cm) for DCS. (c) We varied the second layer’s absorption a total of eight times by adding an ink solution to the second layer. For the 2L phantom experiment, we also varied the flow in the second layer by varying the flow output of the peristaltic pump after each ink addition.

For our phantom experiments, we used a hybrid diffuse optical system developed at the University of Campinas, which is described elsewhere.^38,39 Briefly, the system combines a commercial FD-DOS system (Imagent, ISS) and a homemade DCS system. The FD-DOS system contains 4 detectors and 32 lasers equally split among 4 wavelengths (685, 705, 750, and 830 nm). The DCS system contains a single-laser emitting at 785 nm (CrystaLaser) and 16 single-photon counters (SPCM-AQ4C, Excelitas). We used two detectors and 16 sources from the FD-DOS module for the experimental phantom measurements and one source and twelve detectors from the DCS module (two detectors at 0.7 cm, four at 1.5 cm, three at 2 cm, and four at 2.5 cm). Our phantom [Fig. 3(b)] allowed for FD-DOS measurements with eight different source–detector separations ($\rho =0.7$, 1.2, 1.5, 1.7, 2.0, 2.2, 2.5, and 3.0 cm), and DCS measurements with four source–detector separations ($\rho =0.7$, 1.5, 2.0, and 2.5 cm). Prior to each experiment, the FD-DOS measurements were calibrated using three different solid phantoms, as described elsewhere.^39,40

Fig. 3

2L fitting scheme. The scheme assumes a priori knowledge of the extracerebral layer thickness ($\ell$). First, we fit the multidistance FD-DOS amplitude (${\mathrm{AC}}_{\mathrm{meas}}({\rho }_i)$) and phase (${\theta }_{\mathrm{meas}}({\rho }_i)$) data to the two-layer (2L) cylindrical frequency-domain Green’s function solution to recover the extracerebral and cerebral layer absorption coefficients (${\mu }_{a,1}$, ${\mu }_{a,2}$), and the reduced scattering coefficient (${\mu }_s^{\prime }$; the scheme assumes ${\mu }_{s,1}^{\prime }={\mu }_{s,2}^{\prime }={\mu }_s^{\prime }$). Next, using the recovered ${\mu }_{a,1}$ and ${\mu }_s^{\prime }$ as inputs, we fit the DCS measurement at the short source–detector separation ($g_2^{(\mathrm{meas},T)}({\rho }_S,\tau )$) to the SI correlation diffusion Green’s function solution to recover the extracerebral flow index ($F_1$), assuming knowledge of the $\beta$ factor from Siegert’s relation (${\beta }_s=0.5$). Finally, using ${\mu }_{a,1}$, ${\mu }_{a,2}$, ${\mu }_s^{\prime }$, and $F_1$ as inputs, we fit the DCS measurements at a long source–detector separation ($g_2^{(\mathrm{meas},T)}({\rho }_L,\tau )$) to the 2L cylindrical Green’s function solution to recover the cerebral flow index ($F_2$), assuming ${\beta }_L=0.5$.

Using the setup described above, we performed two separate phantom experiments. In the first experiment, we measured the phantom in a homogeneous geometry (i.e., without the thin plastic film separating the layers). We started with a solution consisting of 4.5 L of distilled water, 200 mL of Liponfundin 20%, and 2 mL of an ink solution. The ink solution was made with 0.5 mL of ink (Nankin, Acrilex) and 25 mL of distilled water. After measuring the initial solution, we increased the homogeneous medium’s absorption coefficient by repeatedly adding 0.5 mL of the ink solution to the phantom; in total, we simulated eight different absorption coefficients. After each ink addition, we waited at least 6 min for the solution to reach equilibrium before moving to the next step. The data from this homogeneous experiment were used to estimate the true optical properties for different volume ratios of the liquid phantom.

In the second experiment, we started with a solution with the same ink and intralipid concentrations used in the start of the homogenous phantom experiment described above. However, before any ink additions, we placed a thin plastic film at 1.2 cm from the phantom wall to simulate a 2L geometry; data collection started after placing the thin film. Here we varied the absorption coefficient of the second layer via multiple additions of 0.5 mL of the ink solution. After each ink addition, we also changed the second layer’s flow by changing the flow output of a peristaltic pump between 1.5 and $3\mathrm{L}/\min$ [Fig. 3(c)]. Similar to the homogeneous phantom experiment, we waited at least 6 min after each absorption and flow change step, and we repeated this process to simulate eight distinct second-layer absorption coefficients.

2.5.

Two-Layer Fitting Algorithm

The 2L fitting scheme is depicted in Fig. 3. The scheme assumes that homogeneous tissue reduced scattering (i.e., ${\mu }_{s,1}^{\prime }={\mu }_{s,2}^{\prime }={\mu }_s^{\prime }$; we performed an analysis to justify this assumption, which we discuss below) and that the extracerebral layer thickness is known a priori. We first used a nonlinear constrained global optimizer implemented in MATLAB R2020a (fmincon, Mathworks, Natick, Massachusetts, United States) to obtain estimates of ${\mu }_{a,1}$, ${\mu }_{a,2}$, and ${\mu }_s^{\prime }$ by fitting multidistance FD-DOS data (i.e., ${\mathrm{AC}}_{\mathrm{meas}}({\rho }_i)$, ${\theta }_{\mathrm{meas}}({\rho }_i)$) to the 2L analytical frequency-domain Green’s function solution (i.e., ${\mathrm{AC}}_{\mathrm{theo},2\mathrm{L}}({\rho }_i)$ and ${\theta }_{\mathrm{theo},2L}({\rho }_i)$, see Sec. 2.2). Specifically, we used fmincon to find the set of ${\mu }_{a,1}$, ${\mu }_{a,2}$, and ${\mu }_s^{\prime }$ parameters that minimize the cost function ${\chi }_{\mathrm{FD}}=\sqrt{{\chi }_{\mathrm{AC}}^2}+\sqrt{{\chi }_{\mathrm{F}\mathrm{D}}^2}$, where ${\chi }_{\mathrm{AC}}$ and ${\chi }_{\theta }$ are defined as

Here $N=8$ is the total number of source–detector distances. The minimization was also subject to the following constraints: $0.005\le {\mu }_{a,1}({\mathrm{cm}}^{-1})<0.6$, $0.005\le {\mu }_{a,2}({\mathrm{cm}}^{-1})<0.6$, and $4\le {\mu }_{s,i}^{\prime }({\mathrm{cm}}^{-1})<20$. These constraints were based on an adult head’s expected ranges of optical properties.^27,28 The known extracerebral layer thickness, an index of refraction ($n=1.4$), and a cylindrical radius ($a=30\mathrm{cm}$) were used as inputs in the minimization, and the initial guesses used for ${\mu }_{a,1}$, ${\mu }_{a,2}$, and ${\mu }_s^{\prime }$ in the minimization were 0.1, 0.1, and $10{\mathrm{cm}}^{-1}$, respectively. The normalization of the amplitude and phase data by their values at the shortest source–detector distance in the cost functions removes the need to fit for additional amplitude and phase scaling factors. Additionally, because the amplitude decreases exponentially with increasing the source–detector distance, we used the logarithm of the amplitude to minimize bias to the shorter distances (i.e., such that fitting errors at each distance are weighted approximately evenly in the cost function).

In the next step, we fit the short-separation DCS data (i.e., $g_2^{(\mathrm{meas},T)}({\rho }_s,\tau )$, see Secs. 2.2 and 2.3) to the SI correlation diffusion solution to obtain the extracerebral flow index $F_1$ given the ${\mu }_{a,1}$ and ${\mu }_s^{\prime }$ inputs determined from FD-DOS. The short separation (i.e., ${\rho }_s=0.8\mathrm{cm}$) was chosen such that the detected light is predominantly confined to the extracerebral layer for the adult head geometry.^10,41 The use of a homogeneous SI tissue model for the short-separation data is thus reasonable. We employed the same nonlinear optimizer (fmincon) to find an $F_1$ value that minimizes the cost function ${\chi }_{\mathrm{DCS},{\rho }_S}$, which is defined as

In the third and final step, we fit the long-separation DCS data [i.e., $g_2^{(\mathrm{meas},T)}({\rho }_L,\tau )$, see Secs. 2.2 and 2.3] to the 2L correlation diffusion solution to obtain the cerebral flow index $F_2$ given the inputs of ${\mu }_{a,1}$, ${\mu }_{a,2}$, ${\mu }_s^{\prime }$, and $F_1$ from the first two steps. Additional inputs in the fit were the extracerebral layer thickness, index of refraction, DCS wavelength ($\lambda =785\mathrm{nm}$), the cylindrical radius of $a=30\mathrm{cm}$, and ${\beta }_L=0.5$. We used fmincon to find the $F_2$ value that minimizes the cost function ${\chi }_{\mathrm{DCS},{\rho }_L}$, which is defined as

2.6.

Data Analysis

3.1.

Accuracy of the Two-Layer and Homogenous Approaches

The first step of our 2L fitting algorithm was the recovery of the optical properties of each layer from the FD-DOS measures of ${\mathrm{AC}}_{\mathrm{meas}}$ and ${\theta }_{\mathrm{meas}}$. With our algorithm, we recovered the tissue absorption and reduced scattering coefficients with excellent agreement between the recovered and actual values for the forward-model, 2L cube, and 2L realistic head simulations (see Fig. 4 and Table 1; exemplar ${\mathrm{AC}}_{\mathrm{meas}}$ and ${\theta }_{\mathrm{meas}}$ fits shown in Appendix B). In these geometries, median errors were $<8\%$. The best-fit linear regression lines for the comparison of ${\mu }_{a,2}$ and ${\mu }_{a,2,\mathrm{act}}$ approached the unity line. However, the agreement for the SI analysis was not as good; the slope of the linear best-fit line between ${\mu }_{a,\mathrm{SI}}$ and ${\mu }_{a,2,\mathrm{act}}$ (i.e., 0.5) deviated from the unity line.

Fig. 4

Recovered versus actual tissue optical properties. (a)–(c) The recovered tissue absorption (${\mu }_{a,i}$) and tissue reduced scattering coefficient (${\mu }_{s,i}^{\prime }$) for the extracerebral (green lines) and cerebral (red lines) layers are plotted against the actual values of the second-layer absorption (${\mu }_{a,2,\mathrm{act}}$) for the (a), (d) forward-model, (b), (e) cube, and (c), (f) realistic head simulations (circles denote the medians of the recovered values across all simulations run for each actual value; shaded areas represent the IQR). The corresponding recovered tissue absorption (${\mu }_{a,\mathrm{SI}}$) and reduced scattering coefficients from the SI model are also plotted against the actual cerebral absorption values (blue diamonds). Dashed lines represent the actual relationships between each parameter and the cerebral absorption coefficient.

Table 1

MAPE of the optical properties (μa,1,μa,2, and μs′) and flow indices (F1 and F2) recovered with the 2L approach and with the SI approach (μa,SI, μs,SI′, and FSI) for all datasets generated. The linear best-fit relations between the recovered and actual values for the second layer are also reported.

For assessing the reconstructed $F_1$ and $F_2$ accuracy, we used the simulated datasets with a DCS integration time of $T=10\mathrm{s}$. For the forward-model and 2L cube simulations, we were able to accurately recover $F_1$ with median errors below 3% [Table 1, Figs. 5(a) and 5(b); the exemplar intensity autocorrelation function fits are shown in Appendix B]. For the 2L realistic head simulations, our method of using an SI model to recover $F_1$ from a short DCS source–detector separation was modestly less accurate, with errors around 13% [Table 1 and Fig. 5(c)].

Fig. 5

Recovered versus actual flow indices. The (a)–(c) recovered extracerebral flow ($F_1$, green lines) and (d)–(f) cerebral flow ($F_2$, red lines) indices are plotted against the actual values of the second-layer flow ($F_{2,\mathrm{act}}$) for the [(b), (e)] cube and [(c), (f)] realistic head simulations. For the forward-model simulations [(a), (d)], we plot the extracerebral and CBF values against their actual values. The corresponding recovered flow index ($F_{\mathrm{SI}}$) from the SI model are also plotted against the actual values (blue diamonds). Dashed lines represent the actual relationships between each parameter and the $x$ axis. In all cases, circles denote the medians of the recovered values across all simulations (with varying noise, flow indices, and varying absorption), and shaded areas represent the IQR.

We observed excellent agreement between the recovered and actual $F_2$ for the forward-model and 2L cube simulations [Table 1 and Figs. 5(d) and 5I]. For both datasets, the errors were $<10\%$ on average, and the best-fit linear regression lines approached the unity line (Table 1). In the 2L realistic head simulations, however, the recovered $F_2$ systematically underestimated the true value by a median error of 34% [Table 1 and Fig. 5(f)]. The small IQRs of the recovered flow values demonstrate robustness against noise and optical absorption changes.

When neglecting the extracerebral layer using an SI model to estimate $F_{2,\mathrm{act}}$, the systematic errors (i.e., $\mathrm{MAPE}>69\%$) were larger than the errors recovered with our 2L approach in all simulated datasets ($p<0.001$). The SI homogeneous model recovered the correct directional trends for the 2L cube and realistic head simulations, where the first-layer flow was held constant. However, for the forward-model simulations, the recovered $F_{\mathrm{SI}}$ values were highly sensitive to variations in first layer flow [$F_{1,\mathrm{act}}$, Fig. 5(a)]. Note that the small IQRs for $F_2$ across variations in extracerebral blood flow indicate minimal cross talk between extracerebral and CBF [Fig. 5(d)].

Finally, for conditions wherein scalp and skull blood flow are substantially different (i.e., the three-layer realistic head geometry), the recovered errors with the constrained 2L algorithm were substantially larger than those for the 2L head simulations, but they were still smaller than the errors in the SI estimates (see Appendix C).

3.2.

Sensitivity of the FD-DOS Two-Layer Green’s Function to Tissue Optical Property Changes

We found that the 2L cylindrical FD-DOS Green’s function solution is minimally sensitive to changes in ${\mu }_{s,2}^{\prime }$ for source–detector distances ($\rho$) up to 5 cm (Fig. 6). Variations in the solutions for FD-DOS amplitude (${\mathrm{AC}}_{\mathrm{theo},2\mathrm{L}}(\rho )$) and phase (${\theta }_{\mathrm{theo},2\mathrm{L}}(\rho )$) by variations in ${\mu }_{s,2}^{\prime }$ from 5 to $15{\mathrm{cm}}^{-1}$ are smaller or on the same order of the expected noise [Figs. 4(a) and 4(d)]. Here we fixed the other tissue parameters at $\ell =1.2\mathrm{cm}$, ${\mu }_{a,1}=0.13{\mathrm{cm}}^{-1}$, ${\mu }_{a,2}=0.13{\mathrm{cm}}^{-1}$, and ${\mu }_{s,1}^{\prime }=10.5{\mathrm{cm}}^{-1}$. The sensitivities of the FD-DOS amplitude and phase to extracerebral and cerebral layer optical properties are also plotted versus source–detector distance in Fig. 6. The sensitivities are defined by the evaluation of the derivatives $\partial \log {\mathrm{AC}}_{\mathrm{theo},2\mathrm{L}}({\rho }_i)/\partial x_i$ and $\partial {\theta }_{\mathrm{theo},2\mathrm{L}}({\rho }_i)/\partial x_i$ at ${\mu }_{a,1}={\mu }_{a,2}=0.13{\mathrm{cm}}^{-1}$, ${\mu }_{s,1}^{\prime }={\mu }_{s,2}^{\prime }=10.5{\mathrm{cm}}^{-1}$, and $\ell =1.2\mathrm{cm}$ ($x_i$ denotes ${\mu }_{a,1},{\mu }_{a,2},{\mu }_{s,1}^{\prime }$, or ${\mu }_{s,2}^{\prime }$). Note that the sensitivities to ${\mu }_{s,2}^{\prime }$ are much lower than those for the other optical properties. For example, at $\rho =4\mathrm{cm}$, the sensitivities of the FD-DOS amplitude and phase to ${\mu }_{s,2}^{\prime }$ are 5% and $-3\%$ of the corresponding sensitivities to ${\mu }_{s,1}^{\prime }$ and $<0.5\%$ of the sensitivities to ${\mu }_{a,1}$ and ${\mu }_{a,2}$. Given its minimal sensitivity to the FD-DOS measurements, ${\mu }_{s,2}^{\prime }$ is not a good fitting parameter. These results justify the need to assume homogeneous tissue reduced scattering.

Fig. 6

Sensitivity of the 2L FD-DOS Green’s function solution to changes in tissue optical properties. The (a) amplitude logarithm ($\log {\mathrm{AC}}_{\mathrm{theo},2\mathrm{L}}$) and (d) phase $({\theta }_{\mathrm{theo},2\mathrm{L}})$ of the cylindrical 2L FD-DOS Green’s function solution are plotted against source–detector distance ($\rho$) for a wide range of cerebral tissue reduced scattering coefficients (${\mu }_{s,2}^{\prime }$) between 5 (blue) and $15{\mathrm{cm}}^{-1}$ (red). For each ${\mu }_{s,2}^{\prime }$ evaluation, the extracerebral and cerebral tissue absorption coefficients (${\mu }_{a,1}$ and ${\mu }_{a,2}$) were both fixed at $0.13{\mathrm{cm}}^{-1}$, the extracerebral tissue reduced scattering coefficient (${\mu }_{s,1}^{\prime }$) was fixed at $10.5{\mathrm{cm}}^{-1}$, and the extracerebral layer thickness was fixed at 1.2 cm. The sensitivities of the [(b), (c)] amplitude logarithm [$\partial \log {\mathrm{AC}}_{\mathrm{theo},2\mathrm{L}}/\partial x_i$] and [(e), (f)] phase [$\partial {\theta }_{\mathrm{theo},2\mathrm{L}}/\partial x_i$] to tissue optical property changes are also plotted against $\rho$; $x_i$ refers to ${\mu }_{a,1}$, ${\mu }_{s,1}^{\prime }$, ${\mu }_{a,2}$, and ${\mu }_{s,2}^{\prime }$. All derivatives were evaluated at the same optical properties used for (a) and (d) with ${\mu }_{s,2}^{\prime }$ fixed at $10.5{\mathrm{cm}}^{-1}$.

3.3.

Error Arising from the Assumption of Homogeneous Tissue Reduced Scattering

Although the assumption of homogeneous tissue reduced scattering is justified for 2L fitting, inhomogeneous tissue reduced scattering remains a source of error in the recovery of $F_2$ and ${\mu }_{a,2}$. As a pilot characterization of this issue, we visualized the median percent error of the recovered $F_2$ and ${\mu }_{a,2}$ values as a function of the ${\mu }_{s,2,\mathrm{act}}^{\prime }/{\mu }_{s,1,\mathrm{act}}^{\prime }$ ratio in the forward-model simulation dataset. Specifically, the median and IQR of the percent errors for the simulations run at each ${\mu }_{s,2,\mathrm{act}}^{\prime }/{\mu }_{s,1,\mathrm{act}}^{\prime }$ ratio are plotted against ${\mu }_{s,2,\mathrm{act}}^{\prime }/{\mu }_{s,1,\mathrm{act}}^{\prime }$ for different extracerebral layer thicknesses (Fig. 7). The magnitudes of the errors from inhomogeneous scattering in the recovered $F_2$ are smaller than those in the recovered ${\mu }_{a,2}$. For example, at ${\mu }_{s,2,\mathrm{act}}^{\prime }/{\mu }_{s,1,\mathrm{act}}^{\prime }=0.8$ and $\ell =1.2\mathrm{cm}$, the median percent errors in the recovered ${\mu }_{a,2}$ and $F_2$ are 15.8% and $-8.8\%$, respectively. Note also that, as $\ell$ increases, the percent errors are more variable across all ${\mu }_{s,2,\mathrm{act}}^{\prime }/{\mu }_{s,1,\mathrm{act}}^{\prime }$ ratios (i.e., the IQR of the errors is larger, see Fig. 7). This larger variability arises because of lower brain sensitivities (fitting of parameters with lower brain sensitivities are more prone to crosstalk caused by measurement noise). Additionally, as the extracerebral thickness increases, the sensitivity of the recovered ${\mu }_{a,2}$ error to ${\mu }_{s,2,\mathrm{act}}^{\prime }/{\mu }_{s,1,\mathrm{act}}^{\prime }$ decreases.

Fig. 7

Errors arising from inhomogeneous tissue reduced scattering. Median (circles) and IQRs (shaded area) of the percent errors in the 2L recovered (a)–(c) CBF and (d)–(f) cerebral absorption coefficient plotted against the cerebral to extracerebral ratio of the actual tissue reduced scattering coefficients (the median and IQR are across all simulations ran for each ratio). These plots were generated from the forward-model simulations of three extracerebral thicknesses [(a), (d) $\ell =1.0\mathrm{cm}$; (b), (e) $\ell =1.2\mathrm{cm}$; and (c), (f) $\ell =1.6\mathrm{cm}$].

3.4.

Error Arising from the Assumption of Homogeneous Tissue Absorption

To test the effects of assuming homogeneous tissue absorption as an additional constraint in the recovery of flow indices, we reanalyzed our NIRFASTer simulations using homogeneous optical properties (i.e., using an SI model for FD-DOS, in which we assume ${\mu }_{a,1}={\mu }_{a,2}={\mu }_{a,\mathrm{SI}}$ and ${\mu }_{s,1}^{\prime }={\mu }_{s,2}^{\prime }={\mu }_{s,\mathrm{SI},}^{\prime }$), but a 2L model for DCS (to separately recover $F_1$ and $F_2$). Here we focused on cases in which $F_{1,\mathrm{act}}$, $F_{2,\mathrm{act}}$, ${\mu }_{a,1,\mathrm{act}}$, ${\mu }_{s,1,\mathrm{act}}^{\prime }$, and ${\mu }_{s,2,\mathrm{act}}^{\prime }$ were fixed at ${10}^{-8}{\mathrm{cm}}^2/\mathrm{s}$, $6\times {10}^{-8}{\mathrm{cm}}^2/\mathrm{s}$, $0.1{\mathrm{cm}}^{-1}$, $10{\mathrm{cm}}^{-1}$, and $10{\mathrm{cm}}^{-1}$, respectively, whereas ${\mu }_{a,2,\mathrm{act}}$ varied between 0.08 and $0.18{\mathrm{cm}}^{-1}$. When homogeneous tissue absorption was assumed, the actual changes in ${\mu }_{a,2}$ translated to erroneous changes in $F_2$ in both geometries (Fig. 8).

Fig. 8

Errors in the recovery of CBF if homogeneous tissue absorption is assumed. The 2L recovered CBF ($F_2$) for the (a) cube and (b) 2L realistic head simulations plotted against the actual cerebral tissue absorption coefficient (${\mu }_{a,2,\mathrm{act}}$; all other actual tissue flow and optical property parameters were held fixed [see Sec. 2.5)]. $F_2$ was recovered using: (1) homogeneous absorption and scattering derived from the SI model (purple squares) and (2) the primary 2L fitting scheme depicted in Fig. 2 (red circles). The purple square and red circles denote the medians across all simulations ran for each ${\mu }_{a,2,\mathrm{act}}$, and the shaded areas denote the IQRs.

3.5.

Error Arising from DCS Correlation Noise

For our fourth secondary analysis, we examined the influence of DCS correlation noise on the accuracy of CBF measurements recovered with the 2L scheme ($F_2$) and with the SI scheme ($F_{\mathrm{SI}}$). DCS correlation noise increases with decreasing DCS integration time ($T$). Correlation noise did not substantially influence the accuracy of the recovered $F_{\mathrm{SI}}$ values in any of the 2L geometries (Fig. 9). The MAPE and IQR of the absolute percent errors were comparable for all three $T$. For the 2L scheme, however, there was a strong effect (Fig. 9). The MAPE and IQR of the absolute percent errors in $F_2$ are noticeably larger at $T=0.1\mathrm{s}$ than at $T=10\mathrm{s}$ for all three geometries. As expected, the 2L model reconstruction is thus more susceptible to correlation noise than the SI model. Note, however, that the reconstruction errors at $T=1\mathrm{s}$ are only minimally to modestly larger than those at $T=10\mathrm{s}$, depending on the geometry.

Fig. 9

Correlation noise effects on the recovered CBF accuracy. The MAPE (circles) and the IQR of the absolute percent errors (dashed lines) in the SI ($F_{\mathrm{SI}}$) and 2L ($F_2$) CBF reconstructions are plotted against DCS integration time ($T$). The median and IQR were computed across all (a) forward-model simulations, (b) 2L cube simulations, and (c) 2L realistic head simulations ran for each $T$.

Finally, the CNR estimated from the realistic head simulations (see Sec. 2.6.5) for $F_2$ at $T=0.1$, 1 and 10 s were 0.7, 2.5, and 6.0, respectively, and the corresponding CNR for $F_{\mathrm{SI}}$ were 0.7, 1.3, and 1.4, respectively. Note that, although these SI CNR levels were obtained from fitting the upper half of the autocorrelation curve ($g_2\ge 1.25$), the SI CNR levels were the same if the entire $g_2$ curve was fit for instead ($g_2>1$).

3.6.

Error Arising from Inaccurate Extracerebral Thickness

We also used the realistic head simulations to investigate the influence of errors in the extracerebral layer thickness on the recovery of $F_2$ and ${\mu }_{a,2}$. The influence was considerable for the $F_2$ recovery but more modest for the ${\mu }_{a,2}$ recovery (Fig. 10). For example, $\pm 0.2$ cm errors in $\ell$ resulted in median errors of up to 15% in ${\mu }_{a,2}$ and up to 60% in $F_2$. Surprisingly, the minimum error for the recovery of $F_2$ occurred when the extracerebral layer thickness was overestimated by $\mathrm{\Delta }\ell \approx 0.1\mathrm{cm}$.

Fig. 10

Reconstruction errors in CBF and tissue absorption arising from inaccurate extracerebral layer thickness ($\ell$). Errors in the recovered (a) CBF ($F_2$) and (b) cerebral tissue absorption coefficient (${\mu }_{a,2}$) plotted against errors in the extracerebral layer thickness used for the 2L fits of the realistic head synthetic data. $\mathrm{\Delta }\ell$ is the difference between the extracerebral layer thickness used in the fits and the actual extracerebral layer thickness (1.22 cm). The circles and dashed lines denote the median and IQR of the absolute percent errors across all simulations with actual CBF larger than actual extracerebral blood flow.

3.7.

Accuracy of This Methodology in a Two-Layer Liquid Phatom

Figure 11 and Table 1 show the comparison of the absolute values recovered with the proposed 2L fitting algorithm and the expected values from the phantom experiment at the 705-nm wavelength. Using our 2L fitting algorithm, we obtained MAPE [IQR] of 10 [9, 11]% for ${\mu }_{a,1}$, 8 [5, 12]% for ${\mu }_{a,2}$, and 10.2 [10.1, 10.4]% for ${\mu }_s^{\prime }$. Similar to our simulations, the standard SI homogeneous approximation recovered the optical properties with larger errors (MAPE [IQR] equal to 19 [9, 27]% for ${\mu }_{a,\mathrm{SI}}$ and 14 [11, 16]% for ${\mu }_{s,\mathrm{SI}}^{\prime }$). Note that the accuracies for the 685-nm wavelength were very similar to these results (data not shown); the data at the 750 and 830 nm wavelengths were unusable due to calibration and SNR issues.

Fig. 11

Recovered optical properties and flow indices for the 2L phantom experiment. The 2L and SI (a) recovered absorption coefficients and (b) reduced scattering coefficients are plotted against the actual second-layer values (the actual second-layer values were estimated based on a homogeneous phantom measurement that contained the same volume ratio of ink and intralipid). (c) The 2L and SI recovered flow indices are plotted against the second-layer pump flow used in the experiment. The green lines denote the recovered absorption coefficient (${\mu }_{a,1}$) and flow index ($F_1$) for the first layer from the 2L fitting algorithm. The red lines denote the recovered absorption coefficient (${\mu }_{a,2}$) and flow index ($F_2$) for the second layer, as well as the recovered homogenous reduced scattering coefficient (${\mu }_s^{\prime }$), from the 2L fitting algorithm. Blue lines denote the values recovered using the SI model (i.e., ${\mu }_{a,\mathrm{SI}}$, ${\mu }_{\mathrm{s},\mathrm{SI}}^{\prime }$, and $F_{\mathrm{SI}}$). Each point denotes the median of the recovered values across all measurements [i.e., across either varying pump flows for (a) and (b) or varying second-layer absorption for (c); see text]; shaded areas represent the IQR. The IQRs for (a) and (b) were smaller than 0.002 and $0.03{\mathrm{cm}}^{-1}$, respectively, and thus are not apparent.

As expected, the flow indices recovered for the first layer were relatively independent of the second layer’s changes in absorption and flow [Fig. 11(c)]. Additionally, flow indices for the second layer were not entangled with the second layer’s absorption changes [as seen by the small IQR in Fig. 11(c)]. When using an SI homogeneous model to estimate changes in the second layer, we were able to recover the correct trend in the flow changes. However, the absolute values recovered using the SI approximation were significantly lower than those retrieved with our 2L approach. Of note, for the DCS measurements, we recovered $\beta =0.45$ for the short SDS (0.7 cm) and 0.49 for the long SDS (2.5 cm).

5.1.

Semi-Infinite Media

To obtain the solution for a homogenous SI media, we consider the source, $S(\overrightarrow{r},\omega )$ to be a punctual source located at $z_0={({\mu }_a+{\mu }_s^{\prime })}^{-1}$.^1,3 Using the extrapolated zero boundary condition, we arrive at the SI solution as:

Here $r_1=\sqrt{(z-z_0{)}^2+{\rho }^2}$ and $r_2=\sqrt{(z+2z_b+z_0{)}^2+{\rho }^2}$, and $z_b$ is defined as

5.2.

Two-Layer Media

In this study, we opted to use the solution of the photon diffusion equation for a 2L cylinder as it is more computationally robust than the standard 2L solution for a laterally unbounded medium.^16,23 Specifically, we modeled the tissue as an infinitely thick cylinder with radius $a$, composed of two layers: the first layer, with thickness $\ell$, represents the extracerebral tissues; the second layer is infinitely thick and means the cerebral tissues [Fig. 1(a)]. Although we restricted our discussion to the FD-DOS and CW-DCS solutions, our results can be extended to time-domain measurements by applying a Fourier-transform to the solution presented below.

By solving the FD-DOS diffusion equation [Eq. (5) for a 2L cylinder, it is possible to show that the fluence in the $k$’th layer can be written as^16,23