2.1.

EAC Computation

EAC estimation relies on the known exponential decay in optical fluence rate as a function of distance from the illumination point on the tissue surface. The slope (i.e., the exponential factor linking optical fluence to source–detector distance) of this decay function is determined by the optical properties of the tissue (and, more precisely, by a combination of the absorption and reduced scattering coefficients of the tissue, EAC), as well as by its geometric properties and boundary conditions. Specifically, for a semi-infinite, homogeneous medium with zero boundary conditions, the CW reflectance $I(r)$ recorded at a distance $r$ from a source is linked to the EAC through the following equation (derived from Ref. 21):

2.2.

Light Modeling

NIR light propagation in biological tissue is well approximated by the diffusion equation.³² In order to test our assumptions and to run comparisons between simulations and real results, we used both analytical solutions³³ and FEM approaches applied to the diffusion equation. ³⁴

For FEM simulations, high-resolution meshes (maximum tetrahedral $\text{volume}=2{\mathrm{mm}}^3$) were generated using the MATLAB^® function iso2mesh³⁵ for homogeneous slab geometry, homogeneous spherical geometry, and heterogeneous head geometry. The heterogeneous head geometries were based on segmented T1w MR images. Segmentation of subjects’ heads into skull and scalp, CSF, white matter and gray matter [an example can be seen in Fig. 11(b)] was performed using statistical parametric mapping functions applied to T1w images.³⁶ Baseline optical properties (${\mu }_{\mathrm{a}}$, ${\mu }_{\mathrm{s}}^{\prime }$, and refraction index $\eta$) of the tissues at the relevant wavelength were taken from previous papers.^30,37 A combination of T1w and T2w segmentation is required to differentiate skull and scalp. Here, for practical purposes, and given the similarities of the optical properties of skull and scalp, we only used T1w images and attributed the skull’s values to both structures. Optical properties vary considerably more among skull, CSF, and brain. The optical values at the wavelength of interest (830 nm) were set to the following values: scalp and skull: ${\mu }_{\mathrm{a}}=0.014{\mathrm{mm}}^{-1}$, ${\mu }_{\mathrm{s}}^{\prime }=0.84{\mathrm{mm}}^{-1}$; CSF: ${\mu }_{\mathrm{a}}=0.004{\mathrm{mm}}^{-1}$, ${\mu }_{\mathrm{s}}^{\prime }=0.3{\mathrm{mm}}^{-1}$; gray matter: ${\mu }_{\mathrm{a}}=0.019{\mathrm{mm}}^{-1}$, ${\mu }_{\mathrm{s}}^{\prime }=0.673{\mathrm{mm}}^{-1}$; white matter: ${\mu }_{\mathrm{a}}=0.021{\mathrm{mm}}^{-1}$, ${\mu }_{\mathrm{s}}^{\prime }=1.01{\mathrm{mm}}^{-1}$. The refraction index was set at $\eta =1.4$ for all compartments, simulations, and mediums considered. In order to improve the heterogeneous model, veins were added to the T1w segmented structures in some simulations. Veins were identified using thresholded, high-pass filtered MRI venograms (Gaussian filter, full-width half-maximum: 15 mm, threshold: 15% of maximum value). These venograms were then aligned to the T1w image and, when included in the model, the optical properties of the veins were set to ${\mu }_{\mathrm{a}}=0.3{\mathrm{mm}}^{-1}$, ${\mu }_{\mathrm{s}}^{\prime }=1{\mathrm{mm}}^{-1}$ (obtained from Ref. 38). When included in the model, the vein compartment substituted some of the CSF layer.

The FEM software NIRFAST³⁹ was used to model the light propagation through the investigated medium, both for homogeneous finite geometries and heterogeneous heads. NIRFAST was used to compute the boundary data for a given optode distribution when simulations were employed.

2.3.

EAC Values

EAC values, defined in Eq. (1), were computed using a more general and accurate solution of the diffusion equation³³ as a function of optical parameters when a simple semi-infinite geometry is employed (considering source–detector distances from 20 up to 60 mm). The linear assumption was tested by performing regression analyses on simulated data based on a set of different baseline optical properties covering most of the range commonly observed in human head tissues.

EAC is computed based on interoptode Euclidean distances, but its estimated value can also be affected by surface geometry. To test the effect of a curved surface, EAC was estimated using FEM. A spherical medium was generated with constant optical properties (${\mu }_{\mathrm{a}}=0.01{\mathrm{mm}}^{-1}$, ${\mu }_{\mathrm{s}}^{\prime }=1{\mathrm{mm}}^{-1}$) and different curvatures (simulated by sphere radiuses ranging from 50 to 150 mm). The bulk absorption and reduced scattering coefficients were chosen to represent commonly accepted average values of optical properties of the head in the NIR range.³⁸ The range of sphere radiuses was chosen to approximate adult human head curvatures. Different light intensities were computed by moving a source relative to a fixed detector on the sphere surface. The linear assumption was tested by computing regression analyses for simulations based on a set of surface curvatures covering the range likely to be observed in normal adult human heads.

2.4.

Construction of Topographic EAC Maps from Two-Dimensional Optode Array

In order to quantify local variations in EAC using a topographic two-dimensional (2-D) approach, we developed an algorithm that computed EAC for each channel (source–detector pair) using subsets of channels in close proximity (neighbors) to the channel being estimated. A flowchart of the procedure is reported in Fig. 2; a more formal description of the algorithm is as follows:

Fig. 2

Flowchart of the procedures used for EAC mapping.

The number of channels considered for each EAC estimate can be adjusted by modifying the neighborhood radius $\rho$, thereby affecting (in opposite direction) the accuracy of EAC estimation and the resolution of the topographic image. The effect of varying $\rho$ was investigated in some of the analyses.

The performance of this image reconstruction process was tested using simulated FEM data. Simulations were based on a homogeneous slab medium ($200\mathrm{mm}\times 200\mathrm{mm}\times 200\mathrm{mm}$, ${\mu }_{\mathrm{a}}=0.01{\mathrm{mm}}^{-1}$, ${\mu }_{\mathrm{s}}^{\prime }=1{\mathrm{mm}}^{-1}$) with a cylindrical inhomogeneity (10 mm diameter parallel to the investigation surface, up to 40 mm depth, ${\mu }_{\mathrm{a}}=0.04{\mathrm{mm}}^{-1}$) at the center of the cube face. One face of the cube was covered by a square grid of optodes [alternating source and detector fibers with a 15-mm minimum interfiber distance, see Fig. 3(c)]. Only channels with source–detector distances $>20$ and $<60\mathrm{mm}$ were used for the analysis. The accuracy and resolution capabilities of the algorithm were estimated using different neighborhood radiuses.

Fig. 3

(a) Analytic solutions for light reflectance [multiplied by the $r^2$ term, left-hand side of Eq. (1)] as a function of source–detector distance for a semi-infinite homogeneous medium (semilogarithmic plot). Different curves represent different absorption coefficients of the medium (from ${\mu }_{\mathrm{a}}=0.005{\mathrm{mm}}^{-1}$ to ${\mu }_{\mathrm{a}}=0.04{\mathrm{mm}}^{-1}$, in steps of $0.001{\mathrm{mm}}^{-1}$, arrow depicting ascending values). The reduced scattering coefficient was fixed at ${\mu }_{\mathrm{s}}^{\prime }=1{\mathrm{mm}}^{-1}$. (b) Analytic solutions for light reflectance [multiplied by the $r^2$ term, left-hand side of Eq. (1)] as a function of source–detector distance for a semi-infinite homogeneous medium (semilogarithmic plot). Different curves represent different reduced scattering coefficients of the medium (from ${\mu }_{\mathrm{s}}^{\prime }=0.2{\mathrm{mm}}^{-1}$ to ${\mu }_{\mathrm{s}}^{\prime }=2{\mathrm{mm}}^{-1}$, in steps of $0.1{\mathrm{mm}}^{-1}$, arrow depicting ascending values). The absorption coefficient was fixed at ${\mu }_{\mathrm{a}}=0.01{\mathrm{mm}}^{-1}$. (c) EAC estimates obtained using a multidistance regression approach using interoptode distances between 20 and 60 mm. The EAC estimates are reported as a function of the absorption and reduced scattering coefficients (from ${\mu }_{\mathrm{a}}=0.005{\mathrm{mm}}^{-1}$ to ${\mu }_{\mathrm{a}}=0.04{\mathrm{mm}}^{-1}$ and from ${\mu }_{\mathrm{s}}^{\prime }=0.2{\mathrm{mm}}^{-1}$ to ${\mu }_{\mathrm{s}}^{\prime }=2{\mathrm{mm}}^{-1}$).

2.5.

EAC Estimation in Recorded Data

3.1.

EAC Values

A simulation was run to determine that Eq. (1) is valid across a wide range of optical parameters as applied to a homogeneous semi-infinite medium and to estimate the possible EAC values found in vivo.

Figures 3(a) and 3(b) report the analytic solutions for light reflectance [multiplied by the $r^2$ term, left-hand side of Eq. (1)] as a function of distance between sources and detectors located on the medium surface. Figure 3(a) reports results obtained with a fixed reduced scattering parameter (${\mu }_{\mathrm{s}}^{\prime }=1{\mathrm{mm}}^{-1}$) and variable absorption coefficients (from ${\mu }_{\mathrm{a}}=0.005{\mathrm{mm}}^{-1}$ to ${\mu }_{\mathrm{a}}=0.041{\mathrm{mm}}^{-1}$, with $0.002{\mathrm{mm}}^{-1}$ increments). Figure 3(b) reports results obtained with a fixed absorption parameter (${\mu }_{\mathrm{a}}=0.01{\mathrm{mm}}^{-1}$) and variable reduced scattering coefficients (from ${\mu }_{\mathrm{s}}^{\prime }=0.2{\mathrm{mm}}^{-1}$ to ${\mu }_{\mathrm{s}}^{\prime }=2{\mathrm{mm}}^{-1}$, with $0.1{\mathrm{mm}}^{-1}$ increments).

Figure 3(c) reports the EAC estimates obtained by regressing the log light reflectance (multiplied by the $r^2$ term) as a function of interoptode distance, based on interoptode distance values between 20 and 60 mm. A minimum distance of 20 mm was chosen in order to make the assumption $r\gg \frac{1}{{\mu }_{\mathrm{eff}}}$ hold for all the optical properties considered. The EAC estimates are reported as a function of the absorption and reduced scattering coefficients (from ${\mu }_{\mathrm{a}}=0.005{\mathrm{mm}}^{-1}$ to ${\mu }_{\mathrm{a}}=0.04{\mathrm{mm}}^{-1}$ and from ${\mu }_{\mathrm{s}}^{\prime }=0.2{\mathrm{mm}}^{-1}$ to ${\mu }_{\mathrm{s}}^{\prime }=2{\mathrm{mm}}^{-1}$). EAC values can vary in vivo between the values of 0.05 and $0.5{\mathrm{mm}}^{-1}$. The maximum error of the estimated EAC based on Eqs. (1) and (2) across the optical values considered was $\sim 1\%$ with an intercept standard deviation of 0.2. In order to test the linearity of the log light intensity as a function of distance, correlations were also computed between interoptode distances (20 to 60 mm) and the diffuse optical reflectance (adjusted for the $r^2$ term) for all the combinations of optical parameters considered. The average absolute correlation obtained was 0.999 ($\mathrm{SD}=0.0005$). The lowest absolute correlation obtained was 0.997. These results strongly support the linearity of the relationship between source–detector distance and the log of the intensity of the light reaching the detector [corrected by the $r^2$ term, Eq. (1)]. These data also support the claim that Eq. (1) is valid for semi-infinite, homogeneous media across the source–detector, absorption, and reduced scattering values typically found in diffuse optical imaging applications.

Note that, although the head structure is optically not homogeneous, it has been extensively proven using the FDMD method that the linearity assumption holds locally when in vivo head measurements are performed, both in infants and adults.^16,44,45

A second simulation was performed in order to estimate the effect of a curved, finite space (such as the head) on EAC estimates. Figure 4(a) reports the reflectance (multiplied by the $r^2$ term) derived from FEM simulations as a function of Euclidean interoptode distance in a spherical geometry (${\mu }_{\mathrm{a}}=0.01{\mathrm{mm}}^{-1}$ and ${\mu }_{\mathrm{s}}^{\prime }=1{\mathrm{mm}}^{-1}$). Different colors represent reflectance values obtained with different sphere radiuses (from $r=50\mathrm{mm}$ to $r=170\mathrm{mm}$, with 20 mm increments). The continuous line depicts the reflectance obtained in a flat geometry. Figure 4(b) reports the percent distance in EAC estimates (compared to flat-surface estimates) as a function of the sphere radius. Smaller radiuses (high curvature) decrease the estimated EAC up to a few percentage points. The average absolute correlation obtained between the reflectance [corrected $r^2$ term, Eq. (1)] and the Euclidean interoptode distance across all sphere radiuses was 0.9965 ($\mathrm{SD}=0.00055$). The weakest absolute correlation obtained was 0.9955. No clear effect of sphere radius on the linearity of the problem was found (within FEM noise level).

Fig. 4

(a) Simulated reflectance [multiplied by the $r^2$ term, left-hand side of Eq. (1)] as a function of interoptode distance obtained using FEM in a spherical geometry (${\mu }_{\mathrm{a}}=0.01{\mathrm{mm}}^{-1}$, ${\mu }_{\mathrm{s}}^{\prime }=1{\mathrm{mm}}^{-1}$, semilogarithmic plot). Different colors represent different sphere radiuses (from $r=50\mathrm{mm}$ to $r=170\mathrm{mm}$, in steps of 20 mm, arrow depicting ascending sphere radius). The continuous line depicts the reflectance obtained with a flat geometry. (b) Percent distance in EAC estimates as a function of sphere radius.

These data indicate that the presence of a curved surface produces a small but systematic bias in the EAC estimates obtained from Eqs. (1) and (2), with a tendency to underestimate EAC for regions of higher curvature (smaller radius). This effect, however, is less than 10%, even for the spheres with the highest curvature. Note that the radius of an adult human head is of the order of 80 to 90 mm whereas that of a newborn is of the order of 45 to 50 mm. An in-depth evaluation of curvature effects as a function of the population considered (from neonates to adults) can be found in Ref. 46.

3.2.

Topographic Two-Dimensional Image Reconstruction

Figure 5(a) reports the topographic image of EAC obtained with a neighborhood radius of 20 mm, based on simulated data. Figure 5(b) reports the same image with a neighborhood radius of 40 mm. The white circles represent the original inhomogeneity position projected onto the medium’s surface. The procedure used to derive these topographic images was able to correctly retrieve the position of the inhomogeneity, with a resolution dependent on the neighborhood sphere radius. Note that the spatial localization power was tested in these images. Consider that this particular simulation was affected by the partial volume effect, because of the small inhomogeneity used and the low resolution of the procedure. This is evident when comparing the two images obtained with different neighborhood radiuses [Fig. 5(b) presents a less intense reconstructed inhomogeneity when compared to Fig. 5(a)]. A bigger inhomogeneity is needed to test the accuracy in reconstructed EAC values (refer to Fig. 8 presented later in this paper). Figure 5(d) reports the full width half maximum (FWHM) resolution estimate of the image reconstruction procedure as a function of the neighborhood radius: the resolution of the image is positively correlated with the neighborhood radius. Figure 5(e) reports the average number of channels and related standard deviation that were involved in the computation of EAC as a function of neighborhood radius. Note that the standard deviation bars are nonzero because of the finite geometry employed. In a semi-infinite medium with an infinite, regularly spaced optode montage on the surface, there would be no variability in the number of channels for each EAC computation at a given neighborhood radius. This number is positively correlated with (the square of) the neighborhood radius (when a fixed optode density is considered), and it is an important parameter to take into account when applying the algorithm to recorded data. For recorded data applications, including more channels for each EAC computation leads to better EAC estimates (as it will be discussed later, see Fig. 7), generally losing spatial resolution because of the increased neighborhood radius. Another way of increasing the number of channels considered in each EAC computation is to increase the optode’s density, which will increase precision without sacrificing spatial resolution. Consider that for this simulation, the optode density chosen was particularly high, and as such it can be considered a best-case scenario. Note, finally, that the closest source–detector pairs (15 mm distance) were not used for EAC computation.

Fig. 5

(a and b) Topographic images of EAC obtained with simulated FEM data on a cubic medium ($200\mathrm{mm}\times 200\mathrm{mm}\times 200\mathrm{mm}$, 15 mm minimum interoptode distance grid, ${\mu }_{\mathrm{a}}=0.01{\mathrm{mm}}^{-1}$, ${\mu }_{\mathrm{s}}^{\prime }=1{\mathrm{mm}}^{-1}$, with a cylindrical 10-mm-diameter absorption inhomogeneity, ${\mu }_{\mathrm{a}}=0.04{\mathrm{mm}}^{-1}$, up to 40 mm depth). (a) Topographic image obtained with a neighborhood radius of 20 mm. The white circle indicates the actual position and shape of the inhomogeneity. (b) Topographic image obtained with a neighborhood radius of 40 mm. The white circle indicates the actual position and shape of the inhomogeneity. (c) Optode array used for the simulations (15 mm minimum interoptode distance grid). (d) FWHM of the image reconstruction procedure as a function of the neighborhood radius. (e) Average number of channels, and related standard deviation, involved in the computation of each channel’s EAC as a function of the neighborhood radius.

3.3.

Estimating EAC from Recorded Data

Two main assumptions are considered when computing EAC from recorded data. First, the SNR of the recorded data is mainly influenced by quantum/shot noise and can provide a good estimate of light intensity in the absence of quantitative measures of intensity. Second, the noise introduced by O-S coupling can be diminished by considering a sufficient number of channels for each EAC computation and by statistical correction procedures.

Figure 6(a) reports a typical spectrum of the logarithm of the SNR for signals recorded in vivo. The spectrum can be divided into two regions. In the first region (lower frequency ranges, below 10 Hz), spectrum peaks can be related to physiological fluctuations of the signal (such as hemodynamics, respiration, Mayer waves, and pulse waves). Above 10 Hz, the spectrum is practically flat, and its power is mainly related to quantum/shot noise if a sufficient number of photons are detected. In fact, Fig. 6(b) reports the ${\mathrm{SNR}}^2$ (computed above 10 Hz) for a particular channel during an in vivo measurement where the integration time of the detector was manipulated (between 1 and 5 ms). At a constant source power, the integration time increases the number of photons detected, and the high linearity of the ${\mathrm{SNR}}^2$ versus integration time ($r=0.99$, $p<0.001$) supports the assumption of shot noise prevalence in the frequency band considered [Eq. (3)].

Fig. 6

(a) A typical spectrum of the logarithm of the SNR for signals recorded in vivo. The spectrum is practically flat above 10 Hz, and its power is mainly related to quantum/shot noise. (b) ${\mathrm{SNR}}^2$ (computed above 10 Hz) for a particular channel during an in vivo measurement, as a function of the integration time of the detector. (c) Logarithm of the SNR (multiplied by the $r^2$ term) for each channel as a function of interoptode distance, obtained from a set of data recorded in vivo. The vertical line shows the maximum distance used in the computation. (d) Logarithm of the SNR (multiplied by the $r^2$ term) as a function of interoptode distance obtained from good channels (below 60 mm interoptode distance). Note that the SNR values do not lay on a straight line for both tissue inhomogeneities and O-S coupling. The inset histogram reports the relative frequency of the residuals of a linear interpolation considering all the good channels. (e) Logarithm of the SNR (multiplied by the $r^2$ term) as a function of interoptode distance obtained from good channels (below 60 mm interoptode distance) after O-S coupling correction.

Figure 6(c) reports the logarithm of the SNR (multiplied by the $r^2$ term) for each channel as a function of interoptode distance, obtained from a set of data recorded from one subject. The vertical line shows the maximum distance used in our computation. The maximum distance separates the scatterplot into two regions. In the left region, the SNR decreases linearly with distance (with variability due to the heterogeneity of the human head and O-S coupling). In the right region, the SNR plateaus because light intensity is very low and the shot noise is small compared to other sources of noise. Equation (3) is assumed to be valid for channels with distances below the maximum distance chosen. Therefore, only these channels are included in the computation of EAC. An outlier identification procedure further identified channels with very low levels of light and removed them from the computation algorithm. Specifically, after the identification of the maximum usable interoptode distance, the data were linearly interpolated and the residuals were examined. Channels that had a residual value greater than five times the SD of the residual distribution were excluded from further computations. Figure 6(d) reports the logarithm of the SNR as a function of interoptode distance obtained from good channels (below 60 mm interoptode distance). Note that the SNR values do not lie on a straight line due to both tissue inhomogeneities and O-S coupling. However, the data provide good homoscedasticity, further supporting the implemented statistical procedure for noise dampening. Assuming O-S coupling as the only source of noise, we estimated the O-S coupling for the five in vivo measurements that presented similar results as the ones reported in Fig. 6(d). O-S coupling was estimated from the residual of a linear interpolation of the data [as shown in Fig. 6(d)] and their Gaussian properties were tested with a Kolmogorov–Smirnov test (average noise $\mathrm{SD}=1.5$, all $p’\mathrm{s}<0.01$). This procedure can be considered to represent a “worst case scenario” since it yields overestimation of both intensity and departure from normality of O-S coupling.

Figure 6(e) reports the logarithm of the SNR as a function of interoptode distance obtained after O-S coupling correction. As expected, a reduction in the variance of the data is present in this figure. Figure 7(a) reports the standard error of the slope as a function of the number of independent samples used with varying residual sample noise [Eq. (4)]. Since the standard error in Eq. (4) depends also on the independent variable values (interoptode distances in this particular case), we considered a homogeneous distribution of distances between 20 and 60 mm in this analysis. The thicker line in the graph represents a residual sample noise of $\mathrm{SD}=1.5$, which is similar to the level of noise found in vivo.

Fig. 7

(a) Standard error of the slope as a function of the number of independent sampled used with varying residual sample noise [Eq. (4)]. Since the standard error in Eq. (4) depends also on the independent variable values (interoptode distances in this particular case), a homogeneous distribution of distances between 20 and 60 mm was used for this analysis. The thicker line represents a residual sample noise of $\mathrm{SD}=1.5$, which is similar to the level of noise found in vivo. (b) EAC error (%) as a function of the number of channels averaged in the process, separately for data uncorrected and corrected for optode effects.

Note that the standard error of the slope represents the error of the implemented EAC mapping. This error needs to be compared to the EAC values generally measured in vivo (the actual slope of the regression), and it is generally at least an order of magnitude smaller when a sufficient number of independent samples (channels) are employed. However, the measurements are not completely independent (each optode can be involved in multiple channels), and a further optode efficiency correction procedure can be introduced [Eqs. (5)–(8)]. To quantify the effects of channel numerosity and of the correction procedure for suppressing the effect of O-S coupling on the accuracy of EAC computation, we ran a simulation where Gaussian noise was added to the log intensity data of each source and detector. Figure 7(b) shows the EAC error as a function of the number of channels averaged in the process, separately for uncorrected and corrected data. As the figure shows, the EAC estimate error was reduced to 2% to 3% when a large number ($\sim 20$ or more) of neighboring channels and the correction algorithm were employed in the EAC computation. As stated before, in order for the O-S coupling correction procedure to work, it is important that the parameter $k_0$ [Eq. (5)] does not change as a function of tissue optical properties. We determined that, for a homogeneous medium, the parameter $k_0$ across a wide range of optical properties had a $\mathrm{SD}=0.2$. Note that this variability is lower than that of the O-S coupling introduced in the simulations (which we estimated from in vivo data, $\mathrm{SD}=1.5$). Thus, the correction procedure should improve the EAC estimation reliability at a given number of channels when a flat homogeneous medium is investigated.

Both different medium geometries and tissue inhomogeneities can affect the variability of the parameter $k_0$. This aspect is difficult to address for in vivo head studies since the actual properties of the different head tissues are not exactly known. However, as a first step we simulated light intensity signals and estimated the intercept parameter $k_0$ on heterogeneous head models using the optode montages used in vivo [Fig. 10(b)], segmentations of the T1w images of the five subjects that underwent the in vivo experiment, and a priori optical parameters. For these particular simulations, the parameter $k_0$ was found to have $\mathrm{SD}=0.8$. Note that this variability is still lower than the level of variability found in vivo.

For practical purposes, the use of multiple channels for each EAC computation (with a good combination of sources and detectors) is mandatory in order to suppress O-S coupling variability and to obtain a reliable EAC estimate. The O-S coupling correction can be applied to further improve the reliability of the EAC estimate when the intercept factor $k$ can be accurately estimated for each optode and the optical parameter changes over space are assumed to be within reasonable ranges.

We then investigated the validity of using the data SNR (to estimate light intensity), multiple combinations of sources and detectors (using a regression approach), and O-S coupling corrections to estimate EAC in phantom data. Specifically, we compared the EAC maps computed from data recorded from the surface of a milk tank containing a known inhomogeneity to maps obtained from FEM simulated data (keeping the montage, medium and inhomogeneity constant for the two approaches). Figure 8(a) (left image) reports the optode montage employed in the simulated and real phantom data, together with the channel centers used to compute EAC maps. Figure 8(a) (right image) reports the optode montage employed, highlighting the sources (reported in yellow), the detectors (reported in blue), the channels (represented as magenta lines), and the center of the channels (reported as small black x’s) involved in the computation of the EAC in the position marked by the large gray “X,” using a neighborhood radius of 30 mm. Twenty-four channels were employed in the computation of the EAC in that particular point. Figure 8(b), top graph, reports the log of the intensity (multiplied by the $r^2$ term), computed using FEM simulations, as a function of interoptode distance for the channels involved in the computation of EAC in the black $x$ point [Fig. 8(a), right image]. Figure 8(b), bottom graph, reports the log of the ${\mathrm{SNR}}^2$ (multiplied for the $r^2$ term), for real data in the same location. Although the real data showed an increased level of noise, the EAC estimates for the real and simulated data in the location of interest were identical ($\mathrm{EAC}=0.24$). The estimated standard error for the real measurement was $\mathrm{SE}=0.02$ [based on Eq. (4)].

Fig. 8

(a) Left image: Optode montage used for the simulated and real phantom data, together with the channel centers (midpoints). Right image: Optode montage with sources (yellow), detectors (blue), channels (represented as magenta lines), and the centers of each channel, reported as black x’s, involved in the computation of the EAC at the position indicated by the thick gray X, using a neighboring radius of 30 mm. (b) Top image: Logarithm of the intensity (multiplied by the $r^2$ term), computed using FEM simulations, as a function of interoptode distance, for the channels involved in the computation of EAC at the gray X point [Fig. 8(a), right image]. Bottom image: Logarithm of the ${\mathrm{SNR}}^2$ (multiplied for the $r^2$ term), for real data at the same location. Although increased noise was found for in vivo data, the EAC estimates at the location of interest were identical ($\mathrm{EAC}=0.24$). The estimated standard error for the real measurement was $\mathrm{SE}=0.02$ [based on Eq. (4)]. (c) Top image: EAC maps obtained using light intensity in FEM simulated data. Bottom image: Topographic map of EAC obtained using the signal-to-shot noise ratio of recorded data. The white rectangle represents the original inhomogeneity position and shape (50 mm thick, 40 mm wide, and 100 mm long). Both images were obtained using a neighborhood sphere radius of 30 mm (average number of channels for each $\text{EAC computation}=30$, $\mathrm{SD}=5$).

Figure 8(c) (top) reports a topographic map of EAC obtained using light intensity in FEM simulated data, whereas Fig. 8(c) (bottom) presents the topographic map of EAC obtained using the signal-to-shot noise ratio of recorded data. The white rectangle represents the shape and position of the inhomogeneity. Both images were obtained using a neighborhood sphere radius of 30 mm. The similarities are clear both qualitatively and quantitatively. The spatial correlation between the two images, computed across all voxels, was $r=0.9814$ ($p<0.0001$). Note that this analysis also demonstrated the ability of the procedure to correctly quantify the EAC across space (milk: ${\mu }_{\mathrm{a}}=0.005{\mathrm{mm}}^{-1}$ and ${\mu }_{\mathrm{s}}^{\prime }=1{\mathrm{mm}}^{-1}$, $\mathrm{EAC}=0.122$; phantom: ${\mu }_{\mathrm{a}}=0.02{\mathrm{mm}}^{-1}$ and ${\mu }_{\mathrm{s}}^{\prime }=1{\mathrm{mm}}^{-1}$, $\mathrm{EAC}=0.24$, at 830 nm wavelength).

Figure 9 explores the different image reconstruction capabilities of the procedure as a function of the neighborhood radius employed (from 20 to 60 mm, 10 mm steps from top to bottom), and the type of data (from left to right, FEM simulated, real data uncorrected for the optode effect, and real data corrected for the optode effect). Clear similarities can be noted between the simulated FEM data and the real data (particularly for the corrected ones), even for the smaller neighborhood radiuses. The location and intensity of the optical inhomogeneity are correctly identified using real data. Moreover, a partial volume effect can be noted when the neighborhood radius increases.

Fig. 9

Different images reconstructed using the developed EAC estimation procedure. From top to bottom images differed for the neighborhood radiuses employed (from 20 to 60 mm, in 10 mm steps). From left to right columns indicate different data types (FEM simulated data, real data, and real data corrected for optode effects, respectively).

3.4.

Full-Head EAC Maps from Human Subjects

Figure 10(a) reports topographic EAC images (based on 830 nm light) projected to the axial surface of structural MRI images of five subjects. Very similar images were obtained using 690 nm light (average image correlation $r=0.85$, $\mathrm{SD}=0.10$) with an average small increase in EAC values (5% compared to the 830 nm). These images were obtained using an optode array covering most of the scalp surface [Fig. 10(b)] and a neighborhood radius of 30 mm (average number of channels for each EAC computation: $N=28$, $\mathrm{SD}=6$). Reliability tests were performed on each subject by generating the EAC images twice, employing half of the optodes available for each image computation (with optodes randomly chosen, and no optode shared between the two images). Although the split-half data showed decreased SNR (because of a decreased number of optodes), the root mean square error between the two images across the five subjects was 3%, $\mathrm{SD}=1\%$. For each subject we also computed the correlation between the EAC values for each of the two images (the correlation was computed across voxels). The average value of this correlation across subjects was $r=0.86$, $\mathrm{SD}=0.06$. Note that this correlation was based on images computed using different optodes. Nevertheless the intrinsic low resolution of the images and the high number of channels employed allowed for this procedure to provide a lower-bound estimator of the reliability of the EAC mapping method. There are clear similarities in the patterns of EAC across participants. In particular, dorsal and occipital regions appear to have higher EAC values compared to lateral frontal and temporal regions. Intersubject consistency analyses were performed on the images by correlating each subject’s image with the average of the other four subjects (a single correlation was measured for each subject using all the available voxels). The rationale of this procedure was to show spatial similarities between EAC maps. The correlations ranged from 0.4 to 0.85 (all $p’\mathrm{s}<0.001$) and are reported in Fig. 10(c). The images were coregistered based on their center. Due to small differences in the shapes of the topographic EAC maps, only pixels where EAC was computed for all the images could be used for the spatial correlation analyses. No warping procedure was applied. In fact, on average, only 6% of the pixels (3% SD) were removed from each image during this procedure. Given the relatively low spatial resolution of the method, we would not expect major changes in the spatial correlations if a warping procedure were applied.

Fig. 10

(a) Topographic EAC maps obtained from five adult subjects (18 to 30 years), with a neighborhood radius of 30 mm. The results were obtained using 830-nm light. (b) Example of the optode montage used for optical recording, overlaid over a subject’s structural MRI. (c) Image correlations between each subject’s EAC map and the average map of the other four subjects.

We hypothesized that the major changes in EAC across the head might be due to the presence of large veins in superficial intracranial areas (venous sinuses). Figures 11(a)–11(c) show EAC maps (projected to the axial surface) based on (a) measurements recorded from a participant; (b) simulated measurements based on segmentation of the same subject’s T1w image; and (c) simulated measurements where veins (obtained from a segmented venogram from the same subject) were added to the T1w segmentation. Figure 11(d) shows the spatial correlation of the two simulated results with the actual measurements for each of the five subjects. Positive and statistically significant ($p<0.01$) spatial correlations were obtained for each subject when veins were added to the model, whereas negative spatial correlations were obtained when only skull/scalp, CSF, and white and gray matter were considered. Thus, considering only the tissue layers available with a standard T1w image results in largely inaccurate estimates (based on a spatial correlation metric) of local variations in EAC across the head in our dataset. This problem is greatly ameliorated by the inclusion of the vein compartment in the forward model.

Fig. 11

(a) Recording helmet seen from the top (inset); EAC map (top view) obtained from one participant (S1); (b) MR T1w segmentation of the same subject (inset; purple: scalp and skull; turquoise: CSF; red: gray matter; yellow: white matter); simulated EAC map (top view) based on a standard four-layer model of optical tissue properties, derived from the segmented T1w image of the same subject; (c) MR T1w segmentation combined with the MR-derived venogram of the same subject, superimposed over a gray matter rendering (inset); Simulated EAC map (top view) in which the veins have been added to the four-layer model. (d) Spatial correlations between recorded and simulated EAC maps for each subject. Positive and statistically significant ($p<0.001$) correlations were obtained for each subject when veins were added to the model. Negative correlations were instead obtained when only the T1w-segmented structures were considered.