Measurements of dynamic near-infrared (NIR) light attenuation across the human head together with model-based image reconstruction algorithms allow the recovery of three-dimensional spatial brain activation maps. Previous studies using high-density diffuse optical tomography (HD-DOT) systems have reported improved image quality over sparse arrays. These HD-DOT systems incorporated multidistance overlapping continuous wave measurements that only recover differential intensity attenuation. We investigate the potential improvement in reconstructed image quality due to the additional incorporation of phase shift measurements, which reflect the time-of-flight of the measured NIR light, within the tomographic reconstruction from high-density measurements. To evaluate image reconstruction with and without the additional phase information, we simulated point spread functions across a whole-scalp field of view in 24 subject-specific anatomical models using an experimentally derived noise model. The addition of phase information improves the image quality by reducing localization error by up to 59% and effective resolution by up to 21% as compared to using the intensity attenuation measurements alone. Furthermore, we demonstrate that the phase data enable images to be resolved at deeper brain regions where intensity data fail, which is further supported by utilizing experimental data from a single subject measurement during a retinotopic experiment.

## 1.

## Introduction

Functional neuroimaging is a valuable medical tool that is employed in a wide area of medical applications, with several noninvasive or minimally invasive neuromonitoring techniques for examining functional brain activity available for clinical use. Functional near-infrared spectroscopy (fNIRS)-based optical imaging is a noninvasive, relatively inexpensive technology applied in numerous applications where neuroimaging is crucial, such as functional brain mapping,^{1} psychology studies,^{2}^{,}^{3} intensive care unit patient monitoring,^{4} mental disease monitoring,^{5} and early dementia diagnosis.^{6} Recent advances in system design have also paved the way to promise wireless wearable systems, allowing neuroimaging to be performed in real-life contexts where it is otherwise not possible.^{7}

fNIRS techniques rely on injecting light into the tissue, then measuring the emerging light at nearby locations. Light propagation in biological tissue is dominated by scattering events in the near-infrared (NIR) window, where low absorption allows measurements on the surface of the head to sample the brain cortex.^{8} Differential measurements at two NIR wavelengths, known as NIR spectroscopy (NIRS), can recover oxy (${\mathrm{HbO}}_{2}$) and deoxy (HbR) hemoglobin concentration changes that convey information about functional brain vascular events. When more wavelengths are used, it is possible to additionally recover concentration levels of water, lipids, and even measures of metabolism such as cytochrome c.^{9}

Due to the strong scattering of NIR in tissue, utilizing sparse optical measurements for neuroimaging leads to low spatial specificity as compared to fMRI. Imaging quality and quantitative accuracy can be significantly improved by arranging the NIR optical source and detector elements in a dense grid to provide overlapping measurements that allow tomographic reconstructions, a procedure known as diffuse optical tomography (DOT). Using a high-density (HD) DOT configuration that provides access to overlapping measurements at multiple source–detector (SD) separation distances has been shown to improve image quality^{10} to a point comparable to fMRI.^{1}^{,}^{11} In addition, instruments with a large dynamic range can utilize distant measurement neighborhoods of the HD array, which improves reconstruction quality and depth sensitivity,^{12} as sampling depth increases with SD distance.^{13}

When continuous wave (CW) NIR light is used, measurements provide access to only light intensity attenuation from the tissue. By contrast, amplitude-modulated NIR light, typically in the range of 100 to 200 MHz and known as frequency-domain (FD) NIRS, enables additional measurements of a phase shift of the NIR light.^{14} Measurements of the phase shift are related to the average photon path length and can be used to directly estimate rather than assume scattering parameters within the tissue, which enables more accurate optical parameter recovery.^{15} FD-NIRS has been used in breast optical tomography for characterization and monitoring of lesions,^{16} joint imaging for arthritis detection and treatment monitoring,^{17} infant brain development monitoring,^{18} and brain trauma assessment.^{19} In previous studies utilizing FD NIRS, the phase measurements have typically been disregarded when recovering hemoglobin concentrations.^{20}^{–}^{22} This is largely due to the fact that the accurate absolute measurement of the phase is not trivial and requires complicated and often time-consuming instrument calibration procedures, as the phase can be affected by a number of factors such as fiber length, fiber bending, and attachment angle, each of which can present long-term drifts and requires recalibration. The frequency-domain multidistance approach has been shown to provide increased sensitivity to cerebral oxygenation using intensity and phase data;^{23}^{,}^{24} however, the multidistance approach has been extremely challenging to adopt with an HD-DOT tomographic setup due to the strong demands on the dynamic range of the high bandwidth detector instrumentation.

Alternatively, measurements of picosecond modulated light (time-resolved or TR) provide a direct measure of the temporal point spread function (PSF) of the detected photons to account for both tissue absorption and scattering. However, systems with enough channels to provide overlapping measurements, thus achieving DOT, have a low temporal resolution (75 s for a full imaging cycle^{14}), though various mathematical models have been proposed to increase the temporal resolution.^{25}^{,}^{26} Although TR systems are considered the golden standard within fNIRS, typically only spectroscopic measures are used to provide bulk tissue optical measurements, and these systems have not been widely adopted for HD-DOT.

The CW, FD, and TR measurements are all sensitive to hemodynamics and are therefore indirect measures of variations in brain function. An alternative method, known as the event-related optical signal or fast optical signal, uses just the phase measurements of FD systems and has been reported to be sensitive to neuronal activity.^{27} It has been suggested that this technique does not suffer from the hemodynamic delay effect and may be sensitive to scatter-related changes within firing neurons. However, in practice, this technique requires extensive stimulus repetitions to allow adequate signal averaging to obtain required signal-to-noise, thus canceling the real-time monitoring potential and resulting in experimental sessions of extremely long duration ($\sim 60\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{min}$ for one stimulus^{28}) that are impractical for many applications.

To investigate hemodynamic changes in a practical setting, measurements of differential phase changes are much easier to apply and have shown high correlation with the fMRI blood oxygen level dependent signal, suggesting that phase changes convey information about functional hemodynamic phenomena.^{29} This is an expected result as changes in the absorption of a medium will not only affect the total number of detected photons (the measured intensity) but also the average photon diffusion path-length (the measured phase). How a change in absorption within the volume leads to differential changes in both (logarithmic) intensity and (linear) phase can be understood when considering the distribution of the time of flight of the photons. Early arriving photons travel a shorter distance and stochastically sample more superficial tissues while the later arriving photons travel longer distances and therefore tend to both stochastically sample deeper tissues and have a higher probability of undergoing absorption events.^{30} Therefore, the phase measurement-sensitivity profile to changes in absorption is deeper than the sensitivity of the intensity attenuation measurement because early and late photons contribute equally to the mean time of flight while the total intensity is dominated by early arriving photons.^{31}

To date, FD measurements have not been used in the context of tomographic optical functional brain imaging and the effects of incorporating phase information together with intensity in HD-DOT has yet to be investigated. This work provides a comparison of imaging resolution and localization accuracy between CW and FD HD-DOT by analysis of simulated PSF over 24 subject-specific models including an experimentally measured noise model. The underlying theory is provided to demonstrate the origins of FD measurements. The benefits of FD HD-DOT as compared to CW HD-DOT are also demonstrated through an example of a focal activation map from a retinotopic experiment from a healthy subject.

## 2.

## Methods

## 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^{32} modeling and image reconstruction software package

## Eq. (1)

$$-\nabla \kappa (r)\nabla \mathrm{\Phi}(r,\omega )+[{\mu}_{a}(r)+\frac{i\omega}{{c}_{m}(r)}]\mathrm{\Phi}(r,\omega )=q(r,\omega ),$$^{25}To overcome the computational cost, GPU parallelized linear solvers have been deployed that facilitate the efficient calculation of the forward problem in FD simulations; when using high-resolution meshes with $\sim \mathrm{45,000}$ nodes, with 158 sources and 160 detectors, the forward model calculation time is less than 3 min.

^{33}

## 2.2.

### Inverse Problem

The sensitivity matrix, $\mathbf{J}$, also known as the Jacobian or weight matrix is calculated using the Adjoint method^{34} 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:

## Eq. (3)

$$\mathrm{ln}\left(\frac{y}{{y}_{0}}\right)=\mathrm{ln}\left(\frac{{y}_{\mathrm{I}}{e}^{i{y}_{\mathrm{\Theta}}}}{{y}_{\mathrm{I}0}{e}^{i{y}_{\mathrm{\Theta}0}}}\right)=\mathrm{ln}\left(\frac{{y}_{\mathrm{I}}}{{y}_{\mathrm{I}0}}\right)+i({y}_{\mathrm{\Theta}}-{y}_{\mathrm{\Theta}0}).$$Therefore, the frequency-domain Jacobian describes changes in logarithmic intensity, $\partial \mathrm{ln}\text{\hspace{0.17em}}{\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:

## Eq. (4)

$$\left[\begin{array}{c}\partial \text{\hspace{0.17em}}\mathrm{ln}\text{\hspace{0.17em}}{\mathbf{y}}_{\mathrm{I}}\\ \partial {\mathbf{y}}_{\mathrm{\Theta}}\end{array}\right]=\left[\begin{array}{cc}\frac{\partial \text{\hspace{0.17em}}\mathrm{ln}\text{\hspace{0.17em}}{\mathbf{y}}_{\mathrm{I}}}{\partial {\mu}_{a}}& \frac{\partial \text{\hspace{0.17em}}\mathrm{ln}\text{\hspace{0.17em}}{\mathbf{y}}_{\mathrm{I}}}{\partial \kappa}\\ \frac{\partial {\mathbf{y}}_{\mathrm{\Theta}}}{\partial {\mu}_{a}}& \frac{\partial {\mathbf{y}}_{\mathrm{\Theta}}}{\partial \kappa}\end{array}\right]\left[\begin{array}{c}\partial {\mu}_{a}\\ \partial \kappa \end{array}\right].$$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

## Eq. (5)

$$\left[\begin{array}{c}\partial \text{\hspace{0.17em}}\mathrm{ln}\text{\hspace{0.17em}}{\mathbf{y}}_{\mathrm{I}}\\ \partial {\mathbf{y}}_{\mathrm{\Theta}}\end{array}\right]=\left[\begin{array}{c}\frac{\partial \text{\hspace{0.17em}}\mathrm{ln}\text{\hspace{0.17em}}{\mathbf{y}}_{\mathrm{I}}}{\partial {\mu}_{a}}\\ \frac{\partial {\mathbf{y}}_{\mathrm{\Theta}}}{\partial {\mu}_{a}}\end{array}\right][\partial {\mu}_{a}],$$## Eq. (6)

$$\left[\begin{array}{c}\partial {\mathrm{HbO}}_{2}\\ \partial \mathrm{HbR}\end{array}\right]=\left[\begin{array}{cc}{\epsilon}_{{\mathrm{HbO}}_{2},690}& {\epsilon}_{\mathrm{HbR},690}\\ {\epsilon}_{{\mathrm{HbO}}_{2},830}& {\epsilon}_{\mathrm{HbR},830}\end{array}\right]\left[\begin{array}{c}\partial {\mu}_{a,690}\\ \partial {\mu}_{a,830}\end{array}\right],$$## Eq. (7)

$$\left[\begin{array}{c}\partial \mathrm{ln}{\mathbf{y}}_{\mathrm{I}690}\\ \partial {\mathbf{y}}_{\mathrm{\Theta}690}\\ \partial \mathrm{ln}{\mathbf{y}}_{\mathrm{I}830}\\ \partial {\mathbf{y}}_{\mathrm{\Theta}830}\end{array}\right]=\left[\begin{array}{cc}\frac{\partial \text{\hspace{0.17em}}\mathrm{ln}\text{\hspace{0.17em}}{\mathbf{y}}_{\mathrm{I}690}}{\partial {\mu}_{a,690}}\xb7{\epsilon}_{{\mathrm{HbO}}_{2},690}& \frac{\partial \text{\hspace{0.17em}}\mathrm{ln}{\text{\hspace{0.17em}}\mathbf{y}}_{\mathrm{I}690}}{\partial {\mu}_{a,690}}\xb7{\epsilon}_{\mathrm{HbR},690}\\ \frac{\partial {\mathbf{y}}_{\mathrm{\Theta}690}}{\partial {\mu}_{a,690}}\xb7{\epsilon}_{{\mathrm{HbO}}_{2},690}& \frac{\partial {\mathbf{y}}_{\mathrm{\Theta}690}}{\partial {\mu}_{a,690}}\xb7{\epsilon}_{\mathrm{HbR},690}\\ \frac{\partial \mathrm{ln}\text{\hspace{0.17em}}{\mathbf{y}}_{\mathrm{I}830}}{\partial {\mu}_{a,830}}\xb7{\epsilon}_{{\mathrm{HbO}}_{2},830}& \frac{\partial \text{\hspace{0.17em}}\mathrm{ln}\text{\hspace{0.17em}}{\mathbf{y}}_{\mathrm{I}830}}{\partial {\mu}_{a,830}}\xb7{\epsilon}_{\mathrm{HbR},830}\\ \frac{\partial {\mathbf{y}}_{\mathrm{\Theta}830}}{\partial {\mu}_{a,830}}\xb7{\epsilon}_{{\mathrm{HbO}}_{2},830}& \frac{\partial {\mathbf{y}}_{\mathrm{\Theta}830}}{\partial {\mu}_{a,830}}\xb7{\epsilon}_{\mathrm{HbR},830}\end{array}\right]\left[\begin{array}{c}\partial \mathrm{HbO}\\ \partial \mathrm{HbR}\end{array}\right].$$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.

## Table 1

Baseline optical properties in the head model.35

Region | Oxyhemoglobin (mM) | Deoxyhemoglobin (mM) | Scatter power | Scatter amplitude |
---|---|---|---|---|

Skin | 0.057 | 0.031 | 1.16 | 0.53 |

Skull | 0.044 | 0.019 | 0.89 | 0.73 |

CSF | 0.011 | 0.008 | 0 | 0.3 |

Gray matter | 0.056 | 0.035 | 1.74 | 0.51 |

White matter | 0.068 | 0.027 | 1.31 | 0.82 |

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:

## Eq. (8)

$${\mathbf{J}}^{-1}\approx {\mathbf{J}}^{\#}={\mathbf{J}}^{\mathrm{T}}{({\mathbf{JJ}}^{\mathrm{T}}+\alpha \mathbf{I})}^{-1},$$^{36}when running simulations without noise, fine-tuning of this parameter can improve image resolution beyond an experimentally and physiologically meaningful range. In practice, this parameter should be defined with respect to the noise of the measurements,

^{10}therefore all simulations in this work use the same regularization weight of $a=0.01$, following previously reported methods.

^{1}

For the FD Jacobian, Eq. (5), the regularization is typically applied separately for each of the kernels of log intensity and phase measurements.^{32} 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^{37} 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\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}$, $1\times 1\times 1\text{\hspace{0.17em}\hspace{0.17em}}{\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\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{mm}}^{3}$, resulting in $\sim \mathrm{450,000}$ nodes for each mesh using NIRFAST.^{39} A sparse-sensitivity matrix, $J$, (sensitivity threshold set to 0.1%^{40}), 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).

## 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.^{5} 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.^{41} 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.^{1} 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.

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.^{12} 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.

Amplitude noise model | Phase noise model | |||
---|---|---|---|---|

Coefficients | 830 nm | 690 nm | 830 nm | 690 nm |

$a$ | 0.6019 | 0.2502 | 1.917e-10 | 3.933e-11 |

$b$ | 0.01052 | 0.02913 | 0.3708 | 0.4161 |

$c$ | 9.685e-05 | 4.625e-06 | 0.03573 | 0.0105 |

$d$ | 0.1382 | 0.2128 | 0.02002 | 0.05585 |

## Table 3

The estimated noise levels for each measurement neighborhood.

Amplitude noise (%) | Phase noise (deg) | |||
---|---|---|---|---|

NN | 830 nm | 690 nm | 830 nm | 690 nm |

NN1 | 0.69 | 0.36 | 0.044 | 0.021 |

NN2 | 0.81 | 0.56 | 0.062 | 0.05 |

NN3 | 0.91 | 0.77 | 0.076 | 0.088 |

NN4 | 1 | 0.95 | 0.088 | 0.126 |

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\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{M}$ while for HbR was $-1.8\text{\hspace{0.17em}\hspace{0.17em}}\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:^{10} (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.^{42} 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,^{1} 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.^{43} 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\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{Hz}$. Prior to data analysis, measurements with high noise (standard deviation of the logmean of the data larger than 7.5%^{1}) 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.

## Results

## 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.

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,^{42} 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.

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.

## Table 4

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

Depth (mm) | FWHM (mm) | Localization error (mm) | Effective resolution (mm) | ||||||
---|---|---|---|---|---|---|---|---|---|

NN2 | NN3 | NN4 | NN2 | NN3 | NN4 | NN2 | NN3 | NN4 | |

5 to 10 | ${11.6}_{-10.3}^{+13.7}$ | ${11}_{-9.9}^{+12.4}$ | ${10.6}_{-9.5}^{+12.1}$ | ${2.2}_{-1.7}^{+2.9}$ | ${1.4}_{-1}^{+1.8}$ | ${1.1}_{-0.8}^{+1.5}$ | ${14.9}_{-13.2}^{+17.7}$ | ${13.2}_{-12}^{+15}$ | ${12.5}_{-11.2}^{+14.2}$ |

10 to 15 | ${13.3}_{-11.7}^{+15.7}$ | ${12.9}_{-11.6}^{+14.5}$ | ${12.5}_{-11.2}^{+14.1}$ | ${4.2}_{-3.3}^{+5}$ | ${2.6}_{-2}^{+3.2}$ | ${1.5}_{-1}^{+2}$ | ${19.6}_{-17.3}^{+22.5}$ | ${16.8}_{-15}^{+19.1}$ | ${14.9}_{-13.3}^{+17}$ |

15 to 20 | ${14.6}_{-12.6}^{+17.2}$ | ${15.2}_{-13.4}^{+17.5}$ | ${14.6}_{-12.8}^{+16.8}$ | ${6.4}_{-5.6}^{+7.2}$ | ${4.4}_{-3.6}^{+5.4}$ | ${2.8}_{-2.2}^{+3.6}$ | ${24.5}_{-22.3}^{+27.1}$ | ${22}_{-19.8}^{+25.3}$ | ${19}_{-16.8}^{+22}$ |

20 to 25 | ${14.5}_{-12.7}^{+17.1}$ | ${15.2}_{-12.8}^{+17.8}$ | ${15.3}_{-13}^{+17.8}$ | ${7.1}_{-6.1}^{+7.6}$ | ${5.8}_{-4.7}^{+6.7}$ | ${4.7}_{-3.7}^{+5.8}$ | ${25.1}_{-23.2}^{+28}$ | ${24}_{-21.6}^{+27.2}$ | ${23}_{-20.3}^{+26.1}$ |

## Table 5

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

Depth (mm) | FWHM (mm) | Localization error (mm) | Effective resolution (mm) | ||||||
---|---|---|---|---|---|---|---|---|---|

NN2 | NN3 | NN4 | NN2 | NN3 | NN4 | NN2 | NN3 | NN4 | |

5 to 10 | ${10.8}_{-9.8}^{+12.2}$ | ${10.2}_{-9.3}^{+11.5}$ | ${9.8}_{-8.9}^{+11.1}$ | ${1.1}_{-0.78}^{+1.4}$ | ${0.94}_{-0.68}^{+1.2}$ | ${0.87}_{-0.62}^{+1.2}$ | ${12.5}_{-11.3}^{+14.2}$ | ${11.7}_{-10.6}^{+13.3}$ | ${11.3}_{-10.2}^{+12.9}$ |

10 to 15 | ${13.1}_{-11.8}^{+14.6}$ | ${12.3}_{-11.1}^{+13.7}$ | ${11.3}_{-10.1}^{+12.9}$ | ${1.7}_{-1.2}^{+2.2}$ | ${1.4}_{-1}^{+1.8}$ | ${1.2}_{-0.83}^{+1.6}$ | ${15.6}_{-14}^{+17.5}$ | ${14.5}_{-13}^{+16.3}$ | ${13.3}_{-11.9}^{+15.3}$ |

15 to 20 | ${15.1}_{-13.7}^{+16.8}$ | ${14.9}_{-13.4}^{+16.6}$ | ${13.2}_{-11.5}^{+15.4}$ | ${3.4}_{-2.6}^{+4.2}$ | ${2.4}_{-1.8}^{+3.1}$ | ${1.7}_{-1.3}^{+2.3}$ | ${20.1}_{-18.2}^{+22.3}$ | ${18.7}_{-16.9}^{+21}$ | ${16.3}_{-14.3}^{+19.1}$ |

20 to 25 | ${15}_{-12.9}^{+17}$ | ${15.3}_{-13.4}^{+17.4}$ | ${14.8}_{-12.7}^{+17.4}$ | ${5.4}_{-4.3}^{+6.5}$ | ${4}_{-3.1}^{+5}$ | ${2.7}_{-2}^{+3.5}$ | ${23.1}_{-21.2}^{+25.4}$ | ${21.9}_{-19.6}^{+24.5}$ | ${20}_{-17.1}^{+23.4}$ |

## Table 6

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

FWHM (%) | Localization error (%) | Effective resolution (%) | |||||||
---|---|---|---|---|---|---|---|---|---|

Depth (mm) | NN2 | NN3 | NN4 | NN2 | NN3 | NN4 | NN2 | NN3 | NN4 |

5 to 10 | 7 | 7 | 7.5 | 51.5 | 32.2 | 20.7 | 16.1 | 11 | 9.5 |

10 to 15 | 1.9 | 4.8 | 9.3 | 58.9 | 47.3 | 21.1 | 20.6 | 13.7 | 10.5 |

15 to 20 | $-\mathbf{4}$ | 2.4 | 9.6 | 47.2 | 44.8 | 37.6 | 17.9 | 15.7 | 14 |

20 to 25 | $-\mathbf{4}$ | $-\mathbf{1.3}$ | 3.5 | 23.5 | 30.5 | 43.3 | 7.9 | 9.1 | 13 |

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.

## 4.

## Discussion

To date, most HD-DOT imaging studies have concentrated on the utilization of CW data. Although some studies have reported the use of FD NIR systems for spectroscopic mapping of brain functions, these have been mainly limited to utilization of phase and amplitude data for the separation of optical scattering components from absorption. We have recently demonstrated the information content of phase data from a retinotopic experiment showing that the phase data do contain rich information, potentially arriving from deeper regions within the brain and therefore minimizing the superficial tissue contamination.^{44} In this work, we have undertaken a detailed study of the potential improvements to image quality gained through incorporating measurements of the phase along with intensity in the reconstruction process.

Utilizing a simple layered slab model, the Jacobean (i.e., the sensitivity of both log amplitude and absolute phase to absorption changes) has been calculated and shown in Fig. 1. Several features are evident: (1) both log-amplitude and phase data show a decrease in signal due to a change to absorption; (2) while the sensitivity for both log-amplitude and phase data increases in depth of sampling with SD separation, the phase data demonstrates deeper sensitivity for a given SD separation; and (3) the hypersensitivity in superficial tissues observed with log amplitude data is less pronounced for phase data. This supports our earlier findings that phase data are able to detect changes deeper within the brain and are less prone to be contaminated to changes within the superficial tissue that may mainly be attributed to systemic cardiac and respiratory signal.

Utilizing a set of 24 subject specific models, Fig. 2 highlights that the surface of the cortex can be from 5 to 30 mm deep with respect to the surface of the head, further highlighting the need for the utilization of data-types that will allow deeper sampling of the brain tissue. Although this can, in theory, be achieved through the use of larger SD separations, in reality, this is not possible due to system design and low SNR at measurement distances above 45 mm.^{1}^{,}^{43} To demonstrate that FD data-types (log intensity and phase) can provide a better sampling of deeper tissue, the subject-specific models have been utilized to simulate a set of data corresponding to focal activations at all cortex surface nodes lying underneath the imaging array, highlighted as pink regions, Fig. 2(a). In order to ensure that these simulated data are realistic, the noise characteristics of a commercially available FD system, ISS Imagent have been measured, as a function of SD-NN measurements (Fig. 3) and were added to the simulated data (Table 2). To investigate the benefits of using differing NN (NN2-NN4) measurements for activations at different depths, three sets of representative focal activation recovery maps are shown in Fig. 5. In all cases as the number NN measurements increases, the depth recovery improves. In each case, the presence of noise is shown to introduce artifacts within the recovered maps. Importantly, the utilization of FD data is shown to improve all metrics of image quality relative to using CW. In almost all cases, the localization accuracy of utilizing FD data is shown to be much better than those of CW, even including an activation at 20.5 mm where the CW data failed to recover a meaningful map.

To provide an overview of these benefits, images of focal activations for all noise-added simulated cortical activations have been reconstructed and analyzed to provide group average metrics of localization accuracy (Fig. 7) full-width-half-max (Fig. 8) and effective resolution (Fig. 9) with data summarized in Tables 4–6. For all metrics in the noise-free case, the accuracy decreases as a function of activation depth, with FD data outperforming CW in all cases. To achieve meaningful statistical comparisons, a localization error limit of 8 mm has been imposed as the upper-bound for statistical analysis, excluding all recoveries that are solely products of noise. For all NN measurements, the FD data provide a significant improvement in almost all recovered activations for all reported metrics (summary in Table 6). For example, assuming a set of NN3 measurements, it is shown that as compared to CW, the FD data show up to 59% improvement in localizing accuracy and up to 21% improvement in effective resolution. In addition, the FD data using NN4 are observed to successfully recover activations as deep as 20 to 25 mm, beyond the range accessible to CW.

Finally, to demonstrate the improvements using experimental data, a single focal activation from a set of retinotopic data for a single subject is shown in Fig. 10. Utilizing the same reconstruction algorithms, with the same regularization parameters, it is demonstrated that the inclusion of phase data from an FD system is capable of recovering the focal activation within the visual cortex with much more realistic location accuracy.

It is important to highlight that utilizing a time-resolved system for NIRS, utilization of mean-time of photon travel has been previously reported,^{45} which is analogous to the phase for an FD system. However in this work, the utilization of both log intensity and phase for tomographic reconstructions have been highlighted to demonstrate benefits. Although the sensitivity of the log intensity is shown to be more weighted toward the superficial tissues, it is this variation of sensitivity between log intensity and phase data that is providing the additional accuracy as demonstrated.

## 5.

## Conclusions

This work has examined the effect and potential benefits of employing frequency-domain measurements for functional DOT, using a high-density measurement array. Simulated data using 24 subject-specific models utilizing experimentally measured noise have been employed to provide a set of metrics for characterization of recovered image quality, as well as an *in-vivo* retinotopy experiment from a single subject. The results reveal that, as expected, resolution in the depth is related to the SD distance and level of the noise present, which are in agreement with previous studies.^{6} In addition, using the phase information significantly improves image quality when employing the same number of NN measurements with or without noise. While the use of the phase significantly improves image quality throughout the volume, the improvement is larger in the deeper regions, as phase samples deeper than amplitude for the same measurement. As such, the FD measurements can be used to resolve activations in deeper regions, where the CW completely fails. Overall, this study suggests that the use of FD measurements in the context of HD-fDOT expands the potential for improvements in image quality beyond the current state of the art.

## Acknowledgments

This work has been funded by the National Institutes of Health (NIH) Grant Nos. R01EB009233-2, RO1-CA132750, NIMH-R21-109775, and NIHMH-K01-103594 and in part by Marie Sklodowska-Curie Innovative Training Networks (ITN-ETN) programme, under Grant Agreement No. 675332 (BitMap).

## References

## Biography

**Matthaios Doulgerakis** is received his BEng degree in electronic engineering from Alexander Technological Educational Institute, Thessaloniki, Greece, in 2013, and his MRes in medical imaging from Kingston University, London, United Kingdom. He is a PhD candidate developing real-time algorithms for functional diffuse optical tomography (DOT) at the University of Birmingham, United Kingdom. His research interests are inverse problems and optimization algorithms in medical imaging, particularly in diffuse optics.

**Adam T. Eggebrecht** received his PhD in physics in 2009 and is currently an assistant professor in the Mallinckrodt Institute of Radiology at Washington University School of Medicine. He research is focused on designing and applying diffuse optical tomography (DOT) imaging methods in the human brain to understand mechanisms of typical and atypical functional brain development. He directs development of the NeuroDOT software package for modeling and analyses of diffuse optical measurements of brain function. He is currently funded by an NIH career development award to optimize and apply DOT to the study of altered brain function in children with autism. He is also funded to further develop DOT for studies of brain development in neonates, infants, and toddlers, as well as patients with implanted deep brain stimulators.

**Hamid Dehghani** is a professor of medical imaging in the School of Computer Science at the University of Birmingham, United Kingdom, and an OSA fellow. He has published over 100 peer-reviewed journal papers in the area of image reconstruction and numerical modeling and has a long and established track record in the development of biophotonics based techniques with specific applications in *in-vivo* optical imaging.