High-density functional diffuse optical tomography based on frequency-domain measurements improves image quality and spatial resolution

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


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. 6Recent 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. 7NIRS 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. 8Differential measurements at two NIR wavelengths, known as NIR spectroscopy (NIRS), can recover oxy (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 sourcedetector (SD) separation distances has been shown to improve image quality 10 to a point comparable to fMRI. 1,11In 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. 13hen 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. 14Measurements 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. 15FD-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][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,26lthough 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 eventrelated optical signal or fast optical signal, uses just the phase measurements of FD systems and has been reported to be sensitive to neuronal activity. 27It 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 (∼60 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. 29This 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. 30Therefore, 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. 31o 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.

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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 3 2 6 ; 5 4 0 Φðr; ωÞ ¼ qðr; ωÞ; (1) where r defines a spatial location, ω is the modulation frequency, the diffusion coefficient being the reduced scattering coefficient, Φ is the photon fluence, μ a is the absorption coefficient, c m is the speed of light in the medium, and q denotes a light source.High-resolution meshes (mean mesh element volume <1 mm 3 ) have been used throughout this work, which improve the accuracy of the forward model at a significant increase in the computational cost. 25To 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 ∼45;000 nodes, with 158 sources and detectors, the forward model calculation time is less than 3 min. 33

Inverse Problem
The sensitivity matrix, J, also known as the Jacobian or weight matrix is calculated using the Adjoint method 34 that relates changes in boundary fluence measurements, ∂y, with respect to small changes in underlying tissue optical parameters, ∂x: ; t e m p : i n t r a l i n k -; e 0 0 2 ; 3 2 6 ; 2 6 5 From Eq. (1), in the frequency domain, the measured data are complex numbers and can be decomposed in two parts: amplitude y I , and phase y Θ .Using the Rytov approximation, differential measurements with respect to a baseline y 0 are described as a logarithmic ratio for amplitude and a difference in phase: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 3 2 6 ; 1 7 9 ln y y 0 ¼ ln y I e iy Θ y I0 e iy Θ0 ¼ ln Therefore, the frequency-domain Jacobian describes changes in logarithmic intensity, ∂ ln y I , and linear differences in phase ∂y θ , with respect to changes in the absorption, ∂μ a , and the diffusion, ∂κ 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, ∂μ a , alone, reducing Eq. ( 4) to E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 6 3 ; 6 4 and using measurements at two wavelengths of 690 and 830 nm, the dynamic changes of oxy (HbO 2 ) and deoxy (HbR) hemoglobin concentration can be calculated using 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.
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), J # , to perform a single step linear recovery of absorption: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 6 3 ; 3 0 2 ∂x ¼ J # ∂y; where the regularization term is defined as α ¼ a Ã max½diagðJJ T Þ.The weight of this regularization parameter, a, provides a tuning between high spatial frequency noise and smoothness in the image.Due to this association with image resolution, 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. 32In this work, the absorption coefficient is recovered for each wavelength then spectrally unmixed to recover 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.

Subject-Specific Models
The subject-specific models are based on 14 female and 10 male individual subjects with a mean age of 26 (AE4), with T1weighted MPRAGE [echo time = 3.13 ms, repetition time (TR) = 2400 ms, flip angle ¼ 8 deg, 1 × 1 × 1 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 ∼1 mm 3 , resulting in ∼450;000 nodes for each mesh using NIRFAST. 39A 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 μ 0 s ¼ s a λ −s p , where s a is the scatter amplitude, s p is the scatter power, and λ 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 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).

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. 5The 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. 41The 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. 1Briefly, 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Þ ¼ ae br þ ce 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.
To simulate realistic measured data, perturbations of HbO 2 and HbR values were induced using Eq. ( 7), using perturbation values scaled to minimize the difference between simulated and 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.measured data.The perturbation for HbO 2 was scaled to be 3.8 μM while for HbR was −1.8 μ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.

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. 42This 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.

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. 43Briefly, 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 ∼40 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 HbO 2 and HbR were performed as defined above, using NIRFAST for the calculation of the Jacobian matrix for both CW and FD data.

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 HbO 2 and HbR have been reconstructed and as both demonstrated similar results, for conciseness, only those of 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.
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 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 Wilcoxonsigned 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.

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

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. 44In 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,43To 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.

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. 6In 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.

Fig. 1
Fig. 1 Frequency-domain Jacobian for Log intensity (Joules∕mm 2 ∕mM) and phase (radians/mM) for 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 Table1.

Fig. 2
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).

Fig. 4
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. 6
Fig. 6 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.

Fig. 7
Fig. 7 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
Fig. 8 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
Fig. 9 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 2
The estimated coefficients for the two-term exponential fit of the noise model.

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).
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.Doulgerakis, Eggebrecht, and Dehghani: High-density functional diffuse optical. . .