PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptNIH Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Neuroimage. Author manuscript; available in PMC Feb 1, 2013.
Published in final edited form as:
PMCID: PMC3254723
NIHMSID: NIHMS330962
Short separation channel location impacts the performance of short channel regression in NIRS
Louis Gagnon,abc Robert J. Cooper,a Meryem A. Yücel,a Katherine L. Perdue,d Douglas N. Greve,a and David A. Boasab
aAthinoula A. Martinos Center for Biomedical Imaging, Department of Radiology, Massachusetts General Hospital, Harvard Medical School, Charlestown, MA, USA
bHarvard-MIT Division of Health Sciences and Technology, Cambridge, MA, USA
cDepartment of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA, USA
dThayer School of Engineering, Dartmouth College, Hanover, NH, USA
Near-Infrared Spectroscopy (NIRS) allows the recovery of cortical oxy-and deoxyhemoglobin changes associated with evoked brain activity. NIRS is a back-reflection measurement making it very sensitive to the superficial layers of the head, i.e. the skin and the skull, where systemic interference occurs. As a result, the NIRS signal is strongly contaminated with systemic interference of superficial origin. A recent approach to overcome this problem has been the use of additional short source-detector separation optodes as regressors. Since these additional measurements are mainly sensitive to superficial layers in adult humans, they can be used to remove the systemic interference present in longer separation measurements, improving the recovery of the cortical hemodynamic response function (HRF). One question that remains to answer is whether or not a short separation measurement is required in close proximity to each long separation NIRS channel. Here, we show that the systemic interference occurring in the superficial layers of the human head is inhomogeneous across the surface of the scalp. As a result, the improvement obtained by using a short separation optode decreases as the relative distance between the short and the long measurement is increased. NIRS data was acquired on 6 human subjects both at rest and during a motor task consisting of finger tapping. The effect of distance between the short and the long channel was first quantified by recovering a synthetic hemodynamic response added over the resting-state data. The effect was also observed in the functional data collected during the finger tapping task. Together, these results suggest that the short separation measurement must be located as close as 1.5 cm from the standard NIRS channel in order to provide an improvement which is of practical use. In this case, the improvement in Contrast-to-Noise Ratio (CNR) compared to a standard General Linear Model (GLM) procedure without using any small separation optode reached 50 % for HbO and 100 % for HbR. Using small separations located farther than 2 cm away resulted in mild or negligible improvements only.
Keywords: Near-Infrared Spectroscopy, Systemic Interference, Short Optode Separations, Kalman filtering
Over the past 15 years, Near-Infrared Spectroscopy (NIRS) (Villringer et al., 1993; Obrig and Villringer, 2003; Gibson et al., 2005; Hillman, 2007; Hoshi, 2007) has emerged as a complement to functional Magnetic Resonance Imaging (fMRI) for mapping the hemodynamic response associated with cerebral activity. NIRS non-invasively measures the temporal variations of the two dominant chromophores in the near-infrared window: oxygenated hemoglobin (HbO) and deoxygenated or reduced hemoglobin (HbR).
The advantages of NIRS for the investigation of brain activity include the measurement of both HbO and HbR concentrations, its low cost, and its portability. The portability of NIRS enables long-term monitoring of the hemodynamic reponse associated with, for instance, epilectic activity at the bedside (Roche-Labarbe et al., 2008). Disadvantages of NIRS include modest spatial resolution of the order of one to three centimeters and limited penetration depth (Boas et al., 2004).
A common problem with NIRS recordings is the presence of strong physiology-based systemic interference in the signal which reduces the accuracy of NIRS for detecting brain activation. This interference arises from cardiac activity, respiration and other homeostatic processes (Obrig et al., 2000; Tonorov et al., 2000; Diamond et al., 2009; Payne et al., 2009). The contribution of this interference in the NIRS signal is amplified because the light is both introduced and collected at the surface of the scalp. This back-reflection geometry makes NIRS very sensitive to the superficial layers of the head which contain no brain signal but exhibit strong systemic fluctuations. As such, the NIRS signal is often dominated by systemic interference occurring in the superficial layers of the head including the scalp and the skull.
Several methods have been described in the literature to remove the systemic interference from NIRS measurements. Some post-processing algorithms include low pass filtering (Franceschini et al., 2003; Jasdzewski et al., 2003), principal component analysis (Zhang et al., 2005; Franceschini et al., 2006) and wavelet filtering (Lina et al., 2008; Matteau-Pelletier et al., 2009; Jang et al., 2009; Lina et al., 2010). Multi-distance NIRS measurements with layered models and path length weighted methods have also being investigated (Umeyama and Yamada, 2009a,b; Yamada et al., 2009). Other methods include the subtraction of a NIRS channel acquired in a non-activated region of the brain from the signal of interest to reduce the systemic interference (Franceschini et al., 2003).
A more refined version of this method is to simultaneously collect additional NIRS measurements using short source-detector (SD) separation channels (generally shorter than 1 cm), which are sensitive to superficial layers only (Saager and Berger, 2005). Assuming that the signal collected with these additional short separation measurements is dominated by the same systemic interference present in the longer SD channels, the small separation signals can be used as regressors to filter the systemic interference from the longer SD measurements. Several algorithms have been developed to perform the regression of the small separation measurements. These include linear minimum mean square estimation (LMMSE) (Saager and Berger, 2005, 2008; Gregg et al., 2010; Saager et al., 2011), adaptive filtering (Zhang et al., 2007a,b, 2009) and state-space modeling with Kalman filter estimation (Gagnon et al., 2011).
An important question which was not addressed in these previous papers is the impact that the relative location of the short and long SD channels has on the performance of the short separation method. If good performance is obtained using a short separation channel located far away from the standard long SD channel, then a single short separation channel can be used as a regressor for all longer SD channels on the probe. On the other hand, the performance of the short separation method potentially worsen as the relative distance between the short and the long SD channel increases. In this case several short separation channels would be required and only those closest to the long SD channels would be suitable for regression.
The main contribution of this paper is to quantify the performance of the short separation method as a function of the relative distance between long SD NIRS channels (3 cm) containing the brain signal and short separation (1 cm) channels used as regressors. We investigated this relationship with both simulations and real functional data. NIRS measurements including several short separation channels spread across the probe were acquired on 6 human subjects. The simulations were performed by adding a synthetic hemodynamic response to the resting-state NIRS data. NIRS signals were also collected during a series of finger tapping blocks for each of the 6 subjects. In both cases, the performance of the short separation regression was characterized for different short SD regressors located at different distances from the standard 3 cm channel.
2.1. Experimental data
For this study, 6 healthy adult subjects were recruited. The Massachusetts General Hospital Institutional Review Board approved the study and all subjects gave written informed consent. Data were collected using a TechEn CW6 system operating at 690 and 830 nm. The NIRS probe contained 5 sources and 12 detectors as shown in Fig. 1a. This source-detector geometry resulted in 14 long SD measurements (3 cm) and 7 short SD measurements (1 cm). A set of 200 μm-core fibers was used for the short separation detector optodes to avoid saturation of the photodiode. These fibers are illustrated in orange in Fig. 1. An alternative to avoid photodiode saturation could be the use of standard NIRS fibers with optical filters at the tip of the probe to attenuate light intensity. The probe was secured over the left motor region of each subject as illustrated in Fig. 1b. One of the short separation measurements was acquired over the forehead. In this probe, the relative distances center-to-center between the short and the long channels take the values 1.4, 1.7, 2.4, 3.3, 4.2, 5.2 or 6.2 cm. Examples are given for each case in Fig. 1c. The forehead short separation channel was located more than 10 cm away from any 3 cm channel.
Figure 1
Figure 1
(A) Geometry of the optical probe. Two different SD separations were used: 1 cm and 3 cm. (B) Location of the probe on the subjects. The probe was secured over the motor region. (C) Examples of short and long channel pairs. With this probe arrangement, (more ...)
During the experiment, subjects were sitting in a comfortable chair in front of a computer screen with a black background. The functional runs were divided as shown in Fig. 2. Each run lasted 390 seconds and contained six blocks of 30 s finger tapping interleaved with 30 s resting blocks. Three functional runs were acquired for each subject. During the resting blocks, a small 0.5-by-0.5 cm white square located at the middle of the screen appeared and the subjects were asked to fixate on this square. During the finger tapping blocks, the instruction “tap your fingers” was displayed in white characters on the computer screen using the Psychophysics toolbox in Matlab (Brainard, 1997). At that time, the subjects were asked to touch their right thumb with each of the fingers of their right hand alternately at a rate of 3 Hz. Following the three functional runs, three baseline runs of 5 minutes each were acquired. During the baseline runs, the subjects were asked to simply close their eyes and remain still.
Figure 2
Figure 2
Overview of the finger tapping protocol. A run consisted of 6 blocks of 30 seconds of finger tapping interleaved with 30 seconds of rest. Each runs started and ended with a 30 second resting period. 3 functional runs were acquired for each of the 6 subjects. (more ...)
2.2. Data processing
An overview of the procedure is shown in Fig. 3. Both the short and long SD measurements were bandpass filtered at 0.01–1.25 Hz. Even though the data will be further low pass filtered at 0.5 Hz in the processing stream, it is important to keep the 0.5–1.25 Hz frequency band here, since most of the cardiac oscillations are contained in this frequency band. These cardiac oscillations are strongly present in both the short and the long SD measurements and this increases the baseline correlation between the short and the long separation channel. These cardiac oscillations guide the dynamic estimation of the superficial contamination to more accurately estimate the HRF. As weve shown recently (Gagnon et al., 2011), prefiltering the cardiac oscillations reduces the performance of the dynamic estimation and results in a poorer estimate of the HRF. For both the simulations and the real functional data analysis, the Kalman filter algorithm was used to regress the short separation measurement and recover the hemodynamic response simultaneously. This algorithm was described in detail previously (Gagnon et al., 2011) and only the salient points are reviewed here.
Figure 3
Figure 3
Schematic of the NIRS data analysis. The NIRS data from both the 1 cm and the 3 cm separation channels were first converted to HbO and HbR time courses and bandpass filtered. HbO and HbR were analyzed separately. The 1 cm and 3 cm bandpass time courses (more ...)
The hemodynamic response was modeled by
equation M1
(1)
where bi [n] are normalized Gaussian functions with a standard deviation of 0.5 s and their means separated by 0.5 s. Nw is the number of Gaussian functions used to model the hemodynamic response and was set to 15 for our simulations (section 2.3) and 79 for our finger tapping data to recover the HRF over 0–8 sec and 0–40 sec respectively. The signal in the 3 cm separation channel y3 [n] was modelled by a linear combination of the 1 cm separation signal y1 [n] and the brain response yb [n]. The expression for the 3 cm signal is given by
equation M2
(2)
with
equation M3
(3)
and where u [n] is the onset vector which is a binary vector taking the value 1 when n corresponds to a time when the stimulus was presented and 0 otherwise. It is to note that u [n] is equal to 1 only at the onset of the stimulus and not throughout the duration of the stimulus.
The variable a is the dynamic weight used to model the superficial signal in the 3 cm separation channel as a linear combination of the short separation signal. Only a single time delay was taken from the short separation channel to model the superficial signal in the 3 cm channel since this has been shown to result in a better performance in our previous paper (Gagnon et al., 2011). The states to be estimated by the Kalman filter were the weight of the superficial contribution a and the weights of the temporal bases wi. All these weights were assumed to be time-varying. Eqs. (1), (2) and (3) can be re-written in state-space form:
equation M4
(4)
equation M5
(5)
where w [n] and v [n] are the process and the measurement noise respectively. x [n] is the nth instance of x given by
equation M6
(6)
The quantity I is an Nw + 1 by Nw + 1 identity matrix and C [n] is a 1 by Nw + 1 vector given by
equation M7
(7)
where “*” denotes the convolution operator. The estimate x̂ [n] at each sample n is then computed using the Kalman filter (Kalman, 1960) followed by the Rauch–Tung–Striebel smoother (Rauch et al., 1965).
The convergence of the Kalman filter depends on the initial estimate of the state vector x̂ [0]. To overcome this problem, x̂ [0] was set to the values obtained using a static least-squares estimator as in (Gagnon et al., 2011) to ensure a fast convergence. Moreover, to overcome the problem of selecting a good initial guess for the state covariance estimate P [0], the Kalman filter algorithm was run twice and the initial covariance estimate for the second run was set to the final covariance estimate of the first run. This process makes the performance of the filter almost insensitive to the initial covariance estimate. For the first pass of the Kalman filter, we set P [0] to an identity matrix with diagonal entries of 1×10−1 for the temporal basis states and 5×10−4 for the superficial contribution state. The process noise covariance Q only contained nonzero terms on the diagonal elements. Those diagonal terms were set to 2.5×10−6 for the temporal basis states and 5×10−6 for the superficial contribution state. The measurement noise covariance R was set to an identity matrix scaled by 5×10−2. These values were extensively studied in our previous paper (Gagnon et al., 2011) and multiplying or dividing these values by factor of 10 did not significantly affect the performance of our method. The Kalman filter algorithm was then processed with the following prediction-correction recursion (Gelb, 1974):
equation M8
(8)
equation M9
(9)
equation M10
(10)
equation M11
(11)
equation M12
(12)
After the Kalman algorithm was applied twice, the Rauch–Tung–Striebel smoother was applied in the backward direction (Haykin, 2001):
equation M13
(13)
with Nt the number of time points in the data. The complete time course of the filtered brain signal ŷb [n] containing the estimated hemodynamic response ĥ [n] was then reconstructed for each sample time n using the first Nw final state estimates x̂b = [w1 ··· wNw]T and the temporal basis set contained in C [n]
equation M14
(14)
This reconstructed filtered brain signal time course ŷb [n] was further low pass filtered at 0.5 Hz to remove any cardiac fluctuations potentially present in the time course and the final estimate of the hemodynamic response ĥ [n] was obtained either by applying a standard General Linear Model (GLM) procedure (without any cosine bases or short separation regressor) containing the same temporal basis function as in Eq. (1) or by block-averaging ŷb [n]. More details can be found in our previous paper (Gagnon et al., 2011).
2.3. Simulations
For each baseline measurement, the changes in optical density were converted to changes in hemoglobin concentrations using the modified Beer-Lambert relationship (Cope and Deply, 1998; Deply et al., 1988; Boas et al., 2004). A pathlength correction factor of 6 and a partial volume correction factor of 50 were applied (Huppert et al., 2006a,b). The variance in all 252 (6 subjects × 3 runs × 14 pairs) baseline HbO and HbR time courses from the 3 cm measurements were then computed. To ensure a uniform distribution of the noise in our simulations, only the time courses showing a variance below 25 μM2 were kept in the analysis, corresponding to 28.2 % of the data (71 of the 252 baseline time courses). We have also tested and confirmed that our method was working for higher levels of noise. Due to the non uniform level of noise across the probe, this threshold of 25 μM2 was required in order to compare all distances on equal footing.
Ten individual evoked responses were added over all 71 selected 3 cm baseline measurements at random onset times with an inter-stimulus interval taken randomly from a uniform distribution (10–30 sec). This procedure was repeated 30 times for each baseline measurement to create 30 simulated time courses with 30 different onset times and ensure reproducible averaged results. The duration of the synthetic response was 8 seconds. The HbO time course increased by 15 μM at the peak while the HbR time course decreased by 7 μM. The synthetic hemodynamic response was the same used in our previous paper (Gagnon et al., 2011). The resulting 2130 time courses (71 time courses × 30 simulated runs) were then bandpass filtered (0.01–1.25 Hz) and passed to the Kalman filter algorithm (Fig. 3) using each of the seven short separation (1 cm) measurements available as a regressor. The HRF was also recovered using a standard GLM with a set of cosine basis with 64 s period cutoff (Friston et al., 2000) for comparison (no short separation used). This resulted in 17,040 estimated HRFs (2130 time courses × 8 regressors (7 short separations + 1 standard GLM with cosine basis set)). The HbO and HbR responses were recovered independently. For each 1 cm–3 cm combination, the baseline R2 coefficient before adding the synthetic HRF to the 3 cm channel was computed.
For each short separation used, the relative center-to-center distance between the 3 cm and the short-separation channel was computed. With the probe shown in Fig. 1a, the possible relative distances were 1.4, 1.7, 2.4, 3.3, 4.2, 5.2 or 6.2 cm as well as > 10 cm for the forehead channel and are illustrated in Fig. 1c.
The quality of each recovered HRF was quantified by three different metrics: (1) the Pearson correlation coefficient R2 between the true synthetic HRF (tHRF) and the recovered HRF (rHRF), (2) the mean square error (MSE) between tHRF and rHRF and (3) the Contrast-to-noise ratio (CNR) defined as the amplitude of rHRF divided by the root-mean-square (RMS) of the residual of tHRF and rHRF
equation M15
(15)
The average for each of these three metrics across all the recovered HRFs for each specific relative distance was computed and the results were compared to the corresponding averaged metrics obtained from the HRFs recovered with the standard GLM (no short separation) using a two-tailed paired t-test. As in our previous paper (Gagnon et al., 2011), we used a paired t-test to resolve for small systematic differences. For the Pearson R2 metric, the average was taken after applying a Fisher transformation and the resulting average was then back transformed. This comparison was performed for all 8 relative distances (1.4, 1.7, 2.4, 3.3, 4.2, 5.2, 6.2 cm and forehead).
This entire procedure was then repeated after introducing a time-lag in the short separation channel. For each 1 cm–3 cm combination, the cross-correlation function between the two channels before adding the synthetic HRF was also computed and a time-lag corresponding to the maximum of the cross-correlation function was applied to the short separation measurement. This time-lag could be any number in the interval {−Nt, Nt} with Nt the number of time point in the NIRS time course but typical values obtained from our data ranged from −4 to 4 seconds for both HbO and HbR. The values for R2, MSE and CNR obtained by introducing a time-lag were also compared with the zero-lag values with a two-tailed paired t-test.
The cross-correlation function used to identify the optimal time-lag was normalized such that the zero-lag value corresponded to the Pearson R2 coefficient
equation M16
(16)
with σy3 the standard deviation of y3 [n]. The maximum of this normalized cross-correlation function is the equivalent to shifting one of the channels by the optimal time-lag before computing the standard correlation and thus we will refer to this value as the optimal time-lag correlation for the rest of the text. To avoid any confusion, we will refer to the standard R2 correlation as the zero-lag correlation. The zero-lag and the optimal time-lag correlations were also compared using a two-tailed paired t-test of their Fisher transformed values.
2.4. Functional data
The functional data were analysed in the same way as above with the Kalman filter, but the HRFs were recovered from 0 to 40 seconds after the stimulus onsets. Each 3 cm channel was analyzed using each of the seven short separation channels available and also with a standard block-average for comparison.
3.1. Baseline correlation
The correlation (Pearson R2) between the baseline NIRS time courses are shown in Fig. 4. In panels (A) and (B), the correlation between the 3 cm separation and the short separation channels are plotted as a function of their relative distance on the probe. These values are identified by the label “no time-lag” in the legends. The optimal time-lag correlation values are also plotted and identified by the label “optimum time-lag” in the legends. We observed a fast decay of both the zero-lag and the optimal time-lag correlations as the distance between the two channels went from 1.4 to 2.4 cm and the correlation then plateaued from 2.4 to 6.2 cm. This trend was observed for both HbO and HbR. The optimal time-lag correlation values obtained were slightly higher (but significant at the p < 0.05 level, two-tailed paired t-test) than the zero-lag correlation for all relative distances on the probe. It is good to re-emphasize that the statistical test is pairwise such that a small but constant difference across the sample is marked as significant. The optimal time-lags obtained increased with increasing relative distances indicating that this slight improvement in correlation was real and not due to potential artifact in the processing. Finally, the increases in correlation obtained by introducing time-lag were slightly more prominent for HbR than HbO.
Figure 4
Figure 4
Effect of the relative distance on the initial baseline correlation between the channels. The baseline data were bandpass filtered between 0.01 and 1.25 Hz before the R2 correlation was computed. The values labeled “optimal time-lag” were (more ...)
In panels (C) and (D) of Fig. 4, both the zero-lag and the optimal time-lag correlations between two short separation channels are plotted as a function of their relative distance on the probe. A similar fast decay was observed as the relative distance between the two channels increased from 1 to 2 cm on the probe and then the correlation plateaued for longer distances. The values for the optimal-time lag correlations in this case were significantly higher (p < 0.05, two-tailed paired t-test) compared to the zero-lag correlations.
3.2. Simulation results
The results for the synthetic HRF simulations are shown in Figs. 5, ,66 and and77 for the R2, MSE and CNR metric respectively. On panels (A) and (B) of all three figures, the three metrics are plotted as a function of the relative distance between the 3 cm separation and the short separation channel used as a regressor. Values obtained by introducing a time-lag in the short separation channel are also shown as well as the corresponding values obtained using a standard GLM. We observed a fast decrease of the improvement obtained by the Kalman filter as the relative distance between the 3 cm and the short separation channel was increased from 1.4 to 2.4 cm. The performance then plateaued for longer relative distances. Both the R2 (Fig. 5) and the CNR (Fig. 7) decreased as the relative distance between the long- and short separation channels was increased, while the MSE (Fig. 6) increased. Using a short separation channel located 1.4 cm away from the channel containing the synthetic HRF resulted in a mean increase in CNR of 50 % for HbO and 100 % for HbR relative to the GLM method. Using a short separation channel located farther than 2 cm away from the channel containing the synthetic HRF resulted in significant (p < 0.05, two-tailed paired t-test) but negligible improvements of the order of a few percent compared to the standard GLM procedure. Again, we re-emphasize that the statistical test is pairwise such that a small but constant difference across the sample is marked as significant. However the magnitude of the difference is small.
Figure 5
Figure 5
Effect of the relative distance on the correlation between the recovered HRF and the true HRF. Panels (A)–(B) show the recovered R2 as a function of the distance between the short and the long NIRS channel. Panels (C)–(D) show the recovered (more ...)
Figure 6
Figure 6
Effect of the relative distance on the MSE between the recovered HRF and the true HRF. Panels (A)–(B) show the MSE as a function of the distance between the short and the long NIRS channel. Panels (C)–(D) show the MSE as a function of (more ...)
Figure 7
Figure 7
Effect of the relative distance on the Contrast-to-noise ratio (CNR) defined in Eq. 15. Panels (A)–(B) show the CNR as a function of the distance between the short and the long NIRS channel. Panels (C)–(D) show the CNR as a function of (more ...)
On panels (C) and (D) of Figs. 5, ,66 and and7,7, the same R2, MSE and CNR results are plotted as a function of the baseline zero-lag correlation (Pearson R2) between the 3 cm and the short separation channels. Results obtained by introducing a time-lag in the short separation channel are also plotted as a function of the baseline optimal time-lag correlation. We observed a linear relationship between the improvement obtained with the Kalman filter and both the baseline zero-lag correlations and optimal time-lag correlations between the two channels. A baseline correlation greater than 0.8 resulted in a mean improvement in CNR of 50 % and 100 % compared to the standard GLM for HbO and HbR respectively.
3.3. Functional data results
Each run of finger tapping was analysed independently for each subject. The SD pair showing the strongest functional response was selected manually for each subject. To avoid any bias towards the Kalman filter method, both the Kalman filter and the standard block-average results were taken into account independently. The criteria for selecting the responses were a sustained increase in HbO and a sustained decrease in HbR (Cui et al., 2010), as well as a sustained increase in HbT to avoid pial vein washout contamination. Based on these criteria, the selected channel from the block-average matched the one selected from the Kalman filter result for each subject, although the HRF recovered with the block-average showed weak activation for three of the six subjects. For each individual subject, the HRFs recovered from the selected channel are shown for the first run in Fig. 8. Results from a single run are presented to illustrate the power of our method and the high CNR achieved with only 6 individual finger tapping blocks. Results from the second and the third run were very similar. Columns 1–4 illustrate the corresponding HRFs (same 3 cm channels) recovered using: (1) a block-average with no small-separation channel, (2) the Kalman filter with the small-separation channel located in the forehead as a regressor, (3) the Kalman filter with the short separation channel located 2.4 cm away as a regressor and (4) the Kalman filter with the closest short separation channel (located 1.4 cm away) as a regressor.
Figure 8
Figure 8
Finger tapping results from a single run containing 6 individual blocks. For each subject, the trace from the channel containing the strongest functional response is shown. (1) Standard block-average (2) Kalman filter with the forehead short separation (more ...)
For subjects 1,4 and 6, the HRFs recovered without using the closest short separation channel (columns 1–3) showed weak activation patterns only. However, the activation became clear on the same channels using the Kalman filter together with the closest short separation channel (column 4). For subjects 2, 3 and 5, a clear activation was present whether the Kalman filter was used or not. For subjects 2, 3 and 4, a strong artifact present between 15 and 30 s. was removed using the Kalman filter with the closest short separation channel. Removing this artifact made the HRFs more constant during the duration of the stimulus (0–30 s).
4.1. Systemic interference measured by NIRS is inhomogeneous across the scalp
Systemic interference measured in NIRS has been termed “global” interference previously in the literature (Saager and Berger, 2005; Zhang et al., 2007a,b; Umeyama and Yamada, 2009a; Saager et al., 2011). In contradiction, our present results indicate that systemic interference is actually inhomogeneous across the surface of the scalp, that is, the correlation between systemic interference measured at two different locations decreases with the increasing relative distance between the two measurements. Although the short separation channel measurements might contain some cortical signal, Monte Carlo simulations have shown that this contribution is negligible for a SD distance of 1 cm (Zhang et al., 2007a). As such, panels (C) and (D) of Fig. 4 indicate that the origin of this decorrelation is located in the superficial layers of the head and therefore, is not due to autoregulation mechanisms occurring in the brain tissue (Payne et al., 2009).
Although introducing a time-lag re-established one third of the correlation between the two short separation measurements, the other two thirds of the R2 correlation was still lost after introducing time delays. This finding indicates that only part of the correlation decay can be explained by simple transit time effects across different locations on the scalp. The systemic interferences measured in NIRS are oscillatory processes containing three dominant components (Lina et al., 2008). These are the cardiac pulsations around 1 Hz, respiratory oscillations around 0.4 Hz and other low frequency oscillations (including Mayer waves (Julien, 2007)) around 0.1 Hz. Analyses similar to the one presented in Fig. 4 were performed with the NIRS data bandpass filtered at 0.01–0.2 Hz, 0.2–0.5 Hz and at 0.5–3 Hz (see supplementary figures 9, 10 and 11). These frequency bands correspond to the low frequency, respiratory and cardiac oscillations respectively. These analyses revealed a decay in correlation with increasing relative distances in all these three frequency bands. Although the correlation decayed, it never reached zero even for low frequency oscillations, which is in agreement with recent findings by Tong and Frederick (2010). Up to 3/4 and 1/2 of the correlation lost in the 0.2–0.5 Hz band and the 0.5–3 Hz band respectiveley could be re-established by introducing a time-lag. However, introducing a time-lag in the 0.01–0.2 Hz frequency band resulted in only negligible improvements in correlation. These findings are in agreement with a recent paper from Tian et al (2011) indicating that cardiac fluctuations (~1 Hz) are more global while low frequency oscillations (~0.1 Hz) are less spatially coherent.
From our results, one can conclude that (1) slow oscillations are inhomogeneous across the surface of the scalp and (2) a significant proportion of the correlation decay in the higher frequency bands (cardiac and respiration) is attributed to transit time effects across different spatial regions. These phase mismatches of the cardiac and respiratory pulsation over different locations arise potentially from spatial heterogeneity of the vasculature such as blood vessel length, orientation, size, depth and dilation (Zhang et al., 2007a, 2009). However, 1/4 and 1/2 of the correlation lost in the respiration and cardiac frequency band respectively could not be re-established by introducing a time-lag and future studies will be required to investigate the origin of this correlation decay.
4.2. Impact on the short separation method
This fast decrease in correlation with increasing relative distance has an important impact on the performance of the short separation method. Panels (C) and (D) of Figs. 5, ,66 and and77 illustrate that all three metrics used (R2, MSE and CNR) varied linearly with the baseline correlation, whether a time-lag was used or not. Since the baseline correlation decreased with the relative distance as shown in Fig. 4, we expected the performance of the short separation method assessed with these three metrics to decrease with the relative distance. This expectation was confirmed by our simulations and are shown in panels (A) and (B) of Figs. 5, ,66 and and7.7. For relative distances larger than 2 cm, only mild improvements were obtained using our short separation approach compared to the standard GLM method. No decrease in performance was observed using our Kalman filter with any of the available short separation. This is consistent with our previous findings (Gagnon et al., 2011) that the Kalman filter improves or doesn’t change recovery of the hemodynamic response. At worst, the recovered response will be the same as the one recovered with a standard GLM with no small separation used. To obtain larger improvements that are useful in practice, the short separation channel must be located no more than 1.5 cm from the 3 cm channel from which the HRF is to be recovered.
For large NIRS probes containing several long SD measurements spread over several centimeters, our results indicate that multiple short separation channels are required in order to combine each long-separation measurement with a short separation channel within a 1.5 cm radius. In this case, our simulations indicate that the improvement in CNR using the short separation method is of the order of 50 % for HbO and 100 % for HbR, as shown respectively in panels (A) and (B) of Fig. 7. As in our previous paper (Gagnon et al., 2011), we observed an improvement for both HbO and HbR, in contradiction with Zhang et al. (2009) where no improvement was observed for HbR using an adaptive filter method. In the simulations of our previous paper (Gagnon et al., 2011), we also observed a decrease in performance for HbR using an adaptive filter. As we showed, our Kalman filter approach overcomes this pitfall by regressing the short separation measurement and simultaneously estimating the hemodynamic response.
The necessity of using several short separation channels on the probe was also confirmed by our finger tapping experiment. Using the short separation channel located at 2.4 cm from the 3 cm SD measurement, the HRF obtained from subjects 1, 4 and 6 showed a weak activation pattern only, as shown in Fig. 8. On the other hand, the activation became very clear on these same channels using the short separation channel located 1.4 cm away. The baseline correlations between the short separation and the long-separation channels were around 0.3 for the short separation channel located 2.4 cm away and around 0.5 for the short separation channel located 1.4 cm away. This re-emphasizes the fact that the initial baseline correlation between the short separation and the 3 cm channel is an important factor in determining the performance of the Kalman filter algorithm. In practice, this baseline correlation can be computed to predict the impact of using short optode separations. In our previous paper (Gagnon et al., 2011), we showed that the presence of a hemodynamic response in the 3 cm channel does not impact the baseline initial correlation between the short separation and the 3 cm channel. This occurs because the contribution of systemic interference in NIRS largely dominates the contribution of the hemodynamic response. This was confirmed here with our real functional data.
4.3. How many regressors should be used?
In all the simulations presented in this paper, we used a single regressor at a time. We also ran simulations using multiple regressors simultaneously but all of them resulted in a significantly higher MSE and lower CNR compared to using the single short separation channel. This is explained by a potential overfitting of the data. When more than one regressor is used at the same time, we start fitting noise in the data which introduces errors in the estimation of the hemodynamic response. The same thing occurred in our previous paper (Gagnon et al., 2011) when a single regressor was used but multiple time delays from this regressor were used in the regression. To avoid overfitting the data and to obtain the maximum power from the short separation method, we have found that on average a single regressor performs better than two or more.
4.4. Alternatives for modelling the physiological interference
Alternative methods for modelling the physiological interference have been proposed in the literature. Prince et al. (2003) used a set of sine and cosine functions to model the oscillatory behaviour of the systemic physiology. The linear coefficients of these temporal bases were included as additional states in the state-space model. Alternatively, Abdelnour and Huppert (2009) used a set of sine functions only but included the phase as an additional state. These methods were implemented and compared with the short separation method. In both cases, the short separation approach performed better compared to these modelling techniques. These models contain a higher number of degrees of freedom (i.e. larger number of state) which potentially introduces crosstalk between the state corresponding to the hemodynamic response and the state corresponding to the systemic physiology. This phenomenon degrades the estimation of the hemodynamic response. The short separation method reduces the number of degrees of freedom and reduces crosstalk by measuring directly the systemic interference in the superficial layers of the head and therefore, results in a more accurate estimation of the hemodynamic response.
4.5. Future studies
In this work, a small detector fiber was placed in proximity of each source fiber, resulting in a single small separation channel for every long SD channel. An interesting question is whether or not an additional short separation channel located in proximity of the detector fiber would further improve the recovery of the hemodynamic response. Doing so would maximize the overlap between the pathlength of the small separation channels and the longer SD channel, and might result in a further improvement of the small separation method. However, care must be taken to ensure that the additional source fibers do not inadvertently saturate the detectors and to ensure that the fiber optic probe remains flexible for efficient coupling with the scalp. We have begun to investigate the improvement of additional short separation measurements on the long separation source and long separation detector simultaneously. Our preliminary evidence indicates that as expected the results get better. For sure, we can still use a single short separation regressor but choose the one which has a higher correlation with the long separation measurement; either the short separation coincident with the long separation source, or the long separation detector. When we use both short separation regressors we have to ensure that we are not over fitting the data. Sometimes using both short separation regressors will actually degrade our estimate of the HRF, and thus statistical model testing needs to be implemented to determine if one or two regressors should be used.
In this study, we have determined that the position of the short separation NIRS channel relative to the long-separation channel impacts the performance of the short separation regression method to improve the recovery of the hemodynamic response in NIRS. We showed that the relative distance between the channel of interest and the regressor must be less than 1.5 cm to have a meaningful impact on the recovery of the hemodynamic response. In this case, improvements in CNR were of the order of 50 % for HbO and 100 % for HbR compared to the standard GLM approach. When a short separation channel located farther than 2 cm was used as regressor, only minor improvements were obtained compared to the standard GLM method, which are of little practical use. This decrease in performance for longer relative distances results from a decrease in the baseline correlation between the channel of interest and the regressor. Our results indicate that this correlation decay is due to spatially inhomogeneous hemodynamics in the superficial layers of the head.
Figure 9
Figure 9
Low frequency band (Supplementary figure). Effect of the relative distance on the initial baseline correlation between the channels. The baseline data were bandpass filtered between 0.01 and 0.2 Hz (low frequency oscillations) before the R2 correlation (more ...)
Figure 10
Figure 10
Respiratory band (Supplementary figure). Effect of the relative distance on the initial baseline correlation between the channels. The baseline data were bandpass filtered between 0.2 and 0.5 Hz (respiratory oscillations) before the R2 correlation was (more ...)
Figure 11
Figure 11
Cardiac band (Supplementary figure). Effect of the relative distance on the initial baseline correlation between the channels. The baseline data were bandpass filtered between 0.5 and 3 Hz (cardiac oscillations) before the R2 correlation was computed. (more ...)
Research Highlights
  • Short separation channel location impacts the performance of the method
  • Systemic physiology is inhomogeneous across the scalp
  • CNR improvement of 50% (HbO) and 100% (HbR) if regressor located within 1.5cm
  • Only mild improvement if short channel located farther than 1.5 cm
Supplementary Material
Acknowledgments
The authors are grateful to Drs. Evgeniya Kirilina, Yunjie Tong and Blaise deB. Frederick for fruitful discussions about cerebrovascular physiology. This work was supported by NIH grants P41-RR14075 and R01-EB006385. L. Gagnon was supported by the Fonds Quebecois sur la Nature et les Technologies.
Footnotes
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
  • Abdelnour AF, Huppert TJ. Real-time imaging of humain brain function by near-infrared spectroscopy using an adaptive general linear model. NeuroImage. 2009;46:133–143. [PMC free article] [PubMed]
  • Boas D, Dale AM, Franceschini MA. Diffuse optical imaging of brain activation: Approaches to optimizing image sensitivity, resolution and accuracy. Neuroimage. 2004;23:S275–S288. [PubMed]
  • Brainard DH. The psychophysics toolbox. Spatial Vision. 1997;10:433–436. [PubMed]
  • Cope M, Deply DT. System for long-term measurement of cerebral blood flow and tissue oxygenation on newborn infants by infrared transillumination. Med Biol Eng Comput. 1998;26:289–294. [PubMed]
  • Cui X, Bray S, Reiss AL. Functional near infrared spectroscopy (NIRS) signal improvement based on negative correlation between oxygenated and deoxygenated hemoglobin dynamics. NeuroImages. 2010;49:3039–3046. [PMC free article] [PubMed]
  • Deply DT, Cope M, van der Zee P, Arridge S, Wray S, Wyatt J. Estimation of optical pathlength through tissue from direct time of flight measurement. Phys Med Biol. 1988;33:1433–1442. [PubMed]
  • Diamond SG, Perdue KL, Boas DA. A cerebrovascular response model for functional neuroimaging including dynamic cerebral autoregulation. Math Biosci. 2009;220(2):102–117. [PMC free article] [PubMed]
  • Franceschini M, Fantini S, Thompson J, Culver J, Boas D. Hemodynamic evoked response of the sensorimotor cortex measured noninvasively with near-infrared optical imaging. Psychophysiology. 2003;40 (4):548–560. [PMC free article] [PubMed]
  • Franceschini M, Joseph DK, Huppert TJ, Diamond SG, Boas D. Diffuse optical imaging of the whole head. J Biomed Opt. 2006;11(5):054007. [PMC free article] [PubMed]
  • Friston KJ, Josephs O, Zarahn E, Holmes AP, Rouquette S, Poline J. To smooth or not to smooth? bias and efficiency in fMRI time-series analysis. NeuroImage. 2000;12 (2):196–208. [PubMed]
  • Gagnon L, Perdue K, Greve DN, Goldenholz D, Kaskhedikar G, Boas DA. Improved recovery of the hemodynamic response in diffuse optical imaging using short optode separations and state-space modeling. NeuroImage. 2011;56 (3):1362–1371. [PMC free article] [PubMed]
  • Gelb A. Applied Optimal Estimation. MIT Press; 1974.
  • Gibson A, Hebden J, Arridge S. Recent advances in diffuse optical imaging. Phys Med Biol. 2005;50:R1–R43. [PubMed]
  • Gregg NM, White BR, Zeff BW, Berger AJ, Culver JP. Brain specificity of diffuse optical imaging: improvements from superficial signal regression and tomography. Font in Neuroenergetics. 2010;2 (14):1–8. [PMC free article] [PubMed]
  • Haykin S. Kalman filtering and neural Networks. John Wiley & Sons; New York: 2001.
  • Hillman EMC. Optical brain imaging in vivo: techniques and applications from animal to man. J Biomed Opt. 2007;12(5):051402. [PMC free article] [PubMed]
  • Hoshi Y. Functional near-infrared spectroscopy: current status and future prospects. J Biomed Opt. 2007;12(6):062106. [PubMed]
  • Huppert T, Hoge R, Dale AM, Franceschini M, Boas D. Quantitative spatial comparison of diffuse optical imaging with blood oxygen level-dependent and arterial spin labeling-based functional magnetic resonance imaging. J Biomed Opt. 2006a;11(6):064018. [PMC free article] [PubMed]
  • Huppert T, Hoge R, Diamond S, Franceschini M, Boas D. A temporal comparison of BOLD, ASL, and NIRS hemodynamic responses to motor stimuli in adult humans. NeuroImage. 2006b;29:368–382. [PMC free article] [PubMed]
  • Jang KE, Tak S, Jung J, Jang J, Jeong Y, Ye JC. Wavelet minimum description length detrending for near-infrared spectroscopy. J Biomed Opt. 2009;14(3):034004. [PubMed]
  • Jasdzewski G, Strangman G, Wagner J, Kwong KK, Poldorack RA, Boas DA. Differences in the hemodynamic response to event-related motor and visual paradigms as measured by near- infrared spectroscopy. Neuroimage. 2003;20 (1):479–488. [PubMed]
  • Julien C. The enigma of Mayer waves: facts and models. Cardiovascular Research. 2007;70:12–21. [PubMed]
  • Kalman RE. A new approach to linear filtering and prediction problems. Journal of Basic Engineering. 1960;82 (Series D):35–45.
  • Lina JM, Dehaes M, Matteau-Pelletier C, Lesage F. Complex wavelets applied to diffuse optical spectroscopy for brain activity detection. Optics Express. 2008;16 (2):1029–1050. [PubMed]
  • Lina J-M, Matteau-Pelletier C, Dehaes M, Desjardins M, Lesage F. Wavelet-based estimation of the hemodynamic responses in diffuse optical imaging. Med Imag Anal. 2010;14:606–616. [PubMed]
  • Matteau-Pelletier C, Dehaes M, Lesage F, Lina J-M. 1/f noise in diffuse optical imaging and wavelet-based response estimation. IEEE Trans Med Imag. 2009;28:415–422. [PubMed]
  • Obrig H, Neufang M, Wenzel R, Kohl M, Steinbrink J, Einhaupl K, Villringer A. Spontaneous low frequency oscillations of cerebral hemodynamics and metabolism in human adults. Neuroimage. 2000;12 (6):623–639. [PubMed]
  • Obrig Y, Villringer A. Beyond the visible: imaging the human brain with light. J Cereb Blood Flow Metab. 2003;23:1–18. [PubMed]
  • Payne S, Selb J, Boas DA. Effects of autoregulation and CO2 reactivity on cerebral oxygen transport. Ann Biomed Eng. 2009;37(11):2288–2298. [PMC free article] [PubMed]
  • Prince S, Kolehmainen V, Kaipo JP, Franceschini MA, Boas DA, Arridge SR. Time-series estimation of biological factors in optical diffusion tomography. Phys Med Biol. 2003;48:1491–1504. [PubMed]
  • Rauch HE, Tung F, Striebel CT. Maximum likelihood estimates of linear dynamic systems. AIAA Journal. 1965;3 (8):1445–1450.
  • Roche-Labarbe N, Zaaimi B, Berquin P, Nehlig A, Grebe R, Wallois F. NIRS-measured oxy- and deoxyhemoglobin changes associated with EEG spike-and-wave discharges in children. Epilepsia. 2008;49 (11):1871–1880. [PubMed]
  • Saager RB, Berger AJ. Direct characterization and removal of interfering absorption trends in two-layer turbid media. J Opt Soc Am A. 2005;22 (9):1874–1882. [PubMed]
  • Saager RB, Berger AJ. Measurement of layer-like hemodynamic trends in scalp and cortex: implications for physiological baseline suppression in functional near-infrared spectroscopy. J Biomed Opt. 2008;13(3):034017. [PubMed]
  • Saager RB, Telleri NL, Berger AJ. Two-detector corrected near infrared spectroscopy (C-NIRS) detects hemodynamic activation responses more robustly than single-detector NIRS. NeuroImage. 2011;55 (4):1679–1685. [PubMed]
  • Tian F, Niu H, Khan B, Alexandrakis G, Behbehani K, Liu H. Enhanced functional brain imaging by using adaptive filtering and a depth compensation algorithm in Diffuse Optical Tomography. IEEE Trans Medical Imaging. 2011;30 (6):1239–1251. [PubMed]
  • Tong Y, deB Frederick B. Time lag dependent multimodal processing of concurrent fMRI and near-infrared spectroscopy (NIRS) data suggests a global circulatory origin for low-frequency oscillation signals in human brain. NeuroImage. 2010;53 (2):553–564. [PMC free article] [PubMed]
  • Tonorov V, Franceschini M, Filiaci M, Fantini S, Wolf M, Michalos A, Gratton E. Near-infrared study of fluctuations in cerebral hemodynamics during rest and motor stimulation: temporal analysis and spatial mapping. Med Phys. 2000;27(4):801–815. [PubMed]
  • Umeyama S, Yamada T. Monte carlo study of global interference cancellation by multidistance measurement of near-infrared spectroscopy. J Biomed Opt. 2009a;14(6):064025. [PubMed]
  • Umeyama S, Yamada T. New method of estimating wavelength-dependent optical path length ratios for oxy- and deoxyhemoglobin measurement using near-infrared spectroscopy. J Biomed Opt. 2009b;14(5):054038. [PubMed]
  • Villringer A, Hock C, Schleinkofer L, Dirnagl U. Near infrared spectroscopy (NIRS): A new tool to study hemodynamic changes during activation of brain function in human adults. Neurosci Lett. 1993;154:101–104. [PubMed]
  • Yamada T, Umeyama S, Matsuda K. Multidistance probe arrangement to eliminate artifacts in functional near-infrared spectroscopy. J Biomed Opt. 2009;16(4):06434. [PubMed]
  • Zhang Q, Brooks DH, Franceschini MA, Boas DA. Eigenvector-based spatial filtering for reduction of physiological interference in diffuse optical imaging. J Biomed Opt. 2005;10(1):011014. [PubMed]
  • Zhang Q, Brown EN, Strangman GE. Adaptative filtering for global interference cancellation and real-time recovery of evoked brain activity: a Monte Carlo simulation study. J of Biomed Opt. 2007a;12(4):044014. [PubMed]
  • Zhang Q, Brown EN, Strangman GE. Adaptative filtering to reduce global interference in evoked brain activity detection: a human subject case study. J of Biomed Opt. 2007b;12(6):064009. [PubMed]
  • Zhang Q, Strangman G, Ganis G. Adaptive filtering to reduce global interference in non-invasive nirs measures of brain activation: How well and when does it work? Neuroimage. 2009;45:788–794. [PMC free article] [PubMed]