Source localization approach for functional DOT using MUSIC and FDR control

Jin Wook Jung; Ok Kyun Lee; Jong Chul Ye

doi:10.1364/OE.20.006267

1. 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 deoxy-hemoglobin (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–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.

2. Theory

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

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

Δ ϕ (r_{d}, r_{s}; λ, t) = - ln \frac{U (r_{d}, r_{s}; λ, t)}{U_{0} (r_{d}, r_{s}; λ)} \approx \int \frac{U_{0} (r_{d}, r; λ) U_{0} (r, r_{s}; λ)}{U_{0} (r_{d}, r_{s}; λ)} Δ μ_{a} (r; λ, t) d r,

where Δμ_a(r;λ,t) is an absorption variation at r, λ denotes the illumination wavelength, and U₀(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}})}_{i = 1}^{m}

, Eq. (1) can be represented in a matrix form for a given wavelength λ at time t_j:

\begin{matrix} y_{j}^{λ} = A^{λ} u_{j}^{λ} + w_{j}^{λ}, & j = 1, 2, \dots, J, \end{matrix}

where

\begin{array}{l} y_{j}^{λ} = & {[Δ ϕ (r_{d_{1}}, r_{s_{1}}; λ, t_{j}), Δ ϕ (r_{d_{2}}, r_{s_{2}}; λ, t_{j}), \dots, Δ ϕ (r_{d_{m}}, r_{s_{m}}; λ, t_{j})]}^{T} \in ℝ^{m}, \\ u_{j}^{λ} = {[Δ μ_{a} (r_{1}; λ, t_{j}), Δ μ_{a} (r_{2}; λ, t_{j}), \dots, Δ μ_{a} (r_{n}; λ, t_{j})]}^{T} \in ℝ^{n}, \end{array}

where

{r_{i}}_{i = 1}^{n}

denotes voxel positions, and

w_{j}^{λ} \in ℝ^{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}^{λ}}_{i, j = 1}^{m, n}

is given by

A_{i j}^{λ} = \frac{U_{0} (r_{d_{i}}, r_{j}; λ) U_{0} (r_{j}, r_{s_{i}}; λ)}{U_{0} (r_{d_{i}}, r_{s_{i}}; λ)} δ,

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

Δ μ_{a} (r; λ, t) = ɛ_{H b O}^{λ} Δ c_{H b O} (r; t) + ɛ_{H b R}^{λ} Δ c_{H b R} (r; t),

where

ɛ_{H b O}^{λ} [μ M^{- 1} m m^{- 1}]

and

ɛ_{H b R}^{λ} [μ M^{- 1} m m^{- 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 {λ₁, λ₂}, the forward measurement model Eq. (2) can therefore be represented as

[\begin{matrix} y_{j}^{λ_{1}} \\ y_{j}^{λ_{2}} \end{matrix}] = A [\begin{matrix} x_{H b O, j} \\ x_{H b R, j} \end{matrix}] + [\begin{matrix} w_{j}^{λ_{1}} \\ w_{j}^{λ_{2}} \end{matrix}],

where the sensing matrix A is

A = [\begin{matrix} A_{H b O} & A_{H b R} \end{matrix}] = [\begin{matrix} ɛ_{H b O}^{λ_{1}} A^{λ_{1}} & ɛ_{H b R}^{λ_{1}} A^{λ_{1}} \\ ɛ_{H b O}^{λ_{2}} A^{λ_{2}} & ɛ_{H b R}^{λ_{2}} A^{λ_{2}} \end{matrix}],

and

x_{H b O, j} = {[Δ c_{H b O} (r_{1}; t_{j}), Δ c_{H b O} (r_{2}; t_{j}), \dots, Δ c_{H b O} (r_{n}; t_{j})]}^{T} \in ℝ^{n},

x_{H b R, j} = {[Δ c_{H b R} (r_{1}; t_{j}), Δ c_{H b R} (r_{2}; t_{j}), \dots, Δ c_{H b R} (r_{n}; t_{j})]}^{T} \in ℝ^{n} .

2.3. Source-localization formulation

By collecting all temporal series (j = 1,2,··· ,J), Eq. (5) can be formulated as

\begin{matrix} Y = [\begin{matrix} Y^{λ_{1}} \\ Y^{λ_{2}} \end{matrix}] = A [\begin{matrix} X_{H b O} \\ X_{H b R} \end{matrix}] + [\begin{matrix} W^{λ_{1}} \\ W^{λ_{2}} \end{matrix}] = A X + W, & A \in ℝ^{2 m \times 2 n}, & X \in ℝ^{2 n \times J} \end{matrix}

where

Y^{λ_{i}} = [y_{1}^{λ_{i}}, y_{2}^{λ_{i}}, \dots y_{j}^{λ_{i}}] \in ℝ^{m \times J}

, X_HbO = [x_HbO_,1, x_HbO_,2, ··· ,x_HbO_,_J] ∈ 𝕉ⁿ^×^J, X_HbR = [x_HbR_,1, x_HbR_,2, ··· ,x_HbR_,_J] ∈ 𝕉ⁿ^×^J, and

W^{λ_{i}} = [w_{1}^{λ_{i}}, w_{2}^{λ_{i}}, \dots, w_{J}^{λ_{i}}] \in ℝ^{m \times J}

.

From this formulation, Appendix A describes how this problem can be solved using the so-called MUSIC algorithm. Furthermore, under the Gaussian assumption of the sensing matrix, the following MUSIC spectrum can be assumed to have chi-squared distribution,

ν (j) = m {‖ U_{s}^{H} a_{j} ‖}_{2}^{2} .

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

ν (i) = {\begin{array}{l} m {‖ U_{s}^{H} a_{i} ‖}_{2}^{2}, & i = 1, \dots, n & : & ν_{H b O} \\ m {‖ U_{s}^{H} a_{i} ‖}_{2}^{2}, & i = n + 1, \dots, 2 n & : & ν_{H b R} \end{array}

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

\begin{matrix} ν_{ex} (i) = ν (i) + ν (n + i) = m {‖ U_{s}^{H} [a_{i} a_{n + i}] ‖}_{F}^{2}, & i = 1, \dots, n . & : & ν_{H b O & H b R} \end{matrix}

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.

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.

Download Full Size | PDF

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

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.

Download Full Size | PDF

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.

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.

Download Full Size | PDF

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

F D R = \frac{F P}{F P + T P} = \frac{F P}{T_{1}},

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

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

View Table | View all tables in this article

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:

E (FDR) \leq \frac{n_{0}}{n} q \leq 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₁ ≤ p₂ ≤ ... ≤ p_n which correspond to null hypotheses, H₁, H₂,..., H_n. The final step is to evaluate the following inequality in reverse order,

\begin{matrix} p_{i} \leq i \frac{q}{n}, & i = n, n - 1, \dots, 1. \end{matrix}

Let k be the largest i which satisfies the Eq. (15), then we reject all the hypotheses, H₁, H₂,..., 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.

3. Method

3.1. 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 sphere-shaped 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).

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

Download Full Size | PDF

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

Fig. 5 Block paradigm of RFT experiment.

Download Full Size | PDF

3.3. 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 [TR] = 3000ms, echo time [TE] = 35ms, a flip angle of 80°, 35 slices, 4mm slice thickness). Anatomical T1-weighted scans were successively taken to obtain the structural image of each subject based on fluid-attenuated inversion recovery (FLAIR) sequence ([TR] = 2800ms, [TE] = 16ms, a flip angle of 80°, 35 slices, 4mm slice thickness).

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

3.5. 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–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⁷ photons to calculate the Green’s function U₀(r, r_s;λ) and U₀(r_d, r;λ) (the reciprocity theorem [47] can be used to calculate U₀(r_d, r;λ)), and 10⁸ photons to calculate the optical flux measured at detectors U₀(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].

Table 2. The wavelength-dependent optical properties of the segmented head model (absorption coefficient μ_a and reduced scattering coefficient μ′_s) [41–46].

View Table | View all tables in this article

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

4. Results

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

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.

Download Full Size | PDF

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.

Download Full Size | PDF

4.2. 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 somato-sensory 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.

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.

Download Full Size | PDF

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.

Download Full Size | PDF

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.

Fig. 10 COM positions of fMRI, ν_HbO, ν_HbR, and ν_HbO&HbR.

Download Full Size | PDF

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

View Table | View all tables in this article

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

View Table | View all tables in this article

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 similar to the simulation results.

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

View Table | View all tables in this article

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.

View Table | View all tables in this article

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

Appendix A

Assuming that the activation areas are very sparse, the DOT activation detection problem can be formulated as

\begin{array}{l} min & {‖ X ‖}_{0} \\ subject to & Y = A X, \end{array}

where ||X||₀ denotes the number of non-zero rows, i.e. ||X||₀ = |suppX| = k, and where suppX refers to the indices of non-zero rows. This type of problem is often called a multiple measurement vector (MMV) problem, and is an extension of compressed sensing for joint sparse signals [25, 48].

Classically, the MMV problem Eq. (A.1), which was often termed as direction-of-arrival (DOA) or the bearing estimation problem, had been addressed using sensor array signal processing techniques [22]. One of the most popular and successful DOA estimation algorithms is the MUSIC (MUltiple SIgnal Classification) algorithm [27]. Under the assumption of zero-mean i.i.d. Gaussian noise, i.e. $w_{j}^{λ} \sim 𝒩 (0, σ^{2} I)$ in Eq. (9), the outer product of measurement matrix, $\frac{1}{J} Y Y^{T}$ , approaches the correlation matrix of Y given by [22, 27, 28]:

{\hat{R}}_{Y} = \frac{1}{J} Y Y^{T} \overset{J \to \infty}{\to} A {\hat{R}}_{X} A^{T} + σ^{2} I,

where

{\hat{R}}_{X} = \frac{1}{J} X X^{T}

. If k-incoherent sources exist, AR̂_XA^T can be decomposed as

A {\hat{R}}_{X} A^{T} = U_{S} \sum_{S} U_{S}^{T}

where U_s = [u₁,..., u_k] ∈ 𝕉²^m^×^k denotes a signal subspace, U_n ∈ 𝕉²^m^×2^m⁻^k denotes a noise subspace that orthogonal to span(U_s) and

\sum_{S} = [\begin{array}{l} σ_{1}^{2} \\ σ_{2}^{2} \\ ⋱ \\ σ_{k}^{2} \end{array}] \in ℝ^{k \times k} .

Under this decomposition, we can see that

span (U_{s}) = span (A X) = span (A_{I_{k}}),

where I_k denotes the true support.

The main observation in MUSIC is that the signal subspace should be highly correlated with the atom at the correct support, such that [22, 27, 28]

{‖ U_{S}^{H} a_{j} ‖}_{2}^{2} = 1,

or the noise subspace should be orthogonal to the correct support, such that

{‖ U_{n}^{H} a_{j} ‖}_{2}^{2} = 0.

Therefore, we can determine the correct support through either signal or noise subspace version of MUSIC algorithm.

Under the assumption that the elements of A follow Gaussian distribution N(0,1/m) and U_s is independent of a_j, scaled MUSIC spectrum of signal subspace version

ν (j) = m {‖ U_{S}^{H} a_{j} ‖}_{2}^{2},

follows the chi-squared distribution with degrees of freedom (df) equal to the dimension of the signal subspace [49].

Acknowledgments

This research was supported by the Korea Science and Engineering Foundation (KOSEF) grant funded by the Korea government (MEST) (No. 2011-0000855).

References and links

1. A. Villringer and B. Chance, “Non-invasive optical spectroscopy and imaging of human brain function,” Trends Neurosci. 20, 435–442 (1997). [CrossRef] [PubMed]

2. V. Ntziachristos and B. Chance, “Probing physiology and molecular function using optical imaging: applications to breast cancer,” Breast Cancer Res 3, 41–46 (2001). [CrossRef] [PubMed]

3. R. Weissleder and V. Ntziachristos, “Shedding light onto live molecular targets,” Nature Medicine 9, 123–128 (2003). [CrossRef] [PubMed]

4. D. A. Benaron, S. R. Hintz, A. Villringer, D. Boas, A. Kleinschmidt, J. Frahm, C. Hirth, H. Obrig, J. C. van Houten, E. L. Kermit, W. F. Cheong, and D. K. Stevenson, “Noninvasive functional imaging of human brain using light,” J. Cereb. Blood Flow Metab. 20, 469–477 (2000). [CrossRef] [PubMed]

5. S. Ogawa, D. W. Tank, R. Menon, J. M. Ellermann, S. Kim, H. Merkle, and K. Ugurbil, “Intrinsic signal changes accompanying sensory stimulation: Functional brain mapping with magnetic resonance imaging,” Proc. Nati. Acad. Sci. USA 89, 5951–5955 (1992). [CrossRef]

6. J. Steinbrink, A. Villringer, F. Kempf, D. Haux, S. Boden, and H. Obrig, “Illuminating the BOLD signal: Combined fMRI-fNIRS studies,” Magnetic Resonance Imaging 24, 495–505 (2006). [CrossRef] [PubMed]

7. A. P. Gibson, T. Austin, N. L. Everdell, M. Schweiger, S. R. Arridge, J. H. Meek, J. S. Wyatt, D. T. Delpy, and J. C. Hebden, “Three-dimensional whole-head optical tomography of passive motor evoked responses in the neonate,” NeuroImage 30, 521–528 (2006). [CrossRef]

8. A. Y. Bluestone, G. Abdoulaev, C. H. Schmitz, R. L. Barbour, and A. H. Hielscher, “Three-dimensional optical tomography of hemodynamics in the human head,” Opt. Express 9, 272–286 (2001). [CrossRef] [PubMed]

9. D. A. Boas, A. M. Dale, and M. A. Franceschini, “Diffuse optical imaging of brain activation: approaches to optimizing image sensitivity, resolution, and accuracy,” NeuroImage 23, S275–S288 (2004). [CrossRef] [PubMed]

10. H. Dehghani, B. R. White, B. W. Zeff, A. Tizzard, and J. P. Culver, “Depth sensitivity and image reconstruction analysis of dense imaging arrays for mapping brain function with diffuse optical tomography,” Appl. Opt. 48, D137–D143 (2009). [CrossRef] [PubMed]

11. N. M. Gregg, B. R. White, B. W. Zeff, A. J. Berger, and J. P. Culver, “Brain specificity of diffuse optical imaging: improvements from superficial signal regression and tomography,” Front Neuroenergetics 2, 1–8 (2010).

12. B. W. Zeff, B. R. White, H. Dehghani, B. L. Schlaggar, and J. P. Culver, “Retinotopic mapping of adult human visual cortex with high-density diffuse optical tomography,” Proc. Nat. Acad. Sci. USA 104, 12169–12174 (2007). [CrossRef] [PubMed]

13. S. P. Koch, C. Habermehl, J. Mehnert, C. H. Schmitz, S. Holtze, A. Villringer, J. Steinbrink, and H. Obrig, “High-resolution optical functional mapping of the human somatosensory cortex,” Front Neuroenergetics 2, 1–8 (2010).

14. B. R. White, A. Z. Snyder, A. L. Cohen, S. E. Petersen, M. E. Raichle, B. L. S, and J. P. Culver, “Resting-state functional connectivity in the human brain revealed with diffuse optical tomography,” NeuroImage 47, 148–156 (2009). [CrossRef] [PubMed]

15. H. Niu, Z. J. Lin, F. Tian, S. Dhamne, and H. Liu, “Comprehensive investigation of three-dimensional diffuse optical tomography with depth compensation algorithm,” J. Biomed. Opt. 15, 046005 (2010). [CrossRef] [PubMed]

16. F. Abdelnour, C. Genovese, and T. Huppert, “Hierarchical Bayesian regularization of reconstructions for diffuse optical tomography using multiple priors,” Biomed. Opt. Express 1, 1084–1103 (2010). [CrossRef]

17. M. Jacob, Y. Bresler, V. Toronov, X. Zhang, and A. Webb, “Level-set algorithm for the reconstruction of functional activation in near-infrared spectroscopic imaging,” J. Biomed. Opt. 11, 064029 (2006). [CrossRef]

18. D. A. Boas and A. M. Dale, “Simulation study of magnetic resonance imaging-guided cortically constrained diffuse optical tomography of human brain function,” Appl. Opt. 44, 1957–1968 (2005). [CrossRef] [PubMed]

19. V. Kolehmainen, S. Prince, S. R. Arridge, and J. P. Kaipio, “State-estimation approach to the nonstationary optical tomography problem,” J. Opt. Soc. Am. A 20, 876–889 (2003). [CrossRef]

20. S. Prince, V. Kolehmainen, J. P. Kaipio, M. A. Franceschini, D. Boas, and S. R. Arridge, “Time-series estimation of biological factors in optical diffusion tomography,” Phys. Med. Biol. 48, 1491–1504 (2003). [CrossRef] [PubMed]

21. Y. Zhang, D. H. Brooks, and D. A. Boas, “A haemodynamic response function model in spatio-temporal diffuse optical tomography,” Phys. Med. Biol. 50, 4625–4644 (2005). [CrossRef] [PubMed]

22. H. Krim and M. Viberg, “Two decades of array signal processing research,” IEEE Signal Proc. Mag. pp. 67–94 (1996). [CrossRef]

23. D. L. Donoho, “Compressed sensing,” IEEE Trans. Information Theory 52, 1289–1306 (2006). [CrossRef]

24. E. Candes, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory 52, 489–509 (2006). [CrossRef]

25. J. M. Kim, O. K. Lee, and J. C. Ye, “Compressive MUSIC: revisiting the link between compressive sensing and array signal processing,” IEEE Trans. Inf. Theory 58, 278–301 (2012). [CrossRef]

26. O. K. Lee, J. M. Kim, Y. Bresler, and J. C. Ye, “Compressive diffuse optical tomography: non-iterative exact reconstruction using joint sparsity,” IEEE Trans. Med. Imag. 30, 1129–1142 (2011). [CrossRef]

27. R. Schmidt, “Multiple emitter location and signal parameter estimation,” IEEE Trans. Antennas Propag. 34, 276–280 (1986). [CrossRef]

28. P. Stoica and A. Nehorai, “MUSIC, maximum likelihood, and Cramer-Rao bound,” IEEE Trans. Acoust., Speech, Signal Process. 37, 720–741 (1989). [CrossRef]

29. Y. Benjamini and Y. Hochberg, “Controlling the flase discovery rate: a practical and powerful approach to multiple testing,” J. R. Stat. Soc. Ser B (Methodl.) 57, 289–300 (1995).

30. J. Mosher, P. Lewis, and R. Leahy, “Multiple dipole modelling and localization from spatio-temporal MEG data,” IEEE Trans. Biomed. Eng. 39, 541–557 (1992). [CrossRef] [PubMed]

31. C. Genovese, N. Lazar, and T. Nichols, “Thresholding of statistical maps in functional neuroimaging using the false discovery rate,” NeuroImage 15, 870–878 (2002). [CrossRef] [PubMed]

32. A. Singh and I. Dan, “Exploring the false discovery rate in multichannel NIRS,” NeuroImage 33, 542–549 (2006). [CrossRef] [PubMed]

33. T. Li, H. Gong, and Q. Luo, “Visualization of light propagation in visible Chinese human head for functional near-infrared spectroscopy,” J. Biomed. Opt. 16, 045001 (2011). [CrossRef] [PubMed]

34. J. C. Ye, S. Tak, K. E. Jang, J. Jung, and J. Jang, “NIRS-SPM: Statistical parametric mapping for near-infrared spectroscopy,” NeuroImage 44, 428–447 (2009). [CrossRef]

35. T. Takahashi, Y. Takikawa, R. Kawagoe, S. Shibuya, T. Iwano, and S. Kitazawa, “Influence of skin blood flow on near-infrared spectroscopy signals measured on the forehead during a verbal fluency task,” NeuroImage 57, 991–1002 (2011). [CrossRef] [PubMed]

36. K. E. Jang, S. Tak, J. Jung, J. Jang, and J. C. Ye, “Wavelet minimum description length detrending for near-infrared spectroscopy,” J. Biomed. Opt. 14, 034004 (2009). [CrossRef] [PubMed]

37. K. J. Friston, O. Josephs, E. Zarahn, A. P. Holmes, S. Rouquette, and J.-B. Poline, “To smooth or not to smooth? bias and efficiency in fmri time-series analysis,” NeuroImage 12, 196–208 (2000). [CrossRef] [PubMed]

38. K. J. Friston, J. Ashburner, S. Kiebel, T. Nichols, and W. Penny, eds., Statistical parametric mapping: The analysis of functional brain images (Academic Press, Sandiego, CA, USA, 2006).

39. Q. Fang and D. A. Boas, “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express 17, 20178–20190 (2009). [CrossRef] [PubMed]

40. A. Custo, D. A. Boas, D. Tsuzuki, I. Dan, R. Mesquita, B. Fischl, W. Eric, L. Grimson, and W. Wells III, “Anatomical atlas-guided diffuse optical tomography of brain activation,” NeuroImage 49, 561–567 (2010). [CrossRef]

41. N. Okui and E. Okada, “Wavelength dependence of crosstalk in dual-wavelength measurement of oxy- and deoxy-hemoglogin,” J. Biomed. Opt. 10, 011015 (2005). [CrossRef]

42. R. M. P. Doornbos, R. Lang, M. C. Aalders, F. W. Cross, and H. J. C. M. Sterenborg, “The determination of in vivo human tissue optical properties and absolute chromophore concentrations using spatially resolved steady-state diffuse reflectance spectroscopy,” Phys. Med. Biol. pp. 967–981 (1999). [CrossRef] [PubMed]

43. E. Okada and D. T. Delpy, “Near-infrared light propagation in an adult head model. II. Effect of superficial tissue thickness on the sensitivity of the near-infrared spectroscopy signal,” Appl. Opt. 42, 2915–2922 (2003). [CrossRef] [PubMed]

44. S. A. Prahl, “Optical Properties Spectra,” Http://omlc.ogi.edu/spectra/index.html (Oregon Medical Laser Clinic), 2001.

45. W. Cheong, S. A. Prahl, and A. J. Welch, “A review of the optical properties of biological tissues,” IEEE J. Quantum Electron. 26, 2166–2185 (1990). [CrossRef]

46. C. R. Simpson, M. Kohl, M. Essenpreis, and M. Cope, “Near-infrared optical properties of ex vivo human skin and subcutaneous tissues measured using the Monte Carlo inversion technique,” Phys. Med. Biol. 43, 2465–2478 (1998). [CrossRef] [PubMed]

47. W. C. Chew, Waves and Fields in Inhomogeneous Media, 2nd ed. (IEEE Press, New York, 1995).

48. J. Chen and X. Huo, “Theoretical results on sparse representations of multiple measurement vectors,” IEEE Trans. Signal Process. 54, 4634–4643 (2006). [CrossRef]

49. A. Fletcher, S. Rangan, and V. Goyal, “Necessary and sufficient conditions for sparsity pattern recovery,” IEEE Trans. Inf. Theory 55, 5758–5772 (2009). [CrossRef]

	Declared non-significant	Declared significant	Total

Truly non-significant	TN	FP	n₀
Truly significant	FN	TP	n₁

Total declared	T₀	T₁	n

Wavelength [nm]	781		856

Optical parameters [mm⁻¹]	μ_a	μ′_s	μ_a	μ′_s

Superficial tissue (skin/skull)	0.018	1.83	0.0189	1.5556
CSF	0.0044	0.25	0.0056	0.2062
Gray matter	0.036	2.31	0.0428	2.0515
White matter	0.016	9.25	0.0175	8.6325

	unit [mm]	subject1		subject2		subject3		Group
		BA4	BA123	BA4	BA123	BA4	BA123	BA4	BA123
(a)	fMRI_COM(BOLD)	1.8711	1.1543	2.3348	5.1339	7.5629	3.5905	7.6072	4.9646
(b)	ν_COM(HbO)	14.7827	18.6245	10.2632	8.53	29.5933	18.2131	12.7561	5.6856
(c)	ν_COM(HbR)	14.1393	17.505	10.4962	9.2238	25.9842	17.0267	11.0966	7.6136
(d)	ν_COM(HbO&HbR)	14.3324	17.5451	10.7566	7.134	30.0659	18.5857	12.136	6.0928

		subject1			subject2			subject3			Group
	∩	Total	BA4	BA123	Total	BA4	BA123	Total	BA4	BA123	Total	BA4	BA123
(a)	fMRI_(BOLD)	3587	246	373	726	130	69	9315	327	672	17956	251	415
(b)	ν_(HbO)	8074	218	192	4737	200	328	6619	18	274	7809	185	319
(c)	ν_(HbR)	4464	99	88	3604	143	268	3463	0	99	5781	156	317
(d)	ν_(HbO&HbR)	10848	330	407	6877	256	406	9144	94	375	10898	208	337

	unit [mm]	subject 1	subject 2	subject 3	Group
(a)	ν_COM(HbO)	19.8966	28.0584	23.3319	14.7694
(b)	ν_COM(HbR)	19.9905	28.7928	24.0994	17.5162
(c)	ν_COM(HbO&HbR)	18.8452	27.247	23.3463	12.0618

Source localization approach for functional DOT using MUSIC and FDR control

Abstract

1. Introduction

2. Theory

2.1. Notation

2.2. Forward model

2.3. Source-localization formulation

2.4. Signal subspace selection

2.5. Inference using MUSIC spectrum

3. Method

3.1. Simulation setup

3.2. Task protocol for in vivo experiment

3.3. Real data acquisition

3.4. Preprocessing

3.5. Computation of sensing matrix

3.6. Quantitative analysis

4. Results

4.1. Synthetic data

4.2. In vivo real data

5. Conclusion and discussion

Appendix A

Acknowledgments

References and links

Cited By

Figures (10)

Tables (6)

Equations (24)

Optics Express

	unit [mm]	subject1		subject2		subject3		Group
		depth	lateral	depth	lateral	depth	lateral	depth	lateral
(a)	ν_COM(HbO)	4.6477	19.3461	11.684	25.5099	8.2951	21.8075	11.3086	9.5001
(b)	ν_COM(HbR)	1.3603	19.9442	11.9276	26.2061	12.6373	20.5203	13.0308	11.7052
(c)	ν_COM(HbO&HbR)	4.8906	18.1995	11.152	24.8602	6.2693	22.4888	8.4128	8.6436