1.

Introduction

Diffuse optical imaging (DOI) is an effective technique for the non-invasive monitoring of cerebral hemodynamics and oxygenation.^{1 2 3} Compared to functional magnetic resonance imaging (fMRI), it has several advantages, including high temporal resolution, safety, and low cost. In addition, DOI can detect local changes in the concentrations of oxyhemoglobin (Δ[HbO]) and deoxyhemoglobin (Δ[Hb]) in tissue based on differential absorption at multiple wavelengths.⁴ The most important trade-off is spatial resolution, which is inferior to that achieved by fMRI. An additional problem with DOI is that, in addition to being sensitive to the evoked hemodynamic changes, it is also sensitive to cerebral hemodynamic fluctuations of systemic origin. Such systemic fluctuations are associated, for instance, with arterial pulse, respiration, and heart rate fluctuations,^{5 6 7} which appear as interfering signals, obscuring or even overwhelming the activation-evoked response.

What is more, since heart rate may increase in response to the same stimuli that are presented to evoke neuronal response, especially during motor activity, these systemic hemodynamic effects may exhibit a temporal pattern of amplitude change that is highly correlated with the evoked vascular activity of interest. Hemodynamic changes due to heart rate increase may appear over a large spatial area,⁸ making it more difficult to recover the localized hemodynamic changes evoked by stimulation.

In particular, currently identified sources of physiological interference include

Previous attempts to treat this problem used temporal filtering or modeled the interference signals with predefined basis functions. In Ref. 18, for multi-trial experiments, the time courses of the measurements were bandpass filtered to correct for slow drifts and to eliminate arterial pulsations, and then averaged across trials to improve signal-to-noise ratio (SNR) for each channel (source-detector pair). Also, to estimate and correct for the effect of systemic hemodynamics, “inactivated” pixels were manually identified and their hemoglobin concentration changes were assigned to systemic hemodynamics and subtracted from the measurements at all other pixels. In Ref. 19, the temporal measurement of optical intensity was modeled as the sum of three quasi-sinusoidal signals with known frequency but slowly varying random amplitudes and phases. The amplitude and phase parameters became coefficients to be reconstructed. The model frequencies were chosen to correspond to heart rate, respiration rate, and the block cycle rate of the experimental independent variable (e.g., finger tapping), all calculated from auxiliary signal measurements. The amplitude and phase of the sinusoids were then estimated using a fixed-interval extended Kalman smoothing algorithm.

It was shown in Refs. 5, 6, and 7 that during finger tapping, systemic rhythms such as heartbeat, respiration, and heart rate fluctuations exhibited a significant degree of synchronization with the periodic stimulation/rest sequence. This phenomenon, together with the global characteristics of systematic interference, suggest that temporal filtering or the use of predefined temporal basis functions may not take full advantage of available information to separate the evoked cerebral hemodynamics from physiological interferences. In particular, one might expect to be able to exploit “spatial diversity,” that is, the difference in spatial behavior between local desired signals and global “interference” signals, to improve suppression of the systemic effects.

The hypothesis of this work then is that the spatial behavior of the interference across the sensor array can be used to distinguish it from that of the activation. We have observed that the spatial modes of the systemic response will be much broader and of higher amplitude than that of the activation signals, which will in contrast be more focal and have less total energy. We propose to exploit that difference to achieve interference reduction. In particular, we make a mathematical assumption of orthogonality between the spatial interference subspace and the spatial evoked activation subspace.

We note that spatial processing should be essentially independent of and complementary to temporal filtering or modeling. In addition, the difference in spatial behavior should hold true whether the interference is temporally uncorrelated or temporally correlated (e.g., time locked) with the activation signal via heart rate variation or other systemic stimulus-dependent physiological responses. Since the amplitude of the interference is large compared to the activation signal, sometimes even completely obscuring the activation, adequate interference suppression is essential to be able to exploit the advantages of diffuse optical imaging, especially in the case of weak activation signals, and there is a need to take advantage of these complementary approaches to interference reduction.

In array signal processing, a now classic approach to attenuate interference in sensor array data is to exploit the eigenstructure of the data spatial correlation matrix.²⁰ In the problem at hand, since the interference signals are assumed to be much larger than the desired signals in both total energy and spatial spread, we expect that the interference will be the dominant spatial pattern in baseline (i.e., without stimulus) data. We find these spatial patterns via an eigendecomposition of a particular spatial correlation matrix, and then project the stimulus data onto the orthogonal complement of this dominant baseline subspace.

This implied orthogonality assumption is not based on any specific model of physiological processes. But we hypothesize that it is reasonable due to the insight that the spatial behavior is likely to be an independent distinguishing characteristic with reference to temporal behavior to discriminate interference from evoked response, while temporal behavior of all signals is quite similar. On the other hand, this orthogonality condition is not absolutely necessary, because our purpose is to remove the strong interference and get a more localized activation area. Thus, as long as the signal component in the subspace we project onto is larger than the interference component in that subspace, with respect to their respective power in the measurements, we can hope to benefit.

The reason that we estimate the interference subspace from baseline data, rather than directly from the stimulus data, is that during stimulus, the physiological interference can be time-locked with the frequency of stimulation, making it hard to separate the interference subspace from the induced activation subspace. The baseline data, by contrast, contain the spatial patterns of the interference only, which should be the same in both baseline and stimulus data. Thus, we are using a principle components approach to identify an interference subspace from baseline data and then attenuate the interference in stimulus data.

In the next section, we describe our experimental and mathematical methods, and in the following section we present results of testing the method on data acquired from two human subjects. A discussion section and a section presenting our conclusions follow.

2.

Methods

2.1.

Data Acquisition

The data used in this work were previously published in Refs. 8 and 18. In these experiments, a study of the hemodynamic evoked response to a finger opposition task, passive finger tactile stimulation, and electrical median nerve stimulation in the human sensorimotor cortex were performed. Only data from experiments using the first two stimuli were considered in this work.

The instrument used was a multichannel continuous wave (CW) optical imager (CW4 system), developed at Martinos Center, Massachusetts General Hospital, which employs 16 laser diodes and 16 avalanche photodiode detectors (Hamamatsu C5460-01). Eight of the laser diodes emit light at a wavelength of 690 nm and eight at 830 nm. Each source location has sources at both wavelengths, for a total of eight source positions. Sources and detectors are fiber-coupled and placed on the head according to the arrangement shown schematically in Fig. 1. The minimum source-source and source-detector separations are 1.9 and 3.0 cm, respectively. The area imaged on each hemisphere is 5.6×6.0 cm. The signal from each detector was sampled at 25 Hz.¹⁸

Figure 1

Geometrical arrangement of the source and detector fibers on the head.

For the finger opposition task experiment, participants were instructed to touch their thumb with the index and middle finger following a 4 to 5-Hz metronome beep. In the passive stimulation experiment, tactile stimulation was delivered manually by an investigator by touching the first three fingers of the participant at the same rate as in the finger opposition task.¹⁸ The block protocol alternated stimulation periods (20 s long) and rest periods (20 s long). A complete session consisted of an initial reference rest period, in which baseline data were recorded, followed by a block design run of 10 task/rest sequences. Other details about the experiments can be found in Ref. 18. For each subject, the experiment was repeated for four runs by stimulating left and right hands with both finger tapping and tactile stimulation. In this paper, spatial eigenfiltering was performed on measurements from two subjects (denoted A and B).

3.

Results

We used the 48 channels that had either the minimum 3-cm (28/48) or next shortest 4-cm (20/48) source-detector separation, as shown in Fig. 2. Channels with separation more than 4 cm had low SNR values and the measurements were dominated by noise.

Figure 2

Forty-eight channels selected for the data matrix; dotted lines show the 4-cm channels and solid lines the 3-cm channels.

Figure 3 shows the block-averaged time courses of Δ[HbO] response to the left hand finger-tapping stimuli, for subject A. For clarity, for all time course figures in this paper, we plot only the 28 channels with 3-cm separation. In each subplot, the original response (dotted line) and the result after filtering using the baseline interference subspace (thick solid line) are presented together. The horizontal bar on each panel indicates the duration of the 20-s stimulus period. Each subplot is labeled by its source-detector pair, numbered as shown in Figs. 1 and 2. Note that the upper two rows correspond to the channels on right hemisphere (contralateral side for left-hand stimulus), while the lower two rows show the channels on the left hemisphere (ipsilateral side).

Figure 3

Block-averaged Δ[HbO] response in a left hand finger-tapping stimulus of subject A: original response (dotted line) and baseline filtered result (thick solid line). The horizontal bar indicates duration of the 20-s stimulus period. The title of each subplot shows the source-detector pair and the time axis ranges from −10 to 30 s.

For this result and all other results for subject A, two baseline eigenvectors were projected out. The number of eigenvectors to be removed was chosen by observing the spatial pattern of the eigenvectors, which becomes clearer in the following when we present spatial maps of the eigenvectors.

One can see in Fig. 3 that the original data showed a response spread over a large area on the contralateral hemisphere, presumably due to the global systemic interference arising from a stimulus-induced heart rate acceleration.⁵ On the ipsilateral side, although the response had a smaller magnitude, many channels still had a similar shape response pattern. Spatial filtering with the eigenvectors produced by the baseline data significantly reduced the global interference, leaving the remaining response on a more localized area (the channels associated with detectors 6 and 7). Also the filtered results show a deactivation response, i.e., decrease of Δ[HbO] for the ipsilateral channels.

Figure 4 shows back-projected spatial maps of the first three eigenvectors of both baseline and stimulus data sets. The maps are shown in pairs. The left-hand member of each pair corresponds to the imaged area on the left hemisphere and correspondingly for the right-hand member. The blank column in the middle of each pair of panels schematically indicates the separation area between two hemispheres. The percentage given above each pair of panels indicates the proportion of the sum over all eigenvalues that was accounted for by that particular eigenvalue. All the plots have the same color map, ranging from −0.6 to +0.6.

Figure 4

Back-projected spatial maps of the first three spatial eigenvectors of both baseline data and the stimulus data matrix of the left-hand finger-tapping stimulus for subject A. For each eigenvector, the left map corresponds to the imaged area on the left hemisphere and the right map to the right hemisphere. All the plots have the same color map, ranging from −0.6 to +0.6.

One can see that for the baseline data, the first three eigenvectors exhibit both a rather homogeneous distribution, presumably corresponding to the systemic physiological interference, plus some peak values located at noisy channels. Almost all these noisy channels have 4-cm separation. For the stimulus data, the first eigenvector shows a more homogeneous distribution than the baseline data does, presumably corresponding to the summation of systemic physiological interference and evoked activation, which are time-locked to the stimulus, and the corresponding eigenvalue is a much higher proportion of the total compared to the baseline data. These observations hold for all cases studied. The second and third eigenvectors, like their counterpart in the baseline data, result in spatial maps of “noisy” channels, i.e., high magnitudes at noisy channels but small magnitudes elsewhere, and the corresponding eigenvalues take much lower energy.

Figure 5 shows the block-averaged time courses of the Δ[HbO] response to the right-hand finger-tapping stimuli for subject A. The postfiltered response is localized around channel S8-D12, and to some extent S6-D12. Note that before filtering, not only do we have a broad spatial response, but these two channels showed a smaller response than even several ipsilateral channels.

Figure 5

Block-averaged Δ[HbO] response in a right-hand finger-tapping stimulus of subject A.

Looking at Figs. 3 and 5, one observes a spatial asymmetry in the tapping stimulation responses, for both left and right hands. This observation, which also is shown in the response to tactile stimulus, is due to the practical difficulty of placing the probes perfectly symmetrically on the two hemispheres.

Figure 6 shows the Δ[HbO] block-averaged time courses for left-hand passive tactile stimulus for subject A. Similarly to the preceding results, the eigenfiltered response has a more localized response area on the contralateral side, associated with detectors 3 and 8. Again, prefiltered ipsilateral responses of similar magnitude to contralateral responses (see S6-D11, for example) were significantly attenuated by the eigenfiltering. In contrast, the indicated contralateral channels are hardly attenuated at all.

Figure 6

Block-averaged Δ[HbO] response in a left-hand passive finger tactile stimulus of subject A.

For subject B, we present the result of filtering the right-hand tapping stimulus data in Fig. 7. For this result and all other results for subject B, three baseline eigenvectors were projected out. The “double-peak” interference temporal pattern is global and very obvious at most channels. Projection onto the baseline subspace gave a more localized response area focused around channels S6-D14 and S6-D15. Moreover, the very noisy channel S1-D1 is effectively attenuated by eigenfiltering.

Figure 7

Block-averaged Δ[HbO] response in a right-hand tapping stimulus experiment of subject B.

Figures 8 and 9 show back-projected spatial maps of CNR values for all data sets in this study. Results with data from subject A are in Fig. 8 and from subject B in Fig. 9. For clarity, any negative CNR values, which are the result of an apparent deactivation response on some ipsilateral channels, were set to zero. (See Sec. 4 for some comments on this “deactivation” result.) Note also that the amplitude range of the CNR maps vary from map to map, and in general the amplitude range is lower for the filtered results simply due to the decrease in total energy from the filtering. Note that the filtered maps show a much more concentrated region of large CNR in comparison to the CNR maps before filtering. This change is generally more notable in the finger-tapping experiments than it is in the passive stimulation experiments. Results are comparable across both subjects.

Figure 8

Spatial maps of the CNR values of subject A: lhtap, left-hand tapping; rhtap, right hand tapping; lhpas, left-hand passive; and rhpas, right-hand passive. The original response is shown on the left and the baseline filtered results on the right. The arrangement of each pair of maps is as in Fig. 4.

Figure 9

Spatial maps of the CNR values of subject B.

Figures 10 and 11 show the back-projected spatial maps of CC values of all experiments for subjects A and B. We use different color maps for the four original data sets for subject A, but a common color map for subject B, for clarity. For all spatially filtered results for both subjects, the color maps range on [−1,1]. The CC maps of the original response show a high, homogeneous correlation with the stimulus function on both hemispheres. After eigenfiltering, most channels have low CC values, except for a few channels on the contralateral side. Generally, the original response to passive stimulation has smaller CC values than that of tapping stimulation, due to the fact that the latter stimulus will cause more stimulus-locked response.

Figure 10

Spatial maps of the CC values of subject A.

Figure 11

Spatial maps of the CC values of subject B.

4.

Discussion

Although these results are not reported here, projecting onto baseline spatial eigenvectors gives much better results than projections onto the spatial eigenvectors of the stimulus data matrix itself. We speculate that the reason is that during stimulus, the physiological interference is time-locked with the frequency of stimulation, making it hard to separate the interference subspace from the induced-activation subspace, whereas the baseline data contain spatial patterns of the interference only. The spatial map of eigenvectors reinforces that speculation. Looking at Fig. 4, one can see in the baseline data, the first eigenvalue only accounts for 63.96% of the sum of all eigenvalues, whereas for the stimulus data it is 84.11%. This suggests that in the stimulus data, the first eigenvector was a combination of activation and interference, coupled to each other. Similar observations were true for all cases studied. The percentage of the total eigenvalue sum taken by the first three eigenvalues for the baseline data and for averages over both handed stimulus data are presented in Table 1. One can see that in both tapping and passive tactile stimulus, the relative size of the first eigenvalue is much higher than that of the baseline data.

Table 1

Both the baseline and stimulus eigenvectors show spatial maps of “noisy” channels. This was true for all cases studied and thus also gives an indication of how many eigenvectors should be filtered out: we should remove enough eigenvectors so that the global physiological interference is reduced and the noisy channels are removed. A contraindicating consideration is the remaining energy; projecting out too many eigenvectors will leave only a very weak response, which is generally not desirable.

Observing all time course figures, one can see that the filtered results show a deactivation response, i.e., decrease of Δ[HbO] for the ipsilateral channels. This observation was true for all cases studied and may result from nonorthogonality between the interference subspace and the ipsilateral pattern. The physiological reason is not clear yet, but it agrees with some results from a recent fMRI study,²⁶ where negative BOLD responses to motor stimulation in the ipsilateral sensorimotor cortex and subcortical regions were observed.

During the finger-tapping task, an increase in the heart rate synchronous with the stimulus was observed.¹⁸ Because oxyhemoglobin concentration is mostly representative of the arterial compartment, it is more affected by these changes compared to deoxyhemoglobin concentration, which comes mostly from the venous compartment.¹⁸ In tapping stimulation, a more localized activation-induced response in the Δ[Hb] map with respect to the Δ[HbO] map was reported.⁵ Thus, spatial eigenfiltering should show more improvement in localization of Δ[HbO] than for Δ[Hb]. Processing results of all cases studied here verified this supposition.

In the work of this paper, we performed spatial eigenfiltering directly on hemoglobin concentration data, rather than eigenfiltering the absorption coefficient data and then transforming to chromophore concentrations. The reason is that this transformation (i.e., the matrix of extinction coefficients of HbO and Hb at two wavelengths) is not orthogonal, and thus will make the eigenvectors in the absorption coefficient domain lose their orthogonality in the concentration domain. Since we are interested in recovering concentrations, it is more effective to filter in the space domain of the concentration results.

Basically, what we are making in this paper is a kind of stationarity assumption, which enables us to estimate the spatial correlation matrix by averaging over time. Note the temporal behavior of the interference subspace will be variant, but the spatial distribution displays relative stability. A possible way to improve this would be a temporally adaptive (adaptive-filtering-like) method, i.e., one could update the spatial eigenvectors in time. However our experimental experience has been that the baseline behavior is quite stable and repeatable.

Besides the signal processing methods to reduce the interferences, note also that there is in principle the possibility of designing a new hardware system that includes short-distance light channels to estimate the systemic interference, where the superficial variation or nonactivation background signals will dominate the measurement, and then can be removed from other measurements. In this paper, we have only relatively long distance channels in the system, such as 3- and 4-cm separation channels, thus a pure signal processing method is used.

5.

Conclusions and Future Work

We presented a new method using spatial eigenfiltering to separate an activity-evoked response from physiological interference in diffuse optical imaging data. The basic assumption is that the first several spatial eigenvectors of the baseline concentration correlation matrices are dominated by interference patterns and hold most of the total energy. By projecting the stimulus data onto the orthogonal nullspace of these eigenvectors, we can obtain a more localized response, and improved CNR and CC maps.

Future work includes the investigation of optimal criteria for the selection of which source-detector channels to use and the development of combined spatial and temporal filtering. The possible considerations might include making the spatial approach temporally adaptive, which means one can update the eigenvectors in time; or using simultaneously recorded auxiliary signals (cardiac, respiration, blood pressure) to make a temporal model as well. All of this physiologically reasonable prior knowledge in space or time can be put into a joint penalty function.