|
1.IntroductionMeasuring 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 () 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 , 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 layers5 using the short-distance measurements as regressors6 for removing component of signal change due to surface effects, and using analysis of moments of the temporal point spread function (TPSF),9 or principal component (PCA) and independent component analysis (ICA) to separate contributions from the brain and the superficial layers,10 as well as using peripheral measurements such as pulse and blood pressure and include these in a physiological model to predict the optical signal.11 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 method12,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,14 verification with real data,15 marginalization over the unknown inhomogeneous scattering coefficient16 using Born approximation,17 and using diffusion approximation instead of radiative transfer model with low scattering targets.18 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.Methods2.1.Monte Carlo SimulationWe 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.20 Literature-derived values of optical parameters were used for the tissue classes (see Table 1). Table 1Optical 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.
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. 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 . Changes in the blood flow of the scalp were simulated by spatially smoothly varying absorption changes with peak contrast of approximately 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.22 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. 2.2.ReconstructionThe perturbation Monte Carlo (pMC)23–25 method was used for reconstruction. In pMC, the Jacobian matrix that describes the linear change in optical signal due to changes in optical parameters is derived from a Monte Carlo simulation at the linearization point. The reconstruction is obtained by minimizing an objective function: where is the measured data during the brain activation, and is the reference measurement during baseline conditions, both obtained from simulation in this study. is the Jacobian matrix computed at linearization point . We distinguish the subset of as the parameter and region of interest, for which the difference from base state is reconstructed. This may include all or part of . While is the Jacobian matrix for the parameter to be reconstructed, in the region of interest , it depends on the parameters in the whole domain . The covariance matrix of the measurement error is used for weighting the measurements. is the regularization term. In this work, we used , where and are regularization parameters.We constrain the reconstruction to the brain, and, as is common in brain activation imaging, we reconstruct only for the optical absorption. The reconstructed quantity is thus the three-dimensional distribution of change in optical absorption in the cortex, relative to the baseline value . Due to the high sensitivity of the measurement to superficial tissues, the anatomical constraint is necessary in order to obtain a solution within the brain. This method has been shown to provide superior topographical reconstruction26 compared with unconstrained reconstruction, when the observed change in the optical measurement is due to a hemodynamic change within the brain. However, when extra-cerebral hemodynamic changes contribute to the measured change in optical signal, the constraint can lead to erroneous results. The same error can also occur in a nonconstrained reconstruction, but in this case, most of the reconstructed change will be placed in the scalp and skull regardless of its origin. In this study, we explore the potential of the approximation error method to correct for these artifacts. 2.3.Approximation Error MethodApproximation error method can be used to reduce or remove artifacts in the solution of the inverse problem (reconstruction of parameters) that are due to inaccuracy of the forward model.12,14 The general principle is to examine the difference between an ‘exact’ model and its computational approximation . Using this difference, we can define the ‘augmented model’The realization of the approximation error is unknown and depends on the unknown parameters . However, by using a Gaussian approximation for the joint statistics of parameters and approximation error , the reconstruction problem can be premarginalized with respect the approximation error similarly as the (weighted) least squares estimation is premarginalized with respect the measurement errors (for details, see Ref. 16). By a making a further technical approximation that and are mutually uncorrelated, the computation of the maximum a posteriori (MAP) estimate amounts to where is the mean of modeling error, and is a sum of the signal covariance due to measurement noise and modeling error . The mean and covariance of the modeling error can be calculated using a set of samples from an a priori distribution of the parameters .In the present study, we consider the modeling error in the measured optical signal due to changes in blood flow of the scalp during a brain activation imaging, and linear reconstruction. The approximate linear model for the signal change is then and the accurate linear model that includes absorption changes in the scalp is and are the Jacobian matrices corresponding to the brain and the surface structures, and and are the absorption changes in these regions. Thus, the modeling error consists of the part of the accurate model corresponding the superficial structuresUsing the approximation error correction, Eq. (1) becomes where the objective function is minimized with respect to absorption in the brain , while the possibility of absorption changes in the scalp is included in the modeling error .We calculated the approximation error statistics using 500 draws from a smoothness prior model of absorptive changes in the scalp similarly as in Ref. 15. The prior model consisted of a uniform background component and inhomogeneities. The 2 s.t.d. limits for the uniform part and the inhomogeneities were and , respectively. Both the homogeneous and inhomogeneous part were centered on the baseline optical parameters for the scalp. The correlation length of the inhomogeneous part was 10 cm. The correlation length represents a rough prior estimate for the inhomogeneity of hemodynamic effects in the scalp. This corresponds to smoothly varying absorption change that will typically be either negative or positive under or all or most of the measurement optode array (roughly 7 cm by 8 cm). We consider this to be a realistic assumption since changing state of systemic circulation and autonomic nervous system typically do not produce sharply localized effects in the blood flow. However, we consider allowing for spatial variation in the absorption change important to accommodate for variations in local vasculature and its responses to activation of the autonomous nervous system. 3.ResultsWe 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 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 , where the sums go over positions and reconstructed absorption changes of all voxels in . The CNR was calculated as , 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 , obtained by dilating by two millimeters to allow for the diffuse nature of the reconstruction. The L2-error norm was calculated for the difference , 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). Table 2Quantitative 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 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 was overestimated by between 27% and 54%, and the absorption coefficient was overestimated by between 0% and 30% in all tissues expect CSF, where no inaccuracy was simulated. In the second case (WP2), was underestimated by between 8% and 40%, and 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 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.ConclusionsWe 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. AcknowledgmentsThis work was supported by the European Community’s Seventh Framework Programme (FP7/2007-2013) under grant agreement n 201076, and by the Academy of Finland (Projects 119270, 136220, 140731, 140984, 218183, 213476, and 250215 Finnish Centre of Excellence in Inverse Problems Research). ReferencesY. HoshiM. Tamura,
“Dynamic multichannel near-infrared optical imaging of human brain activity,”
J. Appl. Physiol., 75
(4), 1842
–1846
(1993). JAPYAA 0021-8987 Google Scholar
A. VillringerU. Dirnagl,
“Coupling of brain activity and cerebral blood flow: basis of functional neuroimaging,”
Cerebrovasc. Brain Metab. Rev., 7
(3), 240
–276
(1995). 1040-8827 Google Scholar
M. A. FranceschiniD. A. Boas,
“Noninvasive measurement of neuronal activity with near-infrared optical imaging,”
NeuroImage, 21
(1), 372
–386
(2004). http://dx.doi.org/10.1016/j.neuroimage.2003.09.040 NEIMEF 1053-8119 Google Scholar
I. Tachtsidiset al.,
“False positives in functional near-infrared topography,”
Adv. Exp. Med. Biol., 645 307
–314
(2009). http://dx.doi.org/10.1007/978-0-387-85998-9 AEMBAP 0065-2598 Google Scholar
S. UmeyamaT. Yamada,
“Monte Carlo study of global interference cancellation by multidistance measurement of near—infrared spectroscopy,”
J. Biomed. Opt., 14
(6), 064025
(2009). http://dx.doi.org/10.1117/1.3275466 JBOPFO 1083-3668 Google Scholar
L. Gagnonet al.,
“Improved recovery of the hemodynamic response in diffuse optical imaging using short optode separations and state-space modeling,”
NeuroImage, 56
(3), 1362
–1371
(2011). http://dx.doi.org/10.1016/j.neuroimage.2011.03.001 NEIMEF 1053-8119 Google Scholar
T. Takahashi,
“Influence of skin blood flow on near-infrared spectroscopy signals measured on the forehead during a verbal fluency task,”
NeuroImage, 57
(3), 991
–1002
(2011). http://dx.doi.org/10.1016/j.neuroimage.2011.05.012 NEIMEF 1053-8119 Google Scholar
L. Minatiet al.,
“Intra- and extra-cranial effects of transient blood pressure changes on brain near-infrared spectroscopy (NIRS) measurements,”
J. Neurosci. Method., 197
(2), 283
–288
(2011). http://dx.doi.org/10.1016/j.jneumeth.2011.02.029 JNMEDT 0165-0270 Google Scholar
H. Wabnitzet al.,
“Time-resolved near-infrared spectroscopy and imaging of the adult human brain,”
Adv. Exp. Med. Biol., 662 143
–148
(2010). http://dx.doi.org/10.1007/978-1-4419-1241-1 AEMBAP 0065-2598 Google Scholar
J. VirtanenT. NoponenP. Merila¨inen,
“Comparison of principal and independent component analysis in removing extracerebral interference from near-infrared spectroscopy signals,”
J. Biomed. Opt., 14
(5), 054032
(2009). http://dx.doi.org/10.1117/1.3253323 JBOPFO 1083-3668 Google Scholar
S. G Diamondet al.,
“Dynamic physiological modeling for functional diffuse optical tomography,”
NeuroImage, 30
(1), 88
–101
(2006). http://dx.doi.org/10.1016/j.neuroimage.2005.09.016 NEIMEF 1053-8119 Google Scholar
J. KaipioE. Somersalo, Statistical and Computational Inverse Problems, Springer, New York
(2004). Google Scholar
J. KaipioE. Somersalo,
“Statistical inverse problems: discretization, model reduction and inverse crimes,”
J. Comput. Appl. Math., 198
(2), 493
–504
(2007). http://dx.doi.org/10.1016/j.cam.2005.09.027 JCAMDI 0377-0427 Google Scholar
S. R. Arridgeet al.,
“‘Approximation errors and model reduction with an application in optical diffusion tomography,”
Inv. Probl., 22
(1), 175
–195
(2006). http://dx.doi.org/10.1088/0266-5611/22/1/010 INPEEY 0266-5611 Google Scholar
V. Kolehmainenet al.,
“Approximation errors and model reduction in three-dimensional diffuse optical tomography,”
J. Opt. Soc. Am. A., 26
(10), 2257
–2268
(2009). http://dx.doi.org/10.1364/JOSAA.26.002257 JOAOD6 0740-3232 Google Scholar
V. Kolehmainenet al.,
“Marginalization of uninteresting distributed parameters in inverse problems—application to diffuse optical tomography,”
Int. J. Uncertain. Quantification, 1
(1), 1
–17
(2011). http://dx.doi.org/10.1615/Int.J.UncertaintyQuantification.v1.i1 2152-5080 Google Scholar
T. Tarvainenet al.,
“Corrections to linear methods for diffuse optical tomography using approximation error modelling,”
Biomed. Opt. Express, 1
(1), 209
–222
(2010). http://dx.doi.org/10.1364/BOE.1.000209 BOEICL 2156-7085 Google Scholar
T. Tarvainenet al.,
“Approximation error approach for compensating modelling errors between the radiative transfer equation and the diffusion approximation in diffuse optical tomography,”
Inv. Probl., 26
(1), 015005
(2010). http://dx.doi.org/10.1088/0266-5611/26/1/015005 INPEEY 0266-5611 Google Scholar
D. A. BoasA. M. DaleM. A. Franceschini,
“Diffuse optical imaging of brain activation: approaches to optimizing image sensitivity, resolution, and accuracy,”
NeuroImage, 23 S275
–S288
(2004). http://dx.doi.org/10.1016/j.neuroimage.2004.07.011 NEIMEF 1053-8119 Google Scholar
J. Heiskalaet al.,
“Modeling anisotropic light propagation in a realistic model of the human head,”
Appl. Opt., 44
(11), 2049
–2057
(2005). http://dx.doi.org/10.1364/AO.44.002049 APOPAI 0003-6935 Google Scholar
G. StrangmanM. A. FranceschiniD. A. Boas,
“Factors affecting the accuracy of near-infrared spectroscopy calculations for focal changes in oxygenation parameters,”
NeuroImage, 18
(4), 865
–879
(2003). http://dx.doi.org/10.1016/S1053-8119(03)00021-1 NEIMEF 1053-8119 Google Scholar
E. Hillman,
“Experimental and theoretical investigations of near infrared tomographic imaging methods and clinical applications,”
University College London,
(2002). Google Scholar
C. Hayakawa,
“Perturbation Monte Carlo methods to solve inverse photon migration problems in heterogeneous tissues,”
Opt. Lett., 26
(17), 1335
–1337
(2001). http://dx.doi.org/10.1364/OL.26.001335 OPLEDP 0146-9592 Google Scholar
P. KumarR. M. Vasu,
“Reconstruction of optical properties of low-scattering tissue using derivative estimated through perturbation Monte-Carlo method,”
J. Biomed. Opt., 9
(5), 1002
–1012
(2004). http://dx.doi.org/10.1117/1.1778733 JBOPFO 1083-3668 Google Scholar
J. Heiskalaet al.,
“Probabilistic atlas can improve reconstruction from optical imaging of the neonatal brain,”
Opt. Express, 17
(17), 14977
–14992
(2009). http://dx.doi.org/10.1364/OE.17.014977 OPEXFF 1094-4087 Google Scholar
D. A. BoasA. M. Dale,
“Simulation study of magnetic resonance imaging-guided cortically constrained diffuse optical tomography of human brain function,”
Appl. Opt., 44
(10), 1957
–1968
(2005). http://dx.doi.org/10.1364/AO.44.001957 APOPAI 0003-6935 Google Scholar
|