2.1.

Forward and Inverse Problems

The goal of image reconstruction (inverse problem) is the recovery of properties of the investigated object for a given recording pattern. For functional NIR DOT, the property of interest is the optical absorption perturbation $\mathrm{\Delta }{\mu }_{\mathrm{a}}(\stackrel{\rightharpoonup }{r})$ over time at each position in space ($\stackrel{\rightharpoonup }{r}$) using measurements of light intensity from the head surface, from which changes in the concentration of oxy- and deoxyhemoglobin can be computed (provided that multiple wavelengths are used). Previous research has shown that NIR light propagation in highly scattering biological tissue (such as the human head, for which the reduced scattering coefficient is ${\mu }_{\mathrm{s}}^{\prime }\gg {\mu }_{\mathrm{a}}$), is well approximated by the diffusion equation.³¹ According to this equation, the diffusion of light in tissue depends on its baseline optical properties for each location ${\mu }_{\mathrm{a}}(\stackrel{\rightharpoonup }{r})$, ${\mu }_{\mathrm{s}}^{\prime }(\stackrel{\rightharpoonup }{r})$. In complex media such as the head, the diffusion model can be solved using FEM,³² in which continuous variations of optical parameters in space are approximated by discrete changes occurring at finite nodes located at small distances from each other. Importantly, when using FEM to solve the diffusion equation, the optical parameters for each node have to be known. This generates a circularity, because in order to compute the absorption coefficient change at a particular location $\mathrm{\Delta }{\mu }_{\mathrm{a}}(\stackrel{\rightharpoonup }{r})$, we must first know the base value of the absorption coefficient ${\mu }_{\mathrm{a}}(\stackrel{\rightharpoonup }{r})$ for all locations in the medium, including the one being estimated (a similar problem exists for the scattering coefficient). Thus, the problem is not linear:³³ a change in the recorded pattern due to $N$ perturbations occurring at the same time in different positions is not equal to the sum of the changes in the recording pattern due to each perturbation occurring separately. This issue is important when the perturbations are large and require computing absolute optical parameters. In this case, iterative processes need to be employed. Multiple algorithms have been developed that rely on different iterative procedures and regularizations.^33–35 They are generally applied to frequency-domain data and show promising results.

In the case of relative changes due to physiological factors, this apparent conundrum can be solved by considering that the changes in absorption (or scattering) related to brain activity are small with respect to the base value (1% to 5% for oxygenation, with the expected changes in scattering being even smaller) and can therefore be ignored as a first approximation when describing the diffusion of light through tissue. The issue of estimating local changes in $\mathrm{\Delta }{\mu }_{\mathrm{a}}(\stackrel{\rightharpoonup }{r})$ can then be isolated from that of how light diffuses (linearization). This means that we can assume that the baseline optical properties of the tissue are not sufficiently changed as a function of brain activity to require us to update the light diffusion model to include such changes. We can therefore estimate a light diffusion model only once without having to iterate the process to include the effects of the changes due to local perturbations to the basic model. This linearized model is mathematically described by a matrix (sensitivity, or Jacobian, matrix) describing the effects on the light intensity measured by optodes located on the medium surface as a function of unitary absorption perturbations at any location inside the medium (forward model):

2.2.

Computation of the Forward Light Model: Finite-Element Method

In our forward model computation, a fine mesh (maximum tetrahedral $\text{volume}=2{\mathrm{mm}}^3$) was generated for FEM using the software iso2mesh³⁶ for both slab geometry (used for the simulated and actual milk-tank data) and heterogeneous models of the head (used for the simulated and actual brain imaging data). Heterogeneous models require segmentation of a structural MR image of the subject’s head into skull and scalp, cerebrospinal fluid (CSF), white matter, and gray matter, which was performed using statistical parametric mapping (SPM) functions applied to T1 structural images.³⁷ Baseline optical properties (absorption parameter ${\mu }_{\mathrm{a}}$, reduced scattering ${\mu }_{\mathrm{s}}^{\prime }$, and refraction index $\eta$) of the tissues at the relevant wavelengths (830 nm) were taken from Tian and Liu.²⁷ Given the similarities of optical properties of skull and scalp and the difficulty of separating them in T1 images, we decided to attribute the same optical property values to both structures. Specifically, the values were set equal to those reported for the skull (see Table 1 for the actual optical values used).

Table 1

Optical properties for different types of tissues (830 nm wavelength, from Tian and Liu27).

The FEM software NIRFAST³³ was used to model light propagation through the medium. NIRFAST was used to compute the boundary data for a given optode montage when simulations were employed. Finally, NIRFAST was also used to compute the sensitivity (Jacobian) matrices for both simulated and real data using the adjoint method.³⁸

2.3.

Inverse Modeling: Energy, Spatially Variant, and Laplacian Regularized Minimum Norm

The inverse problem requires solving Eq. (1) for $\mathrm{\Delta }{\mu }_{\mathrm{a}}(\stackrel{\rightharpoonup }{r})$ based on a known $\mathrm{\Delta }\mathrm{OD}$ (the change in optical density measured at each channel) and J (the Jacobian matrix). However, ${\mathrm{J}}_{[{\mu }_{\mathrm{a}}(\overrightarrow{r}),{\mu }_{\mathrm{s}}^{\prime }(\overrightarrow{r})]}$ is generally not square, and therefore cannot be inverted (generating the well-known ill-posedness of inverse problems). In order to provide a unique and stable solution to the inverse form of Eq. (1), further constraints need to be made (norms). Least squares estimation and energy regularization are the most common procedures employed (generally using the ${\ell }_2$, or Euclidean, norm).³⁹ In this case, $\mathrm{\Delta }{\mu }_{\mathrm{a}}(\stackrel{\rightharpoonup }{r})$ is estimated by minimizing the function $\mathrm{O}[\mathrm{\Delta }{\mu }_{\mathrm{a}}(\stackrel{\rightharpoonup }{r})]$:

Because the ${\ell }_2$ norm is differentiable and convex, the problem can be solved by differentiating $\mathrm{O}[\mathrm{\Delta }{\mu }_{\mathrm{a}}(\stackrel{\rightharpoonup }{r})]$ with respect to $\mathrm{\Delta }{\mu }_{\mathrm{a}}(\stackrel{\rightharpoonup }{r})$ and setting the equation to zero. The following equation, providing a linearized image reconstruction (ER-DOT), is then obtained:

Note that if $ϵ=0$, the problem simplifies to a classic least squares solution, which would provide a 3-D reconstruction with the highest theoretical resolution. However, in most cases, the matrix resulting from the product ${\mathrm{J}}_{[{\mu }_{\mathrm{a}}(\overrightarrow{r}),{\mu }_{\mathrm{s}}^{\prime }(\overrightarrow{r})]}{{\mathrm{J}}_{[{\mu }_{\mathrm{a}}(\overrightarrow{r}),{\mu }_{\mathrm{s}}^{\prime }(\overrightarrow{r})]}}^{\mathrm{T}}$ is almost singular (i.e., the determinant is very close to zero) and the inversion of this matrix leads to unstable results (this practically leads to very noisy 3-D images)—hence, the need for $ϵ>0$. Thus, the parameter $ϵ$ provides reliable numerical results by trading spatial resolution in data fitting for stability. The increased stability is reached by forcing the unknown to be small.

The optimal value of $ϵ$ can be found using statistical methods^40,41 or empirically. In brain DOT, $ϵ$ is normally found empirically as a value proportional to the maximum of ${\Vert {\mathrm{J}}_{[{\mu }_{\mathrm{a}}(\overrightarrow{r}),{\mu }_{\mathrm{s}}^{\prime }(\overrightarrow{r})]}\Vert }^2$.^12,22,23 This means that $ϵ$ is tuned to a value proportional to the squared Euclidean norm of the Jacobian for the node where it is the biggest. The practical consequence of this is that the algorithm finds solutions that minimize the mean squared $\mathrm{\Delta }{\mu }_{\mathrm{a}}(\stackrel{\rightharpoonup }{r})$ required, accounting for the observed effects (minimum norm).

As mentioned above, the minimum norm procedure tends to favor lower-energy superficial solutions over larger-energy deep solutions. Several procedures have been developed to address this issue. They mainly rely on weighting the Jacobian as a function of space. Probably the most common procedure in brain DOT is the spatially variant regularization (SVR-DOT) introduced by Culver et al.²² Thus, here, we employ the SVR-DOT algorithm, where $\mathrm{\Delta }{\mu }_{\mathrm{a}}(\stackrel{\rightharpoonup }{r})$ is estimated by minimizing the function $\mathrm{O}[\mathrm{\Delta }{\mu }_{\mathrm{a}}(\stackrel{\rightharpoonup }{r})]$:

The linearized image reconstruction (SVR-DOT) is obtained as

In order to dampen the effect of decreased sensitivity to deep phenomena, the regularization needs to be increased when the Euclidean norm for a particular node is high and vice-versa. In a general approach, $W$ is a diagonal matrix equal to the $\Vert {\mathrm{J}}_{[{\mu }_{\mathrm{a}}(\overrightarrow{r}),{\mu }_{\mathrm{s}}^{\prime }(\overrightarrow{r})]}\Vert$ for each node at position ($\stackrel{\rightharpoonup }{r}$). It can be computed as

However, this procedure tends to make the solution unstable. In order to address this issue, a spatial regularization parameter $\beta$ needs to be introduced:

Beta is generally chosen to be proportional to the maximum value of ${\Vert {\mathrm{J}}_{[{\mu }_{\mathrm{a}}(\overrightarrow{r}),{\mu }_{\mathrm{s}}^{\prime }(\overrightarrow{r})]}\Vert }^2$. Tuning the parameter using $\beta$ produces a trade-off between stability of the linear inversion and ability to reconstruct deep phenomena.

As mentioned above, a possible alternative method to the energy regularization procedure that is suited to correct for the shallow-depth bias is the Laplacian regularization method. This method may, however, introduce an opposite bias toward deep regions of activation. Therefore, here, we introduce a combination of both approaches, the ELR-DOT algorithm, in order to overcome the issue of depth sensitivity without generating a depth bias. With ELR-DOT, the reconstructed image $\mathrm{\Delta }{\mu }_{\mathrm{a}}(\stackrel{\rightharpoonup }{r})$ is estimated by minimizing the function $\mathrm{O}[\mathrm{\Delta }{\mu }_{\mathrm{a}}(\stackrel{\rightharpoonup }{r})]$:

In the current paper, this parameter was selected based on simulations on a slab geometry model and then used unchanged for all other datasets. In a 3-D space and for a regular cubic grid of nodes, the discrete Laplacian can be approximated using Dirichlet boundary conditions, such that

To speed up computation and to reduce the number of unknowns in the inversion algorithm, the Jacobians were resampled into a cubic grid with a minimum interposition distance $d$ of 5 mm using NIRFAST. With NIRFAST, the Jacobian matrices can be resampled at any desired isovoxel resolution. A 5-mm resolution for the Jacobian was selected as an optimal value, balancing the need for a high level of detail in the description of the spatial properties of the 3-D model while limiting the illposedness of the problem. In particular, the value chosen (5 mm) provides a spatial frequency that is beyond the Nyquist frequency involved in the DOT reconstruction process (spatial resolution of DOT is around 15 mm).¹² For display purposes, 1-mm-voxel images were reconstructed by cubic spline interpolation of $\mathrm{\Delta }{\mu }_{\mathrm{a}}(\stackrel{\rightharpoonup }{r})$ in node space.

2.4.

Metrics

In order to estimate the performance of the ELR-DOT algorithm compared to ER-DOT and SVR-DOT, several different metrics were used. One metric estimated the stability of the algorithm inversion independently of the recorded data. Looking at Eqs. (3), (5), and (10), the linear inversion consists of multiplying recorded data $\mathrm{\Delta }\mathrm{OD}(n)$ by a matrix that we called $G$:

We decided to estimate the algorithm stability by computing the Frobenius norm of matrix $G$ for different regularization parameters and different algorithms. The Frobenius norm is computed as

The Frobenius norm tests for the magnitude of the values of the matrix $G$ and, for a given forward solution, is a test of the inverse procedure stability. When the solution is more stable, lower values of ${\Vert G\Vert }_F$ are obtained and vice-versa.⁴²

Other metrics considered were related to the reconstructed images. When a single perturbation was considered, three metrics were used. The position error (PE) was computed as the difference between the centroid of the known position of a perturbation and the reconstructed one, namely,

For each $\mathrm{\Delta }{\mu }_{\mathrm{a}}(\stackrel{\rightharpoonup }{r})>{0.5}_{\max }[\mathrm{\Delta }{\mu }_{\mathrm{a}}(\stackrel{\rightharpoonup }{r})]$, where ${\stackrel{\rightharpoonup }{r}}_{\text{known}}$ is the known position of the perturbation, $n$ is the number of the nodes, and $\mathrm{\Delta }{\mu }_{\mathrm{a}}(\stackrel{\rightharpoonup }{r})$ is the reconstructed absorption changes. The constraint $\mathrm{\Delta }{\mu }_{\mathrm{a}}({\stackrel{\rightharpoonup }{r}}_i)>{0.5}_{\max }[\mathrm{\Delta }{\mu }_{\mathrm{a}}(\stackrel{\rightharpoonup }{r})]$ was chosen in order to dampen the effect of small noise in the PE estimation.

The estimated depth (ED) of a reconstructed perturbation was computed as

In addition, the full-width half-maximum (FWHM) of the reconstructed perturbation was computed. The average FWHM along the three main axes is reported as

When two perturbations were considered, a value called the resolution parameter (RP) was computed. RP is the ratio between the absorption value at the midway position between the two perturbations and the average value of absorption at the centers of the perturbations:

2.5.

Simulations: Slab Geometry

In order to estimate the optimal combination of energy and Laplacian regularization parameters, we ran a simulation considering a homogeneous cubic medium of $100\times 100\times 100\mathrm{mm}$, having baseline absorption and reduced scattering coefficients of ${\mu }_{\mathrm{a}}=0.01{\mathrm{mm}}^{-1}$ and ${\mu }_{\mathrm{s}}^{\prime }=1{\mathrm{mm}}^{-1}$, respectively. These absorption and reduced scattering coefficients were chosen to represent commonly accepted average values of optical properties in the NIR range⁴³ obtained with in-vivo measurements of biological tissues. A square array of optodes (minimum source–detector distance of 15 mm) was arranged on a face of the cube (Fig. 1). Only source–detector distances within 60 mm were considered meaningful. A spherical absorber with a 4-mm radius and ${\mu }_{\mathrm{a}}=0.011{\mathrm{mm}}^{-1}$ (10% change from bulk value) was embedded at the center of the investigation surface at different depths (centered from 5 to 35 mm from the recording surface). Gaussian noise amounting to 5% of the maximum intensity change was introduced in the simulation as an initial step to assess the effects of noise on the DOT reconstruction algorithm.³³ Intensities were computed for each meaningful source–detector pair for the different conditions. The optimal energy and Laplacian regularization parameters were estimated by minimizing the maximum PE among the different depths considered (we required the algorithm to be accurate for all those depths). The ELR-DOT regularization algorithm was compared to ER-DOT and SVR-DOT as a function of regularization values using the metrics Frobenius norm, PE, and FWHM.

Fig. 1

Homogeneous cubic mesh ($100{\mathrm{mm}}^3$) used for simulations and real data image reconstruction. A square array of optodes, with alternating sources and detectors, was arranged over parallel lines at distances of 15 mm on a face of the cube. This same array was also used for the real (phantom) data collected in the milk tank. For both the simulation and the phantom experiments, an absorption perturbation was positioned at the center of the investigation surface at different depths.

In order to estimate the resolution capabilities of ELR-DOT, two absorbing perturbations ($\text{radius}=4\mathrm{mm}$, ${\mu }_{\mathrm{a}}=0.011{\mathrm{mm}}^{-1}$) were embedded at a depth of 25 mm, while the distance between them was varied by moving their centers from 10 to 40 mm apart (perturbation distance, PD) in 2-mm steps. The RP metric was evaluated for each PD.

2.6.

Simulations: Real-Head Geometry

The robustness of parameter-tuned ELR-DOT to changes in geometry and baseline optical values were tested using a heterogeneous realistic head geometry. The head’s geometry and optical values were derived from T1-weighted MR images obtained in five subjects who were recruited and run using procedures approved by the University of Illinois Institutional Review Board. The T1-weighted MR images were acquired with a Siemens Trio^® 3T scanner using a 3-D MPRAGE (magnetization prepared rapid gradient echo) sequence [$\text{re}\text{pitition\hspace{0.17em}}\text{time}(\mathrm{TR})=1900\mathrm{ms}$, $\text{\hspace{0.17em}inversion time}(\text{TI})=900\mathrm{ms}$, $\text{echo}\text{\hspace{0.17em}time}(\mathrm{TE})=2.32\mathrm{ms}$, $\text{field of view}=230\times 230\times 172.8{\mathrm{mm}}^3$ (sagittal), $\text{matrix size}=256\times 256\times 192$, $\text{flip angle}=9\mathrm{deg}$, $\text{slice thickness}=0.9\mathrm{mm}$]. Single MRI image slices were concatenated and converted into NIfTI format using the dcm2nii conversion tool⁴⁴ called by a MATLAB^® toolbox. Each image was then placed into a left–right (L-R), posterior–anterior (P-A), inferior–superior (I-S) spatial orientation and resampled into 1-mm isovoxels by cubic spline interpolation of the original image.⁴⁵ Segmentation of the subject’s head into skull and scalp, CSF, white matter, and gray matter was performed using SPM functions applied to T1 structural images.³⁷

A fine mesh (maximum tetrahedral $\text{volume}=2{\mathrm{mm}}^3$) was generated for NIRFAST computing using the software iso2mesh.³⁶ A square array of optodes (minimum source–detector $\text{distance}=15\mathrm{mm}$) was simulated on the scalp over the occipital cortex of the subjects (Fig. 2). Only source–detector distances within 60 mm were considered meaningful. An absorbing sphere with $\text{radius}=4\mathrm{mm}$ and ${\mu }_{\mathrm{a}}=0.011{\mathrm{mm}}^{-1}$ (10% change from bulk value) was embedded in the occipital area at different depths from the scalp (center varying from 3 to 33 mm from the surface). Intensities for the different conditions were computed. Gaussian noise amounting to 5% of the maximum intensity change was introduced in the simulation. ELR-DOT, SVR-DOT, and ER-DOT were applied to the simulated recorded intensities, and ED and FWHM metrics were evaluated using optimum regularization values, determined using the homogeneous slab simulations.

Fig. 2

Example of a heterogeneous head mesh (maximum tetrahedral $\text{volume}=2{\mathrm{mm}}^3$) generated for FEM, computed using iso2mesh. (a) Segmentation divisions of skull and scalp, CSF, white matter, and gray matter, performed on T1 structural images. Different optical properties (see Table 1) were attributed to different tissues. (b) Source and detector locations (same geometry as Fig. 1) on the occipital surface of the head.

2.7.

Phantom Data: Slab Geometry

We also tested the ability of the combination of the FEM forward solver with the ELR-DOT algorithm to reconstruct real changes of local absorption in a homogeneous slab medium using reflectance geometry. Because of the known similarity between the scattering properties of milk and human tissue in the NIR spectral range, homogenized skim milk was used as a bulk substance.⁴⁶ The milk was put in a 4-l black tank and was kept at a constant temperature using ice. India ink was added in order to increase the homogeneous absorption to more closely mimic the overall optical properties of head tissue. Data were acquired with a multichannel frequency-domain commercial NIR spectrometer (ISS Imagent™, Champaign, Illinois) equipped with 128 laser diodes (64 emitting light at 690 nm and 64 at 830 nm) and 24 photomultiplier tubes. Time multiplexing was employed, so that each detector picked up light from 16 different sources at different times within a multiplexing cycle.

A square array of optodes (minimum source–detector $\text{distance}=15\mathrm{mm}$) was arranged on the surface of the milk in the same configuration reported in Fig. 1. Baseline absorption values were estimated using the frequency-domain system in a multidistance configuration.⁴⁷ The estimated absorption and scattering values over a recording epoch of 2 min were ${\mu }_{\mathrm{a}}=0.0095{\mathrm{mm}}^{-1}$ and ${\mu }_{\mathrm{s}}^{\prime }=0.85{\mathrm{mm}}^{-1}$ at 830 nm.

An absorbing sphere with $\text{radius}=3\mathrm{mm}$ and ${\mu }_{\mathrm{a}}=0.02{\mathrm{mm}}^{-1}$, ${\mu }_{\mathrm{s}}^{\prime }=0.95{\mathrm{mm}}^{-1}$ at 830 nm of wavelength was embedded at the center of the investigation surface and moved with a mechanical mechanism to different depths from the milk surface (center varying from 3 to 33 mm from the surface). Data at different depths were acquired for 30 s with a 40-Hz sampling rate. Reference (homogeneous medium) intensities were acquired with the absorber at $\text{depth}=40\mathrm{mm}$. Reference data were taken for 30 s after each depth recording in order to compensate for slow instrumentation drifts. Data acquired using 830-nm wavelength NIR light were used for image reconstruction. Channels that presented a variation of the intensity $>40\%$ for any of the depth recordings when compared to baseline were considered noisy and disregarded during the image reconstruction process. ED and FWHM metrics were estimated for the reconstructed images computed using ELR-DOT, ER-DOT, and SVR-DOT.

2.8.

Real (In-Vivo) Data: Eccentricity Study

In order to quantify the ability of ELR-DOT to localize functional activation in vivo, we compared results obtained with fMRI and ELR-DOT for one subject performing a visual eccentricity task. The participant was recruited and run with approval from the University of Illinois Institutional Review Board. The visual eccentricity task is particularly suitable for evaluating algorithm localization performance. In fact, it is well known that the visual cortex is spatially organized as a function of the eccentricity (and polar angle) of visual stimulation. More importantly for testing the accuracy of depth estimations, the more eccentric the stimulus, the deeper the areas of the primary visual cortex involved in its processing.^48–50

The data reported in this paper are from one participant, who underwent a typical hemodynamic imaging eccentricity study in which a ring comprising a contrast-reversing checkerboard pattern (2 Hz) was presented at the center of the screen and expanded (and in a separate block, contracted) continuously in size, increasing or decreasing the stimulated visual angle.⁵¹ Visual angles from 0 to 10 deg were spanned during the experiment. Within a stimulation period of 48 s, the expanding or contracting cycle was repeated eight times. The stimulus was designed to generate asynchronous responses in regions of the primary visual cortex corresponding to different visual angles, and therefore varying in depth. fMRI data were acquired with a Siemens Trio^® 3T scanner using a BOLD sequence ($\mathrm{TR}=2000\mathrm{ms}$, $\mathrm{TE}=25\mathrm{ms}$, $\text{flip angle}=90\mathrm{deg}$, in-plane $\text{resolution}=2.5\times 2.5\mathrm{mm}$, $\text{slice thickness}=2.5\mathrm{mm}$). fMRI images were coregistered to the structural resampled MRI of the same participant, producing a four-dimensional (4-D) time-course image in structural space.

fNIRS data were acquired with a multichannel frequency-domain NIR spectrometer (ISS Imagent™, Champaign, Illinois). Optical data from the same participant were recorded using an array of eight source pairs (830 and 690 nm) and 16 detectors positioned on the scalp over the occipital cortex [with source–detector distances ranging from 20 to 70 mm; Fig. 11(c)]. The optode array was designed so as to provide coverage of the primary visual cortex, and, most importantly for the purposes of the current paper, to provide a wide variance of source–detector distances over a relatively small surface. As in Ref. 48, this was achieved by grouping all the sources over one hemisphere and all the detectors over the other. Simulations and real-task data indicate that this “grouped” geometry is very effective (compared to a square-grid geometry) in providing a good 3-D reconstruction of perturbations located in between the sources and the detectors. This montage also avoids the problem of having too many short-distance channels, challenging the dynamic range of the detectors. Fiducial markers were placed on the participant’s left and right preauricular points and on the nasion (Na). Fiducials, optodes, and other scalp locations were digitized with a Polhemus FastTrak 3-D digitizer (Colchester, VT; accuracy: 0.8 mm) using a recording stylus and three head-mounted receivers, which allowed for small movements of the head in between measurements. Optode locations and structural MRI data were coregistered using fiducials and a surface-fitting Levenberq–Marquard algorithm.⁴⁵

Only data recorded with the 830-nm wavelength were used for the computation of DOT. Preprocessing involved movement correction,⁵² removal of noisy channels (i.e., channels that presented a variation in intensity greater than 40%), and signal resampling (0.5 Hz). The last step of preprocessing regressed out the shortest channels’ signals from all the good channels in order to dampen the effect of superficial interference.⁵³

After applying ELR-DOT, we obtained an estimate of the absorption changes at 830 nm for each node considered. Nodes with a very small sensitivity were masked (norm equal to ${10}^{-6}$ of the maximum norm estimated). A voxel-based image of absorption changes in the head was obtained by cubic spline interpolation of the absorption values in node space.

Because of the increased oxygenation of the brain in active areas and the hemoglobin extinction coefficients at 830 nm of wavelength, changes in absorption at the wavelength of interest should be proportional to the BOLD response of the fMRI. The same statistical analysis in the same structural space (T1 subject MRI space) was conducted on both fMRI data and fNIRS ELR-DOT data.

The 4-D functional images of the two blocks (expanding and contracting rings) were subtracted from one another in order to eliminate the effect of the lag of the hemodynamic response function. Each voxel in the obtained 4-D images was correlated in time with a sine and a cosine of a period of $T=48\mathrm{s}$. As the stimulus varied continuously in eccentricity as a function of time, for both fMRI and fNIRS data, the phase of the hemodynamic response reflected the tuning of the brain activity to a particular stimulus eccentricity. The phase estimate at the particular oscillation frequency (48 s) was obtained for each voxel as

A $t$-score $t(x,y,z)$ was obtained by correlating the signal in each voxel with a cosine shifted by the amount $\mathrm{Ph}(x,y,z)$, corresponding to each stimulation eccentricity. Knowing the ranges of stimulated visual angles and the time when they occurred, it was possible to linearly relate $\mathrm{Ph}(x,y,z)$ to visual eccentricity, thus obtaining a $\mathrm{VA}(x,y,z)$, where $\mathrm{VA}(x,y,z)$ was a brain map of the visual angle that maximally stimulated a particular voxel.⁵⁴

fMRI results were smoothed using a Gaussian filter with a $\mathrm{FWHM}=20\mathrm{mm}$ in order to account for the reduced spatial resolution of DOT when compared to fMRI. fMRI and ELR-DOT results were compared by correlating statistically significant voxels ($p<0.05$) $\mathrm{VA}(x,y,z)$ over space of the two imaging procedures. Here, we report results obtained for a single subject to show feasibility. This subject had particularly robust effects for both fMRI and DOT analysis. Statistical analyses on a bigger sample of subjects will be reported in a separate paper.

3.1.

Simulations: Slab Geometry

The purpose of the first set of simulations was to estimate the relative weights to be given to the energy and Laplacian regularization matrices ($ϵ$ and $\lambda$) when combining them. In order to better quantify the effects of the regularization parameters $ϵ$ and $\lambda$, the quantity ${\mathrm{J}}_{[{\mu }_{\mathrm{a}}(\overrightarrow{r}),{\mu }_{\mathrm{s}}^{\prime }(\overrightarrow{r})]}{(I+\frac{\lambda }{3}{\mathrm{L}}^{\mathrm{T}}L)}^{-1}{{\mathrm{J}}_{[{\mu }_{\mathrm{a}}(\overrightarrow{r}),{\mu }_{\mathrm{s}}^{\prime }(\overrightarrow{r})]}}^{\mathrm{T}}$ was set to 0.

Figure 3 reports the channel average inversion matrix $G$ for a slice passing through the middle of and perpendicular to the recording surface for different ratios of $ϵ$ to $\lambda$. For a ratio $ϵ/\lambda ={10}^4$, the matrix $G$ is almost identical to the forward solver ${{\mathrm{J}}_{[{\mu }_{\mathrm{a}}(\overrightarrow{r}),{\mu }_{\mathrm{s}}^{\prime }(\overrightarrow{r})]}}^{\mathrm{T}}$. This ratio shows the effect of a strongly energy-weighted regularized solver for which the data are forced toward the surface, the region of higher sensitivity. As the ratio of $ϵ$ to $\lambda$ decreases, the Laplacian regularization increases its effect, clearly moving the preferential solutions to deeper tissues. However, when the ratio is very low (i.e., $={10}^{-4}$), the operator generates physically unlikely (very deep) solutions (e.g., $>5\mathrm{cm}$). Thus, Fig. 3 suggests that an intermediate ratio between the two regularization parameters might provide an optimal solution.

Fig. 3

Channel average inverse operator using cube geometry, presented on a slice perpendicular to the recording surface, passing through the middle of the cube. Different images show a different combination of ratios between $ϵ$ (energy regularization parameter) and $\lambda$ (Laplacian regularization parameter). As the ratio between $ϵ$ and $\lambda$ decreases, the Laplacian regularization increases its effect, clearly moving the preferential solutions to deeper tissues.

The second step of our analysis was to estimate the best $ϵ$ and $\lambda$ for ELR-DOT. For comparison purposes, the best $ϵ$ for ER-DOT and the best $ϵ$ and $\beta$ for SVR-DOT were computed. We estimated the best regularization values running forward simulations with perturbations of baseline optical properties, positioned in the homogeneous media at different depths (see Sec. 2). The goal of the optimization was to find the minimum of the maximum PE for all the depths (we wanted accurate results for all the depths considered) given the parameters $ϵ$, $\lambda$, and $\beta$. The values were estimated normalized to the $\max {(\Vert {\mathrm{J}}_{[{\mu }_{\mathrm{a}}(\overrightarrow{r}),{\mu }_{\mathrm{s}}^{\prime }(\overrightarrow{r})]}\Vert )}^2$:

The analysis was conducted by varying $E$, $\mathrm{\Lambda }$, and $B$ from ${10}^{-7}$ to ${10}^{-1}$. Similarly to previous studies,^12,22,23 we found an optimal energy regularization parameter equal to $E={10}^{-5}$. Thus, in order to focus on the effects of the Laplacian regularization and the SVR, here we present results obtained at a fixed energy regularization equal to $E={10}^{-5}$.

Increasing the value of $\mathrm{\Lambda }$ increases the relative weight given to the Laplacian regularization, whereas decreasing the value of B increases the relative weight given to the spatial regularization. For display purposes, we report a value of spatial regularization $\tilde{B}={10}^{-\log 10(B)-8}$. Increasing $\tilde{B}$ increases the effect of the spatial regularization.

Figure 4(a) reports maximum PEs as a function of Laplacian regularization $\mathrm{\Lambda }$ and SVR $\tilde{B}$. It can be noticed that ELR-DOT provides parameter combinations that result in lower maximum PEs than SVR-DOT, especially with increasing values of $\mathrm{\Lambda }$ and $\tilde{B}$, i.e., because SVR-DOT becomes very unstable for high values of $\tilde{B}$ (low values of $B$). This results in a very noisy reconstructed image even when a low level of noise is added to the simulations, as was the case here. In general, regularization is introduced to provide stability to a generally ill-posed problem, generated by the fact that the matrix ${\mathrm{J}}_{[{\mu }_{\mathrm{a}}(\overrightarrow{r}),{\mu }_{\mathrm{s}}^{\prime }(\overrightarrow{r})]}{{\mathrm{J}}_{[{\mu }_{\mathrm{a}}(\overrightarrow{r}),{\mu }_{\mathrm{s}}^{\prime }(\overrightarrow{r})]}}^{\mathrm{T}}$ is typically almost singular. In order to better address this issue, we tested how the regularization parameters $\mathrm{\Lambda }$ and $\tilde{B}$ affected the stability of the procedure for the given energy regularization. Figure 4(b) reports the computed Frobenius norms of $G$, ${\Vert G\Vert }_F$, for different values of Laplacian regularization $\mathrm{\Lambda }$ and SVR $\tilde{B}$ (smaller values indicate greater stability). Note that the abscissa’s origin can be considered the ER-DOT ${\Vert G\Vert }_F$ with $E={10}^{-5}$.

Fig. 4

(a) Maximum PE as a function of the Laplacian regularization parameter $\mathrm{\Lambda }$ and the SVR parameter $\tilde{B}$ at a fixed energy regularization parameter $E={10}^{-5}$. (b) Frobenius norms of the inverse operator $G$ as a function of the Laplacian regularization parameter $\mathrm{\Lambda }$ and the SVR parameter $\tilde{B}$ at a fixed energy regularization parameter $E={10}^{-5}$. The results were obtained with simulations in the slab geometry configuration.

One important conclusion can be derived from Fig. 4. ELR-DOT stabilizes the problem at a given ER-DOT solution (decreasing the Frobenius norm), whereas SVR-DOT does the opposite (increasing the Frobenius norm). This finding indicates that with respect to the stability criterion, preference should be given to ELR-DOT over SVR-DOT when choosing a method for correcting for shallow-depth bias.

Based on the results depicted in Figs. 4(a) and 4(b), and closely matching values independently reported in the literature for SVR-DOT and ER-DOT,^12,22,23 we chose the optimal regularization parameters as

These parameters were optimized for a coarse mesh with 5-mm internode distance. Therefore, optimal values could differ if a different spatial sampling were used. In fact, since the Laplacian operator is defined over space, its effects will vary as a function of internode distance. However, unless the spatial sampling is drastically altered, we do not expect a substantial (i.e., orders of magnitude) change in optimal regularization parameters. Note also that the Frobenius norm value for the chosen parameter is $\sim 10$ to 15 times lower for ELR-DOT compared to ER-DOT and SVR-DOT, respectively [Fig. 4(b)].

Figure 5 reports the PEs and the FHWMs based on the chosen parameters $E$, $\mathrm{\Lambda }$, and $B$, at different perturbation depths for ELR-DOT, ER-DOT, and SVR-DOT. As can be noticed from this figure, using ER-DOT leads to a low PE in the shallow depth range, with PE increasing linearly for depths $>15\mathrm{mm}$. Using SVR-DOT leads to a very low PE ($\sim 1\mathrm{mm}$) in a range between 11 and 23 mm, but, as for the ER-DOT, the PE increases linearly for depths $>23\mathrm{mm}$ (it also leads to greater errors for depths $<10\mathrm{mm}$). ELR-DOT leads to stable PEs ($\sim 2$ to 3 mm) for a wide range of depths (from 5 to 35 mm, i.e., for the entire depth spectrum explored). The FWHMs for the three methods are shown in Fig. 5(b). At greater depths, FWHM increases, with ELR-DOT plateauing around 25 mm and ER-DOT and SVR-DOT plateauing earlier and at a lower overall FWHM.

Fig. 5

(a) PEs for ELR-DOT, ER-DOT, and SVR-DOT obtained for perturbation depths between 5 and 35 mm. (b) FWHMs for ELR-DOT, ER-DOT, and SVR-DOT obtained for perturbation depths between 5 and 35 mm. The regularization parameters were ELR-DOT: $E={10}^{-5}$, $\mathrm{\Lambda }={10}^{-4}$, ER-DOT: $E={10}^{-5}$, SVR-DOT: $E={10}^{-5}$, $B={10}^{-3}$. The values were obtained during simulations in the slab geometry configuration.

Taken together, Figs. 4 and 5 clearly show how ELR-DOT trades spatial resolution at deep locations for stability and improved depth localization above 20 mm. However, for a gain of $\sim 10$ to 15 times in stability (${\Vert G\Vert }_F$) and a smaller PE across the range of depths (5 to 35 mm), the FWHM of the reconstructed perturbation was increased only by $\sim 30\%$.

In Fig. 6, reconstructed images related to perturbations varying in depths (indicated by white circles) are displayed from left to right and from top to bottom. The minimum depth of the perturbation centroid was 5 mm and the maximum 35 mm, with 2-mm steps. As the figure shows, ELR-DOT follows the perturbation’s true location (white circle) accurately for the whole range of depths considered. Only for the last two slices (33 and 35 mm of centroid depth) do the reconstructed perturbation and the real position seem off center. However, what is off center is the maximum value of the reconstructed image. It is clear that the attenuation of the reconstructed perturbation is not symmetric along depths, being more elongated along the horizontal axis for deeper perturbations.

Fig. 6

A single slice of the reconstructed images obtained during simulations for the cube geometry. The slice is perpendicular to the recording surface and passes through the middle of the cube. Perturbation depth (represented by the white circles) increases (moving from left to right and from top to bottom) from a centroid at 5 mm to one at 35 mm, in 2-mm steps. The color scale represents the change in absorption in arbitrary units.

To assess the resolution (as opposed to localization) of ELR-DOT, we performed a simulation based on the same homogeneous slab geometry, considering two perturbations at different lateral separation distances (PD) and a constant depth of 25 mm. The depth was chosen to maximize the FWHM [Fig. 5(b)]; thus, these results are obtained in a “worst-case scenario” of deep perturbations.

Three reconstructed images over a range of PDs at a constant depth are displayed in Fig. 7(a), with their corresponding profiles reported in Fig. 7(b). Clearly, at a distance of 10 mm, the algorithm fuses the two perturbations, whereas at distances of 24 and 36 mm, the algorithm is able to distinguish between them. The parameter used to assess the ability of the algorithm to separate two close perturbations was RP (see Sec. 2.4). Figure 7(c) reports RP as a function of centroid distance of the two perturbations. When RP is equal to or greater than 1, the algorithm is not able to separate the two perturbations. Not surprisingly, Fig. 7(c) shows that the minimum distance (considered from the centroids) that allows ELR-DOT to separate them is $\sim 22\mathrm{mm}$.

Fig. 7

(a) Reconstructed images (single slice) obtained during simulations with two perturbations at a constant depth (25 mm) in the cube geometry. Perturbations distances (PDs) of 10, 24, and 36 mm are shown. (b) Profiles at 25-mm depth of the reconstructed images for each of the PDs (10, 24, and 36 mm). (c) RP as a function of the PD.

3.2.

Simulations: Realistic Head Geometry

The stability of the regularization parameters of ELR-DOT, ER-DOT, and SVR-DOT to changes in geometry and baseline optical values was tested using a set of heterogeneous head models obtained from anatomical images of real heads. The head’s geometry and optical values were derived from T1-weighted MRIs from five adult subjects. Simulations were conducted by placing optical perturbations at different depths inside the head. A forward solver was computed to estimate the amplitudes of intensity effects using an optical montage positioned over the occipital cortex (Fig. 2). Inverse procedures were performed using tuned regularization parameters derived from the slab geometry simulation.

Figure 8(a) reports examples of reconstructed absorption perturbations using ELR-DOT for a single subject at three depths from the scalp ($\text{blue}=15\mathrm{mm}$, $\text{red}=25\mathrm{mm}$, $\text{green}=33\mathrm{mm}$). The absorption image is thresholded at 50% of its maximum value. The true position of the centroid of the perturbations is displayed by the intersection of the cross. Qualitatively, Fig. 8(a) shows that the ELR-DOT algorithm is able to retrieve the actual position of the perturbation and that the FHWM increases as the depth increases. Figure 8(b) reports the average ED and related standard errors across the five subjects as a function of the actual perturbations centroid depth for ELR-DOT, ER-DOT, and SVR-DOT. No statistically significant discrepancies from the ideal curve were found for ELR-DOT. Significant underestimation of depth reconstruction was found for ER-DOT and SVR-DOT for depths $>20\mathrm{mm}$. Figure 8(c) reports the average FHWM as function of the centroid depth. The FWHM of the reconstructed image increased as the actual depth increased up to a value of $\sim 23\mathrm{mm}$ for ELR-DOT. Similar behavior was found for ER-DOT and SVR-DOT up to a value of $\sim 14$ to 16 mm. These results are in agreement with the slab geometry simulations and show the stability and accuracy of ELR-DOT in a realistic perturbations medium when compared to ER-DOT and SVR-DOT.

Fig. 8

(a) Example of a reconstructed absorption perturbation for a single subject at a depth of 15 (blue), 25 (red), and 33 mm from the scalp (green). The absorption is reported until an attenuation of 50% from its maximum value (FWHM). The true position of the centroid of the perturbations is displayed by the cross-hair intersections. (b) Average ED and related standard errors across the five subjects as a function of the actual perturbation’s centroid depth for ELR-DOT, ER-DOT, and SVR-DOT. (c) Average FHWMs and related standard errors as a function of the actual perturbation’s centroid depth for ELR-DOT, ER-DOT, and SVR-DOT. Results were obtained running simulations on heterogeneous head geometry.

3.3.

Phantom Data: Slab Geometry

A phantom study was performed in order to test the ability of the FEM forward solver and the ELR-DOT algorithm to reconstruct optical images. Figure 9(a) reports the ED as function of the centroid depth (mm) for ELR-DOT, ER-DOT, and SVR-DOT. Consistent with the simulations, ELR_DOT showed a nice fit with the ideal curve (solid line), with an average error of $\sim 2\mathrm{mm}$, and ELR_DOT performed better than ER-DOT and SVR-DOT, especially at greater depths. As expected, based on the simulations, Fig. 9(b) indicates that the FWHM of the reconstructed image is slightly greater for ELR-DOT compared to ER-DOT and SVR-DOT at depths greater than 10 mm.

Fig. 9

(a) EDs as a function of the actual perturbation’s centroid depth for ELR-DOT, ER-DOT, and SVR-DOT. (b) FHWMs as a function of the actual perturbation’s centroid depth for ELR-DOT, ER-DOT, and SVR-DOT. Results were obtained using real (phantom) data in the slab geometry configuration.

Figure 10 displays the reconstructed images resulting from placing the absorbing object at different depths (indicated by a white circle) within the homogenous milk medium. The minimum depth of the perturbations centroid was 3 mm, the maximum 33 mm with 2-mm steps. Although the images were obtained for a $100\times 100\times 100\mathrm{mm}$ cube, the images are shown up to a depth of 50 mm for display purposes. These real images are noisier than the simulated images reported in Fig. 6. However, they show that the ELR-DOT method accurately retrieved the true perturbations’ location (white circle) for all the depths considered (up to 33 mm).

Fig. 10

Slices of the reconstructed image obtained for real (phantom) data collected with the homogenous slab geometry. Each slice is perpendicular to the recording surface and is passing through the middle of the cube. Perturbation depth (represented by the white circles) increases (moving from left to right and from top to bottom) from a centroid at 3 to 33 mm in 2-mm steps. The color scale represents the change in absorption in arbitrary units.

3.4.

Real (In-Vivo) Data: Eccentricity Study

Finally, to test the ability of ELR-DOT to localize functional activation in vivo, we compared the depth localization of brain activity obtained with fMRI and ELR-DOT for a subject undergoing a visual eccentricity task. Figure 11 reports a mapping of the brain response, measured with fMRI [Fig. 11(a)] and ELR-DOT [Fig. 11(b)], as a function of the visual angle of stimulation $\mathrm{VA}(x,y,z)$ on a particular axial slice overlaid onto the subject’s anatomical image. $\mathrm{VA}(x,y,z)$ is displayed only for statistically significant voxels ($p<0.05$, unadjusted for multiple comparisons). Colors represent the visual angle at which the response was maximal, coded according to a color scale presented in Fig. 11(c). The fMRI results [Fig. 11(a)] showed the characteristic spatial organization of the visual cortex as function of the eccentricity of the visual stimulation. This organization consists of a progressively deeper response for more eccentric stimuli. Figure 11(b) reports the same image obtained by applying ELR-DOT to 830-nm optical intensity fNIRS data. The volume of the significantly activated region is smaller for fNIRS when compared to fMRI. This result reflects the lower signal-to-noise ratio, the masking procedure for deep regions, and the limited optode coverage for lateral activations in DOT. The loss of very superficial activations for ELR-DOT when compared to fMRI is caused by the subtraction of short channels (at 20-mm interfiber distance) to compensate for systemic fluctuations. However, qualitatively, the two maps show strong similarities. To quantify them, we first masked the fMRI locations to match the more limited area showing activation with ELR-DOT [see black tracing in Figs. 11(a) and 11(b)]. We then divided the area into four ROIs, the boundaries of which were defined by the visual angle (2.5 deg, 2.5 to 5.0 deg, 5.0 to 7.5 deg, and 7.5 to 10.0 deg), which optimally activated each voxel in the ELR-DOT data (i.e., the ROIs encompassed four receptive field ranges). The centroid (in degrees of angle) of the receptive field was estimated as the average visual angle associated with the maximum response for each voxel in a given ROI (the same operation was computed for ELR-DOT data). Figure 11(d) shows a scatterplot of the two estimates of the receptive field’s centroids for each ROI (estimates based on ELR-DOT on the abscissa and on fMRI on the ordinate). The data show a very good correlation between the centroids of the receptive fields obtained with the two techniques ($r=0.95$, $p=0.04$). Across the four ROIs, the largest mismatch between fMRI- and ELR-DOT–based receptive field centroids was $\sim 2\mathrm{deg}$. This is translated into a maximum mismatch of $\sim 4\mathrm{mm}$ in depth estimation between the two methods. Importantly, the maximum depth of the ELR-DOT activated area from the surface of the scalp was 33 mm. These results support the claim that (1) ELR-DOT provides depth estimates that closely correspond with those obtained with fMRI, the gold-standard for 3-D functional imaging in vivo, and (2) ELR-DOT can accurately localize brain activity up to a depth of 30 to 33 mm from the surface of the scalp.

Fig. 11

(a) fMRI: colors represent activation correlated with different visual angles overlaid on an axial anatomical slice through primary visual cortex (significant voxels, $p<0.05$, uncorrected). Black contours define statistically activated regions for both fMRI and ELR-DOT. (b) ELR-DOT: colors represent activation correlated with different visual angles overlaid on an axial anatomical slice through primary visual cortex (significant voxels, $p<0.05$, uncorrected). Black contours define statistically activated regions for both fMRI and ELR-DOT. (c) Actual digitized montage used for the optical recording (rear view of the head): the sources are indicated in blue and the detectors in red—all combinations of sources and detectors were used. (d) Scatterplot of the average visual angle in four ROIs for ELR-DOT and the average visual angle for the corresponding fMRI regions ($r=0.95$, $p=0.04$). The color coding for the different visual angles is presented in bar graph (c).