2.1.

Forward Problem

To model the light propagation in tissue, commonly known as the forward model, a numerical model based on the finite element method (FEM) has been employed to solve the diffusion approximation through the NIRFAST³² modeling and image reconstruction software package

2.2.

Inverse Problem

The sensitivity matrix, $\mathbf{J}$, also known as the Jacobian or weight matrix is calculated using the Adjoint method³⁴ that relates changes in boundary fluence measurements, $\partial \mathbf{y}$, with respect to small changes in underlying tissue optical parameters, $\partial \mathbf{x}$:

From Eq. (1), in the frequency domain, the measured data are complex numbers and can be decomposed in two parts: amplitude ${\mathbf{y}}_{\mathrm{I}}$, and phase ${\mathbf{y}}_{\mathbf{\Theta }}$. Using the Rytov approximation, differential measurements with respect to a baseline ${\mathbf{y}}_0$ are described as a logarithmic ratio for amplitude and a difference in phase:

Therefore, the frequency-domain Jacobian describes changes in logarithmic intensity, $\partial \ln {\mathbf{y}}_{\mathrm{I}}$, and linear differences in phase $\partial {\mathbf{y}}_{\mathit{\theta }}$, with respect to changes in the absorption, $\partial {\mu }_a$, and the diffusion, $\partial \mathbf{\kappa }$ coefficients:

Assuming that only changes in the absorption coefficient are associated with vascular hemodynamic changes within the brain, this work will focus on modeling and reconstructing changes in absorption, $\partial {\mu }_a$, alone, reducing Eq. (4) to

This equation is used in conjunction with a noise model to produce realistic boundary data by perturbating hemoglobin concentration values.

To demonstrate the difference in the sensitivity distribution and magnitude of both log-amplitude and phase, as a function of source and detector distance, the Jacobian for a three-layered model, at a single wavelength of 830 nm, at a modulation frequency of 140 MHz is shown in Fig. 1. Here, the FD Jacobian has been calculated for four source and detector separations of 10, 20, 30, and 40 mm, for a three-layered two-dimensional model, with a skin/scalp layer of 5 mm, skull layer of 8 mm, and brain layer (combining gray and white matter) with optical properties shown in Table 1.

Fig. 1

Frequency-domain Jacobian for Log intensity ($\text{Joules}/{\mathrm{mm}}^2/\mathrm{mM}$) and phase (radians/mM) for ${\mathrm{HbO}}_2$. Each row corresponds to a different source and detector distances of 10 to 40 mm. The dashed line represents the boundary of each layer, skin/scalp (top), bone (middle), and brain (bottom) tissue, each having different optical properties, as depicted in Table 1.

Table 1

Baseline optical properties in the head model.35

Several features are evident from Fig. 1: (1) the Jacobian for the phase is showing sensitivity weighting at deeper regions as compared to log intensity, highlighting the potential of deeper tissue sampling, (2) the magnitude of change for log intensity measurements is larger than that compared to the phase but, (3) the distribution of the sensitivity for the phase appears more homogeneous across all tissues as compared to log intensity reducing the “hypersensitivity” to superficial signal contamination. These same features can also be demonstrated for HbR and is also true for a wavelength of 690 nm.

Functional DOT (fDOT) is concerned with three-dimensional spatial mapping of dynamic hemodynamic changes that occur during brain activations, using temporal data as measured from the scalp. Assuming that these vascular dynamic changes are small, the inverse model as used for parameter recovery can be considered as a linear problem. A Moore–Penrose pseudoinverse with Tikhonov regularization can be used to find an approximation of the inverse of the Jacobian in Eq. (2), ${\mathbf{J}}^{\#}$, to perform a single step linear recovery of absorption:

For the FD Jacobian, Eq. (5), the regularization is typically applied separately for each of the kernels of log intensity and phase measurements.³² In this work, the absorption coefficient is recovered for each wavelength then spectrally unmixed to recover ${\mathrm{HbO}}_2$ and HbR as described in Eq. (6). In previous work, spatially constrained reconstructions³⁷ or spatially variant regularization schemes^20,38 have been proposed that facilitate reconstructions of perturbations at deeper locations. To allow a direct evaluation of the potential increase in imaging performance afforded by the additional phase-based measurements, no spatially variant regularization methods have been considered in this work.

2.3.

Subject-Specific Models

The subject-specific models are based on 14 female and 10 male individual subjects with a mean age of 26 ($\pm 4$), with T1-weighted MPRAGE [echo time = 3.13 ms, repetition time (TR) = 2400 ms, $\text{flip angle}=8\mathrm{deg}$, $1\times 1\times 1{\mathrm{mm}}^3$ isotropic voxels] scans acquired for each subject. All subjects passed MR screening to ensure their safe participation. Informed consent was obtained and the research was approved by the Human Research Protection Office at Washington University School of Medicine.

High-resolution subject-specific FEM, consisting of five tissue layers were created, with mean element volume of $\sim 1{\mathrm{mm}}^3$, resulting in $\sim \mathrm{450,000}$ nodes for each mesh using NIRFAST.³⁹ A sparse-sensitivity matrix, $J$, (sensitivity threshold set to 0.1%⁴⁰), was then created for each of the 24 subject-specific models with realistic optical properties for a 690- and 830-nm wavelengths, shown in Table 1. The reduced scattering coefficient was derived from the scatter amplitude and power, using ${\mu }_s^{\prime }=s_a{\lambda }^{-s_p}$, where $s_a$ is the scatter amplitude, $s_p$ is the scatter power, and $\lambda$ is the wavelength in micrometers.

In order to generate the boundary data, the nodes on the surface of the brain for each subject laying directly under the imaging array were selected to simulate point perturbations for ${\mathrm{HbO}}_2$ and HbR, Fig. 2(a). The data were then used for the calculation of a PSF to allow an analysis of the quality of the image reconstruction within a reasonable field of view. For each subject, the selected region depth from the surface of the brain varies, as the brain surface changes along the sulci and gyri morphology. To better demonstrate the depth variation of the surface of the brain, the histograms of brain surface depth across all 24 subjects are shown in Fig. 2(b) whereas Fig. 2(c) shows the brain surface depth probability distribution for each subject-specific model. This further highlights the importance and the need for improving the depth resolution of HD-DOT and demonstrates the need for modeling individual surface activation points for the presented analysis. The brain surface depth varies across the subjects, Fig. 2(c), ranging from 5 to 30 mm and brain surface points presented in this study lie on the cortical surface with a median depth of 12.5 mm, Figs. 2(b) and 2(c).

Fig. 2

(a) The modeled HD-DOT system with 158 sources (red) and 166 detectors (cyan), the surface of the brain surface is shown in pink and blue. The pink region denotes cortex locations of point perturbations in the simulations. (b) Histogram of brain surface depth across 24 subjects. (c) Brain surface depth probability distribution for each subject-specific model. The red stars mark the maximum probability for each subject while the red line is the most probable depth across all subjects (12.5 mm).

2.4.

Realistic System Noise Model

The modeled system array adopted in this work is an HD-DOT system, with SD separation distances ranging from 13 to 47 mm as described previously.⁵ The SD pairs are extended over the whole of the head, in order to allow coverage of a large area of the optically accessible cortical surface of the brain that has previously been used as a benchmark for fDOT simulation studies.⁴¹ The modeled HD-DOT array contains a total of 158 light source at two wavelengths of 690 and 830 nm, and 166 detectors, resulting in 3,500 associated SD pairs within 50 mm separation at each wavelength, Fig. 2(a). The SD pairs are grouped into sets of nearest neighbor (NN) measurements (i.e., NN1, NN2, NN3, NN4, for first through fourth NN, respectively) with mean distances of 13, 29, 39, and 47 mm. In this work, the labeling for reconstruction with respect to a set of NN measurements used is the maximum NN; for example, for NN4 includes NN1, NN2, NN3, and NN4 measurements.

To generate a realistic noise model to be used with simulations, data were recorded using a frequency-domain system console (ISS Imagent™, Champaign, Illinois) employing six detectors at distances of 18 to 58 mm from the source, placed over the optical cortex. The source was modulated at 140 MHz and data were sampled at 39.74 Hz. Data were recorded for 120 s while the subject sat still and quietly fixated on a blank screen. Data were then processed following methods that have been previously described.¹ Briefly, the data were band-pass filtered to 0.01 to 0.1 Hz to remove slow drifts, systemic physiology (heartbeat and respiration), downsampled to 1 Hz, and block averaged in 10 blocks of 12 s each. The same processing was applied to the measured amplitude and phase signals. The noise was then estimated as the standard deviation of the log mean intensity signals; therefore, it is a ratio-metric difference for the intensity and a linear phase difference (in degrees) for the phase and is shown in Fig. 3. As the system utilized photomultiplier tubes for signal detection, each detector was at a different voltage bias to increase signal gain and throughout the experiment, all data were ensured to be above the expected noise floor.

Fig. 3

In-vivo measured noise and a two-term exponential fit of noise versus distance from the source for intensity (left) and phase (right) measurements, on 690 nm (red diamonds) and 830 nm (pink squares) wavelengths.

The noise model $n$ was calculated as a function of SD distance, $r$, via a two-term exponential of the form $n(r)=a{e}^{br}+c{e}^{dr}$ and was found to provide the best fit with coefficients shown in Table 2. In a previous work, using the same HD array configuration and similar signal processing but acquiring measurements with a CW system with LEDs at 830 nm, the amplitude noise was reported to be 0.12% for NN1, 0.15% for NN2, 0.41% for NN3, and 1.42% for NN4.¹² The empirically derived noise model generated in this work corresponds to a noise level estimate for each neighborhood as shown in Table 3.

Table 2

The estimated coefficients for the two-term exponential fit of the noise model.

Table 3

The estimated noise levels for each measurement neighborhood.

To simulate realistic measured data, perturbations of ${\mathrm{HbO}}_2$ and HbR values were induced using Eq. (7), using perturbation values scaled to minimize the difference between simulated and measured data. The perturbation for ${\mathrm{HbO}}_2$ was scaled to be $3.8\mu \mathrm{M}$ while for HbR was $-1.8\mu \mathrm{M}$. Evaluations were carried out using noise-free simulations and simulations that included the empirically derived noise model. Stochastic realistic noise was linearly added to each measurement based on sampling from a Gaussian distribution with zero mean and a standard deviation that was a function of SD distance with parameters given in Table 3.

2.5.

Evaluation Metrics

To evaluate the image quality of each PSF, for each modeled focal activation, metrics as previously described in similar work have been utilized:¹⁰ (1) full width at half maximum (FWHM) defined as the maximum separation between all pairs of points in the activation above half of the maximum recovery; (2) localization error defined as the distance between the centroid of the recovery and the known location of the perturbed node; and (3) the effective resolution defined as the diameter of a sphere centered at the known perturbed node that can enclose all recovered nodes.

It is likely that the recovered volume above half of the maximum will not create a single focal recovery, resulting in multiple scattered activations.⁴² This is more evident in the cases with added noise, as the reconstruction often has multiple recovered activations. Therefore, to automate the quantification process, such activations are isolated by a connected component algorithm, using the FEM elements as the connectivity index and considering only the nodes that are associated with each other through at least one connected element. The recovered activations are then selected based on the maximum accumulated recovery value for evaluation without the use of any prior regional or locational knowledge. Therefore, a recovery near the perturbation will not necessarily be selected for evaluation, unless it is the strongest integrated recovery in the volume. Any recovery with localization error greater than 8 mm is considered a product of noise and therefore not included in the statistical analysis. The Wilcoxon signed rank test was used to assess the statistical difference between metrics of CW and FD reconstructions.

2.6.

In-Vivo Experiment

An FD NIRS device (IMAGENT™, ISS Inc., Illinois) was used to obtain the NIRS-based data on a healthy subject. The system consists of 32 sources, modulated at 140 MHz, and 30 detector fibers with each source coupled to laser diodes emitting at 690 nm and 830 nm with the detectors being compact and fast photomultiplier tubes (PMT). Subjects were recruited from the University of Birmingham community, and written informed consent was obtained. Subjects were seated facing an adjustable screen while the imaging pad was attached over the occipital cortex with hook and loop strapping. Visual stimuli consisted of rotating logarithmic checkerboards, reversing intensity at 10 Hz, on a 50% intensity gray background. Each rotation completed in 36 seconds and 10 rotations occurred in the run, resulting in 7 minutes total scanning session time, including 30 second baseline periods at the beginning and the end. The raw signals were processed according to previously proposed pipeline,¹ with all of the analysis being performed using NeuroDot, an extendable and user-friendly environment for analysis of optical data from raw light measurements ( https://github.com/WUSTL-ORL/NeuroDOT_Beta).

Scans were in line with a previously reported CW high-density DOT imaging system setup.⁴³ Briefly, the imaging cap consisted of 24 source positions (with two NIR wavelengths—690 nm and 830 nm—at each position) and 28 detectors interleaved in a high-density array. First- (13 mm, 1NN), through fourth NN (48 mm, 4NN) optode pairs were utilized, giving a total of 348 total possible measurements, at a frame rate of $\sim 40\mathrm{Hz}$. Prior to data analysis, measurements with high noise (standard deviation of the logmean of the data larger than 7.5%¹) were disregarded, resulting in 300 measurements passing the quality control for the scan. For the CW data, this was based on FD amplitude measurements only, whereas FD data contained both amplitude and phase. Image recovery of ${\mathrm{HbO}}_2$ and HbR were performed as defined above, using NIRFAST for the calculation of the Jacobian matrix for both CW and FD data.

3.1.

Simulations Results

In this work, the brain surface under the imaging array was selected in 24 subject specific models and point activations modeled as perturbations were induced for each node in the selected area, allowing the recovery of PSFs both with and without noise. To provide an overview of the spatial distribution characteristics of the recovered PSFs, consider three points on the primary motor cortex within Brodmann Area 4, along the precentral gyrus, as shown in Fig. 4(a). These three focal activations are modeled at depths of 13.3 mm, 17.2 mm, and 20.5 mm, Fig. 4(b). PSF of each individual modeled focal activations for both ${\mathrm{HbO}}_2$ and HbR have been reconstructed and as both demonstrated similar results, for conciseness, only those of ${\mathrm{HbO}}_2$ are shown as in Fig. 5.

Fig. 4

(a) Axial view of the brain, the red line denotes the coronal slice selected to view the primary motor cortex, Brodmann Area 4, along the precentral gyrus. (b) The five layers of the model: skin, skull, cerebral fluids, white matter, and gray matter. Selected cortex points in different depths: 13.3 mm (left/red), 17.2 mm (middle/green), and 20.5 mm (bottom/blue).

Fig. 5

Recovered PSFs with different NN measurements (NN2, NN3, NN4), for activations at different depths. The perturbation location is denoted with the green dot. The FWHM contour (when present) is noted in cyan. Recoveries with (a)–(c) CW without noise, (d)–(f) CW with noise, (g)–(i) FD without noise, and (j)–(l) FD with noise).

For the focal point at 13.3 mm depth, Fig. 5, top row, both CW and FD provide qualitatively similar reconstructions for all SD pairs, but the FD demonstrated a visually better reconstruction in each case particularly for shorter NN pairs (NN1 and NN2) measurements. Specifically, for the noise-added CW recoveries, the localization error is 6.5 mm for NN2, 2.3 mm for NN3, and 1 mm for NN4, whereas for the FD recoveries are 2.2 mm, 1.6 mm, and 0.85 mm respectively. Similarly, the effective resolution for CW is 22.9 mm for NN2, 21.4 mm for NN3, and 11 mm for NN4, whereas for the FD recoveries is 17.7 mm, 16.2 mm and 9.8 mm respectively, demonstrating that for larger NN, the image recovery of the focal activation is improving for both CW and FD, and FD is consistently outperforming CW.

For a perturbation at 17.2 mm depth, Fig. 5, middle row, the noise added CW data fails to reconstruct for NN2 and NN3, whereas the NN4 has a 5.2 mm localization error and 17.9 mm effective resolution. The reconstructions based on FD data are qualitatively more accurate and have a localization error of 6.2 mm for NN2, 3.8 mm for NN3, and 2.2 mm for NN4; and the effective resolution is 21.42 mm, 20.58 mm, and 17.9 mm respectively. This case demonstrates that for some depth range the CW reconstructions may fail if the dynamic range of the detectors cannot provide measurements with a good signal to noise ratio (SNR), and the FD recoveries can perform reasonably well even when using only NN2.

When reconstructing for simulated data without noise, an activation positioned 20.5 mm deep from the surface of the head may be recovered reasonably well using any number of measurement neighborhoods, using CW only or FD data (Fig. 5) though the FD based data provide reconstructions with better image quality for all NN cases. In the presence of noise, recovery based on the NN2 and NN3 completely fail in both CW and FD. The recovery for CW using NN3 appears to provide a recovery near the perturbated point; however, this is not the strongest recovery in the volume, and therefore is considered as noise artifact and as such is not selected for evaluation. However, when using NN4, the FD data recovery demonstrates a better recovery, achieving 6.1 mm localization error with 24.8 mm effective resolution. In the presence of noise, recoveries for such deep points often fail, which leads to large increases in the localization and resolution metrics.

To provide a comprehensive evaluation of the benefits of FD data over CW, images were reconstructed for a set of NN1-NN4 SD pairs, using both noise-free and noise added data. In the presence of noise, the recoveries for deep activation points fail, which leads to significant increases in the calculated error metrics. The most obvious change is seen in the localization error, as beyond certain depths the recovered focal activations are not successful. In accordance to existing literature,⁴² any recovery with localization error greater than 8 mm is considered a product of noise. To identify the percentage of recoveries within the imposed localization limit, a success rate metric is introduced (Fig. 6), where a 50% success rate limit was set to examine the depths at which different approaches perform adequately. Recoveries exceeding this limit are not included in any further results or statistical analysis.

Fig. 6

${\mathrm{HbO}}_2$ success rate for 8 mm localization error limit. The dashed blue and red lines denote the depth where the recoveries reach 50% success rate for CW and FD, respectively.

Each recovered focal activation is assessed using metrics of the FWHM (Fig. 7), localization error (Fig. 8), and effective resolution (Fig. 9). These results, quantized according to depth from the surface of the head, are summarized in Tables 4 and 5 and compared in Table 6. The noise-free evaluation is also presented to provide a benchmark of limitations imposed due to the ill-posed inverse problem, the employed regularization, and to provide a lower-bound on metric quality given further advancements in instrumentation that will lead to hardware with lower noise characteristics.

Fig. 7

${\mathrm{HbO}}_2$ localization error median for CW (blue asterisk line) and FD (red cross line); the colored bands represent the 25th to 75th percentile range. The percentile bands stop at the point where the 50% success rate is reached, while the median line continues until the 10% success rate is reached.

Fig. 8

${\mathrm{HbO}}_2$ FWHM median for CW (blue asterisk line) and FD (red cross line); the colored bands represent the 25th to 75th percentile range. The percentile bands stop at the point where the 50% success rate is reached, while the median line continues until the 10% success rate is reached.

Fig. 9

${\mathrm{HbO}}_2$ effective resolution median for CW (blue asterisk line) and FD (red cross line); the colored bands represent the 25th to 75th percentile range. The percentile bands stop at the point where the 50% success rate is reached, while the median line continues until the 10% success rate is reached.

Table 4

Image quality metrics medians and quartiles for recoveries with CW noise-added model.

Table 5

Image quality metrics medians and quartiles for recoveries with FD noise-added model.

Table 6

FD reconstruction relative difference. The bolded cells are not statistically significant (p>0.01).

For the noise-free cases, the success rate maintains near 100% until it steadily declines at a characteristic distance (Fig. 6). The FD approach consistently performs better than the CW case in deeper regions. For CW recoveries, the 50% success rate is crossed at 17.18 mm for NN2, 22.27 mm for NN3, and 24.97 mm for NN4; and the FD 50% success rate points are at 23.04, 27.53, and 29.38 mm, respectively. For the noise-added cases, the success rate drops under 50% at shallower depths. Specifically, for the CW case, the 50% success rate is crossed at 14.87 mm for NN2, 16.59 mm for NN3, and 19.71 mm for NN4, and the FD 50% success rate points are at 18.43, 19.75, and 23.04 mm, respectively.

The FD data significantly improve the localization throughout the imaging domain for noiseless simulations (Fig. 7). At each depth, the use of FD data demonstrates a statistically significant improvement ($p<0.01$) in localization accuracy. For the noise-added case, the localization error is larger for activations deep within the brain, and the improvement of the FD reconstructions is again statistically significant ($p<0.01$) at each depth interval. Specifically, when comparing a depth range of 10 to 15 mm, the localization error for FD presents a 58.9% improvement for NN2 recoveries, 47.3% for NN3 recoveries, and 21.1% for NN4 (Table 6). These improvements are statistically significant for all modeled focal activations ($p<0.01$) using the Wilcoxon-signed rank test.

The FWHM for the noise-free cases increases with depth, but the FD reconstructions do not always significantly improve the FWHM (e.g., in the NN2 case); however, employing the FD approach results in slightly improved performance and smaller variations within each depth interval (Fig. 8). For the noise-added case, the FWHM still presents slight improvements when employing FD recoveries, however, the improvement is not always statistically significant (Table 6), particularly for NN2 at activation depths of greater than 15 mm.

Finally, the effective resolution with respect to depth from the head surface is shown in Fig. 9. As expected, in the noise-free cases, the effective resolution increases with depth. As the effective resolution is a function of both FWHM and localization, in the noise-added cases, the improvements observed are larger than those of the FWHM but smaller than the localization improvements. However, the observed improvement for FD data is statistically significant for each depth interval. Specifically, when comparing for a depth range of 10 to 15 mm, the effective resolution for FD presents a 20.6% improvement for NN2 recoveries, 13.7% for NN3 recoveries, and 10.5% for NN4 (Table 6). These improvements are all statistically significant for all modeled focal activations ($p<0.01$) using the Wilcoxon-signed rank test.

The results are summarized in Tables 4 and 5 where the reconstruction metrics for recoveries with noise-added data have been binned and analyzed in 5-mm depth intervals and the median, upper, and lower quartile for each interval is reported. Finally, in Table 6, the percentage relative change for reconstruction using FD data as compared to CW is shown.

3.2.

In-Vivo Results

A single snapshot of a recovered focal activation for a single subject is shown in Fig. 10. Here the cross-sectional sagittal and transverse planes of a single focal activation within the left visual cortex are shown as a change in recovered ${\mathrm{HbO}}_2$. The CW reconstruction is shown to recover an activation at a much more superficial depth, mainly appearing in the area of the superficial tissue; whereas, the FD-based reconstruction is able to recover the activation overlapping with the cortex of the brain. It is important to highlight that the same regularization parameter for the inverse model has been used, namely ensuring that the improvement in image recovery and accuracy is primarily due to the inclusion of phase data within the FD-based reconstruction.

Fig. 10

In-vivo visual cortex activation (50% threshold, $\mathrm{\Delta }{\mathrm{HbO}}_2$ [a.u.]) using (a, b) CW data and (c, d) FD data.