1.

Introduction

Measuring the attenuation of near-infrared and visible red light that has passed through tissue can provide physiologically interesting information. In near-infrared spectroscopy (NIRS), changes in optical signal measured at multiple wavelengths are converted into changes in oxygenated (${\mathrm{HbO}}_2$) and deoxygenated hemoglobin (HbR). In diffuse optical imaging (DOI), multiple measurements are used to reconstruct the spatial distribution of optical parameters or changes in these parameters.

During the past two decades, NIRS and DOI have been used for imaging the hemodynamic responses associated with brain activity.^1–3 In comparison with other functional imaging modalities such as functional magnetic resonance imaging (fMRI), positron emission tomography (PET), and electro- and magnetoencephalography (EEG and MEG), the advantages of optical imaging include good temporal resolution for the hemodynamic response, ability to separately measure both HbR and ${\mathrm{HbO}}_2$, and relatively portable and light-weight and low-cost equipment that can be used at the bed side. The disadvantages include modest spatial resolution and limited sensitivity to deep-laying tissues. The last point is especially true when imaging brain activations in the adult, since the size of head allows only measurements using reflection geometry.

The high sensitivity of the measurement to superficial tissues causes optical imaging of brain activations to be highly susceptible to artifacts due to metabolic and hemodynamic effects that occur either systemically or locally in the scalp.^4–8 Strategies to reduce the artifacts have included using multi-distance measurements and a layered anatomical model to estimate contributions of the various tissue layers⁵ using the short-distance measurements as regressors⁶ for removing component of signal change due to surface effects, and using analysis of moments of the temporal point spread function (TPSF),⁹ or principal component (PCA) and independent component analysis (ICA) to separate contributions from the brain and the superficial layers,¹⁰ as well as using peripheral measurements such as pulse and blood pressure and include these in a physiological model to predict the optical signal.¹¹ These methods have achieved some success in reducing the artifacts, but often they either do not completely correct the artifacts or depend on measurements or data that may not be available. In Ref. 7, it was shown that task-related activation of the sympathetic nervous system may be responsible for changes in optical signal that are not correlated with systemic measurements such as pulse and blood pressure. Methods that rely on using the short distance measurements to estimate the surface component in observed signals depend on accurate modelling of the physical situation, and the surface component in the signal is in many cases assumed to be due to uniform change in the scalp. The analysis of moments of the TPSF requires time-domain (TD) data, and while using a suitable optode grid information about the localization of the brain activation can be achieved, the method is best suited to optical spectroscopy. Methods based on ICA and PCA require no explicit physical or physiological model, but only achieve a separation of signal components, and reconstructing the interesting parameters requires further analysis. They are also best suited to spectroscopic imaging.

In this study, we propose to use the Bayesian approximation error method^12,13 to remove the artifacts caused by systemic and superficial effects. Earlier, the Bayesian approximation error approach has been applied to optical tomography in order to recover from different types of uncertainties and modeling errors, including handling pure discretization errors,¹⁴ verification with real data,¹⁵ marginalization over the unknown inhomogeneous scattering coefficient¹⁶ using Born approximation,¹⁷ and using diffusion approximation instead of radiative transfer model with low scattering targets.¹⁸

In contrast to many of the previously proposed methods, while the possibility of hemodynamic changes in the scalp is included in our model, its effect is not explicitly included in the forward model. By posing a certain probability distribution model for changes in the absorption at the superficial structures, we compute estimates for the mean and covariance of the approximation error noise process and then employ these estimated statistics in the reconstruction process for weighting the measurements to compensate for the uncertainty caused by the superficial structures in the reconstructed image of a brain activation. The method can accommodate a smoothly varying change in the scalp blood flow, to allow for local variation in the surface vasculature and effects such as blushing that may not be spatially uniformly distributed. Rather than spectroscopic analysis, the method presented reconstructs spatially resolved images of brain activations. The method is applicable to averaged brain activation data, and no time series data is used. We consider mean time of flight and signal intensity from TD data, and the results are also applicable for the frequency-domain case.

2.

Methods

2.1.

Monte Carlo Simulation

We used a voxel based Monte Carlo model.^19,20 to generate the simulated measurements. We considered TD measurement, and used the signal intensity and mean time of flight of photons.

The anatomical model was based on an MR image of an adult male. The image was segmented into scalp, skull, cerebrospinal fluid (CSF), and gray and white brain matter.²⁰ Literature-derived values of optical parameters were used for the tissue classes (see Table 1).

Table 1

Optical properties of tissue types.21 Given are reduced scattering coefficient (μs′), absorption (μa) coefficient. The ’true parameters,’ used in generating the simulated measurements, are estimates for wavelength λ=830  nm. ’WP1’ and ’WP2’ are two sets of deviated optical parameters that were used to test the behaviour of the reconstruction method when the baseline optical properties are inaccurately known. Refractive index n=1.4 was used for all tissue types. GM=Gray matter, WM=White matter, CSF=Cerebrospinal fluid.

Tissue	True parameters		WP1		WP2
	${\mu }_s^{\prime }({\mathrm{mm}}^{-1})$	${\mu }_{\mathrm{a}}({\mathrm{mm}}^{-1})$	${\mu }_s^{\prime }({\mathrm{mm}}^{-1})$	${\mu }_{\mathrm{a}}({\mathrm{mm}}^{-1})$	${\mu }_s^{\prime }({\mathrm{mm}}^{-1})$	${\mu }_{\mathrm{a}}({\mathrm{mm}}^{-1})$
Scalp	1.4	0.016	2.0	0.020	1.0	0.015
Skull	1.25	0.024	1.8	0.024	1.0	0.020
GM	0.65	0.0186	1.0	0.02	0.6	0.015
WM	1.1	0.014	1.4	0.02	1.0	0.010
CSF	0.025	0.004	0.025	0.004	0.025	0.004

The optode array used in the simulation consisted of 8 sources and 15 detectors, located on the scalp at the back of the head overlying the visual cortex (see Fig. 1). The distance from a source to the closest detectors was approximately 2 cm.

Fig. 1

The optode array used in the simulations. Source locations marked with ‘$X$,’ detector locations marked with ‘$O$’.

We simulated localized activations in the visual cortex, represented by local increases in optical absorption. In our example cases, described in Sec. 3, we used two approximately spherical simulated activations with radii 7 mm (Cases 1 and 2) and 11 mm (Case 3). The activations were constrained to the gray matter, and had homogeneous contrast $\delta {\mu }_a=0.002{\mathrm{mm}}^{-1}$. Changes in the blood flow of the scalp were simulated by spatially smoothly varying absorption changes with peak contrast of approximately $\delta {\mu }_a=0.00022{\mathrm{mm}}^{-1}$ or 11% of the contrast in the brain.

The change in the TD data during a brain activation or change in the scalp blood flow was recorded, and these difference data were used for reconstructing the brain activations. We only considered the difference between the base state and a single activated state, i.e., no time series data were used. Gaussian white noise was added to the simulated difference data. The noise consists of two parts corresponding to stochastic noise due to dark counts and thermal noise which is same for all channels, and shot noise which is proportional to the square root of the number of photons observed. The effect of shot noise on the signal-to-noise ratio is thus proportional to reciprocal of the square root of signal intensity.²² At the shortest measurement distance the noise was set at 0.15% for the intensity measurement, and 0.15 ps for the mean time of flight measurement, with equal contributions from shot noise and stochastic noise. The noise as function of measurement distance is shown in Fig 2.

Fig. 2

The simulated noise level. On the left, standard deviation of noise in the intensity measurement as percentage of measured intensity. On the right, standard deviation of noise in the mean time measurements in picoseconds. The x axis shows distance between source and detector in millimeters.

3.

Results

We tested the correction procedure in the following four cases. These were

The above cases were selected to examine the performance of the approximation error model correction in a range of situations that can occur in optical brain activation imaging. In order to isolate the problem of artifacts due to superficial contamination of the optical signal, this was the only modeling inaccuracy that was considered and compensated using the approximation error method.

In Figs. 3–6, reconstructions from the simulated cases described above are shown. The figures show reconstructed contrast at the surface of the cortex. In each case, we compare the results obtained with and without the aid of the approximation error method. To facilitate the evaluation of the performance of our method, we also computed quantitative metrics of the quality of the reconstruction. These metrics were the localization error of the center-of-mass of the reconstructed absorptive changes, contrast-to-noise ratio (CNR), and the L2-norm of the error of the three-dimensional reconstruction. For the calculation of the center-of-mass and the CNR, we defined the area of the reconstructed brain activation ${\mathbf{R}}_{\mathbf{B}\mathbf{A}}$ as the three-dimensional continuous volume around the peak of the activation, with absorption change at least half of the maximum at the peak (26-connectivity was used). The center-of-mass was calculated as ${\overrightarrow{r}}_{\mathrm{CoM}}=\sum _{r_i\in {\mathbf{R}}_{\mathbf{B}\mathbf{A}}}{\overrightarrow{r}}_i·\delta {\mu }_a({\overrightarrow{r}}_i)/\sum _{r_i\in {\mathbf{R}}_{\mathbf{B}\mathbf{A}}}\delta {\mu }_a({\overrightarrow{r}}_i)$, where the sums go over positions ${\overrightarrow{r}}_i$ and reconstructed absorption changes $\delta {\mu }_a({\overrightarrow{r}}_i)$ of all voxels $i$ in ${\mathbf{R}}_{\mathbf{B}\mathbf{A}}$. The CNR was calculated as $\mathrm{CNR}=\text{mean}[\delta {\mu }_a({\overrightarrow{r}}_i\in {\mathbf{R}}_{\mathbf{B}\mathbf{A}})]/\mathrm{std}[\delta {\mu }_a({\overrightarrow{r}}_i\notin {\mathbf{R}}_{\mathbf{BA},\mathrm{ext}})]$, or the ratio between the mean contrast of the reconstructed activation and the standard deviation of the background. The background is defined as the volume outside the extended activated area ${\mathbf{R}}_{\mathbf{BA},\mathbf{ext}}$, obtained by dilating ${\mathbf{R}}_{\mathbf{BA}}$ by two millimeters to allow for the diffuse nature of the reconstruction. The L2-error norm was calculated for the difference $\delta {\overrightarrow{\mu }}_{a,\mathrm{recon}}-\delta {\overrightarrow{\mu }}_{a,\text{target}}$, where the vectors contain the voxel-by-voxel absorption change in the reconstruction and in the target distribution containing the simulated brain activation in the relevant volume. Since the absorption change is almost exlusively reconstructed near the surface of the cortex, all metrics were computed for the cortical volume under the optodes down to 5 mm depth. For the localization error and L2-error norm, a smaller value means a smaller error in the reconstruction, and is thus better. For the CNR, a larger value means a more reliably reconstructed brain activation, and is better. These numerical data are presented in Table 2. Also presented in the table are results from a case where the effect of inaccurate knowledge of the background optical properties are inaccurately known (Fig. 7, see below).

Fig. 3

Reconstruction results for Case 1 (no absorptive change in the scalp). (a) Simulated target brain activation, radius 7 mm, (b) reconstruction of the brain activation using the approximation error method assuming absorption changes may be present in the scalp, and (c) reconstruction of the brain activation using conventional constrained reconstruction (Color online only).

Fig. 4

Reconstruction results for Case 2 (predominantly positive absorptive change in the scalp). (a) Simulated target brain activation (radius 7 mm), (b) simulated absorptive change in the scalp, (c) reconstruction of the brain activation using the approximation error method, and (d) reconstruction of the brain activation using conventional constrained reconstruction (Color online only).

Fig. 5

Reconstruction results for Case 3 (predominantly negative absorptive change in the scalp). (a) Simulated target brain activation (radius 11 mm), (b) simulated absorptive change in the scalp, (c) reconstruction of the brain activation using the approximation error method, and (d) reconstruction of the brain activation using conventional constrained reconstruction (Color online only).

Fig. 6

Reconstruction results for Case 4 (no brain activation). (a) Simulated absorptive change in the scalp, (b) reconstruction using the approximation error method, and (c) reconstruction using conventional constrained reconstruction (Color online only).

Fig. 7

Reconstruction results for Case 3 (negative absorptive change in the scalp), assuming inaccurate knowledge of the baseline optical properties. Two inaccurate estimates of the optical properties were considered (WP1, WP2, see Table 1). (a) Simulated absorptive change in the brain, (b) absorptive change in the scalp, (c) reconstruction of the brain activation using inaccurate estimate WP1, and the approximation error model correction, (d) reconstruction of the brain activation using inaccurate estimate WP1, without model correction, (e) reconstruction of the brain activation using inaccurate estimate WP2, and the approximation error model correction, and (f) reconstruction of the brain activation using inaccurate estimate WP2 without model correction (Color online only).

Table 2

Quantitative metrics for evaluating the quality of the reconstructions. Cases as explained above. Case 1 involves no effect in the scalp. For cases 2 and 3, positive and negative absorptive changes are present in the scalp, respectively. WP1 and WP2 stand for two sets of inaccurately know background optical properties. Results using the approximation error to aid reconstruction (AE) and results from conventional reconstruction (Conv.) presented for each case.

In Fig. 3, the results for Case 1 (no absorption change in the scalp) are presented. In this case, the corrected and noncorrected reconstruction results are very similar. This similarity is also shown in the similar values of the quantitative metrics (see Table 2) for Case 1. In Case 2 (predominantly positive absorption change in the scalp, see Fig. 4), the approximation error can be shown to improve the reconstruction, revealing the location of the simulated brain activation in the lower right hand side of the occipital cortex, whereas in the noncorrected reconstruction the reconstruction resembles the pattern of absorption change in the scalp. For Case 2, all quantitative metrics are clearly better for the approximation error corrected case. In Case 3 (predominantly negative absorption change in the scalp, see Fig. 5), the approximation error correction allows reasonably accurate localization of the brain activation, whereas without the correction, the reconstruction only reflects the absorption change in the scalp [compare with Figs. 5(b) and 6(c)]. This is also reflected in the quantitative metrics. The very low CNR for the noncorrected reconstruction suggests that the brain activation is not realiably reconstructed. In Fig. 6, the reconstructions are shown from Case 4 (no brain activation). The approximation error corrected reconstruction is fainter and appears noisier, whereas the uncorrected reconstructions shows peak absorption increase bilaterally in the inferior part of the occipital cortex, resembling the pattern of absorption change in the scalp. The result suggests that the approximation error method may reduce the occurrence of false detections of brain activations that may be due to the combination of extra-cerebral blood flow changes and the sensitivity pattern of the measurement.

For performing the reconstruction, and for calculating the approximation error statistic, we need an estimate of the baseline optical properties $x_0$ in the whole domain. This is used as the linearization point. The reconstruction results depend on the accuracy of the estimate. In order to evaluate the sensitivity of the proposed method to inaccurately known values of the baseline optical properties, we considered two inaccurate estimates WP1 and WP2 (see Table 1). In the first case (WP1), the scattering coefficient ${\mu }_s^{\prime }$ was overestimated by between 27% and 54%, and the absorption coefficient ${\mu }_a$ was overestimated by between 0% and 30% in all tissues expect CSF, where no inaccuracy was simulated. In the second case (WP2), ${\mu }_s^{\prime }$ was underestimated by between 8% and 40%, and ${\mu }_a$ was underestimated by between 6% and 29% in all tissues expect CSF. We then re-performed the reconstruction in Case 3, using the inaccurate estimates WP1 and WP2 of baseline optical properties. The inaccurate estimates were used as the linearization point for calculating the Jacobian, and to calculate the approximation error statistics, using the same distribution of ${\mu }_a$ changes in the scalp with respect to baseline as for the accurate baseline optical properties. These were combined with the respective Jacobians to perform approximation error corrected reconstructions. In Fig. 7, the effect of inaccurately known baseline optical properties is shown. The Case 3 shown in Fig. 5 is considered. The approximation error correction is shown to improve the reconstruction also in this case. The quantitative metrics shown in Table 2 show similar relation between the approximation error corrected and conventional reconstructions as in Case 3 using accurate baseline optical properties.

4.

Conclusions

We have shown the feasibility of using the approximation error method to correct for artifacts in optical brain activation imaging, when the artifacts have been caused by changes in systemic and extra-cerebral hemodynamics. The only assumption used was that the absorptive changes in the scalp are spatially smooth. The method has some degree of tolerance for other inaccuracies in the model apart from unknown signal contribution from the scalp, such as inaccurately known baseline optical properties of tissues of the head. The reconstruction results naturally depend on the level of noise in the measurements. In this study, we used a level of noise that can be achievable experimentally.

The important question in optical brain activation imaging is whether a change seen in the measurement signal has been caused by changes in the cortex, or in the scalp. The results of this study suggest that the approximation error method can be used to reduce artifacts of extra-cerebral origin, while also taking into account another important source of artifacts, which is the measurement noise. We expect that combining the approximation error method with physiological modeling can make estimation of the origin of observed signals more reliable.

In the present study, we only considered artifacts due to the presence of unknown perturbations in blood circulation in the scalp. Other inaccuracies in the model were not considered in the model correction. The contribution of this paper is to show by simulations what the performance of the approximation error method in compensating for the above mentioned perturbations may be. Subsequent papers will validate the method using phantom experiments and study of the performance of the method on experimental brain activation data.