Near infrared spectroscopy (NIRS) and diffuse optical imaging (DOI) are noninvasive technologies that can measure functional brain activity 1, 2, 3, 4, 5, 6 and have been used to study a variety of hemodynamic responses including vision, 7, 8, 9, 10 motor, 11, 12, 13, 14 and others, 15, 16, 17, 18 as well as fast neuronal signals.19, 20, 21 To detect evoked brain activity, photons must necessarily pass through superficial tissue layers (especially scalp and skull) when propagating to and from the cortex. Any hemodynamic variations in superficial layers will hence unavoidably affect the measured deoxyhemoglobin (HHb) and oxyhemoglobin concentrations. 12, 22, 23, 24 Cardiac activity is one of the major sources, and we also expect the hemodynamic signal to vary at the respiration frequency, due to phenomenon such as respiratory sinus arrhythmia,25 “chest pump,”26 and the respiratory wave in arterial blood pressure. 25, 27, 28, 29 Researchers have also categorized the spontaneous physiological low-frequency oscillations as low-frequency oscillations (LFO, vasomotor waves or Mayer waves), centered around , and very low frequency oscillations (VFLO), at about .30 In addition, the optical measurement of evoked brain responses is contaminated with the preceding physiological oscillations in the vasculature inside the brain. The spatial origins of these interferences (either from shallow layers or from inside the brain) are often not distinguished, and the interferences from different layers are normally lumped under the heading of “global interference” or “systemic physiological interference.” Several methods have been developed to suppress the interference in the detection of functional brain activity, including filtering, average wave form subtraction, and others. 13, 16, 19, 31, 32, 33, 34, 35
Previously, we have developed a novel way to remove global interferences based on a multiseparation probe configuration and adaptive filtering.36 In this method, an HHb (or ) measurement acquired with close source-detector separation (reflecting hemodynamic changes in the shallow layers) is supplied as a reference channel to an adaptive filter. The adaptive filter converts this reference to an estimate of the global interference; this estimated global interference is then subtracted from the target HHb (or ) from optical measurements using longer source-detector separations. Our Monte Carlo simulation suggests that signals from superficial layers actually comprise the major component of the total interference, thus making it an ideal reference measurement for adaptive filtering–based interference cancellation.36 Our Monte Carlo simulation results suggest that this method is very effective in suppressing such global interference and improving the contrast-to-noise ratio (CNR) in evoked brain activity detection. While simple, the algorithm is fast enough for real-time applications. However, its suitability and performance on human data remains to be evaluated. Here, we describe an example case of evoked visual response detection as a preliminary demonstration of our methodology.
Materials and Methods
The visual stimulation task was performed on a healthy, 36-year-old male, with no known history of neurological, psychiatric, or vision problems. NIRS measurements and simultaneous EKG and respiration recordings were collected during our study. Optical data was collected using NIRS1 (TechEn, Milford, Massachusetts), as used in previous studies.13, 37 We used two source lasers, at 682 and , respectively, and two detectors (Hamamatsu C5460-01). The lasers were intensity modulated by an approximately 5-kHz square wave and an in-phase/quadrature-phase (IQ) detection scheme. The final analog output was low-pass-filtered at using a second-order RC low-pass filter. An EKG was measured using chest leads, and respiration rate was collected using a piezo-electric respiration sensor (Pneumotrace II, UFI, Morro Bay, California) tied at the upper abdomen, and the measurements were digitally low-pass filtered at (eighth-order IIR filter). Heart rate (HR) was calculated from the EKG recordings, done on a beat-to-beat basis from QRS waves. First the derivative of the EKG signal is calculated, and then a threshold is applied to the EKG derivative so that QRS waves can be detected. The threshold is manually selected—it can be either for the rising or falling edge of the QRS wave to get the best result. Errors from this straightforward QRS detection will be corrected manually, and the calculated beat-to-beat HR will then be digitally low-pass filtered at (sixth-order IIR filter). The optical signals, together with the EKG and respiration measurements, were simultaneously sampled at and transferred to a computer using a 16-bit A/D conversion (Measurement Computing, PCM-DAS16D/16); thus, all signal channels were strictly co-registered to enable accurate comparison of timings of event onsets.
The lights from two diode lasers were combined using a bifurcated glass fiber bundle with a nominal core diameter (Fiberoptics Technology, prototype). Detector fibers were glass fiber bundles, also with core diameter. The optical fiber bundles (optodes) were attached to a flexible piece of plastic that was attached to a black Velcro headband that fastened snugly around the head of the subject. Soft black Velcro was also used to absorb stray light beneath the probe. The combined light of two wavelengths were in contact with tissue in one source location [S in Fig. 1a ], and exiting light was collected from two detector locations, and from the source, respectively [D1 and D2 in Fig. 1a].
During the visual task, the multichannel NIRS probe was positioned lateral and parallel to the subject’s midline, over the primary/secondary visual cortex (confirmed by an existing MRI scan), with the source inferior to the inion [Fig. 1a]. MRI images also show that the thickness of the scalp and skull at this location is approximately , and with this probe configuration, photon propagation theory suggests that measurement from S-D1 will be sensitive to scalp and skull changes, while measurements from S-D2 will be sensitive to scalp, skull, CSF, vision cortex, and some white matter regions. The subject sat upright in a comfortable chair in a dark and quiet room facing a laptop screen. The stimulus consisted of a counterphased, 100% contrast radial checkerboard ( luminance reversals) covering approximately of the visual field [Fig. 1b]. This was alternated with a mean-luminance, uniform field with a central fixation cross [Fig. 1c]. One run of 16 blocks was presented, where each block consisted of checkerboard oscillation for , followed by the fixation background for ; thus, the base visual stimulation frequency was .
The raw 200-Hz optical data were offset-corrected and further digitally low-pass filtered, in addition to the instrument filter, using a 1.25-Hz, sixth-order Butterworth filter. The optical measurements for each channel were converted to relative changes in the concentration of HHb and using the modified Beer-Lambert law. 38, 39, 40, 41 Concentration measurements were then bandpass filtered using a sixth-order Butterworth low-pass filter at in series with a sixth order Butterworth high-pass filter at (to remove any slowly drifting signal components) and then down-sampled to . Based on our previous Monte Carlo simulations with common head structure and tissue optical properties,36 we chose DPFs of 5.4 and 7.0 for S-D1 and S-D2 at , and 5.1 and 6.6 for S-D1 and S-D2 at . (We initially set all DPFs to 6.0 for the data analysis, and the results closely mirrored the results reported here, suggesting insensitivity to DPF errors.36 The reported DPFs are based on more realistic Monte Carlo estimates of DPF differences across wavelengths and separations.) The time series of concentration changes acquired from both source-detector pairs (S-D1 and S-D2) were fed into the adaptive filter. (HHb was filtered similarly but separately.) The adaptive filter used a finite impulse response (FIR) and transversal structure (tapped delay line), with 100 nodes (“weights”). The Widrow-Hoff least-mean-squared (LMS) adaptation algorithm was used to optimize the filter coefficients on a sample-by-sample basis. The details can be found in our previous report.36 The convergence control parameter was . In order to speed the convergence, we prenormalized the target and the reference data sets by dividing by their estimated standard deviation so that both time series have a standard deviation of approximately 1. After the adaptive filtering, the quantitative values of the filtered results were recovered by multiplying back the standard deviation. The initial guess of the adaptive filter weights was —in other words, we assume that the shape of the global interference and the S-D1 measurements are identical at the beginning of the adaptive filtering. If this initial guess were to remain constant, it would be similar to Saager and Berger’s approach, which assumes that the entire S-D1 time series can be least-squares fit to and subtracted from the S-D2 measurement.35 Our approach differs in that (1) we assume a linear mapping between S-D1 and S-D2 and that this mapping is dependent on recent history (here, 100 weights instead of the entire time series, making the approach particularly suited to real-time applications), and (2) that this mapping is allowed to change over time, which in principle allows us to handle nonstationary time series (see Sec. 4.2).
As in previous NIRS studies, we performed block averaging using nonoverlapping portions of the data set triggered on stimulus onsets (averaging windows prior to checkerboard onset and after the onset of each block of visual stimulation).
After adaptive filtering, in order to retrieve the quantitative functional hemodynamic changes in the deep cortex layer in our previous Monte Carlo simulation study,36 we calculated a sensitivity correction scaling factor. This is done by introducing a known perturbation into the blood content of the simulated gray matter layer only and then calculating the blood content changes using the surface measurements and the MBLL. The sensitivity correction factor is defined as the ratio of true blood concentration change to measured blood volume change. In the simulation with general physiological parameters, this sensitivity correction factor is calculated to be 27.5 for and 26.4 for HHb, and we assumed the same correction factor in this human study.
Power spectral density (PSD) was used to identify the components and vision stimulation–associated changes in the hemodynamic responses or other physiological signals. PSD was estimated in MATLAB (based on Welch’s averaged, modified periodogram method, with a Hanning window for smoothing). In order to find peaks that are statistically significant, we chose a confidence interval of 0.95 and calculated the upper and lower bound of the power spectrum within this interval. Usually, we say a peak can be identified if the lower bound of the peak is above the averaged surrounding noise floor shown in the PSD.
With PSD, we can also perform a quantitative CNR analysis. For human subject data, we cannot totally separate signal and noise. Thus, we roughly calculate the signal power by integrating the power (in linear scale) from the PSD over the “signal bands,” which are the small frequency segments surrounding the stimulation frequency and its second harmonic. This way, the evoked brain activity peak will normally be included. We integrate the rest of the power spectrum (“noise bands”) to get the noise power. The unavoidable error in this calculation includes real signal power in the noise band and actual noise power in the signal bands. The CNR is then calculated as the square root of the ratio of signal power to noise power. It is worth mentioning that since the PSD is calculated using the whole data set (9674 points at sampling rate, no overlap in the PSD segmentation), the CNR calculated based on this PSD is the result of averaging the whole data set, not the CNR for a single stimulation.
Origins of Global Interference: Comparing Optical, Respiration, and Cardiac Measurements
The respirometry and EKG were recorded simultaneously with optical measurements in order to help us determine the origins and understand the nature of the global interferences. The raw optical measurements (light intensity) at collected from the “far” detector (S-D2 with source-detector separation) and from the “near” detector (S-D1 with source-detector separation) are shown in Figs. 2a and 2c , respectively. In order to more clearly view the details, such as the temporal correlation between different channels, we have shown here only of data. The associated power spectral densities are shown in Figs. 2b and 2d. The simultaneously recorded respiration signal (left axis with arbitrary unit, indicating the movement of the upper abdomen) and heart rate (right axis, calculated from the EKG) are presented in Fig. 2e, and the associated PSD is shown in Fig. 2f. The insets in Figs. 2b, 2d, 2f show the details of PSD at frequencies less than . Each PSD was estimated using all 96,736 data points so that more details can be seen at the low-frequency end. To avoid possible confusion, we need to point out again that the lower curve in Fig. 2e is the heart rate, not the EKG, and its spectrum in Fig. 2f is actually the spectrum for heart rate variation (HRV).
In Figs. 2a and 2c, we see that both optical measurements consist of a slow variation and a fast oscillation. The fast oscillation is obviously due to cardiac activity, it co-registers strictly with the EKG (not shown), and the inset in Fig. 2a shows an enlarged data segment to demonstrate the temporal details of these fast signal fluctuations. The PSD of the optical signals, shown in Figs. 2b and 2d, show that the cardiac activity peaks at , and its second harmonic at about is also very obvious.
From the respiration PSD [the inset in Fig. 2f], we see that the principal respiration frequency is at . However, for both optical measurements [seen in the insets in Figs. 2b and 2d], the major low-frequency peak is (with a peak appearing about down, and among peaks that might be artifacts). The heart rate PSD [inset in Fig. 2f] helps in explaining the origin of this peak in optical channels. The heart rate PSD shows and as the top two principal peaks, with the peak as the highest and above the peak. (The latter is obviously due to the respiratory sinus arrhythmia.) This suggests that the oscillation in the optical channels is physiological rather than artifact. Since is the prototypical LFO frequency, we suspect that this oscillation in the hemodynamics shown in the optical channels is vasomotor waves, and the peak in the heart rate PSD is possibly an associated change in a vasomotor reflex. Full separation of component signals is difficult in this subject, because the subject’s respiration rate is relatively low, and we lack other more direct measurements such as invasive arterial pressure.
Optical measurements from far and near channels are very similar, actually a linear correlation analysis that shows that they are highly correlated ( , ). The fact that the target (far) measurement and reference (near) measurement are highly correlated and that they all demonstrate major PSD peaks at heart beat and LFO frequencies indicates that the optical measurement from S-D2 is heavily contaminated by global interference (measured by S-D1), as was previously derived via Monte Carlo experiments.
We performed a T-test on the 16 averaged heart rates during vision stimulation periods and the 16 averaged heart rates during rest periods . The change was not significant [ , ], suggesting no obvious heart rate increase or associated hemodynamic interference due to the visual stimulation itself.
To summarize, the PSD in the optical and physiological measurements suggests that the systemic interference comes from mixed sources, with cardiac activity and LFO (Mayer waves or vasomotor waves, judging from the frequency) being dominate. The visual stimulation does not cause any heart rate increase.
Changes during Visual Stimulation
Changes in during visual stimulation are shown in Fig. 3 . Here, Figs. 3a and 3b represent the target measurement of calculated from the raw “far” signal and its block average. Figures 3c and 3d show the reference (“near” signal) and its block average. Although the target data set was expected to contain an increase in following stimulus onset (and concomitant decrease following stimulus offset), neither the raw time series nor the block average show any obvious expected signal change. We hypothesized that the global interferences obscured the expected brain-related change. The fact that the signal variations in the target closely match those of the reference (which should contain no visual response) suggested that indeed global interference may dominate the target data set.
The target data set after adaptive filtering is shown in Figs. 3g, 3h, 3i, 3j. Figures 3g and 3h show the raw adaptive filtered result with its block average, and 3(i) and 3(j) show the filtered result following sensitivity correction. Comparing Figs. 3g and 3h [or their sensitivity-corrected versions in 3i and 3j] with 3a and 3b, we see that after adaptive filtering, both LFO and cardiac interference are substantially reduced (actually, 80% of the signal variation was removed), while at the same time, the CNR of the hemodynamic responses to vision stimulation increased dramatically. Both the time series’ and the block averaged results show an expected increase following stimulation onset and return to baseline during rest periods, with temporal changes appropriately associated with the stimulation paradigm. For the detected visual response shown in Fig. 3j, we see a uniphase increase starting from approximately after stimulus onset, which plateaus at about . The response remains steady until after the visual stimulation ends (we suspect that the hump at about is artifact), and begins to decrease at about , recovering to baseline at about . Although unlike in our previous Monte Carlo simulation study, here we do not have a true visual response value to compare with, the quantitative changes match previous reports.9, 10
We also compared the adaptive filtering result with the traditional low-pass filtering result. acquired from the far detector [S-D2, Fig. 3a] is low-pass filtered with an eighth-order Butterworth filter with bandwidth, and the result is shown in Figs. 3e and 3f. From Fig. 3e, we can see that this low-pass filter effectively removes the cardiac oscillations; however, systemic variations due to LFO remain. After block averaging, seen in Fig. 3f, the systemic interference was further suppressed; however, it is still too large to reveal clear evoked brain activity. Comparing with the final adaptive filtering block averaged result [Fig. 3j], we can see that adaptive filtering removes the global interference and recovers the evoked brain activity much more effectively.
The CNR improvement from adaptive filtering can be quantitatively estimated using power spectral density techniques. Figure 4a shows the PSD of the target before adaptive filtering, in which we see two dominating peaks at about and , which correspond to cardiac activity and LFO. Figure 4b shows the details of the PSD of target in 0 to range (solid line), together with its upper and lower bound calculated at a confidence interval of 0.95. In Fig. 4b, using the average local noise level (horizontal dashed line; averaged PSD magnitude in the 0 to range) as the detection threshold (to compare with the lower bound of the PSD), we can see several peaks above this threshold (e.g., at , , , , , and ), the LFO peak being the most significant one, about above. A peak of , the expected visual stimulation frequency, is actually among detected peaks; however, it is the second smallest compared with other detected peaks, only above the threshold, and difficult to distinguish among the other detectable peaks as real or artifact. After adaptive filtering [Fig. 4c], the physiological interference due to LFO and cardiac activity peaks was significantly reduced. From the details presented in Fig. 4d, we can see that the noise floor (horizontal dashed line) was less than the original, while the visual response peak remains [compared with Fig. 4b]. Now only two peaks can be detected ( and ), and the peak at , corresponding to vision stimulation, is the most significant peak and is above the local noise floor.
If we choose the signal bands of the vision response to be roughly 0.028 to (base band) and 0.061 to (second harmonic), and the rest of the spectrum in the 0 to range is the noise bands, our CNR analysis shown that adaptive filtering doubled the CNR (from 40.2% to 80.8%). The very low frequency of was also preserved. It remains to be determined whether signal variation at this frequency range is real or artifact.
HHb Changes during Vision Stimulation
One obvious difference in the HHb data as compared to is that for HHb, the interference was relatively small and the hemodynamic response was visible before adaptive filtering [Figs. 5a and 5b]. In addition, the signal from the near detector (shallow layer hemodynamics) differs more substantially from the far detector signal (deep-layer hemodynamics)—that is, the far detector HHb signal was not dominated by synchronized global interference. We suspect that interferences here have local origins and thus are different at different detector locations. This explanation is supported by the PSD analysis. In the PSD of target HHb measurement [Figs. 6a and 6b], we actually see no obvious LFO, respiration, or cardiac activity peaks. In other words, the major global interference sources in HHb are not significant. In Sec. 5, we discuss why in this subject HHb did not demonstrate obvious global interference peaks, as previously observed in .
Since the amount of interference in the target data set is relatively small and the visual response and the global interference are mostly independent of each other, the adaptive filtering did not substantially alter the signal. As seen in Figs. 5g and 5h, the overall amplitude of variations after adaptive filtering have only slightly changed (17% smaller), unlike the adaptive filtering results, where 80% of the variations in the target data set were removed [Figs. 4g and 4h]. Comparing the results before and after adaptive filtering [Figs. 5a with 5g, respectively], we see that the visual response is clear in both, although adaptive filtering appears to slightly reduce the signal and increase noise. This can be more clearly seen in the PSD analysis in Figs. 5b and 5d, where we see a reduced peak at , and a elevated noise level. The previous CNR analysis shows a 35% reduction of CNR after adaptive filtering.
When comparing Figs. 5i and 5j with the low-pass filtering result shown in Figs. 5e and 5f (which use an eighth-order Butterworth filter with bandwidth), we can see that the low-pass filtering result is smoother. The shape of the response in 5(j) (onset, rise, offset, and fall times) has less low frequency fluctuation and is slightly more temporally consistent with prototypical hemodynamic responses than that in Fig. 5f, but it is unclear if this is because of adaptive filtering.
For the filtered and sensitivity-corrected visual response shown in the block averaged and sensitivity corrected results in Fig. 5j, we see a uniphase decrease starting from approximately after the onset of the stimuli, which continues until about , when it reaches a valley. The signal remains at this value until after the visual stimulation ends, and it begins to increase at about , again recovering to baseline at about . These quantitative HHb changes due to vision stimulation also match those found in previous reports.9, 10
The improvement from adaptive filtering on HHb in this case is not obvious. Although the shape of the filtered results [Fig. 5j] seems to agree more with expectations compared with Fig. 5b without adaptive filtering, it is unclear if this slight improvement is really from adaptive filtering, and the CNR estimated using PSD is actually reduced (from 66.5% to 43%). This suggests that in the cases where hemodynamic signals demonstrate obvious functional brain response and its PSD demonstrate no obvious global interference sources, adaptive filtering may be unnecessary.
Comparison of Shallow Layer Hemodynamics and Estimated Global Interference
The performance of our method depends on whether the shallow layer hemodynamics acquired from S-D1 and global interference in the visual response detection acquired from S-D2 are correlated. When they are well correlated, the global interference can be optimally canceled, and the CNR of the evoked hemodynamic response is improved. A further question is whether the relationship is static or dynamic. Since it is difficult to precisely separate global interference from the raw S-D2 blood content time series, and since the visual response we acquired after adaptive filtering seems quite reasonable, we chose to estimate the global interference by subtracting the visual response [Fig. 3g] from the raw time series with both interference and visual response embedded [Fig. 3a] and use the residual as global interference. Although it is not mathematically strict (the “estimated global interference” is actually the component in the target signal that correlates most with the reference signal, picked up by the adaptive filter), the fact that we subtract a decent visual response makes this estimation acceptable, since generally the target signal is the combination of visual response and interference.
In Fig. 7 , we compare our two raw concentration time series directly by mean-subtraction and dividing each time series by its own (temporal) standard deviation. Only a small segment is shown for clarity. Here, the normalized concentration acquired from the near measurement (S-D1) is shown as the dashed line and the computed global interference [Fig. 3a minus Fig. 3g] is the solid line. As shown, the two time series are highly overlapping. Actually, a linear regression of the two measurements [Fig. 8a ] shows that the correlation coefficient between these two time series is 0.96 . Although the estimated global interference may be different from real global interference, this correlation comparison demonstrates the similarity of the reference signal to the target signal after the visual response is removed and in turn explains why adaptive filtering removes the interference term very effectively.
Our previous analysis suggests that there is reasonably a linear relationship between the shallow layer hemodynamics and global interference; however, this relationship changes with time. In Fig. 8b, one segment of global interference is plotted against the shallow layer hemodynamics (from to ; both are smoothed using a sixth-order Butterworth filter), and the arrow shows how the data moves with time. The data did not strictly follow the line. Thus, we performed short-term linear regressions—that is, we performed linear fits using only 100 data points at a time and slid the data window through the entire data set, with the estimated global interference as the independent time series and the near measurement as the dependent variable. The resulting linear fit slope and intercept values oscillate over time (Fig. 9 ). This dynamic feature can also be viewed in the adaptive filtering, where the node is adjusted at a point-by-point basis and is following the slow nonstationary variations of the time series to optimize the interference cancellation. The update of the first node as a function of time is presented in Fig. 10 .
In this case study, the fact that the shallow layer hemodynamic changes correlated well with the global interference and that adaptive filtering in general follows the nonstationary changes in the time series explains why our method should work well in suppression of the global interference and improve the CNR in the evoked brain activity detection.
This case study provides a preliminary demonstration of the potential effectiveness of our combined adaptive filter/multidistance methodology for global interference reduction. From this study, we see the effect of adaptive filtering both in the situation when global interference dominates (e.g., in ) and when the global interference does not dominate (e.g., in HHb). Clearly, adaptive filtering is substantially more useful in the former case. The key for adaptive filtering to work is that the interference must be common to both the signal channel and the reference channel. In principle, cardiac activity, respiration, and vasomotor waves should be common in the whole circulation system; however, because of the heterogeneity of the vasculature in tissue and the location of the probe (e.g., relative to large vessels), together with local changes at each optode, different optical measurements may be affected by different types of interference to a different degree, and not all optical measurements are dominated by the same type of interferences. To maximize the benefit of adaptive filtering, one must first determine whether a given signal is substantially affected by global interference common to both the signal channel and the reference channel (e.g., via PSD comparison). If the target channel shows interference similar to the reference channel, then apply the adaptive filter. Otherwise, skip adaptive filtering.
It was interesting to observe that HHb (found mostly in venous compartments) exhibited little global interference, whereas the global interference in (present in all compartments) was substantial. This makes sense in light of signal frequencies, which indicate that the interference was largely vasomotor or Mayer waves, which are arterial pressure waves.25 Since this is only a case report and the goal is to demonstrate the methodology, no physiological generalizations should be drawn . For example, the subject in this case report showed an HHb response but not an response after block averaging. We believe this occurred because of excessive physiological interference of the overlying layers and predominantly in , a problem which was substantially resolved by adaptive filtering. In other studies, group averages normally show and HHb responses, with a certain amount of variability. Results from individual subjects, however, are even more variable, sometimes showing clear responses in both and HHb, sometimes only in one Hb species, and sometimes in neither.
The sensitivity correction factor is a constant that converts the evoked brain hemodynamics acquired from surface measurements and MBLL to quantitative blood content changes in the cortex. Based on the photon transport theory, this constant changes with parameters such as the probe configuration and dimensions and the optical properties of different layers of the head. Since not all this information is known in our test, errors in the assumed sensitivity correction number in our study are unavoidable. However, since the focus of our study is to demonstrate the effectiveness of adaptive filtering by looking at the contrast-to-noise changes, and this constant changes the quantitative results of the and but not the CNR (this sensitivity correction factor is canceled in the CNR calculation), the conclusions of this case study should not be affected by the error in the sensitivity correction factor.
The overall performance of this adaptive filtering approach to removing global interference during in vivo applications depends on whether the shallow layer hemodynamics provide a good estimate of the global interference in the target data set. According to photon transport theory, optical measurements are more sensitive to superficial layers, and we have demonstrated using Monte Carlo simulation that hemodynamic changes in the superficial layers is the major component of the global interference in evoked brain activity detection. Furthermore, since the blood vessels for the superficial and deep layers are connected and can be closely traced back to the same origin, namely, the common carotid arteries, one might expect hemodynamics in these layers to be reasonably correlated. Thus, it would seem to be a reasonably general principle that the scalp and skull hemodynamics acquired using optical measurements with short source-detector separation should provide a good estimate of the global interference in measurements of deeper layers. The observed correspondence between shallow layer hemodynamics and global interference in this case study explains why this can work on human subjects when global interference dominates.
In addition, since living human tissue is a very nonstationary system, we expect the relationship between the shallow layer hemodynamics and the global interference to change with time, as observed in this case. The adaptive filtering method assumes that there exists a linear mapping between the reference and the target data set. This mapping is allowed to change, but the change needs to be “slow enough” for adaptive filtering to follow. This slowness depends on the convergence of the adaptive filtering method and the magnitude of the mapping change. We chose one of the simplest available adaptive filtering methods and parameter optimization approaches (LMS). Many other methods could be applied, with better convergence rates, while maintaining good stability and robustness. These approaches remain to be explored, as here we sought only to demonstrate the suitability of the approach for real physiological data.
In many applications, it is desirable to have interferences removed in real time. In some applications, including functional brain area localization and biofeedback, this real-time feature is mandatory. The adaptive filtering method we used is computationally trivial and hence can be easily implemented in real-time manner. Based on the extensive knowledge base surrounding adaptive filtering techniques, several techniques can help improve the convergence of the adaptive filter. First, since the convergence and stability of the LMS algorithm depends on the energy of input signals, normalization (so that both target time series and reference time series have roughly a standard deviation of one) increases the convergence rate and standardization of the choice of initial value of the adaptive filter and the step control parameter. Second, a pretest is needed to train the adaptive filter before the real brain activity detection can be utilized to generate a better initial guess at node weights. From the pretest, we can estimate the amplitude of the target measurement and reference measurement, and these values can be used in the normalization of the two time series. It would also allow the adaptive filter to reach a reasonable value, to be used as the initial filter coefficient for the real test. When really high filtering speed is required, such as in some real-time imaging of brain activity, one may even consider the normalization-subtraction method described in a previous report (Sec. 2, Eq. (5) in Ref. 36), where the amount of calculation is minimized.
In Fig. 3, before filtering, the quantitative amplitude of the concentration changes at the far separation [Figs. 3a and 3b, from S-D2] is only about half of that acquired from the near separation [Figs. 3c and 3d, from S-D1]. This might be due to a partial volume effect and suggests that the spatial origins of the major physiological variations may be limited to a thin, vascularized layer (e.g., the dermis layer of skin). As the sampling volume photon probes increase with source-detector separation, the percentage of area with hemodynamic variation over the whole sampling volume decreases, and thus the amount of variations acquired using MBLL is reduced. In principle, one could compensate for such a layer using sensitivity correction, although this would require knowledge of the location and thickness of that vascularized layer. The result in Fig. 3 also indicates that in our experiment, a direct subtraction of reference signal [Fig. 3c] from the target signal [Fig. 3a] will not work for the purpose of removal of global interference, because quantitatively they are not comparable (due to the difference in the average blood concentration in different tissue area, or partial volume effect), and actually in our case, such subtraction will generate a negative result since S-D1 detects a much larger variation than S-D2.
The fact that our adaptive filtering dramatically removes global interference and improves CNR in in the evoked brain hemodynamic response detection is exciting, and PSD analysis should be practical in helping to judge whether the optical measurements are global interference–dominated and adaptive filtering should be applied. However, the approach needs to be tested in different experimental settings and on larger sample sizes to assess the breadth of the method utility. For example, in this case study, the global interference and the evoked brain activity are not correlated (suggested by the T-test of heart rates at the rest and stimulation periods), and further study is need to explore the results in the situation when the two are partly correlated. The case study described in this paper sets the data analysis platform for future human subject tests.
This work was supported by the NIH (Grant Nos. K25-NS46554, R21-EB02416, R01-DA015644, and R01-EB006589). We would also like to thank Margaret Duff for her help preparing the manuscript.