Source Localization Approach for Functional Dot Using Music and Fdr Control References and Links

In this paper, we formulate diffuse optical tomography (DOT) problems as a source localization problem and propose a MUltiple SIgnal Classification (MUSIC) algorithm for functional brain imaging application. By providing MUSIC spectra for major chromophores such as oxy-hemoglobin (HbO) and deoxy-hemoglobin (HbR), we are able to investigate the spatial distribution of brain activities. Moreover, the false discovery rate (FDR) algorithm can be applied to control the family-wise error in the MUSIC spectra. The minimum distance between the center of mass in DOT and the Montreal Neurological Institute (MNI) coordinates of target regions in experiments was between approximately 6 and 18mm, and the displacement of the center of mass in DOT and fMRI ranged between 12 and 28mm, which demonstrate the legitimacy of the DOT-based imaging. The proposed brain mapping method revealed its potential as an alternative algorithm to monitor the brain activation. 1. A. Villringer and B. Chance, " Non-invasive optical spectroscopy and imaging of human brain function, " Trends Noninvasive functional imaging of human brain using light, " J. Cereb. Intrinsic signal changes accompanying sensory stimulation: Functional brain mapping with magnetic resonance imaging, " Proc. Nati. Three-dimensional whole-head optical tomography of passive motor evoked responses in the neonate, " NeuroImage 30, 521–528 (2006). Three-dimensional optical tomography of hemodynamics in the human head, " Opt. Diffuse optical imaging of brain activation: approaches to optimizing image sensitivity, resolution, and accuracy, " NeuroImage 23, S275–S288 (2004). Depth sensitivity and image reconstruction analysis of dense imaging arrays for mapping brain function with diffuse optical tomography, " Appl. Brain specificity of diffuse optical imaging: improvements from superficial signal regression and tomography, " Front Neuroenergetics 2, 1–8 (2010). Retinotopic mapping of adult human visual cortex with high-density diffuse optical tomography, " Proc. Nat. " High-resolution optical functional mapping of the human somatosensory cortex, " Front Neuroenergetics 2, 1–8 (2010). Resting-state functional connectivity in the human brain revealed with diffuse optical tomography, Comprehensive investigation of three-dimensional diffuse optical tomography with depth compensation algorithm, " J. Hierarchical Bayesian regularization of reconstructions for diffuse optical tomography using multiple priors, " Biomed. Level-set algorithm for the reconstruction of functional activation in near-infrared spectroscopic imaging, " J. Simulation study of magnetic resonance imaging-guided cortically constrained diffuse optical tomography of human brain function, " Appl. Time-series estimation of biological factors in optical diffusion tomography, " Phys. A haemodynamic response function model in spatio-temporal diffuse optical tomography, " Phys. Compressive MUSIC: revisiting …


Introduction
Diffuse optical tomography (DOT) is a non-invasive and low-cost imaging modality that reconstructs optical parameters of cross sections of a highly scattering medium from the measurements of scattered and attenuated optical flux.Due to the relative transparent optical window between 700 and 1000nm [1], near-infrared (NIR) photon can penetrate several centimeters, which makes DOT a complementary tool in bio-medical imaging applications such as breast cancer detection [2], molecular imaging [3], functional brain imaging [4], etc.Recently, many researchers have investigated brain imaging applications of DOT, since it provides rich physiological information such as total hemoglobin (HbT), oxy-hemoglobin (HbO) or deoxyhemoglobin (HbR) for monitoring brain activities, compared to functional magnetic resonance imaging (fMRI) that measures only blood oxygen level dependent (BOLD) signals [1,5].Moreover, DOT has a better temporal resolution than that of fMRI (< 0.5Hz) [6], and the equipment for a DOT system is portable and radioactivity-free.Therefore, it is appropriate for routine bed-side monitoring.
Gibson et al. [7] reconstructed three-dimensional concentration changes in HbT, HbO, and HbR for the passive motor activation in neonate.Bluestone et al. [8] visualized concentration changes in human heads during Valsalva maneuver experiments.Currently, more sophisticated functional mapping has been also demonstrated by using high definition DOT (HD-DOT) systems [9][10][11].HD-DOT uses remote pairs to capture the signal from deep inside the brain, and nearby pairs to remove the signal originating from superficial tissue.Using HD-DOT, Zeff et al. [12] demonstrated the retinotopic mapping of visual cortex in an adult human brain, and Koch et al. [13] successfully differentiated various tasks in the human somatosensory cortex such as finger tapping and vibrotactile stimulus on the first or fifth finger.Furthermore, not only task-based brain imaging, but also resting-state functional connectivity has been studied using DOT [14].
However, the major drawback of DOT is its ill-posedness due to the diffusive nature of light propagation.Therefore, extensive investigations have been performed to resolve this problem.Niu et al. [15] proposed a depth compensation algorithm that exploits singular values of the partitioned forward matrix to compensate for exponentially decreasing sensitivity with depth.Abdelnour et al. [16] have implemented a hierarchical Bayesian regularization, and Jacob et al. [17] proposed a level-set algorithm by imposing additional constraints such as smoothness of the activation pattern and sparseness of the activation area.Boas et al. [18] used cortical constraint to reduce a region of interest and enhance the depth sensitivity.All these methods use spatial constraints to remedy the ill-posedness.
As temporal changes of physiological parameters like HbO or HbR are highly correlated with pre-defined paradigms, we can reduce the noise level of the signal and make the reconstruction more stable by considering this temporal information.Kolehmainen et al. [19] modeled the absorption variations in a multivariate stochastic process and formulated it as a state-estimation problem.They used a Kalman filtering technique to solve the problem and demonstrated the validity in the simulation study.Prince et al. [20] also suggested a similar method by considering the absorption variations as a mixture of quasi-sinusoidal signals and applied it in real experimental data.Zhang et al. [21] proposed a fully spatio-temporal model using general linear model (GLM) constraint in the DOT forward problem to describe the concentration changes in HbO and HbR.
Unlike the existing methods, the main contribution of this paper is to formulate DOT problems as a source localization problem by assuming that neural activation is relatively localized during experiments.For example, in a resting state experiment or event related paradigm, application of GLM as a temporal constraint is not always feasible due to the lack of an accurate paradigm.However, even in this case, we can assume that the activation area determines statistical significance of the experiment.Similarly, in a block paradigm or event related paradigm, target areas are expected to activate regularly according to a given paradigm.Therefore, we can expect a new spatio-temporal regularization scheme by exploiting simple statistics of the temporal data rather than assuming accurate knowledge of the temporal time trace.
Recall that a source localization approach based on sensor array signal processing [22] uses an array of sensors to localize incoming sparse sources by exploiting the second order statistics of the time traces measured through multiple detectors.Since a DOT system usually has multiple detectors that measure time series of the optical signals and the activation area is usually sparsely distributed, sensor array signal processing algorithm can be used.Furthermore, in our recent study, a close link between a source localization problem and compressed sensing [23,24] has been revealed, which demonstrates the array signal processing is a kind of compressed sensing that exploits spatial domain sparsity of the temporal signals [25].Based on this observation, we already demonstrated a sensor array signal processing approach for small animal imaging, which exploits spatial domain diversity by changing source excitation patterns [26].However, brain imaging is a little bit different from small animal imaging since we are mainly interested in relative temporal dynamics of the signals rather than the absolute value at a given time point.Moreover, the number of source-detector pairs is much smaller than small animal imaging and the illumination pattern cannot be changed as done in the study by Lee et al. [26].Therefore, instead of using those particular techniques [26], we use a multiple signal classification (MUSIC) algorithm [27,28] which exploit the eigen decomposition of the temporal covariance matrix to identify the signal and noise subspace.Then, the main idea of MUSIC is that the atoms at the signal location are orthogonal to the noise subspace.Since the time series of DOT is sufficiently long due to the dense temporal sampling, the empirical covariance matrix can approach to the true covariance matrix, hence the estimated noise and signal subspace can converge to the true ones.Hence, MUSIC can identify the true activation consistently.Furthermore, to enable inferences that are required for neuroimaging, we propose a false discovery rate (FDR) controlling scheme [29] by assuming that the resulting signal subspace version of MUSIC spectrum follows the chi-squared distribution.
We carried out experiments with the right finger tapping task based on a block paradigm to investigate the localized activation in the left primary motor cortex (BA4) and somato-sensory cortex (BA123).In order to evaluate the proposed method quantitatively, we computed the minimum distance between the center of mass of each MUSIC spectrum and the Montreal Neurological Institute (MNI) coordinate of the BA4 or BA123, as well as the displacement of the center of mass in activated areas of DOT and fMRI.These comparison procedures were performed using both synthetic and real right finger tapping experimental data so that we elucidate whether the proposed approach can extract the spatial distribution of brain activities.

Notation
The following notation is used in this paper: • m : number of source and detector pairs (channel); • J : number of time points; • λ : illumination wavelength; • n : number of descretized voxels; • r;r i : position vector; position of the i-th voxel; • r s : position vector of the source; • r d : position vector of the detector.

Forward model
An optical flux variation at a detector position r d from a source location at r s can be approximated using the following Rytov approximation: where Δμ a (r; λ ,t) is an absorption variation at r, λ denotes the illumination wavelength, and U 0 (r, r ; λ ) is the optical flux at r generated from a source located at r .By collecting optical density variations for pairs of detectors and sources { r d i , r s i } m i=1 , Eq. ( 1) can be represented in a matrix form for a given wavelength λ at time t j : where where {r i } n i=1 denotes voxel positions, and w λ j ∈ R m is a noise vector for wavelength λ at time t j .Here, the (i, j)-th element of the sensing matrix A λ = {A λ i j } m,n i, j=1 is given by where δ denotes a discrete voxel volume.If we assume that dynamically varying absorption contrast is mainly due to the oxygenation level changes in hemoglobin, we have where ε λ HbO [μM −1 mm −1 ] and ε λ HbR [μM −1 mm −1 ] are the extinction coefficients of HbO and HbR at wavelength of λ , and Δc HbO (r;t) and Δc HbR (r;t) are the concentration changes of HbO and HbR in a voxel r at time t, respectively.Using measurements from two wavelengths {λ 1 , λ 2 }, the forward measurement model Eq. ( 2) can therefore be represented as where the sensing matrix A is and x HbR, j = [Δc HbR (r 1 ;t j ), Δc HbR (r 2 ;t j ),

Source-localization formulation
By collecting all temporal series ( j = 1, 2, • • • , J), Eq. ( 5) can be formulated as where From this formulation, Appendix A describes how this problem can be solved using the socalled MUSIC algorithm.Furthermore, under the Gaussian assumption of the sensing matrix, the following MUSIC spectrum can be assumed to have chi-squared distribution, Therefore, we can plot the spectrum across for all voxel index j and find the appropriate maximum after thresholding to detect the activation.
In applying MUSIC for DOT problems, there are two approaches to detect activation.First, we can assume that HbO and HbR activation are independent of each other.In this case, we need to calculate the separate MUSIC spectrum for HbO and HbR, respectively.Actual implementation of this is relatively simple, since in our forward model in Eq. ( 9), the unknown ambient space dimension is 2n, where one n is for HbO and the other n for HbR.Hence, rather than using voxel space domain as the sparse source location, we plot MUSIC spectrum along the 2n-dimensional index space.Then, the first half of the MUSIC spectrum corresponds to HbO and the remaining half of the MUSIC spectrum is for HbR.More specifically, we have Second, we assume that there are correlations in HbO and HbR in activated areas.This is usually true since neural activation generally increases the HbO level while reducing HbR levels.In this case, our source localization problem has block-sparse structure as shown in Fig. 1.Such a block-sparse problem has been studied in the context of EEG/MEG source localization Fig. 1.Block-sparse MMV model when HbO and HbR are assumed to be correlated with each other.The shaded areas describe the index of simultaneously activated voxels of HbO and HbR, respectively.[30], where the dipole moment direction is assumed unknown and we need to estimate both the directional vector as well as magnitude.For this type of model, the extended version of MUSIC spectrum can be written as where • F denotes the Frobenius norm.By plotting Eq. (11) or Eq. ( 12), we can investigate the spatial distribution of brain activities related with the specific task.

Signal subspace selection
We obtained thirty-four singular values by using DOT data from seventeen channels and two different wavelengths, 781 and 856nm.Among the singular values, the first few dominant singular values are related to signal subspace as shown in Fig. 2. Here, the number of singular values included in signal subspace corresponds to dimension of signal subspace.Also, as described in Appendix A, the dimensions of signal subspace determines the degrees of freedom (df) in statistics (chi-squared distribution) derived from MUSIC spectrum.More specifically, df is equal to the dimension of the signal subspace in Eq. (11).In a case where HbO and HbR are assumed to be correlated with each other in activation region (Eq.( 12)), df is determined as the sum of the dimension of the signal subspace of HbO and HbR, respectively.Under the assumption that the neural activation is sparsely localized during experiments, we can expect the hemoglobin concentration changes from the activated area are highly correlated resulting in very low dimension of the signal subspace (low value of df).Therefore, we selected the dimension of the signal subspace as one or two (df=1,2) in simulation and experimental studies.Figure 3 shows the example of MUSIC spectrum according to various choices of signal subspace dimension (Dim.= 1, 3, 5).As described in the figure, the spatial distribution of MUSIC spectrum is dependent on the choice of signal subspace.

Inference using MUSIC spectrum
One of main goals of statistical analysis in functional neuroimaging is to determine whether the statistics represent evidence of certain effects.A typical approach is one in which the statistics are tested against a null hypothesis (i.e., no activation).By conducting a statistical test against a null hypothesis, we can determine how likely the statistics occur by chance or by meaningful activation.
A recent approach known as false discovery rate (FDR) controls the expected proportion of falsely declared-active detections among the total declared-active hypotheses [29].It has been adopted in functional neuroimaging as an alternative approach for the inference of brain activation [31,32].As shown in Table 1, the total n voxels are classified into one of four types.Then the FDR is defined as follows: which is the proportion of the false-positives (FP) among only the rejected null hypotheses (FP+TP(true-positives)).Hence, controlling the FDR provides higher localizing power than other multiple testing methods such as Bonferroni correction which assess the entire null hypotheses (N) [31,32].Total declared Inference procedure based on FDR consists of the following three steps.First, specify a desired FDR q (0 ≤ q ≤ 1) which ensures that the expected FDR is less than or equal to q: Here, to compute the threshold based on FDR method and associated q-value, one needs uncorrected p-values, i.e. p i = Pr{V ≥ ν(i)}, where V denotes an underlying random variable for MUSIC spectrum.Secondly, sort the uncorrected p-values into ascending order, p 1 ≤ p 2 ≤ . . .≤ p n which correspond to null hypotheses, H 1 , H 2 , . . ., H n .The final step is to evaluate the following inequality in reverse order, Let k be the largest i which satisfies the Eq. ( 15), then we reject all the hypotheses, H 1 , H 2 , . . ., H k .That is, we threshold the MUSIC spectrum at the p k value, and declare the corresponding voxels active.Hence, we investigate the significantly active voxels controlling the expectation of FDR less than q-level.
The FDR algorithm adapts its threshold to the characteristics of the functional neuroimaging data; therefore, it provides consistent inference results across various data without an excessive removal of family-wise error.The intuitive and clear definition of the pre-specified q-value makes the inference procedures simple and easily accepted throughout several studies.Moreover, it gives us higher localizing power than Bonferroni correction while still adjusting the balance with specificity.

Simulation setup
We performed a simulation study for localization of the activation area using the T1 image obtained from a fMRI experiment.Figure 4(a) illustrates the source and detector geometry, which is identical with the one adopted in the real NIRS system (Oxymon MK III, Artinis, Netherlands).The activation area was assumed to be sphere-shaped and placed in the cortical layer.Among the twenty-four channels, we removed the most upper part of seven channels since those were not aligned with fMRI data.The DOT signals from the eliminated channels can be neglected since those channels were placed far from the left primary motor cortex and somato-sensory cortex (The results of using all channels were similar with the removed one).These remaining channels are widely spread over the left hemisphere.Since there are no overlapping and no crossing photon paths for an accurate lateral localization of the activated area, we calculated lateral and depth distance error in simulation and experimental studies to quantify these effects.Since most of the channels are almost placed in parallel with sagittal plane, we considered the distance error on this plane is lateral, and the one on the axis perpendicular to the sagittal plane is depth.Figure 4(b) describes the synthetic concentration changes for HbO and HbR in the activation area, in which Gaussian noise with 30dB signal to noise ratio (SNR) was added.Location of the centers for synthetic activation with a radius of 3mm on the horizontal plane during the simulation are illustrated in Fig. 4(c).In Fig. 4(c), the activation area varies in lateral direction from 0 ∼ 28mm for a given depth, and we changed the depth with 28 ∼ 48mm.The position over 40mm depth seems to deep for NIRS, however, there is a result that a photon can reach even white matter [33] so we want to check the sensitivity of our proposed method in this case.Measurements were obtained from 17 channels and Gaussian noise with 10dB SNR was added at each measurement.To test the inference using FDR control, we put an sphereshaped activation (the big arrow in the Fig. 4(a)) having a radius of 3mm at the same center position with the middle one of the most lateral line in Fig. 4(c).

Task protocol for in vivo experiment
To assess the validity of the proposed approach, right finger tapping tasks were carried out.We used a paradigm based on block design for the experiment.The main target region of these assignments is the left primary motor cortex (BA4) and somato-sensory cortex (BA123), which is located within the penetration depth of near-infrared light.Three, healthy, right-handed subjects were involved in this study, and each subject was instructed to perform a flexion-extension movement of right fingers repeatedly when a word 'Go' appeared until a word 'Stop' was shown.During the rest period, the subjects stared at a fixed point to avoid eye and head movements.As illustrated in Fig. 5, each block consisted of a 15-sec stimulation period followed by a 72-sec control period.This was repeated 4 times for each volunteer, which resulted in a total recording time of 468-sec including the preceding 90-sec and the following 30-sec rest periods.All volunteers were given instructions about the experimental process and all provided written informed consent.No subject had a history of neurological diseases.This study was approved by the Institutional Review Board of the Korea Advanced Institute of Science and Technology (KAIST).

Real data acquisition
Optical density variation was measured by a continuous wave NIRS instrument (Oxymon MK III, Artinis, Netherlands) that offers up to 24 channels with 8 sources and 4 detectors.This  instrument provides only measurements from the first nearest channel where the distance between source and detector is 3.5cm.The second nearest pair is about 7.8cm which is too far to measure the reasonable optical signal.The sampling frequency used for the experiment was 10Hz.Two kinds of wavelengths, 781nm and 856nm laser light, were emitted from each source fiber.A square-shaped holder cap with optodes was attached to the scalp around the primary motor cortex.For simultaneous recordings of DOT and fMRI, 10m optical fiber was used to connect the optodes in the MRI scanner with the NIRS system in the MRI control room.A 3.0T MRI system (ISOL, Republic of Korea) was employed to measure the blood oxygenation level dependent (BOLD) signal.The echo planar imaging (EPI) sequence was applied to acquire the functional images (with repetition time

Preprocessing
In NIRS research, there often occur global drifts of the optical density whose amplitude is comparable to that of the signals induced by brain activation [34].There have been also some results suggesting that in NIRS brain activation imaging, significant artifacts can be caused by changes in systemic and superficial blood circulation [35].Detailed analysis of these effects for activation localization is important and needs thorough study.However, this is beyond the scope of this paper, and it will be reported elsewhere.Instead, we employed the wavelet-MDL detrending method [36] to remove the global drifts.In addition, as suggested in fMRI [37], a low pass filter using the canonical hemodynamic response function (HRF) was also employed to temporally smooth the NIRS time-series, which is required to model the "short-range" temporal correlation that exists between the temporally-neighbored residual signals.We followed the generalized preprocessing procedures to perform MRI data analysis by using the SPM5 version of the statistical parametric mapping (SPM) software (http://www.fil.ion.ucl.ac.uk/spm, [38]).Specifically, the MR time-series were realigned in order to remove the movement artifact of volunteers.The data were then spatially normalized into the standard space of MNI coordinates and spatially smoothed with a Gaussian kernel whose full-width at half maximum (FWHM) was determined as 8mm.

Computation of sensing matrix
The elements of the sensing matrix in Eq. ( 3) were obtained by calculating photon flux using graphics processing unit (GPU) based Monte Carlo simulation [39], which was also adopted in Custo el al. [40].Toward this, we first segmented T1-weighted image into the superficial tissue (skin/skull), cerebrospinal fluid (CSF), gray matter and white matter using SPM5 toolbox in MATLAB (Mathworks, Natick, MA).It was then aligned with the standard MNI coordinates with real optode positions in a NIRS experiment.After the segmentation, we assigned optical parameters at each wavelength, 781nm and 856nm respectively, to each part of the brain as described at Table 2.The optical parameters were obtained from Okui et al. [41].Specifically, the parameters at 856nm were estimated based on the fact that a reduced scattering coefficient follows that of Lorenz-Mie scattering, and the absorption coefficient can be represented by a linear combination of the specific absorption of hemoglobin, oxy-hemoglobin, fat and water [41][42][43][44].We set the anisotropic factor g = 0.9 as the common value for the tissue [45,46], and the refractive index to 1 for MC simulation.We used 10 7 photons to calculate the Green's function U 0 (r, r s ; λ ) and U 0 (r d , r; λ ) (the reciprocity theorem [47] can be used to calculate U 0 (r d , r; λ )), and 10 8 photons to calculate the optical flux measured at detectors U 0 (r d , r s ; λ ).We assumed that neural activation occurs only in the brain cortical area located within 5cm from the optodes and set it as a field of view in DOT.Similar spatial constraint has been also used to improve the resolution of the DOT reconstruction [18].

Quantitative analysis
The minimum distance between the center of mass (COM) in DOT as well as in fMRI and the MNI coordinate of the BA4 or BA123 was calculated, which demonstrates the level of localization error.We also acquired the number of simultaneously activated voxels in DOT and fMRI to investigate the correspondence between them.Moreover, we quantified the displacement of the COM between fMRI and MUSIC spectra in DOT to look into the degree of their inconsistencies.COM is calculated for each activation area of DOT and fMRI that is included in both field of view (FOV).We applied the FDR correction to DOT data with q < 0.01 or 0.05, whereas the SPM5 toolbox was used to threshold fMRI while controlling the family-wise error with p < 0.01 or 0.05.For group analysis, all MUSIC spectra was added up in the overlap region since we used each individual T1-weighted image in MNI coordinate.Group df was determined as the sum of each individual df.

Synthetic data
Figures 6(a)-6(c) illustrate the lateral and depth displacement error, which is defined by the distance between the center of the ground truth activation and the peak of the MUSIC spectra (ν HbO and ν HbR ) from the simulation in Fig. 4(c).The results of ν HbO&HbR are similar with that of ν HbO or ν HbR (results are not shown).The results were averaged after 500 runs.We display the displacement error along the lateral direction when depth is 30mm and 44mm in Fig. 6(a) and 6(b), respectively.In Fig. 6(c), we averaged the entire error for the lateral directions for each depth.From this simulation studies, we confirmed that the difference between lateral and depth displacement error gets bigger when the activations are located deeper.Next, we fixed the location of the activation area at the middle position of the most lateral line in Fig. 4(c) and confirmed the thresholding performance by FDR control with q < 0.05 (df=2).
Figure 7(a) shows the thresholded MUSIC spectrum for HbO, in which the thresholded area is precisely located in the activation area.The cross hair of coronal, sagittal and horizontal section in Fig. 7(b) also demonstrate the declared activation is well localized compared to the ground truth in Fig. 4(a).The results from ν HbR and ν HbO&HbR , like that from ν HbO , also illustrate accurate localization (results are not shown).Fig. 7. Activation area for HbO using MUSIC (q < 0.05, corrected.df=2).The cross hair in coronal, sagittal and horizontal section indicates the position of the peak value of MUSIC spectrum.

In vivo real data
MUSIC spectra of each HbO and HbR, of HbO and HbR combination were acquired from right finger tapping experiment by using Eq. ( 11) and Eq.(12).Figures 8(a) and 9(a) show the inference results of subject 1 and subject 3 based on ν HbO&HbR and FDR control for q < 0.01 (df=1), respectively.In the figures, red colored area corresponds to the simultaneously activated region by fMRI and DOT, which covers around the left primary motor cortex and somatosensory cortex.Likewise, the 2D maps for DOT with coronal, sagittal, and horizontal views in Figs.8(b) and 9(b) shows that the significant voxels are localized tightly in the target regions.
Table 3 describes the minimum distance between the COM of DOT or fMRI and the MNI coordinates of BA4 or BA123.In most cases, the displacement ranges from 7 to 18mm for individual DOT.The COMs from group analysis was visualized as shown in Fig. 10, in which the distance up to BA4 or BA123 from COMs of DOT is comparable with the one from the COM of fMRI.
We quantified the number of activated voxels in DOT or fMRI together with the number of voxels which are placed within BA4 or BA123 among the whole activated voxels from each method (See Table 4).Note that the number of significant voxels from ν HbO&HbR is apparently greater than that from ν HbO or ν HbR , which demonstrates that the sensitivity of ν HbO&HbR is higher than others.
Table 5 illustrates the displacement of the center of mass between DOT and fMRI to investigate the level of disagreement between DOT and fMRI methods.There exists around 18 to 28mm of distance between each center of mass.As depicted in Fig. 10, the COM of group ν HbO&HbR is relatively closer to the COM of group fMRI.Depth and lateral displacement of the COM are also summarized in Table 6.Lateral error is generally bigger than that of depth   3. Minimum distance between the MNI coordinate of left primary motor cortex (BA4) or somato-sensory cortex (BA123) and the COM of (a) fMRI, (b)-(d) MUSIC spectra from HbO, HbR, combination of HbO and HbR, respectively; the degrees of freedom for MUSIC spectrum for each individual case was set as two in (b), (c), and as four in (d).In group analysis, each df was tripled, thereby six for (b), (c), and twelve for (d).

Conclusion and discussion
In this paper, we presented MUSIC algorithm for functional DOT in brain functional imaging.While MUSIC spectrum provides the spatial distribution of brain activities, we need an inference tool to detect activation area.To solve this problem, we applied the false discovery rate to control the family-wise error rate of the statistical map based on the assumption that the resultant MUSIC spectrum follows the chi-squared distribution in a large system limit.FDR algorithm enables us to remove the false positives among the number of voxels that are declared active while maintaining a given FDR q-level.
To validate the proposed method, we performed several quantitative comparison procedures, the minimum distance between the COM in DOT and the MNI coordinate of left primary motor cortex or somato-sensory cortex, and the displacement of the center of mass in DOT and fMRI.Even though DOT itself has relatively poor spatial resolution compared to fMRI, the proposed DOT inverse method based on MUSIC algorithm still shows its potential as an alternative algorithm that maps the brain activation.
In this work, the signal subspace dimension was determined as the number of singular values, so we focused on a few dominant singular values which were likely to be closely related with meaningful signal.While the first two or three singular values seem to be major with a large magnitude as shown in Fig. 2, there still does not exist a standardized criteria for the choice of singular values.Therefore, it is required to establish an objective guideline for the selection of signal subspace dimension, which will be further studied elsewhere.
One of the shortcomings in current study is that it requires simultaneous MRI acquisition to segment out the individual brains.However, an anatomical atlas-guided approach has recently been proposed to solve this problem [40], so in the future, we will incorporate such approach into our framework.

Fig. 2 .
Fig. 2. Singular value distribution: each singular value is included in either signal subspace or noise subspace.The number of singular values in signal subspace is equal to the dimension of signal subspace.

Fig. 3 .
Fig. 3. Example of MUSIC spectrum for HbO; dimension of signal subspace is (a) one, (b) three (c) five.The spectrum distribution depends on the signal subspace dimension.

Fig. 4 .
Fig. 4. (a) Simulation setup for synthetic data, (b) synthetic ΔC HbO (red) and ΔC HbR (blue) at the activation area, and (c) locations of synthetic activation area on the horizontal plane.The big arrow in (a) indicates the activation placed in the middle position of the most lateral line in (c).

Fig. 6 .
Fig. 6.Displacement error [mm] in lateral and depth directions when the activation changes along the lateral direction as shown in Fig. 4(c).The depth is fixed to (a) 30mm, (b) 44mm, respectively.(c) Lateral distance error averaged over depth.

Fig. 8 .
Fig. 8. Activation maps obtained from the right finger tapping task (Subject 1).MUSIC spectrum of DOT is based on combination of HbO and HbR, in which df = 1, q-level for FDR control is given as 0.01.fMRI images is controlled by p < 0.01 (corrected).(a) 3D rendering of the activation map for DOT (blue), fMRI (yellow-green), and intersectional area (red).(b) 2D maps for DOT with coronal, sagittal, and horizontal views.

Fig. 9 .
Fig. 9. Activation maps obtained from the right finger tapping task (Subject 3).MUSIC spectrum of DOT is based on combination of HbO and HbR, in which df = 1, q-level for FDR control is given as 0.01.fMRI images is controlled by p < 0.01 (corrected).(a) 3D rendering of the activation map for DOT (blue), fMRI (yellow-green), and intersectional area (red).(b) 2D maps for DOT with coronal, sagittal, and horizontal views.

Table 1 .
Voxel classification in multiple testing of N hypotheses.

Table 4 .
The number of significant voxels in fMRI or DOT (HbO, HbR, HbO&HbR) and the number of voxels which are located within BA4 or BA123 among the total significant voxels in each fMRI and DOT method.

Table 5 .
Displacement of the center of mass between fMRI and DOT using MUSIC spectra of (a) HbO, (b) HbR, and (c) HbO&HbR.

Table 6 .
Depth and lateral displacement of the center of mass between fMRI and DOT using MUSIC spectra of (a) HbO, (b) HbR, and (c) HbO&HbR.