Extended hierarchical Bayesian diffuse optical tomography for removing scalp artifact

Takeaki Shimokawa; Takashi Kosaka; Okito Yamashita; Nobuo Hiroe; Takashi Amita; Yoshihiro Inoue; Masa-aki Sato

doi:10.1364/BOE.4.002411

1. Introduction

Functional near-infrared spectroscopy (fNIRS) is a technique that non-invasively measures hemodynamic responses to neuronal activation in the cerebral cortex [1–6]. To observe these responses, near-infrared light is emitted into the head from source probes on the scalp, and the light that passes through the head is measured by detector probes on the scalp. Since fNIRS measurement provides such advantages as safety, portability, low cost as well as few physical restrictions, it is expected to extend the application studies of neuroscience as an easy-to-use neuro-monitoring tool.

However, a major problem with fNIRS measurement is the artifact originating from the hemodynamic change in the scalp. The sensitivity of fNIRS measurement is very high near the observation surface, namely the scalp, and it exponentially decreases with increased distance from the scalp. Therefore, a small hemodynamic fluctuation in the scalp results in a large observation signal. Recently, it has been reported that scalp hemodynamic responses show task-related changes, suggesting that the scalp artifact may be misinterpreted as cortical activities in fNIRS measurements [7, 8]. Many methods have been developed for removing the scalp signals. It is empirically known that the hemodynamic change in the superficial layers, including the scalp and the skull, is not localized but instead spatially global [9–13]. To remove or isolate this global artifact, most methods have used temporal information, and often the temporal values of the short-distance measurement channel as a reference of superficial hemodynamic fluctuations for regression [9–15]. However, the inhomogeneity of the temporal patterns in the hemodynamic change of the superficial layers [16, 17] or the temporal correlation between the hemodynamic responses in the scalp and in the cortex [7, 8, 18] negatively impacts such artifact-removal methods based on temporal information. In this study, by contrast, the spatial information of the optical paths of all observation channels, which are not affected by temporal inhomogeneity or correlation, is fully exploited for removing the scalp signals.

We propose a method that separates the signals from the cortex and scalp by estimating the hemodynamic changes in both regions by diffuse optical tomography (DOT), which is an advanced imaging method designed to increase the spatial resolution of fNIRS measurements. Continuous-wave DOT reconstructs the distribution of the changes in optical properties inside the tissue from the changes in measured light intensity using knowledge of the spatial information of the optical paths of measurement channels [6, 19–22]. In our previous study, we proposed a three-dimensional reconstruction method with high spatial resolution in both the depth and horizontal directions [23], where we introduced sparse regularization to improve the depth accuracy and the spatial resolution. In this study, we extend this method and reconstruct the distributions of hemoglobin concentration change in both the cortex and the scalp. Although the introduction of sparse regularization is reasonable for accurate localization of cortical activity, it is not reasonable for estimating the hemodynamic change in the scalp because it is rather global. Therefore, we introduced different types of regularization for the cortex and the scalp. Our proposed method applies sparse regularization to the cortical hemodynamic change and smooth regularization to the scalp’s hemodynamic change. The parameter controlling the smoothness of the scalp regularization is adjusted from the data in the estimation process. Consequently, inhomogeneous scalp hemodynamic change distributions can be appropriately estimated. In this way, our proposed method simultaneously estimates the hemodynamic changes in both the cortex and the scalp.

We conducted two experiments to validate our proposed method. First, we applied it to the data obtained from a two-layer phantom experiment measured with a commercially available near-infrared imaging system. The global absorption change in the scalp and the local absorption change in the cortex were mimicked by a two-layer phantom. Second, we evaluated it by a simulation study using an MRI-based five-layer head model [20]. In both experiments, the proposed method separated the signals from the cortical area and from the scalp area and accurately estimated the three-dimensional position of the local activities in the cortical area.

The following notations are used throughout this paper. Scalars are represented in normal-type letters as x. Vectors are represented in lower-case boldface letters as x. Matrices are represented in upper-case boldface letters as X. The transpose, trace, and determinant of matrix X are denoted as X^T, tr(X), and |X|, respectively. A diagonal matrix with diagonal elements x is denoted as diag(x). A vector consisting of the diagonal elements of matrix X is denoted as diag(X). The L₂ norm of vector x is denoted as ${‖ x ‖}_{2}^{2}$ . The expectation of X with respect to a probability P(X) is denoted as 〈X〉_P₍_X₎. The mean value $\frac{1}{T} Σ_{t = 1}^{T} x_{t}$ is denoted as mean(x). The N-dimensional vector consisting of 1s and the identity matrix of size N are denoted as 1_N and I_N, respectively. N(x; μ, Σ) denotes that a vector random variable x follows a multivariate Gaussian distribution with mean vector μ and covariance matrix Σ. Γ(x; μ, γ) denotes that a scalar random variable x follows a gamma distribution with mean μ and shape parameter γ (note that the scale parameter corresponds to μγ⁻¹).

2. Methods

Here, we introduce our method to reconstruct the three-dimensional distributions of the concentration changes of oxygenated and deoxygenated hemoglobin, ΔHbO₂ and ΔHbR, in both the cortical and scalp areas. First, we reconstructed the three-dimensional distributions of the absorption changes at multiple wavelengths. Then we converted these absorption change distributions to the ΔHbO₂ and ΔHbR distributions.

We used the Bayesian approach to reconstruct the distributions of the absorption changes. The forward model of DOT, which relates the absorption changes inside the biological tissue to the observed light intensity change, is formulated as a conditional probability function. The regularization used to resolve the ill-posed nature of the inverse problem of DOT is formulated as a prior probability function. We obtained the estimate by inverting the conditional probability function with Bayes’ theorem.

Our proposed Bayesian estimation method is an extended version of a previous work [23], where we applied the automatic relevance determination (ARD) sparse prior to regularize the cortical activities. Here, we add the Laplacian smooth prior to regularize the scalp activities.

2.1. Forward model

Assuming that the absorption changes are small, we use the Rytov approximation [24–27], which linearly relates the absorption changes and the measured light intensity changes. Here, the target area of the absorption changes is separated into cortical and scalp areas:

y = A^{c} x^{c} + A^{s} x^{s} + ξ,

where

x^{c} = {(x_{1}^{c}, \dots, x_{N^{c}}^{c})}^{T}

and

x^{s} = {(x_{1}^{s}, \dots, x_{N^{s}}^{s})}^{T}

represent absorption changes Δμ_a in the voxels of the cortical and scalp areas, y = (y₁, ···,y_M)^T represents the measured light intensity changes, A^c and A^s represent the sensitivity matrix of the cortical and scalp areas, and ξ represents the measurement noise. N^c and N^s denote the number of voxels in the cortical and scalp areas, and M denotes the number of measurement logical channels.

In a measurement, we observe multiple wavelength data and a number of time sample data. We summarize them in a matrix form. Observed light intensity changes y_k,t of wavelength k (k = 1,···,K) and time sample t (t = 1,···,T) are summarized as Y = (Y₁, ···, Y_k, ···, Y_K), where Y_k = (y_k_,1, ···, y_k,t, ···, y_k,T). Likewise, the absorption changes are summarized as $X^{r} = (X_{1}^{r}, \dots, X_{k}^{r}, \dots, X_{K}^{r})$ , where $X_{k}^{r} = (x_{k, 1}^{r}, \dots, x_{k, t}^{r}, \dots, x_{k, T}^{r})$ and r ∈ [c, s].

We constructed a probabilistic forward model assuming that ξ is a Gaussian noise. This assumption is validated by the phantom experimental data. We also assume that the noise covariance matrix in the absorption change measurement (task condition) is proportional to the noise covariance matrix under baseline condition (rest condition) Σ_y. Then the conditional probability for observed light-intensity change Y given absorption change X^c, X^s is given by

P (Y | X^{c}, X^{s}, σ) = \prod_{t = 1}^{T} \prod_{k = 1}^{K} N (y_{k, t}; A_{k}^{c} x_{k, t}^{c} + A_{k}^{s} x_{k, t}^{s}, σ^{- 1} Σ_{y, k}),

where σ is the scale parameter inversely proportional to the measurement noise variance. The scale parameter is estimated by introducing a hierarchical prior:

P (σ) \propto σ^{- 1} .

2.2. Regularization model

We used different types of regularization for the cortical and scalp areas. We applied the ARD prior and localized smoothing filter to the cortical area because the cortical activities are localized. In contrast, we applied a smoothing Laplacian prior to the scalp area because the hemo-dynamic responses in the scalp are rather global.

For the cortical area, we assume that solution x^c becomes a mixture of a small number of focal activities. Therefore, we represent the solution as x^c = Wz^c, where W is a spatially localized smoothing filter. The ARD prior is adopted for z^c. The probability distribution of z^c is given by a Gaussian distribution with a zero mean vector and a diagonal covariance matrix whose diagonal elements are adjustable hyperparameters:

P (Z^{c} | λ) = \prod_{t = 1}^{T} \prod_{k = 1}^{K} N (z_{k, t}^{c}; 0, Λ^{- 1}),

where Λ = diag(λ) and λ = (λ₁, ···,λ_N^c)^T represent the hyperparameters that are automatically adjusted from the data. Here, we use the common λ value for all of the wavelength data. Equation (4) can be equivalently rewritten as

P (X^{c} | λ) = \prod_{t = 1}^{T} \prod_{k = 1}^{K} N (x_{k, t}^{c}; 0, W Λ^{- 1} W^{T}) .

Hyperparameter λ is estimated by introducing a hierarchical prior:

P (λ) = \prod_{i = 1}^{N^{c}} Γ (λ_{i}; {\bar{λ}}_{0 i}, γ_{0 i}),

where λ̄_0i and γ_0i represent the mean and reliability (or confidence) of the hierarchical prior, respectively [28, 29].

For the scalp area, we assume solution x^s is spatially smooth. That is, the Laplacian (or 2nd derivative) of solution x^s, Lx^s, becomes small. Therefore, regularization is applied for z^s = Lx^s as

P (Z^{s} | η) = \prod_{t = 1}^{T} \prod_{k = 1}^{K} N (z_{k, t}^{s}; 0, {(η I_{N^{s}})}^{- 1}),

where η is the hyperparameter controlling the degree of smoothness. Equation (7) can be equivalently rewritten as

P (X^{s} | η) = \prod_{t = 1}^{T} \prod_{k = 1}^{K} N (x_{k, t}^{s}; 0, L^{- 1} η^{- 1} L^{- 1^{T}}) .

Hyperparameter η is estimated by introducing a hierarchical prior:

P (η) = Γ (η; {\bar{η}}_{0}, γ_{η_{0}}),

where η̄₀ and γ_η₀ represent the mean and reliability (or confidence) of the hierarchical prior.

2.3. Variational Bayesian method

The posterior probability of unknown X^c, X^s, λ, η, σ given observation data Y is obtained by Bayes’ theorem,

P (X^{c}, X^{s}, λ, η, σ | Y) = \frac{P (Y, X^{c}, X^{s}, λ, η, σ)}{P (Y)},

where P(Y, X^c, X^s, λ, η, σ) = P(Y | X^c, X^s, σ)P(X^c|λ)P(X^s|η)P(λ)P(η)P(σ) and P(Y) = ∫ dXX^cdXX^sdλdηdσP(Y, X^c, X^s, λ, η, σ).

Because analytically calculating the posterior is difficult, we converted this posterior calculation problem to a maximization problem of a free energy function. The free energy for trial distribution Q(X^c, X^s, λ, η, σ), which is the approximation of posterior P(X^c, X^s, λ, η, σ|Y), is defined by

\begin{array}{l} F (Q) & = & \int d X^{c} d X^{s} d λ d η d σ Q (X^{c}, X^{s}, λ, η, σ) log \frac{P (Y, X^{c}, X^{s}, λ, η, σ)}{Q (X^{c}, X^{s}, λ, η, σ)}, \\ = & log P (Y) - KL [Q (X^{c}, X^{s}, λ, η, σ) | | P (X^{c}, X^{s}, λ, η, σ | Y)], \end{array}

where

KL [P (θ) | | Q (θ)] = \int d θ P (θ) log \frac{P (θ)}{Q (θ)}

denotes the Kullback-Leibler distance. Equation (11) implies that the maximization of free energy F(Q) is equivalent to the minimization of the Kullback-Leibler distance between the posterior and the trial distribution. The maximization problem was iteratively solved using a factorization approximation of the trial distribution Q(X^c, X^s, λ, η, σ) (see Appendix). We used the solution of the sensitivity-normalized Tikhonov regularization as an initial value of the iteration algorithm (see Appendix). Finally, we obtained the approximated mean value with respect to the posterior, X̂^c = 〈X^c〉_{Q(X^c,X^s,λ,η,σ)} ≃ 〈X^c〉_{P(X^c,X^s,λ,η,σ|Y)} and X̂^s = 〈X^s〉_{Q(X^c,X^s,λ,η,σ)} ≃ 〈X^s〉_{P(X^c,X^s,λ,η,σ|Y)}.

2.4. Conversion from absorption changes to hemoglobin concentration changes

Estimated absorption changes x̂ = Δμ_a in multiple wavelengths are converted to the concentration changes of the oxygenated and deoxygenated hemoglobin (ΔHbO₂ and ΔHbR) in each voxel. These values are connected by the spectral extinction coefficients:

(\begin{matrix} Δ μ_{a, 1} \\ ⋮ \\ Δ μ_{a, K} \end{matrix}) = (\begin{matrix} ε_{1}^{{HbO}_{2}} & ε_{1}^{HbR} \\ ⋮ & ⋮ \\ ε_{K}^{{HbO}_{2}} & ε_{K}^{HbR} \end{matrix}) (\begin{matrix} Δ {HbO}_{2} \\ Δ HbR \end{matrix}) \equiv E (\begin{matrix} Δ {HbO}_{2} \\ Δ HbR \end{matrix}),

where

ε_{k}^{{HbO}_{2}}

and

ε_{k}^{HbR}

denote the extinction coefficients of the oxygenated and deoxygenated hemoglobin of the k-th wavelength, respectively [30]. The estimated values of ΔHbO₂ and ΔHbR are derived using the least-squares method [10]:

(\begin{matrix} Δ {HbO}_{2} \\ Δ HbR \end{matrix}) = {(E^{T} E)}^{- 1} E^{T} (\begin{matrix} Δ μ_{a, 1} \\ ⋮ \\ Δ μ_{a, K} \end{matrix}) .

3. Validation of method by two-layer phantom experiments

In this section, we describe how we validated our proposed method using a two-layer phantom experiment. Here, only one-wavelength data were used, and the positional accuracy of the estimated absorption change was evaluated. Because this experiment was static, the temporal information was not useful and the spatial information was only useful for separating the artifact signals.

3.1. Experimental conditions

We used a two-layer phantom consisting of liquid and silicone layers. The liquid layer was a substitute for the cortex, and the silicone layer was a substitute for the scalp. The measurement was done under two conditions: rest and task. We used the light intensity changes between these two conditions for our analysis. In the rest condition, a tank was filled with the liquid, and a silicone plate was placed on the bottom of the tank. In the task condition, we replaced the silicone plate with another one whose absorption coefficient was slightly higher than that of the rest condition to mimic the global activity in a scalp layer. We also added absorbers to mimic the local cortical activity (Fig. 1(a)). These black spherical absorbers had 5-mm diameters made of polyacetal. The liquid was made of 1.4% intralipid and ink. The silicone plates were 5-mm thick and made of silicone, kaolin powder, and ink. The kaolin powder and ink were stirred uniformly into a silicone liquid and cured to make a uniform silicone plate. Note that uncontrollable inhomogeneity remains in the silicone plate. Intralipid and kaoline powder were used to regulate the scattering coefficients. The ink regulated the absorption coefficients. The estimated absorption and reduced scattering coefficients (μ_a, μ′_s) of the liquid, the silicone plate for the rest condition, and the silicone plate for the task condition were (0.018 mm⁻¹, 1.1 mm⁻¹), (0.015 mm⁻¹, 0.9 mm⁻¹), and (0.016 mm⁻¹, 0.9 mm⁻¹), respectively. These estimations were based on measurements with a spectral photometer and a near-infrared imaging system.

Fig. 1 Setup of two-layer phantom experiment. (a) Schematic of phantom. The liquid filling the tank represents the cortical area. Silicone plate on the tank’s bottom represents the scalp area. Measurements were done in rest and task conditions. In the task condition, spherical absorbers were inserted in the liquid. The absorption coefficient of the silicone plate in the task condition was slightly higher than that in the rest condition. (b) Arrangement of source and detector probes on bottom of tank. Red circles represent the source, and blue squares represent the detector. x-y positions are superposed where the spherical absorbers were settled in measurements. One absorber was settled at (i) (0, 0, d) in the one-absorber experiment. Two absorbers were settled at (ii) ( $- \frac{l}{2}$ , 0, d) and ( $- \frac{l}{2} + h$ , 0, d) in the two-absorber experiment. Here, d denotes the depth or the vertical distance between the observation surface and the center of the absorber, h denotes the horizontal distance between the centers of the two absorbers, and l denotes the probe interval.

Download Full Size | PDF

We used a continuous-wave near-infrared imaging system (FOIRE-3000, Shimadzu Corporation) to acquire data to reconstruct the absorption change. To achieve a sufficient dynamic range for DOT, we customized the system’s amplifier to reach 256 times gain at maximum, which is 4 times higher than the commercially available FOIRE-3000. Here, we only used 805-nm wavelength data. Source and detector probes were alternately placed on the tank’s bottom in a rectangular arrangement (Fig. 1(b)). The sampling rate was 130 ms. Each time, we measured for 20 seconds in the rest and task conditions and tested the estimation results under two probe-interval conditions: l = 13 mm and 18.4 mm ( $= 13 mm \times \sqrt{2}$ ).

We tested the proposed method by one- and two-absorber experiments. First, we evaluated the depth accuracy using one spherical absorber placed at the x-y center position (i) (0, 0, d) at various depths (Fig. 1(b)). The results are summarized in section 3.3. Second, we evaluated the spatial resolution using two spherical absorbers placed in (ii) ( $- \frac{l}{2}$ , 0, d) and ( $- \frac{l}{2} + h$ , 0, d) at various depths d and various horizontal distances h (Fig. 1(b)). The results are summarized in section 3.4. One example of this experiment is shown in section 3.2 with a three-dimensional image to understand the proposed method.

The voxel size was 2.5 mm × 2.5 mm × 2.5 mm. The numbers of voxels in the cortical and scalp areas were N^c = 6250 = 25 × 25 × 10 and N^s = 1250 = 25 × 25 × 2, respectively. The number of measurement logical channels M was 56 (l = 13 mm, up to 3rd nearest-neighbor (NN) channels) or 48 (l = 18.4 mm, up to 2nd NN channels). We used a 5-mm FWHM Gaussian filter as W. Control constant β of sensitivity-normalized Tikhonov regularization (see Appendix) was set with reference to the maximum value of the sensitivity for a region deeper than 20 mm from the surface: β = max(ρ_i), i ∈ I_{≥20 mm}, where I_{≥20 mm} denotes the set of voxel indices deeper than 20 mm.

3.2. Estimation example

In this section, we show the advantage of the proposed method compared to previous methods by applying it to one example of a 3D reconstruction in the two-layer phantom experiment. We obtained the example data by the phantom experiment with an arrangement of 18.4-mm probe intervals. Here, 5-mm-thick silicone plates were used to mimic the global scalp activity, and two spherical absorbers were used to mimic the local cortical activities. Center positions of the absorbers were (−9.2, 0, 15) [mm] and (5.8, 0, 15) [mm] (Fig. 2(a)). The horizontal distance between two absorbers h was 15 mm, which was shorter than the probe intervals of 18.4 mm. Note that the coordinate system of Fig. 2 is flipped upside down to that of Fig. 1 so that it resembles the functional neuroimaging in which the probes are usually placed on top of the surface.

Fig. 2 Example of 3D reconstruction image. Top 5-mm-thick layer is the silicone layer representing the scalp, and the remaining region deeper than 5 mm is the liquid layer representing the cortex. Three-dimensional image of the absorption change (tomography) is represented in gray scale, and the x-y maximum intensity projection (topography) of the absorption change in the liquid layer is represented in green scale on the z = 30 surface. (a) True distribution of absorption change. Positions of centers of two absorbers are (−9.2, 0, 15) [mm] and (5.8, 0, 15) [mm]. Red cylinders represent the source, and blue square prisms represent the detector. (b) Solution of sensitivity-normalized Tikhonov regularization. (c) Solution of hierarchical Bayesian estimation method when only using the sparse regularization model. (d) Solution of hierarchical Bayesian estimation method using both the sparse and smooth regularization models.

Download Full Size | PDF

Figure 2(b) shows the estimation results of the sensitivity-normalized Tikhonov regularization. Because Tikhonov regularization cannot localize the solution, the two local absorption changes in the liquid layer and the global absorption change in the silicone layer are mixed. Values deeper than 20 mm are suppressed by the introduction of control constant β (see Appendix). The following hierarchical Bayesian estimations used this solution as an initial value of their iteration algorithm.

Figure 2(c) shows the estimation results of the hierarchical Bayesian estimation where the ARD sparse prior was applied to both the liquid and silicone areas. The sparse prior mismatches the global absorption changes in the silicone area. As a result, the local absorption changes in the liquid layer and the global absorption change in the silicone layer could not be completely distinguished. Accordingly, the estimated position of the two absorption changes in the liquid area is different from the true position.

Figure 2(d) shows the estimation result of the hierarchical Bayesian estimation where the ARD sparse prior is applied to the liquid area and the Laplacian smooth prior is applied to the silicone area. The two local absorption changes in the liquid layer and the global absorption change in the silicone layer can be accurately distinguished. The positional error of the two local absorption changes in the liquid area is smaller than the voxel size of 2.5 mm. This result is summarized again in section 3.4 (d = 15 mm, h = 15 mm).

3.3. One-absorber experiment

In this section, we examine whether the depth of a local absorption change can be accurately estimated by our proposed method in spite of the presence of a superficial global absorption change. We tested the proposed method using a one-absorber setting. The spherical absorber mimicking the local cortical activity was placed in the center position (i) (Fig. 1(b)) at various depths: d = 10 mm, 12.5 mm, 15 mm, 17.5 mm, and 20 mm. The silicone plates were also used to mimic the global scalp activity. The analysis was done under two conditions of probe intervals: l = 18.4 mm and 13 mm.

Figure 3(a) shows the estimation results for various depths under the two probe-interval conditions. Here, the 3D stimation results are represented as the 2D x-z projection images to clearly visualize the estimation accuracy. In the figure, the maximum intensity of the three central voxel regions in the y direction (7.5-mm-thick layer in Fig. 3(b)), including the true position of the 5-mm-diameter spherical absorber, was projected onto the x-z plane to check the positional accuracy in the y direction. If a positional error occurs in the y direction, the estimation value in this region becomes nearly zero. The white circles in the projection images in Fig. 3(a) represent the true position of the spherical absorber. We judged the estimation results as successful when the positional error between the center position of the true spherical absorption change and the voxel center of the maximum estimation value in the liquid layer was within one voxel (2.5 mm) in each x, y, z direction and when the maximum estimation value was greater than 0.025 mm⁻¹. The estimation image is surrounded by a bold frame if the estimation result was judged as successful. Figure 3(c) represents the depth of the maximum estimation value in the liquid layer exceeding 0.025 mm⁻¹ as a function of the true depth of the spherical absorption change.

Fig. 3 Estimation results of one-absorber experiment. (a) Estimation results under various depths and probe-interval conditions. From left, results under depth conditions of the spherical absorption change: d = 10 mm, 12.5 mm, 15 mm, 17.5 mm, and 20 mm. From top, results under two probe-interval conditions: l = 18.4 mm, 13 mm. White circles represent the true position of the 5-mm-diameter spherical absorption change in the liquid layer. Estimation result is surrounded by a bold frame if the positional error of the spherical absorption change was within one voxel (2.5 mm) in each x, y, z direction and the maximum estimation value exceeded 0.025 mm⁻¹. (b) Projection region is the central 7.5-mm-thick layer containing the true distribution of the spherical absorption change. Estimation images of (a) were obtained by projecting the maximum estimation value within the projection region in the y direction onto the x-z space. (c) Estimated depth of spherical absorption changes as a function of the true depth.

Download Full Size | PDF

The first and second rows from the top of Fig. 3(a) show the estimation results for the 18.4-and 13-mm probe-interval conditions, respectively. The spherical absorption changes were accurately estimated from d = 10 to 17.5 mm in the 18.4-mm probe-interval condition and from d = 10 to 20 mm in the 13-mm probe-interval condition, respectively, where the positional errors were within one voxel (see also Fig. 3(c)). The amplitudes of the estimated absorption changes were maintained at a roughly constant level. No absorption change deeper than d = 20 and d = 22.5 mm could be detected in the 13- and 18.4-mm probe-interval conditions, respectively. These results show that the depth of the local absorption change can be accurately estimated in spite of the presence of the superficial global absorption change if the local absorption change is within the depth reached by an adequate amount of diffuse light.

3.4. Two-absorber experiment

In this section, we evaluate the accuracy of the spatial resolution with respect to the local absorption changes in the presence of a superficial global absorption change. We tested the proposed estimation method using a two-absorber setting. In this phantom experiment, two spherical absorbers were placed in (ii) ( $- \frac{l}{2}$ , 0, d) and ( $- \frac{l}{2} + h$ , 0, d) (Fig. 1(b)) to mimic the local cortical activities, where h denotes the horizontal distance, d denotes the depth of the two absorbers, and l denotes the probe interval. Silicone plates were also used to mimic the global scalp activity. We evaluated the estimation results under two probe-interval conditions: 18.4 and 13 mm.

Figure 4 shows the estimation results for the 18.4-mm probe-interval condition. The estimations were tested under four conditions of horizontal distance, h = 17.5, 15, 12.5, and 10 mm, and under three depth conditions, d = 12.5, 15, and 17.5 mm. The 3D estimation results are represented as the 2D x-z intensity projection images, as in Fig. 3(a). The white circles represent the true positions of the two absorbers. We judged the estimation results to be successful if the positional errors of both spherical absorbers were within one voxel (2.5 mm) in each x, y, z direction and both maximum estimation values exceeded 0.025 mm⁻¹. The estimation image is surrounded by a bold frame if the estimation result was judged successful.

Fig. 4 Estimation results for 18.4-mm probe-interval condition. From left, results under four conditions of horizontal distance: h = 17.5, 15, 12.5, and 10 mm. From top, results under three depth conditions: d = 12.5, 15, and 17.5 mm. White circles represent the true positions of the spherical absorption change. Estimation image is surrounded by a bold frame if the positional errors of both spherical absorbers were within one voxel (2.5 mm) in each x, y, z direction and both maximum estimation values exceeded 0.025 mm⁻¹.

Download Full Size | PDF

For horizontal distance h = 17.5 mm, which is almost the same as probe interval l = 18.4 mm, the estimation was accurate up to the depth condition of d = 17.5 mm (Fig. 4) (h = 17.5 mm, d = 17.5 mm). For horizontal distance h = 12.5 mm, which is much shorter than the probe interval, the estimation was difficult but successful up to a depth condition of d = 12.5 mm (Fig. 4) (h = 12.5 mm, d = 12.5 mm).

The same analysis was also done under a 13-mm probe-interval condition. For horizontal distances h = 17.5, 15, 12.5, and 10 mm, the estimation was accurate up to depth conditions of d = 17.5, 15, 12.5, and 12.5 mm, respectively. These results show that even if the observation signal is contaminated by the superficial global absorption change, it is still possible to discriminate the two local absorption changes with distances shorter than the probe interval to a certain depth condition.

4. Validation of method by MRI-based head-model simulations

In this section, we explain how we validated our proposed method by a simulation study using an MRI-based head model with a much more complex shape than that used in the two-layer phantom experiment. Time-varying hemoglobin concentration changes were estimated using multiple wavelength data. The scalp hemodynamic change measured in a real experiment with human subjects was used in this simulation study. Note that the results investigated in section 3, such as the depth limits up to which the estimations were successful or the accuracy of the spatial resolution, do not necessarily correspond to those of this simulation study because the optical parameters of the media and the strengths of the observation noise and absorption changes were different.

4.1. Simulation conditions

The head model was constructed based on an anatomical MRI image (Trio 3T, Siemens, Ger-many) of an experimental subject (YO). T1-weighted MPRAGE (TE = 3.06 ms, TR = 2250 ms, flip angle = 9 deg, voxel size = 1 mm×1 mm×1 mm) was segmented to a five-layer head model by FreeSurfer [31] and VBMEG [32], which consists of the scalp, the skull, cerebrospinal fluid (CSF), and the gray and white matter.

We made a DOT measurement of the left motor area during a right-finger-tapping task and put 5 × 5 probes on the left hemisphere to cover the primary motor area (Fig. 5). The probe interval was l = 13 mm. The sensitivity was calculated by MCX [33] based on the head model and the probe arrangement. We set the optical parameters summarized in the paper by Fang [34] to each layer of the head model: The absorption coefficient, scattering coefficient, anisotropy, and refractive index (μ_a, μ_s, g, n) of the scalp, the skull, CSF, and the gray and white matter were (0.019 mm⁻¹, 7.8 mm⁻¹, 0.89, 1.37), (0.019 mm⁻¹, 7.8 mm⁻¹, 0.89, 1.37), (0.004 mm⁻¹, 0.009 mm⁻¹, 0.89, 1.37), (0.02 mm⁻¹, 9.0 mm⁻¹, 0.89, 1.37), and (0.08 mm⁻¹, 40.9 mm⁻¹, 0.84, 1.37), respectively. The finger-tapping task continued from 0 to 15 sec and the measurement time was from −5 to 35 sec.

Fig. 5 Arrangement of source and detector probes. Red points (e.g., labeled ‘T1’) represent the source probe, and the blue points (e.g., labeled ‘R1’) represent the detector probe. Region considered for DOT reconstruction image is illustrated as a cuboid.

Download Full Size | PDF

To reduce the amount of calculation, we restricted the area used for our DOT simulation study to an 80 mm × 80 mm × 60 mm region (Fig. 5). For the same reason, 1 mm × 1 mm × 1 mm voxels were merged to 2 mm × 2 mm × 2 mm voxels.

We tested our proposed method by putting one localized hemoglobin concentration change in the gray matter area and the global hemoglobin concentration change in the scalp area. The localized hemoglobin concentration change was placed at 84 various positions near the 5 × 5 probes: every second voxel of the area within a 35-mm depth from the head surface and within the square pole made by the four normal directions of R4, R5, R9, and R10 to the head surface. First, the total amount of changes in hemoglobin concentration of ΔHbO₂ and ΔHbR was placed in one of these voxels. Then it was smoothed by an 8-mm FWHM Gaussian smoothing filter, W, whose value was set to zero outside the gray matter area. The integral of W was normalized to 1. This W was also used in the estimation process. The peak values of ΔHbO₂ and ΔHbR were set to +6μM and −2μM, respectively, on the condition that the Gaussian distribution was all inside the gray matter area. The actual peak values depended on the shape of the gray matter area. The actual mean values of ΔHbO₂ and ΔHbR for the 84 peaks were +9.71μM and −3.24μM, respectively. The localized hemoglobin concentration changes of ΔHbO₂ and ΔHbR vary in time along with the hemodynamic response function; a gamma function of a = 6 and b = 0.9 sec is convoluted with the finger-tapping task period from 0 to 15 sec [35,36]. The global hemoglobin concentration changes in the scalp area were generated based on the measured data of a human finger-tapping task with the same probe arrangement. For this generation process, we randomly chose one trial of the finger-tapping task from 45 measured trials. The ΔHbO₂ and ΔHbR of the 1st NN channels calculated by the modified Beer-Lambert law were placed in the mid-points of the channels and smoothed on the scalp by the inverse of Laplacian L⁻¹ of the scalp area. The strength of the hemoglobin concentration changes in the scalp area was adjusted so that the strengths of the observation signals from the gray matter and the scalp were the same on average as $\bar{max (\sqrt{Σ_{t} {(A^{c} x^{c})}^{2}})} = \bar{max (\sqrt{Σ_{t} {(A^{s} x^{s})}^{2}})}$ , where the maximum was applied over the channels of all wavelengths and the bar represents the average over all 84 simulations. The observation noise was set to half the strength of the phantom experiment, which corresponds to one-fourth or one-fifth the strength of the human experiments. This means that we used the data from an average of around 20 trials as the observation data in this simulation study. Three values of wavelength data (780, 805, and 830 nm) were used for spectroscopy. The sampling rate was 190 ms. One estimation example of this simulation study is shown in section 4.2. The summarized results of all 84 conditions are shown in section 4.3.

The numbers of voxels in the cortical and scalp areas were N^c = 5475 and N^s = 4119, respectively. The number of measurement logical channels M was 108 (up to 3rd NN channels). Control constant β of the sensitivity-normalized Tikhonov regularization (see Appendix) was set with reference to the mean value of the sensitivity for the 25-mm deep region from the surface near the probes: logβ = mean(logρ_i), i ∈ I_{≃25 mm}, where I_{≃25 mm} denotes the set of voxel indices within a range of 24.5 to 25.5 mm deep and within the square pole made by the four normal directions of R4, R5, R9, and R10.

4.2. Estimation example

We first show one estimation example to understand the simulation conditions. Figure 6(a) shows the true spatial distribution and the true time courses of the hemoglobin concentration changes. Figure 6(b) shows the observation value resulting from the hemoglobin concentration changes shown in Fig. 6(a) and the observation noise. Figure 6(c) shows the spatial distribution and the time courses of the hemoglobin concentration changes estimated from the observation value shown in the top panel of Fig. 6(b).

Fig. 6 Example of estimation. (a) (top) Spatial distribution of the true oxygenated hemoglobin concentration change. Center position of the local hemoglobin concentration change in the cortex was (x, y, z) = (21, 41, 19) [mm]. Cross-section view of y = 41 plane is shown. (middle) Red solid and blue dashed lines represent time courses of true ΔHbO₂ and ΔHbR, respectively, of the peak value’s voxel in the cortex. (bottom) Time course of true value of ΔHbO₂ and ΔHbR in the peak value’s voxel in the scalp. (b) (top) Observation value of all channels. Observation value is the linear summation of the cortical signal, scalp signal, and observation noise. (middle) Cortical signal resulting from the hemoglobin concentration changes in cortical area. (bottom) Scalp signal resulting from hemoglobin concentration changes in scalp area. (c) (top) Spatial distribution of estimated oxygenated hemoglobin concentration change. (middle) Red solid and blue dashed lines represent time courses of estimated ΔHbO₂ and ΔHbR, respectively, of the peak value’s voxel in cortex. (bottom) Time course of estimated value of ΔHbO₂ and ΔHbR in the peak value’s voxel in scalp.

Download Full Size | PDF

The top panel of Fig. 6(a) shows the true spatial distribution of the oxygenated hemoglobin concentration change, whose center position in the cortex was (x, y, z) = (21, 41, 19) [mm]. The cross-section view of the y = 41 plane is shown. The time average of the local oxygenated hemoglobin concentration change in the cortex within the task response period from 5 to 20 sec is displayed in a color map. The square root of the time-averaged squared value of the global change in the scalp within the measurement period from −5 to 35 sec is also displayed in a color map. Since the mean value of the global change in the scalp area becomes very small, we show the square root of the integrated squared value. The time courses of ΔHbO₂ and ΔHbR of the peak voxel in the gray matter are shown in the middle panel of Fig. 6(a) and those in the scalp are shown in the bottom panel of Fig. 6(a).

The time courses of observation value Y of all channels are superimposed in the top panel of Fig. 6(b). The observation value consists of cortical signal A^cx^c, scalp signal A^sx^s, and observation noise ξ. The middle and bottom panels of Fig. 6(b) show the time courses of the cortical and scalp signals, respectively. The strengths of the cortical and scalp signals are at the same level, as explained in section 4.1. The strength of the noise is stronger than the strengths of the cortical and scalp signals.

The top panel of Fig. 6(c) shows the spatial distribution of the estimated oxygenated hemoglobin concentration change. The peak position of the estimated oxygenated hemoglobin concentration change in the cortex was (x, y, z) = (21, 41, 17) [mm], and the localization error was 2 mm. The time average of the estimated oxygenated hemoglobin concentration change in the cortex within the task response period is displayed in a color map. The square root of the time-averaged squared value of the estimated oxygenated hemoglobin concentration change in the scalp within the measurement period is also displayed in a color map. The time courses of the estimated ΔHbO₂ and ΔHbR of the peak voxel in the gray matter are shown in the middle panel of Fig. 6(c) and those in the scalp are shown in the bottom panel of Fig. 6(c). The spatial distribution and the time courses of the hemoglobin concentration changes were accurately estimated by the proposed method.

4.3. Evaluation of estimation errors

To investigate the effect of the complex shape of the five-layer head model on the estimation, we evaluated the estimation errors with 84 various positions of the localized cortical hemoglobin change. In each position, we repeated the simulation five times with different random seeds to generate the observation noise and the global hemoglobin concentration change in the scalp and thus reduce fluctuations in the error evaluations. The following values of indices for evaluating the errors were averaged five times.

Figure 7 summarizes the localization errors defined as the distance between the true and estimated peak voxels of the time-averaged change in oxygenated hemoglobin concentration in the cortex within the task response period. The point’s color in the gray matter area represents the value of the localization error when the center position of the true hemoglobin concentration change was in that point.

Fig. 7 Localization error map. A point’s color in the gray matter area represents the value of the localization error. The point’s position represents the center position of the true hemoglobin concentration change. Red and blue points represent the positions of the source and detector probes, respectively.

Download Full Size | PDF

Figure 7 shows that the localization error mostly depends on the depth from the head surface and that the positional dependency other than the depth is relatively small. Thus, the complex shape does not have a crucial effect on the estimation. Because the observation signal of the cortical hemoglobin concentration change exponentially decreases with increased depth, the depth has a dominant effect on the estimation errors.

Figure 8 summarizes the estimation errors as functions of the true depth of the cortical hemoglobin concentration change from the head surface. Figure 8(a) summarizes the localization errors shown in Fig. 7 as a function of the true depth. Figure 8(b) shows the estimated depth as a function of the true depth. If the depth is accurately estimated, the value is on the dotted line. Figure 8(c) shows the ratio of the estimated peak value to the true peak value of the time-averaged hemoglobin concentration change in the cortex within the task response period. Figure 8(d) shows the root mean squared errors (RMSE) of the hemoglobin concentration change in the scalp. Here, the squared errors of all of the voxels in the scalp and all of the time samples within the measurement time were averaged.

Fig. 8 Estimation errors as functions of true depth of cortical hemoglobin concentration change. Red x and blue + marks represent the error values of ΔHbO₂ and ΔHbR, respectively, of the 84 positions shown in Fig. 7. Red solid and blue dashed lines represent average values of x and + marks, respectively. These average values were calculated for every 1-mm range of true depth. (a) Localization error of cortical hemoglobin concentration change. (b) Depth of peak voxel of estimated cortical hemoglobin concentration change. (c) Ratio of estimated peak value to true peak value of cortical hemoglobin concentration change. (d) Root mean squared error of hemoglobin concentration change in the scalp.

Download Full Size | PDF

Figures 8(a) and (b) show that the estimation was successful and not biased toward the depth direction if the true depth was less than 25–26 mm but that the estimation suddenly became impossible above this depth. The averaged localization errors of ΔHbO₂ and ΔHbR are less than 5 mm up to the 26- and 25-mm depth conditions, respectively, and are much larger than 5 mm in deeper conditions (Fig. 8(a)). The averages of the estimated depths are consistent with the true depths within an error of 1.5 mm up to a 26-mm depth and plateaued beyond it (Fig. 8(b)). This means that the depth estimation became impossible in conditions over a 26-mm depth. Figure 8(c) shows that the hemoglobin concentration changes were underestimated because of the regularization term. The degree of underestimation depends on the depth. In deep conditions, because the regularization effect is more dominant than the observation signal, the estimated hemoglobin concentration changes became nearly zero. Figure 8(d) shows that the estimation errors in the scalp are small and only slightly affected by the position of the hemoglobin concentration change in the gray matter. The RMSE value around 0.05μM is sufficiently small compared to the hemoglobin concentration change in the scalp (see bottom panel of Fig. 6(a)). The subtle increase of RMSE in the shallow conditions represents a slight effect from the cortical signal.

To see the effect of the signal-to-noise ratio on the estimation errors, we ran the same simulation except for the observation noise level, which was halved or doubled compared to the standard noise condition shown in Fig. 8 (explained in section 4.1). Figure 9(a) shows the localization errors of ΔHbO₂ in the cortex. The red solid line is the same as that of Fig. 8(a). The dotted and dashed lines represent the localization errors with the half- and double-strength conditions of the observation noise. The depth limit until which the estimation is successful depends on the signal-to-noise ratio. The localization error is less than 5 mm up to depths of 28, 26, and 23 mm with observation noise levels of half, standard, and double strengths, respectively. Figure 9(b) shows the ratio of the estimated peak value to the true peak value of the time-averaged hemoglobin concentration change in the cortex within the task response period. The dotted, solid, and dashed lines represent the ratio with the half-, standard-, and double-strength conditions of the observation noise. This result shows that the ratio approaches 1 in the high signal-to-noise ratio condition.

Fig. 9 Estimation errors with three conditions of signal-to-noise ratio. Dotted, solid, and dashed lines represent estimation errors of ΔHbO₂ with observation noise levels of half, standard, and double strength, respectively, compared to data shown in Fig. 8. (a) Localization error of cortical hemoglobin concentration change. (b) Ratio of estimated peak value to true peak value of cortical hemoglobin concentration change.

Download Full Size | PDF

Next, we show the advantage of our proposed method over the regression method for removing the scalp artifact by comparing the localization errors. Figure 10 shows the localization errors of ΔHbO₂ in the cortex. The dotted line represents the estimation error without scalp artifact removal, where only the cortical hemodynamic change was estimated by the hierarchical Bayesian method with the sparse regularization model. The dashed line represents the estimation error with the scalp artifact removal by the regression method, where the cortical hemodynamic change was estimated by the hierarchical Bayesian method after removing scalp signals by the regression method. Here, we used the 1st principal component of the time series of the 1st NN channels for regression removal. The solid line represents the estimation error obtained by the proposed method, where both the cortical and scalp hemodynamic changes were estimated for removing the scalp artifact. This red solid line is the same as that of Fig. 8(a). These results show that the proposed method increased the estimation accuracy to a higher level than did the regression method. This is because the regression method removed the globally coherent scalp signal but it could not remove the spatiotemporally inhomogeneous scalp signal, which was better estimated by the proposed method. The inhomogeneous scalp signal made a difference in the estimation accuracy, even though the scalp hemodynamic changes were spatially very smooth in this simulation study as shown in the top panel of Fig. 6(a).

Fig. 10 Comparison of regression and proposed methods for scalp artifact removal. Each line represents localization error of ΔHbO₂ in the cortex. Dotted line represents the estimation error without scalp artifact removal, where only the cortical hemodynamic change was estimated by the hierarchical Bayesian method. Dashed line represents the estimation error with scalp artifact removal by the regression method, where the cortical hemodynamic change was estimated after removing scalp signals by the regression method. Solid line represents the estimation error obtained by the proposed method, where both the cortical and scalp hemodynamic changes were estimated for scalp artifact removal.

Download Full Size | PDF

Finally, we evaluated the robustness of the proposed method, which is important for practical use. We considered two kinds of variabilities. The first one is the optical parameter variability. We tested the case where there was a discrepancy in the optical parameters. Although we generated the data in the same manner described in section 4.1, we estimated the hemodynamic changes using a head model whose optical parameters were different from the true ones. We used the optical parameters summarized in the paper by Okada et al. [37] in the estimation process: The absorption coefficient and reduced scattering coefficient (μ_a, μ′_s(= μ_s(1 − g))) of the scalp, the skull, CSF, and the gray and white matter were (0.018 mm⁻¹, 1.9 mm⁻¹), (0.016 mm⁻¹, 1.6 mm⁻¹), (0.002 mm⁻¹, 0.001 mm⁻¹), (0.036 mm⁻¹, 2.2 mm⁻¹), and (0.014 mm⁻¹, 9.1 mm⁻¹), respectively. We used the same anisotropy g and refractive index n summarized in section 4.1. The second kind of variability involves the anatomical structure. We tested the case where there was a discrepancy in the anatomical structure. Although we generated the data in the same manner described in section 4.1, we estimated the hemodynamic changes using a simplified head model that was different from the true head model. The simplified head model had the same external boundary (scalp) but the internal boundaries were given by the parallel surfaces with the external boundary. They were determined by the shortest distance from the external boundary as follows: scalp (0–3 mm), skull (3–12 mm), CSF (12–16mm), and gray matter (16∼ mm). Consequently, there was no folding structure in the tissue boundaries. Figure 11(a) and (b) show the localization error and the estimated depth, respectively, of ΔHbO₂ in the cortex as a function of the true depth. The red solid line, which is the same as that of Fig. 8(a) and (b), represents the case where there was no discrepancy. The dotted lines represent the estimation results in the case where there was a discrepancy in the optical parameter. The results show only a slight increase in estimation error. The dashed lines represent the estimation result in the case where there was a discrepancy in the anatomical structure. The results show that the estimation was still possible with a relatively low error level despite the large difference in shape between the simplified head model and the true MRI-based head model. The localization error particularly increased in the shallow area (Fig. 11(a)), where the shape of the gray matter widely varied depending on the location. The increase in localization error in the shallow area was mostly due to the estimation error in the depth direction (Fig. 11(b)).

Fig. 11 Estimation errors when there was a discrepancy between the data generation model and the estimation model. Dotted and dashed lines represent the cases where there was discrepancy in the optical parameter and the anatomical structure, respectively. The solid lines represent the case where there was no discrepancy. (a) Localization error of ΔHbO₂ in the cortex. (b) Estimated depth of ΔHbO₂ in the cortex as a function of the true depth.

Download Full Size | PDF

5. Conclusion

We proposed a method to reconstruct the three-dimensional distributions of hemoglobin concentration changes in both cortical and scalp areas by three-dimensional diffuse optical tomography. This method spatially separates the signals from the cortex and the scalp in the estimation process using the spatial information of the optical diffusion process. Therefore, it accurately estimates the position of the cortical hemodynamic change regardless of whether the hemodynamic changes in the cortex and the scalp are temporally correlated. Different types of regularization models were introduced for the cortical and scalp areas depending on the nature of their hemodynamic responses. We applied the sparse regularization model to the cortical area for better localization [23, 38] and a smooth regularization model to the scalp area. Because the parameter controlling the smoothness of scalp regularization is automatically adjusted from the data, our method can estimate the inhomogeneous distribution of the scalp hemodynamic changes. We validated our proposed method by two-layer phantom experiments and MRI-based head-model simulations.

First, we evaluated our proposed method by applying it to the data obtained from the two-layer phantom experiment. The silicone layer of the two-layer phantom mimics the scalp and the liquid layer mimics the cortex. We found from the one-absorber experiment that we could accurately estimate both the horizontal and vertical positions of the local absorption change in the liquid layer even if the observation signal were contaminated by the superficial global absorption change in the silicone layer. Next, we confirmed the method’s high spatial resolution in the two-absorber experiment. The two local absorption changes in the liquid layer with distances shorter than the probe interval could be discriminated in spite of the presence of the superficial global absorption change in the silicone layer.

Second, we evaluated our proposed method by a simulation study using an MRI-based head model to examine the effect of the head model’s complex shape on the estimation. Our results show that the degree of estimation error mostly depends on the depth of the true position and that the positional dependency other than the depth is relatively small. Thus, the complex shape did not have a crucial effect on the estimation. We found that the depth limit up to which the estimation was successful strongly depended on the signal-to-noise ratio. We also confirmed that the proposed method has better accuracy than the regression method for removing the scalp artifact due to its ability to accurately estimate inhomogeneous hemodynamic changes in the scalp.

Finally, we tested the robustness of the proposed method with respect to the two kinds of variabilities. In the first case, we tested the robustness on the optical parameters of the head model because they may vary by individual, and thus the exact values are not known. As a result, we found only a slight increase in estimation error due to discrepancies in the optical parameters. In the second case, we tested the method’s robustness on the anatomical structure. The results showed that estimation was still possible with a relatively low level of error, even though we used a simplified head model instead of the MRI-based head model. This is probably because the diffused light is not so sensitive to the fine anatomical structure. It is also possible to register the head model with the standard MNI space, even when there is no MRI information, by the virtual spatial registration method proposed by Tsuzuki et al. [39]. These suggest that the simplified head model or the standard brain could be used for our proposed DOT method as an alternative to an MRI-based head model if some degree of estimation error were accepted.

Appendix Details of variational Bayes method

The maximization of free energy can be solved using a factorization approximation that restricts the solution space [40, 41]:

Q (X, λ, η, σ) = Q_{X} (X) Q_{λ} (λ, η, σ),

where the absorption changes in the cortical and scalp areas are summarized as X = (X^cT, X^sT)^T to simplify the expression. Under the factorization assumption, maximum free energy is obtained by alternately maximizing the free energy with respect to Q_X and Q_λ. In the first step (X-step), free energy F(Q) is maximized with respect to Q_X while Q_λ is fixed. Solution Q_X (X) becomes a normal distribution:

log Q_{X} (X) = {〈 log P (Y, X, λ, η, σ) 〉}_{Q_{λ} (λ, η, σ)} = Σ_{t = 1}^{T} Σ_{k = 1}^{K} log N (x_{k, t}; {\bar{x}}_{k, t}, Σ_{x, k}),

where x̄_k,t and Σ_x,k denote the expectation value and the covariance matrix of x_k,t, respectively. In the second step (λ-step), free energy F(Q) is maximized with respect to Q_λ while Q_X is fixed. Solution Q_λ (λ, η, σ) becomes a gamma distribution:

\begin{array}{l} log Q_{λ} (λ, η, σ) & = & {〈 log P (Y, X, λ, η, σ) 〉}_{Q_{X} (X)} \\ = & Σ_{i = 1}^{N^{c}} log Γ (λ_{i}; {\bar{λ}}_{i}, γ_{i}) + log Γ (η; \bar{η}, γ_{η}) + log Γ (σ; \bar{σ}, γ_{σ}), \end{array}

where λ̄_i and γ_i denote the expectation and shape parameter values of λ_i, η̄ and γ_η denote the expectation and shape parameter values of η, and σ̄ and γ_σ denote the expectation and shape parameter values of σ. The above X- and λ-steps are repeated until the free energy converges.

In this procedure, the free energy function can be expressed as

\begin{array}{l} F (Q) & = & - \frac{1}{2} Σ_{t = 1}^{T} Σ_{k = 1}^{K} [\bar{σ} {‖ Σ_{y, k}^{- \frac{1}{2}} (y_{k, t} - A_{k} {\bar{x}}_{k, t}) ‖}_{2}^{2} + {‖ {\bar{Λ}}^{- \frac{1}{2}} W^{- 1} {\bar{x}}_{k, t}^{c} ‖}_{2}^{2} + \bar{η} {‖ L {\bar{x}}_{k, t}^{s} ‖}_{2}^{2}] \\ + Σ_{i = 1}^{N^{c}} γ_{0 i} [log ({\bar{λ}}_{i} / {\bar{λ}}_{0 i}) - ({\bar{λ}}_{i} / {\bar{λ}}_{0 i}) + 1] + γ_{η_{0}} [log (\bar{η} / {\bar{η}}_{0}) - (\bar{η} / {\bar{η}}_{0}) + 1] \\ + Σ_{2}^{T} Σ_{k = 1}^{K} [log | Σ_{x, k} | + log | \bar{Λ} | + N^{s} log \bar{η} + log | \bar{σ} Σ_{y, k}^{- 1} |] \\ - \frac{T}{2} Σ_{k = 1}^{K} tr [Σ_{x, k} (\bar{σ} A_{k}^{T} Σ_{y, k}^{- 1} A_{k} + (\begin{matrix} W^{- 1 T} \bar{Λ} W^{- 1} & 0 \\ 0 & L^{T} \bar{η} L \end{matrix}))] + const ., \end{array}

where x̄_k,t, λ̄_i, η̄, and σ̄ represent the expectation values of x_k,t, λ_i, η, and σ with respect to trial distribution Q, respectively. Σ_x,k represents the covariance matrix of x_k,t with respect to trial distribution Q. Here, to simplify the expression, the sensitivity matrix of the cortical and scalp areas is summarized as

A_{k} = (A_{k}^{c}, A_{k}^{s})

.

Sensitivity-normalized Tikhonov regularization

We used a solution of the sensitivity-normalized Tikhonov regularization [23] as an initial value of the variational Bayes algorithm. The solution, X̂_D = min_XC_D(X; λ), is obtained by minimizing the following cost function:

C_{D} (X; λ) = Σ_{t = 1}^{T} Σ_{k = 1}^{K} [{‖ Σ_{y, k}^{- \frac{1}{2}} (y_{k, t} - A_{k} x_{k, t}) ‖}_{2}^{2} + λ {‖ D_{k}^{\frac{1}{2}} x_{k, t} ‖}_{2}^{2}],

where

D_{k} = diag (ρ_{k} + β 1_{N}),

ρ_{k, i} = {(A_{k}^{T} Σ_{y, k}^{- 1} A_{k})}_{i i} .

We determined appropriate λ value by minimizing Akaike’s Bayesian information criterion (or maximizing the marginal likelihood) [42, 43].

Estimation algorithm

The hierarchical prior and initial values are set.
In this study, we set the hierarchical prior and initial values of the cortical regularization model to γ₀ = 0 and ${\bar{λ}}_{init} = {(mean ({\hat{x}}_{D}^{c}) / max (W))}^{- 2}$ . We set the hierarchical prior and initial values of the scalp regularization model as γ_η₀ = 0 and ${\bar{η}}_{init} = {(mean ({(L {\hat{x}}_{D}^{s})}^{2}))}^{- 1}$ . We set the initial value of the noise scale parameter as σ̄_init = 1.
The free energy is maximized by alternately repeating the X-and λ-steps. In the following equations, := denotes the assignment operator; the value of the left-hand-side variable is updated by the current value of the right-hand-side equations.
We repeated these steps until the iteration count reached 2,000 for the estimation in the two-layer phantom experiments. We repeated these steps until the iteration count reached 10,000 for the estimation in the MRI-based head-model simulations.
[X-step] In the X-step, the following values are calculated:
$Σ_{k} : = {\bar{σ}}^{- 1} Σ_{y, k} + A_{k}^{c} W {\bar{Λ}}^{- 1} W^{T} A_{k}^{c^{T}} + A_{k}^{s} L^{- 1} {\bar{η}}^{- 1} L^{- 1^{T}} A_{k}^{s^{T}},$ $h^{c} : = Σ_{k = 1}^{K} diag [{\bar{Λ}}^{- 1} W^{T} A_{k}^{c^{T}} Σ_{k}^{- 1} Y_{k} Y_{k}^{T} Σ_{k}^{- 1} A_{k}^{c} W + T (1_{N^{c}} - {\bar{Λ}}^{- 1} W^{T} A_{k}^{c^{T}} Σ_{k}^{- 1} A_{k}^{c} W)],$ $h^{s} : = Σ_{k = 1}^{K} tr [{\bar{η}}^{- 1} L^{- 1^{T}} A_{k}^{s^{T}} Σ_{k}^{- 1} Y_{k} Y_{k}^{T} Σ_{k}^{- 1} A_{k}^{s} L^{- 1} + T (I_{N^{s}} - {\bar{η}}^{- 1} L^{- 1^{T}} A_{k}^{s^{T}} Σ_{k}^{- 1} A_{k}^{s} L^{- 1})],$ $S : = Σ_{k = 1}^{K} tr [{\bar{σ}}^{- 1} Σ_{y, k} Σ_{k}^{- 1} Y_{k} Y_{k}^{T} Σ_{k}^{- 1} + T ({\bar{Λ}}^{- 1} W^{T} A_{k}^{c^{T}} Σ_{k}^{- 1} A_{k}^{c} W + {\bar{η}}^{- 1} L^{- 1^{T}} A_{k}^{s^{T}} Σ_{k}^{- 1} A_{k}^{s} L^{- 1})] .$ The expectation and covariance of X with respect to Q_X can be calculated using these values as ${\bar{X}}_{k} = (\begin{matrix} {\bar{X}}_{k}^{c} \\ {\bar{X}}_{k}^{s} \end{matrix}) = {(\begin{matrix} W^{- 1 T} \bar{Λ} W^{- 1} & 0 \\ 0 & L^{T} \bar{η} L \end{matrix})}^{- 1} A_{k}^{T} Σ_{k}^{- 1} Y_{k}$ and $Σ_{x, k}^{- 1} = \bar{σ} A_{k}^{T} Σ_{y, k}^{- 1} A_{k} + (\begin{matrix} W^{- 1 T} \bar{Λ} W^{- 1} & 0 \\ 0 & L^{T} \bar{η} L \end{matrix})$ .
[λ-step] In the λ-step, the expectations and shape parameters of Λ, η, and σ with respect to Q_λ are calculated as
$γ_{i} : = γ_{0 i} + \frac{K}{2} T, γ_{i} {\bar{λ}}_{i}^{- 1} : = γ_{0 i} {\bar{λ}}_{0 i}^{- 1} + \frac{1}{2} h_{i}^{c} {\bar{λ}}_{i}^{- 1},$ $γ_{η} : = γ_{η_{0}} + \frac{K}{2} N^{s} T, γ_{η} {\bar{η}}^{- 1} : = γ_{η_{0}} {\bar{η}}_{0}^{- 1} + \frac{1}{2} h^{s} {\bar{η}}^{- 1},$ $γ_{σ} : = \frac{1}{2} M T, γ_{σ} {\bar{σ}}^{- 1} : = \frac{1}{2} S {\bar{σ}}^{- 1} .$
After the convergence criteria are met, the X-step is performed in one more step. Then the estimation result is obtained: ${\hat{X}}_{k}^{c} = W {\bar{Λ}}^{- 1} W^{T} A_{k}^{c^{T}} Σ_{k}^{- 1} Y_{k},$ ${\hat{X}}_{k}^{s} = L^{- 1} η^{- 1} L^{- 1^{T}} A_{k}^{s^{T}} Σ_{k}^{- 1} Y_{k} .$

Acknowledgments

This research was supported by a contract with the National Institute of Information and Communications Technology entitled, 'Development of network dynamics modeling methods for human brain data simulation systems' and 'Multimodal integration for brain imaging measurements'. A part of this research was supported by SRPBS, MEXT.

References and links

1. F. F. Jöbsis, “Noninvasive, infrared monitoring of cerebral and myocardial oxygen sufficiency and circulatory parameters,” Science 198, 1264–1267 (1977). [CrossRef] [PubMed]

2. Y. Hoshi and M. Tamura, “Detection of dynamic changes in cerebral oxygenation coupled to neuronal function during mental work in man,” Neurosci. Lett. 150, 5–8 (1993). [CrossRef] [PubMed]

3. T. Kato, A. Kamei, S. Takashima, and T. Ozaki, “Human visual cortical function during photic stimulation monitoring by means of near-infrared spectroscopy,” J. Cereb. Blood Flow Metab. 13, 516–520 (1993). [CrossRef] [PubMed]

4. A. Villringer, J. Planck, C. Hock, L. Schleinkofer, and U. Dirnagl, “Near infrared spectroscopy (NIRS): a new tool to study hemodynamic changes during activation of brain function in human adults,” Neurosci. Lett. 154, 101–104 (1993). [CrossRef] [PubMed]

5. B. Chance, Z. Zhuang, C. UnAh, C. Alter, and L. Lipton, “Cognition-activated low-frequency modulation of light absorption in human brain,” Proc. Natl. Acad. Sci. U.S.A. 90, 3770–3774 (1993). [CrossRef] [PubMed]

6. M. Ferrari and V. Quaresima, “A brief review on the history of human functional near-infrared spectroscopy (fNIRS) development and fields of application,” NeuroImage 63, 921–935 (2012). [CrossRef] [PubMed]

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

8. E. Kirilina, A. Jelzow, A Heine, M. Niessing, H. Wabnitz, R. Brühl, B. Ittermann, A. M. Jacobs, and I. Tachtsidis, “The physiological origin of task-evoked systemic artefacts in functional near infrared spectroscopy,” NeuroImage 61, 70–81 (2012). [CrossRef] [PubMed]

9. Y. Zhang, D. H. Brooks, M. A. Franceschini, and D. A. Boas, “Eigenvector-based spatial filtering for reduction of physiological interference in diffuse optical imaging,” J. Biomed. Opt. 10, 011014 (2005). [CrossRef]

10. S. Kohno, I. Miyai, A. Seiyama, I. Oda, A. Ishikawa, S. Tsuneishi, T. Amita, and K. Shimizu, “Removal of the skin blood flow artifact in functional near-infrared spectroscopic imaging data through independent component analysis,” J. Biomed. Opt. 12, 062111 (2007). [CrossRef]

11. Q. Zhang, E. N. Brown, and G. E. Strangman, “Adaptive filtering for global interference cancellation and real-time recovery of evoked brain activity: a Monte Carlo simulation study,” J. Biomed. Opt. 12, 044014 (2007). [CrossRef] [PubMed]

12. Q. Zhang, G. E. Strangman, and G. Ganis, “Adaptive filtering to reduce global interference in non-invasive NIRS measures of brain activation: how well and when does it work?” NeuroImage 45, 788–794 (2009). [CrossRef] [PubMed]

13. T. Yamada, S. Umeyama, and K. Matsuda, “Multidistance probe arrangement to eliminate artifacts in functional near-infrared spectroscopy,” J. Biomed. Opt. 14, 064034 (2009). [CrossRef]

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

15. R. B. Saager, N. L. Telleri, and A. J. Berger, “Two-detector Corrected Near Infrared Spectroscopy (C-NIRS) detects hemodynamic activation responses more robustly than single-detector NIRS,” NeuroImage 55, 1679–1685 (2011). [CrossRef] [PubMed]

16. L. Gagnon, R. J. Cooper, M. A. Yücel, K. L. Perdue, D. N. Greve, and D. A. Boas, “Short separation channel location impacts the performance of short channel regression in NIRS,” NeuroImage 59, 2518–2528 (2012). [CrossRef]

17. L. Gagnon, M. A. Yücel, D. A. Boas, and R. J. Cooper, “Further improvement in reducing superficial contamination in NIRS using double short separation measurements,” NeuroImage (in press).

18. T. Funane, H. Atsumori, T. Katura, A. N. Obata, H. Sato, Y. Tanikawa, E. Okada, and M. Kiguchi, “Quantitative evaluation of deep and shallow tissue layers’ contribution to fNIRS signal using multi-distance optodes and independent component analysis,” NeuroImage (in press).

19. A. P. Gibson, J. C. Hebden, and S. R. Arridge, “Recent advances in diffuse optical imaging,” Phys. Med. Biol. 50, R1–R43 (2005). [CrossRef] [PubMed]

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

21. A. T. Eggebrecht, B. R. White, S. L. Ferradal, C. Chen, Y. Zhan, A. Z. Snyder, H. Dehghani, and J. P. Culver, “A quantitative spatial comparison of high-density diffuse optical tomography and fMRI cortical mapping,” NeuroImage 61, 1120–1128 (2012). [CrossRef] [PubMed]

22. C. Habermehl, S. Holtze, J. Steinbrink, S. P. Koch, H. Obrig, J. Mehnert, and C. H. Schmitz, “Somatosensory activation of two fingers can be discriminated with ultrahigh-density diffuse optical tomography,” NeuroImage 59, 3201–3211 (2011). [CrossRef] [PubMed]

23. T. Shimokawa, T. Kosaka, O. Yamashita, N. Hiroe, T. Amita, Y. Inoue, and M. Sato, “Hierarchical Bayesian estimation improves depth accuracy and spatial resolution of diffuse optical tomography,” Opt. Express 20, 20427–20446 (2012). [CrossRef] [PubMed]

24. A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (IEEE Press, New York, 1988).

25. M. A. O’Leary, “Imaging with diffuse photon density waves,” Ph.D. Thesis, University of Pennsylvania (1996).

26. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl. 15, R41–R93 (1999). [CrossRef]

27. A. Kienle, M. S. Patterson, N. Dögnitz, R. Bays, G. Wagnières, and H. van den Bergh, “Noninvasive determination of the optical properties of two-layered turbid media,” Appl. Opt. 37, 779–791 (1998). [CrossRef]

28. C. M. Bishop, “Variational principal components,” In Proc. of ICANN 1, 509–514 (1999).

29. M. Sato, T. Yoshioka, S. Kajihara, K. Toyama, N. Goda, K. Doya, and M. Kawato, “Hierarchical Bayesian estimation for MEG inverse problem,” NeuroImage 23, 806–826 (2004). [CrossRef] [PubMed]

30. S. J. Matcher, C. E. Elwell, C. E. Cooper, M. Cope, and D. T. Delpy, “Performance comparison of several published tissue near-infrared spectroscopy algorithms,” Anal. Biochem. 227, 54–68 (1995). [CrossRef] [PubMed]

31. “FreeSurfer,” http://surfer.nmr.mgh.harvard.edu/.

32. “VBMEG, Variational Bayesian Multimodal EncephaloGraphy,” http://vbmeg.atr.jp/.

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

34. Q. Fang, “Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates,” Biomed. Opt. Express 1, 165–175 (2010). [CrossRef] [PubMed]

35. N. Lange and S. L. Zeger, “Non-linear Fourier time series analysis for human brain mapping by functional magnetic resonance imaging,” Appl. Statist. 46, 1–29 (1997). [CrossRef]

36. K. J. Worsley, C. H. Liao, J. Aston, V. Petre, G. H. Duncan, F. Morales, and A. C. Evans, “A general statistical analysis for fMRI data,” NeuroImage 15, 1–15 (2002). [CrossRef] [PubMed]

37. E. Okada and D. T. Delpy, “Near-infrared light propagation in an adult head model. I. Modeling of low-level scattering in the cerebrospinal fluid layer,” Appl. Opt. 42, 2906–2914 (2003). [CrossRef] [PubMed]

38. V. C. Kavuri, Z. J. Lin, F. Tian, and H. Liu, “Sparsity enhanced spatial resolution and depth localization in diffuse optical tomography,” Biomed. Opt. Express 3, 943–957 (2012). [CrossRef] [PubMed]

39. D. Tsuzuki, V. Jurcak, A. K. Singh, M. Okamoto, E. Watanabe, and I. Dan, “Virtual spatial registration of stand-alone fNIRS data to MNI space,” NeuroImage 34, 1506–1518 (2007). [CrossRef] [PubMed]

40. H. Attias, “Inferring parameters and structure of latent variable models by variational Bayes,” Proc. 15th Conf. on Uncertainty in Artificial Intelligence, Morgan Kaufmann, 21–30 (1999).

41. M. Sato, “Online model selection based on the variational Bayes,” Neural Comput. 13, 1649–1681 (2001). [CrossRef]

42. H. Akaike, “Likelihood and the Bayes procedure,” in Bayesian Statistics,J. M. Bernardo, M. H. De Groot, D. V. Lindley, and A. F. M. Smith, eds. (Univ. Press, Valencia, 1980), 143–166.

43. O. Yamashita, “Dynamical EEG inverse problem and causality analysis of fMRI data,” Ph.D. Thesis, Graduate University for Advanced Studies (SOKENDAI) (2004).

Extended hierarchical Bayesian diffuse optical tomography for removing scalp artifact

Abstract

1. Introduction

2. Methods

2.1. Forward model

2.2. Regularization model

2.3. Variational Bayesian method

2.4. Conversion from absorption changes to hemoglobin concentration changes

3. Validation of method by two-layer phantom experiments

3.1. Experimental conditions

3.2. Estimation example

3.3. One-absorber experiment

3.4. Two-absorber experiment

4. Validation of method by MRI-based head-model simulations

4.1. Simulation conditions

4.2. Estimation example

4.3. Evaluation of estimation errors

5. Conclusion

Appendix Details of variational Bayes method

Sensitivity-normalized Tikhonov regularization

Estimation algorithm

Acknowledgments

References and links

Cited By

Figures (11)

Equations (29)

Biomedical Optics Express